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/VisBuffer2.h>
29 : #include <casacore/images/Images/ImageInterface.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/Containers/Block.h>
32 : #include <casacore/casa/Containers/Record.h>
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Arrays/ArrayMath.h>
35 : #include <casacore/casa/Arrays/Array.h>
36 : #include <casacore/casa/Arrays/MaskedArray.h>
37 : #include <casacore/casa/Arrays/Vector.h>
38 : #include <casacore/casa/Arrays/Matrix.h>
39 : #include <casacore/casa/Arrays/Cube.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 : #include <casacore/casa/Exceptions/Error.h>
43 : #include <casacore/casa/OS/Timer.h>
44 : #include <sstream>
45 :
46 : #include <synthesis/TransformMachines/StokesImageUtil.h>
47 : #include <synthesis/TransformMachines2/VisModelData.h>
48 : #include <casacore/lattices/Lattices/LatticeLocker.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/TransformMachines2/MultiTermFTNew.h>
55 : #include <synthesis/TransformMachines2/MosaicFTNew.h>
56 : #include <synthesis/TransformMachines2/Utils.h>
57 :
58 : // This is the list of FTMachine types supported by MultiTermFTNew
59 : //#include <synthesis/TransformMachines/GridFT.h>
60 : //#include <synthesis/TransformMachines/AWProjectFT.h>
61 : //#include <synthesis/TransformMachines/AWProjectWBFT.h>
62 :
63 : #include<synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 : namespace refim { // namespace for refactor
68 :
69 : using namespace casacore;
70 : using namespace casa;
71 : using namespace casacore;
72 : using namespace casa::refim;
73 : using namespace casacore;
74 : using namespace casa::vi;
75 :
76 : #define PSOURCE false
77 : #define psource (IPosition(4,1536,1536,0,0))
78 :
79 : //----------------------------------------------------------------------
80 : //-------------------- constructors and descructors ----------------------
81 : //----------------------------------------------------------------------
82 0 : MultiTermFTNew::MultiTermFTNew(CountedPtr<FTMachine>&subftm, Int nterms, Bool forward)
83 0 : :FTMachine(), nterms_p(nterms),
84 0 : reffreq_p(0.0), imweights_p(Matrix<Float>(0,0)), machineName_p("MultiTermFTNew")
85 : // donePSF_p(false)
86 : {
87 :
88 0 : this->setBasePrivates(*subftm);
89 0 : canComputeResiduals_p = subftm->canComputeResiduals();
90 :
91 0 : if( forward ) psfnterms_p = nterms_p;
92 0 : else psfnterms_p = 2*nterms_p-1;
93 :
94 0 : subftms_p.resize(psfnterms_p);
95 0 : for(uInt termindex=0;termindex<psfnterms_p;termindex++)
96 : {
97 : // cout << "Creating new FTM of type : " << subftm->name() << endl;
98 0 : if( termindex==0 ){
99 0 : subftms_p[termindex] = subftm;
100 : }
101 : else {
102 0 : subftms_p[termindex] = getNewFTM(subftm);
103 0 : if((subftms_p[termindex]->name())=="MosaicFTNew"){
104 0 : (static_cast<MosaicFTNew *>(subftms_p[termindex].get()))->setConvFunc( (static_cast<MosaicFTNew * >(subftm.get()))->getConvFunc());
105 :
106 : }
107 : }
108 :
109 0 : subftms_p[termindex]->setMiscInfo(termindex);
110 : }
111 :
112 : // printFTTypes();
113 0 : }
114 :
115 : //----------------------------------------------------------------------
116 : // Construct from the input state record
117 0 : MultiTermFTNew::MultiTermFTNew(const RecordInterface& stateRec)
118 0 : : FTMachine()
119 : {
120 0 : String error;
121 0 : if (!fromRecord(error, stateRec)) {
122 0 : throw (AipsError("Failed to create gridder: " + error));
123 : };
124 0 : }
125 :
126 : //----------------------------------------------------------------------
127 : // Copy constructor
128 0 : MultiTermFTNew::MultiTermFTNew(const MultiTermFTNew& other) : FTMachine(), machineName_p("MultiTermFTNew")
129 : {
130 0 : operator=(other);
131 0 : }
132 :
133 0 : MultiTermFTNew& MultiTermFTNew::operator=(const MultiTermFTNew& other)
134 : {
135 :
136 0 : if(this!=&other)
137 : {
138 0 : FTMachine::operator=(other);
139 :
140 : // Copy local privates
141 0 : machineName_p = other.machineName_p;
142 0 : nterms_p = other.nterms_p;
143 0 : psfnterms_p = other.psfnterms_p;
144 0 : reffreq_p = other.reffreq_p;
145 : // donePSF_p = other.donePSF_p;
146 :
147 : // Make the list of subftms
148 0 : subftms_p.resize(other.subftms_p.nelements());
149 0 : for (uInt termindex=0;termindex<other.subftms_p.nelements();termindex++)
150 : {
151 0 : subftms_p[termindex] = getNewFTM( (other.subftms_p[termindex]) );
152 0 : subftms_p[termindex]->setMiscInfo(termindex);
153 : }
154 : // subftms_p[termindex] = getNewFTM( &(*(other.subftms_p[termindex])) );
155 :
156 : // Just checking....
157 0 : AlwaysAssert( subftms_p.nelements()>0 , AipsError );
158 :
159 : // Check if the sub ftm type can calculate its own residuals....
160 0 : canComputeResiduals_p = subftms_p[0]->canComputeResiduals();
161 :
162 : }
163 :
164 : // cout << "Checking FTtypes at the end of MTFT operator= for " << ( (other.subftms_p.nelements() > nterms_p)?String("grid"):String("degrid") ) << endl;
165 : // printFTTypes();
166 :
167 0 : return *this;
168 :
169 : }
170 :
171 0 : CountedPtr<FTMachine> MultiTermFTNew::getNewFTM(const CountedPtr<FTMachine>& ftm)
172 : {
173 :
174 0 : return ftm->cloneFTM();
175 :
176 : /*
177 : if(ftm->name()=="GridFT")
178 : {
179 : return new GridFT(static_cast<const GridFT&>(*ftm));
180 : }
181 : else if(ftm->name()=="AWProjectWBFT")
182 : { return new AWProjectWBFT(static_cast<const AWProjectWBFT&>(*ftm)); }
183 : else
184 : {throw(AipsError("FTMachine "+ftm->name()+" is not supported with MS-MFS")); }
185 :
186 : return NULL;
187 : */
188 :
189 : }
190 :
191 0 : FTMachine* MultiTermFTNew::cloneFTM()
192 : {
193 0 : return new MultiTermFTNew(*this);
194 : }
195 :
196 :
197 : //----------------------------------------------------------------------
198 0 : MultiTermFTNew::~MultiTermFTNew()
199 : {
200 0 : }
201 0 : void MultiTermFTNew::setMovingSource(const String& sourcename, const String& ephemtable){
202 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
203 0 : (subftms_p[k])->setMovingSource(sourcename, ephemtable);
204 0 : }
205 0 : void MultiTermFTNew::setMovingSource(const MDirection& mdir){
206 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
207 0 : (subftms_p[k])->setMovingSource(mdir);
208 0 : }
209 0 : void MultiTermFTNew::setLocation(const MPosition& mloc){
210 0 : for (uInt k=0; k < subftms_p.nelements(); ++k)
211 0 : (subftms_p[k])->setLocation(mloc);
212 0 : }
213 : //---------------------------------------------------------------------------------------------------
214 : //------------ Multi-Term Specific Functions --------------------------------------------------------
215 : //---------------------------------------------------------------------------------------------------
216 :
217 : // Multiply the imaging weights by Taylor functions - in place
218 : // This function MUST be called in ascending Taylor-term order
219 : // NOTE : Add checks to ensure this.
220 0 : Bool MultiTermFTNew::modifyVisWeights(VisBuffer2& vb,uInt thisterm)
221 : {
222 : {
223 : ///There is no setImagingWeight in vb2 !
224 0 : Matrix<Float>& imwgt=const_cast<Matrix<Float>& >(vb.imagingWeight());
225 : ////remove the above line when using setImagingWeight
226 :
227 0 : if(imweights_p.shape() != vb.imagingWeight().shape())
228 0 : imweights_p.resize(vb.imagingWeight().shape());
229 0 : imweights_p = vb.imagingWeight();
230 :
231 0 : Float freq=0.0,mulfactor=1.0;
232 0 : Vector<Double> selfreqlist(vb.getFrequencies(0));
233 :
234 0 : for (rownr_t row=0; row<vb.nRows(); row++)
235 0 : for (Int chn=0; chn<vb.nChannels(); chn++)
236 : {
237 0 : freq = selfreqlist(IPosition(1,chn));
238 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
239 0 : (imwgt)(chn,row) *= pow( mulfactor,(Int)thisterm);
240 : // sumwt_p += (vb.imagingWeight())(chn,row);
241 : }
242 0 : vb.setImagingWeight(imwgt);
243 0 : }
244 : /* // For debugging.
245 : else
246 : {
247 : for (Int row=0; row<vb.nRow(); row++)
248 : for (Int chn=0; chn<vb.nChannel(); chn++)
249 : {
250 : sumwt_p += (vb.imagingWeight())(chn,row);
251 : }
252 : }
253 : */
254 0 : return true;
255 : }
256 :
257 0 : void MultiTermFTNew::initMaps(const VisBuffer2& vb){
258 :
259 0 : for (uInt k=0; k < subftms_p.nelements(); ++k){
260 0 : (subftms_p[k])->setFrameValidity( freqFrameValid_p);
261 0 : (subftms_p[k])->freqInterpMethod_p=this->freqInterpMethod_p;
262 0 : (subftms_p[k])->initMaps(vb);
263 : }
264 0 : }
265 : // Reset the imaging weights back to their original values
266 : // to be called just after "put"
267 0 : void MultiTermFTNew::restoreImagingWeights(VisBuffer2 &vb)
268 : {
269 0 : AlwaysAssert( imweights_p.shape() == vb.imagingWeight().shape() ,AipsError);
270 : ///There is no setImagingWeight in vb2 !
271 : // Matrix<Float>& imwgt=const_cast<Matrix<Float>& >(vb.imagingWeight());
272 : // ////remove the above line when using setImagingWeight
273 : // imwgt = imweights_p;
274 0 : vb.setImagingWeight(imweights_p);
275 0 : }
276 :
277 :
278 : // Multiply the model visibilities by the Taylor functions - in place.
279 0 : Bool MultiTermFTNew::modifyModelVis(VisBuffer2& vb, uInt thisterm)
280 : {
281 0 : Float freq=0.0,mulfactor=1.0;
282 0 : Vector<Double> selfreqlist(vb.getFrequencies(0));
283 :
284 : // DComplex modcount=0.0;
285 :
286 0 : Cube<Complex> modelVisCube(vb.visCubeModel().shape());
287 0 : modelVisCube=vb.visCubeModel();
288 0 : for (uInt pol=0; pol< uInt((vb.visCubeModel()).shape()[0]); pol++)
289 0 : for (uInt chn=0; chn< uInt(vb.nChannels()); chn++)
290 0 : for (uInt row=0; row< uInt(vb.nRows()); row++)
291 : {
292 : // modcount += ( vb.modelVisCube())(pol,chn,row);
293 0 : freq = selfreqlist(IPosition(1,chn));
294 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
295 0 : modelVisCube(pol,chn,row) *= pow(mulfactor, (Int) thisterm);
296 : }
297 0 : vb.setVisCubeModel(modelVisCube);
298 : // cout << "field : " << vb.fieldId() << " spw : "
299 : // << vb.spectralWindow() << " --- predicted model before taylor wt mult :"
300 : // << thisterm << " sumvis : " << modcount << endl;
301 :
302 0 : return true;
303 0 : }
304 :
305 :
306 : //---------------------------------------------------------------------------------------------------
307 : //---------------------- Prediction and De-gridding -----------------------------------
308 : //---------------------------------------------------------------------------------------------------
309 :
310 : // 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)
311 :
312 0 : void MultiTermFTNew::initializeToVisNew(const VisBuffer2& vb,
313 : CountedPtr<SIImageStore> imstore)
314 : {
315 :
316 : // Convert Stokes planes to correlation planes..
317 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
318 : {
319 :
320 0 : if(!(imstore->forwardGrid(taylor)).get())
321 0 : throw(AipsError("MultiTermFTNew::InitializeToVisNew error imagestore has no valid grid initialized for taylor term "+String::toString(taylor)));
322 0 : stokesToCorrelation( *(imstore->model(taylor)) , *(imstore->forwardGrid(taylor) ) );
323 :
324 0 : if(vb.polarizationFrame()==MSIter::Linear) {
325 0 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ), StokesImageUtil::LINEAR);
326 : } else {
327 0 : StokesImageUtil::changeCStokesRep( *(imstore->forwardGrid(taylor) ) , StokesImageUtil::CIRCULAR);
328 : }
329 : }
330 :
331 0 : reffreq_p = imstore->getReferenceFrequency();
332 :
333 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
334 : {
335 0 : subftms_p[taylor]->initializeToVis(*(imstore->forwardGrid(taylor)),vb);
336 : }
337 :
338 0 : }// end of initializeToVis
339 :
340 :
341 :
342 0 : void MultiTermFTNew::get(VisBuffer2& vb, Int row)
343 : {
344 :
345 : // De-grid the model for the zeroth order Taylor term
346 0 : subftms_p[0]->get(vb,row);
347 : // Save the model visibilities in a local cube
348 0 : modviscube_p.assign( vb.visCubeModel() );
349 :
350 0 : for(uInt tix=1;tix<nterms_p;tix++) // Only nterms.... not 2nterms-1
351 : {
352 : // Reset the model visibilities to zero
353 0 : vb.setVisCubeModel(Complex(0.0,0.0));
354 : // De-grid the model onto the modelviscube (other Taylor terms)
355 0 : subftms_p[tix]->get(vb,row);
356 : // Multiply visibilities by taylor-weights
357 0 : modifyModelVis(vb,tix);
358 : // Accumulate model visibilities across Taylor terms
359 0 : modviscube_p += vb.visCubeModel();
360 : }
361 : // Set the vb.modelviscube to what has been accumulated
362 0 : vb.setVisCubeModel(modviscube_p);
363 0 : }
364 :
365 0 : void MultiTermFTNew::finalizeToVis()
366 : {
367 0 : AlwaysAssert(subftms_p.nelements() >= nterms_p , AipsError);
368 0 : for(uInt taylor=0;taylor<nterms_p;taylor++) subftms_p[taylor]->finalizeToVis();
369 0 : }
370 :
371 : //---------------------------------------------------------------------------------------------------
372 : //---------------------- Calculate Residual Visibilities -------------------------------
373 : //---------------------------------------------------------------------------------------------------
374 0 : void MultiTermFTNew::ComputeResiduals(VisBuffer2 &vb, Bool useCorrected)
375 : {
376 :
377 0 : if(subftms_p[0]->canComputeResiduals()) subftms_p[0]->ComputeResiduals(vb,useCorrected);
378 0 : else throw(AipsError("MultiTerm::ComputeResiduals : subftm of MultiTermFTNew cannot compute its own residuals !"));
379 :
380 0 : }
381 :
382 : //---------------------------------------------------------------------------------------------------
383 : //---------------------- Gridding --------------------------------------------------------------
384 : //---------------------------------------------------------------------------------------------------
385 :
386 : /// void MultiTermFTNew::initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageV
387 : //ec, Block<Matrix<Float> >& weightsVec, const VisBuffer& vb, const Bool dopsf)
388 0 : Long MultiTermFTNew::estimateRAM(const CountedPtr<SIImageStore>& imstor){
389 0 : Long mem=0;
390 0 : for(uInt k=0; k < subftms_p.nelements(); ++k){
391 :
392 0 : mem+=subftms_p[k]->estimateRAM(imstor);
393 : }
394 :
395 0 : return mem;
396 : }
397 :
398 0 : void MultiTermFTNew::initializeToSkyNew(const Bool dopsf,
399 : const VisBuffer2& vb,
400 : CountedPtr<SIImageStore> imstore)
401 : {
402 :
403 : // If PSF is already done, don't ask again !
404 : // AlwaysAssert( !(donePSF_p && dopsf) , AipsError );
405 :
406 : // The PSF needs to be the first thing made (because of weight images)
407 : // AlwaysAssert( !(dopsf==false && donePSF_p==false) , AipsError);
408 :
409 : // if(donePSF_p==true)
410 0 : if(dopsf==false)
411 : {
412 0 : if( subftms_p.nelements() != nterms_p )
413 : {
414 0 : subftms_p.resize( nterms_p ,true);
415 : // cout << "MTFT::initializeToSky : resizing to " << nterms_p << " terms" << endl;
416 : }
417 : }
418 :
419 : // Make the relevant float grid.
420 : // This is needed mainly for facetting (to set facet shapes), but is harmless for non-facetting.
421 0 : if( dopsf ) { imstore->psf(0); } else { imstore->residual(0); }
422 :
423 0 : reffreq_p = imstore->getReferenceFrequency();
424 :
425 0 : Matrix<Float> sumWeight;
426 0 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p);taylor++)
427 : {
428 :
429 0 : if(! (imstore->backwardGrid(taylor)).get())
430 0 : throw(AipsError("MultiTermFTNew::InitializeToSkyNew error imagestore has no valid grid initialized for taylor term "+String::toString(taylor)));
431 0 : subftms_p[taylor]->initializeToSky(*(imstore->backwardGrid(taylor) ), sumWeight,vb);
432 : }
433 :
434 0 : }// end of initializeToSky
435 :
436 0 : void MultiTermFTNew::initBriggsWeightor(vi::VisibilityIterator2& vi){
437 :
438 :
439 0 : for (uInt k=0; k < subftms_p.nelements(); ++k){
440 0 : subftms_p[k]->initBriggsWeightor(vi);
441 : }
442 :
443 0 : }
444 :
445 0 : void MultiTermFTNew::put(VisBuffer2& vb, Int row, Bool dopsf, refim::FTMachine::Type type)
446 : {
447 :
448 0 : subftms_p[0]->put(vb,row,dopsf,type);
449 0 : if (!dryRun())
450 : {
451 0 : Int gridnterms=nterms_p;
452 0 : if(dopsf==true) // && donePSF_p==false)
453 : {
454 0 : gridnterms=2*nterms_p-1;
455 : }
456 :
457 : //cerr << " Calling put for " << gridnterms << " terms, nelements : " << subftms_p.nelements() << " and dopsf " << dopsf << " reffreq " << reffreq_p << endl;
458 :
459 0 : for(Int tix=1;tix<gridnterms;tix++)
460 : {
461 0 : modifyVisWeights(vb,tix);
462 0 : subftms_p[tix]->put(vb,row,dopsf,type);
463 0 : restoreImagingWeights(vb);
464 : }
465 : }
466 :
467 0 : }// end of put
468 :
469 : //----------------------------------------------------------------------
470 :
471 : ///-------------------
472 0 : void MultiTermFTNew::gridImgWeights(const VisBuffer2& vb)
473 : {
474 :
475 0 : subftms_p[0]->gridImgWeights(vb);
476 0 : if (!dryRun())
477 : {
478 0 : Int gridnterms=2*nterms_p-1;
479 :
480 :
481 : //cerr << " Calling put for " << gridnterms << " terms, nelements : " << subftms_p.nelements() << " and dopsf " << dopsf << " reffreq " << reffreq_p << endl;
482 :
483 0 : for(Int tix=1;tix<gridnterms;tix++)
484 : {
485 0 : modifyVisWeights(const_cast<vi::VisBuffer2&>(vb),tix);
486 0 : subftms_p[tix]->gridImgWeights(vb);
487 0 : restoreImagingWeights(const_cast<vi::VisBuffer2&>(vb));
488 : }
489 : }
490 :
491 0 : }// end of gridImgWeights
492 :
493 :
494 :
495 :
496 : ////---------------------------
497 : // 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*/)
498 :
499 0 : void MultiTermFTNew::finalizeToSkyNew(Bool dopsf,
500 : const VisBuffer2& /*vb*/,
501 : CountedPtr<SIImageStore> imstore )
502 : {
503 :
504 : // Collect images and weights from all FTMs
505 0 : for(uInt taylor=0;taylor< (dopsf ? psfnterms_p : nterms_p) ;taylor++)
506 : {
507 0 : Matrix<Float> sumWeights;
508 0 : subftms_p[taylor]->finalizeToSky();
509 0 : shared_ptr<ImageInterface<Float> > theim=dopsf ? imstore->psf(taylor) : imstore->residual(taylor);
510 0 : { LatticeLocker lock1 (*theim, FileLocker::Write);
511 0 : correlationToStokes( subftms_p[taylor]->getImage(sumWeights, false) , *theim, dopsf);
512 0 : }
513 0 : if( subftms_p[taylor]->useWeightImage() && dopsf ) {
514 0 : LatticeLocker lock1 (*(imstore->weight(taylor)), FileLocker::Write);
515 0 : subftms_p[taylor]->getWeightImage(*(imstore->weight(taylor)), sumWeights);
516 0 : }
517 :
518 : // Take sumWeights from corrToStokes here....
519 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt(taylor))->shape()[2], (imstore->sumwt(taylor))->shape()[3] );
520 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
521 :
522 0 : AlwaysAssert( ( (imstore->sumwt(taylor))->shape()[2] == sumWeightStokes.shape()[0] ) &&
523 : ((imstore->sumwt(taylor))->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
524 0 : LatticeLocker lock1 (*(imstore->sumwt(taylor)), FileLocker::Write);
525 0 : (imstore->sumwt(taylor))->put( sumWeightStokes.reform((imstore->sumwt(taylor))->shape()) );
526 :
527 : // cout << "taylor : " << taylor << " sumwt : " << sumWeights << endl;
528 :
529 0 : }// end for taylor
530 :
531 : // if( dopsf ) donePSF_p = true;
532 :
533 0 : }//end of finalizeToSkyNew
534 :
535 :
536 : //---------------------------------------------------------------------------------------------------
537 :
538 : //////-------------------------------------------------------------------
539 0 : void MultiTermFTNew::finalizeToWeightImage(const VisBuffer2& /*vb*/,
540 : CountedPtr<SIImageStore> imstore ){
541 :
542 0 : for(uInt taylor=0;taylor< psfnterms_p ;taylor++)
543 : {
544 0 : Matrix<Float> sumWeights;
545 0 : subftms_p[taylor]->finalizeToSky();
546 0 : if( subftms_p[taylor]->useWeightImage() )
547 : {
548 0 : LatticeLocker lock1 (*(imstore->weight(taylor)), FileLocker::Write);
549 0 : subftms_p[taylor]->getWeightImage(*(imstore->weight(taylor)), sumWeights);
550 :
551 :
552 : // Take sumWeights from corrToStokes here....
553 0 : Matrix<Float> sumWeightStokes( (imstore->sumwt())->shape()[2], (imstore->sumwt())->shape()[3] );
554 0 : StokesImageUtil::ToStokesSumWt( sumWeightStokes, sumWeights );
555 :
556 0 : AlwaysAssert( ( (imstore->sumwt(taylor))->shape()[2] == sumWeightStokes.shape()[0] ) &&
557 : ((imstore->sumwt(taylor))->shape()[3] == sumWeightStokes.shape()[1] ) , AipsError );
558 0 : LatticeLocker lock2 (*(imstore->sumwt(taylor)), FileLocker::Write);
559 0 : (imstore->sumwt(taylor))->put( sumWeightStokes.reform((imstore->sumwt(taylor))->shape()) );
560 0 : }//useWeight
561 : // cout << "taylor : " << taylor << " sumwt : " << sumWeights << endl;
562 :
563 0 : }// end for taylor
564 :
565 :
566 :
567 :
568 0 : }
569 : ////////---------------------------------------------
570 : //----------------------------- Obtain Images -----------------------------------------------------
571 : //---------------------------------------------------------------------------------------------------
572 : //----------------------------------------------------------------------
573 0 : void MultiTermFTNew::makeImage(refim::FTMachine::Type type, VisibilityIterator2& vi,
574 : ImageInterface<Complex>& theImage, Matrix<Float>& weight)
575 : {
576 : // cout << "MTFT :: makeImage for taylor 0 only "<< endl;
577 0 : subftms_p[0]->makeImage(type, vi, theImage, weight);
578 0 : }
579 0 : void MultiTermFTNew::makeMTImages(refim::FTMachine::Type type,
580 : vi::VisibilityIterator2& vi,
581 : casacore::Vector<casacore::CountedPtr<casacore::ImageInterface<casacore::Complex> > >& theImage,
582 : casacore::Vector<casacore::CountedPtr<casacore::Matrix<casacore::Float> > >& weight){
583 0 : Int ntaylor= (type== refim::FTMachine::PSF) ? psfnterms_p : nterms_p;
584 :
585 :
586 :
587 0 : vi::VisBuffer2* vb=vi.getVisBuffer();
588 :
589 : // Initialize put (i.e. transform to Sky) for this model
590 0 : vi.origin();
591 0 : for(Int taylor=0;taylor< ntaylor ; ++taylor) {
592 0 : if(vb->polarizationFrame()==MSIter::Linear) {
593 0 : StokesImageUtil::changeCStokesRep(*(theImage[taylor]), StokesImageUtil::LINEAR);
594 : }
595 : else {
596 0 : StokesImageUtil::changeCStokesRep(*(theImage[taylor]), StokesImageUtil::CIRCULAR);
597 : }
598 0 : subftms_p[taylor]->initializeToSky(*(theImage[taylor]), *(weight[taylor]),*vb);
599 :
600 : }
601 : {
602 0 : Vector<Double> refpix = (theImage[0]->coordinates().spectralCoordinate()).referencePixel();
603 0 : (theImage[0]->coordinates().spectralCoordinate()).toWorld( reffreq_p, refpix[0] );
604 0 : }
605 0 : Bool useCorrected= !(MSColumns(vi.ms()).correctedData().isNull());
606 0 : if((type==FTMachine::CORRECTED) && (!useCorrected))
607 0 : type=FTMachine::OBSERVED;
608 :
609 :
610 : // Loop over the visibilities, putting VisBuffers
611 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
612 0 : for (vi.origin(); vi.more(); vi.next()) {
613 :
614 0 : switch(type) {
615 0 : case FTMachine::RESIDUAL:
616 0 : vb->setVisCube(vb->visCubeCorrected());
617 0 : vb->setVisCube(vb->visCube()-vb->visCubeModel());
618 0 : put(*vb, -1, false);
619 0 : break;
620 0 : case FTMachine::CORRECTED:
621 0 : put(*vb, -1, false, FTMachine::CORRECTED);
622 0 : break;
623 0 : case FTMachine::PSF:
624 0 : vb->setVisCube(Complex(1.0,0.0));
625 0 : put(*vb, -1, true, FTMachine::PSF);
626 0 : break;
627 0 : case FTMachine::OBSERVED:
628 0 : put(*vb, -1, false, FTMachine::OBSERVED);
629 0 : break;
630 0 : default:
631 0 : throw(AipsError("Cannot make multiterm image of type requested"));
632 : break;
633 : }
634 : }
635 : }
636 : ///
637 0 : for(Int taylor=0;taylor< ntaylor ; ++taylor) {
638 0 : subftms_p[taylor]->finalizeToSky();
639 0 : subftms_p[taylor]->getImage(*(weight[taylor]), false);
640 :
641 : }
642 0 : }
643 : //---------------------------------------------------------------------------------------------------
644 : //------------------------ To / From Records ---------------------------------------------------------
645 : //---------------------------------------------------------------------------------------------------
646 0 : Bool MultiTermFTNew::toRecord(String& error, RecordInterface& outRec, Bool withImage, const String diskimage)
647 : {
648 : // cout << "MTFTNew :: toRecord for " << subftms_p.nelements() << " subftms" << endl;
649 0 : Bool retval = true;
650 0 : outRec.define("name", this->name());
651 0 : outRec.define("nterms",nterms_p);
652 0 : outRec.define("reffreq",reffreq_p);
653 0 : outRec.define("machinename",machineName_p);
654 0 : outRec.define("psfnterms",psfnterms_p);
655 : // outRec.define("donePSF_p",donePSF_p);
656 :
657 0 : outRec.define("numfts", (Int)subftms_p.nelements() ); // Since the forward and reverse ones are different.
658 :
659 0 : for(uInt tix=0;tix<subftms_p.nelements();tix++)
660 : {
661 0 : Record subFTContainer;
662 0 : String elimage="";
663 0 : if(diskimage != ""){
664 0 : elimage=diskimage+String("_term_")+String::toString(tix);
665 : }
666 0 : subftms_p[tix]->toRecord(error, subFTContainer,withImage, elimage);
667 0 : outRec.defineRecord("subftm_"+String::toString(tix),subFTContainer);
668 0 : }
669 :
670 0 : return retval;
671 : }
672 :
673 : //---------------------------------------------------------------------------------------------------
674 0 : Bool MultiTermFTNew::fromRecord(String& error, const RecordInterface& inRec)
675 : {
676 : //cout << "MTFTNew :: fromRecord "<< endl;
677 0 : Bool retval = true;
678 :
679 0 : inRec.get("nterms",nterms_p);
680 0 : inRec.get("reffreq",reffreq_p);
681 0 : inRec.get("machinename",machineName_p);
682 0 : inRec.get("psfnterms",psfnterms_p);
683 : // inRec.get("donePSF_p",donePSF_p);
684 :
685 0 : Int nftms=1;
686 0 : inRec.get("numfts",nftms);
687 :
688 0 : subftms_p.resize(nftms);
689 : //cerr << "number of ft " << nftms << endl;
690 0 : for(Int tix=0;tix<nftms;tix++)
691 : {
692 0 : Record subFTMRec=inRec.asRecord("subftm_"+String::toString(tix));
693 0 : subftms_p[tix]=VisModelData::NEW_FT(subFTMRec);
694 0 : if(subftms_p[tix].null())
695 0 : throw(AipsError("Could not recover from record term "+String::toString(tix)+" ftmachine"));
696 0 : retval = (retval || subftms_p[tix]->fromRecord(error, subFTMRec));
697 0 : if(!retval) throw(AipsError("Could not recover term "+String::toString(tix)+" ftmachine; \n Error being "+error));
698 0 : }
699 :
700 :
701 0 : return retval;
702 : }
703 : //---------------------------------------------------------------------------------------------------
704 :
705 0 : Bool MultiTermFTNew::storeAsImg(String fileName, ImageInterface<Float>& theImg)
706 : {
707 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
708 0 : LatticeExpr<Float> le(theImg);
709 0 : tmp.copyData(le);
710 0 : return true;
711 0 : }
712 : //
713 : // Set the supplied CFCache for all internal FTMs if they do use CFCache mechanism
714 : //
715 0 : void MultiTermFTNew::setCFCache(casacore::CountedPtr<CFCache>& cfc, const casacore::Bool resetCFC)
716 : {
717 0 : for (unsigned int i=0;i<subftms_p.nelements();i++)
718 : {
719 0 : if (subftms_p[i]->isUsingCFCache())
720 0 : subftms_p[i]->setCFCache(cfc,resetCFC);
721 : }
722 0 : }
723 :
724 :
725 : //---------------------------------------------------------------------------------------------------
726 : //---------------------------------------------------------------------------------------------------
727 : //---------------------------------------------------------------------------------------------------
728 : } // end namespace refim
729 : } //# NAMESPACE CASA - END
730 :
|
