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