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        1047 :   SIMinorCycleController::SIMinorCycleController(): 
      38             :                                 itsCycleNiter(0),
      39             :                                 itsCycleThreshold(0.0),
      40             :                                 itsNsigmaThreshold(0.0),
      41             :                                 itsLoopGain(0.1),
      42             :                                 itsIsThresholdReached(false),
      43             :                                 itsUpdatedModelFlag(false),
      44             :                                 itsIterDone(0),
      45             :                                 itsCycleIterDone(0),
      46             :                                 itsTotalIterDone(0),
      47             :                                 itsMaxCycleIterDone(0),
      48             :                                 itsPeakResidual(0),
      49             :                                 itsIntegratedFlux(0),
      50             :                                 itsMaxPsfSidelobe(0),
      51             :                                 itsMinResidual(0),itsMinResidualNoMask(0),
      52             :                                 itsPeakResidualNoMask(0), itsNsigma(0),
      53             :                                 itsMadRMS(0), itsMaskSum(0),
      54             :                                 //itsSummaryMinor(IPosition(2,
      55             :                                 //                            SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
      56        1047 :                                 itsSummaryMinor(IPosition(2, SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
      57             :                                                             0)),
      58        1047 :                                 itsDeconvolverID(0) 
      59        1047 :   {}
      60             : 
      61             : 
      62             : 
      63        1047 :   SIMinorCycleController::~SIMinorCycleController(){}
      64             : 
      65             : 
      66        5968 :   Int SIMinorCycleController::majorCycleRequired(Float currentPeakResidual)
      67             :   {
      68       11936 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
      69             :     
      70        5968 :     Int stopCode=0;
      71             : 
      72             :     // Reached iteration limit
      73        5968 :     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        5968 :     os << LogIO::DEBUG1<< "itsNsigmaThreshoild="<<itsNsigmaThreshold<<" itsCycleThreshold="<<itsCycleThreshold<<" currentPeakRes="<<currentPeakResidual<<LogIO::POST;
      80        5968 :     if (itsCycleThreshold >= itsNsigmaThreshold) {
      81             :       //if( fabs(currentPeakResidual) <= itsCycleThreshold ) { stopCode=2; }
      82        5895 :       if( fabs(currentPeakResidual) <= itsCycleThreshold ) { 
      83             :         //itsNsigmaThreshold = 0.0; // since stopped by gobal threshold, reset itsNsigmaThreshold 
      84        1103 :         stopCode=2; 
      85             :       }
      86             :     }
      87             :     else {
      88          73 :       if( fabs(currentPeakResidual) <= itsNsigmaThreshold && !(itsIterDiff<=0)) { if (itsNsigma!=0.0) stopCode=6; }
      89             :     }
      90             :     // Zero iterations done
      91        5968 :     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        2821 :     if( itsIterDiff>0 &&
      95        8789 :         ( fabs(itsMinResidual) > 0.0 ) && 
      96        2821 :         ( fabs(currentPeakResidual) - fabs(itsMinResidual) )/ fabs(itsMinResidual) >0.1  ) 
      97           4 :       {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       11936 :     return stopCode;
     109             : 
     110             : 
     111             :   }
     112             : 
     113             : 
     114        3741 :   Float SIMinorCycleController::getLoopGain()
     115             :   {
     116        3741 :     return itsLoopGain;
     117             :   }
     118             :   
     119             :   
     120        3141 :   void SIMinorCycleController::setUpdatedModelFlag(Bool updatedmodel)
     121             :   {
     122        3141 :     itsUpdatedModelFlag = updatedmodel;
     123        3141 :   }
     124             : 
     125        3075 :   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        3075 :     itsIterDiff = itersDonePerStep;
     136        3075 :     itsIterDone += itersDonePerStep;
     137        3075 :     itsTotalIterDone += itersDonePerStep;
     138        3075 :     itsCycleIterDone += itersDonePerStep;
     139        3075 :   }
     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        9311 :   void SIMinorCycleController::setPeakResidual(Float peakResidual)
     150             :   {
     151        9311 :     itsPeakResidual = peakResidual;
     152             :     //    cout << "Setting peak res (SIMinorCycleController) : " << itsPeakResidual << endl;
     153             : 
     154        9311 :     if( itsMinResidual > itsPeakResidual )
     155        2870 :       itsMinResidual = itsPeakResidual;
     156             : 
     157        9311 :   }
     158             : 
     159        2431 :   void SIMinorCycleController::setPeakResidualNoMask(Float peakResidual)
     160             :   {
     161        2431 :     itsPeakResidualNoMask = peakResidual;
     162             :     //    cout << "Setting peak res (SIMinorCycleController) : " << itsPeakResidual << endl;
     163             : 
     164        2431 :     if( itsMinResidualNoMask > itsPeakResidualNoMask )
     165           0 :       itsMinResidualNoMask = itsPeakResidualNoMask;
     166             : 
     167        2431 :   }
     168             : 
     169           0 :   void SIMinorCycleController::setMadRMS(Float madRMS)
     170             :   {
     171           0 :     itsMadRMS = madRMS;
     172           0 :   }
     173             : 
     174        3147 :   Float SIMinorCycleController::getNsigma()
     175             :   {
     176        3147 :     return itsNsigma;
     177             :   }
     178             : 
     179           0 :   void SIMinorCycleController::setNsigma(Float nSigma)
     180             :   {
     181           0 :     itsNsigma = nSigma;
     182           0 :   }
     183             : 
     184        2833 :   void SIMinorCycleController::setNsigmaThreshold(Float nsigmaThreshold)
     185             :   {
     186             :     
     187        5666 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     188        2833 :     itsNsigmaThreshold = nsigmaThreshold;
     189        2833 :   }
     190             : 
     191        2431 :   void SIMinorCycleController::setPBMask(Float pbMaskLevel)
     192             :   {
     193        2431 :     itsPBMaskLevel = pbMaskLevel;
     194        2431 :   }
     195             : 
     196        2431 :   void SIMinorCycleController::setMaskSum(Float maskSum)
     197             :   {
     198        2431 :     itsMaskSum = maskSum;
     199        2431 :   }
     200             : 
     201        2431 :   void SIMinorCycleController::setFullSummary(bool fullSummary)
     202             :   {
     203        2431 :     itsFullSummary = fullSummary;
     204        2431 :   }
     205             : 
     206        3147 :   void SIMinorCycleController::resetMinResidual()
     207             :   {
     208        3147 :     itsMinResidual = itsPeakResidual;
     209        3147 :     itsIterDiff=-1;
     210        3147 :   }
     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        2431 :   void SIMinorCycleController::setMaxPsfSidelobe(Float maxPsfSidelobe)
     228             :   {
     229        2431 :     itsMaxPsfSidelobe = maxPsfSidelobe;
     230        2431 :   }
     231             : 
     232       12834 :   Int SIMinorCycleController::getIterDone()
     233             :   {
     234       12834 :     return itsIterDone;
     235             :   }
     236             : 
     237        6362 :   Int SIMinorCycleController::getCycleNiter()
     238             :   {
     239        6362 :     return itsCycleNiter;
     240             :   }
     241             : 
     242       10732 :   Float SIMinorCycleController::getCycleThreshold()
     243             :   {
     244       10732 :     return itsCycleThreshold;
     245             :   }
     246             :  
     247          52 :   Bool SIMinorCycleController::isThresholdReached()
     248             :   {
     249          52 :     return itsIsThresholdReached;
     250             :   }
     251             : 
     252         912 :   Record SIMinorCycleController::getCycleExecutionRecord() {
     253        2736 :     LogIO os( LogOrigin("SISkyModel",__FUNCTION__,WHERE) );
     254         912 :     Record returnRecord;
     255             : 
     256         912 :     returnRecord.define( RecordFieldId("iterdone"),  itsIterDone);
     257         912 :     returnRecord.define(RecordFieldId("peakresidual"), itsPeakResidual);
     258         912 :     returnRecord.define(RecordFieldId("updatedmodelflag"), 
     259         912 :                         itsUpdatedModelFlag);
     260         912 :     returnRecord.define(RecordFieldId("summaryminor"), itsSummaryMinor);
     261         912 :     returnRecord.define(RecordFieldId("updatedmodelflag"),
     262         912 :                         itsUpdatedModelFlag);
     263         912 :     returnRecord.define(RecordFieldId("maxcycleiterdone"),
     264             :                         itsMaxCycleIterDone);
     265         912 :     returnRecord.define( RecordFieldId("peakresidualnomask"), itsPeakResidualNoMask);
     266             : 
     267        1824 :     return returnRecord;
     268             :   }
     269             : 
     270         920 :   Float SIMinorCycleController::getPBMask() {
     271         920 :     return itsPBMaskLevel;
     272             :   }
     273             : 
     274        2431 :   Record SIMinorCycleController::getCycleInitializationRecord() {
     275        7293 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     276             : 
     277        2431 :     Record returnRecord;
     278             : 
     279             :     /* Control Variables */
     280        2431 :     returnRecord.define(RecordFieldId("peakresidual"), itsPeakResidual);
     281        2431 :     returnRecord.define(RecordFieldId("maxpsfsidelobe"), itsMaxPsfSidelobe);
     282        2431 :     returnRecord.define( RecordFieldId("peakresidualnomask"), itsPeakResidualNoMask);
     283        2431 :     returnRecord.define( RecordFieldId("madrms"), itsMadRMS);
     284        2431 :     returnRecord.define( RecordFieldId("masksum"), itsMaskSum);
     285        2431 :     returnRecord.define( RecordFieldId("nsigmathreshold"), itsNsigmaThreshold);
     286        2431 :     returnRecord.define( RecordFieldId("nsigma"), itsNsigma);
     287        2431 :     returnRecord.define( RecordFieldId("fullsummary"), itsFullSummary);
     288             : 
     289             :     /* Reset Counters and summary for the current set of minorcycle iterations */
     290        2431 :     itsIterDone = 0;
     291        2431 :     itsIterDiff = -1;
     292             :     //int nSummaryFields = SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     293        2431 :     int nSummaryFields = !itsFullSummary ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
     294        2431 :     itsSummaryMinor.resize( IPosition( 2, nSummaryFields, 0) , true );
     295             : 
     296        4862 :     return returnRecord;
     297             :   }
     298             : 
     299        1179 :   void SIMinorCycleController::setCycleControls(Record &recordIn) {
     300        2358 :     LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
     301             : 
     302        1179 :     if (recordIn.isDefined("cycleniter"))
     303        1179 :       {recordIn.get(RecordFieldId("cycleniter"), itsCycleNiter);}
     304             :     else
     305           0 :       {throw(AipsError("cycleniter not defined in input minor-cycle controller") );}
     306             : 
     307        1179 :     if (recordIn.isDefined("cyclethreshold")) 
     308        1179 :       {recordIn.get(RecordFieldId("cyclethreshold"),itsCycleThreshold);}
     309             :     else
     310           0 :       {throw(AipsError("cyclethreshold not defined in input minor-cycle controller") );}
     311             :      
     312        1179 :     if (recordIn.isDefined("thresholdreached"))
     313        1179 :       {recordIn.get(RecordFieldId("thresholdreached"), itsIsThresholdReached);}
     314             :     else
     315           0 :       { throw(AipsError("thresholdreached not defined in input minor-cycle controller") );}
     316             : 
     317        1179 :     if (recordIn.isDefined("loopgain")) 
     318        1179 :       {recordIn.get(RecordFieldId("loopgain"), itsLoopGain);}
     319             :     else
     320           0 :       {throw(AipsError("loopgain not defined in input minor-cycle controller") );}
     321             : 
     322        1179 :     if (recordIn.isDefined("nsigma"))
     323        1179 :       {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        1179 :     itsMaxCycleIterDone = 0;
     335        1179 :     itsCycleIterDone = 0;
     336        1179 :     itsUpdatedModelFlag = false;
     337        1179 :   }
     338             : 
     339        3141 :   void SIMinorCycleController::resetCycleIter(){
     340        3141 :     itsMaxCycleIterDone = max(itsCycleIterDone, itsMaxCycleIterDone);
     341        3141 :     itsCycleIterDone = 0;
     342        3141 :   }
     343             : 
     344        3141 :   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        9423 :     LogIO os( LogOrigin("SIMinorCycleController", __FUNCTION__ ,WHERE) );
     349        6282 :     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        3141 :     int nSummaryFields = !fullsummary ? 6 : SIMinorCycleController::nSummaryFields;
     355        3141 :     if( shp.nelements() != 2 && shp[0] != nSummaryFields ) 
     356           0 :       throw(AipsError("Internal error in shape of minor-cycle summary record"));
     357        3141 :      itsSummaryMinor.resize( IPosition( 2, nSummaryFields, shp[1]+1 ) , true );
     358             :      // iterations done
     359             :      // make it non-cummulative for all cases
     360        3141 :      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        3141 :      itsSummaryMinor( IPosition(2, 1, shp[1] ) ) = (Double) peakresidual;
     369             :      // model flux
     370        3141 :      itsSummaryMinor( IPosition(2, 2, shp[1] ) ) = (Double) modelflux;
     371             :      // cycle threshold
     372        3141 :      itsSummaryMinor( IPosition(2, 3, shp[1] ) ) = itsCycleThreshold;
     373             :      // mapper id (or multifield id temporary CAS-13683 workaround)
     374        3141 :      itsSummaryMinor( IPosition(2, 4, shp[1] ) ) = deconvolverid;
     375             :      // channel id
     376        3141 :      itsSummaryMinor( IPosition(2, 5, shp[1] ) ) = chan;
     377             :      //if (!uss) {
     378        3141 :      if (fullsummary) {
     379             :          // polarity id
     380          48 :          itsSummaryMinor( IPosition(2, 6, shp[1] ) ) = pol;
     381             :          // cycle start iterations done (ie earliest iterDone for the entire minor cycle)
     382          48 :          itsSummaryMinor( IPosition(2, 7, shp[1] ) ) = cycleStartIter;
     383             :          // starting iterations done
     384          48 :          itsSummaryMinor( IPosition(2, 8, shp[1] ) ) = startIterDone;
     385             :          // starting peak residual
     386          48 :          itsSummaryMinor( IPosition(2, 9, shp[1] ) ) = (Double) startpeakresidual;
     387             :          // starting model flux
     388          48 :          itsSummaryMinor( IPosition(2, 10, shp[1] ) ) = (Double) startmodelflux;
     389             :          // starting peak residual, not limited to the user's mask
     390          48 :          itsSummaryMinor( IPosition(2, 11, shp[1] ) ) = (Double) startpeakresidualnomask;
     391             :          // peak residual, not limited to the user's mask
     392          48 :          itsSummaryMinor( IPosition(2, 12, shp[1] ) ) = (Double) peakresidualnomask;
     393             :          // number of pixels in the mask
     394          48 :          itsSummaryMinor( IPosition(2, 13, shp[1] ) ) = (Double) masksum;
     395             :          // mpi server
     396          48 :          itsSummaryMinor( IPosition(2, 14, shp[1] ) ) = mpiRank;
     397             :          // outlier field id, to be provided in grpcInteractiveCleanManager::mergeMinorCycleSummary
     398          48 :          itsSummaryMinor( IPosition(2, 15, shp[1] ) ) = 0;
     399             :          // stopcode
     400          48 :          itsSummaryMinor( IPosition(2, 16, shp[1] ) ) = stopCode;
     401             :      }
     402             : 
     403        3141 :   }// 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             :