LCOV - code coverage report
Current view: top level - synthesis/ImagerObjects - SIMinorCycleController.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 187 0.0 %
Date: 2024-10-28 15:53:10 Functions: 0 32 0.0 %

          Line data    Source code
       1             : 
       2             : 
       3             : //# SIISubterBot.cc: This file contains the implementation of the SISubIterBot class.
       4             : //#
       5             : //#  CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
       6             : //#  Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
       7             : //#  Copyright (C) European Southern Observatory, 2011, All rights reserved.
       8             : //#
       9             : //#  This library is free software; you can redistribute it and/or
      10             : //#  modify it under the terms of the GNU Lesser General Public
      11             : //#  License as published by the Free software Foundation; either
      12             : //#  version 2.1 of the License, or (at your option) any later version.
      13             : //#
      14             : //#  This library is distributed in the hope that it will be useful,
      15             : //#  but WITHOUT ANY WARRANTY, without even the implied warranty of
      16             : //#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      17             : //#  Lesser General Public License for more details.
      18             : //#
      19             : //#  You should have received a copy of the GNU Lesser General Public
      20             : //#  License along with this library; if not, write to the Free Software
      21             : //#  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
      22             : //#  MA 02111-1307  USA
      23             : //# $Id: $
      24             : 
      25             : #include <synthesis/ImagerObjects/SIMinorCycleController.h>
      26             : 
      27             : /* Records Interface */
      28             : #include <casacore/casa/Containers/Record.h>
      29             : #include <casacore/casa/BasicMath/Math.h>
      30             : #include <casacore/casa/Utilities/GenSort.h>
      31             : #include <math.h> // For FLT_MAX
      32             : 
      33             : using namespace casacore;
      34             : namespace casa { //# NAMESPACE CASA - BEGIN
      35             :   
      36             :   
      37           0 :   SIMinorCycleController::SIMinorCycleController(): 
      38           0 :                                 itsCycleNiter(0),
      39           0 :                                 itsCycleThreshold(0.0),
      40           0 :                                 itsNsigmaThreshold(0.0),
      41           0 :                                 itsLoopGain(0.1),
      42           0 :                                 itsIsThresholdReached(false),
      43           0 :                                 itsUpdatedModelFlag(false),
      44           0 :                                 itsIterDone(0),
      45           0 :                                 itsCycleIterDone(0),
      46           0 :                                 itsTotalIterDone(0),
      47           0 :                                 itsMaxCycleIterDone(0),
      48           0 :                                 itsPeakResidual(0),
      49           0 :                                 itsIntegratedFlux(0),
      50           0 :                                 itsMaxPsfSidelobe(0),
      51           0 :                                 itsMinResidual(0),itsMinResidualNoMask(0),
      52           0 :                                 itsPeakResidualNoMask(0), itsNsigma(0),
      53           0 :                                 itsMadRMS(0), itsMaskSum(0),
      54             :                                 //itsSummaryMinor(IPosition(2,
      55             :                                 //                            SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
      56           0 :                                 itsSummaryMinor(IPosition(2, SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
      57             :                                                             0)),
      58           0 :                                 itsDeconvolverID(0) 
      59           0 :   {}
      60             : 
      61             : 
      62             : 
      63           0 :   SIMinorCycleController::~SIMinorCycleController(){}
      64             : 
      65             : 
      66           0 :   Int SIMinorCycleController::majorCycleRequired(Float currentPeakResidual)
      67             :   {
      68           0 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
      69             :     
      70           0 :     Int stopCode=0;
      71             : 
      72             :     // Reached iteration limit
      73           0 :     if (itsCycleIterDone >= itsCycleNiter ) {stopCode=1;}
      74             :     // Reached cyclethreshold
      75             :     //if( fabs(currentPeakResidual) <= itsCycleThreshold ) { stopCode=2; }
      76             :     // Reached cyclethreshold or n-sigma threshold
      77             :     //debug (TT)
      78             :     //os << LogIO::DEBUG1<< "itsNsigma="<<itsNsigma<<" itsIterDiff="<<itsIterDiff<<LogIO::POST;
      79           0 :     os << LogIO::DEBUG1<< "itsNsigmaThreshoild="<<itsNsigmaThreshold<<" itsCycleThreshold="<<itsCycleThreshold<<" currentPeakRes="<<currentPeakResidual<<LogIO::POST;
      80           0 :     if (itsCycleThreshold >= itsNsigmaThreshold) {
      81             :       //if( fabs(currentPeakResidual) <= itsCycleThreshold ) { stopCode=2; }
      82           0 :       if( fabs(currentPeakResidual) <= itsCycleThreshold ) { 
      83             :         //itsNsigmaThreshold = 0.0; // since stopped by gobal threshold, reset itsNsigmaThreshold 
      84           0 :         stopCode=2; 
      85             :       }
      86             :     }
      87             :     else {
      88           0 :       if( fabs(currentPeakResidual) <= itsNsigmaThreshold && !(itsIterDiff<=0)) { if (itsNsigma!=0.0) stopCode=6; }
      89             :     }
      90             :     // Zero iterations done
      91           0 :     if( itsIterDiff==0 ) {stopCode=3;}
      92             :     // Diverged : CAS-8767, CAS-8584
      93             :     //cout << " itsIterDiff : " << itsIterDiff << "  itsPeak : " << itsPeakResidual << " currentPeak : " << currentPeakResidual << " itsMin : " << itsMinResidual << " stopcode so far : " << stopCode ;
      94           0 :     if( itsIterDiff>0 &&
      95           0 :         ( fabs(itsMinResidual) > 0.0 ) && 
      96           0 :         ( fabs(currentPeakResidual) - fabs(itsMinResidual) )/ fabs(itsMinResidual) >0.1  ) 
      97           0 :       {stopCode=4;}
      98             : 
      99             :    //  cout << " -> " << stopCode << endl;
     100             : 
     101             :     /*    // Going nowhere
     102             :     if( itsIterDiff > 1500 && 
     103             :         ( fabs(itsPeakResidual) > 0.0 ) && 
     104             :         ( fabs(currentPeakResidual) - fabs(itsPeakResidual) )/ fabs(itsPeakResidual) < itsLoopGain  ) 
     105             :       {stopCode=5;}
     106             :     */
     107             : 
     108           0 :     return stopCode;
     109             : 
     110             : 
     111           0 :   }
     112             : 
     113             : 
     114           0 :   Float SIMinorCycleController::getLoopGain()
     115             :   {
     116           0 :     return itsLoopGain;
     117             :   }
     118             :   
     119             :   
     120           0 :   void SIMinorCycleController::setUpdatedModelFlag(Bool updatedmodel)
     121             :   {
     122           0 :     itsUpdatedModelFlag = updatedmodel;
     123           0 :   }
     124             : 
     125           0 :   void SIMinorCycleController::incrementMinorCycleCount(Int itersDonePerStep)
     126             :   {
     127             :     /*
     128             :     if( itersDonePerStep <= 0 )
     129             :       {
     130             :         LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     131             :         os << LogIO::WARN << "Zero iterations done after " << itsCycleIterDone << LogIO::POST;
     132             :       }
     133             :     */
     134             : 
     135           0 :     itsIterDiff = itersDonePerStep;
     136           0 :     itsIterDone += itersDonePerStep;
     137           0 :     itsTotalIterDone += itersDonePerStep;
     138           0 :     itsCycleIterDone += itersDonePerStep;
     139           0 :   }
     140             : 
     141           0 :   Float SIMinorCycleController::getPeakResidual()
     142             :   {
     143           0 :     return itsPeakResidual;
     144             :   }
     145             : 
     146             :   // This is the max residual across all channels/stokes (subimages).
     147             :   // This is returned in the end-minor-cycle record.
     148             :   /// TODO : Make arrays/lists for peakresidual per 'subimage'. Max over subims gets returned.
     149           0 :   void SIMinorCycleController::setPeakResidual(Float peakResidual)
     150             :   {
     151           0 :     itsPeakResidual = peakResidual;
     152             :     //    cout << "Setting peak res (SIMinorCycleController) : " << itsPeakResidual << endl;
     153             : 
     154           0 :     if( itsMinResidual > itsPeakResidual )
     155           0 :       itsMinResidual = itsPeakResidual;
     156             : 
     157           0 :   }
     158             : 
     159           0 :   void SIMinorCycleController::setPeakResidualNoMask(Float peakResidual)
     160             :   {
     161           0 :     itsPeakResidualNoMask = peakResidual;
     162             :     //    cout << "Setting peak res (SIMinorCycleController) : " << itsPeakResidual << endl;
     163             : 
     164           0 :     if( itsMinResidualNoMask > itsPeakResidualNoMask )
     165           0 :       itsMinResidualNoMask = itsPeakResidualNoMask;
     166             : 
     167           0 :   }
     168             : 
     169           0 :   void SIMinorCycleController::setMadRMS(Float madRMS)
     170             :   {
     171           0 :     itsMadRMS = madRMS;
     172           0 :   }
     173             : 
     174           0 :   Float SIMinorCycleController::getNsigma()
     175             :   {
     176           0 :     return itsNsigma;
     177             :   }
     178             : 
     179           0 :   void SIMinorCycleController::setNsigma(Float nSigma)
     180             :   {
     181           0 :     itsNsigma = nSigma;
     182           0 :   }
     183             : 
     184           0 :   void SIMinorCycleController::setNsigmaThreshold(Float nsigmaThreshold)
     185             :   {
     186             :     
     187           0 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     188           0 :     itsNsigmaThreshold = nsigmaThreshold;
     189           0 :   }
     190             : 
     191           0 :   void SIMinorCycleController::setPBMask(Float pbMaskLevel)
     192             :   {
     193           0 :     itsPBMaskLevel = pbMaskLevel;
     194           0 :   }
     195             : 
     196           0 :   void SIMinorCycleController::setMaskSum(Float maskSum)
     197             :   {
     198           0 :     itsMaskSum = maskSum;
     199           0 :   }
     200             : 
     201           0 :   void SIMinorCycleController::setFullSummary(bool fullSummary)
     202             :   {
     203           0 :     itsFullSummary = fullSummary;
     204           0 :   }
     205             : 
     206           0 :   void SIMinorCycleController::resetMinResidual()
     207             :   {
     208           0 :     itsMinResidual = itsPeakResidual;
     209           0 :     itsIterDiff=-1;
     210           0 :   }
     211             : 
     212           0 :   Float SIMinorCycleController::getIntegratedFlux()
     213             :   {
     214           0 :     return itsIntegratedFlux;
     215             :   }
     216             : 
     217           0 :   void SIMinorCycleController::addIntegratedFlux(Float integratedFlux)
     218             :   {
     219           0 :     itsIntegratedFlux += integratedFlux;
     220           0 :   }
     221             : 
     222           0 :   Float SIMinorCycleController::getMaxPsfSidelobe()
     223             :   {
     224           0 :     return itsMaxPsfSidelobe;
     225             :   }
     226             : 
     227           0 :   void SIMinorCycleController::setMaxPsfSidelobe(Float maxPsfSidelobe)
     228             :   {
     229           0 :     itsMaxPsfSidelobe = maxPsfSidelobe;
     230           0 :   }
     231             : 
     232           0 :   Int SIMinorCycleController::getIterDone()
     233             :   {
     234           0 :     return itsIterDone;
     235             :   }
     236             : 
     237           0 :   Int SIMinorCycleController::getCycleNiter()
     238             :   {
     239           0 :     return itsCycleNiter;
     240             :   }
     241             : 
     242           0 :   Float SIMinorCycleController::getCycleThreshold()
     243             :   {
     244           0 :     return itsCycleThreshold;
     245             :   }
     246             :  
     247           0 :   Bool SIMinorCycleController::isThresholdReached()
     248             :   {
     249           0 :     return itsIsThresholdReached;
     250             :   }
     251             : 
     252           0 :   Record SIMinorCycleController::getCycleExecutionRecord() {
     253           0 :     LogIO os( LogOrigin("SISkyModel",__FUNCTION__,WHERE) );
     254           0 :     Record returnRecord;
     255             : 
     256           0 :     returnRecord.define( RecordFieldId("iterdone"),  itsIterDone);
     257           0 :     returnRecord.define(RecordFieldId("peakresidual"), itsPeakResidual);
     258           0 :     returnRecord.define(RecordFieldId("updatedmodelflag"), 
     259           0 :                         itsUpdatedModelFlag);
     260           0 :     returnRecord.define(RecordFieldId("summaryminor"), itsSummaryMinor);
     261           0 :     returnRecord.define(RecordFieldId("updatedmodelflag"),
     262           0 :                         itsUpdatedModelFlag);
     263           0 :     returnRecord.define(RecordFieldId("maxcycleiterdone"),
     264             :                         itsMaxCycleIterDone);
     265           0 :     returnRecord.define( RecordFieldId("peakresidualnomask"), itsPeakResidualNoMask);
     266             : 
     267           0 :     return returnRecord;
     268           0 :   }
     269             : 
     270           0 :   Float SIMinorCycleController::getPBMask() {
     271           0 :     return itsPBMaskLevel;
     272             :   }
     273             : 
     274           0 :   Record SIMinorCycleController::getCycleInitializationRecord() {
     275           0 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     276             : 
     277           0 :     Record returnRecord;
     278             : 
     279             :     /* Control Variables */
     280           0 :     returnRecord.define(RecordFieldId("peakresidual"), itsPeakResidual);
     281           0 :     returnRecord.define(RecordFieldId("maxpsfsidelobe"), itsMaxPsfSidelobe);
     282           0 :     returnRecord.define( RecordFieldId("peakresidualnomask"), itsPeakResidualNoMask);
     283           0 :     returnRecord.define( RecordFieldId("madrms"), itsMadRMS);
     284           0 :     returnRecord.define( RecordFieldId("masksum"), itsMaskSum);
     285           0 :     returnRecord.define( RecordFieldId("nsigmathreshold"), itsNsigmaThreshold);
     286           0 :     returnRecord.define( RecordFieldId("nsigma"), itsNsigma);
     287           0 :     returnRecord.define( RecordFieldId("fullsummary"), itsFullSummary);
     288             : 
     289             :     /* Reset Counters and summary for the current set of minorcycle iterations */
     290           0 :     itsIterDone = 0;
     291           0 :     itsIterDiff = -1;
     292             :     //int nSummaryFields = SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     293           0 :     int nSummaryFields = !itsFullSummary ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     294           0 :     itsSummaryMinor.resize( IPosition( 2, nSummaryFields, 0) , true );
     295             : 
     296           0 :     return returnRecord;
     297           0 :   }
     298             : 
     299           0 :   void SIMinorCycleController::setCycleControls(Record &recordIn) {
     300           0 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     301             : 
     302           0 :     if (recordIn.isDefined("cycleniter"))
     303           0 :       {recordIn.get(RecordFieldId("cycleniter"), itsCycleNiter);}
     304             :     else
     305           0 :       {throw(AipsError("cycleniter not defined in input minor-cycle controller") );}
     306             : 
     307           0 :     if (recordIn.isDefined("cyclethreshold")) 
     308           0 :       {recordIn.get(RecordFieldId("cyclethreshold"),itsCycleThreshold);}
     309             :     else
     310           0 :       {throw(AipsError("cyclethreshold not defined in input minor-cycle controller") );}
     311             :      
     312           0 :     if (recordIn.isDefined("thresholdreached"))
     313           0 :       {recordIn.get(RecordFieldId("thresholdreached"), itsIsThresholdReached);}
     314             :     else
     315           0 :       { throw(AipsError("thresholdreached not defined in input minor-cycle controller") );}
     316             : 
     317           0 :     if (recordIn.isDefined("loopgain")) 
     318           0 :       {recordIn.get(RecordFieldId("loopgain"), itsLoopGain);}
     319             :     else
     320           0 :       {throw(AipsError("loopgain not defined in input minor-cycle controller") );}
     321             : 
     322           0 :     if (recordIn.isDefined("nsigma"))
     323           0 :       {recordIn.get(RecordFieldId("nsigma"), itsNsigma);}
     324             :     else 
     325           0 :       { throw(AipsError(" nsigma is not defined in input minor-cycle controller ") );}
     326             : 
     327             :     //if (recordIn.isDefined("fullsummary"))
     328             :    //  {recordIn.get(RecordFieldId("fullsummary"), itsFullSummary);}
     329             :     //else 
     330             :     //  { throw(AipsError(" fullsummary is not defined in input minor-cycle controller ") );}
     331             :     //os<<" in setCycleControls done... "<<LogIO::POST;
     332             : 
     333             :     /* Reset the counters for the new cycle */
     334           0 :     itsMaxCycleIterDone = 0;
     335           0 :     itsCycleIterDone = 0;
     336           0 :     itsUpdatedModelFlag = false;
     337           0 :   }
     338             : 
     339           0 :   void SIMinorCycleController::resetCycleIter(){
     340           0 :     itsMaxCycleIterDone = max(itsCycleIterDone, itsMaxCycleIterDone);
     341           0 :     itsCycleIterDone = 0;
     342           0 :   }
     343             : 
     344           0 :   void SIMinorCycleController::addSummaryMinor(uInt deconvolverid, uInt chan, uInt pol,
     345             :                                                Int cycleStartIter, Int startIterDone, Float startmodelflux, Float startpeakresidual, Float startpeakresidualnomask,
     346             :                                                Float modelflux, Float peakresidual, Float peakresidualnomask, Float masksum, Int mpiRank, Int stopCode, bool fullsummary)
     347             :   {
     348           0 :     LogIO os( LogOrigin("SIMinorCycleController", __FUNCTION__ ,WHERE) );
     349           0 :     IPosition shp = itsSummaryMinor.shape();
     350             :     //bool uss = SIMinorCycleController::useSmallSummaryminor(); // temporary CAS-13683 workaround
     351             :     //int nSummaryFields = uss ? 6 : SIMinorCycleController::nSummaryFields;
     352             :     //int nSummaryFields = fullsummary ? 6 : SIMinorCycleController::nSummaryFields;
     353             : 
     354           0 :     int nSummaryFields = !fullsummary ? 6 : SIMinorCycleController::nSummaryFields;
     355           0 :     if( shp.nelements() != 2 && shp[0] != nSummaryFields ) 
     356           0 :       throw(AipsError("Internal error in shape of minor-cycle summary record"));
     357           0 :      itsSummaryMinor.resize( IPosition( 2, nSummaryFields, shp[1]+1 ) , true );
     358             :      // iterations done
     359             :      // make it non-cummulative for all cases
     360           0 :      itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone - startIterDone;
     361             :      //if(!fullsummary) {
     362             :      //    itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone - startIterDone;
     363             :      //}
     364             :      //else {
     365             :      //    itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone;
     366             :      // }
     367             :      // peak residual
     368           0 :      itsSummaryMinor( IPosition(2, 1, shp[1] ) ) = (Double) peakresidual;
     369             :      // model flux
     370           0 :      itsSummaryMinor( IPosition(2, 2, shp[1] ) ) = (Double) modelflux;
     371             :      // cycle threshold
     372           0 :      itsSummaryMinor( IPosition(2, 3, shp[1] ) ) = itsCycleThreshold;
     373             :      // mapper id (or multifield id temporary CAS-13683 workaround)
     374           0 :      itsSummaryMinor( IPosition(2, 4, shp[1] ) ) = deconvolverid;
     375             :      // channel id
     376           0 :      itsSummaryMinor( IPosition(2, 5, shp[1] ) ) = chan;
     377             :      //if (!uss) {
     378           0 :      if (fullsummary) {
     379             :          // polarity id
     380           0 :          itsSummaryMinor( IPosition(2, 6, shp[1] ) ) = pol;
     381             :          // cycle start iterations done (ie earliest iterDone for the entire minor cycle)
     382           0 :          itsSummaryMinor( IPosition(2, 7, shp[1] ) ) = cycleStartIter;
     383             :          // starting iterations done
     384           0 :          itsSummaryMinor( IPosition(2, 8, shp[1] ) ) = startIterDone;
     385             :          // starting peak residual
     386           0 :          itsSummaryMinor( IPosition(2, 9, shp[1] ) ) = (Double) startpeakresidual;
     387             :          // starting model flux
     388           0 :          itsSummaryMinor( IPosition(2, 10, shp[1] ) ) = (Double) startmodelflux;
     389             :          // starting peak residual, not limited to the user's mask
     390           0 :          itsSummaryMinor( IPosition(2, 11, shp[1] ) ) = (Double) startpeakresidualnomask;
     391             :          // peak residual, not limited to the user's mask
     392           0 :          itsSummaryMinor( IPosition(2, 12, shp[1] ) ) = (Double) peakresidualnomask;
     393             :          // number of pixels in the mask
     394           0 :          itsSummaryMinor( IPosition(2, 13, shp[1] ) ) = (Double) masksum;
     395             :          // mpi server
     396           0 :          itsSummaryMinor( IPosition(2, 14, shp[1] ) ) = mpiRank;
     397             :          // outlier field id, to be provided in grpcInteractiveCleanManager::mergeMinorCycleSummary
     398           0 :          itsSummaryMinor( IPosition(2, 15, shp[1] ) ) = 0;
     399             :          // stopcode
     400           0 :          itsSummaryMinor( IPosition(2, 16, shp[1] ) ) = stopCode;
     401             :      }
     402             : 
     403           0 :   }// end of addSummaryMinor
     404             : 
     405             :   // temporary CAS-13683 workaround
     406             :   /***
     407             :   Bool SIMinorCycleController::useSmallSummaryminor()
     408             :   {
     409             :     if (const char* use_small_summaryminor_p = std::getenv("USE_SMALL_SUMMARYMINOR"))
     410             :     {
     411             :         string use_small_summaryminor(use_small_summaryminor_p);
     412             :         if (use_small_summaryminor.compare("TRUE") == 0 || use_small_summaryminor.compare("true") == 0) {
     413             :             return true;
     414             :         }
     415             :     }
     416             :     return false;
     417             :   }
     418             :   ***/
     419             :   
     420             : } //# NAMESPACE CASA - END
     421             : 

Generated by: LCOV version 1.16