Line data Source code
1 : //# NewMultiTermFT.cc: Implementation of NewMultiTermFT 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 <casacore/images/Images/ImageInterface.h>
49 : #include <casacore/images/Images/SubImage.h>
50 :
51 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
52 :
53 : #include <synthesis/TransformMachines/NewMultiTermFT.h>
54 : #include <synthesis/TransformMachines/Utils.h>
55 :
56 : // This is the list of FTMachine types supported by NewMultiTermFT
57 : //#include <synthesis/TransformMachines/GridFT.h>
58 : #include <synthesis/TransformMachines/AWProjectFT.h>
59 : #include <synthesis/TransformMachines/AWProjectWBFT.h>
60 :
61 : using namespace casacore;
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 : #define PSOURCE false
65 : #define psource (IPosition(4,1536,1536,0,0))
66 :
67 : //----------------------------------------------------------------------
68 : //-------------------- constructors and descructors ----------------------
69 : //----------------------------------------------------------------------
70 0 : NewMultiTermFT::NewMultiTermFT(FTMachine *subftm, Int nterms, Double reffreq)
71 0 : :FTMachine(), nterms_p(nterms), donePSF_p(false),doingPSF_p(false),
72 0 : reffreq_p(reffreq), imweights_p(Matrix<Float>(0,0)), machineName_p("NewMultiTermFT"),
73 0 : pblimit_p((Float)1e-04), doWideBandPBCorrection_p(false), cacheDir_p("."),
74 0 : donePBTaylor_p(false), useConjBeams_p(true)
75 : {
76 0 : dbg_p=false;
77 :
78 0 : this->setBasePrivates(*subftm);
79 0 : canComputeResiduals_p = subftm->canComputeResiduals();
80 0 : if(dbg_p) cout << "MTFT :: constructor with subftm : "<< subftm->name() << ". It can compute residuals : " << canComputeResiduals_p << endl;
81 :
82 0 : subftms_p.resize(2*nterms_p-1);
83 0 : for(uInt termindex=0;termindex<2*nterms_p-1;termindex++)
84 : {
85 0 : if(dbg_p) cout << "Creating new FTM of type : " << subftm->name() << endl;
86 0 : subftms_p[termindex] = getNewFTM(subftm);
87 0 : subftms_p[termindex]->setMiscInfo(termindex);
88 : }
89 0 : if(dbg_p)
90 : {
91 0 : cout << "Checking FTtypes at the end of MTFT constructor" << endl;
92 0 : printFTTypes();
93 : }
94 :
95 0 : sumwt_p=0.0;
96 :
97 : /// Make empty lists
98 0 : sensitivitymaps_p.resize(0);
99 0 : sumweights_p.resize(0);
100 0 : pbcoeffs_p.resize(0);
101 :
102 0 : if(subftm->name()=="AWProjectWBFT")
103 : {
104 : // doWideBandPBCorrection_p = ((AWProjectWBFT*)subftm)->getDOPBCorrection();
105 : // useConjBeams_p = ((AWProjectWBFT*)subftm)->getConjBeams();
106 0 : setDOPBCorrection(((AWProjectWBFT*)subftm)->getDOPBCorrection());
107 0 : setConjBeams(((AWProjectWBFT*)subftm)->getConjBeams());
108 :
109 0 : pblimit_p = ((AWProjectWBFT*)subftm)->getPBLimit();
110 0 : cout << "dowidebandpbcor : " << doWideBandPBCorrection_p << " conjbeams : " << useConjBeams_p << endl;
111 : }
112 :
113 0 : if(dbg_p) cout << "Running MTFT with doWideBandPBCorrection : " << doWideBandPBCorrection_p
114 0 : << " and pblimit : " << pblimit_p << endl;
115 :
116 0 : pblimit_p = 1e-07;
117 :
118 0 : cacheDir_p = subftm->getCacheDir();
119 :
120 0 : time_get=0.0;
121 0 : time_put=0.0;
122 0 : time_res=0.0;
123 0 : }
124 :
125 : //----------------------------------------------------------------------
126 : // Construct from the input state record
127 0 : NewMultiTermFT::NewMultiTermFT(const RecordInterface& stateRec)
128 0 : : FTMachine()
129 : {
130 0 : String error;
131 0 : if (!fromRecord(error, stateRec)) {
132 0 : throw (AipsError("Failed to create gridder: " + error));
133 : };
134 0 : }
135 :
136 : //----------------------------------------------------------------------
137 : // Copy constructor
138 0 : NewMultiTermFT::NewMultiTermFT(const NewMultiTermFT& other) : FTMachine(), machineName_p("NewMultiTermFT")
139 : {
140 0 : operator=(other);
141 0 : }
142 :
143 0 : NewMultiTermFT& NewMultiTermFT::operator=(const NewMultiTermFT& other)
144 : {
145 :
146 0 : dbg_p = other.dbg_p;
147 0 : if(dbg_p) cout << "In MTFT operator= " << endl;
148 :
149 0 : if(this!=&other)
150 : {
151 0 : FTMachine::operator=(other);
152 :
153 : // Copy local privates
154 0 : machineName_p = other.machineName_p;
155 0 : nterms_p = other.nterms_p;
156 0 : reffreq_p = other.reffreq_p;
157 0 : sumwt_p = other.sumwt_p;
158 0 : donePSF_p=other.donePSF_p;
159 0 : doingPSF_p=other.doingPSF_p;
160 0 : doWideBandPBCorrection_p = other.doWideBandPBCorrection_p;
161 0 : pblimit_p = other.pblimit_p;
162 0 : cacheDir_p = other.cacheDir_p;
163 0 : donePBTaylor_p = other.donePBTaylor_p;
164 0 : useConjBeams_p = other.useConjBeams_p;
165 :
166 : // Make the list of subftms
167 0 : subftms_p.resize(other.subftms_p.nelements());
168 0 : for (uInt termindex=0;termindex<other.subftms_p.nelements();termindex++)
169 : {
170 0 : subftms_p[termindex] = getNewFTM( &(*(other.subftms_p[termindex])) );
171 0 : subftms_p[termindex]->setMiscInfo(termindex);
172 : }
173 : // subftms_p[termindex] = getNewFTM( &(*(other.subftms_p[termindex])) );
174 :
175 : // Just checking....
176 0 : AlwaysAssert( subftms_p.nelements()>0 , AipsError );
177 :
178 : // Check if the sub ftm type can calculate its own residuals....
179 0 : canComputeResiduals_p = subftms_p[0]->canComputeResiduals();
180 :
181 : }
182 :
183 0 : if(dbg_p)
184 : {
185 0 : cout << "Checking FTtypes at the end of MTFT operator=" << endl;
186 0 : printFTTypes();
187 : }
188 :
189 0 : return *this;
190 :
191 : }
192 :
193 0 : FTMachine* NewMultiTermFT::getNewFTM(const FTMachine *ftm)
194 : {
195 : // if(dbg_p) cout << "MTFT::getNewFTM " << endl;
196 : /*
197 : if(ftm->name()=="GridFT")
198 : {
199 : return new GridFT(static_cast<const GridFT&>(*ftm));
200 : }
201 0 : else */ if(ftm->name()=="AWProjectWBFT")
202 0 : { return new AWProjectWBFT(static_cast<const AWProjectWBFT&>(*ftm)); }
203 : else
204 0 : {throw(AipsError("FTMachine "+ftm->name()+" is not supported with MS-MFS")); }
205 :
206 : return NULL;
207 : }
208 :
209 : //----------------------------------------------------------------------
210 0 : NewMultiTermFT::~NewMultiTermFT()
211 : {
212 0 : if(dbg_p) cout << "MTFT :: destructor - assumes automatic deletion of subftm " << endl;
213 0 : }
214 :
215 :
216 : //---------------------------------------------------------------------------------------------------
217 : //------------ Multi-Term Specific Functions --------------------------------------------------------
218 : //---------------------------------------------------------------------------------------------------
219 :
220 : // Multiply the imaging weights by Taylor functions - in place
221 : // This function MUST be called in ascending Taylor-term order
222 : // NOTE : Add checks to ensure this.
223 0 : Bool NewMultiTermFT::modifyVisWeights(VisBuffer& vb,uInt thisterm)
224 : {
225 : {
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.frequency());
233 :
234 0 : for (Int row=0; row<vb.nRow(); row++)
235 0 : for (Int chn=0; chn<vb.nChannel(); chn++)
236 : {
237 0 : freq = selfreqlist(IPosition(1,chn));
238 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
239 0 : (vb.imagingWeight())(chn,row) *= pow( mulfactor,(Int)thisterm);
240 : // sumwt_p += (vb.imagingWeight())(chn,row);
241 : }
242 0 : }
243 : /* // For debugging.
244 : else
245 : {
246 : for (Int row=0; row<vb.nRow(); row++)
247 : for (Int chn=0; chn<vb.nChannel(); chn++)
248 : {
249 : sumwt_p += (vb.imagingWeight())(chn,row);
250 : }
251 : }
252 : */
253 0 : return true;
254 : }
255 :
256 : // Reset the imaging weights back to their original values
257 : // to be called just after "put"
258 0 : void NewMultiTermFT::restoreImagingWeights(VisBuffer &vb)
259 : {
260 0 : AlwaysAssert( imweights_p.shape() == vb.imagingWeight().shape() ,AipsError);
261 0 : vb.imagingWeight() = imweights_p;
262 0 : }
263 :
264 :
265 : // Multiply the model visibilities by the Taylor functions - in place.
266 0 : Bool NewMultiTermFT::modifyModelVis(VisBuffer& vb, uInt thisterm)
267 : {
268 0 : Float freq=0.0,mulfactor=1.0;
269 0 : Vector<Double> selfreqlist(vb.frequency());
270 :
271 : // DComplex modcount=0.0;
272 :
273 0 : for (uInt pol=0; pol< uInt((vb.modelVisCube()).shape()[0]); pol++)
274 0 : for (uInt chn=0; chn< uInt(vb.nChannel()); chn++)
275 0 : for (uInt row=0; row< uInt(vb.nRow()); row++)
276 : {
277 : // modcount += ( vb.modelVisCube())(pol,chn,row);
278 0 : freq = selfreqlist(IPosition(1,chn));
279 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
280 0 : (vb.modelVisCube())(pol,chn,row) *= pow(mulfactor, (Int) thisterm);
281 : }
282 :
283 : // cout << "field : " << vb.fieldId() << " spw : "
284 : // << vb.spectralWindow() << " --- predicted model before taylor wt mult :"
285 : // << thisterm << " sumvis : " << modcount << endl;
286 :
287 0 : return true;
288 0 : }
289 :
290 :
291 : //---------------------------------------------------------------------------------------------------
292 : //---------------------- Prediction and De-gridding -----------------------------------
293 : //---------------------------------------------------------------------------------------------------
294 :
295 0 : void NewMultiTermFT::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)
296 : {
297 0 : if(dbg_p) cout << "MTFT::initializeToVis " << endl;
298 0 : AlwaysAssert(compImageVec.nelements()==nterms_p, AipsError);
299 0 : AlwaysAssert(modelImageVec.nelements()==nterms_p, AipsError);
300 0 : AlwaysAssert(weightImageVec.nelements()==nterms_p, AipsError);
301 0 : AlwaysAssert(fluxScaleVec.nelements()==nterms_p, AipsError);
302 0 : AlwaysAssert(weightsVec.nelements()==nterms_p, AipsError);
303 0 : Matrix<Float> tempWts;
304 :
305 : // Use doWideBandPBCorrection_p to trigger whether or not to do a wideband PB correction before prediction, for the first go ( i.e. simulation )
306 :
307 :
308 : // This is to make sure weight images and avgPBs are ready.
309 : //AlwaysAssert( donePSF_p == true && donePBTaylor_p == true, AipsError )
310 :
311 : if(PSOURCE) cout << "------ model, before de-gridding norm : " << modelImageVec[0]->getAt(psource) << "," << modelImageVec[1]->getAt(psource) << endl;
312 :
313 :
314 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
315 : {
316 : // Make the sensitivity Image if applicable
317 0 : subftms_p[taylor]->findConvFunction(*(compImageVec[taylor]), vb);
318 : // Get the sensitivity Image
319 0 : tempWts.resize();
320 0 : subftms_p[taylor]->getWeightImage(*(weightImageVec[taylor]), tempWts);
321 : }
322 :
323 : /////////////////////////////
324 : //UU normAvgPBs(weightImageVec);
325 :
326 : // Pre-prediction correction of the model, to rid it of the Primary Beam.
327 : if(PSOURCE) cout << "Divide the models by the PB before prediction" << endl;
328 0 : if( useConjBeams_p == true )
329 : {
330 : // Model contains only avgPB scaling, but no PB frequency dependence
331 : // Divide all terms of the model image by the sensitivity image only (from Taylor0)
332 0 : if (doWideBandPBCorrection_p)
333 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
334 : {
335 : // Divide by PB ////// PBWeight
336 : //normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)1);// normtype 1 divides by weightImageVec and ignores wegithsVec
337 : // Divide by sqrt(PB) ////// PBSQWeight
338 0 : normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)4);
339 : }
340 : }
341 : else
342 : {
343 0 : cout << "NewMultiTermFT::inittoVis --> Prediction normalization of model will be wrong....... " << endl;
344 : // The model contains avgPB scaling and twice the frequency-dependence of the PB in it.
345 : // Multiply all terms of the model by pbcoeffs_0. Divide TWICE by the wideband PB.
346 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
347 : {
348 0 : normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)3); // normtype 3 multiplies the model image with the pb
349 : }
350 : // Connect pbcoeffs_p to fluxScaleVec. This is where PB Taylor coefficients were put in by the 'iftm'.
351 0 : if(pbcoeffs_p.nelements() != nterms_p) pbcoeffs_p.resize(nterms_p);
352 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
353 : {
354 0 : pbcoeffs_p[taylor] = &(*(fluxScaleVec[taylor]));
355 : }
356 :
357 0 : applyWideBandPB( String("divide") , modelImageVec );
358 0 : applyWideBandPB( String("divide") , modelImageVec );
359 : }
360 :
361 : // for(uInt taylor=0;taylor<nterms_p;taylor++)
362 : // storeAsImg("flattenedModel_"+String::toString(taylor) , *(modelImageVec[taylor]) );
363 :
364 : if(PSOURCE) cout << "------ model, after de-gridding norm : " << modelImageVec[0]->getAt(psource) << "," << modelImageVec[1]->getAt(psource) << endl;
365 :
366 :
367 :
368 : // Convert Stokes planes to correlation planes..
369 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
370 : {
371 0 : stokesToCorrelation( *(modelImageVec[taylor]) , *(compImageVec[taylor]) );
372 : }
373 :
374 : // Call initializeToVis for all sub ftms..
375 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
376 : {
377 0 : subftms_p[taylor]->initializeToVis(*(compImageVec[taylor]),vb);
378 : }
379 0 : time_get=0.0;
380 :
381 : /// Multiply the model with the avgPB again, so that it's ready for the minor cycle incremental accumulation
382 : if(PSOURCE) cout << "Multiplying the models by the weightimage to reset it to flat-noise for the minor cycle" << endl;
383 :
384 0 : if( useConjBeams_p == true )
385 : {
386 :
387 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
388 : {
389 : //Mulitply by PB ///// PBWeight
390 : //normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)3); // normtype 3 multiplies the model image with the pb
391 : //Mulitply by sqrt(PB) //// PBSQWeight
392 0 : normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)5);
393 :
394 : }
395 : }
396 : else
397 : {
398 :
399 0 : applyWideBandPB( String("multiply") , modelImageVec );
400 0 : applyWideBandPB( String("multiply") , modelImageVec );
401 :
402 : // Divide by the avg PB 0
403 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
404 : {
405 0 : normalizeImage( *(modelImageVec[taylor]) , weightsVec[0], *(weightImageVec[0]) , false, (Float)pblimit_p, (Int)1);
406 : }
407 :
408 :
409 : }
410 :
411 : if(PSOURCE) cout << "------ model, gone back : " << modelImageVec[0]->getAt(psource) << "," << modelImageVec[1]->getAt(psource) << endl;
412 :
413 :
414 0 : }// end of initializeToVis
415 :
416 :
417 :
418 0 : void NewMultiTermFT::get(VisBuffer& vb, Int row)
419 : {
420 :
421 : // De-grid the model for the zeroth order Taylor term
422 0 : subftms_p[0]->get(vb,row);
423 : // Save the model visibilities in a local cube
424 0 : modviscube_p.assign( vb.modelVisCube() );
425 :
426 0 : for(uInt tix=1;tix<nterms_p;tix++) // Only nterms.... not 2nterms-1
427 : {
428 : // Reset the model visibilities to zero
429 0 : vb.setModelVisCube(Complex(0.0,0.0));
430 : // De-grid the model onto the modelviscube (other Taylor terms)
431 0 : subftms_p[tix]->get(vb,row);
432 : // Multiply visibilities by taylor-weights
433 0 : modifyModelVis(vb,tix);
434 : // Accumulate model visibilities across Taylor terms
435 0 : modviscube_p += vb.modelVisCube();
436 : }
437 : // Set the vb.modelviscube to what has been accumulated
438 0 : vb.setModelVisCube(modviscube_p);
439 0 : }
440 :
441 0 : void NewMultiTermFT::finalizeToVis()
442 : {
443 0 : if(dbg_p) cout << "MTFT::finalizeToVis for " << nterms_p << " terms "<<endl;
444 0 : AlwaysAssert(subftms_p.nelements() >= nterms_p , AipsError);
445 0 : for(uInt taylor=0;taylor<nterms_p;taylor++) subftms_p[taylor]->finalizeToVis();
446 0 : }
447 :
448 :
449 : //---------------------------------------------------------------------------------------------------
450 : //---------------------- Calculate Residual Visibilities -------------------------------
451 : //---------------------------------------------------------------------------------------------------
452 0 : void NewMultiTermFT::ComputeResiduals(VisBuffer &vb, Bool useCorrected)
453 : {
454 :
455 0 : if(subftms_p[0]->canComputeResiduals()) subftms_p[0]->ComputeResiduals(vb,useCorrected);
456 0 : else throw(AipsError("MultiTerm::ComputeResiduals : subftm of NewMultiTermFT cannot compute its own residuals !"));
457 :
458 0 : }
459 :
460 : //---------------------------------------------------------------------------------------------------
461 : //---------------------- Gridding --------------------------------------------------------------
462 : //---------------------------------------------------------------------------------------------------
463 :
464 0 : void NewMultiTermFT::initializeToSky(Block<CountedPtr<ImageInterface<Complex> > > & compImageVec, Block<Matrix<Float> >& weightsVec, const VisBuffer& vb, const Bool dopsf)
465 : {
466 0 : if(dbg_p) cout << "MTFT::initializeToSky with dopsf : " << dopsf << endl;
467 :
468 : // If PSF is already done, don't ask again !
469 0 : AlwaysAssert( !(donePSF_p && dopsf) , AipsError );
470 :
471 : // The PSF needs to be the first thing made (because of weight images)
472 0 : AlwaysAssert( !(dopsf==false && donePSF_p==false) , AipsError);
473 :
474 0 : doingPSF_p=false;
475 0 : if(donePSF_p==true) // TODO : Check if we can clean up the extra ftmachines, or if the weightImageVecs from extra ones are still used later
476 : {
477 0 : if( subftms_p.nelements() != nterms_p )
478 : {
479 0 : subftms_p.resize( nterms_p ,true);
480 0 : if(dbg_p) cout << "MTFT::initializeToSky : resizing to " << nterms_p << " terms" << endl;
481 : }
482 : }
483 0 : uInt gridnterms=nterms_p;
484 0 : if(dopsf){ gridnterms=2*nterms_p-1; }
485 :
486 0 : AlwaysAssert(gridnterms <= subftms_p.nelements() , AipsError);
487 0 : AlwaysAssert(compImageVec.nelements()==gridnterms, AipsError);
488 0 : AlwaysAssert(weightsVec.nelements()==gridnterms, AipsError);
489 :
490 :
491 0 : if(dbg_p) cout << "MTFT : Calling subft initializeToSky for " << gridnterms << " terms " << endl;
492 0 : for(uInt taylor=0;taylor<gridnterms;taylor++)
493 : {
494 0 : subftms_p[taylor]->initializeToSky(*(compImageVec[taylor]),weightsVec[taylor],vb);
495 : }
496 :
497 0 : }// end of initializeToSky
498 :
499 :
500 :
501 0 : void NewMultiTermFT::put(VisBuffer& vb, Int row, Bool dopsf, FTMachine::Type type)
502 : {
503 :
504 0 : subftms_p[0]->put(vb,row,dopsf,type);
505 :
506 0 : Int gridnterms=nterms_p;
507 0 : if(dopsf==true && donePSF_p==false)
508 : {
509 0 : gridnterms=2*nterms_p-1;
510 0 : doingPSF_p=true;
511 : }
512 :
513 : //cout << " Calling put for " << gridnterms << " terms, nelements : " << subftms_p.nelements() << " and dopsf " << dopsf << endl;
514 :
515 0 : for(Int tix=1;tix<gridnterms;tix++)
516 : {
517 0 : modifyVisWeights(vb,tix);
518 0 : subftms_p[tix]->put(vb,row,dopsf,type);
519 0 : restoreImagingWeights(vb);
520 : }
521 :
522 :
523 0 : }// end of put
524 :
525 : //----------------------------------------------------------------------
526 0 : void NewMultiTermFT::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*/)
527 : {
528 0 : uInt gridnterms=nterms_p;
529 0 : if(dopsf==true) { gridnterms=2*nterms_p-1; }
530 :
531 0 : if(dbg_p) cout << "MTFT : finalizeToSky for " << gridnterms << " terms and dopsf : " << dopsf << endl;
532 :
533 0 : AlwaysAssert(gridnterms <= subftms_p.nelements() , AipsError);
534 0 : AlwaysAssert(compImageVec.nelements()==gridnterms, AipsError);
535 0 : AlwaysAssert(resImageVec.nelements()==gridnterms, AipsError);
536 0 : AlwaysAssert(weightImageVec.nelements()==gridnterms, AipsError);
537 0 : AlwaysAssert(fluxScaleVec.nelements()==gridnterms, AipsError);
538 0 : AlwaysAssert(weightsVec.nelements()==gridnterms, AipsError);
539 :
540 : // Collect images and weights from all FTMs
541 0 : for(uInt taylor=0;taylor<gridnterms;taylor++)
542 : {
543 : // Call finalizeToSky for subftm
544 0 : subftms_p[taylor]->finalizeToSky();
545 : // Get the gridded image
546 0 : (*(compImageVec[taylor])).copyData(subftms_p[taylor]->getImage(weightsVec[taylor],false));
547 : // Convert to Stokes planes
548 0 : correlationToStokes((*(compImageVec[taylor])) , (*(resImageVec[taylor])) , dopsf );
549 : // Get the weight image.
550 0 : subftms_p[taylor]->getWeightImage(*(weightImageVec[taylor]), weightsVec[taylor]);
551 : }// end for taylor
552 :
553 : //Normalize the weight images by the peak of the zero'th order weightImage.
554 : //UU normAvgPBs(weightImageVec);
555 :
556 : // Norm by sumwt for PSFs and Residuals.
557 0 : for(uInt taylor=0;taylor<gridnterms;taylor++)
558 : {
559 0 : normalizeImage( *(resImageVec[taylor]) , weightsVec[0] , *(weightImageVec[0]) ,
560 0 : dopsf , (Float)pblimit_p, (Int)0);
561 :
562 : }// end for taylor
563 :
564 :
565 : if(PSOURCE) cout << "------ residual, before normalization : " << resImageVec[0]->getAt(psource) << "," << resImageVec[1]->getAt(psource) << endl;
566 :
567 : // Use sumwts to make Hessian, invert, apply to weight images to fill in pbcoeffs
568 : // TODO : Clean up this list of state variables !
569 0 : if(dopsf==true && doingPSF_p==true && donePSF_p==false)
570 : {
571 :
572 0 : if( donePBTaylor_p == false )
573 : {
574 :
575 0 : cout << "MTFT::finalizeToSky for PSF and Weights : Calculating PB coefficients" << endl;
576 :
577 : // Gather normalized sumweights for the Hessian matrix.
578 0 : sumweights_p.resize(gridnterms);
579 0 : for(uInt taylor=0;taylor<gridnterms;taylor++)
580 0 : { sumweights_p[taylor] = weightsVec[taylor]/weightsVec[0]; }
581 :
582 : // Connect pbcoeffs_p to fluxScaleVec. This is where PB Taylor coefficients will end up.
583 0 : if(pbcoeffs_p.nelements() != nterms_p) pbcoeffs_p.resize(nterms_p);
584 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
585 : {
586 0 : pbcoeffs_p[taylor] = &(*(fluxScaleVec[taylor]));
587 : }
588 :
589 : // Do the calculation
590 : // UUU Taken out ( Sep2013 ) because sqrt(PB) is being accumulated and this is not OK
591 : // for calculating Coeff PBs. Need to eventually.....
592 : // ......store PBSQ, Calc coeffs for PB squared, Do a polynomial square root.
593 : // Until then, 'coeffPB' cannot be used for post-deconv corrections.
594 : // Use sensitivityPB, and remember when it's a square and when not.
595 : //
596 : //calculateTaylorPBs(weightImageVec);
597 :
598 0 : donePBTaylor_p = true;
599 :
600 :
601 : }
602 :
603 : // Show the value....
604 : if(PSOURCE) cout << " PB 0 : " << pbcoeffs_p[0]->getAt(psource) << " PB 1 : " << pbcoeffs_p[1]->getAt(psource) << endl;
605 :
606 :
607 : }// if dopsf
608 :
609 :
610 : /*
611 : // Normalize all by the Taylor0 weights
612 : AlwaysAssert( sensitivitymaps_p.nelements() > 0 , AipsError );
613 : for(uInt taylor=0;taylor<gridnterms;taylor++)
614 : {
615 : normalizeImage( *(resImageVec[taylor]) , weightsVec[0] , *(weightImageVec[0]) , dopsf , (Float)pblimit_p, (Int)1); //// use locally-normalized avgPB0.
616 : }// end for taylor
617 : */
618 :
619 : if(PSOURCE) cout << "------ residual, after normalization : " << resImageVec[0]->getAt(psource) << "," << resImageVec[1]->getAt(psource) << endl;
620 :
621 :
622 :
623 0 : if(doingPSF_p==true)
624 0 : {doingPSF_p=false; donePSF_p=true; if(dbg_p) cout << "Setting donePSF to true" << endl;}
625 :
626 :
627 0 : }//end of finalizeToSky
628 :
629 : //---------------------------------------------------------------------------------------------------
630 : //----------------------------- Obtain Images -----------------------------------------------------
631 : //---------------------------------------------------------------------------------------------------
632 : //----------------------------------------------------------------------
633 0 : void NewMultiTermFT::makeImage(FTMachine::Type type, VisSet& vs,
634 : ImageInterface<Complex>& theImage, Matrix<Float>& weight)
635 : {
636 0 : if(dbg_p) cout << "MTFT :: makeImage for taylor 0 only "<< endl;
637 0 : subftms_p[0]->makeImage(type, vs, theImage, weight);
638 0 : }
639 :
640 : // Connect sensitivitymaps_p to weightImages, and normalize all by peak of 0th one.
641 0 : void NewMultiTermFT::normAvgPBs(PtrBlock<SubImage<Float> *>& weightImageVec)
642 : {
643 0 : AlwaysAssert( weightImageVec.nelements() >= nterms_p , AipsError );
644 0 : Matrix<Float> tempMat;
645 0 : Array<Float> tempArr;
646 0 : ( *(weightImageVec[0]) ).get(tempArr,true);
647 0 : tempMat.reference(tempArr);
648 0 : Float maxval = max(tempMat);
649 :
650 : /*
651 : if(sensitivitymaps_p.nelements()==0 && weightImageVec.nelements() == 2*nterms_p-1)
652 : {
653 : sensitivitymaps_p.resize(2*nterms_p-1);
654 : for(uInt taylor=0;taylor<2*nterms_p-1;taylor++)
655 : {
656 : Float rmaxval = maxval;
657 : cout << "Normalizing pb : " << taylor << " by peak of zeroth : " << rmaxval << endl;
658 :
659 : sensitivitymaps_p[taylor] = new PagedImage<Float>( (weightImageVec[taylor])->shape() , (weightImageVec[taylor])->coordinates() , cacheDir_p+"/sensitivityPB_"+String::toString(taylor) );
660 : sensitivitymaps_p[taylor]->copyData( (LatticeExpr<Float>) ( (*(weightImageVec[taylor]))/rmaxval ) );
661 : }
662 : }
663 : */
664 :
665 : // Normalize weightImageVecs in-place
666 0 : for(uInt taylor=0;taylor<weightImageVec.nelements();taylor++)
667 : {
668 : //cout << "MTFT :: Normalizing wtimg : " << taylor << " by peak of zeroth : " << maxval << endl;
669 0 : weightImageVec[taylor]->copyData( (LatticeExpr<Float>) ( (*(weightImageVec[taylor]))/maxval ) );
670 : }
671 :
672 0 : }
673 :
674 :
675 : //---------------------------------------------------------------------------------------------------
676 : //------------------------ To / From Records ---------------------------------------------------------
677 : //---------------------------------------------------------------------------------------------------
678 0 : Bool NewMultiTermFT::toRecord(String& error, RecordInterface& outRec, Bool withImage, const String diskimage)
679 : {
680 0 : if(dbg_p) cout << "MTFT :: toRecord for " << nterms_p << endl;
681 0 : Bool retval = true;
682 :
683 0 : for(uInt tix=0;tix<nterms_p;tix++)
684 : {
685 0 : Record subFTContainer;
686 0 : String elimage="";
687 0 : if(diskimage != ""){
688 0 : elimage=diskimage+String("_term_")+String::toString(tix);
689 : }
690 0 : subftms_p[tix]->toRecord(error, subFTContainer,withImage, elimage);
691 0 : outRec.defineRecord("subftm_"+String::toString(tix),subFTContainer);
692 0 : }
693 :
694 : // Record subFTContainer;
695 : // subftm_p->toRecord(error, subFTContainer,withImage);
696 : // outRec.defineRecord("subftm",subFTContainer);
697 :
698 0 : outRec.define("nterms",nterms_p);
699 0 : outRec.define("reffreq",reffreq_p);
700 :
701 0 : return retval;
702 : }
703 :
704 : //---------------------------------------------------------------------------------------------------
705 0 : Bool NewMultiTermFT::fromRecord(String& error, const RecordInterface& inRec)
706 : {
707 0 : if(dbg_p) cout << "MTFT :: fromRecord "<< endl;
708 0 : Bool retval = true;
709 :
710 : // Record subFTMRec=inRec.asRecord("subftm");
711 : // retval = (retval || subftm_p->fromRecord(error, subFTMRec));
712 :
713 0 : inRec.get("nterms",nterms_p);
714 0 : inRec.get("reffreq",reffreq_p);
715 :
716 0 : subftms_p.resize(nterms_p);
717 0 : for(uInt tix=0;tix<nterms_p;tix++)
718 : {
719 0 : Record subFTMRec=inRec.asRecord("subftm_"+String::toString(tix));
720 0 : retval = (retval || subftms_p[tix]->fromRecord(error, subFTMRec));
721 0 : }
722 :
723 :
724 0 : return retval;
725 : }
726 : //---------------------------------------------------------------------------------------------------
727 :
728 0 : Bool NewMultiTermFT::storeAsImg(String fileName, ImageInterface<Float>& theImg)
729 : {
730 0 : PagedImage<Float> tmp(theImg.shape(), theImg.coordinates(), fileName);
731 0 : LatticeExpr<Float> le(theImg);
732 0 : tmp.copyData(le);
733 0 : return true;
734 0 : }
735 :
736 :
737 : //---------------------------------------------------------------------------------------------------
738 : // Make pixel-by-pixel matrices from pbcoeffs, (invert), and multiply with the given vector (in place)
739 0 : void NewMultiTermFT::applyWideBandPB(String action, PtrBlock<SubImage<Float> *> &imageVec)
740 : {
741 : //readAvgPBs(); // Should read only once.
742 0 : AlwaysAssert( pbcoeffs_p.nelements()==nterms_p , AipsError );
743 0 : AlwaysAssert( imageVec.nelements()==nterms_p , AipsError );
744 :
745 0 : Int nX=imageVec[0]->shape()(0);
746 0 : Int nY=imageVec[0]->shape()(1);
747 0 : Int npol=imageVec[0]->shape()(2);
748 0 : Int nchan=imageVec[0]->shape()(3);
749 :
750 0 : AlwaysAssert(nchan==1,AipsError);
751 0 : AlwaysAssert(npol==1,AipsError);
752 :
753 0 : Double deter=0.0;
754 0 : Matrix<Double> mat( IPosition(2,nterms_p,nterms_p) );
755 0 : Matrix<Double> invmat( IPosition(2,nterms_p,nterms_p) );
756 :
757 0 : mat.set(0.0);
758 0 : invmat.set(0.0);
759 :
760 : // Go over all pixels for which coeffPB_0 is above pblimit_p
761 0 : Vector<Float> pbvec(nterms_p), invec(nterms_p), outvec(nterms_p);
762 0 : IPosition tip(4,0,0,0,0);
763 0 : for(tip[0]=0;tip[0]<nX;tip[0]++)
764 : {
765 0 : for(tip[1]=0;tip[1]<nY;tip[1]++)
766 : {
767 : // Normalize only if pb > limit
768 0 : if(fabs(pbcoeffs_p[0]->getAt(tip)) > pblimit_p )
769 : {
770 : // Fill in the single-pixel Hessian, and RHS and LHS vectors
771 0 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
772 : {
773 0 : for(uInt taylor2=0;taylor2<=taylor1;taylor2++)
774 : {
775 0 : mat(taylor1,taylor2) = (pbcoeffs_p[ taylor1 - taylor2 ])->getAt(tip);
776 : }
777 0 : invec[taylor1] = imageVec[taylor1]->getAt(tip);
778 0 : outvec[taylor1]=0.0;
779 : }
780 :
781 : // Invert hess_p into invhess_p;
782 : try
783 : {
784 : //invertSymPosDef((invmat),deter,(mat));
785 0 : invert((invmat),deter,(mat));
786 : }
787 0 : catch(AipsError &x)
788 : {
789 0 : cout << "The non-invertible Matrix is : " << (mat) << endl;
790 0 : throw( AipsError("Cannot Invert matrix for PB application: " + x.getMesg() ) );
791 0 : }
792 :
793 : if(PSOURCE)
794 : {
795 : if( tip[0] == psource[0] && tip[1] == psource[1] )
796 : {
797 : cout << "PB mat : " << mat << endl;
798 : cout << "InvPB mat : " << invmat << endl;
799 : }
800 : }
801 :
802 : // Multiply invhess_p by RHS and fill in LHS vector
803 0 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
804 : {
805 0 : outvec[taylor1]=0.0;
806 0 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
807 : {
808 0 : if( action == String("divide") )
809 : {
810 0 : outvec[taylor1] += invmat(taylor1,taylor2) * (invec[taylor2]) ;
811 : }
812 : else
813 : {
814 0 : outvec[taylor1] += mat(taylor1,taylor2) * (invec[taylor2]) ;
815 : }
816 : } // for taylor2
817 : }// for taylor1
818 :
819 : /*
820 : if(tip==psource)
821 : {
822 : cout << "------ normalizeWideBandPB2 : before : " << imageVec[0]->getAt(tip) << "," << imageVec[1]->getAt(tip) << " after : " << outvec[0] << "," << outvec[1] << " hess : " << hess_p << endl;
823 : }
824 : */
825 :
826 : // Put the solution into the residual images.
827 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
828 0 : imageVec[taylor]->putAt(outvec[taylor], tip);
829 :
830 : }// if larger than pblimit
831 : else // if smaller than pblimit
832 : {
833 : // Put 0.0 into the residual images for un-normalizable parts
834 0 : for(uInt taylor=0;taylor<nterms_p;taylor++)
835 0 : imageVec[taylor]->putAt((Float)0.0, tip);
836 : }
837 :
838 : }//end of for y
839 : }// end of for x
840 :
841 0 : }// end of applyWideBandPB
842 :
843 :
844 : /*
845 :
846 : //---------------------------------------------------------------------------------------------------
847 : // Make pixel-by-pixel matrices from pbcoeffs, (invert), and multiply with the given vector (in place)
848 : void NewMultiTermFT::applyWideBandPB(String action, PtrBlock<SubImage<Float> *> &imageVec)
849 : {
850 : //readAvgPBs(); // Should read only once.
851 : AlwaysAssert( sensitivitymaps_p.nelements()==2*nterms_p-1 , AipsError );
852 : AlwaysAssert( imageVec.nelements()==nterms_p , AipsError );
853 :
854 : Int nX=imageVec[0]->shape()(0);
855 : Int nY=sensitivitymaps_p[0]->shape()(1);
856 : Int npol=sensitivitymaps_p[0]->shape()(2);
857 : Int nchan=sensitivitymaps_p[0]->shape()(3);
858 :
859 : AlwaysAssert(nchan==1,AipsError);
860 : AlwaysAssert(npol==1,AipsError);
861 :
862 : Double deter=0.0;
863 : hess_p.resize(IPosition(2,nterms_p,nterms_p) );
864 : invhess_p.resize(IPosition(2,nterms_p,nterms_p) );
865 :
866 : // IPosition psource(4,nX/2+22,nY/2,0,0);
867 :
868 : // cout << "Source Position : " << psource << endl;
869 :
870 : // Go over all pixels for which avgPB_tt0 is above pblimit_p
871 : Vector<Float> pbvec(nterms_p), resvec(nterms_p), normedvec(nterms_p);
872 : IPosition tip(4,0,0,0,0);
873 : for(tip[0]=0;tip[0]<nX;tip[0]++)
874 : {
875 : for(tip[1]=0;tip[1]<nY;tip[1]++)
876 : {
877 : // Normalize only if pb > limit
878 : if(1) // fabs(sensitivitymaps_p[0]->getAt(tip)) > pblimit_p )
879 : {
880 :
881 : // Fill in the single-pixel Hessian, and RHS and LHS vectors
882 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
883 : {
884 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
885 : {
886 : hess_p(taylor1,taylor2) = (sensitivitymaps_p[taylor1+taylor2])->getAt(tip);
887 : }
888 : resvec[taylor1] = imageVec[taylor1]->getAt(tip);
889 : normedvec[taylor1]=0.0;
890 : }
891 :
892 : // Invert hess_p into invhess_p;
893 : try
894 : {
895 : invertSymPosDef((invhess_p),deter,(hess_p));
896 : }
897 : catch(AipsError &x)
898 : {
899 : cout << "The non-invertible Hessian is : " << (hess_p) << endl;
900 : throw( AipsError("Cannot Invert matrix : " + x.getMesg() ) );
901 : }
902 :
903 : // Multiply invhess_p by RHS and fill in LHS vector
904 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
905 : {
906 : normedvec[taylor1]=0.0;
907 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
908 : {
909 : normedvec[taylor1] += invhess_p(taylor1,taylor2) * (resvec[taylor2]) ;
910 : } // for taylor2
911 : }// for taylor1
912 :
913 : // Put the solution into the residual images.
914 : for(uInt taylor=0;taylor<nterms_p;taylor++)
915 : imageVec[taylor]->putAt(normedvec[taylor], tip);
916 :
917 : }// if larger than pblimit
918 : else // if smaller than pblimit
919 : {
920 : // Put 0.0 into the residual images for un-normalizable parts
921 : for(uInt taylor=0;taylor<nterms_p;taylor++)
922 : imageVec[taylor]->putAt((Float)0.0, tip);
923 : }
924 :
925 : }//end of for y
926 : }// end of for x
927 :
928 : }// end of normalizeWideBandPB2
929 :
930 : */
931 :
932 : //---------------------------------------------------------------------------------------------------
933 : // Use sumwts to make a Hessian, invert it, apply to weight images, fill in pbcoeffs_p
934 : // This should get called only once, while making PSFs, when there are 2n-1 terms.
935 0 : void NewMultiTermFT::calculateTaylorPBs(PtrBlock<SubImage<Float> *> & weightImageVec)
936 : {
937 :
938 0 : AlwaysAssert( weightImageVec.nelements() == 2*nterms_p-1 , AipsError );
939 :
940 0 : for(uInt taylor=0;taylor<2*nterms_p-1;taylor++)
941 : {
942 0 : storeImg( cacheDir_p+"/sensitivityPB_"+String::toString(taylor) , *weightImageVec[taylor] );
943 : }
944 :
945 : /// Just read from weightimage, instead of sensitivity images...
946 :
947 0 : AlwaysAssert( pbcoeffs_p.nelements()==nterms_p, AipsError );
948 0 : AlwaysAssert( sumweights_p.nelements()==2*nterms_p-1, AipsError );
949 :
950 0 : Int npol=weightImageVec[0]->shape()(2);
951 0 : Int nchan=weightImageVec[0]->shape()(3);
952 :
953 0 : AlwaysAssert(nchan==1,AipsError);
954 0 : AlwaysAssert(npol==1,AipsError);
955 :
956 0 : Double deter=0.0;
957 0 : hess_p.resize(IPosition(2,nterms_p,nterms_p) );
958 0 : invhess_p.resize(IPosition(2,nterms_p,nterms_p) );
959 :
960 : // Fill hess_p from sumweights_p;
961 0 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
962 : {
963 0 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
964 : {
965 0 : hess_p(taylor1,taylor2) = (Double)( (sumweights_p[taylor1+taylor2])(0,0) );
966 : }// for taylor2
967 : }// for taylor1
968 0 : cout << "Hessian : " << hess_p << endl;
969 :
970 : // Invert hess_p into invhess_p;
971 : try
972 : {
973 0 : invertSymPosDef((invhess_p),deter,(hess_p));
974 : }
975 0 : catch(AipsError &x)
976 : {
977 0 : cout << "The non-invertible Hessian is : " << (hess_p) << endl;
978 0 : throw( AipsError("Cannot Invert matrix : " + x.getMesg() ) );
979 0 : }
980 0 : cout << "Inverse Hessian : " << invhess_p << endl;
981 :
982 0 : multiplyHMatrix( invhess_p, weightImageVec, pbcoeffs_p, String("/coeffPB_") );
983 :
984 : /*
985 : // Multiply invhess_p by sumweights_p and fill in pbcoeffs_p (Fig 7.3)
986 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
987 : {
988 : LatticeExprNode len(0.0);
989 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
990 : {
991 : len = len + LatticeExprNode(invhess_p(taylor1,taylor2) * (*(weightImageVec[taylor2])) );
992 : } // for taylor2
993 :
994 : pbcoeffs_p[taylor1]->copyData(LatticeExpr<Float> (iif( abs(len)<pblimit_p*pblimit_p , 0.0 , len )) );
995 : ///pbcoeffs_p[taylor1]->copyData(LatticeExpr<Float> (len) );
996 : storeImg(cacheDir_p+"/coeffPB_"+String::toString(taylor1) , *(pbcoeffs_p[taylor1]) );
997 : }// for taylor1
998 : */
999 :
1000 0 : return;
1001 : }
1002 :
1003 0 : void NewMultiTermFT::multiplyHMatrix( Matrix<Double> &hmat,
1004 : PtrBlock<SubImage<Float>* > &invec,
1005 : PtrBlock<SubImage<Float>* > &outvec,
1006 : String saveImagePrefix )
1007 : {
1008 0 : AlwaysAssert( hmat.shape().nelements()==2 &&
1009 : hmat.shape()[0] == nterms_p &&
1010 : hmat.shape()[1] == nterms_p , AipsError );
1011 0 : AlwaysAssert( invec.nelements() >= nterms_p , AipsError );
1012 0 : AlwaysAssert( outvec.nelements() >= nterms_p, AipsError );
1013 :
1014 0 : for(uInt taylor1=0;taylor1<nterms_p;taylor1++)
1015 : {
1016 0 : LatticeExprNode len(0.0);
1017 0 : for(uInt taylor2=0;taylor2<nterms_p;taylor2++)
1018 : {
1019 0 : len = len + LatticeExprNode(hmat(taylor1,taylor2) * (*(invec[taylor2])) );
1020 : } // for taylor2
1021 :
1022 0 : outvec[taylor1]->copyData(LatticeExpr<Float> (iif( abs(len)<pblimit_p*pblimit_p , 0.0 , len )) );
1023 :
1024 0 : if( saveImagePrefix.length() > 0 )
1025 : {
1026 0 : storeImg(cacheDir_p+saveImagePrefix+String::toString(taylor1) , *(outvec[taylor1]) );
1027 : }
1028 0 : }// for taylor1
1029 :
1030 0 : }// end of multiplyHMatrix
1031 :
1032 :
1033 : /*
1034 :
1035 : //---------------------------------------------------------------------------------------------------
1036 : // Apply inverse of single Hessian to residuals
1037 : void NewMultiTermFT::normalizeWideBandPB(PtrBlock<SubImage<Float> *> &resImageVec, PtrBlock<SubImage<Float> *>& scratchImageVec)
1038 : {
1039 : AlwaysAssert( scratchImageVec.nelements()==nterms_p , AipsError );
1040 : AlwaysAssert( pbcoeffs_p.nelements()==nterms_p , AipsError );
1041 : AlwaysAssert( resImageVec.nelements()>=nterms_p , AipsError );
1042 :
1043 :
1044 : cout << "MTFT::normWideBandPB : normalizing with " << nterms_p << " terms and pblim : " << pblimit_p << endl;
1045 : // Do a polynomial division of the PB from the residual images (using pbcoeffs_p) (Fig 7.4)
1046 :
1047 :
1048 : switch(nterms_p)
1049 : {
1050 : case 1:
1051 : {
1052 : LatticeExprNode deter1( (LatticeExpr<Float>) ( (*(pbcoeffs_p[0])) ) );
1053 : ///(resImageVec[0])->copyData( (LatticeExpr<Float>)( (*(resImageVec[0]))/(*(pbcoeffs_p[0])) ) );
1054 : (resImageVec[0])->copyData( (LatticeExpr<Float>)
1055 : ( iif( abs(deter1)<pblimit_p , 0.0 , (*(resImageVec[0]))/deter1 ) ) );
1056 : break;
1057 : }// end case 1
1058 :
1059 : case 2:
1060 : {
1061 : LatticeExprNode deter2( (LatticeExpr<Float>) ( (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[0])) ) );
1062 :
1063 : (scratchImageVec[0])->copyData( (LatticeExpr<Float>)(
1064 : (*(pbcoeffs_p[0]))*(*(resImageVec[0])) ) );
1065 : (scratchImageVec[1])->copyData( (LatticeExpr<Float>) (
1066 : ( ( -1.0* (*(pbcoeffs_p[1]))*(*(resImageVec[0])) ) +
1067 : ( (*(pbcoeffs_p[0]))*(*(resImageVec[1])) )
1068 : ) ) );
1069 :
1070 : for (uInt tay=0;tay<nterms_p;tay++)
1071 : {
1072 : (resImageVec[tay])->copyData( (LatticeExpr<Float>)
1073 : ( iif( abs(deter2)<pblimit_p*pblimit_p , 0.0 , (*(scratchImageVec[tay]))/deter2 ) ) );
1074 : }
1075 : break;
1076 : }// end case 2
1077 :
1078 : case 3:
1079 : {
1080 : LatticeExprNode deter3( (LatticeExpr<Float>) ( (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[0]) ) ) );
1081 : (scratchImageVec[0])->copyData( (LatticeExpr<Float>) (
1082 : (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[0])) * (*(resImageVec[0])) ) );
1083 : (scratchImageVec[1])->copyData( (LatticeExpr<Float>) ( (
1084 : ( -1.0 * (*(pbcoeffs_p[0]) )*(*(pbcoeffs_p[1]))*(*(resImageVec[0])) ) +
1085 : ( (*(pbcoeffs_p[0]))*(*(pbcoeffs_p[0]))*(*(resImageVec[1])) )
1086 : ) ) );
1087 : (scratchImageVec[2])->copyData((LatticeExpr<Float>) ( (
1088 : (
1089 : ( (*(pbcoeffs_p[1]))*(*(pbcoeffs_p[1])) - (*(pbcoeffs_p[0]))*(*(pbcoeffs_p[2]))
1090 : ) * (*(resImageVec[0]))
1091 : ) +
1092 : ( -1.0 * (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[1])) * (*(resImageVec[1])) ) +
1093 : ( (*(pbcoeffs_p[0])) * (*(pbcoeffs_p[0])) * (*(resImageVec[2])) )
1094 : ) ) );
1095 : for (uInt tay=0;tay<nterms_p;tay++)
1096 : {
1097 : (resImageVec[tay])->copyData( (LatticeExpr<Float>)
1098 : ( iif( abs(deter3)<pblimit_p*pblimit_p*pblimit_p , 0.0 , (*(scratchImageVec[tay]))/deter3 ) ) );
1099 : }
1100 : break;
1101 : }// end case 3
1102 :
1103 : default:
1104 : throw( AipsError("Cannot compute PB Coefficients for more than nterms=3 right now...") );
1105 :
1106 : }// end of switch nterms
1107 :
1108 : /////// Note : To test this, apply it to the pbcoeffs_p and write to disk. They should be ones.
1109 :
1110 : return;
1111 : }
1112 :
1113 : */
1114 :
1115 : /*
1116 : // invert, apply to residuals
1117 : void NewMultiTermFT::applyWideBandPB(PtrBlock<SubImage<Float> *> &modelImageVec)
1118 : {
1119 : AlwaysAssert( pbcoeffs_p.nelements()==nterms_p , AipsError );
1120 : AlwaysAssert( modelImageVec.nelements()>=nterms_p , AipsError );
1121 :
1122 : multiplyHMatrix( hess_p, modelImageVec, pbcoeffs_p, String("/coeffPB_") );
1123 : }
1124 : */
1125 :
1126 :
1127 : //---------------------------------------------------------------------------------------------------
1128 : //---------------------------------------------------------------------------------------------------
1129 : //---------------------------------------------------------------------------------------------------
1130 :
1131 : } //# NAMESPACE CASA - END
1132 :
|
