Line data Source code
1 : //# MultiTermFTNew.cc: Implementation of MultiTermFTNew class
2 : //# Copyright (C) 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: casa-feedback@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 <msvis/MSVis/VisBuffer.h>
29 : #include <msvis/MSVis/VisSet.h>
30 : #include <casacore/images/Images/ImageInterface.h>
31 : #include <casacore/images/Images/PagedImage.h>
32 : #include <casacore/casa/Containers/Block.h>
33 : #include <casacore/casa/Containers/Record.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 : #include <casacore/casa/Arrays/ArrayMath.h>
36 : #include <casacore/casa/Arrays/Array.h>
37 : #include <casacore/casa/Arrays/MaskedArray.h>
38 : #include <casacore/casa/Arrays/Vector.h>
39 : #include <casacore/casa/Arrays/Matrix.h>
40 : #include <casacore/casa/Arrays/Cube.h>
41 : #include <casacore/casa/BasicSL/String.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 : #include <casacore/casa/OS/Timer.h>
45 : #include <sstream>
46 :
47 : #include <synthesis/TransformMachines/StokesImageUtil.h>
48 : #include <synthesis/TransformMachines/VisModelData.h>
49 : #include <casacore/images/Images/ImageInterface.h>
50 : #include <casacore/images/Images/SubImage.h>
51 :
52 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
53 :
54 : #include <synthesis/TransformMachines/MultiTermFTNew.h>
55 : #include <synthesis/TransformMachines/Utils.h>
56 :
57 : // This is the list of FTMachine types supported by MultiTermFTNew
58 : //#include <synthesis/TransformMachines/GridFT.h>
59 : #include <synthesis/TransformMachines/AWProjectFT.h>
60 : #include <synthesis/TransformMachines/AWProjectWBFT.h>
61 :
62 : #include<synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
63 :
64 : using namespace casacore;
65 : namespace casa { //# NAMESPACE CASA - BEGIN
66 :
67 : #define PSOURCE false
68 : #define psource (IPosition(4,1536,1536,0,0))
69 :
70 : //----------------------------------------------------------------------
71 : //-------------------- constructors and descructors ----------------------
72 : //----------------------------------------------------------------------
73 0 : MultiTermFTNew::MultiTermFTNew(CountedPtr<FTMachine>&subftm, Int nterms, Bool forward)
74 0 : :FTMachine(), nterms_p(nterms),
75 0 : reffreq_p(0.0), imweights_p(Matrix<Float>(0,0)), machineName_p("MultiTermFTNew")
76 : // donePSF_p(false)
77 : {
78 :
79 0 : this->setBasePrivates(*subftm);
80 0 : canComputeResiduals_p = subftm->canComputeResiduals();
81 :
82 0 : if( forward ) psfnterms_p = nterms_p;
83 0 : else psfnterms_p = 2*nterms_p-1;
84 :
85 0 : subftms_p.resize(psfnterms_p);
86 0 : for(uInt termindex=0;termindex<psfnterms_p;termindex++)
87 : {
88 : // cout << "Creating new FTM of type : " << subftm->name() << endl;
89 0 : if( termindex==0 ){ subftms_p[termindex] = subftm; }
90 0 : else { subftms_p[termindex] = getNewFTM(subftm); }
91 :
92 0 : subftms_p[termindex]->setMiscInfo(termindex);
93 : }
94 :
95 : // printFTTypes();
96 :
97 0 : }
98 :
99 : //----------------------------------------------------------------------
100 : // Construct from the input state record
101 0 : MultiTermFTNew::MultiTermFTNew(const RecordInterface& stateRec)
102 0 : : FTMachine()
103 : {
104 0 : String error;
105 0 : if (!fromRecord(error, stateRec)) {
106 0 : throw (AipsError("Failed to create gridder: " + error));
107 : };
108 0 : }
109 :
110 : //----------------------------------------------------------------------
111 : // Copy constructor
112 0 : MultiTermFTNew::MultiTermFTNew(const MultiTermFTNew& other) : FTMachine(), machineName_p("MultiTermFTNew")
113 : {
114 0 : operator=(other);
115 0 : }
116 :
117 0 : MultiTermFTNew& MultiTermFTNew::operator=(const MultiTermFTNew& other)
118 : {
119 :
120 0 : if(this!=&other)
121 : {
122 0 : FTMachine::operator=(other);
123 :
124 : // Copy local privates
125 0 : machineName_p = other.machineName_p;
126 0 : nterms_p = other.nterms_p;
127 0 : psfnterms_p = other.psfnterms_p;
128 0 : reffreq_p = other.reffreq_p;
129 : // donePSF_p = other.donePSF_p;
130 :
131 : // Make the list of subftms
132 0 : subftms_p.resize(other.subftms_p.nelements());
133 0 : for (uInt termindex=0;termindex<other.subftms_p.nelements();termindex++)
134 : {
135 0 : subftms_p[termindex] = getNewFTM( (other.subftms_p[termindex]) );
136 0 : subftms_p[termindex]->setMiscInfo(termindex);
137 : }
138 : // subftms_p[termindex] = getNewFTM( &(*(other.subftms_p[termindex])) );
139 :
140 : // Just checking....
141 0 : AlwaysAssert( subftms_p.nelements()>0 , AipsError );
142 :
143 : // Check if the sub ftm type can calculate its own residuals....
144 0 : canComputeResiduals_p = subftms_p[0]->canComputeResiduals();
145 :
146 : }
147 :
148 : // cout << "Checking FTtypes at the end of MTFT operator= for " << ( (other.subftms_p.nelements() > nterms_p)?String("grid"):String("degrid") ) << endl;
149 : // printFTTypes();
150 :
151 0 : return *this;
152 :
153 : }
154 :
155 0 : CountedPtr<FTMachine> MultiTermFTNew::getNewFTM(const CountedPtr<FTMachine>& ftm)
156 : {
157 :
158 0 : return ftm->cloneFTM();
159 :
160 : /*
161 : if(ftm->name()=="GridFT")
162 : {
163 : return new GridFT(static_cast<const GridFT&>(*ftm));
164 : }
165 : else if(ftm->name()=="AWProjectWBFT")
166 : { return new AWProjectWBFT(static_cast<const AWProjectWBFT&>(*ftm)); }
167 : else
168 : {throw(AipsError("FTMachine "+ftm->name()+" is not supported with MS-MFS")); }
169 :
170 : return NULL;
171 : */
172 :
173 : }
174 :
175 0 : FTMachine* MultiTermFTNew::cloneFTM()
176 : {
177 0 : return new MultiTermFTNew(*this);
178 : }
179 :
180 :
181 : //----------------------------------------------------------------------
182 0 : MultiTermFTNew::~MultiTermFTNew()
183 : {
184 0 : }
185 :
186 :
187 : //---------------------------------------------------------------------------------------------------
188 : //------------ Multi-Term Specific Functions --------------------------------------------------------
189 : //---------------------------------------------------------------------------------------------------
190 :
191 : // Multiply the imaging weights by Taylor functions - in place
192 : // This function MUST be called in ascending Taylor-term order
193 : // NOTE : Add checks to ensure this.
194 0 : Bool MultiTermFTNew::modifyVisWeights(VisBuffer& vb,uInt thisterm)
195 : {
196 : {
197 :
198 0 : if(imweights_p.shape() != vb.imagingWeight().shape())
199 0 : imweights_p.resize(vb.imagingWeight().shape());
200 0 : imweights_p = vb.imagingWeight();
201 :
202 0 : Float freq=0.0,mulfactor=1.0;
203 0 : Vector<Double> selfreqlist(vb.frequency());
204 :
205 0 : for (Int row=0; row<vb.nRow(); row++)
206 0 : for (Int chn=0; chn<vb.nChannel(); chn++)
207 : {
208 0 : freq = selfreqlist(IPosition(1,chn));
209 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
210 0 : (vb.imagingWeight())(chn,row) *= pow( mulfactor,(Int)thisterm);
211 : // sumwt_p += (vb.imagingWeight())(chn,row);
212 : }
213 0 : }
214 : /* // For debugging.
215 : else
216 : {
217 : for (Int row=0; row<vb.nRow(); row++)
218 : for (Int chn=0; chn<vb.nChannel(); chn++)
219 : {
220 : sumwt_p += (vb.imagingWeight())(chn,row);
221 : }
222 : }
223 : */
224 0 : return true;
225 : }
226 :
227 0 : void MultiTermFTNew::initMaps(const VisBuffer& vb){
228 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
229 0 : (subftms_p[k])->initMaps(vb);
230 0 : }
231 : // Reset the imaging weights back to their original values
232 : // to be called just after "put"
233 0 : void MultiTermFTNew::restoreImagingWeights(VisBuffer &vb)
234 : {
235 0 : AlwaysAssert( imweights_p.shape() == vb.imagingWeight().shape() ,AipsError);
236 0 : vb.imagingWeight() = imweights_p;
237 0 : }
238 :
239 :
240 : // Multiply the model visibilities by the Taylor functions - in place.
241 0 : Bool MultiTermFTNew::modifyModelVis(VisBuffer& vb, uInt thisterm)
242 : {
243 0 : Float freq=0.0,mulfactor=1.0;
244 0 : Vector<Double> selfreqlist(vb.frequency());
245 :
246 : // DComplex modcount=0.0;
247 :
248 0 : for (uInt pol=0; pol< uInt((vb.modelVisCube()).shape()[0]); pol++)
249 0 : for (uInt chn=0; chn< uInt(vb.nChannel()); chn++)
250 0 : for (uInt row=0; row< uInt(vb.nRow()); row++)
251 : {
252 : // modcount += ( vb.modelVisCube())(pol,chn,row);
253 0 : freq = selfreqlist(IPosition(1,chn));
254 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
255 0 : (vb.modelVisCube())(pol,chn,row) *= pow(mulfactor, (Int) thisterm);
256 : }
257 :
258 : // cout << "field : " << vb.fieldId() << " spw : "
259 : // << vb.spectralWindow() << " --- predicted model before taylor wt mult :"
260 : // << thisterm << " sumvis : " << modcount << endl;
261 :
262 0 : return true;
263 0 : }
264 :
265 :
266 : //---------------------------------------------------------------------------------------------------
267 : //---------------------- Prediction and De-gridding -----------------------------------
268 : //---------------------------------------------------------------------------------------------------
269 :
270 : // void MultiTermFTNew::initializeToVis(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec,PtrBlock<SubImage<Float> *> & modelImageVec, PtrBlock<SubImage<Float> *>& weightImageVec, PtrBlock<SubImage<Float> *>& fluxScaleVec,Block<Matrix<Float> >& weightsVec, const VisBuffer& vb)
271 :
272 0 : void MultiTermFTNew::initializeToVisNew(const VisBuffer& vb,
273 : CountedPtr<SIImageStore> imstore)
274 : {
275 :
276 : // Convert Stokes planes to correlation planes..
277 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
278 : {
279 0 : stokesToCorrelation( *(imstore->model(taylor)) , *(imstore->forwardGrid(taylor) ) );
280 :
281 0 : if(vb.polFrame()==MSIter::Linear) {
282 0 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ), StokesImageUtil::LINEAR);
283 : } else {
284 0 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ) , StokesImageUtil::CIRCULAR);
285 : }
286 : }
287 :
288 0 : reffreq_p = imstore->getReferenceFrequency();
289 :
290 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
291 : {
292 0 : subftms_p[taylor]->initializeToVis(*(imstore->forwardGrid(taylor)),vb);
293 : }
294 :
295 0 : }// end of initializeToVis
296 :
297 :
298 :
299 0 : void MultiTermFTNew::get(VisBuffer& vb, Int row)
300 : {
301 :
302 : // De-grid the model for the zeroth order Taylor term
303 0 : subftms_p[0]->get(vb,row);
304 : // Save the model visibilities in a local cube
305 0 : modviscube_p.assign( vb.modelVisCube() );
306 :
307 0 : for(uInt tix=1;tix<nterms_p;tix++) // Only nterms.... not 2nterms-1
308 : {
309 : // Reset the model visibilities to zero
310 0 : vb.setModelVisCube(Complex(0.0,0.0));
311 : // De-grid the model onto the modelviscube (other Taylor terms)
312 0 : subftms_p[tix]->get(vb,row);
313 : // Multiply visibilities by taylor-weights
314 0 : modifyModelVis(vb,tix);
315 : // Accumulate model visibilities across Taylor terms
316 0 : modviscube_p += vb.modelVisCube();
317 : }
318 : // Set the vb.modelviscube to what has been accumulated
319 0 : vb.setModelVisCube(modviscube_p);
320 0 : }
321 :
322 0 : void MultiTermFTNew::finalizeToVis()
323 : {
324 0 : AlwaysAssert(subftms_p.nelements() >= nterms_p , AipsError);
325 0 : for(uInt taylor=0;taylor<nterms_p;taylor++) subftms_p[taylor]->finalizeToVis();
326 0 : }
327 :
328 : //---------------------------------------------------------------------------------------------------
329 : //---------------------- Calculate Residual Visibilities -------------------------------
330 : //---------------------------------------------------------------------------------------------------
331 0 : void MultiTermFTNew::ComputeResiduals(VisBuffer &vb, Bool useCorrected)
332 : {
333 :
334 0 : if(subftms_p[0]->canComputeResiduals()) subftms_p[0]->ComputeResiduals(vb,useCorrected);
335 0 : else throw(AipsError("MultiTerm::ComputeResiduals : subftm of MultiTermFTNew cannot compute its own residuals !"));
336 :
337 0 : }
338 :
339 : //---------------------------------------------------------------------------------------------------
340 : //---------------------- Gridding --------------------------------------------------------------
341 : //---------------------------------------------------------------------------------------------------
342 :
343 : /// void MultiTermFTNew::initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageV
344 : //ec, Block<Matrix<Float> >& weightsVec, const VisBuffer& vb, const Bool dopsf)
345 :
346 0 : void MultiTermFTNew::initializeToSkyNew(const Bool dopsf,
347 : const VisBuffer& vb,
348 : CountedPtr<SIImageStore> imstore)
349 : {
350 :
351 : // If PSF is already done, don't ask again !
352 : // AlwaysAssert( !(donePSF_p && dopsf) , AipsError );
353 :
354 : // The PSF needs to be the first thing made (because of weight images)
355 : // AlwaysAssert( !(dopsf==false && donePSF_p==false) , AipsError);
356 :
357 : // if(donePSF_p==true)
358 0 : if(dopsf==false)
359 : {
360 0 : if( subftms_p.nelements() != nterms_p )
361 : {
362 0 : subftms_p.resize( nterms_p ,true);
363 : // cout << "MTFT::initializeToSky : resizing to " << nterms_p << " terms" << endl;
364 : }
365 : }
366 :
367 : // Make the relevant float grid.
368 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
369 0 : if( dopsf ) { imstore->psf(0); } else { imstore->residual(0); }
370 :
371 0 : reffreq_p = imstore->getReferenceFrequency();
372 :
373 0 : Matrix<Float> sumWeight;
374 0 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p);taylor++)
375 : {
376 0 : subftms_p[taylor]->initializeToSky(*(imstore->backwardGrid(taylor) ), sumWeight,vb);
377 : }
378 :
379 0 : }// end of initializeToSky
380 :
381 :
382 :
383 0 : void MultiTermFTNew::put(VisBuffer& vb, Int row, Bool dopsf, FTMachine::Type type)
384 : {
385 :
386 0 : subftms_p[0]->put(vb,row,dopsf,type);
387 :
388 0 : if (!dryRun())
389 : {
390 0 : Int gridnterms=nterms_p;
391 0 : if(dopsf==true) // && donePSF_p==false)
392 : {
393 0 : gridnterms=2*nterms_p-1;
394 : }
395 :
396 : //cout << " Calling put for " << gridnterms << " terms, nelements : " << subftms_p.nelements() << " and dopsf " << dopsf << endl;
397 :
398 0 : for(Int tix=1;tix<gridnterms;tix++)
399 : {
400 0 : modifyVisWeights(vb,tix);
401 0 : subftms_p[tix]->put(vb,row,dopsf,type);
402 0 : restoreImagingWeights(vb);
403 : }
404 : }
405 :
406 :
407 0 : }// end of put
408 :
409 : //----------------------------------------------------------------------
410 :
411 : // void MultiTermFTNew::finalizeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec, PtrBlock<SubImage<Float> *> & resImageVec, PtrBlock<SubImage<Float> *>& weightImageVec, PtrBlock<SubImage<Float> *>& fluxScaleVec, Bool dopsf, Block<Matrix<Float> >& weightsVec, const VisBuffer& /*vb*/)
412 :
413 0 : void MultiTermFTNew::finalizeToSkyNew(Bool dopsf,
414 : const VisBuffer& /*vb*/,
415 : CountedPtr<SIImageStore> imstore )
416 : {
417 :
418 : // Collect images and weights from all FTMs
419 0 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p) ;taylor++)
420 : {
421 0 : Matrix<Float> sumWeights;
422 0 : subftms_p[taylor]->finalizeToSky();
423 0 : correlationToStokes( subftms_p[taylor]->getImage(sumWeights, false) , ( dopsf ? *(imstore->psf(taylor)) : *(imstore->residual(taylor)) ), dopsf);
424 :
425 0 : if( subftms_p[taylor]->useWeightImage() && dopsf ) {
426 0 : subftms_p[taylor]->getWeightImage(*(imstore->weight(taylor)), sumWeights);
427 : }
428 :
429 : // Take sumWeights from corrToStokes here....
430 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
431 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
432 :
433 0 : AlwaysAssert( ( (imstore->sumwt(taylor))->shape()[2] == sumWeightStokes.shape()[0] ) &&
434 : ((imstore->sumwt(taylor))->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
435 :
436 0 : (imstore->sumwt(taylor))->put( sumWeightStokes.reform((imstore->sumwt(taylor))->shape()) );
437 :
438 : // cout << "taylor : " << taylor << " sumwt : " << sumWeights << endl;
439 :
440 0 : }// end for taylor
441 :
442 : // if( dopsf ) donePSF_p = true;
443 :
444 0 : }//end of finalizeToSkyNew
445 :
446 :
447 : //---------------------------------------------------------------------------------------------------
448 : //----------------------------- Obtain Images -----------------------------------------------------
449 : //---------------------------------------------------------------------------------------------------
450 : //----------------------------------------------------------------------
451 0 : void MultiTermFTNew::makeImage(FTMachine::Type type, VisSet& vs,
452 : ImageInterface<Complex>& theImage, Matrix<Float>& weight)
453 : {
454 : // cout << "MTFT :: makeImage for taylor 0 only "<< endl;
455 0 : subftms_p[0]->makeImage(type, vs, theImage, weight);
456 0 : }
457 :
458 :
459 : //---------------------------------------------------------------------------------------------------
460 : //------------------------ To / From Records ---------------------------------------------------------
461 : //---------------------------------------------------------------------------------------------------
462 0 : Bool MultiTermFTNew::toRecord(String& error, RecordInterface& outRec, Bool withImage, const String diskimage)
463 : {
464 : // cout << "MTFTNew :: toRecord for " << subftms_p.nelements() << " subftms" << endl;
465 0 : Bool retval = true;
466 0 : outRec.define("name", this->name());
467 0 : outRec.define("nterms",nterms_p);
468 0 : outRec.define("reffreq",reffreq_p);
469 0 : outRec.define("machinename",machineName_p);
470 0 : outRec.define("psfnterms",psfnterms_p);
471 : // outRec.define("donePSF_p",donePSF_p);
472 :
473 0 : outRec.define("numfts", (Int)subftms_p.nelements() ); // Since the forward and reverse ones are different.
474 :
475 0 : for(uInt tix=0;tix<subftms_p.nelements();tix++)
476 : {
477 0 : Record subFTContainer;
478 0 : String elimage="";
479 0 : if(diskimage != ""){
480 0 : elimage=diskimage+String("_term_")+String::toString(tix);
481 : }
482 0 : subftms_p[tix]->toRecord(error, subFTContainer,withImage, elimage);
483 0 : outRec.defineRecord("subftm_"+String::toString(tix),subFTContainer);
484 0 : }
485 :
486 0 : return retval;
487 : }
488 :
489 : //---------------------------------------------------------------------------------------------------
490 0 : Bool MultiTermFTNew::fromRecord(String& error, const RecordInterface& inRec)
491 : {
492 : // cout << "MTFTNew :: fromRecord "<< endl;
493 0 : Bool retval = true;
494 :
495 0 : inRec.get("nterms",nterms_p);
496 0 : inRec.get("reffreq",reffreq_p);
497 0 : inRec.get("machinename",machineName_p);
498 0 : inRec.get("psfnterms",psfnterms_p);
499 : // inRec.get("donePSF_p",donePSF_p);
500 :
501 0 : Int nftms=1;
502 0 : inRec.get("numfts",nftms);
503 :
504 0 : subftms_p.resize(nftms);
505 0 : for(Int tix=0;tix<nftms;tix++)
506 : {
507 0 : Record subFTMRec=inRec.asRecord("subftm_"+String::toString(tix));
508 0 : subftms_p[tix]=VisModelData::NEW_FT(subFTMRec);
509 0 : retval = (retval || subftms_p[tix]->fromRecord(error, subFTMRec));
510 0 : }
511 :
512 :
513 0 : return retval;
514 : }
515 : //---------------------------------------------------------------------------------------------------
516 :
517 0 : Bool MultiTermFTNew::storeAsImg(String fileName, ImageInterface<Float>& theImg)
518 : {
519 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
520 0 : LatticeExpr<Float> le(theImg);
521 0 : tmp.copyData(le);
522 0 : return true;
523 0 : }
524 :
525 :
526 :
527 : //---------------------------------------------------------------------------------------------------
528 : //---------------------------------------------------------------------------------------------------
529 : //---------------------------------------------------------------------------------------------------
530 :
531 : } //# NAMESPACE CASA - END
532 :
|
