LCOV - code coverage report
Current view: top level - synthesis/ImagerObjects - SynthesisDeconvolver.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 521 639 81.5 %
Date: 2024-12-11 20:54:31 Functions: 27 31 87.1 %

          Line data    Source code
       1             : //# SynthesisDeconvolver.cc: Implementation of Imager.h
       2             : //# Copyright (C) 1997-2020
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This program is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU General Public License as published by the Free
       7             : //# Software Foundation; either version 2 of the License, or (at your option)
       8             : //# any later version.
       9             : //#
      10             : //# This program is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
      13             : //# more details.
      14             : //#
      15             : //# You should have received a copy of the GNU General Public License along
      16             : //# with this program; if not, write to the Free Software Foundation, Inc.,
      17             : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: casa-feedback@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id$
      27             : 
      28             : #include <casacore/casa/Exceptions/Error.h>
      29             : #include <iostream>
      30             : #include <sstream>
      31             : 
      32             : #include <casacore/casa/Arrays/Matrix.h>
      33             : #include <casacore/casa/Arrays/ArrayMath.h>
      34             : #include <casacore/casa/Arrays/ArrayLogical.h>
      35             : 
      36             : #include <casacore/casa/Logging.h>
      37             : #include <casacore/casa/Logging/LogIO.h>
      38             : #include <casacore/casa/Logging/LogMessage.h>
      39             : #include <casacore/casa/Logging/LogSink.h>
      40             : #include <casacore/casa/Logging/LogMessage.h>
      41             : 
      42             : #include <casacore/casa/OS/DirectoryIterator.h>
      43             : #include <casacore/casa/OS/File.h>
      44             : #include <casacore/casa/OS/Path.h>
      45             : 
      46             : #include <casacore/casa/OS/HostInfo.h>
      47             : 
      48             : #include <casacore/images/Images/TempImage.h>
      49             : #include <casacore/images/Images/SubImage.h>
      50             : #include <casacore/images/Regions/ImageRegion.h>
      51             : #include <casacore/lattices/Lattices/LatticeLocker.h>
      52             : #include <synthesis/ImagerObjects/CubeMinorCycleAlgorithm.h>
      53             : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
      54             : #include <synthesis/ImagerObjects/SynthesisDeconvolver.h>
      55             : #include <synthesis/ImagerObjects/SIMinorCycleController.h>
      56             : 
      57             : #include <sys/types.h>
      58             : #include <unistd.h>
      59             : using namespace std;
      60             : 
      61             : using namespace casacore;
      62             : extern casa::Applicator casa::applicator;
      63             : namespace casa { //# NAMESPACE CASA - BEGIN
      64             : 
      65        1063 :   SynthesisDeconvolver::SynthesisDeconvolver() :
      66        1063 :                                        itsDeconvolver(nullptr),
      67        1063 :                                        itsMaskHandler(nullptr ),
      68        1063 :                itsImages(nullptr),
      69        1063 :                itsImageName(""),
      70             :                                        //itsPartImageNames(Vector<String>(0)),
      71        1063 :                                        itsBeam(0.0),
      72        1063 :                                        itsDeconvolverId(0),
      73        1063 :                                        itsScales(Vector<Float>()),
      74        1063 :                itsMaskType(""),
      75        1063 :                itsPBMask(0.0),
      76             :                                        //itsMaskString(String("")),
      77        1063 :                itsIterDone(0.0),
      78        1063 :                itsChanFlag(Vector<Bool>(0)),
      79        1063 :                itsRobustStats(Record()),
      80        1063 :                initializeChanMaskFlag(false),
      81        1063 :                itsPosMask(nullptr),
      82        1063 :                                        itsIsMaskLoaded(false),
      83        1063 :                                        itsMaskSum(-1e+9),
      84        1063 :                                        itsPreviousFutureRes(0.0),
      85        5315 :                                        itsPreviousIterBotRec_p(Record())
      86             :   {
      87        1063 :   }
      88             : 
      89        1063 :   SynthesisDeconvolver::~SynthesisDeconvolver()
      90             :   {
      91        2126 :         LogIO os( LogOrigin("SynthesisDeconvolver","descructor",WHERE) );
      92        1063 :         os << LogIO::DEBUG1 << "SynthesisDeconvolver destroyed" << LogIO::POST;
      93        1063 :         SynthesisUtilMethods::getResource("End SynthesisDeconvolver");
      94             : 
      95        1063 :   }
      96             : 
      97        1063 :   void SynthesisDeconvolver::setupDeconvolution(const SynthesisParamsDeconv& decpars)
      98             :   {
      99        2126 :     LogIO os( LogOrigin("SynthesisDeconvolver","setupDeconvolution",WHERE) );
     100             : 
     101             :     //Copy this decpars into a private variable that can be used elsewhere
     102             :     //there is no proper copy operator (as public casa::Arrays members = operator fails)
     103        1063 :     itsDecPars.fromRecord(decpars.toRecord());
     104        1063 :     itsImageName = decpars.imageName;
     105        1063 :     itsStartingModelNames = decpars.startModel;
     106        1063 :     itsDeconvolverId = decpars.deconvolverId;
     107             : 
     108        1063 :     os << "Set Deconvolution Options for [" << itsImageName << "] : " << decpars.algorithm ;
     109             : 
     110        1063 :     if( itsStartingModelNames.nelements()>0 && itsStartingModelNames[0].length() > 0 )
     111          26 :       os << " , starting from model : " << itsStartingModelNames;
     112        1063 :     os << LogIO::POST;
     113             : 
     114             :     try
     115             :       {
     116        1063 :         if(decpars.algorithm==String("hogbom"))
     117             :           {
     118         800 :             itsDeconvolver.reset(new SDAlgorithmHogbomClean());
     119             :           }
     120         263 :         else if(decpars.algorithm==String("mtmfs"))
     121             :           {
     122         148 :             itsDeconvolver.reset(new SDAlgorithmMSMFS( decpars.nTaylorTerms, decpars.scales, decpars.scalebias ));
     123             :           }
     124         115 :         else if(decpars.algorithm==String("clark_exp"))
     125             :           {
     126           0 :             itsDeconvolver.reset(new SDAlgorithmClarkClean("clark"));
     127             :           }
     128         115 :         else if(decpars.algorithm==String("clarkstokes_exp"))
     129             :           {
     130           0 :             itsDeconvolver.reset(new SDAlgorithmClarkClean("clarkstokes"));
     131             :           }
     132         115 :         else if(decpars.algorithm==String("clark"))
     133             :           {
     134          64 :             itsDeconvolver.reset(new SDAlgorithmClarkClean2("clark"));
     135             :           }
     136          51 :         else if(decpars.algorithm==String("clarkstokes"))
     137             :           {
     138           5 :             itsDeconvolver.reset(new SDAlgorithmClarkClean2("clarkstokes"));
     139             :           }
     140          46 :         else if(decpars.algorithm==String("multiscale"))
     141             :           {
     142          34 :             itsDeconvolver.reset(new SDAlgorithmMSClean( decpars.scales, decpars.scalebias ));
     143             :           }
     144          12 :         else if(decpars.algorithm==String("mem"))
     145             :           {
     146           1 :             itsDeconvolver.reset(new SDAlgorithmMEM( "entropy" ));
     147             :           }
     148          11 :         else if (decpars.algorithm==String("asp"))
     149             :           {
     150          11 :       bool isSingle = false;
     151          11 :       if (decpars.specmode == String("mfs"))
     152           0 :         isSingle = true;
     153             : 
     154          11 :             itsDeconvolver.reset(new SDAlgorithmAAspClean(decpars.fusedThreshold, isSingle, decpars.largestscale));
     155             :           }
     156             :         else
     157             :           {
     158           0 :             throw( AipsError("Un-known algorithm : "+decpars.algorithm) );
     159             :           }
     160             : 
     161             :         // Set restoring beam options
     162        1063 :         itsDeconvolver->setRestoringBeam( decpars.restoringbeam, decpars.usebeam );
     163        1063 :         itsUseBeam = decpars.usebeam;// store this info also here.
     164             : 
     165             :         // Set Masking options
     166             :         //      itsDeconvolver->setMaskOptions( decpars.maskType );
     167        1063 :         itsMaskHandler.reset(new SDMaskHandler());
     168             :         //itsMaskString = decpars.maskString;
     169        1063 :         itsMaskType = decpars.maskType;
     170        1063 :         if(itsMaskType=="auto-thresh")
     171             :           {
     172           0 :             itsAutoMaskAlgorithm="thresh";
     173             :           }
     174        1063 :         else if(itsMaskType=="auto-thresh2")
     175             :           {
     176           0 :             itsAutoMaskAlgorithm="thresh2";
     177             :           }
     178        1063 :         else if(itsMaskType=="auto-multithresh")
     179             :           {
     180          46 :             itsAutoMaskAlgorithm="multithresh";
     181             :           }
     182        1017 :         else if(itsMaskType=="auto-onebox")
     183             :           {
     184           0 :             itsAutoMaskAlgorithm="onebox";
     185             :           }
     186             :         else {
     187        1017 :             itsAutoMaskAlgorithm="";
     188             :         }
     189        1063 :         itsPBMask = decpars.pbMask;
     190        1063 :         itsMaskString = decpars.maskString;
     191        2663 :         if(decpars.maskList.nelements()==0 ||
     192        1600 :             (decpars.maskList.nelements()==1 && decpars.maskList[0]==""))
     193             :           {
     194        1063 :             itsMaskList.resize(1);
     195        1063 :             itsMaskList[0] = itsMaskString;
     196             :           }
     197             :         else
     198             :           {
     199           0 :             itsMaskList = decpars.maskList;
     200             :           }
     201        1063 :         itsMaskThreshold = decpars.maskThreshold;
     202        1063 :         itsFracOfPeak = decpars.fracOfPeak;
     203        1063 :         itsMaskResolution = decpars.maskResolution;
     204        1063 :         itsMaskResByBeam = decpars.maskResByBeam;
     205        1063 :         itsNMask = decpars.nMask;
     206             :         //itsAutoAdjust = decpars.autoAdjust;
     207             :         //desable autoadjust
     208        1063 :         itsAutoAdjust = false;
     209        1063 :         itsSidelobeThreshold = decpars.sidelobeThreshold;
     210        1063 :         itsNoiseThreshold = decpars.noiseThreshold;
     211        1063 :         itsLowNoiseThreshold = decpars.lowNoiseThreshold;
     212        1063 :         itsNegativeThreshold = decpars.negativeThreshold;
     213        1063 :         itsSmoothFactor = decpars.smoothFactor;
     214        1063 :         itsMinBeamFrac = decpars.minBeamFrac;
     215        1063 :         itsCutThreshold = decpars.cutThreshold;
     216        1063 :         itsGrowIterations = decpars.growIterations;
     217        1063 :         itsDoGrowPrune = decpars.doGrowPrune;
     218        1063 :         itsMinPercentChange = decpars.minPercentChange;
     219        1063 :         itsVerbose = decpars.verbose;
     220        1063 :         itsFastNoise = decpars.fastnoise;
     221        1063 :               itsIsInteractive = decpars.interactive;
     222        1063 :         itsNsigma = decpars.nsigma;
     223        1063 :         itsNoRequireSumwt = true; //decpars.noRequireSumwt;
     224        1063 :         itsFullSummary = decpars.fullsummary;
     225             :       }
     226           0 :     catch(AipsError &x)
     227             :       {
     228           0 :         throw( AipsError("Error in constructing a Deconvolver : "+x.getMesg()) );
     229           0 :       }
     230             : 
     231        1063 :     itsAddedModel=false;
     232        1063 :   }
     233             : 
     234         771 :   Long SynthesisDeconvolver::estimateRAM(const vector<int>& imsize){
     235             : 
     236         771 :     Long mem=0;
     237             :     /* This does not work
     238             :     if( ! itsImages )
     239             :       {
     240             :         itsImages = makeImageStore( itsImageName );
     241             :       }
     242             :     */
     243         771 :     if(itsDeconvolver)
     244         771 :       mem=itsDeconvolver->estimateRAM(imsize);
     245         771 :     return mem;
     246             :   }
     247             : 
     248        2205 :   Record SynthesisDeconvolver::initMinorCycle() {
     249             :     /////IMPORTANT initMinorCycle has to be called before setupMask...that order has to be kept !
     250             : 
     251        2205 :     if(!itsImages)
     252         219 :       itsImages=makeImageStore(itsImageName, itsNoRequireSumwt);
     253             :     //For cubes as we are not doing a post major cycle residual automasking
     254             :     //Force recalculation of robust stats to update nsigmathreshold with
     255             :     //most recent residual
     256             : 
     257        2205 :     if(itsAutoMaskAlgorithm=="multithresh" && itsImages->residual()->shape()[3] >1 && itsNsigma > 0.0){
     258           0 :       Record retval;
     259           0 :       Record backupRobustStats=itsRobustStats;
     260           0 :       itsRobustStats=Record();
     261           0 :       retval=initMinorCycle(itsImages);
     262           0 :       itsRobustStats=backupRobustStats;
     263           0 :       return retval;
     264           0 :     }
     265             : 
     266             :     /* else if (itsAutoMaskAlgorithm=="multithresh" && itsImages->residual()->shape()[3]){
     267             :       ///As the automask for cubes pre-CAS-9386...
     268             :       /// was tuned to look for threshold in future mask
     269             :       ///It is as good as somewhere in between no mask and mask
     270             :       //      Record backupRobustStats=itsRobustStats;
     271             :       Record retval=initMinorCycle(itsImages);
     272             :       //cerr << "INITMINOR " << itsRobustStats << endl;
     273             :       //itsRobustStats=backupRobustStats;
     274             :       if(retval.isDefined("peakresidualnomask")){
     275             :         Float futureRes=Float(retval.asFloat("peakresidualnomask")-(retval.asFloat("peakresidualnomask")-retval.asFloat("peakresidual"))/1000.0);
     276             :         if(futureRes != itsPreviousFutureRes){
     277             :           //itsLoopController.setPeakResidual(retval.asFloat("peakresidualnomask"));
     278             :           retval.define("peakresidual", futureRes);
     279             :           itsPreviousFutureRes=futureRes;
     280             :         }
     281             :       }
     282             :       return retval;
     283             :       }
     284             :     */
     285        2211 :     Record retval= initMinorCycle(itsImages);
     286             :     //    cerr << "INITMINOR retval" << retval << endl;
     287             : 
     288        2199 :     return retval;
     289        2199 :   }
     290        2468 :   Record SynthesisDeconvolver::initMinorCycle(std::shared_ptr<SIImageStore> imstor )
     291             :   {
     292        4936 :     LogIO os( LogOrigin("SynthesisDeconvolver","initMinorCycle",WHERE) );
     293        2468 :     Record returnRecord;
     294        2468 :     Timer timer;
     295        2468 :     Timer tim;
     296        2468 :     tim.mark();
     297             :     try {
     298             : 
     299             :       //os << "---------------------------------------------------- Init (?) Minor Cycles ---------------------------------------------" << LogIO::POST;
     300             : 
     301             :       Int nSubChans, nSubPols;
     302        2468 :       nSubChans = imstor->getShape()(3);
     303        2468 :       nSubPols = imstor->getShape()(2);
     304        2468 :       itsImages = imstor;
     305             : 
     306             :       // If a starting model exists, this will initialize the ImageStore with it. Will do this only once.
     307        2468 :       setStartingModel();
     308             : 
     309             :       //itsIterDone is currently only used by automask code so move this to inside setAutomask
     310             :       //itsIterDone += itsLoopController.getIterDone();
     311             : 
     312             :       //      setupMask();
     313             :       //
     314             :       Float masksum;
     315        2466 :       if( ! itsImages->hasMask() ) // i.e. if there is no existing mask to re-use...
     316             :           {
     317         599 :           masksum = -1.0;
     318             :       }
     319             :       else
     320             :         {
     321        1867 :           masksum = itsImages->getMaskSum();
     322        1863 :           itsImages->mask()->unlock();
     323             :         }
     324        2462 :       Bool validMask = ( masksum > 0 );
     325             : 
     326             :       //    os << LogIO::NORMAL3 << "****INITMINOR Masksum stuff "<< tim.real() << LogIO::POST;
     327             :       // tim.mark();
     328             : 
     329             :       // Calculate Peak Residual and Max Psf Sidelobe, and fill into SubIterBot.
     330        2462 :       Float peakresnomask = itsImages->getPeakResidual();
     331             :       // CAS-14201 : if the mask is all zero, then the peakresinmask is 0
     332        2462 :       Float peakresinmask= validMask ? itsImages->getPeakResidualWithinMask() : 0;
     333             :       //os << LogIO::NORMAL3 << "****INITMINOR residual peak "<< tim.real() << LogIO::POST;
     334             :       //tim.mark();
     335        2462 :       itsLoopController.setPeakResidual( validMask ? peakresinmask : peakresnomask );
     336             :       //os << LogIO::NORMAL3 << "****INITMINOR OTHER residual peak "<< tim.real() << LogIO::POST;
     337             :       //tim.mark();
     338        2462 :       itsLoopController.setPeakResidualNoMask( peakresnomask );
     339        2462 :       itsLoopController.setMaxPsfSidelobe( itsImages->getPSFSidelobeLevel() );
     340             : 
     341             :       //re-calculate current nsigma threhold
     342             :       //os<<"Calling calcRobustRMS ....syndeconv."<<LogIO::POST;
     343        2462 :       Float nsigmathresh = 0.0;
     344        2462 :       Bool useautomask = ( itsAutoMaskAlgorithm=="multithresh" ? true : false);
     345        2462 :       Int iterdone = itsLoopController.getIterDone();
     346             : 
     347             :       //cerr << "INIT automask " << useautomask << " alg " << itsAutoMaskAlgorithm << " sigma " << itsNsigma  << endl;
     348        2462 :       if ( itsNsigma >0.0) {
     349         143 :         itsMaskHandler->setPBMaskLevel(itsPBMask);
     350         143 :         Array<Double> medians, robustrms;
     351             :         // 2 cases to use existing stats.
     352             :         // 1. automask has run and so the image statistics record has filled
     353             :         // or
     354             :         // 2. no automask but for the first cycle but already initial calcRMS has ran to avoid duplicate
     355             :         //
     356             : 
     357             :         //cerr << "useauto " << useautomask << " nfields " << itsRobustStats.nfields() << " iterdone " << iterdone << endl;
     358         278 :         if ((useautomask && itsRobustStats.nfields()) ||
     359         135 :             (!useautomask && iterdone==0 && itsRobustStats.nfields()) ) {
     360          30 :            os <<LogIO::DEBUG1<<"automask on: check the current stats"<<LogIO::POST;
     361             :            //os<< "itsRobustStats nfield="<< itsRobustStats.nfields() << LogIO::POST;;
     362          30 :            if (itsRobustStats.isDefined("medabsdevmed")) {
     363           2 :              Array<Double> mads;
     364           2 :              itsRobustStats.get(RecordFieldId("medabsdevmed"), mads);
     365           2 :              os<<LogIO::DEBUG1<<"Using robust rms from automask ="<< mads*1.4826 <<LogIO::POST;
     366           2 :              robustrms = mads*1.4826;
     367           2 :            }
     368          28 :            else if(itsRobustStats.isDefined("robustrms")) {
     369          28 :              itsRobustStats.get(RecordFieldId("robustrms"), robustrms);
     370             :            }
     371          30 :            itsRobustStats.get(RecordFieldId("median"), medians);
     372             : 
     373             :         }
     374             :        else { // do own stats calculation
     375         113 :           timer.mark();
     376             : 
     377         113 :           os<<LogIO::DEBUG1<<"Calling calcRobustRMS .. "<<LogIO::POST;
     378         113 :           robustrms = itsImages->calcRobustRMS(medians, itsPBMask, itsFastNoise);
     379             :           os<< LogIO::NORMAL << "time for calcRobustRMS:  real "<< timer.real() << "s ( user " << timer.user()
     380         113 :              <<"s, system "<< timer.system() << "s)" << LogIO::POST;
     381             :           //reset itsRobustStats
     382             :           //cerr << "medians " << medians << " pbmask " << itsPBMask << endl;
     383             :           try {
     384             :             //os<<"current content of itsRobustStats nfields=="<<itsRobustStats.nfields()<<LogIO::POST;
     385         113 :             itsRobustStats.define(RecordFieldId("robustrms"), robustrms);
     386         113 :             itsRobustStats.define(RecordFieldId("median"), medians);
     387             :           }
     388           0 :           catch(AipsError &x) {
     389           0 :             throw( AipsError("Error in storing the robust image statistics") );
     390           0 :           }
     391             : 
     392             :                 //cerr << this << " DOING robust " << itsRobustStats << endl;
     393             :        }
     394             : 
     395             :         /***
     396             :         Array<Double> robustrms =kitsImages->calcRobustRMS(medians, itsPBMask, itsFastNoise);
     397             :         // Before the first iteration the iteration parameters have not been
     398             :         // set in SIMinorCycleController
     399             :         // Use nsigma pass to SynthesisDeconvolver directly for now...
     400             :         //Float nsigma = itsLoopController.getNsigma();
     401             :         ***/
     402             :         Double minval, maxval;
     403         143 :         IPosition minpos, maxpos;
     404             :         //Double maxrobustrms = max(robustrms);
     405         143 :         if(robustrms.empty())
     406           0 :           throw(AipsError("No valid values to deconvolve"));
     407             : 
     408         143 :         minMax(minval, maxval, minpos, maxpos, robustrms);
     409             : 
     410             :         //Float nsigmathresh = nsigma * (Float)robustrms(IPosition(1,0));
     411             :         //nsigmathresh = itsNsigma * (Float)maxrobustrms;
     412         143 :         String msg="";
     413         143 :         if (useautomask) {
     414          12 :           nsigmathresh = (Float)(medians(maxpos)) + itsNsigma * (Float)maxval;
     415          12 :           msg+=" (nsigma*rms + median)";
     416             :         }
     417             :         else {
     418         131 :           nsigmathresh = itsNsigma * (Float)maxval;
     419             :         }
     420         143 :         os << "Current nsigma threshold (maximum along spectral channels ) ="<<nsigmathresh<< msg <<LogIO::POST;
     421         143 :       }
     422             : 
     423        2462 :       itsLoopController.setNsigmaThreshold(nsigmathresh);
     424        2462 :       itsLoopController.setPBMask(itsPBMask);
     425        2462 :       itsLoopController.setFullSummary(itsFullSummary);
     426             : 
     427        2462 :       if ( itsAutoMaskAlgorithm=="multithresh" && !initializeChanMaskFlag ) {
     428          46 :         IPosition maskshp = itsImages->mask()->shape();
     429          46 :         Int nchan = maskshp(3);
     430          46 :         itsChanFlag=Vector<Bool>(nchan,False);
     431          46 :         initializeChanMaskFlag=True;
     432             :         // also initialize posmask, which tracks only positive (emission)
     433             : 
     434          46 :         if(!itsPosMask){
     435             :           //itsPosMask = TempImage<Float> (maskshp, itsImages->mask()->coordinates(),SDMaskHandler::memoryToUse());
     436          46 :           itsPosMask=itsImages->tempworkimage();
     437             :           //you don't want to modify this here...
     438             :           //It is set to 0.0 in SIImageStore first time it is created.
     439             :           //itsPosMask->set(0);
     440          46 :           itsPosMask->unlock();
     441             :         }
     442          46 :       }
     443        2462 :       os<<LogIO::DEBUG1<<"itsChanFlag.shape="<<itsChanFlag.shape()<<LogIO::POST;
     444             : 
     445             :       /*
     446             :       Array<Double> rmss = itsImages->calcRobustRMS();
     447             :       AlwaysAssert( rmss.shape()[0]>0 , AipsError);
     448             :       cout  << "madRMS : " << rmss << "  with shape : " << rmss.shape() << endl;
     449             :       //itsLoopController.setMadRMS( rmss[0] );
     450             :       */
     451             : 
     452        2462 :       if( itsMaskSum != masksum || masksum == 0.0 ) // i.e. mask has changed.
     453             :         {
     454        1515 :           itsMaskSum = masksum;
     455        1515 :           itsLoopController.setMaskSum( masksum );
     456             :         }
     457             :       else // mask has not changed...
     458             :         {
     459         947 :           itsLoopController.setMaskSum( -1.0 );
     460             :         }
     461        2462 :       float psfsidelobelevel = itsImages->getPSFSidelobeLevel();
     462             : 
     463        2462 :       returnRecord = itsLoopController.getCycleInitializationRecord();
     464             : 
     465             :       // cerr << "INIT record " << returnRecord << endl;
     466             : 
     467             :       //      itsImages->printImageStats();
     468        2462 :       os << " Absolute Peak residual within mask : " << peakresinmask << ", over full image : " << peakresnomask  << LogIO::POST;
     469        2462 :       itsImages->releaseLocks();
     470             : 
     471        2462 :       os << LogIO::DEBUG2 << "Initialized minor cycle. Returning returnRec" << LogIO::POST;
     472             : 
     473           6 :     } catch(AipsError &x) {
     474           6 :       throw( AipsError("Error initializing the Minor Cycle for "  + itsImageName + " : "+x.getMesg()) );
     475           6 :     }
     476             : 
     477        4924 :     return returnRecord;
     478        2474 :   }
     479          15 :   void SynthesisDeconvolver::setChanFlag(const Vector<Bool>& chanflag){
     480             :     //ignore if it has not been given a size yet in initminorcycle
     481          15 :     if(itsChanFlag.nelements()==0)
     482           0 :       return;
     483          15 :     if(itsChanFlag.nelements() != chanflag.nelements())
     484           0 :       throw(AipsError("cannot set chan flags for different number of channels"));
     485          15 :     itsChanFlag =chanflag;
     486             : 
     487             :   }
     488         263 :   Vector<Bool> SynthesisDeconvolver::getChanFlag(){
     489         263 :     return itsChanFlag;
     490             :   }
     491          15 :   void SynthesisDeconvolver::setRobustStats(const Record& rec){
     492          15 :     itsRobustStats=Record();
     493          15 :     itsRobustStats=rec;
     494             : 
     495          15 :   }
     496         263 :   Record SynthesisDeconvolver::getRobustStats(){
     497         263 :     return itsRobustStats;
     498             :   }
     499             : 
     500           0 :   Record SynthesisDeconvolver::interactiveGUI(Record& iterRec)
     501             :   {
     502           0 :     LogIO os( LogOrigin("SynthesisDeconvolver","interactiveGUI",WHERE) );
     503           0 :     Record returnRecord;
     504             : 
     505             :     try {
     506             : 
     507             :       // Check that required parameters are present in the iterRec.
     508           0 :       Int niter=0,cycleniter=0,iterdone;
     509           0 :       Float threshold=0.0, cyclethreshold=0.0;
     510           0 :       if( iterRec.isDefined("niter") &&
     511           0 :           iterRec.isDefined("cycleniter") &&
     512           0 :           iterRec.isDefined("threshold") &&
     513           0 :           iterRec.isDefined("cyclethreshold") &&
     514           0 :           iterRec.isDefined("iterdone") )
     515             :         {
     516           0 :           iterRec.get("niter", niter);
     517           0 :           iterRec.get("cycleniter", cycleniter);
     518           0 :           iterRec.get("threshold", threshold);
     519           0 :           iterRec.get("cyclethreshold", cyclethreshold);
     520           0 :           iterRec.get("iterdone",iterdone);
     521             :         }
     522           0 :       else throw(AipsError("SD::interactiveGui() needs valid niter, cycleniter, threshold to start up."));
     523             : 
     524           0 :       if( ! itsImages ) itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
     525             : 
     526             :       //      SDMaskHandler masker;
     527           0 :       String strthresh = String::toString(threshold)+"Jy";
     528           0 :       String strcycthresh = String::toString(cyclethreshold)+"Jy";
     529             : 
     530           0 :       if( itsMaskString.length()>0 ) {
     531           0 :           itsMaskHandler->fillMask( itsImages, itsMaskString );
     532             :       }
     533             : 
     534           0 :       Int iterleft = niter - iterdone;
     535           0 :       if( iterleft<0 ) iterleft=0;
     536             : 
     537           0 :       Int stopcode = itsMaskHandler->makeInteractiveMask( itsImages, iterleft, cycleniter, strthresh, strcycthresh );
     538             : 
     539           0 :       Quantity qa;
     540           0 :       casacore::Quantity::read(qa,strthresh);
     541           0 :       threshold = qa.getValue(Unit("Jy"));
     542           0 :       casacore::Quantity::read(qa,strcycthresh);
     543           0 :       cyclethreshold = qa.getValue(Unit("Jy"));
     544             : 
     545           0 :       itsIsMaskLoaded=true;
     546             : 
     547           0 :       returnRecord.define( RecordFieldId("actioncode"), stopcode );
     548           0 :       returnRecord.define( RecordFieldId("niter"), iterdone + iterleft );
     549           0 :       returnRecord.define( RecordFieldId("cycleniter"), cycleniter );
     550           0 :       returnRecord.define( RecordFieldId("threshold"), threshold );
     551           0 :       returnRecord.define( RecordFieldId("cyclethreshold"), cyclethreshold );
     552             : 
     553           0 :     } catch(AipsError &x) {
     554           0 :       throw( AipsError("Error in Interactive GUI : "+x.getMesg()) );
     555           0 :     }
     556           0 :     return returnRecord;
     557           0 :   }
     558             : 
     559             : 
     560           5 :   void SynthesisDeconvolver::setMinorCycleControl(const Record& minorCycleControlRec){
     561             :     //Don't know what itsloopcontroller does not need a const record;
     562           5 :     Record lala=minorCycleControlRec;
     563           5 :     itsLoopController.setCycleControls(lala);
     564             : 
     565           5 :   }
     566             : 
     567         934 :   Record SynthesisDeconvolver::executeMinorCycle(Record& minorCycleControlRec)
     568             :   {
     569             :     // LogIO os( LogOrigin("SynthesisDeconvolver","executeMinorCycle",WHERE) );
     570             : 
     571             : 
     572             :     //    itsImages->printImageStats();
     573         934 :     SynthesisUtilMethods::getResource("Start Deconvolver");
     574             :     ///if cube execute cube deconvolution...check on residual shape as itsimagestore return 0 shape sometimes
     575         934 :     if(!itsImages)
     576           0 :       throw(AipsError("Initminor Cycle has not been called yet"));
     577         934 :     if(itsImages->residual()->shape()[3]> 1){
     578         253 :      return  executeCubeMinorCycle(minorCycleControlRec);
     579             :     }
     580             : 
     581             :     //  os << "---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------" << LogIO::POST;
     582         681 :     return executeCoreMinorCycle(minorCycleControlRec);
     583             :     SynthesisUtilMethods::getResource("End Deconvolver");
     584             :   }
     585         934 :   Record SynthesisDeconvolver::executeCoreMinorCycle(Record& minorCycleControlRec)
     586             :   {
     587             : 
     588         934 :     Record returnRecord;
     589             :     try {
     590             :       //      if ( !itsIsInteractive ) setAutoMask();
     591             :       //cerr << "MINORCYCLE control Rec " << minorCycleControlRec << endl;
     592         934 :       itsLoopController.setCycleControls(minorCycleControlRec);
     593         934 :       bool automaskon (false);
     594         934 :       if (itsAutoMaskAlgorithm=="multithresh") {
     595          34 :         automaskon=true;
     596             :       }
     597             :       //itsDeconvolver->deconvolve( itsLoopController, itsImages, itsDeconvolverId, automaskon, itsFastNoise );
     598             :       // include robust stats rec
     599         942 :       itsDeconvolver->deconvolve( itsLoopController, itsImages, itsDeconvolverId, automaskon, itsFastNoise, itsRobustStats, itsFullSummary );
     600             : 
     601         926 :       returnRecord = itsLoopController.getCycleExecutionRecord();
     602             : 
     603             :       //scatterModel(); // This is a no-op for the single-node case.
     604             : 
     605         926 :       itsImages->releaseLocks();
     606             : 
     607           8 :     } catch(AipsError &x) {
     608           8 :       throw( AipsError("Error in running Minor Cycle : "+x.getMesg()) );
     609           8 :     }
     610             : 
     611             : 
     612             : 
     613         926 :     return returnRecord;
     614           8 :   }
     615         263 :   Record SynthesisDeconvolver::executeCubeMinorCycle(Record& minorCycleControlRec, const Int AutoMaskFlag)
     616             :   {
     617         526 :         LogIO os( LogOrigin("SynthesisDeconvolver","executeCubeMinorCycle",WHERE) );
     618         263 :     Record returnRecord;
     619         263 :     Int doAutoMask=AutoMaskFlag;
     620             : 
     621         263 :     SynthesisUtilMethods::getResource("Start Deconvolver");
     622             : 
     623             :     try {
     624             :       //      if ( !itsIsInteractive ) setAutoMask();
     625         263 :       if(doAutoMask < 1){
     626         253 :         itsLoopController.setCycleControls(minorCycleControlRec);
     627             :       }
     628          10 :       else if(doAutoMask==1){
     629          10 :         minorCycleControlRec=itsPreviousIterBotRec_p;
     630             :       }
     631         263 :       CubeMinorCycleAlgorithm cmc;
     632             :       //casa::applicator.defineAlgorithm(cmc);
     633             :       ///argv and argc are needed just to callthe right overloaded init
     634         263 :       Int argc=1;
     635         263 :       char **argv=nullptr;
     636         263 :       applicator.init(argc, argv);
     637         263 :       if(applicator.isController()){
     638         263 :         os << ((AutoMaskFlag != 1) ? "---------------------------------------------------- Run Minor Cycle Iterations  ---------------------------------------------" : "---------------------------------------------------- Run Automask  ---------------------------------------------" )<< LogIO::POST;
     639             :         /*{///TO BE REMOVED
     640             :           LatticeExprNode le( sum( *(itsImages->mask()) ) );
     641             :           os << LogIO::WARN << "#####Sum of mask BEFORE minor cycle " << le.getFloat() << endl;
     642             :             }
     643             :         */
     644         263 :         Timer tim;
     645         263 :         tim.mark();
     646             :         //itsImages->printImageStats();
     647             :         // Add itsIterdone to be sent to child processes ...needed for automask
     648             :         //cerr << "before record " << itsIterDone << " loopcontroller " << itsLoopController.getIterDone() << endl;
     649         263 :         minorCycleControlRec.define("iterdone", itsIterDone);
     650         263 :         if(doAutoMask < 0) // && itsPreviousIterBotRec_p.nfields() >0)
     651         253 :           doAutoMask=0;
     652         263 :         minorCycleControlRec.define("onlyautomask",doAutoMask);
     653         263 :         if(itsPosMask){
     654          15 :           minorCycleControlRec.define("posmaskname", itsPosMask->name());
     655             :         }
     656             :         //Int numprocs = applicator.numProcs();
     657             :         //cerr << "Number of procs: " << numprocs << endl;
     658             : 
     659         263 :           Int numchan=itsImages->residual()->shape()[3];
     660         263 :           Vector<Int> startchans;
     661         263 :           Vector<Int> endchans;
     662         263 :           Int numblocks=numblockchans(startchans, endchans);
     663         263 :           String psfname=itsImages->psf()->name();
     664             : 
     665         263 :           Float psfsidelobelevel=itsImages->getPSFSidelobeLevel();
     666         263 :           String residualname=itsImages->residual()->name();
     667         263 :           String maskname=itsImages->mask()->name();
     668         263 :           String modelname=itsImages->model()->name();
     669             :           ////Need the pb too as calcrobustrms in SynthesisDeconvolver uses it
     670         263 :           String pbname="";
     671         263 :           if(itsImages->hasPB())
     672         263 :             pbname=itsImages->pb()->name();
     673         263 :           itsImages->releaseLocks();
     674             :           ///in lieu of = operator go via record
     675             :           // need to create a proper = operator for SynthesisParamsDeconv
     676         263 :           SynthesisParamsDeconv decpars;
     677             :           ///Will have to create a = operator...right now kludging
     678             :           ///from record has a check that has to be bypassed for just the
     679             :           /// usage as a = operator
     680             :           {
     681         263 :             String tempMaskString= itsDecPars.maskString;
     682         263 :             itsDecPars.maskString="";
     683         263 :             decpars.fromRecord(itsDecPars.toRecord());
     684             :             //return itsDecPars back to its original state
     685         263 :             itsDecPars.maskString=tempMaskString;
     686         263 :           }
     687             :           ///remove starting model as already dealt with in this deconvolver
     688         263 :           decpars.startModel="";
     689             :           ///masking is dealt already by this deconvolver so mask image
     690             :           //is all that is needed which is sent as maskname to subdeconvolver
     691         263 :           decpars.maskString="";
     692         263 :           (decpars.maskList).resize();
     693         263 :           Record decParsRec = decpars.toRecord();
     694             : 
     695             :           /////Now we loop over channels and deconvolve each
     696             :           ///If we do want to do block of channels rather than 1 channel
     697             :           ///at a time chanRange can take that and the trigger to call this
     698             :           ///function in executeMinorCycle has to change.
     699         263 :           Int rank(0);
     700             :           Bool assigned;
     701         263 :           Bool allDone(false);
     702         263 :           Vector<Int> chanRange(2);
     703             :           //Record beamsetRec;
     704         263 :           Vector<Bool> retvals(numblocks, False);
     705         263 :           Vector<Bool> chanFlag(0);
     706         263 :           if(itsChanFlag.nelements()==0){
     707         212 :             itsChanFlag.resize(numchan);
     708         212 :             itsChanFlag.set(False);
     709             :           }
     710         263 :           Record chanflagRec;
     711         263 :           Int indexofretval=0;
     712         526 :           for (Int k=0; k < numblocks; ++k) {
     713             :             //os << LogIO::DEBUG1 << "deconvolving channel "<< k << LogIO::POST;
     714         263 :             assigned=casa::applicator.nextAvailProcess(cmc, rank);
     715             :             //cerr << "assigned "<< assigned << endl;
     716         263 :             while(!assigned) {
     717             :               //cerr << "SErial ? " << casa::applicator.isSerial() << endl;
     718           0 :               rank = casa::applicator.nextProcessDone(cmc, allDone);
     719             :               //cerr << "while rank " << rank << endl;
     720             :               //receiving output of CubeMinorCycleAlgorithm::put
     721             :               //#1
     722           0 :               Vector<Int> chanRangeProcessed;
     723           0 :               casa::applicator.get(chanRangeProcessed);
     724             :               //#2
     725             : 
     726           0 :               Record chanflagRec;
     727           0 :               casa::applicator.get(chanflagRec);
     728             : 
     729             :               //#3
     730           0 :               Record retval;
     731           0 :               casa::applicator.get(retval);
     732             : 
     733           0 :               Vector<Bool> retchanflag;
     734           0 :               chanflagRec.get("chanflag", retchanflag);
     735           0 :               if(retchanflag.nelements() >0)
     736           0 :                 itsChanFlag(Slice(chanRangeProcessed[0], chanRangeProcessed[1]-chanRangeProcessed[0]+1))=retchanflag;
     737           0 :               Record substats=chanflagRec.asRecord("statsrec");
     738           0 :               setSubsetRobustStats(substats, chanRangeProcessed[0], chanRangeProcessed[1], numchan);
     739             : 
     740           0 :               retvals(indexofretval)=(retval.nfields() > 0);
     741           0 :               ++indexofretval;
     742             :               ///might need to merge these retval
     743           0 :               if(doAutoMask <1)
     744           0 :                 mergeReturnRecord(retval, returnRecord, chanRangeProcessed[0]);
     745             :               /*if(retval.nfields())
     746             :                 cerr << k << "deconv rank " << rank << " successful " << endl;
     747             :                else
     748             :                 cerr << k << "deconv rank " << rank << " failed " << endl;
     749             :               */
     750             :               //cerr <<"rank " << rank << " return rec "<< retval << endl;
     751           0 :               assigned = casa::applicator.nextAvailProcess(cmc, rank);
     752             : 
     753           0 :             }
     754             : 
     755             :             ///send process info
     756             :             // put dec sel params #1
     757         263 :             applicator.put(decParsRec);
     758             :             // put itercontrol  params #2
     759         263 :             applicator.put(minorCycleControlRec);
     760             :             // put which channel to process #3
     761         263 :             chanRange[0]=startchans[k];  chanRange[1]=endchans[k];
     762         263 :             applicator.put(chanRange);
     763             :             // psf  #4
     764         263 :             applicator.put(psfname);
     765             :             // residual #5
     766         263 :             applicator.put(residualname);
     767             :             // model #6
     768         263 :             applicator.put(modelname);
     769             :             // mask #7
     770         263 :             applicator.put(maskname);
     771             :             //pb #8
     772         263 :             applicator.put(pbname);
     773             :             //#9 psf side lobe
     774         263 :             applicator.put(psfsidelobelevel);
     775             :             //# put chanflag
     776         263 :             chanFlag.resize();
     777         263 :             chanFlag=itsChanFlag(IPosition(1, chanRange[0]), IPosition(1, chanRange[1]));
     778             : 
     779         263 :             chanflagRec.define("chanflag", chanFlag);
     780         263 :             Record statrec=getSubsetRobustStats(chanRange[0], chanRange[1]);
     781         263 :             chanflagRec.defineRecord("statsrec", statrec);
     782         263 :             applicator.put(chanflagRec);
     783             :             /// Tell worker to process it
     784         263 :             applicator.apply(cmc);
     785             : 
     786         263 :           }
     787             :           // Wait for all outstanding processes to return
     788         263 :           rank = casa::applicator.nextProcessDone(cmc, allDone);
     789         526 :           while (!allDone) {
     790             : 
     791         263 :             Vector<Int> chanRangeProcessed;
     792         263 :             casa::applicator.get(chanRangeProcessed);
     793         263 :             Record chanflagRec;
     794         263 :             casa::applicator.get(chanflagRec);
     795         263 :             Record retval;
     796         263 :             casa::applicator.get(retval);
     797         263 :             Vector<Bool> retchanflag;
     798         263 :             chanflagRec.get("chanflag", retchanflag);
     799         263 :             if(retchanflag.nelements() >0)
     800          15 :               itsChanFlag(Slice(chanRangeProcessed[0], chanRangeProcessed[1]-chanRangeProcessed[0]+1))=retchanflag;
     801         263 :             Record substats=chanflagRec.asRecord("statsrec");
     802         263 :             setSubsetRobustStats(substats, chanRangeProcessed[0], chanRangeProcessed[1], numchan);
     803         263 :             retvals(indexofretval)=(retval.nfields() > 0);
     804         263 :             ++indexofretval;
     805         263 :             if(doAutoMask < 1)
     806         253 :               mergeReturnRecord(retval, returnRecord, chanRangeProcessed[0]);
     807         263 :             if(retval.nfields() >0)
     808             :               //cerr << "deconv remainder rank " << rank << " successful " << endl;
     809         263 :               cerr << "";
     810             :             else
     811           0 :               cerr << "deconv remainder rank " << rank << " failed " << endl;
     812             : 
     813         263 :             rank = casa::applicator.nextProcessDone(cmc, allDone);
     814         263 :             if(casa::applicator.isSerial())
     815         263 :               allDone=true;
     816             : 
     817         263 :           }
     818             : 
     819         263 :           if(anyEQ(retvals, False))
     820           0 :             throw(AipsError("one of more section of the cube failed in deconvolution"));
     821         263 :           if(doAutoMask < 1){
     822         253 :             itsLoopController.incrementMinorCycleCount(returnRecord.asInt("iterdone"));
     823         253 :             itsIterDone+=returnRecord.asInt("iterdone");
     824             :           }
     825         263 :           itsPreviousIterBotRec_p=Record();
     826         263 :           itsPreviousIterBotRec_p=minorCycleControlRec;
     827             :           /*{///TO BE REMOVED
     828             :           LatticeExprNode le( sum( *(itsImages->mask()) ) );
     829             :           os << LogIO::WARN << "#####Sum of mask AFTER minor cycle " << le.getFloat()  << "loopcontroller iterdeconv " << itsLoopController.getIterDone() << endl;
     830             :           }*/
     831             : 
     832         263 :       }///end of if controller
     833             :       /////////////////////////////////////////////////
     834             : 
     835             :       //scatterModel(); // This is a no-op for the single-node case.
     836             : 
     837         263 :       itsImages->releaseLocks();
     838             : 
     839         263 :     } catch(AipsError &x) {
     840             :       //MPI_Abort(MPI_COMM_WORLD, 6);
     841           0 :       throw( AipsError("Error in running Minor Cycle : "+x.getMesg()) );
     842           0 :     }
     843             : 
     844         263 :     SynthesisUtilMethods::getResource("End Deconvolver");
     845             : 
     846         526 :     return returnRecord;
     847         263 :   }
     848             : 
     849             :   // Restore Image.
     850         734 :   void SynthesisDeconvolver::restore()
     851             :   {
     852        1468 :     LogIO os( LogOrigin("SynthesisDeconvolver","restoreImage",WHERE) );
     853             : 
     854         734 :     if( ! itsImages )
     855             :       {
     856           7 :         itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
     857             :       }
     858             : 
     859         734 :     SynthesisUtilMethods::getResource("Restoration");
     860             : 
     861         734 :     itsDeconvolver->restore(itsImages);
     862         734 :     itsImages->releaseLocks();
     863             : 
     864         734 :   }
     865             : 
     866         253 :   void SynthesisDeconvolver::mergeReturnRecord(const Record& inRec, Record& outRec, const Int chan){
     867             : 
     868             :     ///Something has to be done about what is done in SIIterBot_state::mergeMinorCycleSummary if it is needed
     869             :     //int nSummaryFields = SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     870             :     //int nSummaryFields = SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     871         253 :     int nSummaryFields = !itsFullSummary ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     872             : 
     873         253 :     Matrix<Double> summaryminor(nSummaryFields,0);
     874         253 :     if(outRec.isDefined("summaryminor"))
     875           0 :       summaryminor=Matrix<Double>(outRec.asArrayDouble("summaryminor"));
     876         253 :     Matrix<Double> subsummaryminor;
     877         253 :     if(inRec.isDefined("summaryminor"))
     878         253 :       subsummaryminor=Matrix<Double>(inRec.asArrayDouble("summaryminor"));
     879         253 :     if(subsummaryminor.nelements() !=0){
     880             :       ///The 6th element is supposed to be the subimage id
     881         253 :       subsummaryminor.row(5)= subsummaryminor.row(5)+(Double(chan));
     882         253 :       Matrix<Double> newsummary(nSummaryFields, summaryminor.shape()[1]+subsummaryminor.shape()[1]);
     883         253 :       Int ocol=0;
     884         253 :       for (Int col=0; col< summaryminor.shape()[1]; ++col, ++ocol)
     885           0 :         newsummary.column(ocol)=summaryminor.column(col);
     886        2678 :       for (Int col=0; col< subsummaryminor.shape()[1]; ++col, ++ocol)
     887        2425 :         newsummary.column(ocol)=subsummaryminor.column(col);
     888         253 :       summaryminor.resize(newsummary.shape());
     889         253 :       summaryminor=newsummary;
     890         253 :     }
     891         253 :     outRec.define("summaryminor", summaryminor);
     892             :     //cerr << "inRec summ minor " << inRec.asArrayDouble("summaryminor") << endl;
     893             :     //cerr << "outRec summ minor " << summaryminor << endl;
     894         253 :     outRec.define("iterdone", Int(inRec.asInt("iterdone")+ (outRec.isDefined("iterdone") ? outRec.asInt("iterdone"): Int(0))));
     895         253 :     outRec.define("maxcycleiterdone", outRec.isDefined("maxcycleiterdone") ? max(inRec.asInt("maxcycleiterdone"), outRec.asInt("maxcycleiterdone")) :inRec.asInt("maxcycleiterdone")) ;
     896             : 
     897         253 :     outRec.define("peakresidual", outRec.isDefined("peakresidual") ? max(inRec.asFloat("peakresidual"), outRec.asFloat("peakresidual")) :inRec.asFloat("peakresidual")) ;
     898             : 
     899             :     ///is not necessarily defined it seems
     900         253 :     Bool updatedmodelflag=False;
     901         253 :     if(inRec.isDefined("updatedmodelflag"))
     902         253 :       inRec.get("updatedmodelflag", updatedmodelflag);
     903         253 :     outRec.define("updatedmodelflag", outRec.isDefined("updatedmodelflag") ? updatedmodelflag || outRec.asBool("updatedmodelflag") : updatedmodelflag) ;
     904             : 
     905             : 
     906             : 
     907         253 :   }
     908             :   // get channel blocks
     909         263 :   Int SynthesisDeconvolver::numblockchans(Vector<Int>& startchans, Vector<Int>& endchans){
     910         263 :     Int nchan=itsImages->residual()->shape()[3];
     911             :     //roughly 8e6 pixel to deconvolve per lock/process is a  minimum
     912         263 :     Int optchan= 8e6/(itsImages->residual()->shape()[0])/(itsImages->residual()->shape()[1]);
     913             :     // cerr << "OPTCHAN" << optchan  << endl;
     914         263 :     if(optchan < 10) optchan=10;
     915         263 :     Int nproc= applicator.numProcs() < 2 ? 1 : applicator.numProcs()-1;
     916             :     /*if(nproc==1){
     917             :       startchans.resize(1);
     918             :       endchans.resize(1);
     919             :       startchans[0]=0;
     920             :       endchans[0]=nchan-1;
     921             :       return 1;
     922             :       }
     923             :     */
     924         263 :     Int blksize= nchan/nproc > optchan ? optchan : Int( std::floor(Float(nchan)/Float(nproc)));
     925         263 :     if(blksize< 1) blksize=1;
     926         263 :     Int nblk=Int(nchan/blksize);
     927         263 :     startchans.resize(nblk);
     928         263 :     endchans.resize(nblk);
     929         526 :     for (Int k=0; k < nblk; ++k){
     930         263 :       startchans[k]= k*blksize;
     931         263 :       endchans[k]=(k+1)*blksize-1;
     932             :     }
     933         263 :     if(endchans[nblk-1] < (nchan-1)){
     934           0 :       startchans.resize(nblk+1,True);
     935           0 :       startchans[nblk]=endchans[nblk-1]+1;
     936           0 :       endchans.resize(nblk+1,True);
     937           0 :       endchans[nblk]=nchan-1;
     938           0 :       ++nblk;
     939             :     }
     940             :     //cerr << "nblk " << nblk << " beg " << startchans << " end " << endchans << endl;
     941         263 :     return nblk;
     942             :   }
     943             : 
     944             :   // pbcor Image.
     945          43 :   void SynthesisDeconvolver::pbcor()
     946             :   {
     947          86 :     LogIO os( LogOrigin("SynthesisDeconvolver","pbcor",WHERE) );
     948             : 
     949          43 :     if( ! itsImages )
     950             :       {
     951           0 :         itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
     952             :       }
     953             : 
     954          43 :     itsDeconvolver->pbcor(itsImages);
     955          43 :     itsImages->releaseLocks();
     956             : 
     957          43 :   }
     958             : 
     959             : 
     960             : 
     961             :   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     962             :   ////    Internal Functions start here.  These are not visible to the tool layer.
     963             :   /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     964             : 
     965         798 :   std::shared_ptr<SIImageStore> SynthesisDeconvolver::makeImageStore( String imagename, Bool noRequireSumwt )
     966             :   {
     967         798 :     std::shared_ptr<SIImageStore> imstore;
     968         798 :     if( itsDeconvolver->getAlgorithmName() == "mtmfs" )
     969         148 :       {  imstore.reset( new SIImageStoreMultiTerm( imagename, itsDeconvolver->getNTaylorTerms(), true, noRequireSumwt ) ); }
     970             :     else
     971         650 :       {  imstore.reset( new SIImageStore( imagename, true, noRequireSumwt ) ); }
     972             : 
     973         798 :     return imstore;
     974             : 
     975           0 :   }
     976             : 
     977             : 
     978             :   // #############################################
     979             :   // #############################################
     980             :   // #############################################
     981             :   // #############################################
     982             : 
     983             :   // Set a starting model.
     984        2468 :   void SynthesisDeconvolver::setStartingModel()
     985             :   {
     986        4936 :     LogIO os( LogOrigin("SynthesisDeconvolver","setStartingModel",WHERE) );
     987             : 
     988        2468 :     if(itsAddedModel==true) {return;}
     989             : 
     990             :     try
     991             :       {
     992             : 
     993         908 :         if( itsStartingModelNames.nelements()>0 && itsImages )
     994             :           {
     995             :             //      os << "Setting " << itsStartingModelNames << " as starting model for deconvolution " << LogIO::POST;
     996          11 :             itsImages->setModelImage( itsStartingModelNames );
     997             :           }
     998             : 
     999         906 :         itsAddedModel=true;
    1000             : 
    1001             :       }
    1002           2 :     catch(AipsError &x)
    1003             :       {
    1004           2 :         throw( AipsError("Error in setting  starting model for deconvolution: "+x.getMesg()) );
    1005           2 :       }
    1006             : 
    1007        2468 :   }
    1008             : 
    1009             :   // Set mask
    1010        1414 :   Bool SynthesisDeconvolver::setupMask()
    1011             :   {
    1012             : 
    1013             :     ////Remembet this has to be called only after initMinorCycle
    1014        2828 :     LogIO os( LogOrigin("SynthesisDeconvolver","setupMask",WHERE) );
    1015        1414 :     if(!itsImages)
    1016           0 :       throw(AipsError("Initminor Cycle has not been called yet"));
    1017             : 
    1018        1414 :     Bool maskchanged=False;
    1019             :     //debug
    1020        1414 :     if( itsIsMaskLoaded==false ) {
    1021             :       // use mask(s)
    1022         660 :       if(  itsMaskList[0] != "" || itsMaskType == "pb" || itsAutoMaskAlgorithm != "" ) {
    1023             :         // Skip automask for non-interactive mode.
    1024         128 :         if ( itsAutoMaskAlgorithm != "") { // && itsIsInteractive) {
    1025          52 :           os << "[" << itsImages->getName() << "] Setting up an auto-mask"<<  ((itsPBMask>0.0)?" within PB mask limit ":"") << LogIO::POST;
    1026             :           ////For Cubes this is done in CubeMinorCycle
    1027             :           //cerr << this << "SETUP mask " << itsRobustStats << endl;
    1028          52 :           if(itsImages->residual()->shape()[3] ==1)
    1029          42 :             setAutoMask();
    1030          10 :           else if((itsImages->residual()->shape()[3] >1)){
    1031          10 :             Record dummy;
    1032          10 :             executeCubeMinorCycle(dummy, 1);
    1033          10 :           }
    1034             :           /***
    1035             :           if ( itsPBMask > 0.0 ) {
    1036             :             itsMaskHandler->autoMaskWithinPB( itsImages, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsPBMask);
    1037             :           }
    1038             :           else {
    1039             :             cerr<<"automask by automask.."<<endl;
    1040             :             itsMaskHandler->autoMask( itsImages, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam);
    1041             :           }
    1042             :           ***/
    1043             :         }
    1044             : 
    1045             :         //else if( itsMaskType=="user" && itsMaskList[0] != "" ) {
    1046         128 :         if( itsMaskType=="user" && itsMaskList[0] != "" ) {
    1047          63 :           os << "[" << itsImages->getName() << "] Setting up a mask from " << itsMaskList  <<  ((itsPBMask>0.0)?" within PB mask limit ":"") << LogIO::POST;
    1048             : 
    1049          63 :           itsMaskHandler->fillMask( itsImages, itsMaskList);
    1050          63 :           if( itsPBMask > 0.0 ) {
    1051           7 :             itsMaskHandler->makePBMask(itsImages, itsPBMask, True);
    1052             :           }
    1053             :         }
    1054          65 :         else if( itsMaskType=="pb") {
    1055          13 :           os << "[" << itsImages->getName() << "] Setting up a PB based mask" << LogIO::POST;
    1056          18 :           itsMaskHandler->makePBMask(itsImages, itsPBMask);
    1057             :         }
    1058             : 
    1059         123 :         os << "----------------------------------------------------------------------------------------------------------------------------------------" << LogIO::POST;
    1060             : 
    1061             :         }
    1062             :         else {
    1063             :             // new im statistics creates an empty mask and need to take care of that case
    1064         532 :             Bool emptyMask(False);
    1065         532 :             if( itsImages->hasMask() )
    1066             :             {
    1067             :                 // CAS-14203 - Check if mask is empty AND user didn't specify an empty mask
    1068          36 :                 if (itsImages->getMaskSum()==0.0 && itsMaskList[0] != "") {
    1069           0 :                     emptyMask=True;
    1070             :                 }
    1071             :             }
    1072         532 :             if( ! itsImages->hasMask() || emptyMask ) // i.e. if there is no existing mask to re-use...
    1073             :             {
    1074         496 :                 LatticeLocker lock1 (*(itsImages->mask()), FileLocker::Write);
    1075         496 :                 if( itsIsInteractive ) itsImages->mask()->set(0.0);
    1076         496 :                 else itsImages->mask()->set(1.0);
    1077         496 :                 os << "[" << itsImages->getName() << "] Initializing new mask to " << (itsIsInteractive?"0.0 for interactive drawing":"1.0 for the full image") << LogIO::POST;
    1078         496 :             }
    1079             :             else {
    1080          36 :                 os << "[" << itsImages->getName() << "] Initializing to existing mask" << LogIO::POST;
    1081             :             }
    1082             :         }
    1083             : 
    1084             :       // If anything other than automasking, don't re-make the mask here.
    1085         655 :       if ( itsAutoMaskAlgorithm == "" )
    1086         603 :         {       itsIsMaskLoaded=true; }
    1087             : 
    1088             :       // Get the number of mask pixels (sum) and send to the logger.
    1089         655 :       Float masksum = itsImages->getMaskSum();
    1090         655 :       Float npix = (itsImages->getShape()).product();
    1091             : 
    1092             :       //Int npix2 = 20000*20000*16000*4;
    1093             :       //Float npix2f = 20000*20000*16000*4;
    1094             : 
    1095             :       //cout << " bigval : " << npix2 << " and " << npix2f << endl;
    1096             : 
    1097         655 :       os << "[" << itsImages->getName() << "] Number of pixels in the clean mask : " << masksum << " out of a total of " << npix << " pixels. [ " << 100.0 * masksum/npix << " % ]" << LogIO::POST;
    1098             : 
    1099         655 :       maskchanged=True;
    1100             :     }
    1101             :     else {
    1102             :     }
    1103             : 
    1104        1409 :     itsImages->mask()->unlock();
    1105        1409 :     return maskchanged;
    1106        1414 : }
    1107             : 
    1108          52 :   void SynthesisDeconvolver::setAutoMask()
    1109             :   {
    1110             :      //modify mask using automask otherwise no-op
    1111          52 :      if ( itsAutoMaskAlgorithm != "" )  {
    1112          52 :        itsIterDone += itsLoopController.getIterDone();
    1113             : 
    1114             : 
    1115             : 
    1116          52 :        Bool isThresholdReached = itsLoopController.isThresholdReached();
    1117             :              //cerr << this << " setAuto " << itsRobustStats << endl;
    1118         104 :        LogIO os( LogOrigin("SynthesisDeconvolver","setAutoMask",WHERE) );
    1119          52 :        os << "Generating AutoMask" << LogIO::POST;
    1120             :        //os << LogIO::WARN << "#####ItsIterDone value " << itsIterDone << endl;
    1121             : 
    1122             :        //os << "itsMinPercentChnage = " << itsMinPercentChange<< LogIO::POST;
    1123             :        //cerr << "SUMa of chanFlag before " << ntrue(itsChanFlag) << endl;
    1124          52 :        if ( itsPBMask > 0.0 ) {
    1125             :          //itsMaskHandler->autoMaskWithinPB( itsImages, itsPosMask, itsIterDone, itsChanFlag, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust,  itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold,itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached, itsPBMask);
    1126             :          //pass robust stats
    1127           0 :          itsMaskHandler->autoMaskWithinPB( itsImages, *itsPosMask, itsIterDone, itsChanFlag, itsRobustStats, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust,  itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold,itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached, itsPBMask);
    1128             :        }
    1129             :        else {
    1130             :          //itsMaskHandler->autoMask( itsImages, itsPosMask, itsIterDone, itsChanFlag,itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold, itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached );
    1131             : 
    1132             :         // pass robust stats
    1133          52 :         itsMaskHandler->autoMask( itsImages, *itsPosMask, itsIterDone, itsChanFlag, itsRobustStats, itsAutoMaskAlgorithm, itsMaskThreshold, itsFracOfPeak, itsMaskResolution, itsMaskResByBeam, itsNMask, itsAutoAdjust, itsSidelobeThreshold, itsNoiseThreshold, itsLowNoiseThreshold, itsNegativeThreshold, itsCutThreshold, itsSmoothFactor, itsMinBeamFrac, itsGrowIterations, itsDoGrowPrune, itsMinPercentChange, itsVerbose, itsFastNoise, isThresholdReached );
    1134             :        }
    1135             :        //cerr <<this << " SETAutoMask " << itsRobustStats << endl;
    1136             :        //cerr << "SUM of chanFlag AFTER " << ntrue(itsChanFlag) << endl;
    1137          52 :      }
    1138          52 :   }
    1139             : 
    1140             :   // check if restoring beam is reasonable
    1141         572 :   void SynthesisDeconvolver::checkRestoringBeam()
    1142             :   {
    1143        1144 :     LogIO os( LogOrigin("SynthesisDeconvolver","checkRestoringBeam",WHERE) );
    1144             :     //check for a bad restoring beam
    1145         572 :     GaussianBeam beam;
    1146             : 
    1147         572 :     if( ! itsImages ) itsImages = makeImageStore( itsImageName, itsNoRequireSumwt );
    1148         572 :     ImageInfo psfInfo = itsImages->psf()->imageInfo();
    1149         572 :     if (psfInfo.hasSingleBeam()) {
    1150         346 :       beam = psfInfo.restoringBeam();
    1151             :     }
    1152         226 :     else if (psfInfo.hasMultipleBeams() && itsUseBeam=="common") {
    1153           9 :       beam = CasaImageBeamSet(psfInfo.getBeamSet()).getCommonBeam();
    1154             :     }
    1155         572 :     Double beammaj = beam.getMajor(Unit("arcsec"));
    1156         572 :     Double beammin = beam.getMinor(Unit("arcsec"));
    1157         572 :     if (std::isinf(beammaj) || std::isinf(beammin)) {
    1158           0 :       String msg;
    1159           0 :       if (itsUseBeam=="common") {
    1160           0 :         msg+="Bad restoring beam using the common beam (at least one of the beam axes has an infinite size)  ";
    1161           0 :         throw(AipsError(msg));
    1162             :       }
    1163           0 :     }
    1164         572 :     itsImages->releaseLocks();
    1165         572 :   }
    1166             : 
    1167             :   // This is for interactive-clean.
    1168           0 :   void SynthesisDeconvolver::getCopyOfResidualAndMask( TempImage<Float> &/*residual*/,
    1169             :                                            TempImage<Float> &/*mask*/ )
    1170             :   {
    1171             :     // Actually all I think we need here are filenames JSK 12/12/12
    1172             :     // resize/shape and copy the residual image and mask image to these in/out variables.
    1173             :     // Allocate Memory here.
    1174           0 :   }
    1175           0 :   void SynthesisDeconvolver::setMask( TempImage<Float> &/*mask*/ )
    1176             :   {
    1177             :     // Here we will just pass in the new names
    1178             :     // Copy the input mask to the local main image mask
    1179           0 :   }
    1180          15 :   void SynthesisDeconvolver::setIterDone(const Int iterdone){
    1181             :     //cerr << "SETITERDONE iterdone " << iterdone << endl;
    1182             :     ///this get lost in initMinorCycle
    1183             :     //itsLoopController.incrementMinorCycleCount(iterdone);
    1184          15 :     itsIterDone=iterdone;
    1185             : 
    1186          15 :   }
    1187           0 :   void SynthesisDeconvolver::setPosMask(std::shared_ptr<ImageInterface<Float> > posmask){
    1188           0 :     itsPosMask=posmask;
    1189             : 
    1190           0 :   }
    1191             : 
    1192         409 :   auto key2Mat = [](const Record& rec, const String& key, const Int npol) {
    1193         409 :      Matrix<Double> tmp;
    1194             :      //cerr << "KEY2mat " << key <<"  "<< rec.asArrayDouble(key).shape() << endl;
    1195         409 :      if(rec.asArrayDouble(key).shape().nelements()==1){
    1196         274 :        if(rec.asArrayDouble(key).shape()[0] != npol){
    1197         274 :          tmp.resize(1,rec.asArrayDouble(key).shape()[0]);
    1198         274 :          Vector<Double>tmpvec=rec.asArrayDouble(key);
    1199         274 :          tmp.row(0)=tmpvec;
    1200         274 :        }
    1201             :        else{
    1202           0 :          tmp.resize(rec.asArrayDouble(key).shape()[0],1);
    1203           0 :          Vector<Double>tmpvec=rec.asArrayDouble(key);
    1204           0 :          tmp.column(0)=tmpvec;
    1205           0 :        }
    1206             : 
    1207             :      }
    1208             :      else{
    1209         135 :        tmp=rec.asArrayDouble(key);
    1210             :      }
    1211         409 :      return tmp;
    1212           0 :    };
    1213             : 
    1214         263 :   Record SynthesisDeconvolver::getSubsetRobustStats(const Int chanBeg, const Int chanEnd){
    1215         263 :     Record outRec;
    1216             :     //cerr << "getSUB " << itsRobustStats << endl;
    1217         263 :     if(itsRobustStats.nfields()==0)
    1218         214 :       return outRec;
    1219          49 :     Matrix<Double> tmp;
    1220         441 :     std::vector<String> keys={"min", "max", "rms", "medabsdevmed", "med", "robustrms", "median"};
    1221         392 :     for (auto it = keys.begin() ; it != keys.end(); ++it){
    1222         343 :       if(itsRobustStats.isDefined(*it)){
    1223         128 :         tmp.resize();
    1224         128 :         tmp=key2Mat(itsRobustStats, *it, itsImages->residual()->shape()[2]);
    1225             :         /*
    1226             :         cerr << "size of " << *it << "   " << itsRobustStats.asArrayDouble(*it).shape() << endl;
    1227             :         if(itsRobustStats.asArrayDouble(*it).shape().nelements()==1){
    1228             :           tmp.resize(1, itsRobustStats.asArrayDouble(*it).shape()[0]);
    1229             :           Vector<Double>tmpvec=itsRobustStats.asArrayDouble(*it);
    1230             :           tmp.row(0)=tmpvec;
    1231             : 
    1232             :         }
    1233             :         else
    1234             :           tmp=itsRobustStats.asArrayDouble(*it);
    1235             :         */
    1236             :         //      cerr << std::setprecision(12) << tmp[chanBeg] << " bool " <<(tmp[chanBeg]> (C::dbl_max-(C::dbl_max*1e-15))) << endl;
    1237         128 :         if(tmp(0,chanBeg)> (C::dbl_max-(C::dbl_max*1e-15)))
    1238           0 :           return Record();
    1239             :         //cerr << "GETSUB blc "<< IPosition(2, 0, chanBeg)<<  " trc " << IPosition(2, tmp.shape()[0]-1, chanEnd) << " shape " << tmp.shape() << endl;
    1240         128 :         outRec.define(*it, tmp(IPosition(2, 0, chanBeg), IPosition(2, tmp.shape()[0]-1, chanEnd)));
    1241             :       }
    1242             :     }
    1243             : 
    1244             :     // When itsImages->residual() is called, it (sometimes) creates a read-lock. Release that lock.
    1245          49 :     shared_ptr<ImageInterface<Float> > resimg = itsImages->residual();
    1246          49 :     itsImages->releaseImage( resimg );
    1247             : 
    1248             :     //cerr <<"chanbeg " << chanBeg << " chanend " << chanEnd << endl;
    1249             :     //cerr << "GETSUB " << outRec << endl;
    1250          49 :     return outRec;
    1251         263 :   }
    1252             : 
    1253         263 :   void SynthesisDeconvolver::setSubsetRobustStats(const Record& inrec, const Int chanBeg, const Int chanEnd, const Int numchan){
    1254         263 :     if(inrec.nfields()==0)
    1255         209 :       return ;
    1256          54 :     Matrix<Double> tmp;
    1257         486 :     std::vector<String> keys={"min", "max", "rms", "medabsdevmed", "med", "robustrms", "median"};
    1258             : 
    1259         432 :     for (auto it = keys.begin() ; it != keys.end(); ++it){
    1260         378 :       if(inrec.isDefined(*it)){
    1261         153 :         tmp.resize();
    1262         153 :         tmp=key2Mat(inrec, *it,itsImages->residual()->shape()[2] );
    1263         153 :         Matrix<Double> outvec;
    1264         153 :         if(itsRobustStats.isDefined(*it)){
    1265         128 :           outvec=key2Mat(itsRobustStats, *it, itsImages->residual()->shape()[2]);
    1266             :         }
    1267             :         else{
    1268          25 :           outvec.resize(itsImages->residual()->shape()[2], numchan);
    1269          25 :           outvec.set(C::dbl_max);
    1270             :         }
    1271             : 
    1272             : 
    1273         153 :         outvec(IPosition(2, 0, chanBeg), IPosition(2,outvec.shape()[0]-1, chanEnd))=tmp;
    1274         153 :         itsRobustStats.define(*it, outvec);
    1275         153 :       }
    1276             :     }
    1277             : 
    1278             :     // When itsImages->residual() is called, it (sometimes) creates a read-lock. Release that lock.
    1279          54 :     shared_ptr<ImageInterface<Float> > resimg = itsImages->residual();
    1280          54 :     itsImages->releaseImage( resimg );
    1281             : 
    1282             :     //cerr << "SETT " << itsRobustStats << endl;
    1283             :     //cerr << "SETT::ItsRobustStats " << Vector<Double>(itsRobustStats.asArrayDouble("min")) << endl;
    1284          54 :   }
    1285             : } //# NAMESPACE CASA - END
    1286             : 

Generated by: LCOV version 1.16