Line data Source code
1 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU Library General Public License as published by
6 : //# the Free Software Foundation; either version 2 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
12 : //# License for more details.
13 : //#
14 : //# You should have received a copy of the GNU Library General Public License
15 : //# along with this library; if not, write to the Free Software Foundation,
16 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: casa-feedback@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 : //# $Id: MultiTermMatrixCleaner.cc 13656 2010-12-04 02:08:02 UrvashiRV$
26 : //
27 : // Implementation of the paper "A multi-scale multi-frequency deconvolution algorithm
28 : // for synthesis imaging in radio interferometry" by Rau and Cornwell
29 : // https://www.aanda.org/articles/aa/pdf/2011/08/aa17104-11.pdf
30 :
31 : // Same include list as in MatrixCleaner.cc
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/Cube.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/MatrixMath.h>
36 : #include <casacore/casa/IO/ArrayIO.h>
37 : #include <casacore/casa/BasicMath/Math.h>
38 : #include <casacore/casa/BasicSL/Complex.h>
39 : #include <casacore/casa/Logging/LogIO.h>
40 : #include <casacore/casa/OS/File.h>
41 : #include <casacore/casa/Containers/Record.h>
42 :
43 : #include <casacore/lattices/LRegions/LCBox.h>
44 : #include <casacore/casa/Arrays/Slicer.h>
45 : #include <casacore/scimath/Mathematics/FFTServer.h>
46 : #include <casacore/casa/OS/HostInfo.h>
47 : #include <casacore/casa/Arrays/ArrayError.h>
48 : #include <casacore/casa/Arrays/ArrayIter.h>
49 : #include <casacore/casa/Arrays/VectorIter.h>
50 :
51 :
52 : #include <casacore/casa/Utilities/GenSort.h>
53 : #include <casacore/casa/BasicSL/String.h>
54 : #include <casacore/casa/Utilities/Assert.h>
55 : #include <casacore/casa/Utilities/Fallible.h>
56 :
57 : #include <casacore/casa/BasicSL/Constants.h>
58 :
59 : #include <casacore/casa/Logging/LogSink.h>
60 : #include <casacore/casa/Logging/LogMessage.h>
61 :
62 : #include <synthesis/MeasurementEquations/MatrixCleaner.h>
63 : #include <casacore/coordinates/Coordinates/TabularCoordinate.h>
64 :
65 : // Additional include files
66 : #include <casacore/lattices/Lattices/SubLattice.h>
67 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
68 : #include <casacore/images/Images/PagedImage.h>
69 :
70 : #include<synthesis/MeasurementEquations/MultiTermMatrixCleaner.h>
71 :
72 : using namespace casacore;
73 : namespace casa { //# NAMESPACE CASA - BEGIN
74 :
75 : #define MIN(a,b) ((a)<=(b) ? (a) : (b))
76 : #define MAX(a,b) ((a)>=(b) ? (a) : (b))
77 :
78 0 : MultiTermMatrixCleaner::MultiTermMatrixCleaner():
79 : MatrixCleaner(),
80 0 : ntaylor_p(0),psfntaylor_p(0),nscales_p(0),nx_p(0),ny_p(0),totalIters_p(0),
81 0 : maxscaleindex_p(0), globalmaxpos_p(IPosition(0)),
82 0 : /*donePSF_p(false),*/donePSP_p(false),doneCONV_p(false),memoryMB_p(0),
83 0 : adbg(false)
84 0 : { }
85 :
86 : /*
87 : MultiTermMatrixCleaner::MultiTermMatrixCleaner(const MultiTermMatrixCleaner & other):
88 : MatrixCleaner(),ntaylor_p(other.ntaylor_p) //And others... minus some...
89 : {
90 : }
91 :
92 : MultiTermMatrixCleaner &MultiTermMatrixCleaner:: operator=(const MultiTermMatrixCleaner & other)
93 : {
94 : operator=(other);
95 : return *this;
96 : }
97 : */
98 :
99 0 : MultiTermMatrixCleaner::~MultiTermMatrixCleaner()
100 : {
101 : //cout << "In MultiTermMatrixCleaner destructor" << endl;
102 0 : }
103 :
104 :
105 0 : Bool MultiTermMatrixCleaner::setscales(const Vector<Float> & scales)
106 : {
107 0 : nscales_p = scales.nelements();
108 0 : scaleSizes_p.resize();
109 0 : scaleSizes_p = scales;
110 0 : totalScaleFlux_p.resize(nscales_p);
111 0 : totalScaleFlux_p.set(0.0);
112 0 : maxScaleVal_p.resize(nscales_p);
113 0 : maxScaleVal_p.set(0.0);
114 0 : maxScalePos_p.resize(nscales_p);
115 0 : globalmaxval_p=-1e+10;
116 0 : globalmaxpos_p=IPosition(2,0);
117 0 : maxscaleindex_p=0;
118 0 : return true;
119 : }
120 :
121 :
122 0 : Bool MultiTermMatrixCleaner::setntaylorterms(const int & nterms)
123 : {
124 0 : ntaylor_p = nterms;
125 0 : psfntaylor_p = 2*nterms-1;
126 0 : totalTaylorFlux_p.resize(ntaylor_p);
127 0 : totalTaylorFlux_p.set(0.0);
128 0 : return true;
129 : }
130 :
131 : // Allocate memory, based on nscales and ntaylor
132 0 : Bool MultiTermMatrixCleaner::initialise(Int nx, Int ny)
133 : {
134 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "initialise()", WHERE));
135 :
136 : /* Verify Image Shapes */
137 : // nx_p = model.shape(0);
138 0 : nx_p = nx;
139 0 : ny_p = ny;
140 :
141 0 : if(adbg) os << "Checking shapes" << LogIO::POST;
142 :
143 : /* Verify nscales_p and ntaylor_p */
144 0 : AlwaysAssert(nscales_p>0, AipsError);
145 0 : AlwaysAssert(ntaylor_p>0, AipsError);
146 0 : AlwaysAssert(scaleSizes_p.nelements()==static_cast<uInt>(nscales_p), AipsError);
147 :
148 0 : if(adbg) os << "Verify Scale sizes" << LogIO::POST;
149 0 : verifyScaleSizes();
150 :
151 : /* Calculate PSF/scale support - after verifying scale sizes and before allocating memory */
152 :
153 : // PSF support : main lobe width x how many times this width the patch should cover
154 0 : Float maxscalesize = scaleSizes_p[nscales_p-1];
155 :
156 : //// N times the size of the main lobe of the PSF at the max scale size
157 0 : Float psfbeam = 4.0;
158 0 : Float nbeams = 20.0;
159 0 : Int psupport = findBeamPatch(maxscalesize, nx_p, ny_p, psfbeam, nbeams);
160 : // //----------------------Encapsulated in the findBeamPatch() method from here...-----------
161 : // Int psupport = (Int) ( sqrt( psfbeam*psfbeam + maxscalesize*maxscalesize ) * nbeams );
162 :
163 : // // At least this big...
164 : // if(psupport < psfbeam*nbeams ) psupport = static_cast<Int>(psfbeam*nbeams);
165 :
166 : // // Not too big...
167 : // if(psupport > nx_p || psupport > ny_p) psupport = MIN(nx_p,ny_p);
168 :
169 : // // Make it even.
170 : // if (psupport%2 != 0) psupport -= 1;
171 : // //----------------------...till here------------------------------------
172 :
173 :
174 : // Full image
175 : //psfsupport_p = IPosition(2,MIN(nx_p,ny_p),MIN(nx_p,ny_p));
176 : //psfpeak_p = IPosition(2,nx_p/2, ny_p/2);
177 :
178 : // Inner quater image
179 : //psfsupport_p = IPosition(2,MIN(nx_p/2,ny_p/2),MIN(nx_p/2,ny_p/2));
180 : //psfpeak_p = IPosition(2,nx_p/4, ny_p/4);
181 :
182 : // Generic support-size.
183 0 : psfsupport_p = IPosition(2,psupport,psupport);
184 0 : psfpeak_p = IPosition(2,psupport/2, psupport/2);
185 :
186 0 : os << "Using a PSF patch of " << psupport << " pixels on each side for minor-cycle updates." << endl;
187 :
188 : /* Force the scale images to be */
189 :
190 :
191 0 : if(adbg) os << "Start allocating mem" << LogIO::POST;
192 :
193 : /* Allocate memory for many many TempMatrices. */
194 : /* Check on a return value of -1 - for insuffucient memory */
195 0 : if( allocateMemory() == -1 ) return false;
196 :
197 : /* FFTServer */
198 0 : fftcomplex = FFTServer<Float,Complex>(IPosition(2,nx_p,ny_p));
199 :
200 : /* Create the scale functions and their FTs */
201 0 : setupScaleFunctions();
202 :
203 0 : totalIters_p=0;
204 0 : prev_max_p = 1e+10;
205 0 : min_max_p = 1e+10;
206 0 : rmaxval_p = -1.0;
207 :
208 0 : if(adbg) os << "Finished initializing MultiTermMatrixCleaner" << LogIO::POST;
209 0 : return true;
210 0 : }
211 :
212 :
213 0 : Bool MultiTermMatrixCleaner::buildImagePatches()
214 : {
215 :
216 0 : gip = IPosition(2,nx_p,ny_p);
217 :
218 : /* The update region. */
219 0 : IPosition inc(2,1,1);
220 0 : blc_p = IPosition(globalmaxpos_p - psfsupport_p/2);
221 0 : trc_p = IPosition(globalmaxpos_p + psfsupport_p/2 - IPosition(2,1,1));
222 : //cout << "residual box 1 : " << blc << trc << endl;
223 0 : LCBox::verify(blc_p, trc_p, inc, gip);
224 : //cout << "residual box 2 : " << blc << trc << endl;
225 :
226 :
227 : /* Shifted region, with the psf at the globalmaxpos_p. */
228 0 : blcPsf_p = IPosition(blc_p + psfpeak_p - globalmaxpos_p); // OLD
229 0 : trcPsf_p = IPosition(trc_p + psfpeak_p - globalmaxpos_p); // OLD
230 : //cout << "psf box 1 : " << blcPsf << trcPsf << endl;
231 0 : LCBox::verify(blcPsf_p, trcPsf_p, inc, psfsupport_p); // NEW
232 : //cout << "psf box 2 : " << blcPsf << trcPsf << endl;
233 :
234 : /* Reconcile box sizes/locations with the image size */
235 0 : makeBoxesSameSize(blc_p,trc_p,blcPsf_p,trcPsf_p);
236 :
237 : //cout << "maxpos : " << globalmaxpos_p << " blc : " << blc_p << " trc : " << trc_p << " blcPsf : " << blcPsf_p << " trcPsf : " << trcPsf_p << endl;
238 :
239 0 : return true;
240 0 : }
241 :
242 :
243 0 : Bool MultiTermMatrixCleaner::setpsf(int order, Matrix<Float> & psf)
244 : {
245 0 : AlwaysAssert((order>=(int)0 && order<(int)vecPsfFT_p.nelements()), AipsError);
246 0 : if(order==0) AlwaysAssert(validatePsf(psf), AipsError);
247 : // Directly store the FFT of the PSFs.
248 : // The FT'd matrix is reshaped automatically.
249 0 : fftcomplex.fft0(vecPsfFT_p[order],psf,false);
250 : //fftcomplex.flip(vecPsfFT_p[order],false,false);
251 : // Matrix<Float> temp(real(vecPsfFT_p[order]));
252 : // writeMatrixToDisk("psfFT_"+String::toString(order), temp);
253 0 : return true;
254 : }
255 :
256 : /* Input : Dirty Images */
257 0 : Bool MultiTermMatrixCleaner::setresidual(int order, Matrix<Float> & dirty)
258 : {
259 0 : AlwaysAssert((order>=(int)0 && order<(int)vecDirty_p.nelements()), AipsError);
260 0 : vecDirty_p[order].reference(dirty);
261 0 : return true;
262 : }
263 :
264 : /* Input : Model Component Image */
265 : //TODO : This is an extra copy of the model image. Why not hold it by reference ???
266 0 : Bool MultiTermMatrixCleaner::setmodel(int order, Matrix<Float> & model)
267 : {
268 0 : AlwaysAssert((order>=(int)0 && order<(int)vecModel_p.nelements()), AipsError);
269 0 : vecModel_p[order].assign(model);
270 0 : vecInitialModel_p[order].assign(model);
271 0 : totalTaylorFlux_p[order] = (sum(vecModel_p[order]));
272 0 : return true;
273 : }
274 :
275 : /* Input : Mask */
276 :
277 0 : Bool MultiTermMatrixCleaner::setmask(Matrix<Float> & mask)
278 : {
279 0 : if(itsMask.null())
280 0 : itsMask = new Matrix<Float>(mask); // it's a counted ptr
281 : else
282 : {
283 0 : AlwaysAssert(itsMask->shape()==mask.shape(), AipsError);
284 0 : (*itsMask).assign(mask);
285 : }
286 0 : return true;
287 : }
288 :
289 : /* Output : Model Component Image */
290 :
291 0 : Bool MultiTermMatrixCleaner::getmodel(int order, Matrix<Float> & model)
292 : {
293 0 : AlwaysAssert((order>=(int)0 && order<(int)vecModel_p.nelements()), AipsError);
294 0 : model.assign(vecModel_p[order]);
295 0 : return true;
296 : }
297 :
298 : /* Output Residual Image */
299 :
300 :
301 0 : Bool MultiTermMatrixCleaner::getresidual(int order, Matrix<Float> & residual)
302 : {
303 0 : AlwaysAssert((order>=(int)0 && order<(int)vecDirty_p.nelements()), AipsError);
304 : //AlwaysAssert(residual, AipsError);
305 0 : residual.assign(vecDirty_p[order]);
306 : ///residual.assign( matR_p[IND2(order,0)] ); // This is the one that gets updated during iters.
307 :
308 0 : return true;
309 : }
310 :
311 :
312 :
313 : /* Output Hessian matrix */
314 0 : Bool MultiTermMatrixCleaner::getinvhessian(Matrix<Double> & invhessian)
315 : {
316 0 : invhessian.resize((invMatA_p[0]).shape());
317 0 : invhessian = (invMatA_p[0]);
318 0 : return true;
319 : }
320 :
321 :
322 : /* Do the deconvolution */
323 0 : Int MultiTermMatrixCleaner::mtclean(Int maxniter, Float stopfraction, Float inputgain, Float userthreshold)
324 : {
325 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
326 0 : if(adbg)os << "SOLVER for Multi-Frequency Synthesis deconvolution" << LogIO::POST;
327 :
328 : /* Initialize some variables */
329 :
330 0 : maxniter_p = maxniter;
331 0 : stopfraction_p= stopfraction;
332 0 : inputgain_p=inputgain;
333 0 : userthreshold_p=userthreshold;
334 :
335 : //cout << "MTMC : User threshold : " << userthreshold_p << endl;
336 :
337 :
338 0 : Int convergedflag = 0;
339 0 : Float fluxlimit = -1;
340 0 : Float loopgain = 0.5;
341 0 : if(inputgain>(Float)0.0) loopgain=inputgain;
342 0 : Int criterion=5; // This chooses among a list of 'penalty functions'
343 0 : itercount_p=0;
344 0 : Int iterdone = 0;
345 :
346 : /* Compute all convolutions and the matrix A */
347 : /* If matrix is not invertible, return ! */
348 : /* This is done only once - for the first major cycle. Null operation otherwise */
349 0 : if( doneCONV_p == false )
350 : {
351 0 : if( computeHessianPeak() == -2 )
352 0 : return -2;
353 : }
354 :
355 : /* Set up the Mask image */
356 : /* This must be after the Hessian is computed.. */
357 0 : setupUserMask();
358 :
359 : /* Compute the convolutions of the current residual with all PSFs and scales */
360 0 : os << "Calculating convolutions of residual images with scales " << LogIO::POST;
361 0 : computeRHS();
362 :
363 : /* Check if peak RHS is already less than the threshold */
364 : /* This step computes the flux limit also - when called here... */
365 0 : convergedflag = checkConvergence(criterion,fluxlimit,loopgain);
366 0 : if(convergedflag==1) return 0;
367 :
368 : /* Compute the flux limits that determine the depth of the minor cycles. */
369 : // computeFluxLimit(fluxlimit,thresh);
370 :
371 : /********************** START MINOR CYCLE ITERATIONS ***********************/
372 : //os << "Doing deconvolution iterations..." << LogIO::POST;
373 :
374 0 : for(itercount_p=0;itercount_p<maxniter_p;itercount_p++)
375 : {
376 0 : globalmaxval_p=-1e+10;
377 :
378 : // if itercount_p==0, set the blc/trc to the full image size.
379 : // For later iterations, it will get updated according to globalmaxpos_p and psfsupport_p
380 : // when buildImagePatches() is called.
381 0 : if(itercount_p==0)
382 : {
383 0 : blc_p = IPosition(2,0,0);
384 0 : trc_p = IPosition(2,nx_p-1,ny_p-1);
385 : }
386 :
387 0 : Int scale=0;
388 0 : Int ntaylor=ntaylor_p;
389 0 : IPosition blc(blc_p), trc(trc_p);
390 : //OMP// #pragma omp parallel default(shared) private(scale) firstprivate(ntaylor,criterion,blc,trc)
391 : {
392 : //OMP// #pragma omp for
393 0 : for(scale=0;scale<nscales_p;scale++)
394 : {
395 : /* Solve the matrix eqn for all pixels */
396 0 : solveMatrixEqn(ntaylor,scale,blc,trc);
397 : /* Choose a component from the list of solutions. Record the location for each scale.*/
398 : // Calculate penalty function
399 : // Find max across all pixels.
400 0 : chooseComponent(ntaylor, scale,criterion,blc,trc);
401 : }
402 : }//end pragma omp
403 :
404 : /* Find the best component over all scales */
405 0 : for(Int scale=0;scale<nscales_p;scale++)
406 : {
407 0 : if((maxScaleVal_p[scale]*scaleBias_p[scale]) > globalmaxval_p)
408 : {
409 0 : globalmaxval_p = maxScaleVal_p[scale]*scaleBias_p[scale];
410 0 : globalmaxpos_p = maxScalePos_p[scale];
411 0 : maxscaleindex_p = scale;
412 : }
413 : }
414 :
415 : /* Update the image and psf patch sizes according to the
416 : current globalmaxval and psfsupport.
417 : This patch is over which the next-iteration's solveMatrixEqn is computed */
418 0 : buildImagePatches();
419 :
420 :
421 : /* Update the current solution by this chosen step */
422 0 : updateModelAndRHS(loopgain);
423 :
424 : /* Increment iteration count */
425 0 : totalIters_p++;
426 0 : iterdone++;
427 :
428 : /* Check for convergence */
429 0 : convergedflag = checkConvergence(criterion,fluxlimit,loopgain);
430 :
431 : /* Reached stopping threshold for this major cycle */
432 : /* Break out of minor-cycle loop, but signal that major cycles should continue */
433 0 : if(convergedflag)
434 : {
435 0 : break;
436 : }
437 :
438 0 : }//end minor cycle
439 0 : if(convergedflag==0)
440 : {
441 0 : os << "Reached max number of iterations for this minor cycle " << LogIO::POST;
442 : }
443 :
444 : /* returning because of non-convergence */
445 0 : if(convergedflag == -1)
446 : {
447 0 : os << "Stopping minor cycle iterations because of possible divergence." << LogIO::POST;
448 : //return(-1);
449 : }
450 :
451 : /********************** END MINOR CYCLE ITERATIONS ***********************/
452 :
453 : /* Print out flux counts so far */
454 : //if(adbg)
455 : {
456 0 : os << "Total flux by scale :";
457 0 : for(Int scale=0;scale<nscales_p;scale++) os << " [" << scaleSizes_p[scale] << "]: " << totalScaleFlux_p[scale] ;
458 0 : os << " (in this run) " << LogIO::POST;
459 0 : os << "Total flux by Taylor coefficient :";
460 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++) os << " [" << taylor << "]: " << totalTaylorFlux_p[taylor] ;
461 0 : os << LogIO::POST;
462 : // for(Int taylor=0;taylor<ntaylor_p;taylor++) os << "Taylor " << taylor << " has total flux = " << totalTaylorFlux_p[taylor] << LogIO::POST;
463 : }
464 :
465 : /* Write out the multiscale model images - to disk */
466 : // for(Int scale=0;scale<nscales_p;scale++)
467 : // writeMatrixToDisk("modelimage_scale_"+String::toString(scale) , vecScaleModel_p[scale]);
468 :
469 :
470 : // At this stage, vecDirty_p contains the original residual images from the start of minorcycle iterations,
471 : // Here, we need to update vecDirty_p with image-domain residuals that account for the model
472 : // components that have been found in this call to mtclean().
473 : // Math : Update residual : Ires_t = Ires_t - sum_q [ Ipsf_tq * Im_q ]
474 : // where Ipsf_tq = sum_nu ( w^{t+q} Ipsf_nu = Ipsf_{00, 11, 12/21}
475 0 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
476 : {
477 0 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
478 : {
479 0 : Matrix<Float>newModel;
480 0 : newModel = vecModel_p[taylor2] - vecInitialModel_p[taylor2];
481 : // Convolve model with psf
482 0 : Matrix<Complex> modelFT;
483 0 : fftcomplex.fft0(modelFT, newModel, false);
484 0 : modelFT= ( vecPsfFT_p[taylor1+taylor2] ) * modelFT;
485 0 : Matrix<Float> smoothMod(vecDirty_p[taylor1].shape());
486 0 : fftcomplex.fft0(smoothMod, modelFT, false);
487 0 : fftcomplex.flip(smoothMod, false, false);
488 : // Subtract from initial residual (vecDirty_p)
489 0 : vecDirty_p[taylor1] = vecDirty_p[taylor1] - smoothMod;
490 0 : }
491 : }
492 : // Note : In the restore step, vecDirty_p is over-written and reused in computeprincipalsolution.
493 :
494 : /* Return the number of minor-cycle iterations completed in this run */
495 0 : return iterdone;
496 0 : }
497 :
498 : /* Indexing Rules... */
499 0 : Int MultiTermMatrixCleaner::IND2(Int taylor, Int scale)
500 : {
501 0 : return taylor * nscales_p + scale;
502 : }
503 :
504 0 : Int MultiTermMatrixCleaner::IND4(Int taylor1, Int taylor2, Int scale1, Int scale2)
505 : {
506 0 : Int tt1=taylor1;
507 0 : Int tt2=taylor2;
508 0 : Int ts1=scale1;
509 0 : Int ts2=scale2;
510 0 : scale1 = MAX(ts1,ts2);
511 0 : scale2 = MIN(ts1,ts2);
512 0 : taylor1 = MAX(tt1,tt2);
513 0 : taylor2 = MIN(tt1,tt2);
514 0 : Int totscale = nscales_p*(nscales_p+1)/2;
515 0 : return ((taylor1*(taylor1+1)/2)+taylor2)*totscale + ((scale1*(scale1+1)/2)+scale2);
516 : }
517 :
518 :
519 : /* Check if scale sizes are appropriate to the image size
520 : If some scales are too big, ignore them.
521 : Reset nscales_p BEFORE all arrays are allocated */
522 0 : Int MultiTermMatrixCleaner::verifyScaleSizes()
523 : {
524 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "verifyScaleSizes()", WHERE));
525 0 : Vector<Int> scaleflags(nscales_p); scaleflags=0;
526 0 : for(Int scale=0;scale<nscales_p;scale++)
527 : {
528 0 : if(scaleSizes_p[scale] > nx_p/2 || scaleSizes_p[scale] > ny_p/2)
529 : {
530 0 : scaleflags[scale]=1;
531 0 : os << LogIO::WARN << "Scale size of " << scaleSizes_p[scale] << " pixels is too big for an image size of " << nx_p << " x " << ny_p << " pixels. This scale size will be ignored. " << LogIO::POST;
532 : }
533 : }
534 0 : if(sum(scaleflags)>0) // prune the scalevector and change nscales_p
535 : {
536 0 : Vector<Float> newscalesizes(nscales_p - sum(scaleflags));
537 0 : uInt i=0;
538 0 : for(Int scale=0;scale<nscales_p;scale++)
539 : {
540 0 : if(!scaleflags[scale])
541 : {
542 0 : AlwaysAssert(i<newscalesizes.nelements(),AipsError);
543 0 : newscalesizes[i]=scaleSizes_p[scale];
544 0 : i++;
545 : }
546 : }
547 0 : AlwaysAssert(i==newscalesizes.nelements(), AipsError);
548 0 : scaleSizes_p.assign(newscalesizes);
549 0 : nscales_p = newscalesizes.nelements();
550 0 : totalScaleFlux_p.resize(nscales_p,true);
551 0 : maxScaleVal_p.resize(nscales_p,true);
552 0 : maxScalePos_p.resize(nscales_p,true);
553 :
554 0 : }
555 0 : os << "Scale sizes to be used for deconvolution : " << scaleSizes_p << LogIO::POST;
556 :
557 0 : return scaleSizes_p.shape()(0);
558 0 : }
559 :
560 : /*************************************
561 : * Allocate Memory
562 : *************************************/
563 :
564 0 : Int MultiTermMatrixCleaner::allocateMemory()
565 : {
566 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "allocateMemory()", WHERE));
567 : // Define max memory usage. (half of available);
568 :
569 0 : Int ntotal4d = (nscales_p*(nscales_p+1)/2) * (ntaylor_p*(ntaylor_p+1)/2);
570 :
571 0 : Int nHess = ntotal4d; // cubeA
572 0 : Int ntempfull = 1 // scratch
573 0 : + nscales_p * 3 // vecScales, vecScaleMasks, vecWork_p, ////vecScaleModel_p
574 0 : + ntaylor_p // vecModel (vecDirty is a ref)
575 0 : + 2 * nscales_p * ntaylor_p; // matR and matCoeff
576 :
577 0 : Int ntemphalf = 2 // scratch
578 0 : + nscales_p // vecSCalesFT
579 0 : + psfntaylor_p; // vecPsfFT
580 :
581 0 : Int numMB = static_cast<Int>(( nx_p*ny_p*4*(ntempfull + ntemphalf/2.0) + psfsupport_p[0]*psfsupport_p[1]*nHess )/(1024*1024));
582 0 : memoryMB_p = Double(HostInfo::memoryTotal()/1024);
583 :
584 0 : if(adbg) os << "This algorithm needs to allocate " << numMB << " MBytes." << LogIO::POST;
585 0 : if (numMB > 0.75*memoryMB_p)
586 : {
587 0 : os << LogIO::WARN << "This algorithm needs to allocate " << numMB << " MBytes for " << ntempfull + ntemphalf << " Matrices. " << LogIO::POST;
588 0 : os << LogIO::WARN << "Available memory for this process is " << memoryMB_p << " MB. Please reduce imsize/nscales/nterms and try again " << LogIO::POST;
589 0 : return -1;
590 : }
591 :
592 : // shape of all matrices
593 0 : gip = IPosition(2,nx_p,ny_p);
594 0 : IPosition tgip(2,ntaylor_p,ntaylor_p);
595 :
596 : // Small HessianPeak matrix to be inverted for each point..
597 0 : matA_p.resize(nscales_p); invMatA_p.resize(nscales_p);
598 0 : for(Int i=0;i<nscales_p;i++)
599 : {
600 0 : matA_p[i].resize(tgip);
601 0 : invMatA_p[i].resize(tgip);
602 : }
603 :
604 : // I_D
605 0 : dirtyFT_p.resize();
606 :
607 : // Temporary work-holder
608 0 : cWork_p.resize();
609 : //tWork_p.resize(gip);
610 :
611 : // Scales
612 0 : vecScales_p.resize(nscales_p);
613 0 : vecScalesFT_p.resize(nscales_p);
614 0 : vecScaleMasks_p.resize(nscales_p);
615 0 : vecWork_p.resize(nscales_p);
616 : // vecScaleModel_p.resize(nscales_p);
617 0 : for(Int i=0;i<nscales_p;i++)
618 : {
619 : //vecScales_p[i].resize(gip); // OLD
620 0 : vecScales_p[i].resize(psfsupport_p)
621 : ;
622 0 : vecScalesFT_p[i].resize();
623 0 : vecScaleMasks_p[i].resize(gip);
624 0 : vecWork_p[i].resize(gip);
625 : // vecScaleModel_p[i].resize(gip);
626 : // vecScaleModel_p[i] = 0.0;
627 : }
628 :
629 : // Psfs and Models
630 0 : vecPsfFT_p.resize(psfntaylor_p);
631 0 : for(Int i=0;i<psfntaylor_p;i++)
632 : {
633 0 : vecPsfFT_p[i].resize();
634 : }
635 :
636 : // Dirty/Residual Images
637 0 : vecDirty_p.resize(ntaylor_p);
638 0 : vecModel_p.resize(ntaylor_p);
639 0 : vecInitialModel_p.resize(ntaylor_p);
640 0 : for(Int i=0;i<ntaylor_p;i++)
641 : {
642 0 : vecDirty_p[i].resize();
643 0 : vecModel_p[i].resize(gip);
644 0 : vecInitialModel_p[i].resize(gip);
645 : }
646 :
647 : // (Psf * Scales) * (Psf * Scales)
648 0 : cubeA_p.resize(ntotal4d);
649 :
650 0 : for(Int i=0;i<ntotal4d;i++)
651 : {
652 : //cubeA_p[i].resize(gip); // OLD
653 0 : cubeA_p[i].resize(psfsupport_p); // NEW
654 : }
655 :
656 : // I_D * (Psf * Scales)
657 0 : matR_p.resize(ntaylor_p*nscales_p);
658 : // Coefficients to be solved for.
659 0 : matCoeffs_p.resize(ntaylor_p*nscales_p);
660 :
661 0 : for(Int i=0;i<ntaylor_p*nscales_p;i++)
662 : {
663 0 : matR_p[i].resize(gip);
664 0 : matCoeffs_p[i].resize(gip);
665 0 : matCoeffs_p[i] = 0.0;
666 : }
667 :
668 0 : return 0;
669 0 : }
670 :
671 : /* Make masks for each scale size */
672 : /* At the end, force a scale-dependent border */
673 0 : Int MultiTermMatrixCleaner::setupUserMask()
674 : {
675 0 : if(itsMask.null())
676 : {
677 : /* Make a mask the full size, for all scales */
678 0 : for(Int scale=0;scale<nscales_p;scale++)
679 : {
680 0 : vecScaleMasks_p[scale] = 1.0;
681 : }
682 : }
683 : else
684 : {
685 : /* Compute the convolutions of the input mask with each scale size */
686 0 : if(adbg) os << "Convolving user mask with scales, and using 0.1 as a threshold." << LogIO::POST;
687 0 : Matrix<Complex> maskft;
688 0 : fftcomplex.fft0(maskft , *itsMask , false);
689 :
690 0 : for(Int scale=0;scale<nscales_p;scale++)
691 : {
692 0 : AlwaysAssert(maskft.shape() == vecScalesFT_p[scale].shape(),AipsError);
693 0 : cWork_p.assign(maskft * vecScalesFT_p[scale]);
694 0 : fftcomplex.fft0( vecScaleMasks_p[scale] , cWork_p , false );
695 0 : fftcomplex.flip( vecScaleMasks_p[scale] , false , false);
696 :
697 0 : for (Int j=0 ; j < (itsMask->shape())(1); ++j)
698 0 : for (Int k =0 ; k < (itsMask->shape())(0); ++k)
699 : {
700 0 : if(itsMaskThreshold > (Float)0.0) // if negative, use the continuous mask
701 : {
702 0 : (vecScaleMasks_p[scale])(k,j) = (vecScaleMasks_p[scale])(k,j) > (Float)0.1 ? 1.0 : 0.0;
703 : }
704 : }
705 :
706 : // writeMatrixToDisk("scalemask."+String::toString(scale), vecScaleMasks_p[scale]);
707 : }// end of for scale
708 :
709 : /* TO DO... maybe :
710 : Map some local variables to those in the parant MatrixCleaner class
711 : so that its functions can be used. Later, clean this up by directly using
712 : the MatrixCleaner variables. Currently not possible because the
713 : function that creates scales and populates itsScaleXfrs (setscales),
714 : does a lot of extra stuff that we don't need here. */
715 :
716 : /* Must reuse itsNscales, itsScaleSizes, itsScales, itsScaleXfrs, itsScaleMasks, itsScalesValid
717 : and eliminate : vecScales_p, scaleSizes_p, vecScalesFT_p, vecScaleMasks_p */
718 : /*
719 : itsScaleXfrs.resize(nscales_p);
720 : itsNscales = nscales_p;
721 : for(Int scale=0;scale<nscales_p;scale++)
722 : {
723 : itsScaleXfrs[scale].reference(vecScalesFT_p[scale]);
724 : }
725 : itsScalesValid=true;
726 :
727 : makeScaleMasks();
728 : AlwaysAssert(nscales_p==itsScaleMasks.nelements(),AipsError);
729 :
730 : for(Int scale=0;scale<nscales_p;scale++)
731 : {
732 : cout << "Shapes for scale : " << scale << itsScaleMasks[scale] << endl;
733 : vecScaleMasks_p[scale].reference(itsScaleMasks[scale]);
734 : }
735 : */
736 :
737 0 : }// end of if else
738 :
739 :
740 : /* Set the edges of the masks according to the scale size */
741 : // Set the values OUTSIDE the box to zero....
742 : // "The vecScaleMasks_p are used when finding the peakres for each iteration. If we subtract
743 : // the (tt/scale convolved) psf from too close to the edge of the image, we could add a hard edge
744 : // to the image, which would result in ringing in an fft. Therefore, we use a border to ensure
745 : // that we don't pick a peakres position too close to the edge." - approximately what Urvashi said
746 : // Start at 1: don't add a 1-pixel border to vecScaleMasks_p[0], because that mask is used to find
747 : // the peakres in checkConvergence(), which could result in the minor cycle disagreeing with the
748 : // major cycle about what the peakres value is if the peakres is at the edge.
749 0 : for(Int scale=1;scale<nscales_p;scale++)
750 : {
751 0 : Int border = (Int)(scaleSizes_p[scale]*1.5);
752 : // bottom
753 0 : IPosition blc1(2, 0 , 0 );
754 0 : IPosition trc1(2,nx_p-1, border );
755 0 : IPosition inc1(2, 1);
756 0 : LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
757 0 : (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
758 : // top
759 0 : blc1[0]=0; blc1[1]=ny_p-border-1;
760 0 : trc1[0]=nx_p-1; trc1[1]=ny_p-1;
761 0 : LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
762 0 : (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
763 : // left
764 0 : blc1[0]=0; blc1[1]=border;
765 0 : trc1[0]=border; trc1[1]=ny_p-border-1;
766 0 : LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
767 0 : (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
768 : // right
769 0 : blc1[0]=nx_p-border-1; blc1[1]=border;
770 0 : trc1[0]=nx_p; trc1[1]=ny_p-border-1;
771 0 : LCBox::verify(blc1,trc1,inc1,vecScaleMasks_p[scale].shape());
772 0 : (vecScaleMasks_p[scale])(blc1,trc1) = 0.0;
773 :
774 : // writeMatrixToDisk("maskforscale_"+String::toString(scale),vecScaleMasks_p[scale]);
775 0 : }
776 :
777 0 : return 0;
778 : }/* end of setupUserMask() */
779 :
780 :
781 : /***************************************
782 : * Set up the Blobs of various scales.
783 : ****************************************/
784 :
785 0 : Int MultiTermMatrixCleaner::setupScaleFunctions()
786 : {
787 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "setupScaleFunctions", WHERE));
788 : // Set the scale sizes
789 0 : if(scaleSizes_p.nelements()==0)
790 : {
791 0 : scaleSizes_p.resize(nscales_p);
792 0 : Float scaleInc = 2.0;
793 0 : scaleSizes_p[0] = 0.0;
794 0 : for (Int scale=1; scale<nscales_p;scale++)
795 : {
796 0 : scaleSizes_p[scale] = scaleInc * pow(10.0, (Float(scale)-2.0)/2.0) ;
797 : }
798 : }
799 :
800 : /* SCALE BIAS : Can get rid of this entirely. Leaving it in, just in case scalebiases are needed later */
801 0 : scaleBias_p.resize(nscales_p);
802 0 : totalScaleFlux_p.resize(nscales_p);
803 : //Float prefScale=2.0;
804 : //Float fac=6.0;
805 0 : if(nscales_p>1)
806 : {
807 0 : for(Int scale=0;scale<nscales_p;scale++)
808 : {
809 0 : scaleBias_p[scale] = 1 - itsSmallScaleBias * scaleSizes_p[scale]/scaleSizes_p(nscales_p-1);
810 : //scaleBias_p[scale] = 1 - 0.4 * scaleSizes_p[scale]/scaleSizes_p(nscales_p-1);
811 : //scaleBias_p[scale] = 1.0;
812 : //////scaleBias_p[scale] = pow((Float)scale/fac,prefScale)*exp(-1.0*scale/fac)/(pow(prefScale/fac,prefScale)*exp(-1.0*prefScale/fac));
813 : //scaleBias_p[scale] = pow((Float)(scale+1)/fac,prefScale)*exp(-1.0*(scale+1)/fac);
814 0 : os << "scale " << scale+1 << " = " << scaleSizes_p(scale) << " pixels with bias = " << scaleBias_p[scale] << LogIO::POST;
815 0 : totalScaleFlux_p[scale]=0.0;
816 : }
817 : }
818 0 : else scaleBias_p[0]=1.0;
819 :
820 : // Compute the scaled blobs - prolate spheroids with tapering and truncation
821 : // vecScales_p, scaleSizes_p, vecScalesFT_p
822 0 : if(!donePSP_p)
823 : {
824 : // Compute h(s1), h(s2),... depending on the number of scales chosen.
825 : // NSCALES = 1;
826 0 : if(adbg) os << "Calculating scales and their FTs " << LogIO::POST;
827 :
828 0 : for (Int scale=0; scale<nscales_p;scale++)
829 : {
830 : /*
831 : // First make the scale
832 : makeScale(vecScales_p[scale], scaleSizes_p(scale));
833 : // Now store the XFR (shape of FT is set automatically)
834 : fftcomplex.fft0(vecScalesFT_p[scale] , vecScales_p[scale] , false);
835 : //fftcomplex.flip(vecScalesFT_p[scale] , false , false);
836 : */
837 : // First make the scale
838 0 : makeScale(vecWork_p[0], scaleSizes_p(scale));
839 : // Now store the XFR (shape of FT is set automatically)
840 0 : fftcomplex.fft0(vecScalesFT_p[scale] , vecWork_p[0] , false);
841 : // Copy the scale onto the smaller vecScales image.
842 0 : IPosition immid(2,nx_p/2, ny_p/2);
843 0 : Matrix<Float> psfpatch = ( vecWork_p[0] ) (immid-psfsupport_p/2,immid+psfsupport_p/2-IPosition(2,1,1));
844 0 : vecScales_p[scale] = psfpatch;
845 :
846 0 : }
847 0 : donePSP_p=true;
848 : }
849 :
850 0 : return 0;
851 0 : }/* end of setupBlobs() */
852 :
853 :
854 :
855 : /***************************************
856 : * Compute convolutions and the A matrix.
857 : ****************************************/
858 :
859 0 : Int MultiTermMatrixCleaner::computeHessianPeak()
860 : {
861 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeHessianPeak", WHERE));
862 0 : gip = IPosition(2,nx_p,ny_p);
863 :
864 0 : if(!doneCONV_p)
865 : {
866 : // Compute the convolutions of the smoothed psfs with each other.
867 : // Compute Assxx
868 : // Compute A100, A101, A102
869 : // A110, A111, A112
870 : // A120, A121, A122 for h(s1)
871 : // Compute A200, A201, A202
872 : // A210, A211, A212
873 : // A220, A221, A222 for h(s2)
874 : //... depending on the number of scales chosen
875 :
876 : // (PSF * scale) * (PSF * scale) -> cubeA_p [nx_p,ny_p,ntaylor,ntaylor,nscales]
877 0 : os << "Calculating PSF and Scale convolutions " << LogIO::POST;
878 0 : for (Int taylor1=0; taylor1<ntaylor_p;taylor1++)
879 0 : for (Int taylor2=0; taylor2<=taylor1;taylor2++)
880 0 : for (Int scale1=0; scale1<nscales_p;scale1++)
881 0 : for (Int scale2=0; scale2<=scale1;scale2++)
882 : {
883 0 : Int ttay1 = taylor1+taylor2;
884 :
885 : // CALC Hess : Calculate PSF_(t1+t2) * scale_1 * scale 2
886 :
887 0 : cWork_p.assign( (vecPsfFT_p[ttay1]) *(vecScalesFT_p[scale1])*(vecScalesFT_p[scale2]) );
888 :
889 0 : Bool usepatch=true;
890 0 : if(usepatch)
891 : {
892 0 : fftcomplex.fft0( vecWork_p[0] , cWork_p , false );
893 0 : Matrix<Float> psfpatch = ( vecWork_p[0] ) (itsPositionPeakPsf-psfsupport_p/2,itsPositionPeakPsf+psfsupport_p/2-IPosition(2,1,1));
894 0 : cubeA_p[IND4(taylor1,taylor2,scale1,scale2)] = psfpatch;
895 0 : }
896 : else
897 : {
898 0 : fftcomplex.fft0( cubeA_p[IND4(taylor1,taylor2,scale1,scale2)] , cWork_p , false );
899 : }
900 :
901 : //writeMatrixToDisk("psfconv_t_"+String::toString(taylor1)+"-"+String::toString(taylor2)+"_s_"+String::toString(scale1)+"-"+String::toString(scale2)+".im", cubeA_p[IND4(taylor1,taylor2,scale1,scale2)] );
902 : }
903 :
904 : // Construct A, invA for each scale.
905 0 : if(itsPositionPeakPsf != IPosition(2,(nx_p/2),(ny_p/2)))
906 0 : os << LogIO::WARN << "The PSF peak is not at the center of the image..." << LogIO::POST;
907 :
908 0 : Int stopnow=false;
909 0 : for (Int scale=0; scale<nscales_p;scale++)
910 : {
911 : // Fill up A
912 0 : for (Int taylor1=0; taylor1<ntaylor_p;taylor1++)
913 0 : for (Int taylor2=0; taylor2<ntaylor_p;taylor2++)
914 : {
915 : //(matA_p[scale])(taylor1,taylor2) = (cubeA_p[IND4(taylor1,taylor2,scale,scale)])(itsPositionPeakPsf); // OLD
916 0 : (matA_p[scale])(taylor1,taylor2) = (cubeA_p[IND4(taylor1,taylor2,scale,scale)])(psfpeak_p); // NEW
917 : /* Check for exact zeros ON MAIN DIAGONAL. Usually indicative of error */
918 0 : if( taylor1==taylor2 &&
919 0 : fabs( (matA_p[scale])(taylor1,taylor2) ) == 0.0 ) stopnow = true;
920 : }
921 :
922 0 : if(stopnow)
923 : {
924 0 : os << "Multi-Term Hessian has exact zeros on its main diagonal. Not proceeding further." << LogIO::WARN << endl;
925 0 : os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
926 0 : return -2;
927 : }
928 :
929 : /* If all elements are non-zero, check by brute-force if the rows/cols
930 : are nearly linearly dependant, making the matrix nearly singular. */
931 0 : Vector<Float> ratios(ntaylor_p);
932 0 : Float tsum=0.0;
933 0 : for(Int taylor1=0; taylor1<ntaylor_p-1; taylor1++)
934 : {
935 0 : for(Int taylor2=0; taylor2<ntaylor_p; taylor2++)
936 0 : ratios[taylor2] = (matA_p[scale])(taylor1,taylor2) / (matA_p[scale])(taylor1+1,taylor2);
937 0 : tsum=0.0;
938 0 : for(Int taylor2=0; taylor2<ntaylor_p-1; taylor2++)
939 0 : tsum += fabs(ratios[taylor2] - ratios[taylor2+1]);
940 0 : tsum /= ntaylor_p-1;
941 0 : if(tsum < 1e-04) stopnow=true;
942 :
943 : //cout << "Ratios for row " << taylor1 << " are " << ratios << endl;
944 : //cout << "tsum : " << tsum << endl;
945 : }
946 :
947 0 : if(stopnow)
948 : {
949 0 : os << "Multi-Term Hessian has linearly-dependent rows/cols. Not proceeding further." << LogIO::WARN << endl;
950 0 : os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
951 0 : return -2;
952 : }
953 :
954 : /* Try to invert the matrix. If it fails, return.
955 : The invertSymPosDef will check if it's positive definite or not.
956 : By construction, it should be pos-def. */
957 : // Compute inv(A)
958 : // Use MatrixMathLA::invert
959 : // or Use invertSymPosDef...
960 : //
961 : try
962 : {
963 0 : Double deter=0.0;
964 : //invert((*invMatA_p[scale]),deter,(*matA_p[scale]));
965 : //os << "Matrix Inverse : inv(A) : " << (*invMatA_p[scale]) << LogIO::POST;
966 :
967 : //if(adbg)
968 0 : os << "The Matrix [H] for " << scaleSizes_p[scale] << " pixel scale is : " << (matA_p[scale]) << LogIO::POST;
969 0 : invertSymPosDef((invMatA_p[scale]),deter,(matA_p[scale]));
970 0 : if(adbg) os << "Lapack Cholesky Decomposition Inverse of [A] is : " << (invMatA_p[scale]) << LogIO::POST;
971 : //if(adbg)os << "A matrix determinant : " << deter << LogIO::POST;
972 : //if(fabs(deter) < 1.0e-08) os << "SINGULAR MATRIX !! STOP!! " << LogIO::EXCEPTION;
973 : }
974 0 : catch(AipsError &x)
975 : {
976 0 : os << "The Matrix [A] is : " << (matA_p[scale]) << LogIO::POST;
977 0 : os << "Cannot Invert matrix : " << x.getMesg() << LogIO::WARN;
978 0 : return -2;
979 0 : }
980 :
981 : /* Add more checks based on the condition number of these matrices */
982 :
983 0 : }// end for scales
984 :
985 0 : doneCONV_p=true;
986 : }
987 :
988 0 : return 0;
989 0 : }/* end of computeHessianPeak() */
990 :
991 :
992 : /***************************************
993 : * Compute convolutions of the residual images ( all specs ) with scales.
994 : * --> the Right-Hand-Side of the matrix equation.
995 : ****************************************/
996 :
997 0 : Int MultiTermMatrixCleaner::computeRHS()
998 : {
999 : // LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeRHS()", WHERE));
1000 0 : IPosition blc1(2,0,0);
1001 0 : IPosition trc1(2,nx_p,ny_p);
1002 0 : IPosition inc1(2, 1);
1003 :
1004 : /* Compute R10 = I_D*B10, R11 = I_D*B11, R12 = I_D*B12
1005 : * Compute R20 = I_D*B20, R21 = I_D*B21, R22 = I_D*B22
1006 : * ... depending on the number of scales chosen.
1007 : */
1008 :
1009 : /* I_D * (PSF * scale) -> matR_p [nx_p,ny_p,ntaylor,nscales] */
1010 0 : for (Int taylor=0; taylor<ntaylor_p;taylor++)
1011 : {
1012 : /* Compute FT of dirty image */
1013 0 : fftcomplex.fft0( dirtyFT_p , vecDirty_p[taylor] , false );
1014 :
1015 0 : for (Int scale=0; scale<nscales_p;scale++)
1016 : {
1017 :
1018 : // CALC RHS : Calculate Dirty_t * scale_s
1019 :
1020 : // Let cWork get resized if needed. Force matR to have the right shape
1021 : //// cWork_p.assign( (dirtyFT_p)*(vecPsfFT_p[0])*(vecScalesFT_p[scale]) );
1022 0 : cWork_p.assign( (dirtyFT_p)*(vecScalesFT_p[scale]) );
1023 0 : fftcomplex.fft0( matR_p[IND2(taylor,scale)] , cWork_p , false );
1024 0 : fftcomplex.flip( matR_p[IND2(taylor,scale)] , false , false );
1025 : }
1026 : // writeMatrixToDisk("resid_"+String::toString(taylor)+".im", matR_p[IND2(taylor,0)] );
1027 : }
1028 :
1029 0 : return 0;
1030 0 : }/* end of computeRHS() */
1031 :
1032 : /***************************************
1033 : * Compute flux limit for minor cycles
1034 : ****************************************/
1035 :
1036 0 : Int MultiTermMatrixCleaner::computeFluxLimit(Float &fluxlimit, Float threshold)
1037 : {
1038 : // LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeFluxLimit", WHERE));
1039 :
1040 : // For now, since this calculation is being done outside.....
1041 0 : fluxlimit = threshold;
1042 :
1043 : /// fluxlimit = MAX(threshold , itsGain * max....
1044 :
1045 : // Later, implement the equivalent of an iteration-based cycle-speedup - if required.
1046 :
1047 0 : return 0;
1048 : }/* end of computeFluxLimit() */
1049 :
1050 :
1051 :
1052 : /***************************************
1053 : * Solve the matrix eqn for each point in the lattice.
1054 : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe
1055 : ****************************************/
1056 0 : Int MultiTermMatrixCleaner::solveMatrixEqn(Int ntaylor, Int scale, IPosition blc, IPosition trc)
1057 : {
1058 :
1059 0 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1060 : {
1061 0 : Matrix<Float> coeffs = (matCoeffs_p[IND2(taylor1,scale)])(blc,trc);
1062 0 : coeffs = 0.0;
1063 0 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1064 : {
1065 0 : Matrix<Float> rhs = (matR_p[IND2(taylor2,scale)])(blc,trc);
1066 0 : coeffs = coeffs + ((Float)(invMatA_p[scale])(taylor1,taylor2))* rhs;
1067 0 : }
1068 0 : }
1069 :
1070 : /* Solve for the coefficients, one scale at at time*/
1071 : /*
1072 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1073 : {
1074 : (matCoeffs_p[IND2(taylor1,scale)]) = 0.0;
1075 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1076 : {
1077 : matCoeffs_p[IND2(taylor1,scale)] = matCoeffs_p[IND2(taylor1,scale)] + ((Float)(invMatA_p[scale])(taylor1,taylor2))*(matR_p[IND2(taylor2,scale)]);
1078 : }
1079 : }
1080 : */
1081 0 : return 0;
1082 : }/* end of solveMatrixEqn() */
1083 :
1084 : /***************************************
1085 : * Find the 'peak' and its location. Fill this into maxScaleVal_p, maxScalePos_p
1086 : * Compute the update direction
1087 : * Options are
1088 : * (1) Original implementation : something related to the derivative of Chi-Sq. Not sure why it ever worked.
1089 : * (2) Peak residual in the 00 image
1090 : * (3) Peak t_0 component (for each scale)
1091 : * (4) Derivative of chi-square
1092 : * (5) Similar to (1), but with the PSF replaced by the peak of the PSF only.
1093 : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe
1094 : ****************************************/
1095 :
1096 0 : Int MultiTermMatrixCleaner::chooseComponent(Int ntaylor, Int scale, Int /*criterion*/, IPosition blc, IPosition trc)
1097 : {
1098 :
1099 :
1100 0 : Matrix<Float> work = ( vecWork_p[scale] )(blc,trc);
1101 :
1102 0 : work = 0.0;
1103 0 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1104 : {
1105 0 : Matrix<Float> coeffs1 = (matCoeffs_p[IND2(taylor1,scale)])(blc,trc);
1106 0 : Matrix<Float> resid = (matR_p[IND2(taylor1,scale)])(blc,trc);
1107 :
1108 : /*
1109 : work = work + (Float)2.0 * coeffs1 * resid;
1110 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1111 : {
1112 : Matrix<Float> coeffs2 = (matCoeffs_p[IND2(taylor2,scale)])(blc,trc);
1113 : work = work - (Float)((matA_p[scale])(taylor1,taylor2)) * coeffs1 * coeffs2;
1114 : }
1115 : */
1116 :
1117 0 : work = work + coeffs1 * resid;
1118 :
1119 0 : }
1120 0 : findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1121 :
1122 : /*
1123 :
1124 : vecWork_p[scale] = 0.0;
1125 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1126 : {
1127 : vecWork_p[scale] = vecWork_p[scale] + (Float)2.0 * ( (matCoeffs_p[IND2(taylor1,scale)]) * (matR_p[IND2(taylor1,scale)]) );
1128 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1129 : vecWork_p[scale] = vecWork_p[scale] - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
1130 : }
1131 : findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1132 : */
1133 :
1134 :
1135 : // switch(criterion)
1136 : // {
1137 : // case 1 : /* For each scale, find the maximum chi-square derivative (maybe) */
1138 : // {
1139 : /*
1140 : /// Code block using a private matrix
1141 : Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
1142 : ttWork_p = 0.0;
1143 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
1144 : {
1145 : ttWork_p = ttWork_p + (Float)2.0 * ( (matCoeffs_p[IND2(taylor1,scale)]) * (matR_p[IND2(taylor1,scale)]) );
1146 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
1147 : ttWork_p = ttWork_p - matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] * cubeA_p[IND4(taylor1,taylor2,scale,scale)];
1148 : }
1149 : findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1150 : */
1151 : /*
1152 : /// Code block using pre-allocated matrices
1153 : vecWork_p[scale] = 0.0;
1154 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
1155 : {
1156 : vecWork_p[scale] = vecWork_p[scale] + (Float)2.0 * ( (matCoeffs_p[IND2(taylor1,scale)]) * (matR_p[IND2(taylor1,scale)]) );
1157 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
1158 : vecWork_p[scale] = vecWork_p[scale] - matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] * cubeA_p[IND4(taylor1,taylor2,scale,scale)];
1159 : }
1160 : findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1161 : */
1162 : // }
1163 : // break;
1164 : /*
1165 : case 2 : // For each scale, find the peak residual
1166 : {
1167 : Float norm = sqrt((1.0/(matA_p[scale])(0,0)));
1168 : findMaxAbsMask( norm*matR_p[IND2(0,scale)] , vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1169 : }
1170 : break;
1171 : case 3: // For each scale, find the peak t_0 coefficnent
1172 : {
1173 : findMaxAbsMask( matCoeffs_p[IND2(0,scale)] , vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1174 : }
1175 : break;
1176 : case 4: // For each scale find the max chisq derivative (diag approx for all hessians) - BAD
1177 : {
1178 : Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
1179 : ttWork_p = 0.0;
1180 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
1181 : {
1182 : ttWork_p = ttWork_p + (Float)(2.0 * (matA_p[scale])(taylor1,taylor1) ) * (matR_p[IND2(taylor1,scale)]) ;
1183 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
1184 : ttWork_p = ttWork_p - (Float)(2.0*( (matA_p[scale])(taylor1,taylor1) * (matA_p[scale])(taylor1,taylor2))) * matCoeffs_p[IND2(taylor2,scale)] ;
1185 : }
1186 : findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1187 : }
1188 : break;
1189 : */
1190 : // case 5 : /* For each scale, same as 1, but use only psf peaks */
1191 : // {
1192 : /*
1193 : /// Code block using a private matrix
1194 : Matrix<Float> ttWork_p((matR_p[IND2(0,0)]).shape());
1195 : ttWork_p = 0.0;
1196 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
1197 : {
1198 : ttWork_p = ttWork_p + (Float)2.0 * ( (matCoeffs_p[IND2(taylor1,scale)]) * (matR_p[IND2(taylor1,scale)]) );
1199 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
1200 : ttWork_p = ttWork_p - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
1201 : }
1202 : findMaxAbsMask(ttWork_p,vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1203 : */
1204 : /*
1205 : /// Code block using pre-allocated matrices
1206 : vecWork_p[scale] = 0.0;
1207 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1208 : {
1209 : vecWork_p[scale] = vecWork_p[scale] + (Float)2.0 * ( (matCoeffs_p[IND2(taylor1,scale)]) * (matR_p[IND2(taylor1,scale)]) );
1210 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1211 : vecWork_p[scale] = vecWork_p[scale] - (Float)((matA_p[scale])(taylor1,taylor2)) * matCoeffs_p[IND2(taylor1,scale)] * matCoeffs_p[IND2(taylor2,scale)] ;
1212 : }
1213 : findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1214 :
1215 :
1216 : }
1217 : break;
1218 : */
1219 : /*
1220 : case 6 : // chi-square for this scale = sum of abs of residual images for taylor terms
1221 : {
1222 :
1223 : /// Code block using pre-allocated matrices
1224 : vecWork_p[scale] = 0.0;
1225 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1226 : {
1227 : vecWork_p[scale] = vecWork_p[scale] + ( abs(matR_p[IND2(taylor1,scale)]) );
1228 : }
1229 : findMaxAbsMask(vecWork_p[scale],vecScaleMasks_p[scale],maxScaleVal_p[scale],maxScalePos_p[scale]);
1230 :
1231 :
1232 : }
1233 : break;
1234 : */
1235 : // default:
1236 : // os << LogIO::SEVERE << "Internal error : Unknown option for type of update direction" << LogIO::POST;
1237 : // }
1238 :
1239 :
1240 0 : return 0;
1241 0 : }/* end of chooseComponent() */
1242 :
1243 : /* Update the RHS vector - Called from 'updateModelandRHS'.
1244 : Note : This function is called within the 'scale' omp/pragma loop. Needs to be thread-safe for scales.
1245 : */
1246 0 : Int MultiTermMatrixCleaner::updateRHS(Int ntaylor, Int scale, Float loopgain, Vector<Float> coeffs, IPosition blc, IPosition trc, IPosition blcPsf, IPosition trcPsf)
1247 : {
1248 0 : for(Int taylor1=0;taylor1<ntaylor;taylor1++)
1249 : {
1250 0 : Matrix<Float> residSub = (matR_p[IND2(taylor1,scale)])(blc,trc);
1251 0 : for(Int taylor2=0;taylor2<ntaylor;taylor2++)
1252 : {
1253 0 : Matrix<Float> smoothSub = (cubeA_p[IND4(taylor1,taylor2,scale,maxscaleindex_p)])(blcPsf,trcPsf);
1254 0 : residSub -= smoothSub * loopgain * coeffs[taylor2];
1255 : // residSub = residSub - smoothSub * loopgain * (matCoeffs_p[IND2(taylor2,maxscaleindex_p)])(globalmaxpos_p);
1256 0 : }
1257 0 : }
1258 : //
1259 0 : return 0;
1260 : }/* end of updateRHS */
1261 :
1262 :
1263 : /***************************************
1264 : * Update the model images and the convolved residuals
1265 : TODO : The current assumption is that only one scale is chosen at a time,
1266 : However the chi-sq derivative can be computed across scales too
1267 : and the update could be all scales at once for the 'best' location.
1268 : ****************************************/
1269 0 : Int MultiTermMatrixCleaner::updateModelAndRHS(Float loopgain)
1270 : {
1271 :
1272 : /* Max support size for all updates is the full image size */
1273 : // blc, trc :model image -> needs to be centred on the component location
1274 : // blcPsf : psf or scale-blob -> centred on the psf image center (peak).
1275 :
1276 : ///// gip = IPosition(2,nx_p,ny_p);
1277 :
1278 : //IPosition support(2,nx_p,ny_p); // OLD
1279 : //IPosition psfpeak(itsPositionPeakPsf);
1280 : /* The update region. */
1281 : //IPosition inc(2,1,1);
1282 : //IPosition blc(globalmaxpos_p - support/2);
1283 : //IPosition trc(globalmaxpos_p + support/2 - IPosition(2,1,1));
1284 : //LCBox::verify(blc, trc, inc, gip);
1285 : /* Shifted region, with the psf at the globalmaxpos_p. */
1286 : //IPosition blcPsf(blc + psfpeak - globalmaxpos_p); // OLD
1287 : //IPosition trcPsf(trc + psfpeak - globalmaxpos_p); // OLD
1288 : ///LCBox::verify(blcPsf, trcPsf, inc, gip); // OLD
1289 :
1290 :
1291 : /* The update region. */
1292 : /////IPosition inc(2,1,1);
1293 : /////IPosition blc(globalmaxpos_p - psfsupport_p/2);
1294 : ///// IPosition trc(globalmaxpos_p + psfsupport_p/2 - IPosition(2,1,1));
1295 : ///// LCBox::verify(blc, trc, inc, gip);
1296 :
1297 :
1298 : /* Shifted region, with the psf at the globalmaxpos_p. */
1299 : ///// IPosition blcPsf(blc + psfpeak_p - globalmaxpos_p); // OLD
1300 : ///// IPosition trcPsf(trc + psfpeak_p - globalmaxpos_p); // OLD
1301 : ///// LCBox::verify(blcPsf, trcPsf, inc, psfsupport_p); // NEW
1302 :
1303 : /* Reconcile box sizes/locations with the image size */
1304 : ///// makeBoxesSameSize(blc,trc,blcPsf,trcPsf);
1305 :
1306 : //buildImagePatches();
1307 :
1308 : //UUU
1309 : /*
1310 : cout << "Source location : " << globalmaxpos_p << endl;
1311 : cout << "region around peak residual : " << blc << trc << endl;
1312 : cout << "around the PSF peak : " << blcPsf << trcPsf << endl;
1313 : cout << "around the Scale blob : " << blcScale << trcScale << endl;
1314 : */
1315 :
1316 :
1317 : /* Update the model images */
1318 : /// Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcPsf,trcPsf); // OLD
1319 : /// Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcScale, trcScale); // NEW
1320 0 : Matrix<Float> scaleSub = (vecScales_p[maxscaleindex_p])(blcPsf_p, trcPsf_p); // NEWER (same size as psf)
1321 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1322 : {
1323 0 : Matrix<Float> modelSub = (vecModel_p[taylor])(blc_p,trc_p);
1324 0 : modelSub += scaleSub * loopgain * (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
1325 0 : }
1326 :
1327 : /* Update the convolved residuals */
1328 0 : Vector<Float> coeffs(ntaylor_p);
1329 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1330 0 : coeffs[taylor] = (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
1331 :
1332 : Int scale;
1333 0 : Int ntaylor=ntaylor_p;
1334 0 : IPosition blc(blc_p), trc(trc_p), blcPsf(blcPsf_p), trcPsf(trcPsf_p);
1335 : //OMP// #pragma omp parallel default(shared) private(scale) firstprivate(ntaylor,loopgain,coeffs,blc,trc,blcPsf,trcPsf)
1336 : {
1337 : //OMP// #pragma omp for
1338 0 : for(scale=0;scale<nscales_p;scale++)
1339 : {
1340 0 : updateRHS(ntaylor,scale, loopgain, coeffs, blc, trc, blcPsf, trcPsf);
1341 : }
1342 : }//End pragma parallel
1343 :
1344 : /* Update flux counters */
1345 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1346 : {
1347 0 : totalTaylorFlux_p[taylor] += loopgain*(matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p);
1348 : }
1349 0 : totalScaleFlux_p[maxscaleindex_p] += loopgain*(matCoeffs_p[IND2(0,maxscaleindex_p)])(globalmaxpos_p);
1350 :
1351 0 : return 0;
1352 0 : }/* end of updateModelAndRHS() */
1353 :
1354 : /* ................ */
1355 0 : Int MultiTermMatrixCleaner::checkConvergence(Int /*criterion*/, Float &fluxlimit, Float &loopgain)
1356 : {
1357 0 : Float rmaxval=0.0;
1358 :
1359 : /* Use the maximum residual (current), to compare against the convergence threshold */
1360 0 : Float maxres=0.0;
1361 0 : IPosition maxrespos;
1362 :
1363 0 : findMaxAbsMask((matR_p[IND2(0,0)]),vecScaleMasks_p[0],maxres,maxrespos);
1364 0 : Float norma = (1.0/(matA_p[0])(0,0));
1365 0 : rmaxval = abs(maxres*norma);
1366 0 : rmaxval_p = fabs(rmaxval);
1367 :
1368 : /* // Calc the max residual across all scales....
1369 : Int maxscale=0;
1370 : Float maxscaleresidual=0.0;
1371 : for (Int scale =0; scale<nscales_p; scale++)
1372 : {
1373 : findMaxAbsMask((matR_p[IND2(0,scale)]),vecScaleMasks_p[scale],maxres,maxrespos);
1374 : if ( maxscaleresidual < maxres )
1375 : {
1376 : maxscaleresidual = maxres;
1377 : maxscale = scale;
1378 : }
1379 : }
1380 : Float norma = (1.0/(matA_p[maxscale])(0,0));
1381 : rmaxval = abs(maxscaleresidual*norma);
1382 : */
1383 :
1384 : /* Check for convergence */
1385 0 : Int convergedflag = 0;
1386 : /// cout << "MTFT::checkconvergence : maxval : " << fabs(rmaxval) << " userthreshold : " << userthreshold_p << " fluxlimit : " << fluxlimit << endl;
1387 0 : if( fabs(rmaxval) < MAX(userthreshold_p,fluxlimit) )
1388 : {
1389 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
1390 0 : os << "Reached stopping threshold at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) ;
1391 0 : if( ! itsMask.null() ){os << " (within mask) " ;}
1392 0 : os << LogIO::POST;
1393 0 : convergedflag=1;
1394 0 : }
1395 :
1396 : /* Levenberg-Macquart-like change in step size */
1397 0 : if(itercount_p>1 && inputgain_p<=(Float)0.0)
1398 : {
1399 :
1400 0 : if (globalmaxval_p < prev_max_p)
1401 0 : loopgain = loopgain * 1.5;
1402 : else
1403 0 : loopgain = loopgain / 1.5;
1404 :
1405 0 : loopgain = MIN((1-stopfraction_p),loopgain);
1406 0 : loopgain = MIN((Float)0.6,loopgain);
1407 :
1408 : /* detect divergence by approximately 10 consecutive increases in maxval */
1409 0 : if(loopgain < (Float)0.01)
1410 : {
1411 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
1412 0 : os << "Not converging any more. May be diverging. Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << LogIO::POST;
1413 0 : convergedflag=-1;
1414 0 : }
1415 :
1416 : /* Stop if there is divergence : 200% increase from the current minimum value */
1417 0 : if( fabs( (min_max_p-globalmaxval_p)/min_max_p ) > (Float)2.0 )
1418 : {
1419 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
1420 0 : os << "Diverging.... Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << " Max: " << globalmaxval_p << LogIO::POST;
1421 0 : convergedflag=-1;
1422 0 : }
1423 :
1424 : // Stop if the maxval has changed by less than 5% in 5 iterations
1425 : // --- this is similar to saying less than 1% change per iteration.... same as a loopgain < 0.01
1426 :
1427 : // if( abs(prev5_max - globalmaxval_p) < 0.01*abs(prev5_max) )
1428 : // {
1429 : // os << "Not converging any more. Flattening out. Stopping minor cycles at iteration " << totalIters_p << ". Peak residual " << fabs(rmaxval) << LogIO::POST;
1430 : // convergedflag=-1;
1431 : // }
1432 :
1433 :
1434 :
1435 : }// end of if(itercount_p>1)
1436 :
1437 : /* Store current max value - to use in the next iteration */
1438 0 : prev_max_p = globalmaxval_p;
1439 : // if(itercount_p%5 == 0)
1440 : // prev5_max=globalmaxval_p;
1441 0 : min_max_p = MIN(min_max_p,abs(globalmaxval_p));
1442 :
1443 : /* Stop, if there are negatives on the largest scale in the Io image */
1444 : //if(nscales_p>1 && maxscaleindex_p == nscales_p-2)
1445 : // if((*matCoeffs_p[IND2(0,maxscaleindex_p)]).getAt(globalmaxpos_p) < 0.0)
1446 : // {converged = false;break;}
1447 :
1448 : /* Detect divergence, and signal it.... */
1449 : // TODO
1450 :
1451 : /* Print out coefficients for a few iterations */
1452 0 : if(convergedflag==0)
1453 : {
1454 0 : if(fluxlimit==-1)
1455 : {
1456 0 : fluxlimit = rmaxval * stopfraction_p;
1457 : // os << "Peak convolved residual : " << rmaxval << " Minor cycle stopping threshold : " << itsThreshold.getValue("Jy") << LogIO::POST;
1458 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "mtclean()", WHERE));
1459 0 : os << "Peak convolved residual" ;
1460 0 : if( ! itsMask.null() ){os << " (within mask) ";}
1461 0 : os << " : " << rmaxval;
1462 0 : if( fluxlimit > 0.0 ){ os << " : Minor cycle stopping threshold : " << fluxlimit;}
1463 0 : os << LogIO::POST;
1464 0 : }
1465 : else
1466 : {
1467 : // if(1)
1468 0 : if( totalIters_p==maxniter_p || (adbg==(Bool)true) || maxniter_p < (int)5 || (totalIters_p%(Int)20==0) )
1469 : {
1470 :
1471 0 : os << "[" << totalIters_p << "] Res: " << rmaxval << " Max: " << globalmaxval_p;
1472 0 : os << " Gain: " << loopgain;
1473 : //// os << "[" << totalIters_p << "] Res: " << rmaxval;
1474 : //os << "[" << totalIters_p << "] Max: " << globalmaxval_p;
1475 0 : os << " Pos: " << globalmaxpos_p << " Scale: " << scaleSizes_p[maxscaleindex_p];
1476 0 : os << " Coeffs: ";
1477 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1478 0 : os << (matCoeffs_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p) << " ";
1479 0 : if(adbg)
1480 : {
1481 0 : os << " OrigRes: ";
1482 0 : for(Int taylor=0;taylor<ntaylor_p;taylor++)
1483 0 : os << (matR_p[IND2(taylor,maxscaleindex_p)])(globalmaxpos_p) << " ";
1484 : }
1485 0 : os << LogIO::POST;
1486 :
1487 : //UUU cout << "minor " << totalIters_p << " " << (vecModel_p[0])( IPosition(2,nx_p/2, ny_p/2) ) << " " << (vecModel_p[1])( IPosition(2,nx_p/2, ny_p/2) ) << " " << (matR_p[0])( IPosition(2,nx_p/2, ny_p/2) ) << endl;
1488 :
1489 : }
1490 : }
1491 :
1492 : }// if converged-flag is still 0
1493 :
1494 :
1495 0 : return convergedflag;
1496 :
1497 0 : }/* end of checkConvergence */
1498 :
1499 :
1500 :
1501 : /* Save a matrix to disk */
1502 0 : Int MultiTermMatrixCleaner::writeMatrixToDisk(String imagename, Matrix<Float>& themat)
1503 : {
1504 0 : TabularCoordinate tab1(0, 1.0, 0, String("deg"), String("axis1"));
1505 0 : TabularCoordinate tab2(0, 1.0, 0, String("deg"), String("axis2"));
1506 0 : CoordinateSystem csys;
1507 0 : csys.addCoordinate(tab1);
1508 0 : csys.addCoordinate(tab2);
1509 0 : PagedImage<Float> limage(themat.shape(), csys, imagename);
1510 0 : limage.put(themat);
1511 : //return value does not seemed to be used so will make compiler happy
1512 0 : return 1;
1513 0 : }
1514 :
1515 :
1516 : /* Compute principal solution in-place on the list of residual images ( vecDirty )
1517 : -- Call MTMC::setresidual() repeatedly to fill in final residuals before computing
1518 : the principal solution. This does the same as solveMatrixEquation(), but
1519 : stores the results in-place in the residual images */
1520 0 : Bool MultiTermMatrixCleaner::computeprincipalsolution()
1521 : {
1522 0 : LogIO os(LogOrigin("MultiTermMatrixCleaner", "computeprincipalsolution()", WHERE));
1523 :
1524 0 : os << "MTMC :: Computing principal solution on residuals" << LogIO::POST;
1525 :
1526 : /// If this is being called with niter=0, the Hessian won't exist. So, make it.
1527 0 : if( doneCONV_p == false )
1528 : {
1529 0 : if( computeHessianPeak() == -2 )
1530 0 : return false;
1531 : }
1532 :
1533 0 : AlwaysAssert((vecDirty_p.nelements()>0), AipsError);
1534 :
1535 : /* Solve for the coefficients at the delta-function scale*/
1536 0 : for(Int taylor1=0;taylor1<ntaylor_p;taylor1++)
1537 : {
1538 0 : (matCoeffs_p[IND2(taylor1,0)]) = 0.0;
1539 0 : for(Int taylor2=0;taylor2<ntaylor_p;taylor2++)
1540 : {
1541 0 : matCoeffs_p[IND2(taylor1,0)] = matCoeffs_p[IND2(taylor1,0)] + ((Float)(invMatA_p[0])(taylor1,taylor2))*(vecDirty_p[taylor2]);
1542 : }
1543 : }
1544 :
1545 : /* Copy this into the residual vector */
1546 0 : for(Int taylor=0; taylor<ntaylor_p;taylor++)
1547 : {
1548 0 : vecDirty_p[taylor].assign(matCoeffs_p[IND2(taylor,0)]);
1549 : }
1550 :
1551 0 : return true;
1552 0 : }
1553 :
1554 :
1555 :
1556 : } //# NAMESPACE CASA - END
|
