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 : 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 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 : }
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 : }
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 : }
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 :
|
