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