Line data Source code
1 : //# SDAlgorithmMSMFS.cc: Implementation of SDAlgorithmMSMFS classes
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,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 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <synthesis/ImagerObjects/SDAlgorithmMSMFS.h>
31 : #include <components/ComponentModels/SkyComponent.h>
32 : #include <components/ComponentModels/ComponentList.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
39 : #include <casacore/lattices/Lattices/LatticeStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <synthesis/TransformMachines/StokesImageUtil.h>
42 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 : #include <casacore/casa/Utilities/Assert.h>
46 : #include <casacore/casa/OS/Directory.h>
47 : #include <casacore/tables/Tables/TableLock.h>
48 : #include <imageanalysis/ImageAnalysis/CasaImageBeamSet.h>
49 :
50 : #include<synthesis/ImagerObjects/SIMinorCycleController.h>
51 :
52 : #include <sstream>
53 :
54 : #include <casacore/casa/Logging/LogMessage.h>
55 : #include <casacore/casa/Logging/LogIO.h>
56 : #include <casacore/casa/Logging/LogSink.h>
57 :
58 : #include <casacore/casa/System/Choice.h>
59 : #include <msvis/MSVis/StokesVector.h>
60 :
61 :
62 : using namespace casacore;
63 : namespace casa { //# NAMESPACE CASA - BEGIN
64 :
65 :
66 141 : SDAlgorithmMSMFS::SDAlgorithmMSMFS( uInt nTaylorTerms, Vector<Float> scalesizes, Float smallscalebias):
67 : SDAlgorithmBase(),
68 : // itsImages(),
69 : itsMatPsfs(), itsMatResiduals(), itsMatModels(),
70 : itsNTerms(nTaylorTerms),
71 : itsScaleSizes(scalesizes),
72 : itsSmallScaleBias(smallscalebias),
73 : itsMTCleaner(),
74 141 : itsMTCsetup(false)
75 : {
76 141 : itsAlgorithmName=String("mtmfs");
77 141 : if( itsScaleSizes.nelements()==0 ){ itsScaleSizes.resize(1); itsScaleSizes[0]=0.0; }
78 141 : }
79 :
80 282 : SDAlgorithmMSMFS::~SDAlgorithmMSMFS()
81 : {
82 :
83 282 : }
84 :
85 :
86 : // void SDAlgorithmMSMFS::initializeDeconvolver( Float &peakresidual, Float &modelflux )
87 129 : void SDAlgorithmMSMFS::initializeDeconvolver()
88 : {
89 387 : LogIO os( LogOrigin("SDAlgorithmMSMFS","initializeDeconvolver",WHERE) );
90 :
91 129 : AlwaysAssert( (bool) itsImages, AipsError );
92 129 : AlwaysAssert( itsNTerms == itsImages->getNTaylorTerms() , AipsError );
93 :
94 129 : itsMatPsfs.resize( 2*itsNTerms-1 );
95 129 : itsMatResiduals.resize( itsNTerms );
96 129 : itsMatModels.resize( itsNTerms );
97 :
98 : //// Why is this needed ? I hope this is by reference.
99 500 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
100 : {
101 377 : if(tix<itsNTerms)
102 : {
103 256 : (itsImages->residual(tix))->get( itsMatResiduals[tix], true );
104 254 : (itsImages->model(tix))->get( itsMatModels[tix], true );
105 : }
106 373 : (itsImages->psf(tix))->get( itsMatPsfs[tix], true );
107 : }
108 :
109 123 : itsImages->mask()->get( itsMatMask, true );
110 :
111 : //// Initialize the MultiTermMatrixCleaner.
112 :
113 : /// ----------- do once ----------
114 123 : if( itsMTCsetup == false)
115 : {
116 : //cout << "Setting up the MT Cleaner once" << endl;
117 : //Vector<Float> scalesizes(1); scalesizes[0]=0.0;
118 89 : itsMTCleaner.setscales( itsScaleSizes );
119 :
120 89 : if(itsSmallScaleBias > 1)
121 : {
122 0 : os << LogIO::WARN << "Acceptable smallscalebias values are [-1,1].Changing smallscalebias from " << itsSmallScaleBias <<" to 1." << LogIO::POST;
123 0 : itsSmallScaleBias = 1;
124 : }
125 :
126 89 : if(itsSmallScaleBias < -1)
127 : {
128 0 : os << LogIO::WARN << "Acceptable smallscalebias values are [-1,1].Changing smallscalebias from " << itsSmallScaleBias <<" to -1." << LogIO::POST;
129 0 : itsSmallScaleBias = -1;
130 : }
131 :
132 :
133 89 : itsMTCleaner.setSmallScaleBias(itsSmallScaleBias);
134 89 : itsMTCleaner.setntaylorterms( itsNTerms );
135 89 : itsMTCleaner.initialise( itsImages->getShape()[0], itsImages->getShape()[1] );
136 :
137 352 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
138 : {
139 526 : Matrix<Float> tempMat;
140 263 : tempMat.reference( itsMatPsfs[tix] );
141 263 : itsMTCleaner.setpsf( tix, tempMat );
142 : /// itsMTCleaner.setpsf( tix, itsMatPsfs[tix] );
143 : }
144 89 : itsMTCsetup=true;
145 : }
146 : /// -----------------------------------------
147 :
148 : /*
149 : /// Find initial max vals..
150 : findMaxAbsMask( itsMatResiduals[0], itsMatMask, itsPeakResidual, itsMaxPos );
151 : itsModelFlux = sum( itsMatModels[0] );
152 :
153 : peakresidual = itsPeakResidual;
154 : modelflux = itsModelFlux;
155 : */
156 :
157 : // Parts to be repeated at each minor cycle start....
158 :
159 246 : Matrix<Float> tempmask(itsMatMask);
160 123 : itsMTCleaner.setmask( tempmask );
161 :
162 367 : for(uInt tix=0; tix<itsNTerms; tix++)
163 : {
164 488 : Matrix<Float> tempMat;
165 244 : tempMat.reference( itsMatResiduals[tix] );
166 244 : itsMTCleaner.setresidual( tix, tempMat );
167 : // itsMTCleaner.setresidual( tix, itsMatResiduals[tix] );
168 :
169 488 : Matrix<Float> tempMat2;
170 244 : tempMat2.reference( itsMatModels[tix] );
171 244 : itsMTCleaner.setmodel( tix, tempMat2 );
172 : // itsMTCleaner.setmodel( tix, itsMatModels[tix] );
173 : }
174 :
175 : // Print some useful info about the frequency setup for this MTMFS run
176 123 : itsImages->calcFractionalBandwidth();
177 123 : }
178 :
179 130 : Long SDAlgorithmMSMFS::estimateRAM(const vector<int>& imsize){
180 : ///taken from imager_resource_predictor.py
181 :
182 130 : Long mem=0;
183 260 : IPosition shp;
184 130 : if(itsImages){
185 0 : shp=itsImages->getShape();
186 : }
187 130 : else if(imsize.size() >1){
188 130 : shp=IPosition(imsize);
189 : }
190 : else
191 0 : return 0;
192 : //throw(AipsError("Deconvolver cannot estimate the memory usage at this point"));
193 :
194 : // psf patches in the Hessian
195 130 : Long nscales=itsScaleSizes.nelements();
196 130 : Long n4d = (nscales * (nscales+1) / 2.0) * (itsNTerms * (itsNTerms+1) / 2);
197 :
198 130 : Long nsupport = Long(Float(100.0/Float(shp(0)))* Float(100.0/Float(shp(1))* Float(n4d + nscales)));
199 :
200 130 : Long nfull = 2 + 2 + 3 * nscales + 3 * itsNTerms + (2 * itsNTerms - 1) + 2 * itsNTerms * nscales;
201 130 : Long nfftserver = 1 + 2*2 ; /// 1 float and 2 complex
202 :
203 130 : Long mystery = 1 + 1; /// TODO
204 :
205 130 : Long ntotal = nsupport + nfull + nfftserver + mystery;
206 130 : mem=sizeof(Float)*(shp(0))*(shp(1))*ntotal/1024;
207 :
208 130 : return mem;
209 : }
210 :
211 123 : void SDAlgorithmMSMFS::takeOneStep( Float loopgain, Int cycleNiter, Float cycleThreshold, Float &peakresidual, Float &modelflux, Int &iterdone)
212 : {
213 :
214 123 : iterdone = itsMTCleaner.mtclean( cycleNiter, 0.0, loopgain, cycleThreshold );
215 :
216 123 : if( iterdone==-2 ) throw(AipsError("MT-Cleaner error : Non-invertible Hessian. Please check if the multi-frequency data selection is appropriate for a polynomial fit of the desired order."));
217 :
218 367 : for(uInt tix=0; tix<itsNTerms; tix++)
219 : {
220 488 : Matrix<Float> tempMat;
221 244 : tempMat.reference( itsMatModels[tix] );
222 :
223 244 : itsMTCleaner.getmodel( tix, tempMat ); //itsMatModels[tix] );
224 : }
225 :
226 : /////////////////
227 : //findMaxAbs( itsMatResiduals[0], itsPeakResidual, itsMaxPos );
228 : //peakresidual = itsPeakResidual;
229 :
230 123 : peakresidual = itsMTCleaner.getpeakresidual();
231 : // cout << "Peak res from matR : " << peakresidual << endl; // Uncomment for debugging
232 :
233 :
234 : // Retrieve residual to be saved to the .residual file in finalizeDeconvolver
235 367 : for(uInt tix=0; tix<itsNTerms; tix++)
236 : {
237 488 : casacore::Matrix<Float> tmp;
238 244 : itsMTCleaner.getresidual(tix, tmp); // possible room for optimization here -> get residual without extra tmp copy? maybe change getResidual to accept an array?
239 244 : itsMatResiduals[tix] = tmp;
240 : }
241 :
242 123 : peakresidual = max(abs(itsMatResiduals[0]*itsMatMask));
243 : // cout << "Peak res from new math : " << peakresidual << endl; // Uncomment for debugging
244 123 : modelflux = sum( itsMatModels[0] ); // Performance hog ?
245 :
246 123 : }
247 :
248 123 : void SDAlgorithmMSMFS::finalizeDeconvolver()
249 : {
250 : // itsResidual.put( itsMatResidual );
251 : // itsModel.put( itsMatModel );
252 :
253 : // Why is this needed ? If the matrices are by reference, then why do we need this ?
254 367 : for(uInt tix=0; tix<itsNTerms; tix++)
255 : {
256 244 : (itsImages->residual(tix))->put( itsMatResiduals[tix] );
257 244 : (itsImages->model(tix))->put( itsMatModels[tix] );
258 : }
259 123 : }
260 :
261 120 : void SDAlgorithmMSMFS::restore(std::shared_ptr<SIImageStore> imagestore)
262 : {
263 :
264 240 : LogIO os( LogOrigin("SDAlgorithmMSMFS","restore",WHERE) );
265 :
266 120 : if( ! imagestore->hasResidualImage() ) return;
267 :
268 : // Compute principal solution ( if it hasn't already been done to this ImageStore...... )
269 : ////// Put some image misc info in here, to say if it has been done or not.
270 :
271 : // Loop over polarization planes here, as MTC knows only about Matrices.
272 : Int nSubChans, nSubPols;
273 :
274 119 : queryDesiredShape(nSubChans, nSubPols, imagestore->getShape());
275 :
276 : // CAS-13401 : Store restoring beam per plane, so it can be set in the final restored image.
277 238 : ImageBeamSet restoringBeams;
278 119 : restoringBeams.resize(nSubChans, nSubPols);
279 :
280 238 : for( Int chanid=0; chanid<nSubChans;chanid++) // redundant since only 1 chan
281 : {
282 :
283 244 : for( Int polid=0; polid<nSubPols;polid++) // one pol plane at a time
284 : {
285 125 : itsImages = imagestore->getSubImageStore( 0, 1, chanid, nSubChans, polid, nSubPols );
286 :
287 : //cout << "Units for chan " << chanid << " and pol " << polid << " are " << itsImages->image()->units().getName() << endl;
288 :
289 : /// ----------- do once if trying to 'only restore' without model ----------
290 125 : if( itsMTCsetup == false)
291 : {
292 :
293 42 : itsMatPsfs.resize( 2*itsNTerms-1 );
294 166 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
295 : {
296 124 : (itsImages->psf(tix))->get( itsMatPsfs[tix], True );
297 : }
298 :
299 : //cout << "Setting up the MT Cleaner once" << endl;
300 : //Vector<Float> scalesizes(1); scalesizes[0]=0.0;
301 42 : itsMTCleaner.setscales( itsScaleSizes );
302 42 : itsMTCleaner.setntaylorterms( itsNTerms );
303 42 : itsMTCleaner.initialise( itsImages->getShape()[0], itsImages->getShape()[1] );
304 :
305 166 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
306 : {
307 248 : Matrix<Float> tempMat;
308 124 : tempMat.reference( itsMatPsfs[tix] );
309 124 : itsMTCleaner.setpsf( tix, tempMat );
310 : /// itsMTCleaner.setpsf( tix, itsMatPsfs[tix] );
311 : }
312 42 : itsMTCsetup=true;
313 : }
314 : /// -----------------------------------------
315 :
316 :
317 250 : Vector<TempImage<Float> > tempResOrig(itsNTerms);
318 :
319 : // Set residual images into mtcleaner
320 372 : for(uInt tix=0; tix<itsNTerms; tix++)
321 : {
322 494 : Array<Float> tempArr;
323 247 : (itsImages->residual(tix))->get( tempArr, True );
324 247 : Matrix<Float> tempMat;
325 247 : tempMat.reference( tempArr );
326 247 : itsMTCleaner.setresidual( tix, tempMat );
327 :
328 : // Also save them temporarily (copies)
329 247 : tempResOrig[tix] = TempImage<Float>(itsImages->getShape(), itsImages->residual(tix)->coordinates());
330 247 : tempResOrig[tix].copyData( LatticeExpr<Float>(* ( itsImages->residual(tix) ) ) );
331 : }
332 :
333 : // Modify the original in place
334 125 : itsMTCleaner.computeprincipalsolution();
335 :
336 372 : for(uInt tix=0; tix<itsNTerms; tix++)
337 : {
338 247 : Matrix<Float> tempRes;
339 247 : itsMTCleaner.getresidual(tix,tempRes);
340 247 : (itsImages->residual(tix))->put( tempRes );
341 : }
342 :
343 : // Calculate restored image and alpha using modified residuals
344 125 : SDAlgorithmBase::restore( itsImages );
345 :
346 : // This is required because imageInfo() only contains the beam for this chan/pol.
347 250 : GaussianBeam thisBeam = itsImages->image(0)->imageInfo().getBeamSet().getBeam(chanid, polid);
348 125 : restoringBeams.setBeam(chanid, polid, thisBeam);
349 :
350 : // Put back original unmodified residuals.o
351 372 : for(uInt tix=0; tix<itsNTerms; tix++) {
352 247 : (itsImages->residual(tix))->copyData( LatticeExpr<Float>( tempResOrig(tix) ) );
353 : }
354 : } // for polid loop
355 : }// for chanid loop
356 :
357 :
358 : // This log message is important. This call of imagestore->image(...) is the first call if there is
359 : // a multi-channel or multi-pol image. This is what will set the units correctly. Ref. CAS-13153
360 119 : os << LogIO::POST << "Restored images : ";
361 354 : for(uInt tix=0; tix<itsNTerms; tix++)
362 : {
363 235 : os << LogIO::POST << imagestore->image(tix)->name() << " (model=" << imagestore->model(tix)->name() << ") " ;
364 :
365 : // CAS-13401 : Set the per-chan per-pol beam info for the restored image.
366 235 : ImageInfo iminf = imagestore->image(tix)->imageInfo();
367 :
368 235 : iminf.removeRestoringBeam();
369 :
370 : // If restoringbeam="common", then calculate and only set one beam
371 235 : if (itsRestoringBeam.isNull() && itsUseBeam == "common"){
372 0 : GaussianBeam cbeam = CasaImageBeamSet(restoringBeams).getCommonBeam();
373 0 : iminf.setRestoringBeam(cbeam);
374 : }
375 235 : else if (! itsRestoringBeam.isNull()) {
376 2 : iminf.setRestoringBeam(itsRestoringBeam);
377 : }
378 : else {
379 233 : iminf.setBeams(restoringBeams);
380 : }
381 235 : imagestore->image(tix)->setImageInfo(iminf);
382 :
383 : }
384 119 : os << LogIO::POST << endl;
385 :
386 : //cout << "Units for " << imagestore->image()->name() << " aaaaaare " << imagestore->image()->units().getName() << endl;
387 : }// ::restore
388 :
389 :
390 : /*
391 : void SDAlgorithmMSMFS::restorePlane()
392 : {
393 :
394 : LogIO os( LogOrigin("SDAlgorithmMSMFS","restorePlane",WHERE) );
395 :
396 : try
397 : {
398 : // Fit a Gaussian to the PSF.
399 : GaussianBeam beam = itsImages->getPSFGaussian();
400 :
401 : os << "Restore with beam : "
402 : << beam.getMajor(Unit("arcmin")) << " arcmin, "
403 : << beam.getMinor(Unit("arcmin"))<< " arcmin, "
404 : << beam.getPA(Unit("deg")) << " deg" << LogIO::POST;
405 :
406 : // Compute principal solution ( if it hasn't already been done to this ImageStore...... )
407 : itsMTCleaner.computeprincipalsolution();
408 : for(uInt tix=0; tix<itsNTerms; tix++)
409 : {
410 : Matrix<Float> tempRes;
411 : itsMTCleaner.getresidual(tix,tempRes);
412 : (itsImages->residual(tix))->put( tempRes );
413 : }
414 :
415 : // Restore all terms
416 : ImageInfo ii = itsImages->image(0)->imageInfo();
417 : ii.setRestoringBeam( beam );
418 :
419 : for(uInt tix=0; tix<itsNTerms; tix++)
420 : {
421 : (itsImages->image(tix))->set(0.0);
422 : (itsImages->image(tix))->copyData( LatticeExpr<Float>(*(itsImages->model(tix))) );
423 : StokesImageUtil::Convolve( *(itsImages->image(tix)) , beam);
424 : (itsImages->image(tix))->copyData( LatticeExpr<Float>
425 : ( *(itsImages->model(tix)) + *(itsImages->residual(tix)) ) );
426 : itsImages->image()->setImageInfo(ii);
427 : }
428 :
429 : // Calculate alpha and beta
430 : LatticeExprNode leMaxRes = max( *( itsImages->residual(0) ) );
431 : Float maxres = leMaxRes.getFloat();
432 : Float specthreshold = maxres/5.0; //////////// do something better here.....
433 :
434 : os << "Calculating spectral parameters for Intensity > peakresidual/5 = " << specthreshold << " Jy/beam" << LogIO::POST;
435 : LatticeExpr<Float> mask1(iif(((*(itsImages->image(0))))>(specthreshold),1.0,0.0));
436 : LatticeExpr<Float> mask0(iif(((*(itsImages->image(0))))>(specthreshold),0.0,1.0));
437 :
438 : /////// Calculate alpha
439 : LatticeExpr<Float> alphacalc( (((*(itsImages->image(1))))*mask1)/(((*(itsImages->image(0))))+(mask0)) );
440 : itsImages->alpha()->copyData(alphacalc);
441 :
442 : // Set the restoring beam for alpha
443 : itsImages->alpha()->setImageInfo(ii);
444 : // itsImages->alpha()->table().unmarkForDelete();
445 :
446 : // Make a mask for the alpha image
447 : LatticeExpr<Bool> lemask(iif(((*(itsImages->image(0))) > specthreshold) , true, false));
448 :
449 : createMask(lemask, *(itsImages->alpha()));
450 :
451 : }
452 : catch(AipsError &x)
453 : {
454 : throw( AipsError("Restoration Error : " + x.getMesg() ) );
455 : }
456 :
457 : }
458 :
459 :
460 : Bool SDAlgorithmMSMFS::createMask(LatticeExpr<Bool> &lemask, ImageInterface<Float> &outimage)
461 : {
462 : ImageRegion outreg = outimage.makeMask("mask0",false,true);
463 : LCRegion& outmask=outreg.asMask();
464 : outmask.copyData(lemask);
465 : outimage.defineRegion("mask0",outreg, RegionHandler::Masks, true);
466 : outimage.setDefaultMask("mask0");
467 : return true;
468 : }
469 :
470 : */
471 :
472 : } //# NAMESPACE CASA - END
473 :
|
