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 1063 : SIMinorCycleController::SIMinorCycleController():
38 1063 : itsCycleNiter(0),
39 1063 : itsCycleThreshold(0.0),
40 1063 : itsNsigmaThreshold(0.0),
41 1063 : itsLoopGain(0.1),
42 1063 : itsIsThresholdReached(false),
43 1063 : itsUpdatedModelFlag(false),
44 1063 : itsIterDone(0),
45 1063 : itsCycleIterDone(0),
46 1063 : itsTotalIterDone(0),
47 1063 : itsMaxCycleIterDone(0),
48 1063 : itsPeakResidual(0),
49 1063 : itsIntegratedFlux(0),
50 1063 : itsMaxPsfSidelobe(0),
51 1063 : itsMinResidual(0),itsMinResidualNoMask(0),
52 1063 : itsPeakResidualNoMask(0), itsNsigma(0),
53 1063 : itsMadRMS(0), itsMaskSum(0),
54 : //itsSummaryMinor(IPosition(2,
55 : // SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
56 1063 : itsSummaryMinor(IPosition(2, SIMinorCycleController::nSummaryFields, // temporary CAS-13683 workaround
57 : 0)),
58 1063 : itsDeconvolverID(0)
59 1063 : {}
60 :
61 :
62 :
63 1063 : 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 1199 : 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 2791 : 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 : /* Reset Counters and summary for the current set of minorcycle iterations */
280 2462 : itsIterDone = 0;
281 2462 : itsIterDiff = -1;
282 : //int nSummaryFields = SIMinorCycleController::useSmallSummaryminor() ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
283 2462 : int nSummaryFields = !itsFullSummary ? 6 : SIMinorCycleController::nSummaryFields; // temporary CAS-13683 workaround
284 2462 : itsSummaryMinor.resize( IPosition( 2, nSummaryFields, 0) , true );
285 :
286 : /* Control Variables */
287 2462 : returnRecord.define(RecordFieldId("peakresidual"), itsPeakResidual);
288 2462 : returnRecord.define(RecordFieldId("maxpsfsidelobe"), itsMaxPsfSidelobe);
289 2462 : returnRecord.define( RecordFieldId("peakresidualnomask"), itsPeakResidualNoMask);
290 2462 : returnRecord.define( RecordFieldId("madrms"), itsMadRMS);
291 2462 : returnRecord.define( RecordFieldId("masksum"), itsMaskSum);
292 2462 : returnRecord.define( RecordFieldId("nsigmathreshold"), itsNsigmaThreshold);
293 2462 : returnRecord.define( RecordFieldId("nsigma"), itsNsigma);
294 2462 : returnRecord.define( RecordFieldId("fullsummary"), itsFullSummary);
295 2462 : returnRecord.define(RecordFieldId("summaryminor"), itsSummaryMinor);
296 :
297 4924 : return returnRecord;
298 2462 : }
299 :
300 1192 : void SIMinorCycleController::setCycleControls(Record &recordIn) {
301 2384 : LogIO os( LogOrigin("SIMinorCycleController",__FUNCTION__,WHERE) );
302 :
303 1192 : if (recordIn.isDefined("cycleniter"))
304 1192 : {recordIn.get(RecordFieldId("cycleniter"), itsCycleNiter);}
305 : else
306 0 : {throw(AipsError("cycleniter not defined in input minor-cycle controller") );}
307 :
308 1192 : if (recordIn.isDefined("cyclethreshold"))
309 1192 : {recordIn.get(RecordFieldId("cyclethreshold"),itsCycleThreshold);}
310 : else
311 0 : {throw(AipsError("cyclethreshold not defined in input minor-cycle controller") );}
312 :
313 1192 : if (recordIn.isDefined("thresholdreached"))
314 1192 : {recordIn.get(RecordFieldId("thresholdreached"), itsIsThresholdReached);}
315 : else
316 0 : { throw(AipsError("thresholdreached not defined in input minor-cycle controller") );}
317 :
318 1192 : if (recordIn.isDefined("loopgain"))
319 1192 : {recordIn.get(RecordFieldId("loopgain"), itsLoopGain);}
320 : else
321 0 : {throw(AipsError("loopgain not defined in input minor-cycle controller") );}
322 :
323 1192 : if (recordIn.isDefined("nsigma"))
324 1192 : {recordIn.get(RecordFieldId("nsigma"), itsNsigma);}
325 : else
326 0 : { throw(AipsError(" nsigma is not defined in input minor-cycle controller ") );}
327 :
328 : //if (recordIn.isDefined("fullsummary"))
329 : // {recordIn.get(RecordFieldId("fullsummary"), itsFullSummary);}
330 : //else
331 : // { throw(AipsError(" fullsummary is not defined in input minor-cycle controller ") );}
332 : //os<<" in setCycleControls done... "<<LogIO::POST;
333 :
334 : /* Reset the counters for the new cycle */
335 1192 : itsMaxCycleIterDone = 0;
336 1192 : itsCycleIterDone = 0;
337 1192 : itsUpdatedModelFlag = false;
338 1192 : }
339 :
340 3151 : void SIMinorCycleController::resetCycleIter(){
341 3151 : itsMaxCycleIterDone = max(itsCycleIterDone, itsMaxCycleIterDone);
342 3151 : itsCycleIterDone = 0;
343 3151 : }
344 :
345 3151 : void SIMinorCycleController::addSummaryMinor(uInt deconvolverid, uInt chan, uInt pol,
346 : Int cycleStartIter, Int startIterDone, Float startmodelflux, Float startpeakresidual, Float startpeakresidualnomask,
347 : Float modelflux, Float peakresidual, Float peakresidualnomask, Float masksum, Int mpiRank, Int stopCode, bool fullsummary)
348 : {
349 6302 : LogIO os( LogOrigin("SIMinorCycleController", __FUNCTION__ ,WHERE) );
350 3151 : IPosition shp = itsSummaryMinor.shape();
351 : //bool uss = SIMinorCycleController::useSmallSummaryminor(); // temporary CAS-13683 workaround
352 : //int nSummaryFields = uss ? 6 : SIMinorCycleController::nSummaryFields;
353 : //int nSummaryFields = fullsummary ? 6 : SIMinorCycleController::nSummaryFields;
354 :
355 3151 : int nSummaryFields = !fullsummary ? 6 : SIMinorCycleController::nSummaryFields;
356 3151 : if( shp.nelements() != 2 && shp[0] != nSummaryFields )
357 0 : throw(AipsError("Internal error in shape of minor-cycle summary record"));
358 3151 : itsSummaryMinor.resize( IPosition( 2, nSummaryFields, shp[1]+1 ) , true );
359 : // iterations done
360 : // make it non-cummulative for all cases
361 3151 : itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone - startIterDone;
362 : //if(!fullsummary) {
363 : // itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone - startIterDone;
364 : //}
365 : //else {
366 : // itsSummaryMinor( IPosition(2, 0, shp[1] ) ) = itsIterDone;
367 : // }
368 : // peak residual
369 3151 : itsSummaryMinor( IPosition(2, 1, shp[1] ) ) = (Double) peakresidual;
370 : // model flux
371 3151 : itsSummaryMinor( IPosition(2, 2, shp[1] ) ) = (Double) modelflux;
372 : // cycle threshold
373 3151 : itsSummaryMinor( IPosition(2, 3, shp[1] ) ) = itsCycleThreshold;
374 : // mapper id (or multifield id temporary CAS-13683 workaround)
375 3151 : itsSummaryMinor( IPosition(2, 4, shp[1] ) ) = deconvolverid;
376 : // channel id
377 3151 : itsSummaryMinor( IPosition(2, 5, shp[1] ) ) = chan;
378 : //if (!uss) {
379 3151 : if (fullsummary) {
380 : // polarity id
381 46 : itsSummaryMinor( IPosition(2, 6, shp[1] ) ) = pol;
382 : // cycle start iterations done (ie earliest iterDone for the entire minor cycle)
383 46 : itsSummaryMinor( IPosition(2, 7, shp[1] ) ) = cycleStartIter;
384 : // starting iterations done
385 46 : itsSummaryMinor( IPosition(2, 8, shp[1] ) ) = startIterDone;
386 : // starting peak residual
387 46 : itsSummaryMinor( IPosition(2, 9, shp[1] ) ) = (Double) startpeakresidual;
388 : // starting model flux
389 46 : itsSummaryMinor( IPosition(2, 10, shp[1] ) ) = (Double) startmodelflux;
390 : // starting peak residual, not limited to the user's mask
391 46 : itsSummaryMinor( IPosition(2, 11, shp[1] ) ) = (Double) startpeakresidualnomask;
392 : // peak residual, not limited to the user's mask
393 46 : itsSummaryMinor( IPosition(2, 12, shp[1] ) ) = (Double) peakresidualnomask;
394 : // number of pixels in the mask
395 46 : itsSummaryMinor( IPosition(2, 13, shp[1] ) ) = (Double) masksum;
396 : // mpi server
397 46 : itsSummaryMinor( IPosition(2, 14, shp[1] ) ) = mpiRank;
398 : // outlier field id, to be provided in grpcInteractiveCleanManager::mergeMinorCycleSummary
399 46 : itsSummaryMinor( IPosition(2, 15, shp[1] ) ) = 0;
400 : // stopcode
401 46 : itsSummaryMinor( IPosition(2, 16, shp[1] ) ) = stopCode;
402 : }
403 :
404 3151 : }// end of addSummaryMinor
405 :
406 : // temporary CAS-13683 workaround
407 : /***
408 : Bool SIMinorCycleController::useSmallSummaryminor()
409 : {
410 : if (const char* use_small_summaryminor_p = std::getenv("USE_SMALL_SUMMARYMINOR"))
411 : {
412 : string use_small_summaryminor(use_small_summaryminor_p);
413 : if (use_small_summaryminor.compare("TRUE") == 0 || use_small_summaryminor.compare("true") == 0) {
414 : return true;
415 : }
416 : }
417 : return false;
418 : }
419 : ***/
420 :
421 : } //# NAMESPACE CASA - END
422 :
|
