Line data Source code
1 : //# SIImageStoreMultiTerm.cc: Implementation of Imager.h
2 : //# Copyright (C) 1997-2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 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/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 : #include <casacore/images/Images/TempImage.h>
48 : #include <casacore/images/Images/PagedImage.h>
49 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
50 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
51 : #include <synthesis/TransformMachines/StokesImageUtil.h>
52 : #include <casacore/images/Images/TempImage.h>
53 : #include <casacore/images/Images/SubImage.h>
54 : #include <casacore/images/Regions/ImageRegion.h>
55 :
56 : #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
57 :
58 : #include <casacore/casa/Arrays/MatrixMath.h>
59 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
60 :
61 : #include <sys/types.h>
62 : #include <unistd.h>
63 : using namespace std;
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 : //////////////////////////////////////////////////////////////////////////////////////////////////////
69 : //////////////////////////////////////////////////////////////////////////////////////////////////////
70 :
71 0 : SIImageStoreMultiTerm::SIImageStoreMultiTerm():SIImageStore()
72 : {
73 0 : itsNTerms=0;
74 :
75 0 : itsPsfs.resize(0);
76 0 : itsModels.resize(0);
77 0 : itsResiduals.resize(0);
78 0 : itsWeights.resize(0);
79 0 : itsImages.resize(0);
80 0 : itsSumWts.resize(0);
81 0 : itsImagePBcors.resize(0);
82 0 : itsPBs.resize(0);
83 :
84 0 : itsForwardGrids.resize(0);
85 0 : itsBackwardGrids.resize(0);
86 :
87 0 : itsUseWeight=false;
88 :
89 0 : init();
90 :
91 0 : validate();
92 :
93 0 : }
94 :
95 : // Used from SynthesisNormalizer::makeImageStore()
96 116 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename,
97 : const CoordinateSystem &imcoordsys,
98 : const IPosition &imshape,
99 : const String &objectname,
100 : const Record &miscinfo,
101 : const Int /*nfacets*/,
102 : const Bool /*overwrite*/,
103 : uInt ntaylorterms,
104 116 : Bool useweightimage)
105 : {
106 348 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open new Images",WHERE) );
107 :
108 116 : itsNTerms = ntaylorterms;
109 :
110 116 : itsPsfs.resize(2 * itsNTerms - 1);
111 116 : itsModels.resize(itsNTerms);
112 116 : itsResiduals.resize(itsNTerms);
113 116 : itsWeights.resize(2 * itsNTerms - 1);
114 116 : itsImages.resize(itsNTerms);
115 116 : itsSumWts.resize(2 * itsNTerms - 1);
116 116 : itsPBs.resize(itsNTerms);
117 116 : itsImagePBcors.resize(itsNTerms);
118 :
119 116 : itsMask.reset( );
120 116 : itsGridWt.reset( );
121 :
122 116 : itsForwardGrids.resize(itsNTerms);
123 116 : itsBackwardGrids.resize(2 * itsNTerms - 1);
124 :
125 : // itsNFacets = nfacets; // So that sumwt shape happens properly, via checkValidity
126 : // itsFacetId = -1;
127 116 : itsNFacets=1;
128 116 : itsFacetId=0;
129 116 : itsNChanChunks = 1;
130 116 : itsChanId = 0;
131 116 : itsNPolChunks = 1;
132 116 : itsPolId = 0;
133 :
134 116 : itsImageName = imagename;
135 116 : itsCoordSys = imcoordsys;
136 116 : itsImageShape = imshape;
137 116 : itsObjectName = objectname;
138 116 : itsMiscInfo = miscinfo;
139 :
140 116 : itsUseWeight = useweightimage;
141 :
142 116 : init();
143 :
144 116 : validate();
145 :
146 116 : }
147 :
148 : // Used from SynthesisNormalizer::makeImageStore()
149 603 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const String &imagename, uInt ntaylorterms,
150 603 : const Bool ignorefacets, const Bool ignoresumwt)
151 : {
152 1809 : LogIO os( LogOrigin("SIImageStoreMultiTerm","Open existing Images",WHERE) );
153 :
154 603 : itsNTerms = ntaylorterms;
155 :
156 16105 : auto fltPtrReset = [](Block<shared_ptr<ImageInterface<Float> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
157 603 : itsPsfs.resize(2 * itsNTerms - 1); fltPtrReset(itsPsfs);
158 603 : itsModels.resize(itsNTerms); fltPtrReset(itsModels);
159 603 : itsResiduals.resize(itsNTerms); fltPtrReset(itsResiduals);
160 603 : itsWeights.resize(2 * itsNTerms - 1); fltPtrReset(itsWeights);
161 603 : itsImages.resize(itsNTerms); fltPtrReset(itsImages);
162 603 : itsPBs.resize(itsNTerms); fltPtrReset(itsPBs);
163 603 : itsSumWts.resize(2 * itsNTerms - 1); fltPtrReset(itsSumWts);
164 603 : itsMask.reset( );
165 603 : itsGridWt.reset( );
166 603 : itsImagePBcors.resize(itsNTerms); fltPtrReset(itsImagePBcors);
167 :
168 :
169 :
170 603 : itsMiscInfo=Record();
171 4173 : auto cmplxPtrReset = [](Block<shared_ptr<ImageInterface<Complex> > >&a) {for(uInt i=0; i < a.nelements(); ++i) a[i].reset(); };
172 603 : itsForwardGrids.resize(itsNTerms); cmplxPtrReset(itsForwardGrids);
173 603 : itsBackwardGrids.resize(2 * itsNTerms - 1); cmplxPtrReset(itsBackwardGrids);
174 :
175 603 : itsImageName = imagename;
176 :
177 603 : itsNFacets=1;
178 603 : itsFacetId=0;
179 603 : itsNChanChunks = 1;
180 603 : itsChanId = 0;
181 603 : itsNPolChunks = 1;
182 603 : itsPolId = 0;
183 :
184 603 : Bool exists=true;
185 603 : Bool sumwtexists=true;
186 603 : Bool modelexists=true;
187 2380 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
188 : {
189 1777 : if( tix<itsNTerms ) {
190 2862 : exists &= ( doesImageExist( itsImageName+String(".residual.tt")+String::toString(tix) ) ||
191 1672 : doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
192 2966 : modelexists &= ( doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) ||
193 1776 : doesImageExist( itsImageName+String(".model.tt")+String::toString(tix) ) );
194 1190 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
195 : }else {
196 587 : exists &= ( doesImageExist( itsImageName+String(".psf.tt")+String::toString(tix) ) );
197 587 : sumwtexists &= ( doesImageExist( itsImageName+String(".sumwt.tt")+String::toString(tix) ) );
198 : }
199 : }
200 :
201 : // The PSF or Residual images must exist. ( or the gridwt image)
202 : // All this is just for the shape and coordinate system.
203 603 : if( exists || modelexists || doesImageExist(itsImageName+String(".gridwt")) )
204 : {
205 0 : std::shared_ptr<ImageInterface<Float> > imptr;
206 603 : if( doesImageExist(itsImageName+String(".psf.tt0")) )
207 : {
208 603 : buildImage( imptr, (itsImageName+String(".psf.tt0")) );
209 :
210 : //cout << "Opening PSF image to read csys" << endl;
211 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".psf.tt0")) );
212 : }
213 0 : else if( doesImageExist(itsImageName+String(".residual.tt0")) )
214 : {
215 0 : buildImage( imptr, (itsImageName+String(".residual.tt0")) );
216 : //cout << "Opening Residual image to read csys" << endl;
217 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".residual.tt0")) );
218 : }
219 0 : else if( doesImageExist(itsImageName+String(".model.tt0")) )
220 : {
221 0 : buildImage( imptr, (itsImageName+String(".model.tt0")) );
222 : //cout << "Opening Model image to read csys" << endl;
223 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".model.tt0")) );
224 : }
225 : else
226 : {
227 : // How can this be right ?
228 : //cout << "Opening Sumwt image to read csys" << endl;
229 0 : buildImage( imptr, (itsImageName+String(".gridwt")) );
230 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".gridwt")) );
231 : }
232 :
233 603 : itsObjectName=imptr->imageInfo().objectName();
234 603 : itsImageShape = imptr->shape();
235 603 : itsCoordSys = imptr->coordinates();
236 603 : itsMiscInfo=imptr->miscInfo();
237 :
238 : }
239 : else
240 : {
241 0 : throw( AipsError( "Multi-term PSF, Residual or Model Images do not exist. Please create one of them." ) );
242 : }
243 :
244 1451 : if( doesImageExist(itsImageName+String(".residual.tt0")) ||
245 848 : doesImageExist(itsImageName+String(".psf.tt0")) )
246 : {
247 603 : if( sumwtexists )
248 : {
249 574 : std::shared_ptr<ImageInterface<Float> > imptr;
250 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
251 574 : buildImage( imptr, (itsImageName+String(".sumwt.tt0")) );
252 574 : itsNFacets = imptr->shape()[0];
253 574 : itsFacetId = 0;
254 574 : itsUseWeight = getUseWeightImage( *imptr );
255 : /////redo this here as psf may have different coordinates
256 574 : itsCoordSys = imptr->coordinates();
257 574 : itsMiscInfo=imptr->miscInfo();
258 574 : if( itsUseWeight && ! doesImageExist(itsImageName+String(".weight.tt0")) )
259 : {
260 0 : throw(AipsError("Internal error : MultiTerm Sumwt has a useweightimage=true but the weight image does not exist."));
261 : }
262 : }
263 : else
264 : {
265 29 : if(!ignoresumwt)
266 0 : {throw( AipsError( "Multi-term SumWt does not exist. Please create PSFs or Residuals." ) );}
267 : else
268 : {
269 29 : os << "SumWt.ttx do not exist. Proceeding only with PSFs" << LogIO::POST;
270 0 : std::shared_ptr<ImageInterface<Float> > imptr;
271 : // imptr.reset( new PagedImage<Float> (itsImageName+String(".sumwt.tt0")) );
272 29 : if( doesImageExist(itsImageName+String(".residual.tt0")) )
273 29 : {buildImage( imptr, (itsImageName+String(".residual.tt0")) );}
274 : else
275 0 : {buildImage( imptr, (itsImageName+String(".psf.tt0")) );}
276 :
277 29 : itsNFacets = 1;
278 29 : itsFacetId = 0;
279 29 : itsUseWeight = False;
280 29 : itsCoordSys = imptr->coordinates();
281 29 : itsMiscInfo=imptr->miscInfo();
282 : }
283 : }
284 : }// if psf0 or res0 exist
285 :
286 603 : if( ignorefacets==true ) itsNFacets=1;
287 :
288 603 : init();
289 603 : validate();
290 :
291 603 : }
292 :
293 : /*
294 : /////////////Constructor with pointers already created else where but taken over here
295 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(Block<std::shared_ptr<ImageInterface<Float> > > modelims,
296 : Block<std::shared_ptr<ImageInterface<Float> > >residims,
297 : Block<std::shared_ptr<ImageInterface<Float> > >psfims,
298 : Block<std::shared_ptr<ImageInterface<Float> > >weightims,
299 : Block<std::shared_ptr<ImageInterface<Float> > >restoredims,
300 : Block<std::shared_ptr<ImageInterface<Float> > >sumwtims,
301 : std::shared_ptr<ImageInterface<Float> > newmask,
302 : std::shared_ptr<ImageInterface<Float> > newalpha,
303 : std::shared_ptr<ImageInterface<Float> > newbeta)
304 : {
305 :
306 : itsPsfs=psfims;
307 : itsModels=modelims;
308 : itsResiduals=residims;
309 : itsWeights=weightims;
310 : itsImages=restoredims;
311 : itsSumWts=sumwtims;
312 : itsMask = newmask;
313 : itsAlpha = newalpha;
314 : itsBeta = newbeta;
315 :
316 : itsNTerms = itsResiduals.nelements();
317 : itsMiscInfo=Record();
318 :
319 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
320 : AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
321 : AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
322 :
323 : itsForwardGrids.resize( itsNTerms );
324 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
325 :
326 : itsImageShape=psfims[0]->shape();
327 : itsCoordSys = psfims[0]->coordinates();
328 : itsMiscInfo = psfims[0]->miscInfo();
329 :
330 : itsNFacets=sumwtims[0]->shape()[0];
331 : itsUseWeight=getUseWeightImage( *(sumwtims[0]) );
332 :
333 : itsImageName = String("reference"); // This is what the access functions use to guard against allocs...
334 :
335 : init();
336 : validate();
337 :
338 : }
339 : */
340 :
341 : //////////////////////////////////////////////////////////////////////////////////////////////////////
342 : /////////////Constructor with pointers already created else where but taken over here
343 : // used from getSubImageStore(), for example when creating the facets list
344 : // this would be safer if it was refactored as a copy constructor of the generic stuff +
345 : // initialization of the facet related parameters
346 286 : SIImageStoreMultiTerm::SIImageStoreMultiTerm(const Block<std::shared_ptr<ImageInterface<Float> > > &modelims,
347 : const Block<std::shared_ptr<ImageInterface<Float> > > &residims,
348 : const Block<std::shared_ptr<ImageInterface<Float> > > &psfims,
349 : const Block<std::shared_ptr<ImageInterface<Float> > > &weightims,
350 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredims,
351 : const Block<std::shared_ptr<ImageInterface<Float> > > &sumwtims,
352 : const Block<std::shared_ptr<ImageInterface<Float> > > &pbims,
353 : const Block<std::shared_ptr<ImageInterface<Float> > > &restoredpbcorims,
354 : const std::shared_ptr<ImageInterface<Float> > &newmask,
355 : const std::shared_ptr<ImageInterface<Float> > &newalpha,
356 : const std::shared_ptr<ImageInterface<Float> > &newbeta,
357 : const std::shared_ptr<ImageInterface<Float> > &newalphaerror,
358 : const std::shared_ptr<ImageInterface<Float> > &newalphapbcor,
359 : const std::shared_ptr<ImageInterface<Float> > &newbetapbcor,
360 : const CoordinateSystem& csys,
361 : const IPosition &imshape,
362 : const String &imagename,
363 : const String &objectname,
364 : const Record &miscinfo,
365 : const Int facet, const Int nfacets,
366 : const Int chan, const Int nchanchunks,
367 286 : const Int pol, const Int npolchunks)
368 : {
369 286 : itsPsfs=psfims;
370 286 : itsModels=modelims;
371 286 : itsResiduals=residims;
372 286 : itsWeights=weightims;
373 286 : itsImages=restoredims;
374 286 : itsSumWts=sumwtims;
375 286 : itsMask = newmask;
376 286 : itsAlpha = newalpha;
377 286 : itsBeta = newbeta;
378 286 : itsAlphaError = newalphaerror;
379 286 : itsAlphaPBcor = newalphapbcor;
380 286 : itsBetaPBcor = newbetapbcor;
381 :
382 286 : itsPBs=pbims;
383 286 : itsImagePBcors=restoredpbcorims;
384 :
385 286 : itsNTerms = itsResiduals.nelements();
386 286 : itsMiscInfo=Record();
387 :
388 286 : AlwaysAssert( itsPsfs.nelements() == 2*itsNTerms-1 , AipsError );
389 : // AlwaysAssert( itsPsfs.nelements()>0 && itsPsfs[0] , AipsError );
390 : // AlwaysAssert( itsSumWts.nelements()>0 && itsSumWts[0] , AipsError );
391 :
392 286 : itsForwardGrids.resize( itsNTerms );
393 286 : itsBackwardGrids.resize( 2 * itsNTerms - 1 );
394 :
395 286 : itsObjectName = objectname;
396 286 : itsMiscInfo = miscinfo;
397 :
398 286 : itsNFacets = nfacets;
399 286 : itsFacetId = facet;
400 286 : itsNChanChunks = nchanchunks;
401 286 : itsChanId = chan;
402 286 : itsNPolChunks = npolchunks;
403 286 : itsPolId = pol;
404 :
405 286 : itsParentImageShape = imshape;
406 286 : itsImageShape = imshape;
407 : ///// itsImageShape = IPosition(4,0,0,0,0);
408 :
409 286 : itsCoordSys = csys; // Hopefully this doesn't change for a reference image
410 286 : itsImageName = imagename;
411 :
412 :
413 : //-----------------------
414 286 : init(); // Connect parent pointers to the images.
415 : //-----------------------
416 :
417 : // Set these to null, to be set later upon first access.
418 : // Setting to null will hopefully set all elements of each array, to NULL.
419 286 : itsPsfs=std::shared_ptr<ImageInterface<Float> >();
420 286 : itsModels=std::shared_ptr<ImageInterface<Float> >();
421 286 : itsResiduals=std::shared_ptr<ImageInterface<Float> >();
422 286 : itsWeights=std::shared_ptr<ImageInterface<Float> >();
423 286 : itsImages=std::shared_ptr<ImageInterface<Float> >();
424 286 : itsSumWts=std::shared_ptr<ImageInterface<Float> >();
425 286 : itsPBs=std::shared_ptr<ImageInterface<Float> >();
426 :
427 286 : itsMask.reset( );
428 :
429 286 : validate();
430 :
431 286 : }
432 :
433 : //////////////////////////////////////////////////////////////////////////////////////////////////////
434 :
435 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
436 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
437 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
438 :
439 203 : uInt SIImageStoreMultiTerm::getNTaylorTerms(Bool dopsf)
440 : {
441 203 : return dopsf ? (2*itsNTerms-1) : itsNTerms;
442 : }
443 :
444 : // Check if images that are asked-for are ready and all have the same shape.
445 : /*
446 : Bool SIImageStoreMultiTerm::checkValidity(const Bool ipsf, const Bool iresidual,
447 : const Bool iweight, const Bool imodel, const Bool irestored,
448 : const Bool imask,const Bool isumwt,
449 : const Bool ialpha, const Bool ibeta)
450 : {
451 :
452 : // cout << "In MT::checkValidity imask is " << imask << endl;
453 :
454 : Bool valid = true;
455 :
456 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
457 : {
458 :
459 : if(ipsf==true)
460 : { psf(tix);
461 : valid = valid & ( itsPsfs[tix] && itsPsfs[tix]->shape()==itsImageShape ); }
462 : if(iweight==true)
463 : { weight(tix);
464 : valid = valid & ( itsWeights[tix] && itsWeights[tix]->shape()==itsImageShape ); }
465 :
466 : if(isumwt==true) {
467 : IPosition useShape(itsImageShape);
468 : useShape[0]=itsNFacets; useShape[1]=itsNFacets;
469 : sumwt(tix);
470 : valid = valid & ( itsSumWts[tix] && itsSumWts[tix]->shape()==useShape );
471 : }
472 :
473 : if( tix< itsNTerms )
474 : {
475 : if(iresidual==true)
476 : { residual(tix);
477 : valid = valid & ( itsResiduals[tix] && itsResiduals[tix]->shape()==itsImageShape ); }
478 : if(imodel==true)
479 : { model(tix);
480 : valid = valid & ( itsModels[tix] && itsModels[tix]->shape()==itsImageShape); }
481 : if(irestored==true)
482 : { image(tix);
483 : valid = valid & ( itsImages[tix] && itsImages[tix]->shape()==itsImageShape); }
484 : }
485 : }
486 :
487 : if(imask==true)
488 : { mask(); valid = valid & ( itsMask && itsMask->shape()==itsImageShape);
489 : // cout << " Mask null ? " << (bool) itsMask << endl;
490 : }
491 : if(ialpha==true)
492 : { alpha(); valid = valid & ( itsAlpha && itsAlpha->shape()==itsImageShape ); }
493 : if(ibeta==true)
494 : { beta(); valid = valid & ( itsBeta && itsBeta->shape()==itsImageShape ); }
495 :
496 : return valid;
497 :
498 : }
499 : */
500 :
501 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
502 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
503 : //////////////////////////////////////////////////////////////////////////////////////////////////////
504 : //////////////////////////////////////////////////////////////////////////////////////////////////////
505 :
506 1608 : SIImageStoreMultiTerm::~SIImageStoreMultiTerm()
507 : {
508 1608 : }
509 :
510 : //////////////////////////////////////////////////////////////////////////////////////////////////////
511 : //////////////////////////////////////////////////////////////////////////////////////////////////////
512 :
513 2110 : Bool SIImageStoreMultiTerm::releaseLocks()
514 : {
515 4220 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseLocks",WHERE) );
516 :
517 8326 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
518 : {
519 6216 : if( itsPsfs[tix] ) releaseImage( itsPsfs[tix] );
520 6216 : if( itsWeights[tix] ) releaseImage( itsWeights[tix] );
521 6216 : if( itsSumWts[tix] ) releaseImage( itsSumWts[tix] );
522 6216 : if( tix < itsNTerms )
523 : {
524 4163 : if( itsModels[tix] ) releaseImage( itsModels[tix] );
525 4163 : if( itsResiduals[tix] ) releaseImage( itsResiduals[tix] );
526 4163 : if( itsImages[tix] ) releaseImage( itsImages[tix] );
527 4163 : if( itsPBs[tix] ) releaseImage( itsPBs[tix] );
528 4163 : if( itsImagePBcors[tix] ) releaseImage( itsImagePBcors[tix] );
529 : }
530 : }
531 2110 : if( itsMask ) releaseImage( itsMask );
532 2110 : if( itsAlpha ) releaseImage( itsAlpha );
533 2110 : if( itsBeta ) releaseImage( itsBeta );
534 2110 : if( itsAlphaError ) releaseImage( itsAlphaError );
535 2110 : if( itsAlphaPBcor ) releaseImage( itsAlphaPBcor );
536 2110 : if( itsBetaPBcor ) releaseImage( itsBetaPBcor );
537 2110 : if( itsGridWt ) releaseImage( itsGridWt );
538 :
539 4220 : return true; // do something more intelligent here.
540 : }
541 :
542 0 : Bool SIImageStoreMultiTerm::releaseComplexGrids()
543 : {
544 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","releaseComplexGrids",WHERE) );
545 :
546 0 : for(uInt tix=0; tix<2*itsNTerms-1; tix++)
547 : {
548 0 : if( itsBackwardGrids[tix] ) releaseImage( itsBackwardGrids[tix] );
549 0 : if( tix < itsNTerms )
550 : {
551 0 : if( itsForwardGrids[tix] ) releaseImage( itsForwardGrids[tix] );
552 : }
553 : }
554 0 : return True; // do something more intelligent here.
555 : }
556 :
557 :
558 594 : Double SIImageStoreMultiTerm::getReferenceFrequency()
559 : {
560 : Double theRefFreq;
561 :
562 594 : Vector<Double> refpix = (itsCoordSys.spectralCoordinate()).referencePixel();
563 594 : AlwaysAssert( refpix.nelements()>0, AipsError );
564 594 : (itsCoordSys.spectralCoordinate()).toWorld( theRefFreq, refpix[0] );
565 : // cout << "Reading ref freq as : " << theRefFreq << endl;
566 1188 : return theRefFreq;
567 : }
568 :
569 : //////////////////////////////////////////////////////////////////////////////////////////////////////
570 : //////////////////////////////////////////////////////////////////////////////////////////////////////
571 0 : Vector<String> SIImageStoreMultiTerm::getModelImageName()
572 : {
573 0 : Vector<String> mods(itsNTerms);
574 0 : for(uInt tix=0;tix<itsNTerms;tix++)
575 0 : {mods[tix]=itsImageName + imageExts(MODEL)+".tt"+String::toString(tix);}
576 0 : return mods;
577 : }
578 :
579 12 : void SIImageStoreMultiTerm::setModelImage( const Vector<String> &modelnames )
580 : {
581 36 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
582 :
583 12 : if( modelnames.nelements() > itsNTerms )
584 0 : { throw(AipsError("We currently cannot support more than nterms images as the starting model. "));
585 : }
586 :
587 12 : if( modelnames.nelements() > 0 && modelnames.nelements() <= itsNTerms )
588 : {
589 30 : for(uInt tix=0;tix<modelnames.nelements();tix++)
590 : {
591 19 : setModelImageOne( modelnames[tix], tix );
592 : }
593 : }
594 :
595 11 : }
596 :
597 : /*
598 : void SIImageStoreMultiTerm::setModelImage( String modelname )
599 : {
600 : LogIO os( LogOrigin("SIImageStoreMultiTerm","setModelImage",WHERE) );
601 :
602 : for(uInt tix=0;tix<itsNTerms;tix++)
603 : {
604 :
605 : Directory immodel( modelname+String(".model.tt")+String::toString(tix) );
606 : if( !immodel.exists() )
607 : {
608 : os << "Starting model image does not exist for term : " << tix << LogIO::POST;
609 : }
610 : else
611 : {
612 : std::shared_ptr<PagedImage<Float> > newmodel( new PagedImage<Float>( modelname+String(".model.tt")+String::toString(tix) ) );
613 : // Check shapes, coordsys with those of other images. If different, try to re-grid here.
614 :
615 : if( newmodel->shape() != model(tix)->shape() )
616 : {
617 : os << "Regridding input model to target coordinate system for term " << tix << LogIO::POST;
618 : regridToModelImage( *newmodel , tix);
619 : // For now, throw an exception.
620 : //throw( AipsError( "Input model image "+modelname+".model.tt"+String::toString(tix)+" is not the same shape as that defined for output in "+ itsImageName + ".model" ) );
621 : }
622 : else
623 : {
624 : os << "Setting " << modelname << " as model for term " << tix << LogIO::POST;
625 : // Then, add its contents to itsModel.
626 : //itsModel->put( itsModel->get() + model->get() );
627 : /////( model(tix) )->put( newmodel->get() );
628 : model(tix)->copyData( LatticeExpr<Float> (*newmodel) );
629 : }
630 : }
631 : }//nterms
632 : }
633 : */
634 :
635 : //////////////////////////////////////////////////////////////////////////////////////////////////////
636 : //////////////////////////////////////////////////////////////////////////////////////////////////////
637 :
638 3404 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::psf(uInt term)
639 : {
640 3404 : AlwaysAssert( itsPsfs.nelements() > term, AipsError );
641 3408 : accessImage( itsPsfs[term], itsParentPsfs[term], imageExts(PSF)+".tt"+String::toString(term) );
642 3402 : return itsPsfs[term];
643 : }
644 7533 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::residual(uInt term)
645 : {
646 7537 : accessImage( itsResiduals[term], itsParentResiduals[term], imageExts(RESIDUAL)+".tt"+String::toString(term) );
647 : // Record mi = itsResiduals[term]->miscInfo(); ostringstream oss;mi.print(oss);cout<<"MiscInfo(res) " << term << " : " << oss.str() << endl;
648 :
649 7531 : return itsResiduals[term];
650 : }
651 1167 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::weight(uInt term)
652 : {
653 1167 : accessImage( itsWeights[term], itsParentWeights[term], imageExts(WEIGHT)+".tt"+String::toString(term) );
654 1167 : return itsWeights[term];
655 : }
656 5616 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::sumwt(uInt term)
657 : {
658 5616 : accessImage( itsSumWts[term], itsParentSumWts[term], imageExts(SUMWT)+".tt"+String::toString(term) );
659 :
660 :
661 5616 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 )
662 276 : {itsUseWeight = getUseWeightImage( *itsParentSumWts[0] );}
663 5616 : setUseWeightImage( *(itsSumWts[term]) , itsUseWeight); // Sets a flag in the SumWt image.
664 :
665 5616 : return itsSumWts[term];
666 : }
667 3671 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::model(uInt term)
668 : {
669 :
670 3677 : accessImage( itsModels[term], itsParentModels[term], imageExts(MODEL)+".tt"+String::toString(term) );
671 :
672 3669 : itsModels[term]->setUnits("Jy/pixel");
673 3669 : return itsModels[term];
674 : }
675 :
676 3488 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::image(uInt term)
677 : {
678 :
679 3488 : accessImage( itsImages[term], itsParentImages[term], imageExts(IMAGE)+".tt"+String::toString(term) );
680 3488 : itsImages[term]->setUnits("Jy/beam");
681 3488 : return itsImages[term];
682 : }
683 :
684 1325 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::pb(uInt term)
685 : {
686 :
687 1329 : accessImage( itsPBs[term], itsParentPBs[term], imageExts(PB)+".tt"+String::toString(term) );
688 1323 : return itsPBs[term];
689 : }
690 :
691 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::imagepbcor(uInt term)
692 : {
693 :
694 0 : accessImage( itsImagePBcors[term], itsParentImagePBcors[term], imageExts(IMAGE)+".tt"+String::toString(term)+ ".pbcor" );
695 0 : itsImagePBcors[term]->setUnits("Jy/beam");
696 0 : return itsImagePBcors[term];
697 : }
698 :
699 616 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::forwardGrid(uInt term){
700 616 : if( itsForwardGrids[term] )// && (itsForwardGrids[term]->shape() == itsImageShape))
701 506 : return itsForwardGrids[term];
702 220 : Vector<Int> whichStokes(0);
703 220 : IPosition cimageShape;
704 110 : cimageShape=itsImageShape;
705 110 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
706 220 : whichStokes, itsDataPolRep);
707 110 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
708 110 : cimageShape(2)=whichStokes.nelements();
709 110 : itsForwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
710 110 : return itsForwardGrids[term];
711 : }
712 1346 : std::shared_ptr<ImageInterface<Complex> > SIImageStoreMultiTerm::backwardGrid(uInt term){
713 :
714 1346 : if( itsBackwardGrids[term] && (itsBackwardGrids[term]->shape() == itsImageShape))
715 1001 : return itsBackwardGrids[term];
716 : // cout << "MT : Making backward grid of shape : " << itsImageShape << endl;
717 690 : Vector<Int> whichStokes(0);
718 690 : IPosition cimageShape;
719 345 : cimageShape=itsImageShape;
720 345 : CoordinateSystem cimageCoord = StokesImageUtil::CStokesCoord( itsCoordSys,
721 690 : whichStokes, itsDataPolRep);
722 345 : cimageCoord.setObsInfo(itsCoordSys.obsInfo());
723 345 : cimageShape(2)=whichStokes.nelements();
724 345 : itsBackwardGrids[term].reset(new TempImage<Complex>(TiledShape(cimageShape, tileShape()), cimageCoord, memoryBeforeLattice()));
725 345 : return itsBackwardGrids[term];
726 : }
727 :
728 452 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alpha()
729 : {
730 452 : if( itsAlpha && itsAlpha->shape() == itsImageShape ) { return itsAlpha; }
731 : // checkRef( itsAlpha , "alpha" );
732 113 : itsAlpha = openImage( itsImageName+String(".alpha"), false );
733 : // itsAlpha->setUnits("Alpha");
734 113 : return itsAlpha;
735 : }
736 :
737 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::beta()
738 : {
739 0 : if( itsBeta && itsBeta->shape() == itsImageShape ) { return itsBeta; }
740 : // checkRef( itsBeta , "beta" );
741 0 : itsBeta = openImage( itsImageName+String(".beta"), false );
742 : // itsBeta->setUnits("Beta");
743 0 : return itsBeta;
744 : }
745 :
746 678 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphaerror()
747 : {
748 678 : if( itsAlphaError && itsAlphaError->shape() == itsImageShape ) { return itsAlphaError; }
749 : // checkRef( itsAlpha , "alpha" );
750 113 : itsAlphaError = openImage( itsImageName+String(".alpha.error"), false );
751 : // itsAlpha->setUnits("Alpha");
752 113 : return itsAlphaError;
753 : }
754 :
755 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::alphapbcor()
756 : {
757 0 : if( itsAlphaPBcor && itsAlphaPBcor->shape() == itsImageShape ) { return itsAlphaPBcor; }
758 : // checkRef( itsAlphaPBcor , "alpha" );
759 0 : itsAlphaPBcor = openImage( itsImageName+String(".alpha.pbcor"), False );
760 : // itsAlphaPBcor->setUnits("Alpha");
761 0 : return itsAlphaPBcor;
762 : }
763 :
764 0 : std::shared_ptr<ImageInterface<Float> > SIImageStoreMultiTerm::betapbcor()
765 : {
766 0 : if( itsBetaPBcor && itsBetaPBcor->shape() == itsImageShape ) { return itsBetaPBcor; }
767 : // checkRef( itsBetaPBcor , "beta" );
768 0 : itsBetaPBcor = openImage( itsImageName+String(".beta.pbcor"), False );
769 : // itsBetaPBcor->setUnits("Beta");
770 0 : return itsBetaPBcor;
771 : }
772 :
773 :
774 :
775 : // TODO : Move to an image-wrapper class ? Same function exists in SynthesisDeconvolver.
776 17133 : Bool SIImageStoreMultiTerm::doesImageExist(String imagename)
777 : {
778 51399 : LogIO os( LogOrigin("SIImageStoreMultiTerm","doesImageExist",WHERE) );
779 34266 : Directory image( imagename );
780 34266 : return image.exists();
781 : }
782 :
783 :
784 0 : void SIImageStoreMultiTerm::resetImages( Bool resetpsf, Bool resetresidual, Bool resetweight )
785 : {
786 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
787 : {
788 0 : if( resetpsf ) psf(tix)->set(0.0);
789 :
790 0 : if( tix < itsNTerms ) {
791 0 : if( resetresidual ) {
792 : //removeMask( residual(tix) );
793 0 : residual(tix)->set(0.0);
794 : }
795 : }
796 0 : if( resetweight && itsWeights[tix] ) weight(tix)->set(0.0);
797 0 : if( resetweight ) sumwt(tix)->set(0.0);
798 : }//nterms
799 0 : }
800 :
801 0 : void SIImageStoreMultiTerm::addImages( std::shared_ptr<SIImageStore> imagestoadd,
802 : Bool addpsf, Bool addresidual, Bool addweight, Bool adddensity)
803 : {
804 0 : for(uInt tix=0;tix<2*itsNTerms-1;tix++)
805 : {
806 :
807 0 : if(addpsf)
808 : {
809 0 : LatticeExpr<Float> adderPsf( *(psf(tix)) + *(imagestoadd->psf(tix)) );
810 0 : psf(tix)->copyData(adderPsf);
811 : }
812 0 : if(addweight)
813 : {
814 :
815 0 : if(getUseWeightImage( *(imagestoadd->sumwt(tix)) ) ) // Access and add weight only if it is needed.
816 : {
817 0 : LatticeExpr<Float> adderWeight( *(weight(tix)) + *(imagestoadd->weight(tix)) );
818 0 : weight(tix)->copyData(adderWeight);
819 : }
820 :
821 0 : LatticeExpr<Float> adderSumWt( *(sumwt(tix)) + *(imagestoadd->sumwt(tix)) );
822 0 : sumwt(tix)->copyData(adderSumWt);
823 0 : setUseWeightImage( *sumwt(tix), getUseWeightImage(*(imagestoadd->sumwt(tix)) ) );
824 : }
825 :
826 0 : if(tix < itsNTerms && addresidual)
827 : {
828 0 : LatticeExpr<Float> adderRes( *(residual(tix)) + *(imagestoadd->residual(tix)) );
829 0 : residual(tix)->copyData(adderRes);
830 : }
831 :
832 0 : if( tix==0 && adddensity )
833 : {
834 0 : LatticeExpr<Float> adderDensity( *(gridwt()) + *(imagestoadd->gridwt()) );
835 0 : gridwt()->copyData(adderDensity);
836 : }
837 :
838 : }
839 0 : }
840 :
841 105 : void SIImageStoreMultiTerm::dividePSFByWeight(const Float /*pblimit*/)
842 : {
843 315 : LogIO os( LogOrigin("SIImageStoreMultiTerm","dividePSFByWeight",WHERE) );
844 :
845 : //// for(uInt tix=0;tix<2*itsNTerms-1;tix++)
846 414 : for(Int tix=2*itsNTerms-1-1;tix>-1;tix--) // AAH go backwards so that zeroth term is normalized last..... sigh sigh sigh.
847 : {
848 :
849 : //cout << "npsfs : " << itsPsfs.nelements() << " and tix : " << tix << endl;
850 :
851 309 : normPSF(tix);
852 :
853 309 : if ( itsUseWeight) {
854 118 : divideImageByWeightVal( *weight(tix) );
855 : }
856 :
857 : }
858 105 : }
859 :
860 105 : void SIImageStoreMultiTerm::normalizePrimaryBeam(const Float pblimit)
861 : {
862 315 : LogIO os( LogOrigin("SIImageStoreMultiTerm","normalizePrimaryBeam",WHERE) );
863 105 : if ( itsUseWeight) {
864 :
865 40 : makePBFromWeight(pblimit);
866 : }
867 65 : else { makePBImage(pblimit); }
868 : // calcSensitivity();
869 105 : }
870 :
871 :
872 : // Make another for the PSF too.
873 180 : void SIImageStoreMultiTerm::divideResidualByWeight(Float pblimit, const String normtype)
874 : {
875 540 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideResidualByWeight",WHERE) );
876 :
877 180 : if( itsUseWeight )
878 : {
879 64 : itsPBScaleFactor = getPbMax();
880 : }
881 :
882 535 : for(uInt tix=0;tix<itsNTerms;tix++)
883 : {
884 :
885 355 : divideImageByWeightVal( *residual(tix) );
886 :
887 : // if(doesImageExist(itsImageName+String(".weight.tt0")) )
888 355 : if( itsUseWeight )
889 : {
890 127 : LatticeExpr<Float> ratio;
891 127 : Float scalepb = fabs(pblimit) * itsPBScaleFactor * itsPBScaleFactor ;
892 127 : if( normtype=="flatnoise"){
893 357 : LatticeExpr<Float> deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) * itsPBScaleFactor );
894 119 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
895 119 : os << " by [ sqrt(weightimage) * " << itsPBScaleFactor ;
896 119 : os << " ] to get flat noise with unit pb peak."<< LogIO::POST;
897 357 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
898 238 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
899 119 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
900 : }
901 8 : else if(normtype=="pbsquare"){
902 8 : Float deno = itsPBScaleFactor*itsPBScaleFactor ;
903 8 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
904 8 : os << deno ;
905 8 : os << " ] to get optimal sig/noise with unit pb peak."<< LogIO::POST;
906 :
907 8 : ratio=LatticeExpr<Float> ( ( *(residual(tix)) ) / ( deno ) );
908 :
909 :
910 :
911 : }
912 0 : else if( normtype=="flatsky") {
913 0 : LatticeExpr<Float> deno( *(weight(0)) );
914 0 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".residual.tt")+String::toString(tix) ;
915 0 : os << " by [ weight ] to get flat sky"<< LogIO::POST;
916 0 : LatticeExpr<Float> mask( iif( (deno) > scalepb , 1.0, 0.0 ) );
917 0 : LatticeExpr<Float> maskinv( iif( (deno) > scalepb , 0.0, 1.0 ) );
918 0 : ratio=LatticeExpr<Float> ( ( (*(residual(tix))) * mask ) / ( deno + maskinv ) );
919 : }
920 : else{
921 0 : throw(AipsError("Don't know how to proceed with normtype "+normtype));
922 : }
923 :
924 :
925 :
926 127 : residual(tix)->copyData(ratio);
927 : }
928 :
929 355 : if( (residual(tix)->getDefaultMask()=="") && hasPB() && pblimit >= 0.0 )
930 177 : {copyMask(pb(),residual(tix));}
931 :
932 355 : if( pblimit <0.0 && (residual(tix)->getDefaultMask()).matches("mask0") ) removeMask( residual(tix) );
933 :
934 : }
935 180 : }
936 :
937 178 : void SIImageStoreMultiTerm::divideModelByWeight(Float pblimit, const String normtype)
938 : {
939 356 : LogIO os( LogOrigin("SIImageStoreMultiTerm","divideModelByWeight",WHERE) );
940 :
941 356 : if( itsUseWeight // only when needed
942 178 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
943 : {
944 :
945 64 : if( normtype=="flatsky") {
946 0 : Array<Float> arrmod;
947 0 : model(0)->get( arrmod, true );
948 :
949 0 : os << "Model is already flat sky with peak flux : " << max(arrmod);
950 0 : os << ". No need to divide before prediction" << LogIO::POST;
951 :
952 0 : return;
953 : }
954 :
955 64 : itsPBScaleFactor = getPbMax();
956 :
957 191 : for(uInt tix=0;tix<itsNTerms;tix++)
958 254 : { LatticeExpr<Float> deno;
959 127 : if(normtype=="flatnoise"){
960 119 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
961 119 : os << " by [ sqrt(weight) / " << itsPBScaleFactor ;
962 119 : os <<" ] to get to flat sky model before prediction" << LogIO::POST;
963 :
964 119 : deno = LatticeExpr<Float> ( sqrt( abs(*(weight(0))) ) / itsPBScaleFactor );
965 : }
966 8 : else if(normtype=="pbsquare"){
967 8 : os << LogIO::NORMAL1 << "Dividing " << itsImageName+String(".model")+String::toString(tix);
968 8 : os << " by [ (weight) / " << itsPBScaleFactor*itsPBScaleFactor ;
969 8 : os <<" ] to get an optimal sig/noise model before prediction" << LogIO::POST;
970 :
971 8 : deno = LatticeExpr<Float> ( abs(*(weight(0))) / (itsPBScaleFactor*itsPBScaleFactor) );
972 : }
973 381 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
974 381 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
975 381 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) / ( deno + maskinv ) );
976 :
977 127 : itsModels[tix]->copyData(ratio);
978 : }
979 : }
980 : }
981 :
982 :
983 178 : void SIImageStoreMultiTerm::multiplyModelByWeight(Float pblimit, const String normtype)
984 : {
985 356 : LogIO os( LogOrigin("SIImageStoreMultiTerm","multiplyModelByWeight",WHERE) );
986 :
987 :
988 356 : if( itsUseWeight // only when needed
989 178 : && weight() )// i.e. only when possible. For an initial starting model, don't need wt anyway.
990 : {
991 :
992 64 : if( normtype=="flatsky") {
993 0 : os << "Model is already flat sky. No need to multiply back after prediction" << LogIO::POST;
994 0 : return;
995 : }
996 :
997 64 : itsPBScaleFactor = getPbMax();
998 :
999 191 : for(uInt tix=0;tix<itsNTerms;tix++)
1000 : {
1001 254 : LatticeExpr<Float> deno;
1002 127 : if( normtype=="flatnoise") {
1003 119 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1004 119 : os << " by [ sqrt(weight) / " << itsPBScaleFactor;
1005 119 : os << " ] to take model back to flat noise with unit pb peak." << LogIO::POST;
1006 :
1007 119 : deno=LatticeExpr<Float> ( sqrt( abs(*(weight(0)) ) ) / itsPBScaleFactor );
1008 : }
1009 8 : else if ( normtype=="pbsquare"){
1010 8 : os << LogIO::NORMAL1 << "Multiplying " << itsImageName+String(".model")+String::toString(tix);
1011 8 : os << " by [ weight / " << itsPBScaleFactor*itsPBScaleFactor;
1012 8 : os << " ] to take model back to optima sig/noise with unit pb peak." << LogIO::POST;
1013 :
1014 8 : deno=LatticeExpr<Float> ( abs(*(weight(0)) ) / (itsPBScaleFactor*itsPBScaleFactor) );
1015 : }
1016 : else{
1017 0 : throw(AipsError("No idea of what to do for "+normtype));
1018 : }
1019 :
1020 381 : LatticeExpr<Float> mask( iif( (deno) > fabs(pblimit) , 1.0, 0.0 ) );
1021 381 : LatticeExpr<Float> maskinv( iif( (deno) > fabs(pblimit) , 0.0, 1.0 ) );
1022 381 : LatticeExpr<Float> ratio( ( (*(model(tix))) * mask ) * ( deno + maskinv ) );
1023 :
1024 127 : itsModels[tix]->copyData(ratio);
1025 : }
1026 : }
1027 : }
1028 :
1029 :
1030 125 : void SIImageStoreMultiTerm::restore(GaussianBeam& rbeam, String& usebeam, uInt /*term*/, Float psfcutoff)
1031 : {
1032 :
1033 250 : LogIO os( LogOrigin("SIImageStoreMultiTerm","restore",WHERE) );
1034 :
1035 372 : for(uInt tix=0; tix<itsNTerms; tix++)
1036 : {
1037 247 : SIImageStore::restore(rbeam, usebeam, tix, psfcutoff);
1038 : }
1039 :
1040 125 : calculateAlphaBeta("image");
1041 :
1042 125 : }// restore
1043 :
1044 125 : void SIImageStoreMultiTerm::calculateAlphaBeta(String imtype)
1045 : {
1046 375 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calculateAlphaBeta",WHERE) );
1047 : try
1048 : {
1049 :
1050 : // Check if this is Stokes I.
1051 125 : Bool isOnlyStokesI=True;
1052 :
1053 250 : Vector<Int> stokes ( (itsParentCoordSys.stokesCoordinate()).stokes() );
1054 125 : AlwaysAssert( stokes.nelements()>0 , AipsError);
1055 125 : if( stokes.nelements()>1 || stokes[0]!=1 ) isOnlyStokesI=False;
1056 :
1057 125 : if( ! isOnlyStokesI )
1058 9 : { os << "Alpha and Beta images will not be computed for images that contain anything other than stokes I in them." << LogIO::POST; }
1059 :
1060 : //Calculate alpha, beta
1061 125 : if( itsNTerms > 1 && isOnlyStokesI)
1062 : {
1063 : // Calculate alpha and beta
1064 339 : LatticeExprNode leMaxRes = max( *( residual(0) ) );
1065 113 : Float maxres = leMaxRes.getFloat();
1066 113 : Float specthreshold = maxres/10.0; //////////// do something better here.....
1067 :
1068 113 : os << "Calculating spectral parameters for Intensity > peakresidual/10 = " << specthreshold << " Jy/beam" << LogIO::POST;
1069 :
1070 : /////////////////////////////////////////////////////////
1071 : /////// Calculate alpha
1072 113 : if(imtype=="image")
1073 : {
1074 339 : LatticeExpr<Float> mask1(iif(((*(image(0))))>(specthreshold),1.0,0.0));
1075 339 : LatticeExpr<Float> mask0(iif(((*(image(0))))>(specthreshold),0.0,1.0));
1076 339 : LatticeExpr<Float> alphacalc( (((*(image(1))))*mask1)/(((*(image(0))))+(mask0)) );
1077 113 : alpha()->copyData(alphacalc);
1078 :
1079 226 : ImageInfo ii = image(0)->imageInfo();
1080 : // Set the restoring beam for alpha
1081 113 : alpha()->setImageInfo(ii);
1082 : //alpha()->table().unmarkForDelete();
1083 :
1084 : // Make a mask for the alpha image
1085 339 : LatticeExpr<Bool> lemask(iif(((*(image(0))) > specthreshold) , True, False));
1086 113 : removeMask( alpha() );
1087 113 : createMask( lemask, (alpha()) );
1088 :
1089 : /////// Calculate alpha error
1090 113 : alphaerror()->set(0.0);
1091 :
1092 339 : LatticeExpr<Float> alphacalcerror( abs(alphacalc) * sqrt( ( (*residual(0)*mask1)/(*image(0)+mask0) )*( (*residual(0)*mask1)/(*image(0)+mask0) ) + ( (*residual(1)*mask1)/(*image(1)+mask0) )*( (*residual(1)*mask1)/(*image(1)+mask0) ) ) );
1093 113 : alphaerror()->copyData(alphacalcerror);
1094 113 : alphaerror()->setImageInfo(ii);
1095 113 : removeMask(alphaerror());
1096 113 : createMask(lemask, alphaerror());
1097 : // alphaerror()->table().unmarkForDelete();
1098 113 : os << "Written Spectral Index Error Image : " << alphaerror()->name() << LogIO::POST;
1099 :
1100 113 : if(itsNTerms>2) // calculate beta too.
1101 : {
1102 0 : beta()->set(0.0);
1103 0 : LatticeExpr<Float> betacalc( (*image(2)*mask1)/((*image(0))+(mask0))-0.5*(*alpha())*((*alpha())-1.0) );
1104 0 : beta()->copyData(betacalc);
1105 0 : beta()->setImageInfo(ii);
1106 : //imbeta.setUnits(Unit("Spectral Curvature"));
1107 0 : removeMask(beta());
1108 0 : createMask(lemask, beta());
1109 0 : os << "Written Spectral Curvature Image : " << beta()->name() << LogIO::POST;
1110 : }
1111 :
1112 : }
1113 : else
1114 : {
1115 0 : LatticeExpr<Float> mask1(iif(((*(imagepbcor(0))))>(specthreshold),1.0,0.0));
1116 0 : LatticeExpr<Float> mask0(iif(((*(imagepbcor(0))))>(specthreshold),0.0,1.0));
1117 0 : LatticeExpr<Float> alphacalc( (((*(imagepbcor(1))))*mask1)/(((*(imagepbcor(0))))+(mask0)) );
1118 0 : alphapbcor()->copyData(alphacalc);
1119 :
1120 0 : ImageInfo ii = image(0)->imageInfo();
1121 : // Set the restoring beam for alpha
1122 0 : alphapbcor()->setImageInfo(ii);
1123 : //alpha()->table().unmarkForDelete();
1124 :
1125 : // Make a mask for the alpha image
1126 0 : LatticeExpr<Bool> lemask(iif(((*(imagepbcor(0))) > specthreshold) , True, False));
1127 0 : removeMask( alphapbcor() );
1128 0 : createMask( lemask, (alphapbcor()) );
1129 :
1130 : /////////////////////////////////////////////////////////
1131 0 : if(itsNTerms>2) // calculate beta too.
1132 : {
1133 0 : betapbcor()->set(0.0);
1134 0 : LatticeExpr<Float> betacalc( (*imagepbcor(2)*mask1)/((*imagepbcor(0))+(mask0))-0.5*(*alphapbcor())*((*alphapbcor())-1.0) );
1135 0 : betapbcor()->copyData(betacalc);
1136 0 : betapbcor()->setImageInfo(ii);
1137 : //imbeta.setUnits(Unit("Spectral Curvature"));
1138 0 : removeMask(betapbcor());
1139 0 : createMask(lemask, betapbcor());
1140 0 : os << "Written Spectral Curvature Image : " << betapbcor()->name() << LogIO::POST;
1141 : }
1142 :
1143 : }// pbcor
1144 :
1145 : }//if nterms>1
1146 : }
1147 0 : catch(AipsError &x)
1148 : {
1149 0 : throw( AipsError("Error in computing Alpha (and Beta) : " + x.getMesg() ) );
1150 : }
1151 :
1152 125 : }// calculateAlphaBeta
1153 :
1154 1 : void SIImageStoreMultiTerm::pbcor()
1155 : {
1156 :
1157 2 : LogIO os( LogOrigin("SIImageStoreMultiTerm","pbcorPlane",WHERE) );
1158 :
1159 : /// Temp Code to prevent this approximate PBCOR from happening for EVLA data
1160 : if(1)
1161 : {
1162 1 : String telescope = itsCoordSys.obsInfo().telescope();
1163 1 : if ( telescope != "ALMA" )
1164 : {
1165 1 : os << LogIO::WARN << "Wideband (multi-term) PB correction is not yet available via tclean in the 4.7 release. Please use the widebandpbcor task instead. "<< LogIO::POST;
1166 1 : return;
1167 : }
1168 : else
1169 : {
1170 0 : os << LogIO::WARN << "Wideband (multi-term) PB Correction is currently only an approximation. It assumes no PB frequency dependence. This code has been added for the 4.7 release to support the current ALMA pipeline, which does not apply corrections for the frequency dependence of the primary beam across small fractional bandwidths. Please look at the help for the 'pbcor' parameter and use the widebandpbcor task if needed. " <<LogIO::POST;
1171 : }
1172 :
1173 :
1174 : }
1175 :
1176 :
1177 : // message saying that it's only stokes I for now...
1178 :
1179 0 : for(uInt tix=0; tix<itsNTerms; tix++)
1180 : {
1181 0 : SIImageStore::pbcor(tix);
1182 : }
1183 :
1184 0 : calculateAlphaBeta("pbcor");
1185 :
1186 : }
1187 :
1188 105 : void SIImageStoreMultiTerm::calcSensitivity()
1189 : {
1190 315 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calcSensitivity",WHERE) );
1191 :
1192 : // Construct Hessian.
1193 :
1194 210 : Matrix<Float> hess(IPosition(2,itsNTerms,itsNTerms));
1195 312 : for(uInt tay1=0;tay1<itsNTerms;tay1++)
1196 618 : for(uInt tay2=0;tay2<itsNTerms;tay2++)
1197 : {
1198 : //uInt taymin = (tay1<=tay2)? tay1 : tay2;
1199 : //uInt taymax = (tay1>=tay2)? tay1 : tay2;
1200 : //uInt ind = (taymax*(taymax+1)/2)+taymin;
1201 :
1202 411 : uInt ind = tay1+tay2;
1203 411 : AlwaysAssert( ind < 2*itsNTerms-1, AipsError );
1204 :
1205 411 : Array<Float> lsumwt;
1206 411 : sumwt( ind )->get( lsumwt, false );
1207 : // cout << "lsumwt shape : " << lsumwt.shape() << endl;
1208 411 : AlwaysAssert( lsumwt.shape().nelements()==4, AipsError );
1209 411 : AlwaysAssert( lsumwt.shape()[0]>0, AipsError );
1210 :
1211 : // hess(tay1,tay2) = lsumwt(IPosition(1,0)); //Always pick sumwt from first facet only.
1212 411 : hess(tay1,tay2) = lsumwt(IPosition(4,0,0,0,0)); //Always pick sumwt from first facet only.
1213 : }
1214 :
1215 105 : os << "Multi-Term Hessian Matrix : " << hess << LogIO::POST;
1216 :
1217 : // Invert Hessian.
1218 : try
1219 : {
1220 105 : Float deter=0.0;
1221 210 : Matrix<Float> invhess;
1222 105 : invertSymPosDef(invhess,deter,hess);
1223 105 : os << "Multi-Term Covariance Matrix : " << invhess << LogIO::POST;
1224 :
1225 : // Just print the sqrt of the diagonal elements.
1226 :
1227 312 : for(uInt tix=0;tix<itsNTerms;tix++)
1228 : {
1229 207 : os << "[" << itsImageName << "][Taylor"<< tix << "] Theoretical sensitivity (Jy/bm):" ;
1230 207 : if( invhess(tix,tix) > 0.0 ) { os << sqrt(invhess(tix,tix)) << LogIO::POST; }
1231 0 : else { os << "none" << LogIO::POST; }
1232 : }
1233 : }
1234 0 : catch(AipsError &x)
1235 : {
1236 0 : os << LogIO::WARN << "Cannot invert Hessian Matrix : " << x.getMesg() << " || Calculating approximate sensitivity " << LogIO::POST;
1237 :
1238 : // Approximate : 1/h^2
1239 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1240 : {
1241 0 : Array<Float> lsumwt;
1242 0 : AlwaysAssert( 2*tix < 2*itsNTerms-1, AipsError );
1243 0 : sumwt(2*tix)->get( lsumwt , false );
1244 :
1245 0 : IPosition shp( lsumwt.shape() );
1246 : //cout << "Sumwt shape : " << shp << " : " << lsumwt << endl;
1247 : //AlwaysAssert( shp.nelements()==4 , AipsError );
1248 :
1249 0 : os << "[" << itsImageName << "][Taylor"<< tix << "] Approx Theoretical sensitivity (Jy/bm):" ;
1250 :
1251 0 : IPosition it(4,0,0,0,0);
1252 0 : for ( it[0]=0; it[0]<shp[0]; it[0]++)
1253 0 : for ( it[1]=0; it[1]<shp[1]; it[1]++)
1254 0 : for ( it[2]=0; it[2]<shp[2]; it[2]++)
1255 0 : for ( it[3]=0; it[3]<shp[3]; it[3]++)
1256 : {
1257 0 : if( shp[0]>1 ){os << "f"<< it[0]+(it[1]*shp[0]) << ":" ;}
1258 0 : if( lsumwt( it ) > 1e-07 )
1259 : {
1260 0 : os << 1.0/(sqrt(lsumwt(it))) << " " ;
1261 : }
1262 : else
1263 : {
1264 0 : os << "none" << " ";
1265 : }
1266 : }
1267 :
1268 0 : os << LogIO::POST;
1269 : }
1270 : }
1271 :
1272 105 : calcFractionalBandwidth();
1273 105 : }
1274 :
1275 228 : Double SIImageStoreMultiTerm::calcFractionalBandwidth()
1276 : {
1277 :
1278 456 : LogIO os( LogOrigin("SIImageStoreMultiTerm","calcFractionalBandwidth",WHERE) );
1279 :
1280 228 : Double fbw=1.0;
1281 :
1282 912 : for(uInt i=0; i<itsCoordSys.nCoordinates(); i++)
1283 : {
1284 684 : if( itsCoordSys.type(i) == Coordinate::SPECTRAL )
1285 : {
1286 456 : SpectralCoordinate speccoord(itsCoordSys.spectralCoordinate(i));
1287 228 : Double startfreq=0.0,startpixel=-0.5;
1288 228 : Double endfreq=0.0,endpixel=+0.5;
1289 228 : speccoord.toWorld(startfreq,startpixel);
1290 228 : speccoord.toWorld(endfreq,endpixel);
1291 228 : Double midfreq = (endfreq+startfreq)/2.0;
1292 228 : fbw = ((endfreq - startfreq)/midfreq) * 100.0;
1293 228 : os << "MFS frequency range : " << startfreq/1.0e+9 << " GHz -> " << endfreq/1.0e+9 << "GHz.";
1294 228 : os << "Fractional Bandwidth : " << fbw << " %.";
1295 228 : os << "Reference Frequency for Taylor Expansion : "<< getReferenceFrequency()/1.0e+9 << "GHz." << LogIO::POST;
1296 : }
1297 : }
1298 456 : return fbw;
1299 : }
1300 :
1301 :
1302 : // Check for non-zero model (this is different from getting model flux, for derived SIIMMT)
1303 427 : Bool SIImageStoreMultiTerm::isModelEmpty()
1304 : {
1305 : /// There MUST be a more efficient way to do this !!!!! I hope.
1306 : /// Maybe save this info and change it anytime model() is accessed....
1307 427 : Bool emptymodel=true;
1308 1273 : for(uInt tix=0;tix<itsNTerms;tix++)
1309 : {
1310 : //if( fabs( getModelFlux(tix) ) > 1e-08 ) emptymodel=false;
1311 1540 : if( doesImageExist(itsImageName+String(".model.tt")+String::toString(tix)) &&
1312 1174 : fabs( getModelFlux(tix) ) > 1e-08 ) emptymodel=false;
1313 : }
1314 427 : return emptymodel;
1315 : }
1316 :
1317 0 : void SIImageStoreMultiTerm::printImageStats()
1318 : {
1319 0 : LogIO os( LogOrigin("SIImageStoreMultiTerm","printImageStats",WHERE) );
1320 : // FIXME minresmask needs to be initialized here, or else compiler complains
1321 0 : Float minresmask = 0, maxresmask, minres, maxres;
1322 0 : ArrayLattice<Bool> pixelmask(residual()->getMask());
1323 :
1324 : // findMinMax( residual()->get(), mask()->get(), minres, maxres, minresmask, maxresmask );
1325 :
1326 0 : if (hasMask())
1327 : {
1328 : // findMinMaxLattice(*residual(), *mask() , maxres,maxresmask, minres, minresmask);
1329 0 : findMinMaxLattice(*residual(), *mask() , pixelmask, maxres,maxresmask, minres, minresmask);
1330 : }
1331 : else
1332 : {
1333 0 : LatticeExpr<Float> reswithpixmask(iif(pixelmask, *residual(), 0));
1334 : //LatticeExprNode pres( max( *residual() ) );
1335 0 : LatticeExprNode pres( max( reswithpixmask ) );
1336 0 : maxres = pres.getFloat();
1337 : //LatticeExprNode pres2( min( *residual() ) );
1338 0 : LatticeExprNode pres2( min( reswithpixmask ) );
1339 0 : minres = pres2.getFloat();
1340 : }
1341 :
1342 0 : os << "[" << itsImageName << "]" ;
1343 0 : os << " Peak residual (max,min) " ;
1344 0 : if( minresmask!=0.0 || maxresmask!=0.0 )
1345 0 : { os << "within mask : (" << maxresmask << "," << minresmask << ") "; }
1346 0 : os << "over full image : (" << maxres << "," << minres << ")" << LogIO::POST;
1347 :
1348 0 : os << "[" << itsImageName << "] Total Model Flux : " ;
1349 0 : for(uInt tix=0;tix<itsNTerms;tix++)
1350 0 : {os << getModelFlux(tix) << "(tt" << String::toString(tix) << ")"; }
1351 0 : os<<LogIO::POST;
1352 :
1353 0 : }
1354 :
1355 286 : std::shared_ptr<SIImageStore> SIImageStoreMultiTerm::getSubImageStore(const Int facet, const Int nfacets,
1356 : const Int chan, const Int nchanchunks,
1357 : const Int pol, const Int npolchunks)
1358 : {
1359 : std::shared_ptr<SIImageStore> multiTermStore =
1360 286 : std::make_shared<SIImageStoreMultiTerm>(itsModels, itsResiduals, itsPsfs, itsWeights, itsImages, itsSumWts, itsPBs, itsImagePBcors, itsMask, itsAlpha, itsBeta, itsAlphaError,itsAlphaPBcor, itsBetaPBcor, itsCoordSys, itsParentImageShape, itsImageName, itsObjectName, itsMiscInfo, facet, nfacets, chan, nchanchunks, pol, npolchunks);
1361 286 : return multiTermStore;
1362 : }
1363 :
1364 :
1365 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1366 :
1367 : //
1368 : //-------------------------------------------------------------
1369 : // Initialize the internals of the object. Perhaps other such
1370 : // initializations in the constructors can be moved here too.
1371 : //
1372 1005 : void SIImageStoreMultiTerm::init()
1373 : {
1374 1005 : imageExts.resize(MAX_IMAGE_IDS);
1375 :
1376 1005 : imageExts(MASK)=".mask";
1377 1005 : imageExts(PSF)=".psf";
1378 1005 : imageExts(MODEL)=".model";
1379 1005 : imageExts(RESIDUAL)=".residual";
1380 1005 : imageExts(WEIGHT)=".weight";
1381 1005 : imageExts(IMAGE)=".image";
1382 1005 : imageExts(SUMWT)=".sumwt";
1383 1005 : imageExts(GRIDWT)=".gridwt";
1384 1005 : imageExts(PB)=".pb";
1385 1005 : imageExts(FORWARDGRID)=".forward";
1386 1005 : imageExts(BACKWARDGRID)=".backward";
1387 : // imageExts(IMAGEPBCOR)=".image.pbcor";
1388 :
1389 1005 : itsParentPsfs.resize(itsPsfs.nelements());
1390 1005 : itsParentModels.resize(itsModels.nelements());
1391 1005 : itsParentResiduals.resize(itsResiduals.nelements());
1392 1005 : itsParentWeights.resize(itsWeights.nelements());
1393 1005 : itsParentImages.resize(itsImages.nelements());
1394 1005 : itsParentSumWts.resize(itsSumWts.nelements());
1395 1005 : itsParentImagePBcors.resize(itsImagePBcors.nelements());
1396 1005 : itsParentPBs.resize(itsPBs.nelements());
1397 :
1398 1005 : itsParentPsfs = itsPsfs;
1399 1005 : itsParentModels=itsModels;
1400 1005 : itsParentResiduals=itsResiduals;
1401 1005 : itsParentWeights=itsWeights;
1402 1005 : itsParentImages=itsImages;
1403 1005 : itsParentSumWts=itsSumWts;
1404 1005 : itsParentImagePBcors=itsImagePBcors;
1405 1005 : itsParentPBs=itsPBs;
1406 :
1407 1005 : itsParentMask=itsMask;
1408 :
1409 1005 : itsParentImageShape = itsImageShape;
1410 1005 : itsParentCoordSys = itsCoordSys;
1411 1005 : if( itsNFacets>1 || itsNChanChunks>1 || itsNPolChunks>1 ) { itsImageShape=IPosition(4,0,0,0,0); }
1412 :
1413 1005 : itsOpened=0;
1414 :
1415 1005 : }
1416 :
1417 : /*
1418 : Bool SIImageStoreMultiTerm::getUseWeightImage()
1419 : { if( itsParentSumWts.nelements()==0 || ! itsParentSumWts[0] )
1420 : {return false;}
1421 : else
1422 : {
1423 : Bool ret = SIImageStore::getUseWeightImage( *(itsParentSumWts[0]) );
1424 : return ret;
1425 : }
1426 : }
1427 : */
1428 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1429 : //////////////////////////////////////////////////////////////////////////////////////////////////////
1430 :
1431 : } //# NAMESPACE CASA - END
1432 :
|
