diff --git a/apps/src/Hummbee/hummbee.cc b/apps/src/Hummbee/hummbee.cc index 168cfe90..a0bceee4 100644 --- a/apps/src/Hummbee/hummbee.cc +++ b/apps/src/Hummbee/hummbee.cc @@ -228,7 +228,7 @@ float Hummbee( string& imageName, string& modelImageName, string& deconvolver, vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain,float& hogbomgain,bool& autohogbom, float& autotrigger,float& autopower, int& autoxmask, int& autoymask,float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, @@ -286,6 +286,21 @@ float Hummbee( decPars_p.interactive=false; decPars_p.autoAdjust=False; //genie decPars_p.fusedThreshold = fusedthreshold; + decPars_p.waveletScales = Vector(waveletscales); + decPars_p.waveletAmps = Vector(waveletamps); + decPars_p.autoThreshold = autothreshold; + decPars_p.autoMaxiter = automaxiter; + decPars_p.autoGain = autogain; + decPars_p.hogbomGain = hogbomgain; + decPars_p.autoHogbom = autohogbom; + decPars_p.autoTrigger = autotrigger; + decPars_p.autoPower = autopower; + decPars_p.autoXMask = autoxmask; + decPars_p.autoYMask = autoymask; + decPars_p.lbfgsEpsF = lbfgsepsf; + decPars_p.lbfgsEpsX = lbfgsepsx; + decPars_p.lbfgsEpsG = lbfgsepsg; + decPars_p.lbfgsMaxit = lbfgsmaxit; decPars_p.specmode=specmode; //deconvolve task does not have this decPars_p.largestscale = largestscale; decPars_p.scalebias = 0.0; @@ -320,9 +335,9 @@ float Hummbee( float minpsffraction = 0.05; float maxpsffraction = 0.8; - float CycleFactor = 1.0; + //float CycleFactor = 1.0; - Float psffraction = MaxPsfSidelobe * CycleFactor; + Float psffraction = MaxPsfSidelobe * cyclefactor;//CycleFactor; psffraction = casacore::max(psffraction, minpsffraction); psffraction = casacore::min(psffraction, maxpsffraction); Float cyclethreshold = PeakResidual * psffraction; diff --git a/apps/src/Hummbee/hummbee.h b/apps/src/Hummbee/hummbee.h index 4d05f7af..9eef58d5 100644 --- a/apps/src/Hummbee/hummbee.h +++ b/apps/src/Hummbee/hummbee.h @@ -45,6 +45,8 @@ #include #include #include +#include +#include #include #include @@ -81,7 +83,7 @@ using namespace casacore; float Hummbee(std::string& imageName, std::string& modelImageName, std::string& deconvolver, std::vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain, float& hogbomgain,bool& autohogbom,float& autotrigger,float& autopower,int& autoxmask,int& autoymask,float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, @@ -95,7 +97,7 @@ void UI(bool restart, int argc, char **argv, bool interactive, string& imageName, string& modelImageName, string& deconvolver, vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain, float& hogbomgain,bool& autohogbom, float& autotrigger,float& autopower,int& autoxmask,int& autoymask,float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, diff --git a/apps/src/Hummbee/hummbee_cl_interface.cc b/apps/src/Hummbee/hummbee_cl_interface.cc index 697ec726..7f049c52 100644 --- a/apps/src/Hummbee/hummbee_cl_interface.cc +++ b/apps/src/Hummbee/hummbee_cl_interface.cc @@ -63,7 +63,7 @@ void UI(bool restart, int argc, char **argv, bool interactive, Float& pbLimit,*/ string& deconvolver, vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain,float& hogbomgain, bool& autohogbom, float& autotrigger,float& autopower, int& autoxmask, int& autoymask, float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, @@ -129,6 +129,21 @@ void UI(bool restart, int argc, char **argv, bool interactive, //InitMap(watchPoints,exposedKeys); exposedKeys.push_back("largestscale"); exposedKeys.push_back("fusedthreshold"); + exposedKeys.push_back("waveletscales"); + exposedKeys.push_back("waveletamps"); + exposedKeys.push_back("autothreshold"); + exposedKeys.push_back("automaxiter"); + exposedKeys.push_back("autogain"); + exposedKeys.push_back("hogbomgain"); + exposedKeys.push_back("autohogbom"); + exposedKeys.push_back("autotrigger"); + exposedKeys.push_back("autopower"); + exposedKeys.push_back("autoxmask"); + exposedKeys.push_back("autoymask"); + exposedKeys.push_back("lbfgsepsf"); + exposedKeys.push_back("lbfgsepsx"); + exposedKeys.push_back("lbfgsepsg"); + exposedKeys.push_back("lbfgsmaxit"); watchPoints["asp"]=exposedKeys; i=1;clgetSValp("deconvolver", deconvolver, i ,watchPoints); @@ -139,6 +154,21 @@ void UI(bool restart, int argc, char **argv, bool interactive, i=1;clgetFValp("largestscale", largestscale,i); i=1;clgetFValp("fusedthreshold", fusedthreshold,i); + N=0; N=clgetNFValp("waveletscales", waveletscales, N); + N=0; N=clgetNFValp("waveletamps", waveletamps, N); + i=1;clgetFValp("autothreshold", autothreshold,i); + i=1;clgetIValp("automaxiter", automaxiter,i); + i=1;clgetFValp("autogain", autogain,i); + i=1;clgetFValp("hogbomgain", hogbomgain,i); + i=1;clgetBValp("autohogbom", autohogbom,i); + i=1;clgetFValp("autotrigger", autotrigger, i); + i=1;clgetFValp("autopower", autopower, i); + i=1;clgetIValp("autoxmask", autoxmask,i); + i=1;clgetIValp("autoymask", autoymask,i); + i=1;clgetFValp("lbfgsepsf", lbfgsepsf,i); + i=1;clgetFValp("lbfgsepsx", lbfgsepsx,i); + i=1;clgetFValp("lbfgsepsg", lbfgsepsg,i); + i=1;clgetIValp("lbfgsmaxit", lbfgsmaxit,i); i=1;clgetIValp("nterms", nterms,i); i=1;clgetFValp("gain", gain,i); @@ -203,6 +233,21 @@ int main(int argc, char **argv) vector scales; float largestscale = -1; float fusedthreshold = 0; + vector waveletscales; + vector waveletamps; + float autothreshold = 0; + int automaxiter = 0; + float autogain = 0; + float hogbomgain = 0; + bool autohogbom = false; + float autotrigger = 0; + float autopower = 1; + int autoxmask = 0; + int autoymask = 0; + float lbfgsepsf = 0.001; + float lbfgsepsx = 0.001; + float lbfgsepsg = 0.001; + int lbfgsmaxit = 5; int nterms=1; float gain=0.1; float threshold=0.0; @@ -225,7 +270,7 @@ int main(int argc, char **argv) ,doPBCorr, conjBeams, pbLimit,*/ deconvolver, scales, - largestscale, fusedthreshold, + largestscale, fusedthreshold,waveletscales, waveletamps, autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower,autoxmask,autoymask,lbfgsepsf,lbfgsepsx,lbfgsepsg,lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -241,7 +286,7 @@ int main(int argc, char **argv) doPBCorr, conjBeams, pbLimit,*/ deconvolver, scales, - largestscale, fusedthreshold, + largestscale, fusedthreshold,waveletscales, waveletamps, autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower,autoxmask,autoymask,lbfgsepsf,lbfgsepsx,lbfgsepsg,lbfgsmaxit, nterms, gain, threshold, nsigma, diff --git a/apps/src/Hummbee/hummbee_py_interface.cc b/apps/src/Hummbee/hummbee_py_interface.cc index 9f41dcae..39796e4e 100644 --- a/apps/src/Hummbee/hummbee_py_interface.cc +++ b/apps/src/Hummbee/hummbee_py_interface.cc @@ -32,6 +32,21 @@ PYBIND11_MODULE(hummbee2py, m) "scales"_a, "largestscale"_a=-1, "fusedthreshold"_a=0, + "waveletscales"_a, + "waveletamps"_a, + "autothreshold"_a=0, + "automaxiter"_a=0, + "autogain"_a=0, + "hogbomgain"_a=0, + "autohogbom"_a=false, + "autotrigger"_a=0, + "autopower"_a=1, + "autoxmask"_a=0, + "autoymask"_a=0, + "lbfgsepsf"_a=0.001, + "lbfgsepsx"_a=0.001, + "lbfgsepsg"_a=0.001, + "lbfgsmaxit"_a=5, "nterms"_a=1, "gain"_a=0.1, "threshold"_a=0.0, diff --git a/apps/src/tests/test_hummbee.cc b/apps/src/tests/test_hummbee.cc index e9512789..0cdf573f 100644 --- a/apps/src/tests/test_hummbee.cc +++ b/apps/src/tests/test_hummbee.cc @@ -52,6 +52,21 @@ TEST(HummbeeTest, AppLevelCubeAsp) { vector scales; float largestscale = -1; float fusedthreshold = 0; + vector waveletscales; + vector waveletamps; + float autothreshold; + int automaxiter; + float autogain; + float hogbomgain; + bool autohogbom; + float autotrigger; + float autopower; + int autoxmask; + int autoymask; + float lbfgsepsf=0.001; + float lbfgsepsx=0.001; + float lbfgsepsg=0.001; + int lbfgsmaxit=5; int nterms=1; float gain=0.1; float threshold=1e-4; @@ -68,6 +83,9 @@ TEST(HummbeeTest, AppLevelCubeAsp) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -134,6 +152,21 @@ TEST(HummbeeTest, AppLevelMfsAsp) { vector scales; float largestscale = -1; float fusedthreshold = 0.007; + vector waveletscales; + vector waveletamps; + float autothreshold; + int automaxiter; + float autogain; + float hogbomgain; + bool autohogbom; + float autotrigger; + float autopower; + int autoxmask; + int autoymask; + float lbfgsepsf=0.001; + float lbfgsepsx=0.001; + float lbfgsepsg=0.001; + int lbfgsmaxit=5; int nterms=1; float gain=0.2; float threshold=2.6e-07; @@ -152,6 +185,9 @@ TEST(HummbeeTest, AppLevelMfsAsp) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -222,6 +258,21 @@ TEST(HummbeeTest, AppLevelWAsp) { vector scales; float largestscale = 5; float fusedthreshold = 0.05; + vector waveletscales; + vector waveletamps; + float autothreshold; + int automaxiter; + float autogain; + float hogbomgain; + bool autohogbom; + float autotrigger; + float autopower; + int autoxmask; + int autoymask; + float lbfgsepsf=0.001; + float lbfgsepsx=0.001; + float lbfgsepsg=0.001; + int lbfgsmaxit=5; int nterms=3; float gain=0.6; float threshold=0.2; @@ -236,6 +287,9 @@ TEST(HummbeeTest, AppLevelWAsp) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -292,6 +346,21 @@ TEST(HummbeeTest, AppLevelMfsRestore) { vector scales; float largestscale = -1; float fusedthreshold = 0.007; + vector waveletscales; + vector waveletamps; + float autothreshold; + int automaxiter; + float autogain; + float hogbomgain; + bool autohogbom; + float autotrigger; + float autopower; + int autoxmask; + int autoymask; + float lbfgsepsf=0.001; + float lbfgsepsx=0.001; + float lbfgsepsg=0.001; + int lbfgsmaxit=5; int nterms=1; float gain=0.2; float threshold=2.6e-07; @@ -311,6 +380,9 @@ TEST(HummbeeTest, AppLevelMfsRestore) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -339,6 +411,9 @@ TEST(HummbeeTest, AppLevelMfsRestore) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -380,6 +455,21 @@ TEST(HummbeeTest, UIFactory) { vector scales; float largestscale = -1; float fusedthreshold = 0; + vector waveletscales; + vector waveletamps; + float autothreshold; + int automaxiter; + float autogain; + float hogbomgain; + bool autohogbom; + float autotrigger; + float autopower; + int autoxmask; + int autoymask; + float lbfgsepsf=0.001; + float lbfgsepsx=0.001; + float lbfgsepsg=0.001; + int lbfgsmaxit=5; int nterms=1; float gain=0.1; float threshold=0.0; @@ -396,6 +486,9 @@ TEST(HummbeeTest, UIFactory) { deconvolver, scales, largestscale, fusedthreshold, + waveletscales, waveletamps, + autothreshold, automaxiter, autogain, hogbomgain, autohogbom, autotrigger, autopower, autoxmask, autoymask, + lbfgsepsf, lbfgsepsx, lbfgsepsg, lbfgsmaxit, nterms, gain, threshold, nsigma, diff --git a/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.cc b/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.cc index e1225670..2fe9624a 100644 --- a/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.cc +++ b/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.cc @@ -61,18 +61,25 @@ using namespace casacore; namespace casa { //# NAMESPACE CASA - BEGIN - SDAlgorithmAAspClean::SDAlgorithmAAspClean(Float fusedThreshold, bool isSingle, Int largestScale, Int stoppointmode): + SDAlgorithmAAspClean::SDAlgorithmAAspClean(Vector scales, Float hogbomGain, Float fusedThreshold, bool isSingle, Int largestScale, Int stoppointmode, Float lbfgsEpsF, Float lbfgsEpsX, Float lbfgsEpsG, Int lbfgsMaxit): SDAlgorithmBase(), itsMatPsf(), itsMatResidual(), itsMatModel(), itsCleaner(), itsStopPointMode(stoppointmode), - itsFusedThreshold(fusedThreshold), - itsUserLargestScale(largestScale), itsMCsetup(true), + itsFusedThreshold(fusedThreshold), + itsScales(scales), + itsHogbomGain(hogbomGain), + itsLbfgsEpsF(lbfgsEpsF), + itsLbfgsEpsX(lbfgsEpsX), + itsLbfgsEpsG(lbfgsEpsG), + itsLbfgsMaxit(lbfgsMaxit), itsPrevPsfWidth(0), - itsIsSingle(isSingle) + itsIsSingle(isSingle), + itsUserLargestScale(largestScale) { itsAlgorithmName = String("asp"); + LogIO os(LogOrigin("SDAlgorithmAAspClean", "constructor", WHERE)); } SDAlgorithmAAspClean::~SDAlgorithmAAspClean() @@ -109,7 +116,10 @@ namespace casa { //# NAMESPACE CASA - BEGIN if (itsPrevPsfWidth != width) { itsPrevPsfWidth = width; - itsCleaner.setInitScaleXfrs(width); + if (itsScales.size() < 1) + itsCleaner.setInitScaleXfrs(width); + else + itsCleaner.loadInitScaleXfrs(itsScales); } itsCleaner.stopPointMode( itsStopPointMode ); @@ -120,10 +130,18 @@ namespace casa { //# NAMESPACE CASA - BEGIN // Not used. Kept for unit test //Matrix tempMat1(itsMatResidual); //itsCleaner.setOrigDirty( tempMat1 ); - + if (itsHogbomGain < 0) + { + os << LogIO::WARN << "Acceptable hogbomgain values are >= 0. Changing hogbomgain from " << itsHogbomGain << " to 0." << LogIO::POST; + itsHogbomGain = 0.0; + } itsCleaner.setFusedThreshold(itsFusedThreshold); + itsCleaner.setHogbomGain(itsHogbomGain); + } + + itsCleaner.setLBFGSControl(itsLbfgsEpsF,itsLbfgsEpsX,itsLbfgsEpsG,itsLbfgsMaxit); // Parts to be repeated at each minor cycle start.... //itsCleaner.setInitScaleMasks(itsMatMask); //casa6 @@ -135,10 +153,11 @@ namespace casa { //# NAMESPACE CASA - BEGIN tempMat1.reference(itsMatResidual); itsCleaner.setDirty( tempMat1 ); // InitScaleXfrs and InitScaleMasks should already be set - itsScaleSizes.clear(); - itsScaleSizes = itsCleaner.getActiveSetAspen(); - itsScaleSizes.push_back(0.0); // put 0 scale - itsCleaner.defineAspScales(itsScaleSizes); + + itsScaleSizes.clear(); + itsScaleSizes = itsCleaner.getActiveSetAspen(); + itsScaleSizes.push_back(0.0); // put 0 scale + itsCleaner.defineAspScales(itsScaleSizes); } @@ -159,28 +178,29 @@ namespace casa { //# NAMESPACE CASA - BEGIN Matrix prevModel; prevModel = itsMatModel; - //cout << "AAspALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl; + //cout << "AAspALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl; + + // retval + // 1 = converged + // 0 = not converged but behaving normally + // -1 = not converged and stopped on cleaning consecutive smallest scale + // -2 = not converged and either large scale hit negative or diverging + // -3 = clean is diverging rather than converging + itsCleaner.startingIteration( 0 ); + Int retval = itsCleaner.aspclean( tempModel ); + iterdone = itsCleaner.numberIterations(); - // retval - // 1 = converged - // 0 = not converged but behaving normally - // -1 = not converged and stopped on cleaning consecutive smallest scale - // -2 = not converged and either large scale hit negative or diverging - // -3 = clean is diverging rather than converging - itsCleaner.startingIteration( 0 ); - Int retval = itsCleaner.aspclean( tempModel ); - iterdone = itsCleaner.numberIterations(); + if( retval==-1 ) {os << LogIO::WARN << "AspClean minor cycle stopped on cleaning consecutive smallest scale" << LogIO::POST; } + if( retval==-2 ) {os << LogIO::WARN << "AspClean minor cycle stopped at large scale negative or diverging" << LogIO::POST;} + if( retval==-3 ) {os << LogIO::WARN << "AspClean minor cycle stopped because it is diverging" << LogIO::POST; } - if( retval==-1 ) {os << LogIO::WARN << "AspClean minor cycle stopped on cleaning consecutive smallest scale" << LogIO::POST; } - if( retval==-2 ) {os << LogIO::WARN << "AspClean minor cycle stopped at large scale negative or diverging" << LogIO::POST;} - if( retval==-3 ) {os << LogIO::WARN << "AspClean minor cycle stopped because it is diverging" << LogIO::POST; } + // update residual - this is critical + itsMatResidual = itsCleaner.getterResidual(); - // update residual - this is critical - itsMatResidual = itsCleaner.getterResidual(); + peakresidual = itsCleaner.getterPeakResidual(); + //cout << "SDAlg: peakres " << peakresidual << endl; + modelflux = sum( itsMatModel ); - peakresidual = itsCleaner.getterPeakResidual(); - //cout << "SDAlg: peakres " << peakresidual << endl; - modelflux = sum( itsMatModel ); } void SDAlgorithmAAspClean::finalizeDeconvolver() diff --git a/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.h b/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.h index b4e33458..ca6fc6c3 100644 --- a/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.h +++ b/src/synthesis/ImagerObjects/SDAlgorithmAAspClean.h @@ -52,7 +52,7 @@ namespace casa { //# NAMESPACE CASA - BEGIN public: // Empty constructor - SDAlgorithmAAspClean(casacore::Float fusedThreshold = 0.0, bool isSingle = true, casacore::Int largestScale = -1, casacore::Int stoppointmode = -1); + SDAlgorithmAAspClean(casacore::Vector scales, casacore::Float hogbomGain, casacore::Float fusedThreshold = 0.0, bool isSingle = true, casacore::Int largestScale = -1, casacore::Int stoppointmode = -1, casacore::Float lbfgsEpsF= 0.001, casacore::Float lbfgsEpsX= 0.001, casacore::Float lbfgsEpsG= 0.001, casacore::Int lbfgsMaxit= 5); virtual ~SDAlgorithmAAspClean(); protected: @@ -83,6 +83,13 @@ namespace casa { //# NAMESPACE CASA - BEGIN casacore::Int itsStopPointMode; casacore::Float itsFusedThreshold; casacore::Int itsUserLargestScale; + + casacore::Vector itsScales; + casacore::Float itsHogbomGain; + casacore::Float itsLbfgsEpsF; + casacore::Float itsLbfgsEpsX; + casacore::Float itsLbfgsEpsG; + casacore::Int itsLbfgsMaxit; /* casacore::IPosition itsMaxPos; diff --git a/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.cc b/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.cc new file mode 100644 index 00000000..0e6e93ea --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.cc @@ -0,0 +1,293 @@ +//# SDAlgorithmMSClean.cc: Implementation of SDAlgorithmMSClean classes +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# $Id$ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//new +#include + +#include + +#include + +#include +#include +#include + +#include +#include + + +using namespace casacore; +namespace casa { //# NAMESPACE CASA - BEGIN + + + SDAlgorithmAutoClean::SDAlgorithmAutoClean(Float autoThreshold, Int autoMaxiter, Float autoGain, Float hogbomGain, Bool autoHogbom, Float autoTrigger, Float autoPower, Int autoXMask, Int autoYMask): + SDAlgorithmBase(), + itsMatPsf(), itsMatResidual(), itsMatModel(), itsMatTemp(), + itsCleaner(), + itsautoThreshold(autoThreshold), itsautoMaxiter(autoMaxiter), itsautoGain(autoGain), itshogbomGain(hogbomGain), itsautoHogbom(autoHogbom), itsautoTrigger(autoTrigger), itsautoPower(autoPower), itsautoXMask(autoXMask), itsautoYMask(autoYMask) + { + itsAlgorithmName=String("autocorrelation"); + } + + + SDAlgorithmAutoClean::~SDAlgorithmAutoClean() + { + + } + + Long SDAlgorithmAutoClean::estimateRAM(const vector& imsize){ + Long mem=0; + IPosition shp; + if(itsImages) + shp=itsImages->getShape(); + else if(imsize.size() >1) + shp=IPosition(imsize); + else + return 0; + //throw(AipsError("Deconvolver cannot estimate the memory usage at this point")); + + //Number of planes in memory + //npsf=nscales+1 + //nresidual=4+nscales + //nmodel=1 + //nmasks=2+nscales + //transfer functions=nscales*2 + Long nplanes=5*1+7; + //in kB + mem=sizeof(Float)*(shp(0))*(shp(1))*nplanes/1024; + + return mem; + } + + // void SDAlgorithmMSClean::initializeDeconvolver( Float &peakresidual, Float &modelflux ) + void SDAlgorithmAutoClean::initializeDeconvolver() + { + LogIO os( LogOrigin("SDAlgorithmAutoClean","initializeDeconvolver",WHERE) ); + AlwaysAssert( (bool) itsImages, AipsError ); + { + LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Read); + LatticeLocker lock2 (*(itsImages->model()), FileLocker::Read); + LatticeLocker lock3 (*(itsImages->psf()), FileLocker::Read); + LatticeLocker lock4 (*(itsImages->mask()), FileLocker::Read); + //LatticeLocker lock5 (*(itsImages->tempworkimage()), FileLocker::Read); + (itsImages->residual())->get( itsMatResidual , true ); + (itsImages->model())->get( itsMatModel , true ); + (itsImages->psf())->get( itsMatPsf , true ); + itsImages->mask()->get( itsMatMask, true ); + + } + //// Initialize the AutoCleaner. + /// ----------- do once ---------- + { + itsCleaner.ignoreCenterBox( true ); // Clean full image + + Matrix tempMat; + tempMat.reference( itsMatPsf ); + itsCleaner.setPsf( tempMat ); + itsCleaner.setscales(); + + } + /// ----------------------------------------- + + /* + /// Find initial max vals.. + findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos ); + itsModelFlux = sum( itsMatModel ); + + peakresidual = itsPeakResidual; + modelflux = itsModelFlux; + */ + + // Parts to be repeated at each minor cycle start.... + + itsCleaner.setcontrol(CleanEnums::MULTISCALE,0,0,0);/// Needs to come before makeDirtyScales + + Matrix tempmask(itsMatMask); + itsCleaner.setMask( tempmask ); + + Matrix tempMat1; + tempMat1.reference( itsMatResidual ); + itsCleaner.setDirty( tempMat1 ); + + //itsCleaner.initializeCorrProducts(); + //itsCleaner.approximateBasisFunction(); + + if (itsautoThreshold < 0) + { + os << LogIO::WARN << "Acceptable autothreshld values are >= 0. Changing autothreshold from " << itsautoThreshold << " to 0." << LogIO::POST; + itsautoThreshold = 0.0; + } + if (itsautoMaxiter < 0) + { + os << LogIO::WARN << "Acceptable automaxiter values are >= 0. Changing automaxiter from " << itsautoMaxiter << " to 0." << LogIO::POST; + itsautoMaxiter = 0.0; + } + if (itsautoGain < 0) + { + os << LogIO::WARN << "Acceptable autogain values are >= 0. Changing autogain from " << itsautoGain << " to 0." << LogIO::POST; + itsautoGain = 0.0; + } + if (itshogbomGain < 0) + { + os << LogIO::WARN << "Acceptable hogbomgain values are >= 0. Changing hogbomgain from " << itshogbomGain << " to 0." << LogIO::POST; + itshogbomGain = 0.0; + } + if (itsautoPower < 0) + { + os << LogIO::WARN << "Acceptable autopower values are >= 0. Changing autopower from " << itsautoPower << " to 0." << LogIO::POST; + itsautoPower = 1.0; + } + + itsCleaner.setAutoControl(itsautoThreshold, itsautoMaxiter, itsautoGain, itshogbomGain, itsautoHogbom, itsautoTrigger, itsautoPower, itsautoXMask, itsautoYMask); + + if(itsautoHogbom==false){ + if(itsImages->doesImageExist(itsImages->getName()+String(".tildeII")) && + itsImages->doesImageExist(itsImages->getName()+String(".cmap")) && + itsImages->doesImageExist(itsImages->getName()+String(".MtildeMB")) && + itsImages->doesImageExist(itsImages->getName()+String(".mod")) && + itsImages->doesImageExist(itsImages->getName()+String(".tildeMI")) && + itsImages->doesImageExist(itsImages->getName()+String(".BI")) && + itsImages->doesImageExist(itsImages->getName()+String(".tildeMB")) && + itsImages->doesImageExist(itsImages->getName()+String(".tildeMBB")) && + itsImages->doesImageExist(itsImages->getName()+String(".tildeMBI")) && + itsImages->doesImageExist(itsImages->getName()+String(".MtildeMBB")) && + itsImages->doesImageExist(itsImages->getName()+String(".BB")) && + itsImages->doesImageExist(itsImages->getName()+String(".window_basis")) + ){ + os << "Autocorrelation product exist, load autocorrelation products" << LogIO::POST; + itsImages->loadmatrix(itsCleaner.tildeII, itsImages->getName()+String(".tildeII")); + itsImages->loadmatrix(itsCleaner.cmap, itsImages->getName()+String(".cmap")); + itsImages->loadmatrix(itsCleaner.MtildeMB, itsImages->getName()+String(".MtildeMB")); + itsImages->loadmatrix(itsCleaner.mod, itsImages->getName()+String(".mod")); + itsImages->loadmatrix(itsCleaner.tildeMI, itsImages->getName()+String(".tildeMI")); + itsImages->loadmatrix(itsCleaner.BI, itsImages->getName()+String(".BI")); + itsImages->loadmatrix(itsCleaner.tildeMB, itsImages->getName()+String(".tildeMB")); + itsImages->loadmatrix(itsCleaner.tildeMBB, itsImages->getName()+String(".tildeMBB")); + itsImages->loadmatrix(itsCleaner.tildeMBI, itsImages->getName()+String(".tildeMBI")); + itsImages->loadmatrix(itsCleaner.MtildeMBB, itsImages->getName()+String(".MtildeMBB")); + itsImages->loadmatrix(itsCleaner.BB, itsImages->getName()+String(".BB")); + itsImages->loadmatrix(itsCleaner.window_basis, itsImages->getName()+String(".window_basis")); + + itsCleaner.reinitializeCorrProducts(); + } + else{ + itsCleaner.initializeCorrProducts(); + itsCleaner.approximateBasisFunction(); + } + } + } + + void SDAlgorithmAutoClean::takeOneStep( Float loopgain, Int cycleNiter, Float cycleThreshold, Float &peakresidual, Float &modelflux, Int &iterdone) + { + LogIO os( LogOrigin("SDAlgorithmAutoClean","takeOneStep",WHERE) ); + + Quantity thresh( cycleThreshold, "Jy" ); + // Quantity ftthresh( 100.0, "Jy" ); /// Look at MFMSCleanImageSkyModel.cc for more. + itsCleaner.setcontrol(CleanEnums::MULTISCALE, cycleNiter, loopgain, thresh); //, ftthresh); + + Matrix tempModel; + tempModel.reference( itsMatModel ); + //save the previous model + Matrix prevModel; + prevModel=itsMatModel; + + //cout << "SDALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl; + + // retval + // 1 = converged + // 0 = not converged but behaving normally + // -1 = not converged and stopped on cleaning consecutive smallest scale + // -2 = not converged and either large scale hit negative or diverging + // -3 = clean is diverging rather than converging + itsCleaner.startingIteration( 0 ); + Int retval = itsCleaner.clean( tempModel ); + iterdone = itsCleaner.numberIterations(); + + // Retrieve residual to be saved to the .residual file in finalizeDeconvolver. + ////This is going to be wrong if there is no 0 scale; + ///Matrix residual(itsCleaner.residual()); + //Matrix residual(itsCleaner.residual(tempModel-prevModel)); + // cout << "Max tempModel : " << max(abs(tempModel)) << " Max prevModel : " << max(abs(prevModel)) << endl; + //itsMatResidual = itsCleaner.residual(tempModel-prevModel); + itsMatResidual = itsCleaner.residual(); + //itsMatTemp = itsCleaner.cmap; + + // account for mask as well + //peakresidual = max(abs(residual*itsMatMask)); + peakresidual = max(abs(itsMatResidual*itsMatMask)); + modelflux = sum( itsMatModel ); // Performance hog ? + } + + void SDAlgorithmAutoClean::finalizeDeconvolver() + { + ///MatrixCleaner does not modify the original residual image matrix + ///so the first line is a dummy. + LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Write); + LatticeLocker lock2 (*(itsImages->model()), FileLocker::Write); + //LatticeLocker lock5 (*(itsImages->tempworkimage()), FileLocker::Write); + (itsImages->residual())->put( itsMatResidual ); + (itsImages->model())->put( itsMatModel ); + //(itsImages->tempworkimage())->put( itsMatTemp ); + itsImages->savematrix(itsCleaner.tildeII, itsImages->getName()+String(".tildeII")); + itsImages->savematrix(itsCleaner.cmap, itsImages->getName()+String(".cmap")); + itsImages->savematrix(itsCleaner.MtildeMB, itsImages->getName()+String(".MtildeMB")); + itsImages->savematrix(itsCleaner.mod, itsImages->getName()+String(".mod")); + itsImages->savematrix(itsCleaner.tildeMI, itsImages->getName()+String(".tildeMI")); + itsImages->savematrix(itsCleaner.BI, itsImages->getName()+String(".BI")); + itsImages->savematrix(itsCleaner.tildeMB, itsImages->getName()+String(".tildeMB")); + itsImages->savematrix(itsCleaner.tildeMBB, itsImages->getName()+String(".tildeMBB")); + itsImages->savematrix(itsCleaner.tildeMBI, itsImages->getName()+String(".tildeMBI")); + itsImages->savematrix(itsCleaner.MtildeMBB, itsImages->getName()+String(".MtildeMBB")); + itsImages->savematrix(itsCleaner.BB, itsImages->getName()+String(".BB")); + itsImages->savematrix(itsCleaner.window_basis, itsImages->getName()+String(".window_basis")); + } + +} //# NAMESPACE CASA - END + diff --git a/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.h b/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.h new file mode 100644 index 00000000..17d43d68 --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmAutoClean.h @@ -0,0 +1,91 @@ +//# SDAlgorithmMSClean.h: Definition for SDAlgorithmMSClean +//# Copyright (C) 1996,1997,1998,1999,2000,2002 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be adressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id$ + +#ifndef SYNTHESIS_SDALGORITHMAutoCLEAN_H +#define SYNTHESIS_SDALGORITHMAutoCLEAN_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + + /* Forware Declaration */ + class SIMinorCycleController; + + + class SDAlgorithmAutoClean : public SDAlgorithmBase + { + public: + + // Empty constructor + SDAlgorithmAutoClean(casacore::Float autoThreshold, casacore::Int autoMaxiter, casacore::Float autoGain, casacore::Float hogbomGain, casacore::Bool autoHogbom, casacore::Float autoTrigger, casacore::Float autoPower, casacore::Int autoXMask, casacore::Int autoYMask); + virtual ~SDAlgorithmAutoClean(); + + //returns the estimate of memory used in kilobytes (kB); + virtual casacore::Long estimateRAM( const std::vector& imsize); + protected: + + // Local functions to be overloaded by various algorithm deconvolvers. + void takeOneStep( casacore::Float loopgain, casacore::Int cycleNiter, casacore::Float cycleThreshold, + casacore::Float &peakresidual, casacore::Float &modelflux, casacore::Int &iterdone ); + void initializeDeconvolver(); + void finalizeDeconvolver(); + + casacore::Array itsMatPsf, itsMatResidual, itsMatModel; + casacore::Array itsMatMask; // Make an array if we eventually use multi-term masks... + casacore::Array itsMatTemp; + + casacore::Float itsautoThreshold; + casacore::Int itsautoMaxiter; + casacore::Float itsautoGain; + casacore::Float itshogbomGain; + casacore::Bool itsautoHogbom; + casacore::Float itsautoTrigger; + casacore::Float itsautoPower; + casacore::Int itsautoXMask; + casacore::Int itsautoYMask; + + AutoCleaner itsCleaner; + + private: + //casacore::Bool itsMCsetup; + + }; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.cc b/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.cc new file mode 100644 index 00000000..1ec65ff4 --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.cc @@ -0,0 +1,202 @@ +//# SDAlgorithmAAspClean.cc: Implementation of SDAlgorithmAAspClean classes +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# $Id$ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +using namespace casacore; +namespace casa { //# NAMESPACE CASA - BEGIN + + SDAlgorithmSpectralAspClean::SDAlgorithmSpectralAspClean(Vector scales, Float hogbomGain, Float fusedThreshold, bool isSingle, Int largestScale, Int stoppointmode, Float lbfgsEpsF, Float lbfgsEpsX, Float lbfgsEpsG, Int lbfgsMaxit): + SDAlgorithmBase(), + itsMatPsf(), itsMatResidual(), itsMatModel(), + itsCleaner(), + itsStopPointMode(stoppointmode), + itsMCsetup(true), + itsFusedThreshold(fusedThreshold), + itsScales(scales), + itsHogbomGain(hogbomGain), + itsLbfgsEpsF(lbfgsEpsF), + itsLbfgsEpsX(lbfgsEpsX), + itsLbfgsEpsG(lbfgsEpsG), + itsLbfgsMaxit(lbfgsMaxit), + itsPrevPsfWidth(0), + itsIsSingle(isSingle), + itsUserLargestScale(largestScale) + { + itsAlgorithmName = String("spectral_asp"); + LogIO os(LogOrigin("SDAlgorithmSpectralAspClean", "constructor", WHERE)); + } + + SDAlgorithmSpectralAspClean::~SDAlgorithmSpectralAspClean() + { + + } + + void SDAlgorithmSpectralAspClean::initializeDeconvolver() + { + LogIO os(LogOrigin("SDAlgorithmSpectralAspClean", "initializeDeconvolver", WHERE)); + AlwaysAssert((bool)itsImages, AipsError); + + itsImages->residual()->get( itsMatResidual, true ); + itsImages->model()->get( itsMatModel, true ); + itsImages->psf()->get( itsMatPsf, true ); + itsImages->mask()->get( itsMatMask, true ); + + // Initialize the AspMatrixCleaner. + // If it's single channel, this only needs to be computed once. + // Otherwise, it needs to be called repeatedly at each minor cycle start to + // get psf for each channel + if (itsMCsetup) + { + Matrix tempMat(itsMatPsf); + itsCleaner.setPsf(tempMat); + // Initial scales are unchanged and only need to be + // computed when psf width is updated + const Float width = itsCleaner.getPsfGaussianWidth(*(itsImages->psf())); + //if user does not provide the largest scale, we calculate it internally. + itsCleaner.setUserLargestScale(itsUserLargestScale); + // we do not use the shortest baseline approach below b/c of reasons in CAS-940 dated around Jan 2022 + // itsCleaner.getLargestScaleSize(*(itsImages->psf())); + + if (itsPrevPsfWidth != width) + { + itsPrevPsfWidth = width; + if (itsScales.size() < 1) + itsCleaner.setInitScaleXfrs(width); + else + itsCleaner.loadInitScaleXfrs(itsScales); + } + + itsCleaner.stopPointMode( itsStopPointMode ); + itsCleaner.ignoreCenterBox( true ); // Clean full image + // If it's single channel, we do not do the expensive set up repeatedly + if (itsIsSingle) + itsMCsetup = false; + // Not used. Kept for unit test + //Matrix tempMat1(itsMatResidual); + //itsCleaner.setOrigDirty( tempMat1 ); + + if (itsFusedThreshold < 0) + { + os << LogIO::WARN << "Acceptable fusedthreshld values are >= 0. Changing fusedthreshold from " << itsFusedThreshold << " to -1." << LogIO::POST; + itsFusedThreshold = -1.; + } + if (itsHogbomGain < 0) + { + os << LogIO::WARN << "Acceptable hogbomgain values are >= 0. Changing hogbomgain from " << itsHogbomGain << " to 0." << LogIO::POST; + itsHogbomGain = 0.0; + } + + itsCleaner.setFusedThreshold(itsFusedThreshold); + itsCleaner.setHogbomGain(itsHogbomGain); + + } + + itsCleaner.setLBFGSControl(itsLbfgsEpsF,itsLbfgsEpsX,itsLbfgsEpsG,itsLbfgsMaxit); + + // Parts to be repeated at each minor cycle start.... + //itsCleaner.setInitScaleMasks(itsMatMask); //casa6 + Matrix tempMatMask(itsMatMask); + itsCleaner.setInitScaleMasks(tempMatMask); + itsCleaner.setaspcontrol(0, 0, 0, Quantity(0.0, "%"));/// Needs to come before the rest + + Matrix tempMat1; + tempMat1.reference(itsMatResidual); + itsCleaner.setDirty( tempMat1 ); + // InitScaleXfrs and InitScaleMasks should already be set + } + + + void SDAlgorithmSpectralAspClean::takeOneStep( Float loopgain, + Int cycleNiter, + Float cycleThreshold, + Float &peakresidual, + Float &modelflux, + Int &iterdone) + { + LogIO os( LogOrigin("SDAlgorithmSpectralAspClean","takeOneStep", WHERE) ); + + Quantity thresh(cycleThreshold, "Jy"); + itsCleaner.setaspcontrol(cycleNiter, loopgain, thresh, Quantity(0.0, "%")); + Matrix tempModel; + tempModel.reference( itsMatModel ); + //save the previous model + Matrix prevModel; + prevModel = itsMatModel; + + itsCleaner.startingIteration( 0 ); + itsCleaner.MFaspclean( tempModel ); + os << "Aspclean finished" << LogIO::POST; + iterdone = itsCleaner.numberIterations(); + + // update residual - this is critical + itsMatResidual = itsCleaner.getterResidual(); + + peakresidual = itsCleaner.getterPeakResidual(); + //cout << "SDAlg: peakres " << peakresidual << endl; + modelflux = sum( itsMatModel ); + + } + + void SDAlgorithmSpectralAspClean::finalizeDeconvolver() + { + (itsImages->residual())->put( itsMatResidual ); + (itsImages->model())->put( itsMatModel ); + } + +} //# NAMESPACE CASA - END diff --git a/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.h b/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.h new file mode 100644 index 00000000..35fca3b6 --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmSpectralAspClean.h @@ -0,0 +1,108 @@ +//# SDAlgorithmAAspClean.h: Definition for SDAlgorithmAAspClean +//# Copyright (C) 1996,1997,1998,1999,2000,2002 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be adressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id$ + +#ifndef SYNTHESIS_SDALGORITHMSPECTRALASPCLEAN_H +#define SYNTHESIS_SDALGORITHMSPECTRALASPCLEAN_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + + /* Forware Declaration */ + class SIMinorCycleController; + + + class SDAlgorithmSpectralAspClean : public SDAlgorithmBase + { + public: + + // Empty constructor + SDAlgorithmSpectralAspClean(casacore::Vector scales, casacore::Float hogbomGain, casacore::Float fusedThreshold = 0.0, bool isSingle = true, casacore::Int largestScale = -1, casacore::Int stoppointmode = -1, casacore::Float lbfgsEpsF= 0.001, casacore::Float lbfgsEpsX= 0.001, casacore::Float lbfgsEpsG= 0.001, casacore::Int lbfgsMaxit= 5); + virtual ~SDAlgorithmSpectralAspClean(); + + protected: + + // Local functions to be overloaded by various algorithm deconvolvers. + virtual void takeOneStep( casacore::Float loopgain, casacore::Int cycleNiter, casacore::Float cycleThreshold, casacore::Float &peakresidual, casacore::Float &modelflux, casacore::Int &iterdone ); + virtual void initializeDeconvolver(); + virtual void finalizeDeconvolver(); + + + /* + void findNextComponent( casacore::Float loopgain ); + void updateModel(); + void updateResidual(); + */ + + /* + casacore::SubImage itsResidual, itsPsf, itsModel, itsImage; + casacore::Float itsComp; + */ + //casacore::SubImage itsResidual, itsPsf, itsModel, itsImage; + + casacore::Array itsMatPsf, itsMatResidual, itsMatModel; + casacore::Array itsMatMask; // Make an array if we eventually use multi-term masks... + + SpectralAspCleaner itsCleaner; + std::vector itsScaleSizes; + casacore::Int itsStopPointMode; + casacore::Float itsFusedThreshold; + casacore::Int itsUserLargestScale; + + casacore::Vector itsScales; + casacore::Float itsHogbomGain; + casacore::Float itsLbfgsEpsF; + casacore::Float itsLbfgsEpsX; + casacore::Float itsLbfgsEpsG; + casacore::Int itsLbfgsMaxit; + + /* + casacore::IPosition itsMaxPos; + casacore::Float itsPeakResidual; + casacore::Float itsModelFlux; + */ + private: + casacore::Bool itsMCsetup; // if we should do setup or not + casacore::Float itsPrevPsfWidth; + bool itsIsSingle; + + }; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.cc b/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.cc new file mode 100644 index 00000000..4a33b033 --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.cc @@ -0,0 +1,220 @@ +//# SDAlgorithmAAspClean.cc: Implementation of SDAlgorithmAAspClean classes +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# $Id$ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +using namespace casacore; +namespace casa { //# NAMESPACE CASA - BEGIN + + SDAlgorithmWaveletAspClean::SDAlgorithmWaveletAspClean(Vector scales, Float hogbomGain, Vector waveletScales, Vector waveletAmps, Float fusedThreshold, bool isSingle, Int largestScale, Int stoppointmode, Float lbfgsEpsF, Float lbfgsEpsX, Float lbfgsEpsG, Int lbfgsMaxit): + SDAlgorithmBase(), + itsMatPsf(), itsMatResidual(), itsMatModel(), + itsCleaner(), + itsStopPointMode(stoppointmode), + itsMCsetup(true), + itsFusedThreshold(fusedThreshold), + itsScales(scales), + itsHogbomGain(hogbomGain), + itsWaveletScales(waveletScales), + itsWaveletAmps(waveletAmps), + itsLbfgsEpsF(lbfgsEpsF), + itsLbfgsEpsX(lbfgsEpsX), + itsLbfgsEpsG(lbfgsEpsG), + itsLbfgsMaxit(lbfgsMaxit), + itsPrevPsfWidth(0), + itsIsSingle(isSingle), + itsUserLargestScale(largestScale) + { + itsAlgorithmName = String("wavelet_asp"); + LogIO os(LogOrigin("SDAlgorithmWaveletAspClean", "constructor", WHERE)); + } + + SDAlgorithmWaveletAspClean::~SDAlgorithmWaveletAspClean() + { + + } + + void SDAlgorithmWaveletAspClean::initializeDeconvolver() + { + LogIO os(LogOrigin("SDAlgorithmWaveletAspClean", "initializeDeconvolver", WHERE)); + AlwaysAssert((bool)itsImages, AipsError); + + itsImages->residual()->get( itsMatResidual, true ); + itsImages->model()->get( itsMatModel, true ); + itsImages->psf()->get( itsMatPsf, true ); + itsImages->mask()->get( itsMatMask, true ); + + // Initialize the AspMatrixCleaner. + // If it's single channel, this only needs to be computed once. + // Otherwise, it needs to be called repeatedly at each minor cycle start to + // get psf for each channel + if (itsMCsetup) + { + Matrix tempMat(itsMatPsf); + itsCleaner.setPsf(tempMat); + // Initial scales are unchanged and only need to be + // computed when psf width is updated + const Float width = itsCleaner.getPsfGaussianWidth(*(itsImages->psf())); + //if user does not provide the largest scale, we calculate it internally. + itsCleaner.setUserLargestScale(itsUserLargestScale); + // we do not use the shortest baseline approach below b/c of reasons in CAS-940 dated around Jan 2022 + // itsCleaner.getLargestScaleSize(*(itsImages->psf())); + + if (itsPrevPsfWidth != width) + { + itsPrevPsfWidth = width; + if (itsScales.size() < 1) + itsCleaner.setInitScaleXfrs(width); + else + itsCleaner.loadInitScaleXfrs(itsScales); + } + + itsCleaner.stopPointMode( itsStopPointMode ); + itsCleaner.ignoreCenterBox( true ); // Clean full image + // If it's single channel, we do not do the expensive set up repeatedly + if (itsIsSingle) + itsMCsetup = false; + // Not used. Kept for unit test + //Matrix tempMat1(itsMatResidual); + //itsCleaner.setOrigDirty( tempMat1 ); + + if (itsFusedThreshold < 0) + { + os << LogIO::WARN << "Acceptable fusedthreshld values are >= 0. Changing fusedthreshold from " << itsFusedThreshold << " to -1." << LogIO::POST; + itsFusedThreshold = -1.; + } + if (itsHogbomGain < 0) + { + os << LogIO::WARN << "Acceptable hogbomgain values are >= 0. Changing hogbomgain from " << itsHogbomGain << " to 0." << LogIO::POST; + itsHogbomGain = 0.0; + } + + itsCleaner.setFusedThreshold(itsFusedThreshold); + itsCleaner.setHogbomGain(itsHogbomGain); + + } + + itsCleaner.setLBFGSControl(itsLbfgsEpsF,itsLbfgsEpsX,itsLbfgsEpsG,itsLbfgsMaxit); + itsCleaner.setWaveletControl(itsWaveletScales, itsWaveletAmps); + + // Parts to be repeated at each minor cycle start.... + //itsCleaner.setInitScaleMasks(itsMatMask); //casa6 + Matrix tempMatMask(itsMatMask); + itsCleaner.setInitScaleMasks(tempMatMask); + itsCleaner.setaspcontrol(0, 0, 0, Quantity(0.0, "%"));/// Needs to come before the rest + + Matrix tempMat1; + tempMat1.reference(itsMatResidual); + itsCleaner.setDirty( tempMat1 ); + // InitScaleXfrs and InitScaleMasks should already be set + itsScaleSizes.clear(); + itsScaleSizes = itsCleaner.getActiveSetAspen(); + itsScaleSizes.push_back(0.0); // put 0 scale + itsCleaner.defineAspScales(itsScaleSizes); + } + + + void SDAlgorithmWaveletAspClean::takeOneStep( Float loopgain, + Int cycleNiter, + Float cycleThreshold, + Float &peakresidual, + Float &modelflux, + Int &iterdone) + { + LogIO os( LogOrigin("SDAlgorithmWaveletAspClean","takeOneStep", WHERE) ); + + Quantity thresh(cycleThreshold, "Jy"); + itsCleaner.setaspcontrol(cycleNiter, loopgain, thresh, Quantity(0.0, "%")); + Matrix tempModel; + tempModel.reference( itsMatModel ); + //save the previous model + Matrix prevModel; + prevModel = itsMatModel; + + //cout << "AAspALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl; + + // retval + // 1 = converged + // 0 = not converged but behaving normally + // -1 = not converged and stopped on cleaning consecutive smallest scale + // -2 = not converged and either large scale hit negative or diverging + // -3 = clean is diverging rather than converging + itsCleaner.startingIteration( 0 ); + Int retval = itsCleaner.aspclean( tempModel ); + iterdone = itsCleaner.numberIterations(); + + if( retval==-1 ) {os << LogIO::WARN << "AspClean minor cycle stopped on cleaning consecutive smallest scale" << LogIO::POST; } + if( retval==-2 ) {os << LogIO::WARN << "AspClean minor cycle stopped at large scale negative or diverging" << LogIO::POST;} + if( retval==-3 ) {os << LogIO::WARN << "AspClean minor cycle stopped because it is diverging" << LogIO::POST; } + + // update residual - this is critical + itsMatResidual = itsCleaner.getterResidual(); + + peakresidual = itsCleaner.getterPeakResidual(); + //cout << "SDAlg: peakres " << peakresidual << endl; + modelflux = sum( itsMatModel ); + + } + + void SDAlgorithmWaveletAspClean::finalizeDeconvolver() + { + (itsImages->residual())->put( itsMatResidual ); + (itsImages->model())->put( itsMatModel ); + } + +} //# NAMESPACE CASA - END diff --git a/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.h b/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.h new file mode 100644 index 00000000..755bc553 --- /dev/null +++ b/src/synthesis/ImagerObjects/SDAlgorithmWaveletAspClean.h @@ -0,0 +1,110 @@ +//# SDAlgorithmAAspClean.h: Definition for SDAlgorithmAAspClean +//# Copyright (C) 1996,1997,1998,1999,2000,2002 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be adressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id$ + +#ifndef SYNTHESIS_SDALGORITHMWAVELETASPCLEAN_H +#define SYNTHESIS_SDALGORITHMWAVELETASPCLEAN_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + + /* Forware Declaration */ + class SIMinorCycleController; + + + class SDAlgorithmWaveletAspClean : public SDAlgorithmBase + { + public: + + // Empty constructor + SDAlgorithmWaveletAspClean(casacore::Vector scales, casacore::Float hogbomGain, casacore::Vector waveletScales, casacore::Vector waveletAmps, casacore::Float fusedThreshold = 0.0, bool isSingle = true, casacore::Int largestScale = -1, casacore::Int stoppointmode = -1, casacore::Float lbfgsEpsF= 0.001, casacore::Float lbfgsEpsX= 0.001, casacore::Float lbfgsEpsG= 0.001, casacore::Int lbfgsMaxit= 5); + virtual ~SDAlgorithmWaveletAspClean(); + + protected: + + // Local functions to be overloaded by various algorithm deconvolvers. + virtual void takeOneStep( casacore::Float loopgain, casacore::Int cycleNiter, casacore::Float cycleThreshold, casacore::Float &peakresidual, casacore::Float &modelflux, casacore::Int &iterdone ); + virtual void initializeDeconvolver(); + virtual void finalizeDeconvolver(); + + + /* + void findNextComponent( casacore::Float loopgain ); + void updateModel(); + void updateResidual(); + */ + + /* + casacore::SubImage itsResidual, itsPsf, itsModel, itsImage; + casacore::Float itsComp; + */ + //casacore::SubImage itsResidual, itsPsf, itsModel, itsImage; + + casacore::Array itsMatPsf, itsMatResidual, itsMatModel; + casacore::Array itsMatMask; // Make an array if we eventually use multi-term masks... + + WaveletAspCleaner itsCleaner; + std::vector itsScaleSizes; + casacore::Int itsStopPointMode; + casacore::Float itsFusedThreshold; + casacore::Int itsUserLargestScale; + + casacore::Vector itsScales; + casacore::Float itsHogbomGain; + casacore::Vector itsWaveletScales; + casacore::Vector itsWaveletAmps; + casacore::Float itsLbfgsEpsF; + casacore::Float itsLbfgsEpsX; + casacore::Float itsLbfgsEpsG; + casacore::Int itsLbfgsMaxit; + + /* + casacore::IPosition itsMaxPos; + casacore::Float itsPeakResidual; + casacore::Float itsModelFlux; + */ + private: + casacore::Bool itsMCsetup; // if we should do setup or not + casacore::Float itsPrevPsfWidth; + bool itsIsSingle; + + }; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/ImagerObjects/SIImageStore.cc b/src/synthesis/ImagerObjects/SIImageStore.cc index d7227037..0e991197 100644 --- a/src/synthesis/ImagerObjects/SIImageStore.cc +++ b/src/synthesis/ImagerObjects/SIImageStore.cc @@ -3488,6 +3488,34 @@ void SIImageStore::regridToModelImage( ImageInterface &inputimage, Int te return True; } + + /////////////// + + void SIImageStore::savematrix(const Matrix& mat, const String& name){ + std::shared_ptr > imPtr; + IPosition useShape( itsParentImageShape ); + buildImage(imPtr, useShape, itsParentCoordSys, name) ; + { + LatticeLocker lock1(*imPtr, FileLocker::Write); + imPtr->set(0.0); + imPtr->flush(); + imPtr->unlock(); + } + ArrayLattice lat(itsPsf->shape()); + lat.put( mat ); + imPtr->copyData( lat ); + } + + void SIImageStore::loadmatrix(Matrix& mat, const String& name){ + std::shared_ptr > imPtr; + IPosition useShape( itsParentImageShape ); + buildImage(imPtr, name) ; + + ArrayLattice lat(itsPsf->shape()); + imPtr->copyDataTo( lat ); + Array values=lat.get(true); + mat.assign_conforming(values); + } ////////////////////////////////////////////////////////////////////////////////////////////////////// diff --git a/src/synthesis/ImagerObjects/SIImageStore.h b/src/synthesis/ImagerObjects/SIImageStore.h index 8eed9ba8..35fc76ca 100644 --- a/src/synthesis/ImagerObjects/SIImageStore.h +++ b/src/synthesis/ImagerObjects/SIImageStore.h @@ -114,6 +114,8 @@ class SIImageStore virtual std::shared_ptr > backwardGrid(casacore::uInt term=0); virtual std::shared_ptr > sumwt(casacore::uInt term=0); + virtual void savematrix(const casacore::Matrix& mat, const casacore::String &name); + virtual void loadmatrix(casacore::Matrix& mat, const casacore::String &name); virtual std::shared_ptr > alpha(){throw(casacore::AipsError("No Alpha for 1 term"));}; virtual std::shared_ptr > beta(){throw(casacore::AipsError("No Beta for 1 term"));}; diff --git a/src/synthesis/ImagerObjects/SynthesisDeconvolver.cc b/src/synthesis/ImagerObjects/SynthesisDeconvolver.cc index a5f85afc..970b433a 100644 --- a/src/synthesis/ImagerObjects/SynthesisDeconvolver.cc +++ b/src/synthesis/ImagerObjects/SynthesisDeconvolver.cc @@ -151,11 +151,26 @@ namespace casa { //# NAMESPACE CASA - BEGIN if (decpars.specmode == String("mfs")) isSingle = true; - if (decpars.nTaylorTerms > 1) // wide-band itsDeconvolver.reset(new SDAlgorithmMTAspClean(decpars.nTaylorTerms, decpars.fusedThreshold, isSingle, decpars.largestscale)); else //narrow-band. WAsp can be use for cube imaging as well but for now let AspClean handle that. - itsDeconvolver.reset(new SDAlgorithmAAspClean(decpars.fusedThreshold, isSingle, decpars.largestscale)); + itsDeconvolver.reset(new SDAlgorithmAAspClean(decpars.scales, decpars.hogbomGain, decpars.fusedThreshold, isSingle, decpars.largestscale, -1, decpars.lbfgsEpsF, decpars.lbfgsEpsX, decpars.lbfgsEpsG, decpars.lbfgsMaxit)); + } + else if (decpars.algorithm==String("waveletasp")) + { + bool isSingle = false; + if (decpars.specmode == String("mfs")) + isSingle = true; + + itsDeconvolver.reset(new SDAlgorithmWaveletAspClean(decpars.scales, decpars.hogbomGain, decpars.waveletScales, decpars.waveletAmps, decpars.fusedThreshold, isSingle, decpars.largestscale, -1, decpars.lbfgsEpsF, decpars.lbfgsEpsX, decpars.lbfgsEpsG, decpars.lbfgsMaxit)); + } + else if (decpars.algorithm==String("spectralasp")) + { + bool isSingle = false; + if (decpars.specmode == String("mfs")) + isSingle = true; + + itsDeconvolver.reset(new SDAlgorithmSpectralAspClean(decpars.scales, decpars.hogbomGain, decpars.fusedThreshold, isSingle, decpars.largestscale, -1, decpars.lbfgsEpsF, decpars.lbfgsEpsX, decpars.lbfgsEpsG, decpars.lbfgsMaxit)); } /*else if (decpars.algorithm==String("asp")) { @@ -166,6 +181,10 @@ namespace casa { //# NAMESPACE CASA - BEGIN //narrow-band itsDeconvolver.reset(new SDAlgorithmAAspClean(decpars.fusedThreshold, isSingle, decpars.largestscale)); }*/ + else if (decpars.algorithm==String("autocorrelation")) + { + itsDeconvolver.reset(new SDAlgorithmAutoClean(decpars.autoThreshold, decpars.autoMaxiter, decpars.autoGain, decpars.hogbomGain, decpars.autoHogbom, decpars.autoTrigger, decpars.autoPower, decpars.autoXMask, decpars.autoYMask)); + } else { throw( AipsError("Un-known algorithm : "+decpars.algorithm) ); diff --git a/src/synthesis/ImagerObjects/SynthesisDeconvolver.h b/src/synthesis/ImagerObjects/SynthesisDeconvolver.h index 270b98ff..e3cc00dd 100644 --- a/src/synthesis/ImagerObjects/SynthesisDeconvolver.h +++ b/src/synthesis/ImagerObjects/SynthesisDeconvolver.h @@ -43,7 +43,10 @@ #include #include #include +#include +#include #include +#include #include #include diff --git a/src/synthesis/ImagerObjects/SynthesisUtilMethods.cc b/src/synthesis/ImagerObjects/SynthesisUtilMethods.cc index d05982ee..a61b63d1 100644 --- a/src/synthesis/ImagerObjects/SynthesisUtilMethods.cc +++ b/src/synthesis/ImagerObjects/SynthesisUtilMethods.cc @@ -4007,6 +4007,76 @@ namespace casa { //# NAMESPACE CASA - BEGIN else err += "fusedthreshold must be a float or double"; } + err += readVal( inrec, String("waveletscales"), waveletScales ); + err += readVal( inrec, String("waveletamps"), waveletAmps ); + //param for the Autocleaner + if (inrec.isDefined("autothreshold")) + { + if (inrec.dataType("autothreshold") == TpFloat || inrec.dataType("autothreshold") == TpFloat) + err += readVal(inrec, String("autothreshold"), autoThreshold); + else + err += "autothreshold must be a float or double"; + } + if (inrec.isDefined("automaxiter")) + { + if (inrec.dataType("automaxiter") == TpInt) + err += readVal(inrec, String("automaxiter"), autoMaxiter); + else + err += "automaxiter must be a integer"; + } + if (inrec.isDefined("autogain")) + { + if (inrec.dataType("autogain") == TpFloat || inrec.dataType("autogain") == TpDouble) + err += readVal(inrec, String("autogain"), autoGain); + else + err += "autogain must be a float or double"; + } + if (inrec.isDefined("hogbomgain")) + { + if (inrec.dataType("hogbomgain") == TpFloat || inrec.dataType("hogbomgain") == TpDouble) + err += readVal(inrec, String("hogbomgain"), hogbomGain); + else + err += "hogbomgain must be a float or double"; + } + if (inrec.isDefined("autohogbom")) + { + if (inrec.dataType("autohogbom") == TpBool) + err += readVal(inrec, String("autohogbom"), autoHogbom); + else + err += "autohogbom must be a bool"; + } + if (inrec.isDefined("autotrigger")) + { + if (inrec.dataType("autotrigger") == TpFloat || inrec.dataType("autotrigger") == TpDouble) + err += readVal(inrec, String("autotrigger"), autoTrigger); + else + err += "autotrigger must be a float or double"; + } + if (inrec.isDefined("autopower")) + { + if (inrec.dataType("autopower") == TpFloat || inrec.dataType("autopower") == TpDouble) + err += readVal(inrec, String("autopower"), autoPower); + else + err += "autopower must be a float or double"; + } + if (inrec.isDefined("autoxmask")) + { + if (inrec.dataType("autoxmask") == TpInt) + err += readVal(inrec, String("autoxmask"), autoXMask); + else + err += "autoxmask must be a integer"; + } + if (inrec.isDefined("autoxmask")) + { + if (inrec.dataType("autoymask") == TpInt) + err += readVal(inrec, String("autoymask"), autoYMask); + else + err += "autoymask must be a integer"; + } + err += readVal( inrec, String("lbfgsepsf"), lbfgsEpsF ); + err += readVal( inrec, String("lbfgsepsx"), lbfgsEpsX ); + err += readVal( inrec, String("lbfgsepsg"), lbfgsEpsG ); + err += readVal( inrec, String("lbfgsMaxit"), lbfgsMaxit ); if (inrec.isDefined("specmode")) { if(inrec.dataType("specmode") == TpString) @@ -4307,6 +4377,21 @@ namespace casa { //# NAMESPACE CASA - BEGIN interactive=false; autoAdjust=False; fusedThreshold = 0.0; + waveletScales.resize(1); waveletScales[0] = 0.0 ; + waveletAmps.resize(1); waveletAmps[0] = 0.0 ; + autoThreshold = 0.0; + autoMaxiter = 0; + autoGain = 0.0; + hogbomGain = 0.0; + autoHogbom = false; + autoTrigger = 1.0; + autoPower = 1.0; + autoXMask = 0; + autoYMask = 0; + lbfgsEpsF = 0.001; + lbfgsEpsX = 0.001; + lbfgsEpsG = 0.001; + lbfgsMaxit = 5; specmode="mfs"; largestscale = -1; } @@ -4324,6 +4409,21 @@ namespace casa { //# NAMESPACE CASA - BEGIN decpar.define("scalebias",scalebias); decpar.define("usemask",maskType); decpar.define("fusedthreshold", fusedThreshold); + decpar.define("waveletscales", waveletScales); + decpar.define("waveletamps", waveletAmps); + decpar.define("autothreshold", autoThreshold); + decpar.define("automaxiter", autoMaxiter); + decpar.define("autogain", autoGain); + decpar.define("hogbomgain", hogbomGain); + decpar.define("autohogbom", autoHogbom); + decpar.define("autotrigger", autoTrigger); + decpar.define("autopower", autoPower); + decpar.define("autoxmask", autoXMask); + decpar.define("autoymask", autoYMask); + decpar.define("lbfgsepsf", lbfgsEpsF); + decpar.define("lbfgsepsx", lbfgsEpsX); + decpar.define("lbfgsepsg", lbfgsEpsG); + decpar.define("lbfgsmaxit", lbfgsMaxit); decpar.define("specmode", specmode); decpar.define("largestscale", largestscale); if( maskList.nelements()==1 && maskList[0]=="") diff --git a/src/synthesis/ImagerObjects/SynthesisUtilMethods.h b/src/synthesis/ImagerObjects/SynthesisUtilMethods.h index fa160881..a2e6611c 100644 --- a/src/synthesis/ImagerObjects/SynthesisUtilMethods.h +++ b/src/synthesis/ImagerObjects/SynthesisUtilMethods.h @@ -460,6 +460,21 @@ class SynthesisParams int nMask; bool autoAdjust; casacore::Float fusedThreshold; + casacore::Vector waveletScales; + casacore::Vector waveletAmps; + casacore::Float autoThreshold; + casacore::Int autoMaxiter; + casacore::Float autoGain; + casacore::Float hogbomGain; + casacore::Bool autoHogbom; + casacore::Float autoTrigger; + casacore::Float autoPower; + casacore::Int autoXMask; + casacore::Int autoYMask; + casacore::Float lbfgsEpsF; + casacore::Float lbfgsEpsX; + casacore::Float lbfgsEpsG; + casacore::Int lbfgsMaxit; casacore::String specmode; casacore::Int largestscale; // task deconvolve needs to tell siimagestore that we don't need to check for the sumwt image diff --git a/src/synthesis/ImagerObjects/test/hummbee_func.h b/src/synthesis/ImagerObjects/test/hummbee_func.h index 3b2a77e9..f9b1492d 100644 --- a/src/synthesis/ImagerObjects/test/hummbee_func.h +++ b/src/synthesis/ImagerObjects/test/hummbee_func.h @@ -46,6 +46,8 @@ #include #include #include +#include +#include #include #include @@ -269,7 +271,7 @@ void Hummbee(bool restartUI, int argc, char **argv, string& MSNBuf, Float& pbLimit, string& deconvolver, vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain, float& hogbomgain,bool& autohogbom, float& autotrigger, float& autopower, int& autoxmask, int& autoymask, float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, @@ -322,6 +324,21 @@ void Hummbee(bool restartUI, int argc, char **argv, string& MSNBuf, decPars_p.interactive=false; decPars_p.autoAdjust=False; //genie decPars_p.fusedThreshold = fusedthreshold; + decPars_p.waveletScales = Vector(waveletscales); + decPars_p.waveletAmps = Vector(waveletamps); + decPars_p.autoThreshold = autothreshold; + decPars_p.autoMaxiter = automaxiter; + decPars_p.autoGain = autogain; + decPars_p.hogbomGain = hogbomgain; + decPars_p.autoHogbom = autohogbom; + decPars_p.autoTrigger = autotrigger; + decPars_p.autoPower = autopower; + decPars_p.autoXMask = autoxmask; + decPars_p.autoYMask = autoymask; + decPars_p.lbfgsEpsF = lbfgsepsf; + decPars_p.lbfgsEpsX = lbfgsepsx; + decPars_p.lbfgsEpsG = lbfgsepsg; + decPars_p.lbfgsMaxit = lbfgsmaxit; decPars_p.specmode="mfs"; //deconvolve task does not have this decPars_p.largestscale = largestscale; decPars_p.scalebias = 0.0; diff --git a/src/synthesis/ImagerObjects/test/tHummbee.cc b/src/synthesis/ImagerObjects/test/tHummbee.cc index 4f9b6691..d3f3ebda 100644 --- a/src/synthesis/ImagerObjects/test/tHummbee.cc +++ b/src/synthesis/ImagerObjects/test/tHummbee.cc @@ -62,7 +62,7 @@ void UI(bool restart, int argc, char **argv, string& MSNBuf, Float& pbLimit, string& deconvolver, vector& scales, - float& largestscale, float& fusedthreshold, + float& largestscale, float& fusedthreshold,vector& waveletscales,vector& waveletamps,float& autothreshold,int& automaxiter,float& autogain, float& hogbomgain,bool& autohogbom,float& autotrigger,float& autopower, int& autoxmask, int& autoymask, float& lbfgsepsf, float& lbfgsepsx, float& lbfgsepsg, int& lbfgsmaxit, int& nterms, float& gain, float& threshold, float& nsigma, @@ -126,6 +126,21 @@ void UI(bool restart, int argc, char **argv, string& MSNBuf, //InitMap(watchPoints,exposedKeys); exposedKeys.push_back("largestscale"); exposedKeys.push_back("fusedthreshold"); + exposedKeys.push_back("waveletscales"); + exposedKeys.push_back("waveletamps"); + exposedKeys.push_back("autothreshold"); + exposedKeys.push_back("automaxiter"); + exposedKeys.push_back("autogain"); + exposedKeys.push_back("hogbomgain"); + exposedKeys.push_back("autohogbom"); + exposedKeys.push_back("autotrigger"); + exposedKeys.push_back("autopower"); + exposedKeys.push_back("autoxmask"); + exposedKeys.push_back("autoymask"); + exposedKeys.push_back("lbfgsepsf"); + exposedKeys.push_back("lbfgsepsx"); + exposedKeys.push_back("lbfgsepsg"); + exposedKeys.push_back("lbfgsmaxit"); watchPoints["asp"]=exposedKeys; i=1;clgetSValp("deconvolver", deconvolver, i ,watchPoints); @@ -136,6 +151,21 @@ void UI(bool restart, int argc, char **argv, string& MSNBuf, i=1;clgetFValp("largestscale", largestscale,i); i=1;clgetFValp("fusedthreshold", fusedthreshold,i); + N=0; N=clgetNFValp("waveletscales", waveletscales, N); + N=0; N=clgetNFValp("waveletamps", waveletamps, N); + i=1;clgetFValp("autothreshold", autothreshold,i); + i=1;clgetIValp("automaxiter", automaxiter,i); + i=1;clgetFValp("autogain", autogain,i); + i=1;clgetFValp("hogbomgain", hogbomgain,i); + i=1;clgetBValp("autohogbom", autohogbom,i); + i=1;clgetFValp("autotrigger", autotrigger, i); + i=1;clgetFValp("autopower", autopower,i); + i=1;clgetIValp("autoxmask", autoxmask,i); + i=1;clgetIValp("autoymask", autoymask,i); + i=1;clgetFValp("lbfgsepsf", lbfgsepsf,i); + i=1;clgetFValp("lbfgsepsx", lbfgsepsx,i); + i=1;clgetFValp("lbfgsepsg", lbfgsepsg,i); + i=1;clgetIValp("lbfgsmaxit", lbfgsmaxit,i); i=1;clgetIValp("nterms", nterms,i); i=1;clgetFValp("gain", gain,i); @@ -218,6 +248,21 @@ int main(int argc, char **argv) vector scales; float largestscale = -1; float fusedthreshold = 0; + vector waveletscales; + vector waveletamps; + float autothreshold = 0; + int automaxiter = 0; + float autogain = 0; + float hogbomgain = 0; + bool autohogbom = 0; + float autotrigger= 0; + float autopower = 1; + int autoxmask = 0; + int autoymask = 0; + float lbfgsepsf = 0.001; + float lbfgsepsx = 0.001; + float lbfgsepsg = 0.001; + int lbfgsmaxit = 5; int nterms=2; float gain=0.1; float threshold=0.0; @@ -233,7 +278,7 @@ int main(int argc, char **argv) ,doPBCorr, conjBeams, pbLimit, deconvolver, scales, - largestscale, fusedthreshold, + largestscale, fusedthreshold,waveletscales, waveletamps,autothreshold,automaxiter,autogain,hogbomgain,autohogbom,autotrigger,autopower,autoxmask,autoymask,lbfgsepsf,lbfgsepsx,lbfgsepsg,lbfgsmaxit, nterms, gain, threshold, nsigma, @@ -251,7 +296,7 @@ int main(int argc, char **argv) doPBCorr, conjBeams, pbLimit, deconvolver, scales, - largestscale, fusedthreshold, + largestscale, fusedthreshold,waveletscales,waveletamps,autothreshold,automaxiter,autogain,hogbomgain,autohogbom,autotrigger,autopower,autoxmask,autoymask,lbfgsepsf,lbfgsepsx,lbfgsepsg,lbfgsmaxit, nterms, gain, threshold, nsigma, diff --git a/src/synthesis/MeasurementEquations/AspMatrixCleaner.cc b/src/synthesis/MeasurementEquations/AspMatrixCleaner.cc index 9eaa5b23..528127eb 100644 --- a/src/synthesis/MeasurementEquations/AspMatrixCleaner.cc +++ b/src/synthesis/MeasurementEquations/AspMatrixCleaner.cc @@ -29,7 +29,7 @@ #include #include #include -//#include +#include #include #include #include @@ -97,6 +97,11 @@ AspMatrixCleaner::AspMatrixCleaner(): itsPrevPeakResidual(0.0), itsOrigDirty( ), itsFusedThreshold(0.0), + itsdimensionsareeven(true), + itsLbfgsEpsF(0.001), + itsLbfgsEpsX(0.001), + itsLbfgsEpsG(0.001), + itsLbfgsMaxit(5), itsNumNoChange(0), itsBinSizeForSumFlux(4), itsUserLargestScale(-1.0) @@ -146,7 +151,9 @@ Int AspMatrixCleaner::aspclean(Matrix& model, LogIO os(LogOrigin("AspMatrixCleaner", "aspclean()", WHERE)); os << LogIO::NORMAL1 << "Asp clean algorithm" << LogIO::POST; - + + if(itsHogbomGain == 0) + itsHogbomGain = itsGain; //Int scale; @@ -255,6 +262,9 @@ Int AspMatrixCleaner::aspclean(Matrix& model, Float initRMSResidual = 1000.0; // float initModelFlux = 0.0; + itsdimensionsareeven = (psfShape_p(0) == 2*(psfShape_p(0)/2)); + Float tempGain; + os < tempScaleSizes; itsIteration = itsStartingIter; // 0 @@ -320,7 +330,7 @@ Int AspMatrixCleaner::aspclean(Matrix& model, rms = rms / num; // make single optimized scale image - os << LogIO::NORMAL3 << "Making optimized scale " << itsOptimumScaleSize << " at location " << itsPositionOptimum << LogIO::POST; + os << "Making optimized scale " << itsOptimumScaleSize << " at location " << itsPositionOptimum << LogIO::POST; if (itsSwitchedToHogbom) { @@ -422,8 +432,14 @@ Int AspMatrixCleaner::aspclean(Matrix& model, itsNumIterNoGoodAspen.push_back(0); } + //which loop gain to use + if (itsOptimumScaleSize == 0.0) + tempGain = itsHogbomGain; + else + tempGain = itsGain; + // Now add to the total flux - totalFlux += (itsStrengthOptimum*itsGain); + totalFlux += (itsStrengthOptimum*tempGain); itsTotalFlux = totalFlux; if(ii == itsStartingIter) @@ -436,7 +452,7 @@ Int AspMatrixCleaner::aspclean(Matrix& model, os < abs(itsPeakResidual)) minMaximumResidual = abs(itsPeakResidual); @@ -550,7 +566,7 @@ Int AspMatrixCleaner::aspclean(Matrix& model, // Update the model image Matrix modelSub = model(blc, trc); Float scaleFactor; - scaleFactor = itsGain * itsStrengthOptimum; + scaleFactor = tempGain * itsStrengthOptimum; Matrix scaleSub = (itsScale)(blcPsf,trcPsf); modelSub += scaleFactor * scaleSub; @@ -561,9 +577,22 @@ Int AspMatrixCleaner::aspclean(Matrix& model, Matrix itsPsfConvScale = Matrix(psfShape_p); fft.fft0(itsPsfConvScale, cWork, false); fft.flip(itsPsfConvScale, false, false); //need this if conv with 1 scale; don't need this if conv with 2 scales + + /*IPosition nullnull(2,0); + Matrix shift(psfShape_p); + shift.assign_conforming(itsPsfConvScale); + if (itsdimensionsareeven){ + Matrix sub = itsPsfConvScale(nullnull+1,support-1); + sub.assign_conforming(shift(nullnull,support-2)); + } + else{ + Matrix sub = itsPsfConvScale(nullnull+2,support-1); + sub.assign_conforming(shift(nullnull,support-3)); + }*/ + Matrix psfSub = (itsPsfConvScale)(blcPsf, trcPsf); Matrix dirtySub=(*itsDirty)(blc,trc); - + /* debug info float maxvalue; IPosition peakpos; @@ -579,7 +608,7 @@ Int AspMatrixCleaner::aspclean(Matrix& model, cout << "itsDirty pos " << peakpos << " maxval " << (*itsDirty)(peakpos) << endl; cout << "itsPositionOptimum " << itsPositionOptimum << endl; cout << " maxPsfSub " << max(fabs(psfSub)) << " maxPsfConvScale " << max(fabs(itsPsfConvScale)) << " itsGain " << itsGain << endl;*/ - os <& model, switchedToHogbom(); }* / }*/ - + // update peakres itsPrevPeakResidual = itsPeakResidual; maxVal=0; @@ -674,7 +703,7 @@ Int AspMatrixCleaner::aspclean(Matrix& model, (fabs(itsPeakResidual - itsPrevPeakResidual) < 1e-4)) //peakres rarely changes itsNumNoChange += 1; //cout << "after: itsDirty optPos " << (*itsDirty)(itsPositionOptimum) << endl; - + // If we switch to hogbom (i.e. only have 0 scale size), // we still need to do the following Aspen update to get the new optimumStrength if (itsSwitchedToHogbom) @@ -716,13 +745,28 @@ Int AspMatrixCleaner::aspclean(Matrix& model, defineAspScales(tempScaleSizes); } // End of iteration - + vector sumFluxByBins(itsBinSizeForSumFlux,0.0); vector rangeFluxByBins(itsBinSizeForSumFlux+1,0.0); getFluxByBins(vecItsOptimumScaleSize,vecItsStrengthOptimum,itsBinSizeForSumFlux,sumFluxByBins,rangeFluxByBins); + //we will need this later for the MF option, comment for now + /*write_array(Optimums, std::string("./strengthoptimum")); + write_array(ScaleSizes, std::string("./scalesizes")); + casacore::Vector xPositions(itsMaxNiter); + casacore::Vector yPositions(itsMaxNiter); + + for (int ii=0; ii < itsMaxNiter; ii++) + { + xPositions(ii) = Positions(ii)(0); + yPositions(ii) = Positions(ii)(1); + } + + //write_array(Positions, std::string("./positions")); + write_array(xPositions, std::string("./xpositions")); + write_array(yPositions, std::string("./ypositions"));*/ os << " The number of bins for collecting the sum of Flux: " << itsBinSizeForSumFlux << endl; @@ -781,6 +825,20 @@ Bool AspMatrixCleaner::destroyInitMasks() return true; } +float AspMatrixCleaner::aspScaleModel(const Int& i, const Int& j, const Float& scaleSize, const Int& refi, const Int& refj) +{ + return (1.0/(sqrt(2*M_PI)*scaleSize))*exp(-(pow(i-refi,2) + pow(j-refj,2))*0.5/pow(scaleSize,2)); +} + +float AspMatrixCleaner::aspPeakNormModel(const Float& width) +{ + return sqrt(2 * M_PI *width); +} + +float AspMatrixCleaner::aspNormalizationModel(const Float& width1, const Float& width2) +{ + return sqrt(2 * M_PI / (pow(1.0/width1, 2) + pow(1.0/width2, 2))); +} float AspMatrixCleaner::getPsfGaussianWidth(ImageInterface& psf) { @@ -854,8 +912,7 @@ void AspMatrixCleaner::makeInitScaleImage(Matrix& iscale, const Float& sc //const int px = i - refi; //const int py = j - refj; //iscale(i,j) = gbeam(px, py); // gbeam with the above def is equivalent to the following - iscale(i,j) = (1.0/(sqrt(2*M_PI)*scaleSize))*exp(-(pow(i-refi,2) + pow(j-refj,2))*0.5/pow(scaleSize,2)); //this is for 1D, but represents Sanjay's and gives good init scale - //iscale(i,j) = (1.0/(2*M_PI*pow(scaleSize,2)))*exp(-(pow(i-refi,2) + pow(j-refj,2))*0.5/pow(scaleSize,2)); // this is for 2D, gives unit area but bad init scale (always picks 0) + iscale(i,j) = aspScaleModel(i, j, scaleSize, refi, refj); } } @@ -902,10 +959,7 @@ void AspMatrixCleaner::makeScaleImage(Matrix& iscale, const Float& scaleS // iscale(i,j) = gbeam(px, py); // this is equivalent to the following with the above gbeam definition // This is for 1D, but represents Sanjay's and gives good init scale // Note that "amp" is not used in the expression - iscale(i,j) = (1.0/(sqrt(2*M_PI)*scaleSize))*exp(-(pow(i-center[0],2) + pow(j-center[1],2))*0.5/pow(scaleSize,2)); - - // This is for 2D, gives unit area but bad init scale (always picks 0) - // iscale(i,j) = (1.0/(2*M_PI*pow(scaleSize,2)))*exp(-(pow(i-center[0],2) + pow(j-center[1],2))*0.5/pow(scaleSize,2)); + iscale(i,j) = aspScaleModel(i, j, scaleSize, center[0], center[1]); } } @@ -1140,6 +1194,42 @@ void AspMatrixCleaner::setInitScaleXfrs(const Float width) } } +void AspMatrixCleaner::loadInitScaleXfrs(const casacore::Vector & scales) +{ + if(itsInitScales.nelements() > 0) + destroyAspScales(); + + if (itsSwitchedToHogbom) + { + itsNInitScales = 1; + itsInitScaleSizes.resize(itsNInitScales, false); + itsInitScaleSizes = {0.0f}; + } + else + { + itsNInitScales = scales.size(); + itsInitScaleSizes.resize(itsNInitScales, false); + Int scale = 0; + while (scale < itsNInitScales) + { + itsInitScaleSizes[scale] = scales(scale); + scale++; + } + } + + itsInitScales.resize(itsNInitScales, false); + itsInitScaleXfrs.resize(itsNInitScales, false); + fft = FFTServer(psfShape_p); + for (int scale = 0; scale < itsNInitScales; scale++) + { + itsInitScales[scale] = Matrix(psfShape_p); + makeInitScaleImage(itsInitScales[scale], itsInitScaleSizes[scale]); + //cout << "made itsInitScales[" << scale << "] = " << itsInitScaleSizes[scale] << endl; + itsInitScaleXfrs[scale] = Matrix (); + fft.fft0(itsInitScaleXfrs[scale], itsInitScales[scale]); + } +} + // calculate the convolutions of the psf with the initial scales void AspMatrixCleaner::setInitScalePsfs() { @@ -1435,9 +1525,7 @@ void AspMatrixCleaner::maxDirtyConvInitScales(float& strengthOptimum, int& optim if (scale > 0) { float normalization; - //normalization = 2 * M_PI / pow(itsInitScaleSizes[scale], 2); //sanjay's - //normalization = 2 * M_PI / pow(itsInitScaleSizes[scale], 1); // this looks good on M31 but bad on G55 - normalization = sqrt(2 * M_PI *itsInitScaleSizes[scale]); //GSL. Need to recover and re-norm later + normalization = aspPeakNormModel(itsInitScaleSizes[scale]); //GSL. Need to recover and re-norm later maxima(scale) /= normalization; } //sanjay's code doesn't normalize peak here. // Norm Method 2 may work fine with GSL with derivatives, but Norm Method 1 is still the preferred approach. @@ -1458,8 +1546,8 @@ void AspMatrixCleaner::maxDirtyConvInitScales(float& strengthOptimum, int& optim if (optimumScale > 0) { - //const float normalization = 2 * M_PI / (pow(1.0/itsPsfWidth, 2) + pow(1.0/itsInitScaleSizes[optimumScale], 2)); // sanjay - const float normalization = sqrt(2 * M_PI / (pow(1.0/itsPsfWidth, 2) + pow(1.0/itsInitScaleSizes[optimumScale], 2))); // this is good. + + float normalization = aspNormalizationModel(itsPsfWidth, itsInitScaleSizes[optimumScale]); // this is good. // norm method 2 recovers the optimal strength and then normalize it to get the init guess if (itsNormMethod == 2) @@ -1500,6 +1588,7 @@ vector AspMatrixCleaner::getActiveSetAspen() fft.fft0(dirtyFT, *itsDirty); itsDirtyConvInitScales.resize(0); itsDirtyConvInitScales.resize(itsNInitScales); // 0, 1width, 2width, 4width and 8width + //cout << "itsInitScaleSizes.size() " << itsInitScaleSizes.size() << " itsInitScales.size() " << itsInitScales.size() << " NInitScales # " << itsNInitScales << endl; for (int scale=0; scale < itsNInitScales; scale++) { @@ -1523,7 +1612,7 @@ vector AspMatrixCleaner::getActiveSetAspen() maxDirtyConvInitScales(strengthOptimum, optimumScale, positionOptimum); - os << LogIO::NORMAL3 << "Peak among the smoothed residual image is " << strengthOptimum << " and initial scale: " << optimumScale << LogIO::POST; + os << "Peak among the smoothed residual image is " << strengthOptimum << " and initial scale: " << optimumScale << LogIO::POST; //cout << "Peak among the smoothed residual image is " << strengthOptimum << " and initial scale: " << optimumScale << endl; // cout << " its itsDirty is " << (*itsDirty)(positionOptimum); // cout << " at location " << positionOptimum[0] << " " << positionOptimum[1] << " " << positionOptimum[2]; @@ -1566,28 +1655,42 @@ vector AspMatrixCleaner::getActiveSetAspen() s[i] = tempx[i]; //amp s[i+1] = tempx[i+1]; //scale } + + + //real_1d_array s = "[1,1]"; + double epsg = itsLbfgsEpsG;//1e-3; + double epsf = itsLbfgsEpsF;//1e-3; + double epsx = itsLbfgsEpsX;//1e-3; + ae_int_t maxits = itsLbfgsMaxit;//5; + //double epsg = 1e-3; + //double epsf = 1e-3; + //double epsx = 1e-3; + //ae_int_t maxits = 5; + minlbfgsstate state; + minlbfgscreate(1, x, state); + minlbfgssetcond(state, epsg, epsf, epsx, maxits); + minlbfgssetscale(state, s); + + //minlbfgssetprecscale(state); + minlbfgsreport rep; + + + ParamAlglibObj optParam(*itsDirty, *itsXfr, activeSetCenter, fft); + ParamAlglibObj *ptrParam; + ptrParam = &optParam; + + alglib::minlbfgsoptimize(state, objfunc_alglib, NULL, (void *) ptrParam); + + minlbfgsresults(state, x, rep); + //double *x1 = x.getcontent(); + //cout << "x1[0] " << x1[0] << " x1[1] " << x1[1] << endl; + + // end alglib bfgs optimization + + + - ParamAlglibObj optParam(*itsDirty, *itsXfr, activeSetCenter, fft); - ParamAlglibObj *ptrParam; - ptrParam = &optParam; - //real_1d_array s = "[1,1]"; - //real_1d_array s = "[0.001,10]"; - double epsg = 1e-3; - double epsf = 1e-3; - double epsx = 1e-3; - ae_int_t maxits = 5; - minlbfgsstate state; - minlbfgscreate(1, x, state); - minlbfgssetcond(state, epsg, epsf, epsx, maxits); - minlbfgssetscale(state, s); - minlbfgsreport rep; - alglib::minlbfgsoptimize(state, objfunc_alglib, NULL, (void *) ptrParam); - minlbfgsresults(state, x, rep); - //double *x1 = x.getcontent(); - //cout << "x1[0] " << x1[0] << " x1[1] " << x1[1] << endl; - - // end alglib bfgs optimization double amp = x[0]; // i double scale = x[1]; // i+1 @@ -1609,7 +1712,7 @@ vector AspMatrixCleaner::getActiveSetAspen() itsGoodAspCenter = activeSetCenter; // debug - os << LogIO::NORMAL3 << "optimized strengthOptimum " << itsStrengthOptimum << " scale size " << itsOptimumScaleSize << LogIO::POST; + //os << "optimized strengthOptimum " << itsStrengthOptimum << " scale size " << itsOptimumScaleSize << LogIO::POST; //cout << "optimized strengthOptimum " << itsStrengthOptimum << " scale size " << itsOptimumScaleSize << " at " << itsPositionOptimum << endl; } // finish bfgs optimization diff --git a/src/synthesis/MeasurementEquations/AspMatrixCleaner.h b/src/synthesis/MeasurementEquations/AspMatrixCleaner.h index 4c4a43cb..ba7a3015 100644 --- a/src/synthesis/MeasurementEquations/AspMatrixCleaner.h +++ b/src/synthesis/MeasurementEquations/AspMatrixCleaner.h @@ -36,6 +36,7 @@ #include #include #include +#include namespace casa { //# NAMESPACE CASA - BEGIN @@ -65,6 +66,11 @@ class AspMatrixCleaner : public MatrixCleaner // helper functions for ASP float getPsfGaussianWidth(casacore::ImageInterface& psf); void getLargestScaleSize(casacore::ImageInterface& psf); + + //Normalization Helpers + virtual float aspScaleModel(const casacore::Int& i, const casacore::Int& j, const casacore::Float& scaleSize, const casacore::Int& refi, const casacore::Int& refj); + virtual float aspPeakNormModel(const casacore::Float& width); + virtual float aspNormalizationModel(const casacore::Float& width1, const casacore::Float& width2); // Make an image of the specified scale by Gaussian void makeInitScaleImage(casacore::Matrix& iscale, const casacore::Float& scaleSize); @@ -72,6 +78,7 @@ class AspMatrixCleaner : public MatrixCleaner void setInitScales(); void setInitScaleXfrs(const casacore::Float width); + void loadInitScaleXfrs(const casacore::Vector & scales); // calculate the convolutions of the psf with the initial scales void setInitScalePsfs(); @@ -82,7 +89,8 @@ class AspMatrixCleaner : public MatrixCleaner void maxDirtyConvInitScales(float& strengthOptimum, int& optimumScale, casacore::IPosition& positionOptimum); // returns the active-set aspen for cleaning - std::vector getActiveSetAspen(); + virtual std::vector getActiveSetAspen(); + //void MFaspclean(casacore::Matrix & model); // Juat define the active-set aspen scales void defineAspScales(std::vector& scaleSizes); @@ -90,6 +98,7 @@ class AspMatrixCleaner : public MatrixCleaner void switchedToHogbom(bool runlong= false); void setOrigDirty(const casacore::Matrix& dirty); void setFusedThreshold(const casacore::Float fusedThreshold = 0.0) { itsFusedThreshold = fusedThreshold; } + void setHogbomGain(const casacore::Float hogbomGain = 0.0) { itsHogbomGain = hogbomGain; } void setUserLargestScale(const casacore::Int largestScale = -1) { itsUserLargestScale = float(largestScale); } // setter/getter @@ -102,6 +111,7 @@ class AspMatrixCleaner : public MatrixCleaner void setBinSizeForSumFlux(const casacore::Int binSize = 4) { itsBinSizeForSumFlux = binSize ; } ; void getFluxByBins(const std::vector& scaleSizes,const std::vector& optimum, casacore::Int binSize, std::vector& sumFluxByBins, std::vector& rangeFluxByBins); + void setLBFGSControl(const casacore::Float LbfgsEpsF, const casacore::Float LbfgsEpsX, const casacore::Float LbfgsEpsG, const casacore::Int LbfgsMaxit) { itsLbfgsEpsF = LbfgsEpsF; itsLbfgsEpsX = LbfgsEpsX; itsLbfgsEpsG = LbfgsEpsG; itsLbfgsMaxit = LbfgsMaxit;} protected: //private: @@ -177,11 +187,18 @@ class AspMatrixCleaner : public MatrixCleaner const casacore::Int itsDefaultNorm = 1; casacore::Int itsNormMethod; casacore::Float itsFusedThreshold; + casacore::Float itsHogbomGain; unsigned int itsNumNoChange; // number of times peakres rarely changes casacore::Int itsBinSizeForSumFlux ; // number of bins for histogram of the sum of Flux float itsLargestInitScale; // estimated largest initial scale float itsUserLargestScale; // user-specified largest initial scale casacore::IPosition blcDirty, trcDirty; + + casacore::Bool itsdimensionsareeven; + casacore::Float itsLbfgsEpsF; + casacore::Float itsLbfgsEpsX; + casacore::Float itsLbfgsEpsG; + casacore::Int itsLbfgsMaxit; }; } //# NAMESPACE CASA - END diff --git a/src/synthesis/MeasurementEquations/AutoCleaner.cc b/src/synthesis/MeasurementEquations/AutoCleaner.cc new file mode 100644 index 00000000..1798dc11 --- /dev/null +++ b/src/synthesis/MeasurementEquations/AutoCleaner.cc @@ -0,0 +1,1624 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#ifdef _OPENMP +#include +#endif +using namespace casacore; +namespace casa { //# NAMESPACE CASA - BEGIN + +Bool AutoCleaner::validatePsf(const Matrix & psf) +{ + LogIO os(LogOrigin("AutoCleaner", "validatePsf()", WHERE)); + + //Find the peak of the raw psf + AlwaysAssert(psf.shape().product() != 0, AipsError); + Float maxPsf=0; + itsPositionPeakPsf=IPosition(psf.shape().nelements(), 0); + Int psfSupport = findBeamPatch(0.0,psf.shape()(0), psf.shape()(1), 4.0, 20.0); + findPSFMaxAbs(psf, maxPsf, itsPositionPeakPsf, psfSupport); + os << "Peak of PSF = " << maxPsf << " at " << itsPositionPeakPsf << LogIO::POST; + return true; +} + +AutoCleaner::AutoCleaner(): + itsMask( ), + itsMaskThreshold(0.9), + itsDirty( ), + itsMaximumResidual(0.0), + itsStrengthOptimum(0.0), + itsTotalFlux(0.0), + itsChoose(true), + itsDoSpeedup(false), + itsIgnoreCenterBox(false), + psfShape_p(0), + noClean_p(false) +{ + itsMemoryMB=Double(HostInfo::memoryTotal()/1024)/16.0; + itsStartingIter = 0; + itsFracThreshold=Quantity(0.0, "%"); +} + +AutoCleaner::AutoCleaner(const Matrix & psf, + const Matrix &dirty): + itsMask( ), + itsMaximumResidual(0.0), + itsStrengthOptimum(0.), + itsTotalFlux(0.0), + itsChoose(true), + itsDoSpeedup(false), + itsIgnoreCenterBox(false), + itsdimensionsareeven(true), + noClean_p(false) +{ + AlwaysAssert(validatePsf(psf), AipsError); + psfShape_p.resize(0, false); + psfShape_p=psf.shape(); + // Check that everything is the same dimension and that none of the + // dimensions is zero length. + AlwaysAssert(psf.shape().nelements() == dirty.shape().nelements(), + AipsError); + AlwaysAssert(dirty.shape().product() != 0, AipsError); + // looks OK so make the convolver + + // We need to guess the memory use. For the moment, we'll assume + // that about 4 scales will be used, giving about 32 TempLattices + // in all. Also we'll try not to take more that half of the memory + + // Ah, but when we are doing a mosaic, its actually worse than this! + // So, we pass it in + itsMemoryMB=Double(HostInfo::memoryTotal()/1024)/16.0; + + Matrix itsDirty; + itsDirty.resize(psf.shape()); + itsDirty.assign(dirty); + setPsf(psf); + itsStartingIter = 0; + itsFracThreshold=Quantity(0.0, "%"); +} + +void AutoCleaner::setPsf(const Matrix& psf){ + LogIO os(LogOrigin("AutoCleaner", "setPsf()", WHERE)); + AlwaysAssert(validatePsf(psf), AipsError); + psfShape_p.resize(0, false); + psfShape_p = psf.shape(); + itsPsf.resize(psf.shape()); + itsPsf.assign(psf); + + cmap.resize(psf.shape()); + mod.resize(psf.shape()); + delta.resize(psf.shape()); + + shifted.resize(psf.shape()); + clean_map.resize(psf.shape()); + + BB.resize(psf.shape()); + tildeII.resize(psf.shape()); + BI.resize(psf.shape()); + basis_function.resize(psf.shape()); + tildeMB.resize(psf.shape()); + mspsf.resize(psf.shape()); + tildeMI.resize(psf.shape()); + tildeMBB.resize(psf.shape()); + MtildeMBB.resize(psf.shape()); + tildeMBI.resize(psf.shape()); + MtildeMB.resize(psf.shape()); + window_basis.resize(psf.shape()); + + cmap=0.0; + mod=0.0; + delta=0.0; + + shifted=0.0; + clean_map=0.0; + + BB=0.0; + tildeII=0.0; + BI=0.0; + basis_function=0.0; + tildeMB=0.0; + mspsf=0.0; + tildeMI=0.0; + tildeMBB=0.0; + MtildeMBB=0.0; + tildeMBI=0.0; + MtildeMB=0.0; + window_basis=0.0; + + normcmap=0.0; + normcmapsq=0.0; + sumcmap=0.0; + blcautomask=itsDirty.shape(); + trcautomask=itsDirty.shape(); + itsAutoThreshold=0.0; + itsAutoMaxiter=0; + itsAutoGain=0.0; + itsHogbomGain=0.0; + itsAutoHogbom=false; + itsAutoTrigger=1.0; + triggerhogbom=0.0; + itsAutoPower=1.0; + itsmaxbeam=1.0; +} + +AutoCleaner::AutoCleaner(const AutoCleaner & other) + +{ + operator=(other); +} + +AutoCleaner & AutoCleaner::operator=(const AutoCleaner & other) { + if (this != &other) { + itsCleanType = other.itsCleanType; + itsMask = other.itsMask; + itsDirty = other.itsDirty; + itsScales = other.itsScales; + itsStartingIter = other.itsStartingIter; + itsMaximumResidual = other.itsMaximumResidual; + itsIgnoreCenterBox = other.itsIgnoreCenterBox; + itsStrengthOptimum = other.itsStrengthOptimum; + itsTotalFlux=other.itsTotalFlux; + itsMaskThreshold = other.itsMaskThreshold; + psfShape_p.resize(0, false); + psfShape_p=other.psfShape_p; + noClean_p=other.noClean_p; + } + return *this; +} + +AutoCleaner::~AutoCleaner() +{ + if(!itsMask.null()) itsMask=0; +} + +void AutoCleaner::update(const Matrix &dirty) +{ + //itsDirty->assign(dirty); + itsDirty.assign(dirty); + + LogIO os(LogOrigin("AutoCleaner", "update()", WHERE)); + +// here we need to reupdate the other matrices + +} + +// add a mask image +void AutoCleaner::setMask(Matrix & mask, const Float& maskThreshold) +{ + itsMaskThreshold = maskThreshold; + IPosition maskShape = mask.shape(); + + //cerr << "Mask Shape " << mask.shape() << endl; + // This is not needed after the first steps + itsMask = new Matrix(mask.shape()); + itsMask->assign(mask); + noClean_p=(max(*itsMask) < itsMaskThreshold) ? true : false; + +} + +Bool AutoCleaner::setcontrol(CleanEnums::CleanType cleanType, + const Int niter, + const Float gain, + const Quantity& threshold) +{ + return setcontrol(cleanType, niter, gain, threshold, Quantity(0.0, "%")); +} + +// Set up the control parameters +Bool AutoCleaner::setcontrol(CleanEnums::CleanType cleanType, + const Int niter, + const Float gain, + const Quantity& aThreshold, + const Quantity& fThreshold) +{ + itsCleanType=cleanType; + itsMaxNiter=niter; + itsGain=gain; + itsThreshold=aThreshold; + itsFracThreshold=fThreshold; + return true; +} + +// Set up speedup parameters +void AutoCleaner::speedup(const Float nDouble) +{ + itsDoSpeedup=true; + itsNDouble = nDouble; +}; + +// Do the clean as set up +Int AutoCleaner::clean(Matrix& model, + Bool /*showProgress*/) +{ + LogIO os(LogOrigin("AutoCleaner", "clean()", WHERE)); + + AlwaysAssert(model.shape()==itsDirty.shape(), AipsError); + itsdimensionsareeven = (model.shape()(0) == 2*(model.shape()(0)/2)); + AlwaysAssert((itsdimensionsareeven==false) || itsAutoHogbom, AipsError); + + Float tmpMaximumResidual=0.0; + + // Define a subregion for the inner quarter + IPosition blcDirty(model.shape().nelements(), 0); + IPosition trcDirty(model.shape()-1); + + // Handle Mask + if(!itsMask.null()){ + os << "Cleaning using given mask" << LogIO::POST; + if (itsMaskThreshold<0) { + os << LogIO::NORMAL + << "Mask thresholding is not used, values are interpreted as weights" + <shape()(0)==nx, AipsError); + AlwaysAssert(itsMask->shape()(1)==ny, AipsError); + Int xbeg=nx-1; + Int ybeg=ny-1; + Int xend=0; + Int yend=0; + for (Int iy=0;iy0.000001) { + xbeg=min(xbeg,ix); + ybeg=min(ybeg,iy); + xend=max(xend,ix); + yend=max(yend,iy); + } + } + } + + if (!itsIgnoreCenterBox) { + if((xend - xbeg)>nx/2) { + xbeg=nx/4-1; //if larger than quarter take inner of mask + os << LogIO::WARN << "Mask span over more than half the x-axis: Considering inner half of the x-axis" << LogIO::POST; + } + if((yend - ybeg)>ny/2) { + ybeg=ny/4-1; + os << LogIO::WARN << "Mask span over more than half the y-axis: Considering inner half of the y-axis" << LogIO::POST; + } + xend=min(xend,xbeg+nx/2-1); + yend=min(yend,ybeg+ny/2-1); + } + + //blcDirty(0)=xbeg> 0 ? xbeg-1 : 0; + //blcDirty(1)=ybeg > 0 ? ybeg-1 : 0; + blcDirty(0)=xbeg; + blcDirty(1)=ybeg; + trcDirty(0)=xend; + trcDirty(1)=yend; + } + else { + if (itsIgnoreCenterBox) { + os << LogIO::NORMAL << "Cleaning entire image" << LogIO::POST; + os << LogIO::NORMAL1 << "as per MF/WF" << LogIO::POST; // ??? + } + else { + os << "Cleaning inner quarter of the image" << LogIO::POST; + for (Int i=0;i vecWork_p; + vecWork_p.resize(gip); + casacore::Matrix vecWork_pms; + vecWork_pms.resize(gip); + Float tempGain; + Bool cleanhogbom=false; + // + + itsIteration = itsStartingIter; + for (Int ii=itsStartingIter; ii < itsMaxNiter; ii++) { + + itsIteration++; + // Find the peak residual + itsStrengthOptimum = 0.0; + + // Find absolute maximum for the dirty image + Matrix work = vecWork_p(blcDirty,trcDirty); + work = 0.0; + work = work + itsDirty(blcDirty,trcDirty); + + if (!itsMask.null()) { + findMaxAbsMask(vecWork_p, *itsMask, MaximumHogbom, posMaximumHogbom); + } else { + findMaxAbs(vecWork_p, MaximumHogbom, posMaximumHogbom); + } + + // Find absolute maximum for the ms dirty image + Matrix workms = vecWork_pms(blcDirty,trcDirty); + workms = 0.0; + workms = workms + tildeMI(blcDirty,trcDirty); + + if (!itsMask.null()) { + findMaxAbsMask(vecWork_pms, *itsMask, itsStrengthOptimum, posMaximum); + } else { + findMaxAbs(vecWork_pms, itsStrengthOptimum, posMaximum); + } + + if(abs(MaximumHogbom) > abs(itsStrengthOptimum)/normcmap || itsAutoHogbom){ + itsStrengthOptimum = MaximumHogbom; + posMaximum = posMaximumHogbom; + tempGain=itsHogbomGain/max(itsPsf); + cleanhogbom=true; + totalFlux += (itsStrengthOptimum*tempGain); + os << "Using Hogbom-CLEAN step"<< LogIO::POST; + } + else{ + tempGain=itsGain/max(tildeMB); + cleanhogbom=false; + if (abs(itsStrengthOptimum)*max(tildeMB)/max(MtildeMB) > abs(itsDirty(posMaximum))){ + itsStrengthOptimum*=max(tildeMB)/max(MtildeMB); + } + else{ + itsStrengthOptimum=itsDirty(posMaximum); + } + //totalFlux += (itsStrengthOptimum*tempGain*sumcmap); + os << "Using MS-CLEAN step"<< LogIO::POST; + } + + //trigger hogbom if we are stuck in a minimum + if(cleanhogbom==false){ + if(posMaximum == posMaximumlast){ + if(abs(abs(itsStrengthOptimum)-abs(StrengthOptimumlast))<0.1*abs(itsStrengthOptimum) && abs(itsStrengthOptimum+StrengthOptimumlast)<0.1*abs(itsStrengthOptimum)){ + os << "We trigger Hogbom step instead, since the MS-CLEAN step got stuck" << LogIO::POST; + itsStrengthOptimum = MaximumHogbom; + posMaximum = posMaximumHogbom; + tempGain=itsHogbomGain/max(itsPsf); + cleanhogbom=true; + totalFlux += (itsStrengthOptimum*tempGain); + os << "Using Hogbom-CLEAN step"<< LogIO::POST; + } + } + } + + if(cleanhogbom==false){ + totalFlux += (itsStrengthOptimum*tempGain*sumcmap); + } + + // Now add to the total flux + itsTotalFlux=totalFlux; + + // + Float scaleFactor; + scaleFactor=tempGain*itsStrengthOptimum; + + // Continuing: subtract the peak that we found from all dirty images + // Define a subregion so that that the peak is centered + IPosition support(model.shape()); + + IPosition inc(model.shape().nelements(), 1); + IPosition blc(posMaximum-support/2); + IPosition trc(posMaximum+support/2); + LCBox::verify(blc, trc, inc, model.shape()); + + IPosition blcPsf(blc+itsPositionPeakPsf-posMaximum); + IPosition trcPsf(trc+itsPositionPeakPsf-posMaximum); + LCBox::verify(blcPsf, trcPsf, inc, model.shape()); + makeBoxesSameSize(blc,trc,blcPsf,trcPsf); + if (itsdimensionsareeven){ + //blc(0) = posMaximum(0)-support(0)/2; + //blc(1) = posMaximum(1)-support(1)/2; + //trc(0) = posMaximum(0)-support(0)/2-1; + //trc(1) = posMaximum(1)-support(1)/2-1; + blc(posMaximum-support/2); + trc(posMaximum-support/2-1); + LCBox::verify(blc, trc, inc, model.shape()); + + //blcPsf(0) = blc(0)+itsPositionPeakPsf(0)-posMaximum(0); + //blcPsf(1) = blc(1)+itsPositionPeakPsf(1)-posMaximum(1); + //trcPsf(0) = trc(0)+itsPositionPeakPsf(0)-posMaximum(0); + //trcPsf(1) = trc(1)+itsPositionPeakPsf(1)-posMaximum(1); + blcPsf(blc+itsPositionPeakPsf-posMaximum); + trcPsf(trc+itsPositionPeakPsf-posMaximum); + LCBox::verify(blcPsf, trcPsf, inc, model.shape()); + makeBoxesSameSize(blc,trc,blcPsf,trcPsf); + } + + // The inverse shift indices + IPosition blcconj(support/2-posMaximum); + IPosition trcconj(support/2+support-posMaximum-1); + LCBox::verify(blcconj, trcconj, inc, model.shape()); + + IPosition blcPsfconj(blcconj+itsPositionPeakPsf+posMaximum-support+1); + IPosition trcPsfconj(trcconj+itsPositionPeakPsf+posMaximum-support+1); + LCBox::verify(blcPsfconj, trcPsfconj, inc, model.shape()); + makeBoxesSameSize(blcconj,trcconj,blcPsfconj,trcPsfconj); + + if(itsAutoHogbom){ + subtractBeam(itsDirty, itsPsf, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + subtractBeam(model, itsScales, blc, trc, blcPsf, trcPsf, scaleFactor, false, true); + } + else{ + if(cleanhogbom){ + triggerhogbom += itsHogbomGain; + Matrix BIold(BI.shape()); + BIold.assign_conforming(BI); + Matrix tildeII_copy1(tildeII.shape()), tildeII_copy2(tildeII.shape()); + tildeII_copy1 = 0; + tildeII_copy2 = 0; + #pragma omp parallel default(shared) num_threads(7) + { + #pragma omp master + { + + #pragma omp task + tildeII += scaleFactor*scaleFactor*BB; + + //Matrix tildeIISubc=tildeII(blcconj, trcconj); + ////Matrix tildeIIscaleSubc=tildeMBI(blcPsfconj,trcPsfconj); + //tildeIISubc -= scaleFactor*tildeMBI(blcPsfconj,trcPsfconj);//tildeIIscaleSubc; + #pragma omp task + subtractBeam(tildeII_copy1, BIold, blcconj, trcconj, blcPsfconj, trcPsfconj, scaleFactor, false, false); + + //Matrix tildeIISub=tildeII(blc, trc); + ////Matrix tildeIIscaleSubr = reverseArray(reverseArray(tildeMBI, 0), 1)(blcPsf,trcPsf); + //tildeIISub -= scaleFactor*reverseArray(reverseArray(tildeMBI, 0), 1)(blcPsf,trcPsf);//tildeIIscaleSubr; + #pragma omp task + subtractBeam(tildeII_copy2, BIold, blc, trc, blcPsf, trcPsf, scaleFactor, true, false); + + #pragma omp taskwait + #pragma omp task + tildeII += tildeII_copy1+tildeII_copy2; + + //Matrix tildeMBISub=tildeMBI(blc, trc); + ////Matrix tildeMBIscaleSub=MtildeMBB(blcPsf,trcPsf); + //tildeMBISub -= scaleFactor*MtildeMBB(blcPsf,trcPsf);//tildeMBIscaleSub; + //#pragma omp taskwait + #pragma omp task + subtractBeam(tildeMBI, tildeMBB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix dirtySub=itsDirty(blc,trc); + ////Matrix psfSub=tildeMB(blcPsf, trcPsf); + //dirtySub -= scaleFactor*tildeMB(blcPsf, trcPsf);//psfSub; + #pragma omp task + subtractBeam(itsDirty, itsPsf, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix tildeMISub=tildeMI(blc, trc); + ////Matrix tildeMIscaleSub=MtildeMB(blcPsf,trcPsf); + //tildeMISub -= scaleFactor*MtildeMB(blcPsf,trcPsf);//tildeMIscaleSub; + #pragma omp task + subtractBeam(tildeMI, tildeMB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix BISub=BI(blc, trc); + ////Matrix BIscaleSub=tildeMBB(blcPsf,trcPsf); + //BISub -= scaleFactor*tildeMBB(blcPsf,trcPsf);//BIscaleSub; + #pragma omp task + subtractBeam(BI, BB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //// Now do the addition of this scale to the model image.... + ////Matrix clean_mapSub=clean_map(blc, trc); + ////Matrix clean_mapscaleSub=mspsf(blcPsf,trcPsf); + ////clean_mapSub += scaleFactor*clean_mapscaleSub; + + //Matrix modelSub=model(blc, trc); + ////Matrix scaleSub=itsScales(blcPsf,trcPsf); //itsScale needs to be defined and is delta function, but there should be faster option here + //modelSub += scaleFactor*itsScales(blcPsf,trcPsf);//scaleSub; + #pragma omp task + subtractBeam(model, itsScales, blc, trc, blcPsf, trcPsf, scaleFactor/itsmaxbeam, false, true); + } + } + } + else{ + triggerhogbom = 0.0; + Matrix tildeMBIold(tildeMBI.shape()); + tildeMBIold.assign_conforming(tildeMBI); + Matrix tildeII_copy1(tildeII.shape()), tildeII_copy2(tildeII.shape()); + tildeII_copy1 = 0; + tildeII_copy2 = 0; + #pragma omp parallel default(shared) num_threads(7) + { + #pragma omp master + { + + #pragma omp task + tildeII += scaleFactor*scaleFactor*MtildeMBB; + + //Matrix tildeIISubc=tildeII(blcconj, trcconj); + ////Matrix tildeIIscaleSubc=tildeMBI(blcPsfconj,trcPsfconj); + //tildeIISubc -= scaleFactor*tildeMBI(blcPsfconj,trcPsfconj);//tildeIIscaleSubc; + #pragma omp task + subtractBeam(tildeII_copy1, tildeMBIold, blcconj, trcconj, blcPsfconj, trcPsfconj, scaleFactor, false, false); + + //Matrix tildeIISub=tildeII(blc, trc); + ////Matrix tildeIIscaleSubr = reverseArray(reverseArray(tildeMBI, 0), 1)(blcPsf,trcPsf); + //tildeIISub -= scaleFactor*reverseArray(reverseArray(tildeMBI, 0), 1)(blcPsf,trcPsf);//tildeIIscaleSubr; + #pragma omp task + subtractBeam(tildeII_copy2, tildeMBIold, blc, trc, blcPsf, trcPsf, scaleFactor, true, false); + + #pragma omp taskwait + #pragma omp task + tildeII += tildeII_copy1+tildeII_copy2; + + //Matrix tildeMBISub=tildeMBI(blc, trc); + ////Matrix tildeMBIscaleSub=MtildeMBB(blcPsf,trcPsf); + //tildeMBISub -= scaleFactor*MtildeMBB(blcPsf,trcPsf);//tildeMBIscaleSub; + #pragma omp task + subtractBeam(tildeMBI, MtildeMBB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix dirtySub=itsDirty(blc,trc); + ////Matrix psfSub=tildeMB(blcPsf, trcPsf); + //dirtySub -= scaleFactor*tildeMB(blcPsf, trcPsf);//psfSub; + #pragma omp task + subtractBeam(itsDirty, tildeMB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix tildeMISub=tildeMI(blc, trc); + ////Matrix tildeMIscaleSub=MtildeMB(blcPsf,trcPsf); + //tildeMISub -= scaleFactor*MtildeMB(blcPsf,trcPsf);//tildeMIscaleSub; + #pragma omp task + subtractBeam(tildeMI, MtildeMB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //Matrix BISub=BI(blc, trc); + ////Matrix BIscaleSub=tildeMBB(blcPsf,trcPsf); + //BISub -= scaleFactor*tildeMBB(blcPsf,trcPsf);//BIscaleSub; + #pragma omp task + subtractBeam(BI, tildeMBB, blc, trc, blcPsf, trcPsf, scaleFactor, false, false); + + //// Now do the addition of this scale to the model image.... + ////Matrix clean_mapSub=clean_map(blc, trc); + ////Matrix clean_mapscaleSub=mspsf(blcPsf,trcPsf); + ////clean_mapSub += scaleFactor*clean_mapscaleSub; + + //Matrix modelSub=model(blc, trc); + ////Matrix scaleSub=itsScales(blcPsf,trcPsf); //itsScale needs to be defined and is delta function, but there should be faster option here + //modelSub += scaleFactor*itsScales(blcPsf,trcPsf);//scaleSub; + #pragma omp task + subtractBeam(model, cmap, blc, trc, blcPsf, trcPsf, scaleFactor/itsmaxbeam, false, true); + } + } + } + + updateBasisFunction(); + + blcDirty = blc; + trcDirty = trc; + + //update the last stored position for multiscale variant + if(cleanhogbom==false){ + posMaximumlast = posMaximum; + StrengthOptimumlast = itsStrengthOptimum; + } + } + + os << itsDirty(posMaximum) << " " << MaximumHogbom-scaleFactor*max(itsPsf) << LogIO::POST; + os << itsIteration << " " << itsStrengthOptimum << " " << totalFlux << LogIO::POST; + if(triggerhogbom>itsAutoTrigger){ + os << "we switch to Hogbom CLEAN permanently" << LogIO::POST; + itsAutoHogbom=true; + } + } + return 1; + // End of iteration +} + + Bool AutoCleaner::findPSFMaxAbs(const Matrix& lattice, + Float& maxAbs, + IPosition& posMaxAbs, + const Int& supportSize) + { + LogIO os(LogOrigin("AutoCleaner", "findPSFMaxAbs()", WHERE)); + posMaxAbs = IPosition(lattice.shape().nelements(), 0); + maxAbs=0.0; + + Float maxVal=0.0; + IPosition psf2DShape(lattice.shape()); + Int blc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 - supportSize/2 : 0; + Int blc1 = (psf2DShape(1) > supportSize) ? psf2DShape(1)/2 - supportSize/2 : 0; + Int trc0 = (psf2DShape(0) > supportSize) ? psf2DShape(0)/2 + supportSize/2 : psf2DShape(0)-1; + + Int trc1 = (psf2DShape(1) > supportSize) ? (psf2DShape(1)/2 + supportSize/2) : (psf2DShape(1)-1) ; + + + // cerr << "####### " << blc0 << " " << blc1 << " " << trc0 << " " << trc1 << endl; + // cerr << "Max of lattice " << max(lattice) << " min " << min(lattice) << endl; + for (Int j=blc1; j < trc1; ++j) + for (Int i=blc0 ; i < trc0; ++i) + if ((maxAbs = abs(lattice(i,j))) > maxVal) + { + maxVal = maxAbs; + posMaxAbs(0)=i; posMaxAbs(1)=j; + } + maxAbs=maxVal; + //cerr << "######## " << posMaxAbs << " " << maxVal << endl; + return true; + } + +Bool AutoCleaner::findMaxAbs(const Matrix& lattice, + Float& maxAbs, + IPosition& posMaxAbs) +{ + + posMaxAbs = IPosition(lattice.shape().nelements(), 0); + maxAbs=0.0; + + Float minVal; + IPosition posmin(lattice.shape().nelements(), 0); + minMax(minVal, maxAbs, posmin, posMaxAbs, lattice); + //cout << "min " << minVal << " " << maxAbs << " " << max(lattice) << endl; + if(abs(minVal) > abs(maxAbs)){ + maxAbs=minVal; + posMaxAbs=posmin; + } + return true; +} + +Bool AutoCleaner::findMaxAbsMask(const Matrix& lattice, + const Matrix& mask, + Float& maxAbs, + IPosition& posMaxAbs) +{ + + posMaxAbs = IPosition(lattice.shape().nelements(), 0); + maxAbs=0.0; + Float minVal; + IPosition posmin(lattice.shape().nelements(), 0); + minMaxMasked(minVal, maxAbs, posmin, posMaxAbs, lattice, mask); + if(abs(minVal) > abs(maxAbs)){ + maxAbs=minVal; + posMaxAbs=posmin; + } + + return true; +} + +void AutoCleaner::initializeCorrProducts() +{ + //BB = correlate(itsPsf, itsPsf); + //tildeII = correlate(itsDirty, itsDirty); + //BI = correlate(itsPsf, itsDirty); + + LogIO os(LogOrigin("AutoCleaner", "initializeCorrProducts()", WHERE)); + //AlwaysAssert(!itsDirty.null(), AipsError); + + Float max_beam = max(itsPsf); + itsPsf /= max_beam; + + FFTServer fft(itsDirty.shape()); + Matrix FTPsf; + Matrix FTDirty; + Matrix cWork; + + fft.fft0(FTPsf, itsPsf); + cWork=((FTPsf)*(conj(FTPsf))); + fft.fft0(BB, cWork, false); + fft.flip(BB, false, false); + + max_beam = max(BB); + BB /= max_beam; + itsPsf /= sqrt(max_beam); + itsmaxbeam = sqrt(max_beam); + + fft.fft0(FTDirty, itsDirty); + cWork=((FTDirty)*(conj(FTDirty))); + fft.fft0(tildeII, cWork, false); + fft.flip(tildeII, false, false); + + cWork=((FTDirty)*(FTPsf)); + fft.fft0(BI, cWork, false); + fft.flip(BI, false, false); + + BI /= sqrt(max_beam); + + os << "Correlation Products initialized" << LogIO::POST; +} + +void AutoCleaner::reinitializeCorrProducts() +{ + //BB = correlate(itsPsf, itsPsf); + //tildeII = correlate(itsDirty, itsDirty); + //BI = correlate(itsPsf, itsDirty); + + LogIO os(LogOrigin("AutoCleaner", "reinitializeCorrProducts()", WHERE)); + //AlwaysAssert(!itsDirty.null(), AipsError); + + Float max_beam = max(itsPsf); + itsPsf /= max_beam; + + FFTServer fft(itsDirty.shape()); + Matrix FTPsf; + Matrix FTDirty; + Matrix FTcmap; + Matrix FTBI; + Matrix cWork; + + fft.fft0(FTPsf, itsPsf); + cWork=((FTPsf)*(conj(FTPsf))); + fft.fft0(BB, cWork, false); + fft.flip(BB, false, false); + + max_beam = max(BB); + BB /= max_beam; + itsPsf /= sqrt(max_beam); + itsmaxbeam = sqrt(max_beam); + + fft.fft0(FTDirty, itsDirty); + cWork=((FTDirty)*(conj(FTDirty))); + fft.fft0(tildeII, cWork, false); + fft.flip(tildeII, false, false); + + cWork=((FTDirty)*(FTPsf)); + fft.fft0(BI, cWork, false); + fft.flip(BI, false, false); + + BI /= sqrt(max_beam); + + fft.fft0(FTcmap, cmap); + cWork=((FTDirty)*(conj(FTcmap))); + fft.fft0(tildeMI, cWork, false); + fft.flip(tildeMI, false, false); + + fft.fft0(FTBI, BI); + cWork=((FTBI)*(conj(FTcmap))); + fft.fft0(tildeMBI, cWork, false); + fft.flip(tildeMBI, false, false); + + //cWork=((FTDirty)*(FTPsf)*(conj(FTcmap))); + //fft.fft0(tildeMBI, cWork, false); + //fft.flip(tildeMBI, false, false); + + //tildeMBI /= sqrt(max_beam); + + casacore::MArray Mcmap(cmap); + normcmap = norm(Mcmap.flatten()); + normcmapsq = normcmap*normcmap; + sumcmap = sum(cmap); + + //Find mask for the autocorrelation by searching for first sidelobe around central peak + Int ix=0; + Int iy=0; + Int ixy=0; + Bool searchxx=true; + Bool searchyy=true; + Bool searchxy=true; + IPosition support(itsDirty.shape()); + + IPosition posxx(support/2); + IPosition posyy(support/2); + IPosition posxy(support/2); + if(itsAutoXMask == 0 && itsAutoYMask == 0){ + for(Int is=0; is < support(0); is++){ + posxx(0)=support(0)/2+is; + posyy(1)=support(1)/2+is; + posxy(0)=support(0)/2+is; + posxy(1)=support(1)/2+is; + if(tildeII(posxx)<0 && searchxx){ + ix=is; + searchxx=false; + } + if(tildeII(posyy)<0 && searchyy){ + iy=is; + searchyy=false; + } + if(tildeII(posxy)<0 && searchxy){ + ixy=is; + searchxy=false; + } + if(searchxx==false && searchyy==false && searchxy==false){ + break; + } + } + + //blcautomask(support); + //trcautomask(support); + if(ix fft(itsDirty.shape()); + Matrix FTPsf; + Matrix FTDirty; + Matrix FTcmap; + Matrix FTBI; + Matrix FTBB; + Matrix FTMB; + Matrix FTMBB; + Matrix cWork; + + fft.fft0(FTPsf, itsPsf); + fft.fft0(FTcmap, cmap); + cWork=((FTPsf)*(conj(FTcmap))); + fft.fft0(tildeMB, cWork, false); + fft.flip(tildeMB, false, false); + + fft.fft0(FTBB, BB); + cWork=((FTBB)*(conj(FTcmap))); + fft.fft0(tildeMBB, cWork, false); + fft.flip(tildeMBB, false, false); + + fft.fft0(FTDirty, itsDirty); + cWork=((FTDirty)*(conj(FTcmap))); + fft.fft0(tildeMI, cWork, false); + fft.flip(tildeMI, false, false); + + fft.fft0(FTBI, BI); + cWork=((FTBI)*(conj(FTcmap))); + fft.fft0(tildeMBI, cWork, false); + fft.flip(tildeMBI, false, false); + + fft.fft0(FTMB, tildeMB); + cWork=((FTMB)*(FTcmap)); + fft.fft0(MtildeMB, cWork, false); + fft.flip(MtildeMB, false, false); + + fft.fft0(FTMBB, tildeMBB); + cWork=((FTMBB)*(FTcmap)); + fft.fft0(MtildeMBB, cWork, false); + fft.flip(MtildeMBB, false, false); + + os << "Correlation Products recomputed" << LogIO::POST; +} + +void AutoCleaner::approximateBasisFunction() +{ + LogIO os(LogOrigin("AutoCleaner", "approximateBasisFunction()", WHERE)); + os << "Start to clean autocorrelation function" << LogIO::POST; + Float gain=itsAutoGain; + Float power=itsAutoPower; + Int niter=itsAutoMaxiter; + Int ii=0; + Matrix dmap(itsPsf.shape()); + dmap = 0.0; + dmap = dmap+tildeII; + + //Find mask for the autocorrelation by searching for first sidelobe around central peak + Int ix=0; + Int iy=0; + Int ixy=0; + Bool searchxx=true; + Bool searchyy=true; + Bool searchxy=true; + IPosition support(dmap.shape()); + + IPosition posxx(support/2); + IPosition posyy(support/2); + IPosition posxy(support/2); + + if(itsAutoXMask == 0 && itsAutoYMask == 0){ + for(Int is=0; is < support(0); is++){ + posxx(0)=support(0)/2+is; + posyy(1)=support(1)/2+is; + posxy(0)=support(0)/2+is; + posxy(1)=support(1)/2+is; + if(tildeII(posxx)<0 && searchxx){ + ix=is; + searchxx=false; + } + if(tildeII(posyy)<0 && searchyy){ + iy=is; + searchyy=false; + } + if(tildeII(posxy)<0 && searchxy){ + ixy=is; + searchxy=false; + } + if(searchxx==false && searchyy==false && searchxy==false){ + break; + } + } + + //blcautomask(support); + //trcautomask(support); + if(ix strengths(niter); + Vector xindex(niter); + Vector yindex(niter); + + strengths=0.0; + xindex=0; + yindex=0; + + Float threshmaxtildeII = itsAutoThreshold*max(abs(tildeII(blcautomask,trcautomask))); + + //for (Int ii=0; ii < niter; ii++) { + while(max(abs(dmap(blcautomask,trcautomask))) > threshmaxtildeII && ii abs(minval) || window_basis(minpos)!=1.0){ + index = maxpos+blcautomask; + strength = maxval; + } + else{ + index = minpos+blcautomask; + strength = minval; + } + + scaleFactor=gain*strength; + scaleFactorpow=sign(scaleFactor)*pow(abs(scaleFactor),power); + + xindex(ii) = index(0); + yindex(ii) = index(1); + strengths(ii) = strength; + + // Continuing: subtract the peak that we found from all dirty images + // Define a subregion so that that the peak is centered + + IPosition inc(dmap.shape().nelements(), 1); + IPosition blc(index-support/2); + IPosition trc(index+support/2); + LCBox::verify(blc, trc, inc, dmap.shape()); + + IPosition blcBB(blc+itsPositionPeakPsf-index); + IPosition trcBB(trc+itsPositionPeakPsf-index); + LCBox::verify(blcBB, trcBB, inc, dmap.shape()); + makeBoxesSameSize(blc,trc,blcBB,trcBB); + + // The inverse shift indices + IPosition blcconj(support/2-index); + IPosition trcconj(support/2+support-index-1); + LCBox::verify(blcconj, trcconj, inc, dmap.shape()); + + IPosition blcBBconj(blcconj+itsPositionPeakPsf+index-support+1); + IPosition trcBBconj(trcconj+itsPositionPeakPsf+index-support+1); + LCBox::verify(blcBBconj, trcBBconj, inc, dmap.shape()); + makeBoxesSameSize(blcconj,trcconj,blcBBconj,trcBBconj); + + normcmapsq -= cmap(index)*cmap(index); + normcmapsq -= cmap(support-index)*cmap(support-index); + sumcmap += 2*scaleFactor; + + #pragma omp parallel default(shared) num_threads(14) + { + #pragma omp master + { + // Substract the beam + #pragma omp task + subtractBeam(dmap, BB, blc, trc, blcBB, trcBB, scaleFactor, false, false); + + // Now do the addition of this scale to the model image.... + #pragma omp task + subtractBeam(cmap, itsScales, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(mod, BB, blc, trc, blcBB, trcBB, scaleFactor, false, true); + + // Build the channel products + #pragma omp task + subtractBeam(tildeMB, itsPsf, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMI, itsDirty, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBB, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBI, BI, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + // Redo the same for the conjugate position + #pragma omp taskwait + #pragma omp task + subtractBeam(dmap, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, false); + + // Add scale to image + #pragma omp task + subtractBeam(cmap, itsScales, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(mod, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, true); + + #pragma omp task + subtractBeam(tildeMB, itsPsf, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMI, itsDirty, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBB, BB, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBI, BI, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + } + + } + window_basis(index) = 1.0; + window_basis(support-index) = 1.0; + normcmapsq += cmap(index)*cmap(index); + normcmapsq += cmap(support-index)*cmap(support-index); + ii += 1; + } + + normcmap = sqrtop(normcmapsq); + niter = ii; + + // continue with the build of channel products that need precomputed quantities + for (Int i=0; i < niter; i++) { + // Read stored positions and strength + index(0)=xindex(i); + index(1)=yindex(i); + strength=strengths(i); + scaleFactor=gain*strength; + scaleFactorpow=sign(scaleFactor)*pow(abs(scaleFactor),power); + + // Define a subregion so that that the peak is centered + IPosition support(dmap.shape()); + + IPosition inc(dmap.shape().nelements(), 1); + IPosition blc(index-support/2); + IPosition trc(index+support/2); + LCBox::verify(blc, trc, inc, dmap.shape()); + + IPosition blcBB(blc+itsPositionPeakPsf-index); + IPosition trcBB(trc+itsPositionPeakPsf-index); + LCBox::verify(blcBB, trcBB, inc, dmap.shape()); + makeBoxesSameSize(blc,trc,blcBB,trcBB); + + // The inverse shift indices + IPosition blcconj(support/2-index); + IPosition trcconj(support/2+support-index-1); + LCBox::verify(blcconj, trcconj, inc, dmap.shape()); + + IPosition blcBBconj(blcconj+itsPositionPeakPsf+index-support+1); + IPosition trcBBconj(trcconj+itsPositionPeakPsf+index-support+1); + LCBox::verify(blcBBconj, trcBBconj, inc, dmap.shape()); + makeBoxesSameSize(blcconj,trcconj,blcBBconj,trcBBconj); + + #pragma omp parallel default(shared) num_threads(4) + { + #pragma omp master + { + // these two blocks need to be executed after the complete iterations have been done on the previous quantities + //Matrix MtildeMBBSub=MtildeMBB(blc, trc); + ////Matrix MtildeMBBscaleSub=tildeMBB(blcBB,trcBB); + //MtildeMBBSub += scaleFactor*tildeMBB(blcBB,trcBB);//MtildeMBBscaleSub; + #pragma omp task + subtractBeam(MtildeMBB, tildeMBB, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + //Matrix MtildeMBSub=MtildeMB(blc, trc); + ////Matrix MtildeMBscaleSub=tildeMB(blcBB,trcBB); + //MtildeMBSub += scaleFactor*tildeMB(blcBB,trcBB);//MtildeMBscaleSub; + #pragma omp task + subtractBeam(MtildeMB, tildeMB, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + // do the conjugate + //Matrix MtildeMBBSubc=MtildeMBB(blcconj, trcconj); + ////Matrix MtildeMBBscaleSubc=tildeMBB(blcBBconj,trcBBconj); + //MtildeMBBSubc += scaleFactor*tildeMBB(blcBBconj,trcBBconj);//MtildeMBBscaleSubc; + #pragma omp taskwait + #pragma omp task + subtractBeam(MtildeMBB, tildeMBB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + //Matrix MtildeMBSubc=MtildeMB(blcconj, trcconj); + ////Matrix MtildeMBscaleSubc=tildeMB(blcBBconj,trcBBconj); + //MtildeMBSubc += scaleFactor*tildeMB(blcBBconj,trcBBconj);//MtildeMBscaleSubc; + #pragma omp task + subtractBeam(MtildeMB, tildeMB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + } + } + } + os << "we used " << niter << " initial iterations" << LogIO::POST; + os << " to achieve at a threshold " << itsAutoThreshold << LogIO::POST; +} + +void AutoCleaner::updateBasisFunction() +{ + LogIO os(LogOrigin("AutoCleaner", "updateBasisFunction()", WHERE)); + os << "Update the Autocorrelation model" << LogIO::POST; + + Float normalize; + normalize = max(abs(tildeII))/max(abs(mod)); + Float normalizepow; + normalizepow=pow(normalize,itsAutoPower); + + cmap *= normalizepow; + mod *= normalize; + tildeMB *= normalizepow; + tildeMI *= normalizepow; + tildeMBB *= normalizepow; + MtildeMBB *= normalizepow*normalizepow; + tildeMBI *= normalizepow; + MtildeMB *= normalizepow*normalizepow; + + normcmapsq *= normalizepow*normalizepow; + normcmap *= normalizepow; + sumcmap *= normalizepow; + + Matrix dmap(itsPsf.shape()); + dmap = 0.0; + dmap = dmap+tildeII-mod; + Float gain=itsAutoGain; + + Float strength; + IPosition index(dmap.shape()); + Float scaleFactor; + Float scaleFactorpow; + + Float threshmaxtildeII = itsAutoThreshold*max(abs(tildeII(blcautomask,trcautomask))); + Int niter=itsAutoMaxiter; + Int ii=0; + + Vector strengths(niter); + Vector xindex(niter); + Vector yindex(niter); + + strengths=0.0; + xindex=0; + yindex=0; + + while(max(abs(dmap(blcautomask,trcautomask))) > threshmaxtildeII && ii abs(minval) || window_basis(minpos)!=1.0){ + index = maxpos+blcautomask; + strength = maxval; + } + else{ + index = minpos+blcautomask; + strength = minval; + } + + scaleFactor=gain*strength; + scaleFactorpow=sign(scaleFactor)*pow(abs(scaleFactor),itsAutoPower); + + // Define a subregion so that that the peak is centered + IPosition support(dmap.shape()); + + IPosition inc(dmap.shape().nelements(), 1); + IPosition blc(index-support/2); + IPosition trc(index+support/2); + LCBox::verify(blc, trc, inc, dmap.shape()); + + IPosition blcBB(blc+itsPositionPeakPsf-index); + IPosition trcBB(trc+itsPositionPeakPsf-index); + LCBox::verify(blcBB, trcBB, inc, dmap.shape()); + makeBoxesSameSize(blc,trc,blcBB,trcBB); + + // The inverse shift indices + IPosition blcconj(support/2-index); + IPosition trcconj(support/2+support-index-1); + LCBox::verify(blcconj, trcconj, inc, dmap.shape()); + + IPosition blcBBconj(blcconj+itsPositionPeakPsf+index-support+1); + IPosition trcBBconj(trcconj+itsPositionPeakPsf+index-support+1); + LCBox::verify(blcBBconj, trcBBconj, inc, dmap.shape()); + makeBoxesSameSize(blcconj,trcconj,blcBBconj,trcBBconj); + + normcmapsq -= cmap(index)*cmap(index); + normcmapsq -= cmap(support-index)*cmap(support-index); + + sumcmap += 2*scaleFactor; + + xindex(ii) = index(0); + yindex(ii) = index(1); + strengths(ii) = strength; + + #pragma omp parallel default(shared) num_threads(3) + { + #pragma omp master + { + #pragma omp task + subtractBeam(dmap, BB, blc, trc, blcBB, trcBB, scaleFactor, false, false); + #pragma omp task + subtractBeam(cmap, itsScales, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + #pragma omp task + subtractBeam(mod, BB, blc, trc, blcBB, trcBB, scaleFactor, false, true); + + #pragma omp taskwait + #pragma omp task + subtractBeam(dmap, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, false); + #pragma omp task + subtractBeam(cmap, itsScales, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + #pragma omp task + subtractBeam(mod, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, true); + } + } + + window_basis(index) = 1.0; + window_basis(support-index) = 1.0; + + normcmapsq += cmap(index)*cmap(index); + normcmapsq += cmap(support-index)*cmap(support-index); + + ii+=1; + } + + normcmap = sqrtop(normcmapsq); + niter = ii; + + //for just a few iterations, the fastest version is to track down the change in the correlation products by an analytic description + // for more iterations, it makes more sense to just redo the convolution operation + if(niter<10){ + for (Int i=0; i < niter; i++) { + // Read stored positions and strength + index(0)=xindex(i); + index(1)=yindex(i); + strength=strengths(i); + scaleFactor=gain*strength; + scaleFactorpow=sign(scaleFactor)*pow(abs(scaleFactor),itsAutoPower); + + // Define a subregion so that that the peak is centered + IPosition support(dmap.shape()); + + IPosition inc(dmap.shape().nelements(), 1); + IPosition blc(index-support/2); + IPosition trc(index+support/2); + LCBox::verify(blc, trc, inc, dmap.shape()); + + IPosition blcBB(blc+itsPositionPeakPsf-index); + IPosition trcBB(trc+itsPositionPeakPsf-index); + LCBox::verify(blcBB, trcBB, inc, dmap.shape()); + makeBoxesSameSize(blc,trc,blcBB,trcBB); + + // The inverse shift indices + IPosition blcconj(support/2-index); + IPosition trcconj(support/2+support-index-1); + LCBox::verify(blcconj, trcconj, inc, dmap.shape()); + + IPosition blcBBconj(blcconj+itsPositionPeakPsf+index-support+1); + IPosition trcBBconj(trcconj+itsPositionPeakPsf+index-support+1); + LCBox::verify(blcBBconj, trcBBconj, inc, dmap.shape()); + makeBoxesSameSize(blcconj,trcconj,blcBBconj,trcBBconj); + + // store old products for parallel computing + Matrix tildeMBold(tildeMB.shape()); + tildeMBold.assign_conforming(tildeMB); + Matrix tildeMBBold(tildeMBB.shape()); + tildeMBBold.assign_conforming(tildeMBB); + + #pragma omp parallel default(shared) num_threads(6) + { + #pragma omp master + { + + //#pragma omp task + //subtractBeam(dmap, BB, blc, trc, blcBB, trcBB, scaleFactor, false, false); + + // Add scale to image + //#pragma omp task + //subtractBeam(cmap, itsScales, blc, trc, blcBB, trcBB, scaleFactor, false, true); + + //#pragma omp task + //subtractBeam(mod, BB, blc, trc, blcBB, trcBB, scaleFactor, false, true); + + //update channel products, first the beams + #pragma omp task + MtildeMB += scaleFactorpow*scaleFactorpow*itsPsf; + #pragma omp task + MtildeMBB += scaleFactorpow*scaleFactorpow*BB; + + //wait for task to finish before writing to MtildeMB again + #pragma omp taskwait + #pragma omp task + subtractBeam(MtildeMB, tildeMBold, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(MtildeMBB, tildeMBBold, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMB, itsPsf, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBB, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + //Now the maps + #pragma omp task + subtractBeam(tildeMI, itsDirty, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBI, BI, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + // Add the conjugate scale to image as well + //#pragma omp task + //subtractBeam(cmap, itsScales, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, true); + + //#pragma omp task + //subtractBeam(mod, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, true); + + //previous tasks have to finish first to avoid overwriting of MtildeMB and MtildeMBB + #pragma omp taskwait + #pragma omp task + subtractBeam(MtildeMB, tildeMBold, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, true, true); + + #pragma omp task + subtractBeam(MtildeMBB, tildeMBBold, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, true, true); + + // update channel products with the conjugate as well + #pragma omp taskwait + #pragma omp task + {tildeMBold.assign_conforming(tildeMB); + tildeMBBold.assign_conforming(tildeMBB);} + + #pragma omp task + MtildeMB += scaleFactorpow*scaleFactorpow*itsPsf; + #pragma omp task + MtildeMBB += scaleFactorpow*scaleFactorpow*BB; + + //These two conjugates cannot be added above with the reverse, since we need the update of tildeMB and tildeMBB first + #pragma omp taskwait + #pragma omp task + subtractBeam(MtildeMB, tildeMBold, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(MtildeMBB, tildeMBBold, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMB, itsPsf, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBB, BB, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + //Now the maps + #pragma omp task + subtractBeam(tildeMI, itsDirty, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + #pragma omp task + subtractBeam(tildeMBI, BI, blc, trc, blcBB, trcBB, scaleFactorpow, false, true); + + //#pragma omp task + //subtractBeam(dmap, BB, blcconj, trcconj, blcBBconj, trcBBconj, scaleFactor, false, false); + + //wait for task to finish to prevent overwriting + #pragma omp taskwait + #pragma omp task + subtractBeam(MtildeMB, tildeMBold, blc, trc, blcBB, trcBB, scaleFactorpow, true, true); + + #pragma omp task + subtractBeam(MtildeMBB, tildeMBBold, blc, trc, blcBB, trcBB, scaleFactorpow, true, true); + } + } + } + } + else{ + recomputeCorrProducts(); + } + basis_function = cmap; + os << "Using " << niter << " iterations for that" << LogIO::POST; +} + +void AutoCleaner::subtractBeam(Matrix &map, Matrix &beam, IPosition blc, IPosition trc, IPosition blcbeam, IPosition trcbeam, Float factor, Bool reverse, Bool add) +{ + Matrix mapSub=map(blc,trc); + if(reverse){ + if(add){ + mapSub += factor*reverseArray(reverseArray(beam, 0), 1)(blcbeam,trcbeam); + } + else{ + mapSub -= factor*reverseArray(reverseArray(beam, 0), 1)(blcbeam,trcbeam); + } + } + else{ + if(add){ + mapSub += factor*beam(blcbeam,trcbeam); + } + else{ + mapSub -= factor*beam(blcbeam,trcbeam); + } + } +} + +Bool AutoCleaner::setscales() +{ + Int nx=itsPsf.shape()(0); + Int ny=itsPsf.shape()(1); + + Double refi=nx/2; + Double refj=ny/2; + + //itsScales=new Matrix(itsDirty.shape()); + Matrix itsScales_c(nx,ny); + //itsScales(itsPsf.shape()); + itsScales_c = 0.0; + + itsScales_c(Int(refi), Int(refj)) = 1.0; + + itsScales.resize(itsPsf.shape()); + itsScales.assign(itsScales_c); + + return true; +} + +void AutoCleaner::setDirty(const Matrix& dirty){ + + //itsDirty=new Matrix(dirty.shape()); + //Matrix itsDirty; + itsDirty.resize(dirty.shape()); + itsDirty.assign(dirty); +} + +//void AutoCleaner::setPsf(const Matrix& psf){ +// //Matrix itsPsf; +// itsPsf(psf.shape()); +// itsPsf.assign(psf); +//} + + +void AutoCleaner::unsetMask() +{ + if(!itsMask.null()) + itsMask=0; + noClean_p=false; +} + +Float AutoCleaner::threshold() const +{ + if (! itsDoSpeedup) { + return max(itsFracThreshold.get("%").getValue() * itsMaximumResidual /100.0, + itsThreshold.get("Jy").getValue()); + } else { + const Float factor = exp( (Float)( itsIteration - itsStartingIter )/ itsNDouble ) + / 2.7182818; + return factor * max(itsFracThreshold.get("%").getValue() * itsMaximumResidual /100.0, + itsThreshold.get("Jy").getValue()); + } +} + +void AutoCleaner::makeBoxesSameSize(IPosition& blc1, IPosition& trc1, + IPosition &blc2, IPosition& trc2) +{ + const IPosition shape1 = trc1 - blc1; + const IPosition shape2 = trc2 - blc2; + + AlwaysAssert(shape1.nelements() == shape2.nelements(), AipsError); + + if (shape1 == shape2) { + return; + } + for (uInt i=0;i=0, AipsError); + //if (minLength % 2 != 0) { + // if the number of pixels is odd, ensure that the centre stays + // the same by making this number even + //--minLength; // this code is a mistake and should be removed + //} + const Int increment1 = shape1[i] - minLength; + const Int increment2 = shape2[i] - minLength; + blc1[i] += increment1/2; + trc1[i] -= increment1/2 + (increment1 % 2 != 0 ? 1 : 0); + blc2[i] += increment2/2; + trc2[i] -= increment2/2 + (increment2 % 2 != 0 ? 1 : 0); + } +} + +Int AutoCleaner::findBeamPatch(const Float maxScaleSize, const Int& nx, const Int& ny, + const Float psfBeam, const Float nBeams) +{ + Int psupport = (Int) ( sqrt( psfBeam*psfBeam + maxScaleSize*maxScaleSize ) * nBeams ); + + // At least this big... + if(psupport < psfBeam*nBeams ) psupport = static_cast(psfBeam*nBeams); + + // Not too big... + if(psupport > nx || psupport > ny) psupport = std::min(nx,ny); + + // Make it even. + if (psupport%2 != 0) psupport -= 1; + + return psupport; +} + + + +} //# NAMESPACE CASA - END diff --git a/src/synthesis/MeasurementEquations/AutoCleaner.h b/src/synthesis/MeasurementEquations/AutoCleaner.h new file mode 100644 index 00000000..c187bf91 --- /dev/null +++ b/src/synthesis/MeasurementEquations/AutoCleaner.h @@ -0,0 +1,370 @@ + +//# Cleaner.h: this defines Cleaner a class for doing convolution +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id: LatticeCleaner.h 20739 2009-09-29 01:15:15Z Malte.Marquarding $ + +#ifndef SYNTHESIS_AUTOCLEANER_H +#define SYNTHESIS_AUTOCLEANER_H + +//# Includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + +//# Forward Declarations +class MatrixNACleaner; +// A copy of casacore::LatticeCleaner but just using 2-D matrices +// It is a mimic of the casacore::LatticeCleaner class but avoid a lot of +// of the lattice to array and back copies and uses openmp in the obvious places +// + +// A class for doing multi-dimensional cleaning + +// + +// +// + +// +//
  • The mathematical concept of deconvolution +// +// +// + +// The MatrixCleaner class will deconvolve 2-D arrays of floats. + +// +// +// +// This class will perform various types of Clean deconvolution +// on Lattices. +// +// +// +// +// +// +// +// +// +// +// +// +//
  • casacore::AipsError: if psf has more dimensions than the model. +// +// +// +// + +class AutoCleaner +{ +public: + + // Create a cleaner : default constructor + AutoCleaner(); + + // Create a cleaner for a specific dirty image and PSF + AutoCleaner(const casacore::Matrix & psf, const casacore::Matrix & dirty); + + // The copy constructor uses reference semantics + AutoCleaner(const AutoCleaner& other); + + // The assignment operator also uses reference semantics + AutoCleaner & operator=(const AutoCleaner & other); + + // The destructor does nothing special. + ~AutoCleaner(); + + //Set the dirty image without calculating convolutions.. + //can be done by calling makeDirtyScales or setscales if one want to redo the + //psfscales too. + void setDirty(const casacore::Matrix& dirty); + // Update the dirty image only (equiv of setDirty + makeDirtyScales) + void update(const casacore::Matrix & dirty); + + casacore::Bool setscales(); + void initializeCorrProducts(); + void reinitializeCorrProducts(); + void recomputeCorrProducts(); + void approximateBasisFunction(); + void updateBasisFunction(); + void subtractBeam(casacore::Matrix &map, casacore::Matrix &beam, casacore::IPosition blc, casacore::IPosition trc, casacore::IPosition blcbeam, casacore::IPosition trcbeam, casacore::Float factor, casacore::Bool reverse, casacore::Bool add); + + + //change the psf + //don't forget to redo the setscales or run makePsfScales, + //followed by makeDirtyScales + void setPsf(const casacore::Matrix& psf); + + // Set up control parameters + // cleanType - type of the cleaning algorithm to use (HOGBOM, MULTISCALE) + // niter - number of iterations + // gain - loop gain used in cleaning (a fraction of the maximum + // subtracted at every iteration) + // aThreshold - absolute threshold to stop iterations + // fThreshold - fractional threshold (i.e. given w.r.t. maximum residual) + // to stop iterations. This parameter is specified as + // casacore::Quantity so it can be given in per cents. + // choose - unused at the moment, specify false. Original meaning is + // to allow interactive decision on whether to continue iterations. + // This method always returns true. + casacore::Bool setcontrol(casacore::CleanEnums::CleanType cleanType, const casacore::Int niter, + const casacore::Float gain, const casacore::Quantity& aThreshold, + const casacore::Quantity& fThreshold); + + // This version of the method disables stopping on fractional threshold + casacore::Bool setcontrol(casacore::CleanEnums::CleanType cleanType, const casacore::Int niter, + const casacore::Float gain, const casacore::Quantity& threshold); + + // return how many iterations we did do + casacore::Int iteration() const { return itsIteration; } + casacore::Int numberIterations() const { return itsIteration; } + + // what iteration number to start on + void startingIteration(const casacore::Int starting = 0) {itsStartingIter = starting; } + + //Total flux accumulated so far + casacore::Float totalFlux() const {return itsTotalFlux;} + + + // Clean an image. + //return value gives you a hint of what's happening + // 1 = converged + // 0 = not converged but behaving normally + // -1 = not converged and stopped on cleaning consecutive smallest scale + // -2 = not converged and either large scale hit negative or diverging + // -3 = clean is diverging rather than converging + casacore::Int clean(casacore::Matrix & model, casacore::Bool doPlotProgress=false); + + // Set the mask + // mask - input mask lattice + // maskThreshold - if positive, the value is treated as a threshold value to determine + // whether a pixel is good (mask value is greater than the threshold) or has to be + // masked (mask value is below the threshold). Negative threshold switches mask clipping + // off. The mask value is used to weight the flux during cleaning. This mode is used + // to implement cleaning based on the signal-to-noise as opposed to the standard cleaning + // based on the flux. The default threshold value is 0.9, which ensures the behavior of the + // code is exactly the same as before this parameter has been introduced. + void setMask(casacore::Matrix & mask, const casacore::Float& maskThreshold = 0.9); + + // remove the mask; + // useful when keeping object and sending a new dirty image to clean + // one can set another mask then + void unsetMask(); + + // Tell the algorithm to NOT clean just the inner quarter + // (This is useful when multiscale clean is being used + // inside a major cycle for MF or WF algorithms) + // if true, the full image deconvolution will be attempted + void ignoreCenterBox(casacore::Bool huh) { itsIgnoreCenterBox = huh; } + + // speedup() will speed the clean iteration by raising the + // threshold. This may be required if the threshold is + // accidentally set too low (ie, lower than can be achieved + // given errors in the approximate PSF). + // + // threshold(iteration) = threshold(0) + // * ( exp( (iteration - startingiteration)/Ndouble )/ 2.718 ) + // If speedup() is NOT invoked, no effect on threshold + void speedup(const casacore::Float Ndouble); + + // Look at what WE think the residuals look like + // Assumes the first scale is zero-sized + casacore::Matrix residual() { return itsDirty; } + //slightly better approximation of the residual: it convolves the given model + //with the psf and remove it from the dirty image put in setdirty + casacore::Matrix residual(const casacore::Matrix& model); + // Method to return threshold, including any speedup factors + casacore::Float threshold() const; + + // Method to return the strength optimum achieved at the last clean iteration + // The output of this method makes sense only if it is called after clean + casacore::Float strengthOptimum() const { return itsStrengthOptimum; } + + // Helper function to optimize adding + //static void addTo(casacore::Matrix& to, const casacore::Matrix & add); + + casacore::Matrix cmap; + casacore::Matrix mod; + casacore::Matrix delta; + + casacore::Matrix shifted; + casacore::Matrix clean_map; + + casacore::Matrix BB; + casacore::Matrix tildeII; + casacore::Matrix BI; + casacore::Matrix basis_function; + casacore::Matrix tildeMB; + casacore::Matrix mspsf; + casacore::Matrix tildeMI; + casacore::Matrix tildeMBB; + casacore::Matrix MtildeMBB; + casacore::Matrix tildeMBI; + casacore::Matrix MtildeMB; + casacore::Matrix window_basis; + + casacore::Float normcmap; + casacore::Float normcmapsq; + casacore::Sqrt sqrtop; + casacore::Float sumcmap; + casacore::Float triggerhogbom; + casacore::Float itsmaxbeam; + + casacore::IPosition blcautomask; + casacore::IPosition trcautomask; + + void setAutoControl(const casacore::Float autoThreshold = 0.0,const casacore::Int autoMaxiter = 0, const casacore::Float autoGain = 0.0, const casacore::Float hogbomGain = 0.0, const casacore::Bool autoHogbom = false, const casacore::Float autoTrigger = 1.0, const casacore::Float autoPower = 1.0, const casacore::Int autoXMask = 0, const casacore::Int autoYMask = 0) { itsAutoThreshold = autoThreshold; itsAutoMaxiter = autoMaxiter; itsAutoGain = autoGain; itsHogbomGain = hogbomGain; itsAutoHogbom = autoHogbom; itsAutoTrigger = autoTrigger; itsAutoPower = autoPower; itsAutoXMask = autoXMask; itsAutoYMask = autoYMask;} + +protected: + //friend class MatrixNACleaner; + // Make sure that the peak of the Psf is within the image + casacore::Bool validatePsf(const casacore::Matrix & psf); + + // Find the Peak of the matrix + static casacore::Bool findMaxAbs(const casacore::Matrix& lattice, + casacore::Float& maxAbs, casacore::IPosition& posMax); + + // This is made static since findMaxAbs is static(!). + // Why is findMaxAbs static??? + // --SB + static casacore::Bool findPSFMaxAbs(const casacore::Matrix& lattice, + casacore::Float& maxAbs, casacore::IPosition& posMax, + const casacore::Int& supportSize=100); + + casacore::Int findBeamPatch(const casacore::Float maxScaleSize, const casacore::Int& nx, const casacore::Int& ny, + const casacore::Float psfBeam=4.0, const casacore::Float nBeams=20.0); + // Find the Peak of the lattice, applying a mask + casacore::Bool findMaxAbsMask(const casacore::Matrix& lattice, const casacore::Matrix& mask, + casacore::Float& maxAbs, casacore::IPosition& posMax); + + // Helper function to reduce the box sizes until the have the same + // size keeping the centers intact + static void makeBoxesSameSize(casacore::IPosition& blc1, casacore::IPosition& trc1, + casacore::IPosition &blc2, casacore::IPosition& trc2); + + + casacore::CleanEnums::CleanType itsCleanType; + casacore::Float itsGain; + casacore::Int itsMaxNiter; // maximum possible number of iterations + casacore::Quantum itsThreshold; + casacore::CountedPtr > itsMask; + casacore::IPosition itsPositionPeakPsf; + casacore::Bool itsScalesValid; + casacore::Int itsNscales; + casacore::Float itsMaskThreshold; + + //# The following functions are used in various places in the code and are + //# documented in the .cc file. Static functions are used when the functions + //# do not modify the object state. They ensure that implicit assumptions + //# about the current state and implicit side-effects are not possible + //# because all information must be supplied in the input arguments + + + //casacore::CountedPtr > itsDirty; + casacore::Matrix itsDirty; + + //casacore::CountedPtr > itsScales; + casacore::Matrix itsScales; + casacore::Matrix itsPsf; + + //casacore::Matrix cmap; + //casacore::Matrix mod; + //casacore::Matrix delta; + + //casacore::Matrix shifted; + //casacore::Matrix clean_map; + + //casacore::Matrix BB; + //casacore::Matrix tildeII; + //casacore::Matrix BI; + //casacore::Matrix basis_function; + //casacore::Matrix tildeMB; + //casacore::Matrix mspsf; + //casacore::Matrix tildeMI; + //casacore::Matrix tildeMBB; + //casacore::Matrix MtildeMBB; + //casacore::Matrix tildeMBI; + //casacore::Matrix MtildeMB; + + casacore::Int itsIteration; // what iteration did we get to? + casacore::Int itsStartingIter; // what iteration did we get to? + casacore::Quantum itsFracThreshold; + + casacore::Float itsMaximumResidual; + casacore::Float itsStrengthOptimum; + + casacore::Vector itsTotalFluxScale; + casacore::Float itsTotalFlux; + + casacore::Bool itsIgnoreCenterBox; + + casacore::IPosition psfShape_p; + casacore::Bool noClean_p; + +private: + + // casacore::Memory to be allocated per TempLattice + casacore::Double itsMemoryMB; + + // Let the user choose whether to stop + casacore::Bool itsChoose; + + // Threshold speedup factors: + casacore::Bool itsDoSpeedup; // if false, threshold does not change with iteration + casacore::Float itsNDouble; + + casacore::Float itsAutoThreshold; + casacore::Int itsAutoMaxiter; + casacore::Float itsAutoGain; + casacore::Float itsHogbomGain; + casacore::Bool itsAutoHogbom; + casacore::Float itsAutoTrigger; + casacore::Float itsAutoPower; + casacore::Bool itsdimensionsareeven; + casacore::Int itsAutoXMask; + casacore::Int itsAutoYMask; + + //# Stop now? + //#// casacore::Bool stopnow(); Removed on 8-Apr-2004 by GvD + + // threshold for masks. If negative, mask values are used as weights and no pixels are + // discarded (although effectively they would be discarded if the mask value is 0.) +}; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/MeasurementEquations/MatrixCleaner.cc b/src/synthesis/MeasurementEquations/MatrixCleaner.cc index ebbd3788..9a9758fc 100644 --- a/src/synthesis/MeasurementEquations/MatrixCleaner.cc +++ b/src/synthesis/MeasurementEquations/MatrixCleaner.cc @@ -679,7 +679,7 @@ Int MatrixCleaner::clean(Matrix& model, Matrix modelSub=model(blc, trc); Matrix scaleSub=(itsScales[optimumScale])(blcPsf,trcPsf); - + // Now do the addition of this scale to the model image.... modelSub += scaleFactor*scaleSub; @@ -720,7 +720,7 @@ Int MatrixCleaner::clean(Matrix& model, if(!converged) { os << "Failed to reach stopping threshold" << LogIO::POST; } - + return converged; } diff --git a/src/synthesis/MeasurementEquations/SpectralAspCleaner.cc b/src/synthesis/MeasurementEquations/SpectralAspCleaner.cc new file mode 100644 index 00000000..a6ec5f1f --- /dev/null +++ b/src/synthesis/MeasurementEquations/SpectralAspCleaner.cc @@ -0,0 +1,235 @@ +//# Copyright (C) 1996-2010 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: casa-feedback@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# $Id: $AspMatrixCleaner.cc + +// Same include list as in MatrixCleaner.cc +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +// for alglib +#include +using namespace alglib; + +using namespace casacore; +using namespace std::chrono; +namespace casa { //# NAMESPACE CASA - BEGIN +SpectralAspCleaner::SpectralAspCleaner(): + AspMatrixCleaner() +{ +} + +SpectralAspCleaner::~SpectralAspCleaner() +{ + destroyAspScales(); + destroyInitMasks(); + destroyInitScales(); + if(!itsMask.null()) + itsMask=0; +} + +// ALGLIB +void SpectralAspCleaner::MFaspclean(Matrix& model) +{ + LogIO os(LogOrigin("SpectralAspCleaner", "MFaspclean", WHERE)); + + // the new aspen is always added to the active-set + casacore::Vector Optimums(itsMaxNiter); + casacore::Vector ScaleSizes(itsMaxNiter); + casacore::Vector Positions(itsMaxNiter); + casacore::Vector xPositions(itsMaxNiter); + casacore::Vector yPositions(itsMaxNiter); + + vector tempx; + vector activeSetCenter; + + //read in amplitudes and scales + read_array(Optimums, std::string("./strengthoptimum")); + read_array(ScaleSizes, std::string("./scalesizes")); + //read_array(Positions, std::string("./positions")); + read_array(xPositions, std::string("./xpositions")); + read_array(yPositions, std::string("./ypositions")); + + for (int ii=0; ii < itsMaxNiter; ii++) + { + IPosition pos(2,0); + pos(0) = xPositions(ii); + pos(1) = yPositions(ii); + Positions(ii) = pos; + } + + itsMaxNiter = Optimums.size(); + + //package information in the correct way + for (int ii=0; ii < itsMaxNiter; ii++){ + tempx.push_back(Optimums(ii)); + if (ScaleSizes(ii) > 0){ + tempx.push_back(ScaleSizes(ii));} + else{ + tempx.push_back(1.0);} + activeSetCenter.push_back(Positions(ii)); + } + + // initialize alglib option + unsigned int length = tempx.size(); + real_1d_array x; + x.setlength(length); + + real_1d_array s; + s.setlength(length); + + // initialize starting point + for (unsigned int i = 0; i < length; i+=2) + { + x[i] = tempx[i]; + x[i+1] = tempx[i+1]; + + s[i] = tempx[i]; + s[i+1] = tempx[i+1]; + } + + //real_1d_array s = "[1,1]"; + double epsg = itsLbfgsEpsG;//1e-3; + double epsf = itsLbfgsEpsF;//1e-3; + double epsx = itsLbfgsEpsX;//1e-3; + ae_int_t maxits = itsLbfgsMaxit;//5; + minlbfgsstate state; + minlbfgscreate(1, x, state); + minlbfgssetcond(state, epsg, epsf, epsx, maxits); + minlbfgssetscale(state, s); + minlbfgssetprecscale(state); + minlbfgsreport rep; + + ParamAlglibObj optParam(*itsDirty, *itsXfr, activeSetCenter, fft); + ParamAlglibObj *ptrParam; + ptrParam = &optParam; + + alglib::minlbfgsoptimize(state, objfunc_alglib, NULL, (void *) ptrParam); + + minlbfgsresults(state, x, rep); + + // end alglib bfgs optimization + for (int ii=0; ii < itsMaxNiter; ii++){ + double amp = x[2*ii]; + double scale = x[2*ii+1]; + itsGoodAspAmplitude.push_back(amp); + itsGoodAspActiveSet.push_back(scale); + itsGoodAspCenter.push_back(activeSetCenter[ii]); + } + + // now perform the CLEANing + // Initialization + Float totalFlux=0.0; + // Define a subregion so that the peak is centered + IPosition support(model.shape()); + support(0) = max(Int(itsInitScaleSizes[itsNInitScales-1] + 0.5), support(0)); + support(1) = max(Int(itsInitScaleSizes[itsNInitScales-1] + 0.5), support(1)); + IPosition inc(model.shape().nelements(), 1); + Matrix itsScale = Matrix(psfShape_p); + MatrixitsScaleXfr = Matrix (); + Float maxVal=0; + IPosition posmin((*itsDirty).shape().nelements(), 0); + Float minVal=0; + IPosition posmax((*itsDirty).shape().nelements(), 0); + + //Loop over Asps + for(int ii=0; ii < itsMaxNiter; ii++){ + //Load the Asp parameters + itsOptimumScaleSize = itsGoodAspActiveSet[ii]; + itsStrengthOptimum = itsGoodAspAmplitude[ii]; + itsPositionOptimum = itsGoodAspCenter[ii]; + + //Gain is one, since we substitute the full Asp + itsGain = 1.0; + + //Make the scale + makeScaleImage(itsScale, itsOptimumScaleSize, itsStrengthOptimum, itsPositionOptimum); + itsScaleXfr.resize(); + fft.fft0(itsScaleXfr, itsScale); + + // Now add to the total flux + totalFlux += (itsStrengthOptimum*itsGain); + itsTotalFlux = totalFlux; + + //The slicing box + IPosition blc(itsPositionOptimum - support/2); + IPosition trc(itsPositionOptimum + support/2 - 1); + + LCBox::verify(blc, trc, inc, model.shape()); + IPosition blcPsf(blc); + IPosition trcPsf(trc); + + // Update the model image + Matrix modelSub = model(blc, trc); + Float scaleFactor; + scaleFactor = itsGain * itsStrengthOptimum; + Matrix scaleSub = (itsScale)(blcPsf,trcPsf); + modelSub += scaleFactor * scaleSub; + + // Now update the residual image + // PSF * model + Matrix cWork; + cWork = ((*itsXfr)*(itsScaleXfr)); //Asp's + Matrix itsPsfConvScale = Matrix(psfShape_p); + fft.fft0(itsPsfConvScale, cWork, false); + fft.flip(itsPsfConvScale, false, false); //need this if conv with 1 scale; don't need this if conv with 2 scales + + IPosition nullnull(2,0); + Matrix shift(psfShape_p); + shift.assign_conforming(itsPsfConvScale); + Matrix sub = itsPsfConvScale(nullnull+2,support-1); + sub.assign_conforming(shift(nullnull,support-3)); + + Matrix psfSub = (itsPsfConvScale)(blcPsf, trcPsf); + Matrix dirtySub=(*itsDirty)(blc,trc); + + os << "itsStrengthOptimum " << itsStrengthOptimum << LogIO::POST; + + // subtract the peak that we found from the dirty image + dirtySub -= scaleFactor * psfSub; + + // update peakres + itsPrevPeakResidual = itsPeakResidual; + maxVal=0; + posmin = 0; + minVal=0; + posmax = 0; + minMaxMasked(minVal, maxVal, posmin, posmax, (*itsDirty), itsInitScaleMasks[0]); + itsPeakResidual = (fabs(maxVal) > fabs(minVal)) ? fabs(maxVal) : fabs(minVal); + os << "current peakres " << itsPeakResidual << LogIO::POST; + } + os << "Aspclean finished" << LogIO::POST; +} + +} //# NAMESPACE CASA - END diff --git a/src/synthesis/MeasurementEquations/SpectralAspCleaner.h b/src/synthesis/MeasurementEquations/SpectralAspCleaner.h new file mode 100644 index 00000000..64fb7455 --- /dev/null +++ b/src/synthesis/MeasurementEquations/SpectralAspCleaner.h @@ -0,0 +1,56 @@ +//# AspMatrixCleaner.h: Minor Cycle for Asp deconvolution +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id: AspMatrixCleaner.h Genie H. 2020-04-06 +#include +#include +#include +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + +class SpectralAspCleaner : public AspMatrixCleaner +{ +public: + // Create a cleaner : default constructor + SpectralAspCleaner(); + + // The destructor does nothing special. + ~SpectralAspCleaner(); + + void MFaspclean(casacore::Matrix & model); +}; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/MeasurementEquations/WaveletAspCleaner.cc b/src/synthesis/MeasurementEquations/WaveletAspCleaner.cc new file mode 100644 index 00000000..15b3ea17 --- /dev/null +++ b/src/synthesis/MeasurementEquations/WaveletAspCleaner.cc @@ -0,0 +1,248 @@ +//# Copyright (C) 1996-2010 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: casa-feedback@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# $Id: $AspMatrixCleaner.cc + +// Same include list as in MatrixCleaner.cc +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +// for alglib +#include +using namespace alglib; + +using namespace casacore; +using namespace std::chrono; +namespace casa { //# NAMESPACE CASA - BEGIN +WaveletAspCleaner::WaveletAspCleaner(): + AspMatrixCleaner() +{ + itsWaveletScales.resize(0); + itsWaveletAmps.resize(0); +} + +WaveletAspCleaner::~WaveletAspCleaner() +{ + destroyAspScales(); + destroyInitMasks(); + destroyInitScales(); + if(!itsMask.null()) + itsMask=0; +} + +float WaveletAspCleaner::aspScaleModel(const Int& i, const Int& j, const Float& scaleSize, const Int& refi, const Int& refj) +{ + return (1.0/(2*M_PI*pow(scaleSize,2)))*exp(-(pow(i-refi,2) + pow(j-refj,2))*0.5/pow(scaleSize,2)); +} + +float WaveletAspCleaner::aspPeakNormModel(const Float& width) +{ + return 2 * M_PI / pow(width,2) ; +} + +float WaveletAspCleaner::aspNormalizationModel(const Float& width1, const Float& width2) +{ + return 2 * M_PI / (pow(1.0/width1, 2) + pow(1.0/width2, 2)); +} + +// ALGLIB - gold - not "log" + +vector WaveletAspCleaner::getActiveSetAspen() +{ + LogIO os(LogOrigin("AspMatrixCleaner", "getActiveSetAspen()", WHERE)); + + if(int(itsInitScaleXfrs.nelements()) == 0) + throw(AipsError("Initial scales for Asp are not defined")); + + if (!itsSwitchedToHogbom && + accumulate(itsNumIterNoGoodAspen.begin(), itsNumIterNoGoodAspen.end(), 0) >= 5) + { + os << "Switched to hogbom because of frequent small components." << LogIO::POST; + switchedToHogbom(); + } + + if (itsSwitchedToHogbom) + itsNInitScales = 1; + else + itsNInitScales = itsInitScaleSizes.size(); + + // Dirty * initial scales + Matrix dirtyFT; + fft.fft0(dirtyFT, *itsDirty); + itsDirtyConvInitScales.resize(0); + itsDirtyConvInitScales.resize(itsNInitScales); // 0, 1width, 2width, 4width and 8width + + //cout << "itsInitScaleSizes.size() " << itsInitScaleSizes.size() << " itsInitScales.size() " << itsInitScales.size() << " NInitScales # " << itsNInitScales << endl; + for (int scale=0; scale < itsNInitScales; scale++) + { + Matrix cWork; + + itsDirtyConvInitScales[scale] = Matrix(itsDirty->shape()); + cWork=((dirtyFT)*(itsInitScaleXfrs[scale])); + fft.fft0((itsDirtyConvInitScales[scale]), cWork, false); + fft.flip((itsDirtyConvInitScales[scale]), false, false); + + //cout << "remake itsDirtyConvInitScales " << scale << " max itsInitScales[" << scale << "] = " << max(fabs(itsInitScales[scale])) << endl; + //cout << " max itsInitScaleXfrs[" << scale << "] = " << max(fabs(itsInitScaleXfrs[scale])) << endl; + } + + float strengthOptimum = 0.0; + int optimumScale = 0; + IPosition positionOptimum(itsDirty->shape().nelements(), 0); + itsGoodAspActiveSet.resize(0); + itsGoodAspAmplitude.resize(0); + itsGoodAspCenter.resize(0); + + maxDirtyConvInitScales(strengthOptimum, optimumScale, positionOptimum); + + os << "Peak among the smoothed residual image is " << strengthOptimum << " and initial scale: " << optimumScale << LogIO::POST; + //cout << "Peak among the smoothed residual image is " << strengthOptimum << " and initial scale: " << optimumScale << endl; + // cout << " its itsDirty is " << (*itsDirty)(positionOptimum); + // cout << " at location " << positionOptimum[0] << " " << positionOptimum[1] << " " << positionOptimum[2]; + + + itsStrengthOptimum = strengthOptimum; + itsPositionOptimum = positionOptimum; + itsOptimumScale = optimumScale; + itsOptimumScaleSize = itsInitScaleSizes[optimumScale]; + + // initial scale size = 0 gives the peak res, so we don't + // need to do the LBFGS optimization for it + if (itsOptimumScale == 0) + return {}; + else + { + // the new aspen is always added to the active-set + vector tempx; + vector activeSetCenter; + + tempx.push_back(strengthOptimum); + tempx.push_back(itsInitScaleSizes[optimumScale]); + activeSetCenter.push_back(positionOptimum); + + // initialize alglib option + unsigned int length = tempx.size(); + real_1d_array x; + x.setlength(length); + + // for G55 ,etc + real_1d_array s; + s.setlength(length); + + // initialize starting point + for (unsigned int i = 0; i < length; i+=2) + { + x[i] = tempx[i]; //amp + x[i+1] = tempx[i+1]; //scale + + s[i] = tempx[i]; //amp + s[i+1] = tempx[i+1]; //scale + } + + + //real_1d_array s = "[1,1]"; + double epsg = itsLbfgsEpsG;//1e-3; + double epsf = itsLbfgsEpsF;//1e-3; + double epsx = itsLbfgsEpsX;//1e-3; + ae_int_t maxits = itsLbfgsMaxit;//5; + //double epsg = 1e-3; + //double epsf = 1e-3; + //double epsx = 1e-3; + //ae_int_t maxits = 5; + minlbfgsstate state; + minlbfgscreate(1, x, state); + minlbfgssetcond(state, epsg, epsf, epsx, maxits); + minlbfgssetscale(state, s); + + //minlbfgssetprecscale(state); + minlbfgsreport rep; + + ParamAlglibObjWavelet optParam(*itsDirty, activeSetCenter, itsWaveletScales, itsWaveletAmps); + ParamAlglibObjWavelet *ptrParam; + ptrParam = &optParam; + + alglib::minlbfgsoptimize(state, objfunc_alglib_wavelet, NULL, (void *) ptrParam); + + minlbfgsresults(state, x, rep); + //double *x1 = x.getcontent(); + //cout << "x1[0] " << x1[0] << " x1[1] " << x1[1] << endl; + + // end alglib bfgs optimization + + + + + + + double amp = x[0]; // i + double scale = x[1]; // i+1 + + if (fabs(scale) < 0.4) + { + scale = 0; + amp = (*itsDirty)(itsPositionOptimum); // This is to avoid divergence due to amp being too large. + // amp=strengthOptimum gives similar results + } + else + scale = fabs(scale); + + itsGoodAspAmplitude.push_back(amp); // active-set amplitude + itsGoodAspActiveSet.push_back(scale); // active-set + + itsStrengthOptimum = amp; + itsOptimumScaleSize = scale; + itsGoodAspCenter = activeSetCenter; + + // debug + //os << "optimized strengthOptimum " << itsStrengthOptimum << " scale size " << itsOptimumScaleSize << LogIO::POST; + //cout << "optimized strengthOptimum " << itsStrengthOptimum << " scale size " << itsOptimumScaleSize << " at " << itsPositionOptimum << endl; + + } // finish bfgs optimization + + AlwaysAssert(itsGoodAspCenter.size() == itsGoodAspActiveSet.size(), AipsError); + AlwaysAssert(itsGoodAspAmplitude.size() == itsGoodAspActiveSet.size(), AipsError); + + // debug info + /*for (unsigned int i = 0; i < itsAspAmplitude.size(); i++) + { + //cout << "After opt AspApm[" << i << "] = " << itsAspAmplitude[i] << endl; + //cout << "After opt AspScale[" << i << "] = " << itsAspScaleSizes[i] << endl; + //cout << "After opt AspCenter[" << i << "] = " << itsAspCenter[i] << endl; + cout << "AspScale[ " << i << " ] = " << itsAspScaleSizes[i] << " center " << itsAspCenter[i] << endl; + }*/ + + return itsGoodAspActiveSet; // return optimized scale +} + +} //# NAMESPACE CASA - END diff --git a/src/synthesis/MeasurementEquations/WaveletAspCleaner.h b/src/synthesis/MeasurementEquations/WaveletAspCleaner.h new file mode 100644 index 00000000..74baa903 --- /dev/null +++ b/src/synthesis/MeasurementEquations/WaveletAspCleaner.h @@ -0,0 +1,70 @@ +//# AspMatrixCleaner.h: Minor Cycle for Asp deconvolution +//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003 +//# Associated Universities, Inc. Washington DC, USA. +//# +//# This library is free software; you can redistribute it and/or modify it +//# under the terms of the GNU Library General Public License as published by +//# the Free Software Foundation; either version 2 of the License, or (at your +//# option) any later version. +//# +//# This library is distributed in the hope that it will be useful, but WITHOUT +//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public +//# License for more details. +//# +//# You should have received a copy of the GNU Library General Public License +//# along with this library; if not, write to the Free Software Foundation, +//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. +//# +//# Correspondence concerning AIPS++ should be addressed as follows: +//# Internet email: aips2-request@nrao.edu. +//# Postal address: AIPS++ Project Office +//# National Radio Astronomy Observatory +//# 520 Edgemont Road +//# Charlottesville, VA 22903-2475 USA +//# +//# +//# $Id: AspMatrixCleaner.h Genie H. 2020-04-06 +#include +#include +#include +#include +#include + +namespace casa { //# NAMESPACE CASA - BEGIN + +class WaveletAspCleaner : public AspMatrixCleaner +{ +public: + // Create a cleaner : default constructor + WaveletAspCleaner(); + + // The destructor does nothing special. + ~WaveletAspCleaner(); + + // Make an image of the specified scale by Gaussian + std::vector getActiveSetAspen() override; + + //Normalization Helpers + float aspScaleModel(const casacore::Int& i, const casacore::Int& j, const casacore::Float& scaleSize, const casacore::Int& refi, const casacore::Int& refj) override; + float aspPeakNormModel(const casacore::Float& width) override; + float aspNormalizationModel(const casacore::Float& width1, const casacore::Float& width2) override; + + void setWaveletControl(const casacore::Vector waveletScales, const casacore::Vector waveletAmps) { itsWaveletScales = waveletScales; itsWaveletAmps=waveletAmps;} + +protected: +//private: + + casacore::Vector itsWaveletScales; + casacore::Vector itsWaveletAmps; +}; + +} //# NAMESPACE CASA - END + +#endif diff --git a/src/synthesis/MeasurementEquations/objfunc_alglib_wavelets.h b/src/synthesis/MeasurementEquations/objfunc_alglib_wavelets.h new file mode 100644 index 00000000..aef5284a --- /dev/null +++ b/src/synthesis/MeasurementEquations/objfunc_alglib_wavelets.h @@ -0,0 +1,185 @@ +#ifndef SYNTHESIS_OBJFUNCALGLIB_WAVELETS_H +#define SYNTHESIS_OBJFUNCALGLIB_WAVELETS_H + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "lbfgs/optimization.h" + +#ifndef isnan +#define isnan(x) std::isnan(x) +#endif + +namespace casa { //# NAMESPACE CASA - BEGIN + +class ParamAlglibObjWavelet +{ +private: + int nX; + int nY; + unsigned int AspLen; + casacore::Matrix itsMatDirty; + std::vector center; + casacore::Matrix newResidual; + casacore::Matrix AspConvPsf; + casacore::Matrix dAspConvPsf; + casacore::Vector itsWaveletScales; + casacore::Vector itsWaveletAmps; + +public: + ParamAlglibObjWavelet(const casacore::Matrix& dirty, + const std::vector& positionOptimum, const casacore::Vector& waveletScales, const casacore::Vector& waveletAmps): + itsMatDirty(dirty), itsWaveletScales(waveletScales), itsWaveletAmps(waveletAmps), + center(positionOptimum) + { + nX = itsMatDirty.shape()(0); + nY = itsMatDirty.shape()(1); + AspLen = center.size(); + newResidual.resize(nX, nY); + AspConvPsf.resize(nX, nY); + dAspConvPsf.resize(nX, nY); + } + + ~ParamAlglibObjWavelet() = default; + + casacore::Matrix getterDirty() { return itsMatDirty; } + std::vector getterCenter() {return center;} + unsigned int getterAspLen() { return AspLen; } + int getterNX() { return nX; } + int getterNY() { return nY; } + casacore::Matrix getterRes() { return newResidual; } + void setterRes(const casacore::Matrix& res) { newResidual = res; } + casacore::Matrix getterAspConvPsf() { return AspConvPsf; } + void setterAspConvPsf(const casacore::Matrix& m) { AspConvPsf = m; } + casacore::Matrix getterDAspConvPsf() { return dAspConvPsf; } + casacore::Vector getterScales() { return itsWaveletScales; } + casacore::Vector getterAmps() { return itsWaveletAmps; } +}; + +inline void objfunc_alglib_wavelet(const alglib::real_1d_array &x, double &func, alglib::real_1d_array &grad, void *ptr) +{ + // retrieve params for GSL bfgs optimization + casa::ParamAlglibObjWavelet *MyP = (casa::ParamAlglibObjWavelet *) ptr; //re-cast back to ParamAlglibObj to retrieve images + + casacore::Matrix itsMatDirty(MyP->getterDirty()); + std::vector center = MyP->getterCenter(); + const unsigned int AspLen = MyP->getterAspLen(); + const int nX = MyP->getterNX(); + const int nY = MyP->getterNY(); + casacore::Matrix newResidual(MyP->getterRes()); + casacore::Matrix AspConvPsf(MyP->getterAspConvPsf()); + casacore::Matrix dAspConvPsf(MyP->getterDAspConvPsf()); + + casacore::Vector itsWaveletScales(MyP->getterScales()); + casacore::Vector itsWaveletAmps(MyP->getterAmps()); + + int nscales; + itsWaveletScales.shape(nscales); + + func = 0; + double amp = 1; + + const int refi = nX/2; + const int refj = nY/2; + + int minX = nX - 1; + int maxX = 0; + int minY = nY - 1; + int maxY = 0; + + // First, get the amp * AspenConvPsf for each Aspen to update the residual + for (unsigned int k = 0; k < AspLen; k ++) + { + amp = x[2*k]; + double scale = x[2*k+1]; + //std::cout << "f: amp " << amp << " scale " << scale << std::endl; + + if (isnan(amp) || scale < 0.4) // GSL scale < 0 + { + //std::cout << "nan? " << amp << " neg scale? " << scale << std::endl; + // If scale is small (<0.4), make it 0 scale to utilize Hogbom and save time + scale = (scale = fabs(scale)) < 0.4 ? 0 : scale; + //std::cout << "reset neg scale to " << scale << std::endl; + + if (scale <= 0) + return; + } + + AspConvPsf = 0.0; + dAspConvPsf = 0.0; + + for (int l = 0; l < nscales; l++){ + float temp_scale = sqrt(scale*scale+itsWaveletScales(l)*itsWaveletScales(l)); + float toadd; + + const double sigma5 = 5 * temp_scale / 2; + const int minI = std::max(0, (int)(center[k][0] - sigma5)); + const int maxI = std::min(nX-1, (int)(center[k][0] + sigma5)); + const int minJ = std::max(0, (int)(center[k][1] - sigma5)); + const int maxJ = std::min(nY-1, (int)(center[k][1] + sigma5)); + + if (minI < minX) + minX = minI; + if (maxI > maxX) + maxX = maxI; + if (minJ < minY) + minY = minJ; + if (maxJ > maxY) + maxY = maxJ; + + for (int j = minJ; j <= maxJ; j++) + { + for (int i = minI; i <= maxI; i++) + { + const int px = i; + const int py = j; + + toadd = itsWaveletAmps(l) * 1.0/(2*M_PI*pow(temp_scale,2))*exp(-(pow(i-center[k][0],2) + pow(j-center[k][1],2))*0.5/pow(temp_scale,2)); + //toadd = itsWaveletAmps(l) * 1.0/(sqrt(2*M_PI)*temp_scale)*exp(-(pow(i-center[k][0],2) + pow(j-center[k][1],2))*0.5/pow(temp_scale,2)); + AspConvPsf(i,j) += toadd; + dAspConvPsf(i,j)+= toadd * ((pow(i-center[k][0],2) + pow(j-center[k][1],2)) / pow(temp_scale,2) - 2) * fabs(scale) / pow(temp_scale, 2); + //dAspConvPsf(i,j)+= toadd * ((pow(i-center[k][0],2) + pow(j-center[k][1],2)) / pow(temp_scale,2) - 1) * fabs(scale) / pow(temp_scale, 2); + } + } + } + + } // end get AspenConvPsf + + // reset grad to 0. This is important to get the correct optimization. + double dA = 0.0; + double dS = 0.0; + + // Update the residual using the current residual image and the latest Aspen. + // Sanjay used, Res = OrigDirty - active-set aspen * Psf, in 2004, instead. + // Both works but the current approach is simpler and performs well too. + for (int j = minY; j < maxY; ++j) + { + for(int i = minX; i < maxX; ++i) + { + newResidual(i, j) = itsMatDirty(i, j) - amp * AspConvPsf(i, j); + func = func + double(pow(newResidual(i, j), 2)); + + // derivatives of amplitude + dA += double((-2) * newResidual(i,j) * AspConvPsf(i,j)); + // derivative of scale + dS += double((-2) * amp * newResidual(i,j) * dAspConvPsf(i,j)); + } + } + //std::cout << "after f " << func << std::endl; + + grad[0] = dA; + grad[1] = dS; +} + + + +} // end namespace casa + +#endif // SYNTHESIS_OBJFUNCALGLIB_WAVELETS_H