Line data Source code
1 : //# ClarkCleanProgress.cc: Abstract base class to monitor progress in lattice operations
2 : //# Copyright (C) 1997,1998,1999,2000,2001,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 :
29 : //# Includes
30 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
31 : #include <casacore/images/Images/TempImage.h>
32 : #include <casacore/casa/System/PGPlotter.h>
33 : #include <casacore/casa/BasicMath/Math.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 :
36 : #include <sstream>
37 :
38 :
39 : using namespace casacore;
40 : namespace casa { //# NAMESPACE CASA - BEGIN
41 :
42 0 : ClarkCleanProgress::ClarkCleanProgress(PGPlotter* pgplotter, Int inc)
43 : : itsPgplotter(pgplotter),
44 : currentIndex(0),
45 : currentTotalIterations(0),
46 : currentFluxScale(0.0),
47 : currentMaxResidual(0.0),
48 : currentMinResidual(0.0),
49 : plottingIncrement(inc),
50 0 : baseFlux(0.0)
51 : {
52 0 : iterationNumber.resize(100);
53 0 : maxResiduals.resize(100);
54 0 : posResiduals.resize(100);
55 0 : negResiduals.resize(100);
56 0 : totalFluxes.resize(100);
57 :
58 : // to prevent trailing vector elements from being plotted
59 0 : iterationNumber = 99999999;
60 0 : }
61 :
62 :
63 0 : ClarkCleanProgress::~ClarkCleanProgress()
64 : {
65 0 : }
66 :
67 :
68 : // Call back function
69 0 : Bool ClarkCleanProgress::info(const Bool lastcall,
70 : const Int iteration,
71 : const Int numberIterations,
72 : const Float& maxResid,
73 : const IPosition& posMaxResid,
74 : const Float& totalFlux0,
75 : const Bool majorIteration,
76 : const Bool resetBaseFlux)
77 : {
78 0 : LogIO os(LogOrigin("ClarkCleanProgress", "info()", WHERE));
79 :
80 0 : if(currentIndex == 0) {
81 0 : if(hasPGPlotter()) {
82 0 : currentFluxScale = 1.5*maxResid;
83 0 : currentMinFluxScale = 0.0;
84 0 : currentMaxResidual = 1.5*maxResid;
85 0 : currentMinResidual = currentMaxResidual/100;
86 0 : currentTotalIterations = numberIterations;
87 0 : basicSetUp(numberIterations);
88 : }
89 : }
90 :
91 0 : if (resetBaseFlux && currentIndex != 0) {
92 0 : baseFlux = totalFluxes(currentIndex-1);
93 : }
94 0 : Float totalFlux = totalFlux0 + baseFlux;
95 :
96 : // resize if required
97 0 : if (currentIndex >= totalFluxes.nelements() ) {
98 0 : uInt nn = totalFluxes.nelements();
99 0 : totalFluxes.resize(2*nn+1, true);
100 0 : maxResiduals.resize(2*nn+1, true);
101 0 : negResiduals.resize(2*nn+1, true);
102 0 : posResiduals.resize(2*nn+1, true);
103 0 : iterationNumber.resize(2*nn+1, true);
104 : // Do this so trailing (ie, unfilled) elements don't get plotted
105 0 : for (uInt i=nn; i < 2*nn+1; i++) {
106 0 : iterationNumber(i) = 9999999;
107 : }
108 : }
109 :
110 0 : totalFluxes(currentIndex) = totalFlux;
111 0 : maxResiduals(currentIndex) = maxResid;
112 0 : iterationNumber(currentIndex) = iteration+1;
113 :
114 0 : if (maxResid > 0) {
115 0 : posResiduals(currentIndex) = log10(maxResid);
116 0 : negResiduals(currentIndex) = -20; // ie, out of range
117 0 : } else if (maxResid < 0) {
118 0 : negResiduals(currentIndex) = log10(abs(maxResid));
119 0 : posResiduals(currentIndex) = -20;
120 : }
121 0 : Float myMinFlux = min (0.0, totalFlux);
122 :
123 0 : currentIndex++;
124 :
125 0 : if(hasPGPlotter()) {
126 0 : if ( lastcall || ( (iteration) % plottingIncrement) == 0) {
127 : // Then we get to do a plot
128 :
129 0 : Bool redrawBox = false;
130 :
131 0 : if ( totalFlux > (0.9*currentFluxScale)) {
132 0 : currentFluxScale = abs(3.0 * currentFluxScale);
133 0 : redrawBox = true;
134 : }
135 0 : if (myMinFlux < currentMinFluxScale) {
136 0 : currentMinFluxScale = -abs( 3.0 * myMinFlux);
137 0 : redrawBox = true;
138 : }
139 0 : if (abs(maxResid) < (1.2*currentMinResidual)) {
140 0 : currentMinResidual /= 3.0;
141 0 : redrawBox = true;
142 : }
143 0 : if ( numberIterations > (Int)currentTotalIterations) {
144 0 : currentTotalIterations = numberIterations;
145 0 : redrawBox = true;
146 : }
147 :
148 0 : if (redrawBox || lastcall || resetBaseFlux) {
149 0 : basicSetUp(true);
150 0 : plotVectors();
151 : } else {
152 0 : plotOne(iteration+1, maxResid, totalFlux);
153 : }
154 : }
155 : }
156 :
157 0 : if (majorIteration) {
158 0 : os << "Max Resid = " << maxResid << " at "
159 0 : << posMaxResid << endl;
160 : os << "Iteration " << iteration+1 <<
161 0 : " flux [Jy] = " << totalFlux << endl;
162 : }
163 :
164 0 : return true;
165 : };
166 :
167 : // Call back function
168 0 : Bool ClarkCleanProgress::finalize()
169 :
170 : {
171 0 : LogIO os(LogOrigin("ClarkCleanProgress", "info()", WHERE));
172 :
173 0 : basicSetUp(true);
174 0 : plotVectors();
175 :
176 0 : os << "Max Resid = " << maxResiduals(currentIndex-1)<< endl;
177 0 : os << "Iteration " << iterationNumber(currentIndex-1) <<
178 0 : " flux [Jy] = " << totalFluxes(currentIndex-1) << endl;
179 :
180 0 : return true;
181 : };
182 :
183 : // checks to see if the pointer is non null;
184 : // also checks to see if it is attached to a plot, and
185 : // checks to see if the PGPlotter pointer is dangling
186 0 : Bool ClarkCleanProgress::hasPGPlotter()
187 : {
188 0 : if (itsPgplotter == 0) {
189 0 : return false;
190 : } else {
191 : try {
192 0 : if (itsPgplotter->isAttached()) {
193 0 : return true;
194 : } else {
195 0 : return false;
196 : }
197 0 : } catch (AipsError x) {
198 0 : return false;
199 : }
200 : }
201 : return false;
202 : };
203 :
204 :
205 0 : void ClarkCleanProgress::basicSetUp(Bool doPlot)
206 : {
207 :
208 0 : Float logMinRes = log10(abs(currentMinResidual));
209 0 : Float logMaxRes = log10(abs(currentMaxResidual));
210 0 : Float deltaY = abs(logMaxRes - logMinRes);
211 0 : logMaxRes += 0.05*deltaY;
212 0 : logMinRes -= 0.05*deltaY;
213 0 : Float xMax = Float(currentTotalIterations)*1.15;
214 0 : Float xMin = -0.05*Float(currentTotalIterations);
215 :
216 0 : if(itsPgplotter){
217 :
218 0 : itsPgplotter->sch(0.6);
219 0 : itsPgplotter->sci(1);
220 0 : itsPgplotter->page();
221 0 : itsPgplotter->svp(0.06, 0.94, 0.64, 0.92);
222 0 : itsPgplotter->swin(xMin, xMax, logMinRes, logMaxRes);
223 0 : itsPgplotter->box("BCST", 0, 0, "BCNLST", 0, 0);
224 0 : itsPgplotter->lab(" ", "+ Peak Resid (Jy)", "Components subtracted");
225 0 : itsPgplotter->iden();
226 : {
227 0 : itsPgplotter->sci(1);
228 0 : ostringstream oos;
229 0 : oos << "MaxRes ";
230 0 : itsPgplotter->text(0.85*xMax,
231 0 : (logMaxRes - 0.1*deltaY),
232 0 : oos);
233 : }
234 0 : if (doPlot) {
235 0 : itsPgplotter->pt(iterationNumber, posResiduals, 2);
236 : }
237 :
238 :
239 : // middle graph
240 0 : itsPgplotter->sci(1);
241 0 : itsPgplotter->svp(0.06, 0.94, 0.36, 0.64);
242 0 : itsPgplotter->swin(xMin, xMax, logMaxRes, logMinRes);
243 0 : itsPgplotter->box("BCST", 0, 0, "BCNLST", 0, 0);
244 0 : itsPgplotter->lab(" ", "- Peak Resid (Jy)", " ");
245 : {
246 0 : ostringstream oos;
247 0 : oos << "-MaxNegRes ";
248 0 : itsPgplotter->text(0.85*xMax,
249 0 : (logMaxRes - 0.1*deltaY),
250 0 : oos);
251 : }
252 0 : if (doPlot) {
253 0 : itsPgplotter->pt(iterationNumber, negResiduals, 2);
254 : }
255 :
256 : // lower graph
257 0 : itsPgplotter->sci(3);
258 0 : itsPgplotter->svp(0.06, 0.94, 0.09, 0.36);
259 0 : itsPgplotter->swin(xMin, xMax, currentMinFluxScale, currentFluxScale);
260 0 : itsPgplotter->box("BCNST", 0, 0, "BCNST", 0, 0);
261 0 : itsPgplotter->lab("Number of iterations", "Total Flux", " ");
262 0 : if (doPlot) {
263 0 : itsPgplotter->pt(iterationNumber, totalFluxes, 2);
264 : }
265 :
266 : }
267 :
268 0 : };
269 :
270 0 : void ClarkCleanProgress::plotVectors()
271 : {
272 0 : Float logMinRes = log10(abs(currentMinResidual));
273 0 : Float logMaxRes = log10(abs(currentMaxResidual));
274 0 : Float deltaY = abs(logMaxRes - logMinRes);
275 0 : logMaxRes += 0.05*deltaY;
276 0 : logMinRes -= 0.05*deltaY;
277 0 : Float xMax = Float(currentTotalIterations)*1.15;
278 0 : Float xMin = -0.05*Float(currentTotalIterations);
279 :
280 0 : if(itsPgplotter){
281 0 : itsPgplotter->sch(0.6);
282 :
283 : // top graph
284 0 : itsPgplotter->sci(1);
285 0 : itsPgplotter->svp(0.06, 0.94, 0.64, 0.92);
286 0 : itsPgplotter->swin(xMin, xMax, logMinRes, logMaxRes);
287 0 : itsPgplotter->pt(iterationNumber, posResiduals, 2);
288 :
289 : // middle graph
290 0 : itsPgplotter->sci(1);
291 0 : itsPgplotter->svp(0.06, 0.94, 0.36, 0.64);
292 0 : itsPgplotter->swin(xMin, xMax, logMaxRes, logMinRes);
293 0 : itsPgplotter->pt(iterationNumber, negResiduals, 2);
294 :
295 :
296 : // lower graph
297 0 : itsPgplotter->sci(3);
298 0 : itsPgplotter->svp(0.06, 0.94, 0.09, 0.36);
299 0 : itsPgplotter->swin(xMin, xMax, 0.0, currentFluxScale);
300 0 : itsPgplotter->pt(iterationNumber, totalFluxes, 2);
301 : }
302 0 : };
303 :
304 :
305 :
306 0 : void ClarkCleanProgress::plotOne(const Int iteration,
307 : const Float resid, const Float flux)
308 : {
309 :
310 0 : Float logMinRes = log10(abs(currentMinResidual));
311 0 : Float logMaxRes = log10(abs(currentMaxResidual));
312 0 : Float deltaY = abs(logMaxRes - logMinRes);
313 0 : logMaxRes += 0.05*deltaY;
314 0 : logMinRes -= 0.05*deltaY;
315 0 : Float xMax = Float(currentTotalIterations)*1.15;
316 0 : Float xMin = -0.05*Float(currentTotalIterations);
317 :
318 0 : Vector<Float> x(1);
319 0 : Vector<Float> y(1);
320 0 : x(0) = iteration;
321 0 : if(itsPgplotter){
322 0 : itsPgplotter->sch(0.6);
323 0 : if (resid > 0) {
324 : // top graph
325 0 : itsPgplotter->sci(1);
326 0 : itsPgplotter->svp(0.06, 0.94, 0.64, 0.92);
327 0 : itsPgplotter->swin(xMin, xMax, logMinRes, logMaxRes);
328 0 : y(0) = log10(resid);
329 0 : itsPgplotter->pt(x,y,2);
330 0 : } else if (resid < 0) {
331 : // middle graph
332 0 : itsPgplotter->sci(1);
333 0 : itsPgplotter->svp(0.06, 0.94, 0.36, 0.64);
334 0 : itsPgplotter->swin(xMin, xMax, logMaxRes, logMinRes);
335 0 : y(0) = log10(abs(resid));
336 0 : itsPgplotter->pt(x,y,2);
337 : }
338 :
339 : // lower graph
340 0 : itsPgplotter->sci(3);
341 0 : itsPgplotter->svp(0.06, 0.94, 0.09, 0.36);
342 0 : itsPgplotter->swin(xMin, xMax, 0.0, currentFluxScale);
343 0 : y(0) = flux;
344 0 : itsPgplotter->pt(x,y,2);
345 : }
346 0 : };
347 :
348 :
349 : } //# NAMESPACE CASA - END
350 :
|
