Line data Source code
1 : //# MultiTermFT.cc: Implementation of MultiTermFT 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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <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 <lattices/Lattices/ArrayLattice.h>
45 : //#include <measures/Measures/UVWMachine.h>
46 : //#include <lattices/Lattices/SubLattice.h>
47 : //#include <lattices/LRegions/LCBox.h>
48 : //#include <lattices/Lattices/LatticeCache.h>
49 : //#include <lattices/LatticeMath/LatticeFFT.h>
50 : //#include <lattices/Lattices/LatticeIterator.h>
51 : //#include <lattices/Lattices/LatticeStepper.h>
52 : //#include <scimath/Mathematics/ConvolveGridder.h>
53 : //#include <casa/Utilities/CompositeNumber.h>
54 : #include <casacore/casa/OS/Timer.h>
55 : #include <sstream>
56 :
57 : #include <synthesis/TransformMachines/MultiTermFT.h>
58 : #include <synthesis/TransformMachines/VisModelData.h>
59 :
60 : // This is the list of FTMachine types supported by MultiTermFT
61 : #include <synthesis/TransformMachines/GridFT.h>
62 : #include <synthesis/TransformMachines/WProjectFT.h>
63 : //#include <synthesis/MeasurementComponents/AWProjectWBFT.h>
64 : //#include <synthesis/MeasurementComponents/rGridFT.h>
65 :
66 : using namespace casacore;
67 : namespace casa { //# NAMESPACE CASA - BEGIN
68 : //----------------------------------------------------------------------
69 : //-------------------- constructors and descructors ----------------------
70 : //----------------------------------------------------------------------
71 0 : MultiTermFT::MultiTermFT(FTMachine *subftm, String subFTMname, Int nterms, Double reffreq)
72 : :FTMachine(), subftm_p(subftm), subFTMname_p(subFTMname), nterms_p(nterms),
73 0 : thisterm_p(0), reffreq_p(reffreq), imweights_p(Matrix<Float>(0,0)), machineName_p("MultiTermFT")
74 : {
75 0 : dbg_p=false;
76 0 : dotime_p=false;
77 0 : this->setBasePrivates(*subftm_p);
78 0 : canComputeResiduals_p = subftm_p->canComputeResiduals();
79 0 : if(dbg_p) cout << "MTFT :: constructor with subftm : "<< subFTMname_p << endl;
80 0 : if(dbg_p) cout << "can compute residuals : " << canComputeResiduals_p << endl;
81 :
82 0 : sumwt_p=0.0;
83 :
84 0 : time_get=0.0;
85 0 : time_put=0.0;
86 0 : time_res=0.0;
87 0 : }
88 :
89 : //----------------------------------------------------------------------
90 : // Construct from the input state record
91 0 : MultiTermFT::MultiTermFT(const RecordInterface& stateRec)
92 0 : : FTMachine()
93 : {
94 0 : String error;
95 0 : if (!fromRecord(error, stateRec)) {
96 0 : throw (AipsError("Failed to create gridder: " + error));
97 : };
98 0 : }
99 :
100 : //----------------------------------------------------------------------
101 : // Copy constructor
102 0 : MultiTermFT::MultiTermFT(const MultiTermFT& other) : FTMachine(), machineName_p("MultiTermFT")
103 : {
104 0 : operator=(other);
105 0 : }
106 :
107 0 : MultiTermFT& MultiTermFT::operator=(const MultiTermFT& other)
108 : {
109 :
110 0 : dbg_p = other.dbg_p;
111 0 : dotime_p = other.dotime_p;
112 0 : if(dbg_p) cout << "In MTFT operator= for " << other.subftm_p->name() << endl;
113 :
114 0 : if(this!=&other)
115 : {
116 0 : FTMachine::operator=(other);
117 :
118 : //if(other.subftm_p==0) throw(AipsError("Internal Error : Empty subFTMachine"));
119 :
120 : // The operator= would have copied only the pointer, and not made a new instance of subftm_p
121 : // Make the new instance of subftm_p here.
122 : // Ideally, this should call "clone" of that ftm, so that the if/else stuff can all go away.
123 0 : if(other.subFTMname_p=="GridFT")
124 0 : { subftm_p = new GridFT(static_cast<const GridFT&>(*other.subftm_p)); }
125 0 : else if(other.subFTMname_p=="WProjectFT")
126 0 : { subftm_p = new WProjectFT(static_cast<const WProjectFT&>(*other.subftm_p)); }
127 : /* else if(other.subFTMname_p=="wbawp")
128 : { subftm_p = new AWProjectWBFT(static_cast<const AWProjectWBFT&>(*other.subftm_p)); }
129 : else if(other.subFTMname_p=="nift")
130 : {
131 : subftm_p = ((rGridFT*)&(*other.subftm_p))->clone();
132 : //cout << "MTFT : copy constructor : newft->visresampler_p : " << &(*(subftm_p->visResampler_p)) << endl;
133 : }
134 : */
135 : else
136 0 : {throw(AipsError("FTMachine "+other.subFTMname_p+" is not supported with MS-MFS")); }
137 :
138 : // subftm_p->setBasePrivates(*this);
139 :
140 : // Copy local privates
141 0 : subFTMname_p = other.subFTMname_p;
142 0 : machineName_p = other.machineName_p;
143 0 : nterms_p = other.nterms_p;
144 0 : thisterm_p = other.thisterm_p;
145 0 : reffreq_p = other.reffreq_p;
146 0 : sumwt_p = other.sumwt_p;
147 :
148 : // Check if the sub ftm can calculate its own residuals....
149 0 : canComputeResiduals_p = subftm_p->canComputeResiduals();
150 : }
151 :
152 0 : return * this;
153 : }
154 :
155 : //----------------------------------------------------------------------
156 0 : MultiTermFT::~MultiTermFT()
157 : {
158 : // if(dbg_p) cerr << "MTFT :: destructor for term " << thisterm_p << " - deletes subftm explicitly " << endl;
159 0 : if(dbg_p) cout << "MTFT :: destructor for term " << thisterm_p << " - assumes automatic deletion of subftm " << endl;
160 : //if(subftm_p) { delete subftm_p; subftm_p=0; }
161 0 : }
162 :
163 :
164 : //---------------------------------------------------------------------------------------------------
165 : //------------ Multi-Term Specific Functions --------------------------------------------------------
166 : //---------------------------------------------------------------------------------------------------
167 :
168 : // Multiply the imaging weights by Taylor functions - in place
169 : // This function MUST be called in ascending Taylor-term order
170 : // NOTE : Add checks to ensure this.
171 0 : Bool MultiTermFT::modifyVisWeights(VisBuffer& vb)
172 : {
173 0 : if( thisterm_p > 0 )
174 : {
175 :
176 0 : if(imweights_p.shape() != vb.imagingWeight().shape())
177 0 : imweights_p.resize(vb.imagingWeight().shape());
178 0 : imweights_p = vb.imagingWeight();
179 :
180 0 : Float freq=0.0,mulfactor=1.0;
181 0 : Vector<Double> selfreqlist(vb.frequency());
182 :
183 0 : for (Int row=0; row<vb.nRow(); row++)
184 0 : for (Int chn=0; chn<vb.nChannel(); chn++)
185 : {
186 0 : freq = selfreqlist(IPosition(1,chn));
187 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
188 0 : (vb.imagingWeight())(chn,row) *= pow( mulfactor, thisterm_p );
189 : // sumwt_p += (vb.imagingWeight())(chn,row);
190 : }
191 : }
192 : /*
193 : else
194 : {
195 : for (Int row=0; row<vb.nRow(); row++)
196 : for (Int chn=0; chn<vb.nChannel(); chn++)
197 : {
198 : sumwt_p += (vb.imagingWeight())(chn,row);
199 : }
200 : }
201 : */
202 0 : return true;
203 : }
204 :
205 : // Reset the imaging weights back to their original values
206 : // to be called just after "put"
207 0 : Bool MultiTermFT::restoreImagingWeights(VisBuffer &vb)
208 : {
209 0 : if(thisterm_p>0)
210 : {
211 0 : AlwaysAssert( imweights_p.shape() == vb.imagingWeight().shape() ,AipsError);
212 0 : vb.imagingWeight() = imweights_p;
213 : }
214 0 : return true;
215 : }
216 :
217 :
218 : // Multiply the model visibilities by the Taylor functions - in place.
219 0 : Bool MultiTermFT::modifyModelVis(VisBuffer& vb)
220 : {
221 0 : if( thisterm_p > 0 )
222 : {
223 0 : Float freq=0.0,mulfactor=1.0;
224 0 : Vector<Double> selfreqlist(vb.frequency());
225 :
226 0 : for (uInt pol=0; pol< uInt((vb.modelVisCube()).shape()[0]); pol++)
227 0 : for (uInt chn=0; chn< uInt(vb.nChannel()); chn++)
228 0 : for (uInt row=0; row< uInt(vb.nRow()); row++)
229 : {
230 0 : freq = selfreqlist(IPosition(1,chn));
231 0 : mulfactor = ((freq-reffreq_p)/reffreq_p);
232 0 : (vb.modelVisCube())(pol,chn,row) *= pow(mulfactor,thisterm_p);
233 : }
234 : }
235 :
236 0 : return true;
237 : }
238 :
239 :
240 : //---------------------------------------------------------------------------------------------------
241 : //---------------------- Prediction and De-gridding -----------------------------------
242 : //---------------------------------------------------------------------------------------------------
243 0 : void MultiTermFT::initializeToVis(ImageInterface<Complex>& iimage,
244 : const VisBuffer& vb)
245 : {
246 0 : if(dbg_p) cout << "MTFT::initializeToVis for term " << thisterm_p << endl;
247 0 : subftm_p->initializeToVis(iimage,vb);
248 0 : time_get=0.0;
249 0 : }
250 :
251 0 : void MultiTermFT::initMaps(const VisBuffer& vb){
252 0 : subftm_p->initMaps(vb);
253 0 : }
254 :
255 0 : void MultiTermFT::get(VisBuffer& vb, Int row)
256 : {
257 0 : if(dotime_p) tmr_p.mark();
258 :
259 0 : subftm_p->get(vb,row);
260 0 : modifyModelVis(vb);
261 :
262 0 : if(dotime_p) time_get += tmr_p.real();
263 0 : }
264 :
265 0 : void MultiTermFT::finalizeToVis()
266 : {
267 0 : if(dbg_p) cout << "MTFT::finalizeToVis for term " << thisterm_p <<endl;
268 0 : subftm_p->finalizeToVis();
269 0 : if(dotime_p) cout << " taylor " << thisterm_p << "*************** get time : " << time_get << endl;
270 0 : }
271 :
272 :
273 : //---------------------------------------------------------------------------------------------------
274 : //---------------------- Calculate Residual Visibilities -------------------------------
275 : //---------------------------------------------------------------------------------------------------
276 0 : void MultiTermFT::ComputeResiduals(VisBuffer &vb, Bool useCorrected)
277 : {
278 0 : if(dotime_p) tmr_p.mark();
279 :
280 0 : if(subftm_p->canComputeResiduals()) subftm_p->ComputeResiduals(vb,useCorrected);
281 0 : else throw(AipsError("MultiTerm::ComputeResiduals : subftm of MultiTermFT cannot compute its own residuals !"));
282 :
283 0 : if(dotime_p) time_res += tmr_p.real();
284 0 : }
285 :
286 : //---------------------------------------------------------------------------------------------------
287 : //---------------------- Gridding --------------------------------------------------------------
288 : //---------------------------------------------------------------------------------------------------
289 0 : void MultiTermFT::initializeToSky(ImageInterface<Complex>& iimage,
290 : Matrix<Float>& weight, const VisBuffer& vb)
291 : {
292 0 : if(dbg_p) cout << "MTFT::initializeToSky for term " << thisterm_p << endl;
293 :
294 0 : subftm_p->initializeToSky(iimage,weight,vb);
295 :
296 : // sumwt_p=0.0;
297 :
298 0 : if(dotime_p) {time_put=0.0; time_res=0.0;}
299 0 : }
300 :
301 0 : void MultiTermFT::put(VisBuffer& vb, Int row, Bool dopsf,
302 : FTMachine::Type type) //, const Matrix<Float>& imwght) // Last parameter is ignored.
303 : {
304 0 : if(dotime_p) tmr_p.mark();
305 :
306 0 : modifyVisWeights(vb);
307 0 : subftm_p->put(vb,row,dopsf,type);
308 0 : restoreImagingWeights(vb);
309 :
310 0 : if(dotime_p) time_put += tmr_p.real();
311 0 : }
312 :
313 0 : void MultiTermFT::finalizeToSky()
314 : {
315 0 : if(dbg_p) cout << "MTFT::finalizeToSky for term " << thisterm_p << endl;
316 :
317 0 : subftm_p->finalizeToSky();
318 :
319 : // cout << "***** sumwt : " << sumwt_p << endl;
320 :
321 0 : if(dotime_p) cout << " taylor " << thisterm_p << "*************** can compute residual " << canComputeResiduals_p << " res time : " << time_res << " put time :" << time_put << endl;
322 0 : }
323 :
324 : //---------------------------------------------------------------------------------------------------
325 : //----------------------------- Obtain Images -----------------------------------------------------
326 : //---------------------------------------------------------------------------------------------------
327 0 : ImageInterface<Complex>& MultiTermFT::getImage(Matrix<Float>& weights, Bool normalize)
328 : {
329 0 : if(dbg_p) cout << "MTFT :: getImage for term " << thisterm_p << endl;
330 0 : return subftm_p->getImage(weights,normalize);
331 : }
332 :
333 : //----------------------------------------------------------------------
334 0 : void MultiTermFT::getWeightImage(ImageInterface<Float>& weightImage, Matrix<Float>& weights)
335 : {
336 0 : if(dbg_p) cout << "MTFT :: getWeightImage for term " << thisterm_p << endl;
337 0 : subftm_p->getWeightImage(weightImage, weights);
338 0 : }
339 :
340 : //----------------------------------------------------------------------
341 0 : void MultiTermFT::makeImage(FTMachine::Type type, VisSet& vs,
342 : ImageInterface<Complex>& theImage, Matrix<Float>& weight)
343 : {
344 0 : if(dbg_p) cout << "MTFT :: makeImage for term " << thisterm_p << endl;
345 0 : subftm_p->makeImage(type, vs, theImage, weight);
346 0 : }
347 :
348 : //---------------------------------------------------------------------------------------------------
349 : //------------------------ To / From Records ---------------------------------------------------------
350 : //---------------------------------------------------------------------------------------------------
351 0 : Bool MultiTermFT::toRecord(String& error, RecordInterface& outRec, Bool withImage, const String diskimage)
352 : {
353 0 : if(dbg_p) cout << "MTFT :: toRecord for term " << thisterm_p << endl;
354 0 : Bool retval = true;
355 : //no image is held in this machine so no image need to be saved
356 : // the subftm holds the image
357 0 : if(!FTMachine::toRecord(error, outRec, false))
358 0 : return false;
359 :
360 0 : Record subFTContainer;
361 0 : String elimage="";
362 0 : if(diskimage != ""){
363 0 : elimage=diskimage+String("_")+ String::toString(rand());
364 0 : while(Table::isReadable(elimage))
365 0 : elimage=diskimage+String("_")+ String::toString(rand());
366 : }
367 0 : subftm_p->toRecord(error, subFTContainer,withImage, elimage);
368 :
369 0 : outRec.defineRecord("subftm",subFTContainer);
370 0 : outRec.define("subftname", subFTMname_p);
371 0 : outRec.define("nterms",nterms_p);
372 0 : outRec.define("thisterm",thisterm_p);
373 0 : outRec.define("reffreq",reffreq_p);
374 0 : outRec.define("dbg", dbg_p);
375 0 : outRec.define("dotime", dotime_p);
376 0 : outRec.define("time_get", time_get);
377 0 : outRec.define("time_put", time_put);
378 0 : outRec.define("time_res", time_res);
379 :
380 0 : return retval;
381 : }
382 :
383 : //---------------------------------------------------------------------------------------------------
384 0 : Bool MultiTermFT::fromRecord(String& error, const RecordInterface& inRec)
385 : {
386 :
387 0 : Bool retval = true;
388 :
389 0 : if(!FTMachine::fromRecord(error, inRec))
390 0 : return false;
391 :
392 0 : Record subFTMRec=inRec.asRecord("subftm");
393 0 : subftm_p=VisModelData::NEW_FT(subFTMRec);
394 0 : if (subftm_p.null())
395 0 : return false;
396 :
397 :
398 0 : inRec.get("subftname",subFTMname_p);
399 0 : inRec.get("nterms",nterms_p);
400 0 : inRec.get("thisterm",thisterm_p);
401 : //if(dbg_p) cout << "MTFT :: fromRecord for term " << thisterm_p << endl;
402 0 : inRec.get("reffreq",reffreq_p);
403 0 : inRec.get("dbg", dbg_p);
404 0 : inRec.get("dotime", dotime_p);
405 0 : inRec.get("time_get", time_get);
406 0 : inRec.get("time_put", time_put);
407 0 : inRec.get("time_res", time_res);
408 0 : machineName_p="MultiTermFT";
409 0 : return retval;
410 : }
411 : //---------------------------------------------------------------------------------------------------
412 :
413 : } //# NAMESPACE CASA - END
414 :
|
