Line data Source code
1 : //# CubeSkyEquation.cc: Implementation of Cube Optimized Sky Equation classes
2 : //# Copyright (C) 2007
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 <iostream>
29 : #include <casacore/casa/Exceptions/Error.h>
30 : #include <casacore/casa/Utilities/Assert.h>
31 : #include <casacore/casa/BasicSL/Complex.h>
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/MatrixMath.h>
35 : #include <casacore/casa/OS/HostInfo.h>
36 : #include <casacore/casa/System/ProgressMeter.h>
37 : #include <casacore/casa/Utilities/CountedPtr.h>
38 :
39 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
40 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
41 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
42 : #include <casacore/lattices/LEL/LatticeExpr.h>
43 :
44 : #include <synthesis/MeasurementComponents/SkyModel.h>
45 : #include <synthesis/MeasurementEquations/CubeSkyEquation.h>
46 : #include <synthesis/TransformMachines/SkyJones.h>
47 : #include <synthesis/TransformMachines/FTMachine.h>
48 : #include <synthesis/TransformMachines/SetJyGridFT.h>
49 : #include <synthesis/TransformMachines/GridFT.h>
50 : #include <synthesis/TransformMachines/MosaicFT.h>
51 : #include <synthesis/TransformMachines/MultiTermFT.h>
52 : #include <synthesis/TransformMachines/NewMultiTermFT.h>
53 : #include <synthesis/MeasurementComponents/GridBoth.h>
54 : #include <synthesis/TransformMachines/WProjectFT.h>
55 : #include <synthesis/MeasurementComponents/nPBWProjectFT.h>
56 : #include <synthesis/TransformMachines/AWProjectFT.h>
57 : #include <synthesis/TransformMachines/AWProjectWBFT.h>
58 : #include <synthesis/MeasurementComponents/PBMosaicFT.h>
59 : #include <synthesis/TransformMachines/WPConvFunc.h>
60 : #include <synthesis/TransformMachines/SimplePBConvFunc.h>
61 : #include <synthesis/TransformMachines/ComponentFTMachine.h>
62 : #include <synthesis/TransformMachines/SynthesisError.h>
63 : #include <synthesis/TransformMachines/StokesImageUtil.h>
64 : #include <synthesis/TransformMachines/Utils.h>
65 : #include <synthesis/Utilities/SigHandler.h>
66 : #include <casacore/images/Images/ImageInterface.h>
67 : #include <casacore/images/Images/SubImage.h>
68 :
69 : #include <msvis/MSVis/StokesVector.h>
70 : #include <msvis/MSVis/VisBufferUtil.h>
71 : #include <msvis/MSVis/VisSet.h>
72 : #include <msvis/MSVis/VisibilityIterator.h>
73 : #include <msvis/MSVis/VisBuffer.h>
74 : #ifdef _OPENMP
75 : #include <omp.h>
76 : #endif
77 :
78 : using namespace casacore;
79 : namespace casa { //# NAMESPACE CASA - BEGIN
80 :
81 0 : CubeSkyEquation::CubeSkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft,
82 0 : ComponentFTMachine& cft, Bool noModelCol)
83 : : SkyEquation(sm, vs, ft, cft, noModelCol),
84 0 : destroyVisibilityIterator_p (false),
85 0 : internalChangesPut_p(false),
86 0 : internalChangesGet_p(false),
87 0 : firstOneChangesPut_p(false),
88 0 : firstOneChangesGet_p(false)
89 : {
90 :
91 0 : init(ft);
92 :
93 0 : }
94 :
95 130 : CubeSkyEquation::CubeSkyEquation(SkyModel& sm, ROVisibilityIterator& vi, FTMachine& ft,
96 130 : ComponentFTMachine& cft, Bool noModelCol)
97 : : SkyEquation(sm, vi, ft, cft, noModelCol),
98 130 : destroyVisibilityIterator_p (false),
99 130 : internalChangesPut_p(false),
100 130 : internalChangesGet_p(false),
101 130 : firstOneChangesPut_p(false),
102 130 : firstOneChangesGet_p(false)
103 : {
104 130 : init(ft);
105 130 : }
106 :
107 130 : void CubeSkyEquation::init(FTMachine& ft){
108 130 : Int nmod=sm_->numberOfModels()/sm_->numberOfTaylorTerms();
109 :
110 130 : doflat_p=false;
111 :
112 : /// if(sm_->numberOfTaylorTerms()>1)
113 130 : if( ft.name()=="MultiTermFT" )
114 : {
115 0 : nmod = (sm_->numberOfModels()/sm_->numberOfTaylorTerms()) * (2 * sm_->numberOfTaylorTerms() - 1);
116 : }
117 :
118 : //case of component ft only
119 130 : if(nmod==0)
120 97 : nmod=1;
121 :
122 130 : ftm_p.resize(nmod, true);
123 130 : iftm_p.resize(nmod, true);
124 :
125 : //make a distinct ift_ as gridding and degridding can occur simultaneously
126 130 : if(ft.name() == "MosaicFT"){
127 0 : ft_=new MosaicFT(static_cast<MosaicFT &>(ft));
128 0 : ift_=new MosaicFT(static_cast<MosaicFT &>(ft));
129 0 : ftm_p[0]=ft_;
130 0 : iftm_p[0]=ift_;
131 : //For mosaic ...outlier fields get normal GridFT's
132 :
133 0 : MPosition loc=ift_->getLocation();
134 0 : for (Int k=1; k < (nmod); ++k){
135 0 : ftm_p[k]=new GridFT(1000000, 16, "SF", loc, 1.0, false);
136 0 : iftm_p[k]=new GridFT(1000000, 16, "SF", loc, 1.0, false);
137 : }
138 0 : }
139 130 : else if(ft.name()== "WProjectFT"){
140 0 : ft_=new WProjectFT(static_cast<WProjectFT &>(ft));
141 0 : ift_=new WProjectFT(static_cast<WProjectFT &>(ft));
142 0 : ftm_p[0]=ft_;
143 0 : iftm_p[0]=ift_;
144 0 : CountedPtr<WPConvFunc> sharedconvFunc= new WPConvFunc( *(static_cast<WProjectFT &>(ft).getConvFunc()));
145 0 : static_cast<WProjectFT &>(ft).setConvFunc(sharedconvFunc);
146 0 : static_cast<WProjectFT &>(*ftm_p[0]).setConvFunc(sharedconvFunc);
147 0 : static_cast<WProjectFT &>(*iftm_p[0]).setConvFunc(sharedconvFunc);
148 : // For now have all the fields have WProjectFt machines....
149 : //but should be seperated between GridFT's for the outliers and
150 : //WProject for the facets.
151 0 : for (Int k=1; k < (nmod); ++k){
152 0 : ftm_p[k]=new WProjectFT(static_cast<WProjectFT &>(*ft_));
153 0 : iftm_p[k]=new WProjectFT(static_cast<WProjectFT &>(*ift_));
154 : // Give each pair of FTMachine a convolution function set to share
155 0 : static_cast<WProjectFT &>(*ftm_p[k]).setConvFunc(sharedconvFunc);
156 0 : static_cast<WProjectFT &>(*iftm_p[k]).setConvFunc(sharedconvFunc);
157 :
158 : }
159 0 : }
160 130 : else if(ft.name()== "GridBoth"){
161 0 : ft_=new GridBoth(static_cast<GridBoth &>(ft));
162 0 : ift_=new GridBoth(static_cast<GridBoth &>(ft));
163 0 : ftm_p[0]=ft_;
164 0 : iftm_p[0]=ift_;
165 0 : if(nmod > 1){
166 0 : throw(AipsError("No multifield with joint gridding allowed"));
167 : }
168 : }
169 130 : else if(ft.name()== "PBWProjectFT"){
170 0 : ft_=new nPBWProjectFT(static_cast<nPBWProjectFT &>(ft));
171 0 : ift_=new nPBWProjectFT(static_cast<nPBWProjectFT &>(ft));
172 0 : ftm_p[0]=ft_;
173 0 : iftm_p[0]=ift_;
174 0 : if(nmod != (2 * sm_->numberOfTaylorTerms() - 1)) /* MFS */
175 0 : throw(AipsError("No multifield with pb-projection allowed"));
176 0 : for (Int k=1; k < (nmod); ++k){
177 0 : ftm_p[k]=new nPBWProjectFT(static_cast<nPBWProjectFT &>(*ft_));
178 0 : iftm_p[k]=new nPBWProjectFT(static_cast<nPBWProjectFT &>(*ift_));
179 : }
180 : }
181 130 : else if(ft.name()== "AWProjectFT"){
182 0 : ft_=new AWProjectFT(static_cast<AWProjectFT &>(ft));
183 0 : ift_=new AWProjectFT(static_cast<AWProjectFT &>(ft));
184 : // ift_=ft_;
185 0 : ftm_p[0]=ft_;
186 0 : iftm_p[0]=ift_;
187 0 : if(nmod != (2 * sm_->numberOfTaylorTerms() - 1)) /* MFS */
188 0 : throw(AipsError("No multifield with a-projection allowed"));
189 0 : for (Int k=1; k < (nmod); ++k){
190 0 : ftm_p[k]=new AWProjectFT(static_cast<AWProjectFT &>(*ft_));
191 0 : iftm_p[k]=new AWProjectFT(static_cast<AWProjectFT &>(*ift_));
192 : // iftm_p[k]=ftm_p[k];
193 : }
194 : }
195 130 : else if(ft.name()== "AWProjectWBFT"){
196 0 : ft_=new AWProjectWBFT(static_cast<AWProjectWBFT &>(ft));
197 0 : ift_=new AWProjectWBFT(static_cast<AWProjectWBFT &>(ft));
198 : // ift_=ft_;
199 0 : ftm_p[0]=ft_;
200 0 : iftm_p[0]=ift_;
201 : // if(nmod != (2 * sm_->numberOfTaylorTerms() - 1)) /* MFS */
202 : // throw(AipsError("No multifield with a-projection allowed"));
203 0 : for (Int k=1; k < (nmod); ++k){
204 0 : ftm_p[k]=new AWProjectWBFT(static_cast<AWProjectWBFT &>(*ft_));
205 0 : iftm_p[k]=new AWProjectWBFT(static_cast<AWProjectWBFT &>(*ift_));
206 0 : if(sm_->numberOfTaylorTerms()>1)
207 : {
208 0 : for (Int model=0; model < (sm_->numberOfModels()) ; ++model)
209 : {
210 0 : ftm_p[model]->setMiscInfo(sm_->getTaylorIndex(model));
211 0 : iftm_p[model]->setMiscInfo(sm_->getTaylorIndex(model));
212 : }
213 : }
214 : // iftm_p[k]=ftm_p[k];
215 : // if (rvi_p != NULL) cerr << rvi_p->getMSSelectionObj(0).getChanList();
216 : // if (wvi_p != NULL) cerr << rvi_p->getMSSelectionObj(0).getChanList();
217 : }
218 : }
219 130 : else if(ft.name()== "PBMosaicFT"){
220 0 : ft_=new PBMosaicFT(static_cast<PBMosaicFT &>(ft));
221 0 : ift_=new PBMosaicFT(static_cast<PBMosaicFT &>(ft));
222 0 : ftm_p[0]=ft_;
223 0 : iftm_p[0]=ift_;
224 0 : if(nmod != (2 * sm_->numberOfTaylorTerms() - 1)) /* MFS */
225 0 : throw(AipsError("No multifield with pb-mosaic allowed"));
226 0 : for (Int k=1; k < (nmod); ++k){
227 0 : ftm_p[k]=new PBMosaicFT(static_cast<PBMosaicFT &>(*ft_));
228 0 : iftm_p[k]=new PBMosaicFT(static_cast<PBMosaicFT &>(*ift_));
229 : }
230 : }
231 130 : else if (ft.name() == "SetJyGridFT") {
232 4 : ft_=new SetJyGridFT(static_cast<SetJyGridFT &>(ft));
233 4 : ift_=new SetJyGridFT(static_cast<SetJyGridFT &>(ft));
234 : // ftm_p[0]=CountedPtr<FTMachine>(ft_, false);
235 4 : ftm_p[0]=ft_;
236 4 : iftm_p[0]=ift_;
237 4 : for (Int k=1; k < (nmod); ++k){
238 0 : ftm_p[k]=new SetJyGridFT(static_cast<SetJyGridFT &>(*ft_));
239 0 : iftm_p[k]=new SetJyGridFT(static_cast<SetJyGridFT &>(*ift_));
240 : }
241 : }
242 126 : else if (ft.name() == "MultiTermFT") {
243 0 : ft_=new MultiTermFT(static_cast<MultiTermFT &>(ft));
244 0 : ift_=new MultiTermFT(static_cast<MultiTermFT &>(ft));
245 0 : ftm_p[0]=ft_;
246 0 : iftm_p[0]=ift_;
247 0 : for (Int k=1; k < (nmod); ++k){
248 0 : ftm_p[k]=new MultiTermFT(static_cast<MultiTermFT &>(*ft_));
249 0 : iftm_p[k]=new MultiTermFT(static_cast<MultiTermFT &>(*ift_));
250 : }
251 0 : for (Int k=0; k < (nmod); ++k){
252 0 : uInt tayindex = k/(sm_->numberOfModels()/sm_->numberOfTaylorTerms());
253 : //cout << "CubeSkyEqn : model : " << k << " : setting taylor index : " << tayindex << endl;
254 0 : ftm_p[k]->setMiscInfo(tayindex);
255 0 : iftm_p[k]->setMiscInfo(tayindex);
256 : }
257 : }
258 126 : else if (ft.name() == "NewMultiTermFT") {
259 0 : ft_=new NewMultiTermFT(static_cast<NewMultiTermFT &>(ft));
260 0 : ift_=new NewMultiTermFT(static_cast<NewMultiTermFT &>(ft));
261 0 : ftm_p[0]=ft_;
262 0 : iftm_p[0]=ift_;
263 0 : for (Int k=1; k < (nmod); ++k){
264 0 : ftm_p[k]=new NewMultiTermFT(static_cast<NewMultiTermFT &>(*ft_));
265 0 : iftm_p[k]=new NewMultiTermFT(static_cast<NewMultiTermFT &>(*ift_));
266 : }
267 : }
268 126 : else if (ft.name() == "SDGrid") {
269 0 : ft_=new SDGrid(static_cast<SDGrid &>(ft));
270 0 : ift_=new SDGrid(static_cast<SDGrid &>(ft));
271 0 : ftm_p[0]=ft_;
272 0 : iftm_p[0]=ift_;
273 0 : for (Int k=1; k < (nmod); ++k){
274 0 : ftm_p[k]=new SDGrid(static_cast<SDGrid &>(*ft_));
275 0 : iftm_p[k]=new SDGrid(static_cast<SDGrid &>(*ift_));
276 : }
277 : }
278 : else {
279 126 : ft_=new GridFT(static_cast<GridFT &>(ft));
280 126 : ift_=new GridFT(static_cast<GridFT &>(ft));
281 : // ftm_p[0]=CountedPtr<FTMachine>(ft_, false);
282 126 : ftm_p[0]=ft_;
283 126 : iftm_p[0]=ift_;
284 127 : for (Int k=1; k < (nmod); ++k){
285 1 : ftm_p[k]=new GridFT(static_cast<GridFT &>(*ft_));
286 1 : iftm_p[k]=new GridFT(static_cast<GridFT &>(*ift_));
287 : }
288 : }
289 :
290 : Int nmod2;
291 130 : if( ft.name() == "MultiTermFT" || ft.name() == "NewMultiTermFT" )
292 : {
293 0 : nmod2 = (sm_->numberOfModels()/sm_->numberOfTaylorTerms()) * (2 * sm_->numberOfTaylorTerms() - 1);
294 : }
295 : else
296 : {
297 130 : nmod2=nmod;
298 : }
299 :
300 130 : imGetSlice_p.resize(nmod2, true, false);
301 130 : imPutSlice_p.resize(nmod2, true, false);
302 130 : weightSlice_p.resize(nmod2, true, false);
303 :
304 130 : }
305 :
306 260 : CubeSkyEquation::~CubeSkyEquation(){
307 : //As we make an explicit ift_ in the constructor we need
308 : //to take care of it here...
309 : //if(ift_ && (ift_ != ft_))
310 : // delete ift_;
311 :
312 130 : if (destroyVisibilityIterator_p){
313 0 : delete rvi_p;
314 0 : rvi_p = NULL;
315 0 : delete (vb_p.release ()); // free up the associated VisBuffer
316 : }
317 130 : SigHandler::resetSignalHandlers();
318 260 : }
319 :
320 110 : void CubeSkyEquation::setPhaseCenterTime(const Double time){
321 110 : SkyEquation::setPhaseCenterTime(time);
322 124 : for(Int model=0; model < sm_->numberOfModels(); ++model){
323 14 : ftm_p[model]->setPhaseCenterTime(time);
324 14 : iftm_p[model]->setPhaseCenterTime(time);
325 : }
326 :
327 110 : }
328 :
329 116 : void CubeSkyEquation::predict(Bool incremental, MS::PredefinedColumns col) {
330 :
331 116 : VisibilityIterator::DataColumn visCol=VisibilityIterator::Model;
332 116 : if(col==MS::DATA){
333 0 : visCol=VisibilityIterator::Observed;
334 : }
335 116 : if(col==MS::CORRECTED_DATA){
336 0 : visCol=VisibilityIterator::Corrected;
337 : }
338 116 : AlwaysAssert(cft_, AipsError);
339 116 : AlwaysAssert(sm_, AipsError);
340 : //AlwaysAssert(vs_, AipsError);
341 116 : if(sm_->numberOfModels()!= 0) AlwaysAssert(ok(),AipsError);
342 : // if(noModelCol_p)
343 : // throw(AipsError("Cannot predict visibilities without using scratch columns yet"));
344 : // Initialize
345 116 : if(wvi_p==NULL)
346 0 : throw(AipsError("Cannot save model in non-writable ms"));
347 116 : VisIter& vi=*wvi_p;
348 : //Lets get the channel selection for later use
349 116 : vi.getChannelSelection(blockNumChanGroup_p, blockChanStart_p,
350 116 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
351 116 : checkVisIterNumRows(vi);
352 116 : VisBufferAutoPtr vb (vi); // uses write VI so no ROVIA conversion
353 116 : Bool changedVI=false;
354 : // Reset the visibilities only if this is not an incremental
355 : // change to the model
356 116 : Bool initialized=false;
357 116 : predictComponents(incremental, initialized);
358 : //set to zero then loop over model...check for size...subimage then loop over subimages
359 :
360 :
361 116 : Bool isEmpty=true;
362 136 : for (Int model=0; model < (sm_->numberOfModels());++model){
363 20 : isEmpty=isEmpty && (sm_->isEmpty(model));
364 :
365 : }
366 : ////if people want to use model but it isn't there..we'll ignore you
367 116 : if(!noModelCol_p)
368 90 : noModelCol_p=rvi_p->msColumns().modelData().isNull();
369 :
370 :
371 116 : if( (sm_->numberOfModels() >0) && isEmpty && !initialized && !incremental){
372 : // We are at the begining with an empty model as starting point
373 37 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
374 913 : for (vi.origin(); vi.more(); vi++) {
375 883 : if(!noModelCol_p){
376 703 : vb->setModelVisCube(Complex(0.0,0.0));
377 703 : vi.setVis(vb->modelVisCube(),visCol);
378 : }
379 : }
380 : }
381 : }
382 :
383 : //If all model is zero...no need to continue
384 116 : if(isEmpty)
385 104 : return;
386 : //TODO if nomodel set flat to 0.0
387 :
388 :
389 : // Now do the images
390 24 : for (Int model=0; model < (sm_->numberOfModels());++model){
391 : // Change the model polarization frame
392 12 : if(vb->polFrame()==MSIter::Linear) {
393 2 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
394 : StokesImageUtil::LINEAR);
395 : }
396 : else {
397 10 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
398 : StokesImageUtil::CIRCULAR);
399 : }
400 : //UUU/// ft_=&(*ftm_p[model]);
401 12 : scaleImage(model, incremental);
402 : }
403 12 : ft_=&(*ftm_p[0]);
404 : // Reset the various SkyJones
405 12 : resetSkyJones();
406 12 : Int nCubeSlice=1;
407 12 : isLargeCube(sm_->cImage(0), nCubeSlice);
408 24 : for (Int cubeSlice=0; cubeSlice< nCubeSlice; ++cubeSlice){
409 12 : changedVI= getFreqRange(vi, sm_->cImage(0).coordinates(),
410 12 : cubeSlice, nCubeSlice) || changedVI;
411 12 : vi.originChunks();
412 12 : vi.origin();
413 12 : initializeGetSlice(* vb, 0,incremental, cubeSlice, nCubeSlice);
414 : ///vb->newMS does not seem to be set to true with originchunks
415 : //so have to monitor msid
416 12 : Int oldmsid=-1;
417 34 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
418 650 : for (vi.origin(); vi.more(); vi++) {
419 628 : if(noModelCol_p){
420 : //cerr << "noModelCol_p " << noModelCol_p << " newms " << vb->newMS() << " nummod " << (sm_->numberOfModels()) << endl;
421 360 : if(vb->msId() != oldmsid){
422 2 : oldmsid=vb->msId();
423 4 : for (Int model=0; model < (sm_->numberOfModels());++model){
424 2 : Record ftrec;
425 2 : String error;
426 2 : String modImage=vi.ms().getPartNames()[0];
427 2 : if(!(vi.ms().source().isNull()))
428 2 : modImage=(vi.ms()).source().tableName();
429 2 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
430 : //cerr << "in ftrec saving" << endl;
431 2 : if(!(ftm_p[model]->toRecord(error, ftrec, true, modImage)))
432 0 : throw(AipsError("Error in record saving: "+error));
433 2 : vi.putModel(ftrec, false, ((model>0) || incremental || (cubeSlice > 0)));
434 2 : }
435 : }
436 : }
437 : else{
438 268 : if(!incremental&&!initialized && (cubeSlice==0)) {
439 232 : vb->setModelVisCube(Complex(0.0,0.0));
440 : }
441 : // get the model visibility and write it to the model MS
442 268 : getSlice(* vb,incremental, cubeSlice, nCubeSlice);
443 268 : vi.setVis(vb->modelVisCube(),visCol);
444 : }
445 : }
446 : }
447 12 : finalizeGetSlice();
448 12 : if(!incremental&&!initialized) initialized=true;
449 : }
450 :
451 24 : for(Int model=0; model < sm_->numberOfModels(); ++model){
452 : //For now unscale test on name of ft_
453 12 : ft_=&(*ftm_p[model]);
454 12 : unScaleImage(model, incremental);
455 : }
456 12 : ft_=&(*ftm_p[0]);
457 :
458 : //lets return original selection back to iterator
459 12 : if(changedVI)
460 0 : vi.selectChannel(blockNumChanGroup_p, blockChanStart_p,
461 0 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
462 :
463 116 : }
464 :
465 14 : void CubeSkyEquation::makeApproxPSF(PtrBlock<ImageInterface<Float> * >& psfs)
466 : {
467 :
468 14 : if(iftm_p[0]->name()=="MosaicFT")
469 0 : makeMosaicPSF(psfs);
470 : else
471 14 : makeSimplePSF(psfs);
472 :
473 14 : }
474 0 : void CubeSkyEquation::makeMosaicPSF(PtrBlock<ImageInterface<Float> * >& psfs){
475 : //lets try to make the psf directly
476 0 : LogIO os(LogOrigin("SkyEquation", "makeMosaicPSF"));
477 0 : Bool centered=true;
478 : try{
479 0 : makeSimplePSF(psfs);
480 : }
481 0 : catch(...){
482 0 : centered=false;
483 0 : }
484 : Int xpos;
485 : Int ypos;
486 0 : Matrix<Float> psfplane;
487 : Float peak;
488 0 : StokesImageUtil::locatePeakPSF(*(psfs[0]), xpos, ypos, peak, psfplane);
489 0 : Int nx=psfplane.shape()(0);
490 0 : Int ny=psfplane.shape()(1);
491 :
492 : //cerr << "xpos " << xpos << " " << ypos << " " << peak << endl;
493 : // lets ignore misers who made 10x10 pixel images
494 0 : centered=centered && (abs(xpos-nx/2) <=5) && (abs(ypos-ny/2) <=5);
495 :
496 :
497 0 : if(centered){
498 : //for cubes some of the planes may not have a central peak
499 0 : Int nchana= (psfs[0])->shape()(3);
500 0 : if(nchana > 1){
501 0 : IPosition blc(4,nx, ny, 0, nchana);
502 0 : IPosition trc(4, nx, ny, 0, nchana);
503 0 : blc(0)=0; blc(1)=0; trc(0)=nx-1; trc(1)=ny-1;
504 0 : Array<Float> goodplane(psfplane.reform(IPosition(4, nx,ny,1,1)));
505 0 : for (Int k=0; k < nchana ; ++k){
506 0 : blc(3)=k; trc(3)=k;
507 0 : Slicer sl(blc, trc, Slicer::endIsLast);
508 0 : SubImage<Float> psfSub(*(psfs[0]), sl, true);
509 : Float planeMax;
510 0 : LatticeExprNode LEN = max( psfSub );
511 0 : planeMax = LEN.getFloat();
512 0 : if( (planeMax >0.0) && (planeMax < 0.8 *peak)){
513 0 : psfSub.put(goodplane);
514 :
515 : }
516 0 : }
517 0 : }
518 0 : return;
519 : }
520 : //lets back up the ftmachines
521 0 : MosaicFT *ft_back= new MosaicFT(static_cast<MosaicFT &>(*ftm_p[0]));
522 0 : MosaicFT *ift_back = new MosaicFT(static_cast<MosaicFT &>(*iftm_p[0]));
523 : os << LogIO::WARN << "Mosaic psf is off. \nCould be no pointing in center of image \n"
524 : << "Will retry to make an approximate one without primary beam "
525 0 : << LogIO::POST;
526 0 : MPosition loc=iftm_p[0]->getLocation();
527 0 : ftm_p[0]=new GridFT(1000000, 16, "SF", loc, 1.0, false);
528 0 : iftm_p[0]=new GridFT(1000000, 16, "SF", loc, 1.0, false);
529 0 : ft_=&(*ftm_p[0]);
530 0 : ift_=&(*iftm_p[0]);
531 :
532 0 : ft_->setBasePrivates(*ft_back);
533 0 : ift_->setBasePrivates(*ift_back);
534 : // try again with simple ftmachines
535 0 : makeSimplePSF(psfs);
536 : //that 's the best psf you'd get
537 : //restore back MosaicFT machinas
538 0 : ftm_p[0]=ft_back;
539 0 : ft_=ft_back;
540 0 : iftm_p[0]=ift_back;
541 0 : ift_=ift_back;
542 0 : }
543 :
544 14 : void CubeSkyEquation::makeSimplePSF(PtrBlock<ImageInterface<Float> * >& psfs) {
545 :
546 14 : Int nmodels=psfs.nelements();
547 28 : LogIO os(LogOrigin("CubeSkyEquation", "makeSimplePSF"));
548 14 : SigHandler myStopSig;
549 14 : ft_->setNoPadding(noModelCol_p);
550 14 : isPSFWork_p= true; // avoid PB correction etc for PSF estimation
551 14 : Bool doPSF=true;
552 14 : Bool changedVI=false;
553 : // Initialize the gradients
554 14 : sm_->initializeGradients();
555 14 : ROVisIter& vi(*rvi_p);
556 : //Lets get the channel selection for later use
557 14 : vi.getChannelSelection(blockNumChanGroup_p, blockChanStart_p,
558 14 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
559 : // Reset the various SkyJones
560 14 : resetSkyJones();
561 14 : checkVisIterNumRows(vi);
562 : // Loop over all visibilities and pixels
563 14 : VisBufferAutoPtr vb (vi);
564 14 : vi.originChunks();
565 14 : vi.origin();
566 : // Change the model polarization frame
567 28 : for (Int model=0; model < nmodels; ++model){
568 14 : if(vb->polFrame()==MSIter::Linear) {
569 14 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
570 : StokesImageUtil::LINEAR);
571 : }
572 : else {
573 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
574 : StokesImageUtil::CIRCULAR);
575 : }
576 : }
577 :
578 :
579 14 : Int nCubeSlice=1;
580 14 : isLargeCube(sm_->cImage(0), nCubeSlice);
581 28 : for (Int cubeSlice=0; cubeSlice< nCubeSlice; ++cubeSlice){
582 14 : changedVI= getFreqRange(vi, sm_->cImage(0).coordinates(),
583 14 : cubeSlice, nCubeSlice) || changedVI;
584 14 : vi.originChunks();
585 14 : vi.origin();
586 14 : vb->invalidate();
587 14 : Int cohDone=0;
588 14 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
589 : "Gridding weights for PSF",
590 42 : "", "", "", true);
591 :
592 14 : initializePutSlice(* vb, doPSF, cubeSlice, nCubeSlice);
593 :
594 77 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
595 326 : for (vi.origin(); vi.more(); vi++) {
596 263 : if(myStopSig.gotStopSignal())
597 0 : throw(AipsError("Terminating...."));
598 263 : if(noModelCol_p) {
599 : //This here forces the modelVisCube shape and prevents reading model column
600 263 : vb->setModelVisCube(Complex(0.0,0.0));
601 : }
602 263 : putSlice(* vb, doPSF, FTMachine::MODEL, cubeSlice, nCubeSlice);
603 263 : cohDone+=vb->nRow();
604 263 : pm.update(Double(cohDone));
605 :
606 : }
607 : }
608 14 : finalizePutSlice(* vb, doPSF, cubeSlice, nCubeSlice);
609 14 : }
610 :
611 : //Don't need these for now
612 28 : for(Int model=0; model < nmodels; ++model){
613 :
614 14 : sm_->work(model).clearCache();
615 14 : sm_->cImage(model).clearCache();
616 : }
617 :
618 : //lets return original selection back to iterator
619 :
620 :
621 14 : if(changedVI)
622 0 : vi.selectChannel(blockNumChanGroup_p, blockChanStart_p,
623 0 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
624 14 : sm_->finalizeGradients();
625 14 : fixImageScale();
626 28 : for(Int model=0; model < nmodels; ++model){
627 : {
628 : //Normalize the gS image
629 14 : Int nXX=sm_->ggS(model).shape()(0);
630 14 : Int nYY=sm_->ggS(model).shape()(1);
631 14 : Int npola= sm_->ggS(model).shape()(2);
632 14 : Int nchana= sm_->ggS(model).shape()(3);
633 14 : IPosition blc(4,nXX, nYY, npola, nchana);
634 14 : IPosition trc(4, nXX, nYY, npola, nchana);
635 14 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
636 : //max weights per plane
637 28 : for (Int j=0; j < npola; ++j){
638 28 : for (Int k=0; k < nchana ; ++k){
639 :
640 14 : blc(2)=j; trc(2)=j;
641 14 : blc(3)=k; trc(3)=k;
642 14 : Slicer sl(blc, trc, Slicer::endIsLast);
643 14 : SubImage<Float> gSSub(sm_->gS(model), sl, false);
644 14 : SubImage<Float> ggSSub(sm_->ggS(model), sl, false);
645 14 : SubImage<Float> psfSub(*(psfs[model]), sl, true);
646 : Float planeMax;
647 14 : LatticeExprNode LEN = max( ggSSub );
648 14 : planeMax = LEN.getFloat();
649 14 : if(planeMax !=0){
650 14 : psfSub.copyData( (LatticeExpr<Float>)
651 56 : (iif(ggSSub > (0.0),
652 28 : (gSSub/planeMax),0.0)));
653 : }
654 : else{
655 0 : psfSub.set(0.0);
656 : }
657 14 : }
658 : }
659 : //
660 14 : }
661 :
662 : /*
663 : if(0){
664 : PagedImage<Float> thisScreen(psfs[model]->shape(), psfs[model]->coordinates(), String("ELPSF).psf"));
665 : LatticeExpr<Float> le(*psfs[model]);
666 : thisScreen.copyData(le);
667 : }
668 : */
669 14 : LatticeExprNode maxPSF=max(*psfs[model]);
670 14 : Float maxpsf=maxPSF.getFloat();
671 14 : if(abs(maxpsf-1.0) > 1e-3) {
672 : os << "Maximum of approximate PSF for field " << model << " = "
673 0 : << maxpsf << " : renormalizing to unity" << LogIO::POST;
674 : }
675 14 : if(maxpsf > 0.0 ){
676 28 : LatticeExpr<Float> len((*psfs[model])/maxpsf);
677 14 : psfs[model]->copyData(len);
678 14 : }
679 : else{
680 0 : if(sm_->numberOfTaylorTerms()>1) { /* MFS */
681 0 : os << "PSF calculation resulted in a PSF with its peak being 0 or less. This is ok for MS-MFS." << LogIO::POST;
682 : }
683 : else{
684 0 : throw(PSFZero("SkyEquation:: PSF calculation resulted in a PSF with its peak being 0 or less!"));
685 : }
686 : }
687 :
688 14 : sm_->PSF(model).clearCache();
689 14 : sm_->gS(model).clearCache();
690 14 : sm_->ggS(model).clearCache();
691 14 : }
692 :
693 14 : isPSFWork_p=false; // resetting this flag so that subsequent calculation uses
694 : // the right SkyJones correction;
695 14 : }
696 :
697 0 : void CubeSkyEquation::gradientsChiSquared(Bool /*incr*/, Bool commitModel){
698 0 : AlwaysAssert(cft_, AipsError);
699 0 : AlwaysAssert(sm_, AipsError);
700 : //AlwaysAssert(vs_, AipsError);
701 0 : Bool initialized=false;
702 0 : Bool changedVI=false;
703 :
704 : //For now we don't deal with incremental especially when having multi fields
705 0 : Bool incremental=false;
706 :
707 0 : predictComponents(incremental, initialized);
708 0 : Bool predictedComp=initialized;
709 :
710 0 : SigHandler myStopSig;
711 : ////if people want to use model but it isn't there
712 :
713 0 : if(!noModelCol_p)
714 0 : noModelCol_p=rvi_p->msColumns().modelData().isNull();
715 :
716 0 : ROVisibilityIterator * oldRvi = NULL;
717 0 : VisibilityIterator * oldWvi = NULL;
718 :
719 0 : configureAsyncIo (oldRvi, oldWvi);
720 :
721 :
722 : // Timers tInitGrad=Timers::getTime();
723 0 : sm_->initializeGradients();
724 : // Initialize
725 : //ROVisIter& vi=*rvi_p;
726 : //Lets get the channel selection for later use
727 : // Timers tGetChanSel=Timers::getTime();
728 0 : rvi_p->getChannelSelection(blockNumChanGroup_p, blockChanStart_p,
729 0 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
730 : // Timers tCheckVisRows=Timers::getTime();
731 0 : checkVisIterNumRows(*rvi_p);
732 :
733 0 : VisBufferAutoPtr vb (rvi_p);
734 :
735 : // Timers tVisAutoPtr=Timers::getTime();
736 :
737 : /**** Do we need to do this
738 : if( (sm_->isEmpty(0)) && !initialized && !incremental){
739 : // We are at the begining with an empty model as starting point
740 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
741 : for (vi.origin(); vi.more(); vi++) {
742 : vb.setModelVisCube(Complex(0.0,0.0));
743 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
744 : }
745 : }
746 : }
747 : */
748 0 : Bool isEmpty=true;
749 0 : for (Int model=0; model < (sm_->numberOfModels());++model){
750 0 : isEmpty=isEmpty && sm_->isEmpty(model);
751 : }
752 : // Now do the images
753 :
754 0 : for (Int model=0;model< (sm_->numberOfModels()); ++model) {
755 : // Don't bother with empty images
756 : // Change the model polarization frame
757 0 : if(vb->polFrame()==MSIter::Linear) {
758 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
759 : StokesImageUtil::LINEAR);
760 : }
761 : else {
762 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
763 : StokesImageUtil::CIRCULAR);
764 : }
765 : //scaleImage(model, incremental);
766 : ///UUU/// ft_=&(*ftm_p[model]);
767 0 : scaleImage(model);
768 : // Reset the various SkyJones
769 : }
770 : // Timers tChangeStokes=Timers::getTime();
771 :
772 0 : ft_=&(*ftm_p[0]);
773 0 : resetSkyJones();
774 0 : firstOneChangesPut_p=false;
775 0 : firstOneChangesGet_p=false;
776 :
777 0 : Int nCubeSlice=1;
778 :
779 0 : isLargeCube(sm_->cImage(0), nCubeSlice);
780 :
781 : // aInitGrad += tGetChanSel - tInitGrad;
782 : // aGetChanSel += tCheckVisRows - tGetChanSel;
783 : // aCheckVisRows += tVisAutoPtr - tCheckVisRows;
784 : // aChangeStokes += tChangeStokes - tVisAutoPtr;
785 :
786 0 : for (Int cubeSlice=0; cubeSlice< nCubeSlice; ++cubeSlice){
787 :
788 : // vi.originChunks();
789 : // vi.origin();
790 : //sliceCube(imGetSlice_p, model, cubeSlice, nCubeSlice, 1);
791 : //Redo the channel selection in case of chunked cube to match
792 : //data needed for gridding.
793 : // Timers tGetFreqRange=Timers::getTime();
794 0 : changedVI= getFreqRange(*rvi_p, sm_->cImage(0).coordinates(),
795 0 : cubeSlice, nCubeSlice) || changedVI;
796 :
797 : // Timers tOrigChunks=Timers::getTime();
798 0 : rvi_p->originChunks();
799 0 : rvi_p->origin();
800 0 : Bool useCorrected= !(vb->msColumns().correctedData().isNull());
801 :
802 : // Timers tVBInValid=Timers::getTime();
803 : //vb->invalidate();
804 :
805 : // Timers tInitGetSlice=Timers::getTime();
806 : // if(!isEmpty)
807 0 : if( ! isNewFTM() || ! isEmpty )
808 : {
809 0 : initializeGetSlice(* vb, 0, false, cubeSlice, nCubeSlice);
810 : }
811 : // Timers tInitPutSlice=Timers::getTime();
812 0 : initializePutSlice(* vb, false, cubeSlice, nCubeSlice);
813 : // Timers tDonePutSlice=Timers::getTime();
814 0 : Int cohDone=0;
815 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
816 : "Gridding residual",
817 0 : "", "", "", true);
818 : // aGetFreq += tOrigChunks - tGetFreqRange;
819 : // aOrigChunks += tVBInValid - tOrigChunks;
820 : // aVBInValid += tInitGetSlice - tVBInValid;
821 : // aInitGetSlice += tInitPutSlice - tInitGetSlice;
822 : // aInitPutSlice += tDonePutSlice - tInitPutSlice;
823 :
824 0 : Int oldmsid=-1;
825 0 : for (rvi_p->originChunks();rvi_p->moreChunks();rvi_p->nextChunk()) {
826 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++) {
827 : //cerr << "nrow " << vb->nRow() << " ";
828 : // Bool anychanIn=false;
829 : //for(Int model=0; model < (sm_->numberOfModels()) ; ++model)
830 : //anychanIn=anychanIn || iftm_p[model]->matchChannel(vb->spectralWindow(), *vb);
831 : //if(anychanIn)
832 : {
833 0 : if(myStopSig.gotStopSignal())
834 0 : throw(AipsError("Terminating..."));
835 : // Timers tInitModel=Timers::getTime();
836 0 : if(!incremental && !predictedComp && (cubeSlice==0 || (noModelCol_p) )) {
837 : //This here forces the modelVisCube shape and prevents reading model column
838 0 : vb->setModelVisCube(Complex(0.0,0.0));
839 : }
840 : // get the model visibility and write it to the model MS
841 : // Timers tGetSlice=Timers::getTime();
842 :
843 : // Timers tgetSlice=Timers::getTime();
844 0 : if( vb->newMS() )
845 : {
846 0 : useCorrected= !(vb->msColumns().correctedData().isNull());
847 : }
848 :
849 : //#pragma omp parallel default(shared)
850 : {
851 : // cerr << "num_threads " << omp_get_num_threads() << endl;
852 : //#pragma omp sections nowait
853 : {
854 : //#pragma omp section
855 0 : if(!isEmpty)
856 0 : getSlice(* vb, (predictedComp || incremental), cubeSlice, nCubeSlice);
857 : //#pragma omp section
858 : //if(!useCorrected)
859 : // vb->visCube();
860 : //else
861 : // vb->correctedVisCube();
862 : }
863 :
864 : } //saving the model for self-cal most probably
865 : // Timers tSetModel=Timers::getTime();
866 : // Timers tsetModel=Timers::getTime();
867 0 : if(commitModel && wvi_p != NULL){
868 : ///Bow commit model to disk or to the record keyword
869 0 : if(!noModelCol_p)
870 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
871 : else{
872 0 : if(vb->msId() != oldmsid){
873 0 : oldmsid=vb->msId();
874 0 : for (Int model=0; model < (sm_->numberOfModels());++model){
875 0 : Record ftrec;
876 0 : String error;
877 0 : String modImage=wvi_p->ms().getPartNames()[0];
878 0 : if(!(wvi_p->ms().source().isNull()))
879 0 : modImage=(wvi_p->ms()).source().tableName();
880 0 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
881 :
882 : //cerr << "in ftrec saving" << endl;
883 0 : if(!(ftm_p[model]->toRecord(error, ftrec, true, modImage)))
884 0 : throw(AipsError("Error in saving model; "+error));
885 0 : wvi_p->putModel(ftrec, false, ((model>0) || predictedComp || incremental || (cubeSlice >0)));
886 0 : }
887 : }
888 : }
889 : }
890 : // Now lets grid the -ve of residual
891 : // use visCube if there is no correctedData
892 : // Timers tGetRes=Timers::getTime();
893 0 : if (!iftm_p[0]->canComputeResiduals())
894 0 : if(!useCorrected) vb->modelVisCube()-=vb->visCube();
895 0 : else vb->modelVisCube()-=vb->correctedVisCube();
896 : else
897 0 : iftm_p[0]->ComputeResiduals(*vb,useCorrected);
898 :
899 :
900 : // Timers tPutSlice = Timers::getTime();
901 0 : putSlice(* vb, false, FTMachine::MODEL, cubeSlice, nCubeSlice);
902 0 : cohDone+=vb->nRow();
903 0 : pm.update(Double(cohDone));
904 : // Timers tDoneGridding=Timers::getTime();
905 : // aInitModel += tgetSlice - tInitModel;
906 : // aGetSlice += tsetModel - tgetSlice;
907 : // aSetModel += tGetRes - tsetModel;
908 : // aGetRes += tPutSlice - tGetRes;
909 : // aPutSlice += tDoneGridding - tPutSlice;
910 : // aExtra += tDoneGridding - tInitModel;
911 : }
912 : }
913 : }
914 :
915 : // Timers tFinalizeGetSlice=Timers::getTime();
916 0 : finalizeGetSlice();
917 0 : if(!incremental&&!initialized) initialized=true;
918 : // Timers tFinalizePutSlice=Timers::getTime();
919 0 : finalizePutSlice(* vb, false, cubeSlice, nCubeSlice);
920 : // Timers tDoneFinalizePutSlice=Timers::getTime();
921 :
922 : // aFinalizeGetSlice += tFinalizePutSlice - tFinalizeGetSlice;
923 : // aFinalizePutSlice += tDoneFinalizePutSlice - tFinalizePutSlice;
924 0 : }
925 :
926 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
927 : //unScaleImage(model, incremental);
928 0 : sm_->cImage(model).clearCache();
929 0 : sm_->gS(model).clearCache();
930 0 : sm_->ggS(model).clearCache();
931 0 : sm_->work(model).clearCache();
932 :
933 : //UUU// ft_=&(*ftm_p[model]);
934 0 : unScaleImage(model);
935 :
936 : }
937 0 : ft_=&(*ftm_p[0]);
938 :
939 0 : this->fixImageScale();
940 : //lets return original selection back to iterator
941 : // storeImg(String("stokesNormed.im"),sm_->gS(0));
942 0 : if(changedVI)
943 0 : rvi_p->selectChannel(blockNumChanGroup_p, blockChanStart_p,
944 0 : blockChanWidth_p, blockChanInc_p, blockSpw_p);
945 :
946 : // If using async, put things back the way they were.
947 :
948 0 : if (oldRvi != NULL){
949 :
950 0 : delete vb.release(); // get rid of local attached to Vi
951 : //vb_p.set (oldRvi); // reattach vb_p to the old vi
952 :
953 0 : if (oldWvi != NULL){
954 0 : delete wvi_p;
955 0 : wvi_p = oldWvi;
956 : }
957 : else{
958 0 : delete rvi_p; // kill the new vi
959 : }
960 :
961 0 : rvi_p = oldRvi; // make the old vi the current vi
962 0 : rvi_p->originChunks(); // reset it
963 0 : vb_p->attachToVisIter(* rvi_p);
964 : }
965 :
966 : // for (Int model=0;model< (sm_->numberOfModels()); ++model)
967 : // if (!isNewFTM(&(*ftm_p[model])))
968 : // {
969 : // Bool dopsf=false;
970 : // tmpWBNormalizeImage(dopsf);
971 : // }
972 :
973 :
974 : // cerr << "gradChiSq: "
975 : // << "InitGrad = " << aInitGrad.formatAverage().c_str() << " "
976 : // << "GetChanSel = " << aGetChanSel.formatAverage().c_str() << " "
977 : // << "ChangeStokes = " << aChangeStokes.formatAverage().c_str() << " "
978 : // << "CheckVisRows = " << aCheckVisRows.formatAverage().c_str() << " "
979 : // << "GetFreq = " << aGetFreq.formatAverage().c_str() << " "
980 : // << "OrigChunks = " << aOrigChunks.formatAverage().c_str() << " "
981 : // << "VBInValid = " << aVBInValid.formatAverage().c_str() << " "
982 : // << "InitGetSlice = " << aInitGetSlice.formatAverage().c_str() << " "
983 : // << "InitPutSlice = " << aInitPutSlice.formatAverage().c_str() << " "
984 : // << "PutSlice = " << aPutSlice.formatAverage().c_str() << " "
985 : // << "FinalGetSlice = " << aFinalizeGetSlice.formatAverage().c_str() << " "
986 : // << "FinalPutSlice = " << aFinalizePutSlice.formatAverage().c_str() << " "
987 : // << endl;
988 :
989 : // cerr << "VB loop: "
990 : // << "InitModel = " << aInitModel.formatAverage().c_str() << " "
991 : // << "GetSlice = " << aGetSlice.formatAverage().c_str() << " "
992 : // << "SetModel = " << aSetModel.formatAverage().c_str() << " "
993 : // << "GetRes = " << aGetRes.formatAverage().c_str() << " "
994 : // << "PutSlice = " << aPutSlice.formatAverage().c_str() << " "
995 : // << "Extra = " << aExtra.formatAverage().c_str() << " "
996 : // << endl;
997 0 : }
998 :
999 : void
1000 0 : CubeSkyEquation::configureAsyncIo (ROVisibilityIterator * & oldRvi, VisibilityIterator * & oldWvi)
1001 : {
1002 :
1003 : using namespace casacore;
1004 : using namespace casa::asyncio;
1005 :
1006 0 : oldRvi = NULL;
1007 0 : oldWvi = NULL;
1008 :
1009 : Bool isEnabled;
1010 0 : AipsrcValue<Bool>::find (isEnabled, "Imager.asyncio", false);
1011 : // Bool foundSetting = AipsrcValue<Bool>::find (isEnabled, "Imager.asyncio", false);
1012 :
1013 : //isEnabled = ! foundSetting || isEnabled; // let global flag call shots if setting not present
1014 : // For now (release 3.4) make imaging be explicitly turned on
1015 :
1016 0 : if (! (isEnabled && ROVisibilityIterator::isAsynchronousIoEnabled())){
1017 0 : return; // async i/o is not going to be used this time around
1018 : }
1019 :
1020 : // Async I/O is enabled globally and for imaging so prepare to replace the
1021 : // existing VIs with async implmentations.
1022 :
1023 : PrefetchColumns prefetchColumns =
1024 : PrefetchColumns::prefetchColumns (VisBufferComponents::Ant1,
1025 : VisBufferComponents::Ant2,
1026 : VisBufferComponents::ArrayId,
1027 : VisBufferComponents::DataDescriptionId,
1028 : VisBufferComponents::Direction1,
1029 : VisBufferComponents::Direction2,
1030 : VisBufferComponents::Feed1,
1031 : VisBufferComponents::Feed1_pa,
1032 : VisBufferComponents::Feed2,
1033 : VisBufferComponents::Feed2_pa,
1034 : VisBufferComponents::FieldId,
1035 : VisBufferComponents::FlagCube,
1036 : VisBufferComponents::Flag,
1037 : VisBufferComponents::FlagRow,
1038 : VisBufferComponents::Freq,
1039 : VisBufferComponents::NChannel,
1040 : VisBufferComponents::NCorr,
1041 : VisBufferComponents::NRow,
1042 : VisBufferComponents::ObservedCube,
1043 : VisBufferComponents::PhaseCenter,
1044 : VisBufferComponents::PolFrame,
1045 : VisBufferComponents::SpW,
1046 : VisBufferComponents::Time,
1047 : VisBufferComponents::Uvw,
1048 : VisBufferComponents::UvwMat,
1049 : VisBufferComponents::Weight,
1050 0 : -1);
1051 :
1052 0 : Bool addCorrectedVisCube = !(rvi_p->msColumns().correctedData().isNull());
1053 :
1054 0 : if (addCorrectedVisCube){
1055 0 : prefetchColumns.insert (VisBufferComponents::CorrectedCube);
1056 : // This can cause an error if a multi-MS has a mixture of MSs with corrected
1057 : // and without corrected data columns.
1058 : }
1059 :
1060 : // Replace the existing VIs with an async version. Keep pointers to the old
1061 : // ones around (these are kept by the caller) so they can be swapped back afterwards.
1062 :
1063 0 : vb_p->detachFromVisIter();
1064 0 : oldRvi = rvi_p;
1065 :
1066 0 : if (wvi_p != NULL){
1067 0 : oldWvi = wvi_p;
1068 0 : wvi_p = new VisibilityIterator (& prefetchColumns, * wvi_p);
1069 0 : rvi_p = wvi_p;
1070 : }
1071 : else{
1072 0 : rvi_p = new ROVisibilityIterator (& prefetchColumns, * rvi_p);
1073 : }
1074 0 : }
1075 :
1076 26 : void CubeSkyEquation::isLargeCube(ImageInterface<Complex>& theIm,
1077 : Int& nslice) {
1078 :
1079 : //non-cube
1080 26 : if(theIm.shape()[3]==1){
1081 22 : nslice=1;
1082 : }
1083 : else{
1084 8 : LogIO os(LogOrigin("CubeSkyEquation", "isLargeCube"));
1085 4 : Long npix=theIm.shape().product();
1086 : // use memory size denfined in aisprc if exists
1087 4 : Long memtot=HostInfo::memoryTotal(true); // Use aipsrc/casarc
1088 : //check for 32 bit OS and limit it to 2Gbyte
1089 : if( sizeof(void*) == 4){
1090 : if(memtot > 2000000)
1091 : memtot=2000000;
1092 : }
1093 4 : if(memtot < 512000){
1094 0 : ostringstream oss;
1095 0 : oss << "The amount of memory reported " << memtot << " kB is too small to work with" << endl;
1096 0 : throw(AipsError(String(oss)));
1097 :
1098 0 : }
1099 4 : Long pixInMem=(memtot/8)*1024;
1100 : // cerr << "CSE: " << memtot << " " << pixInMem << endl;
1101 : // cerr << npix << " " << pixInMem/8 << endl;
1102 4 : nslice=1;
1103 : //There are roughly 13 float images worth held in memory
1104 : //plus some extra
1105 4 : if(npix > (pixInMem/25)){
1106 : //Lets slice it so grid is at most 1/25th of memory
1107 0 : pixInMem=pixInMem/25;
1108 : //One plane is
1109 0 : npix=theIm.shape()(0)*theIm.shape()(1)*theIm.shape()(2);
1110 0 : nchanPerSlice_p=Int(floor(pixInMem/npix));
1111 : // cerr << "Nchan " << nchanPerSlice_p << " " << pixInMem << " " << npix << " " << pixInMem/npix << endl;
1112 0 : if (nchanPerSlice_p==0){
1113 0 : nchanPerSlice_p=1;
1114 : }
1115 0 : nslice=theIm.shape()(3)/nchanPerSlice_p;
1116 0 : if( (theIm.shape()(3) % nchanPerSlice_p) > 0)
1117 0 : ++nslice;
1118 : }
1119 4 : if(nslice>1) os << "Detected a Large Cube. Making and using " << nslice << " slices/chunks." << LogIO::POST;
1120 4 : }
1121 26 : }
1122 :
1123 14 : void CubeSkyEquation::initializePutSlice(const VisBuffer& vb, Bool dopsf,
1124 : Int cubeSlice, Int nCubeSlice) {
1125 14 : AlwaysAssert(ok(),AipsError);
1126 :
1127 28 : LogIO os(LogOrigin("CubeSkyEquation", "initializePutSlice"));
1128 :
1129 14 : Bool newFTM=false;
1130 14 : newFTM = isNewFTM(&(*ftm_p[0]));
1131 14 : if (newFTM) newInitializePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1132 14 : else oldInitializePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1133 :
1134 : // newInitializePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1135 14 : }
1136 :
1137 :
1138 14 : void CubeSkyEquation::oldInitializePutSlice(const VisBuffer& vb, Bool /*dopsf*/,
1139 : Int cubeSlice, Int nCubeSlice) {
1140 14 : AlwaysAssert(ok(),AipsError);
1141 14 : Bool dirDep= (ej_ != NULL);
1142 28 : for(Int model=0; model < (sm_->numberOfModels()) ; ++model){
1143 14 : sliceCube(imPutSlice_p[model], model, cubeSlice, nCubeSlice, 0);
1144 14 : weightSlice_p[model].resize();
1145 14 : if(nCubeSlice>1){
1146 0 : iftm_p[model]->reset();
1147 : }
1148 14 : iftm_p[model]->initializeToSky(*(imPutSlice_p[model]),weightSlice_p[model],
1149 : vb);
1150 14 : dirDep= dirDep || (ftm_p[model]->name() == "MosaicFT");
1151 : }
1152 14 : assertSkyJones(vb, -1);
1153 : //vb_p is used to finalize things if vb has changed propoerties
1154 14 : vb_p->assign(vb, false);
1155 14 : vb_p->updateCoordInfo(& vb, dirDep);
1156 14 : }
1157 :
1158 0 : void CubeSkyEquation::newInitializePutSlice(const VisBuffer& vb, Bool dopsf,
1159 : Int cubeSlice, Int nCubeSlice) {
1160 0 : AlwaysAssert(ok(),AipsError);
1161 0 : Bool dirDep= (ej_ != NULL);
1162 :
1163 0 : Int ntaylors=sm_->numberOfTaylorTerms(),
1164 0 : nfields = sm_->numberOfModels()/ntaylors;
1165 :
1166 0 : for(Int field=0; field < nfields ; ++field){
1167 :
1168 0 : Int ntaylors = sm_->numberOfTaylorTerms();
1169 0 : if(dopsf) ntaylors = 2 * sm_->numberOfTaylorTerms() - 1;
1170 :
1171 0 : Block<CountedPtr<ImageInterface<Complex> > > imPutSliceVec(ntaylors);
1172 0 : Block<Matrix<Float> > weightSliceVec(ntaylors);
1173 0 : for(Int taylor=0; taylor < ntaylors ; ++taylor)
1174 : {
1175 0 : Int model = sm_->getModelIndex(field,taylor);
1176 0 : sliceCube(imPutSlice_p[model], model, cubeSlice, nCubeSlice, 0);
1177 0 : weightSlice_p[model].resize();
1178 0 : imPutSliceVec[taylor] = imPutSlice_p[model];
1179 0 : weightSliceVec[taylor] = weightSlice_p[model];
1180 : }
1181 :
1182 0 : if(nCubeSlice>1){
1183 0 : iftm_p[field]->reset();
1184 : }
1185 : //U// cout << "CubeSkyEqn :: Calling new initializeToSky with dopsf " << dopsf << endl;
1186 0 : iftm_p[field]->initializeToSky(imPutSliceVec, weightSliceVec,vb,dopsf);
1187 0 : dirDep= dirDep || (ftm_p[field]->name() == "MosaicFT");
1188 0 : }// end of field
1189 0 : assertSkyJones(vb, -1);
1190 : //vb_p is used to finalize things if vb has changed propoerties
1191 0 : vb_p->assign(vb, false);
1192 0 : vb_p->updateCoordInfo(& vb, dirDep);
1193 0 : }
1194 :
1195 :
1196 :
1197 :
1198 0 : void CubeSkyEquation::getCoverageImage(Int model, ImageInterface<Float>& im){
1199 0 : if ((sm_->doFluxScale(model)) && (ftm_p.nelements() > uInt(model))){
1200 0 : ftm_p[model]->getFluxImage(im);
1201 : }
1202 :
1203 0 : }
1204 :
1205 0 : void CubeSkyEquation::getWeightImage(Int model, ImageInterface<Float>& im){
1206 0 : if (iftm_p.nelements() > uInt(model)){
1207 0 : Matrix<Float> weights;
1208 0 : iftm_p[model]->getWeightImage(im, weights);
1209 0 : }
1210 0 : }
1211 :
1212 : void
1213 263 : CubeSkyEquation::putSlice(VisBuffer & vb, Bool dopsf, FTMachine::Type col, Int cubeSlice, Int nCubeSlice) {
1214 :
1215 263 : AlwaysAssert(ok(),AipsError);
1216 263 : Int nRow=vb.nRow();
1217 263 : internalChangesPut_p=false; // Does this VB change inside itself?
1218 263 : firstOneChangesPut_p=false; // Has this VB changed from the previous one?
1219 1052 : if((ftm_p[0]->name() != "MosaicFT") && (ftm_p[0]->name() != "PBWProjectFT") &&
1220 1052 : (ftm_p[0]->name() != "AWProjectFT") && (ftm_p[0]->name() != "AWProjectWBFT")) {
1221 263 : changedSkyJonesLogic(vb, firstOneChangesPut_p, internalChangesPut_p);
1222 : }
1223 : //First ft machine change should be indicative
1224 : //anyways right now we are allowing only 1 ftmachine for GridBoth
1225 263 : Bool IFTChanged=iftm_p[0]->changed(vb);
1226 :
1227 :
1228 : // we might need to recompute the "sky" for every single row, but we
1229 : // avoid this if possible.
1230 :
1231 :
1232 263 : if(internalChangesPut_p || internalChangesGet_p) {
1233 :
1234 0 : if(internalChangesGet_p)
1235 0 : internalChangesGet_p=false;
1236 :
1237 : // Yes there are changes: go row by row.
1238 :
1239 0 : for (Int row=0; row<nRow; row++) {
1240 :
1241 0 : if(IFTChanged||changedSkyJones(vb,row)) {
1242 : // Need to apply the SkyJones from the previous row
1243 : // and finish off before starting with this row
1244 0 : finalizePutSlice(* vb_p, dopsf, cubeSlice, nCubeSlice);
1245 0 : initializePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1246 : }
1247 :
1248 0 : for (Int model=0; model<sm_->numberOfModels()/( isNewFTM()? sm_->numberOfTaylorTerms() : 1 ); ++model){
1249 0 : iftm_p[model]->put(vb, row, dopsf, col);
1250 : }
1251 : }
1252 0 : }
1253 263 : else if (IFTChanged || firstOneChangesPut_p || firstOneChangesGet_p) {
1254 :
1255 0 : if(firstOneChangesGet_p)
1256 0 : firstOneChangesGet_p=false;
1257 :
1258 0 : if(!isBeginingOfSkyJonesCache_p){
1259 0 : finalizePutSlice(*vb_p, dopsf, cubeSlice, nCubeSlice);
1260 : }
1261 0 : initializePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1262 0 : isBeginingOfSkyJonesCache_p=false;
1263 0 : for (Int model=0; model<sm_->numberOfModels()/( isNewFTM()? sm_->numberOfTaylorTerms() : 1 ); ++model){
1264 0 : iftm_p[model]->put(vb, -1, dopsf, col);
1265 : }
1266 0 : }
1267 : else {
1268 526 : for (Int model=0; model<sm_->numberOfModels()/( isNewFTM()? sm_->numberOfTaylorTerms() : 1 ); ++model){
1269 263 : iftm_p[model]->put(vb, -1, dopsf, col);
1270 : }
1271 : }
1272 :
1273 263 : isBeginingOfSkyJonesCache_p=false;
1274 :
1275 263 : }
1276 :
1277 862 : Bool CubeSkyEquation::isNewFTM(FTMachine* ftm)
1278 : {
1279 : return (
1280 862 : (ftm->name() == "AWProjectFT")
1281 1724 : || (ftm->name() == "AWProjectWBFT")
1282 1724 : || (ftm->name() == "PBWProjectFT")
1283 2586 : || (ftm->name() == "NewMultiTermFT")
1284 : // || (ftm->name() == "GridFT")
1285 1724 : );
1286 : }
1287 :
1288 822 : Bool CubeSkyEquation::isNewFTM()
1289 : {
1290 822 : AlwaysAssert( &(*ftm_p[0]) ,AipsError );
1291 822 : return isNewFTM(&(*ftm_p[0]));
1292 : }
1293 :
1294 14 : void CubeSkyEquation::finalizePutSlice(const VisBuffer& vb, Bool dopsf,
1295 : Int cubeSlice, Int nCubeSlice)
1296 : {
1297 : //============================================================================
1298 : // NEW CODE BEGINS
1299 : // Iterate across fields
1300 28 : LogIO os(LogOrigin("CubeSkyEquation", "finalizePutSlice"));
1301 :
1302 14 : Bool newFTM=false;
1303 : /*
1304 : for (Int field=0; field < sm_->numberOfModels(); ++field)
1305 : {
1306 : newFTM = isNewFTM(&(*ftm_p[field]));
1307 :
1308 : if (newFTM) newFinalizePutSlice(vb, dopsf, cubeSlice, nCubeSlice, field);
1309 : else oldFinalizePutSlice(vb, dopsf, cubeSlice, nCubeSlice, field);
1310 : }
1311 : */
1312 14 : newFTM = isNewFTM(&(*ftm_p[0]));
1313 :
1314 14 : if (newFTM) newFinalizePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1315 14 : else oldFinalizePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1316 :
1317 : //newFinalizePutSlice(vb, dopsf, cubeSlice, nCubeSlice);
1318 : // if (!newFTM)
1319 : // tmpWBNormalizeImage(dopsf);
1320 14 : }
1321 :
1322 0 : void CubeSkyEquation::newFinalizePutSlice(const VisBuffer& vb, Bool dopsf,
1323 : Int cubeSlice, Int nCubeSlice)
1324 : {
1325 : //============================================================================
1326 : // NEW CODE BEGINS
1327 : // Iterate across fields
1328 0 : LogIO os(LogOrigin("CubeSkyEquation", "newFinalizePutSlice"));
1329 : {
1330 : // os << "WARNING!!! NEW R&D CODE IN USE...DISABLE IT BEFORE CHECK-IN "
1331 : // << ftm_p[field]->name()
1332 : // << LogIO::WARN << LogIO::POST;
1333 :
1334 0 : sm_->setImageNormalization(true);
1335 0 : for (Int field=0; field < sm_->numberOfModels()/sm_->numberOfTaylorTerms(); ++field)
1336 : {
1337 0 : ft_=&(*ftm_p[field]);
1338 0 : ift_=&(*iftm_p[field]);
1339 :
1340 : // Number of Taylor terms per field
1341 0 : Int ntaylors = sm_->numberOfTaylorTerms();
1342 0 : if(dopsf) ntaylors = 2 * sm_->numberOfTaylorTerms() - 1;
1343 :
1344 : // Build a list of reference images to send into FTMachine
1345 0 : PtrBlock<SubImage<Float> *> gSSliceVec(ntaylors);
1346 0 : PtrBlock<SubImage<Float> *> ggSSliceVec(ntaylors);
1347 0 : PtrBlock<SubImage<Float> *> fluxScaleVec(ntaylors);
1348 0 : Block<CountedPtr<ImageInterface<Complex> > > imPutSliceVec(ntaylors);
1349 0 : Block<Matrix<Float> > weightSliceVec(ntaylors); // this is by value
1350 0 : for (Int taylor=0; taylor < ntaylors; ++taylor)
1351 : {
1352 0 : Int model = sm_->getModelIndex(field,taylor);
1353 0 : sliceCube(gSSliceVec[taylor], sm_->gS(model), cubeSlice, nCubeSlice);
1354 0 : sliceCube(ggSSliceVec[taylor], sm_->ggS(model), cubeSlice, nCubeSlice);
1355 0 : sliceCube(fluxScaleVec[taylor], sm_->fluxScale(model), cubeSlice, nCubeSlice);
1356 0 : imPutSliceVec[taylor] = imPutSlice_p[model];
1357 0 : weightSliceVec[taylor] = weightSlice_p[model];
1358 : }// end of taylor
1359 :
1360 : // Call finalizeToSky for this field.
1361 : // -- calls getImage, getWeightImage, does Stokes conversion, and gS/ggS normalization
1362 : //U// cout << "CubeSkyEqn :: calling new finalizeToSky with dopsf " << dopsf << endl;
1363 0 : iftm_p[field]->finalizeToSky( imPutSliceVec , gSSliceVec , ggSSliceVec , fluxScaleVec, dopsf , weightSliceVec, vb );
1364 :
1365 : // Clean up temporary reference images
1366 0 : for (Int taylor=0; taylor < ntaylors; ++taylor)
1367 : {
1368 0 : Int model = sm_->getModelIndex(field,taylor);
1369 0 : weightSlice_p[model] = weightSliceVec[taylor]; // because this is by value...
1370 0 : delete gSSliceVec[taylor];
1371 0 : delete ggSSliceVec[taylor];
1372 0 : delete fluxScaleVec[taylor];
1373 : }
1374 0 : }// end of field
1375 :
1376 : // storeImg(String("stokesNormed1.im"), *(gSSliceVec[0]));
1377 0 : tmpWBNormalizeImage(dopsf,ft_->getPBLimit());
1378 : // storeImg(String("stokesNormed2.im"), *(gSSliceVec[0]));
1379 :
1380 0 : ft_=&(*ftm_p[0]);
1381 0 : ift_=&(*iftm_p[0]);
1382 : // 4. Finally, we add the statistics
1383 0 : sm_->addStatistics(sumwt, chisq);
1384 : }
1385 0 : }
1386 :
1387 14 : void CubeSkyEquation::oldFinalizePutSlice(const VisBuffer& vb, Bool /*dopsf*/,
1388 : Int cubeSlice, Int nCubeSlice)
1389 : {
1390 : // cerr << "### Using old code: " << ftm_p[model]->name() << endl;
1391 : {
1392 28 : for (Int model=0; model < sm_->numberOfModels(); ++model)
1393 : {
1394 : //the different apply...jones use ft_ and ift_
1395 14 : ft_=&(*ftm_p[model]);
1396 14 : ift_=&(*iftm_p[model]);
1397 : // Actually do the transform. Update weights as we do so.
1398 14 : iftm_p[model]->finalizeToSky();
1399 : // 1. Now get the (unnormalized) image and add the
1400 : // weight to the summed weight
1401 14 : Matrix<Float> delta;
1402 : //imPutSlice_p[model]->copyData(iftm_p[model]->getImage(delta, false));
1403 14 : iftm_p[model]->getImage(delta, false);
1404 : //iftm_p[field]->finalizeToSky( imPutSliceVec , gSSliceVec , ggSSliceVec , fluxScaleVec, dopsf , weightSliceVec );
1405 14 : weightSlice_p[model]+=delta;
1406 :
1407 : // 2. Apply the SkyJones and add to grad chisquared
1408 : SubImage<Float> *workSlice;
1409 : SubImage<Float> *gSSlice;
1410 14 : sliceCube(workSlice, sm_->work(model), cubeSlice, nCubeSlice);
1411 14 : sliceCube(gSSlice, sm_->gS(model), cubeSlice, nCubeSlice);
1412 :
1413 14 : applySkyJonesInv(vb, -1, *(imPutSlice_p[model]), *workSlice,
1414 14 : *gSSlice);
1415 : SubImage<Float> *ggSSlice;
1416 14 : sliceCube(ggSSlice, sm_->ggS(model), cubeSlice, nCubeSlice);
1417 :
1418 : // 3. Apply the square of the SkyJones and add this to gradgrad chisquared
1419 14 : applySkyJonesSquare(vb, -1, weightSlice_p[model], *workSlice,
1420 14 : *ggSSlice);
1421 :
1422 14 : (imPutSlice_p[model])->clearCache();
1423 : //imPutSlice_p[model]->clearCache();
1424 14 : delete workSlice;
1425 14 : delete gSSlice;
1426 14 : delete ggSSlice;
1427 14 : }
1428 14 : ft_=&(*ftm_p[0]);
1429 14 : ift_=&(*iftm_p[0]);
1430 : // 4. Finally, we add the statistics
1431 14 : sm_->addStatistics(sumwt, chisq);
1432 : }
1433 14 : sm_->setImageNormalization(false);
1434 14 : }
1435 :
1436 12 : void CubeSkyEquation::initializeGetSlice(const VisBuffer& vb,
1437 : Int row,
1438 : Bool incremental, Int cubeSlice,
1439 : Int nCubeSlice)
1440 : {
1441 24 : LogIO os(LogOrigin("CubeSkyEquation", "initializeGetSlice"));
1442 :
1443 : // oldInitializeGetSlice(vb, row, incremental, cubeSlice, nCubeSlice);
1444 :
1445 12 : Bool newFTM=false;
1446 : /*
1447 : for (Int field=0; field < sm_->numberOfModels(); ++field)
1448 : {
1449 : newFTM = isNewFTM(&(*ftm_p[field]));
1450 :
1451 : if (newFTM) newInitializeGetSlice(vb, row, incremental,cubeSlice, nCubeSlice);
1452 : else oldInitializeGetSlice(vb, row, incremental,cubeSlice, nCubeSlice);
1453 : }*/
1454 :
1455 12 : newFTM = isNewFTM(&(*ftm_p[0]));
1456 :
1457 12 : if (newFTM) newInitializeGetSlice(vb, row, incremental,cubeSlice, nCubeSlice);
1458 12 : else oldInitializeGetSlice(vb, row, incremental,cubeSlice, nCubeSlice);
1459 :
1460 : //newInitializeGetSlice(vb, row, incremental,cubeSlice, nCubeSlice);
1461 : // if (!newFTM)
1462 : // tmpWBNormalizeImage(dopsf);
1463 12 : }
1464 :
1465 0 : void CubeSkyEquation::newInitializeGetSlice(const VisBuffer& vb,
1466 : Int /*row*/,
1467 : Bool incremental, Int cubeSlice,
1468 : Int nCubeSlice)
1469 : {
1470 : // imGetSlice_p.resize(sm_->numberOfModels(), true, false);
1471 : // for(Int field=0; field < sm_->numberOfFields(); ++field){
1472 0 : sm_->setImageNormalization(true);
1473 0 : imGetSlice_p.resize(sm_->numberOfModels(), true, false);
1474 0 : for(Int model=0; model < sm_->numberOfModels()/sm_->numberOfTaylorTerms(); ++model)
1475 : {
1476 0 : if(nCubeSlice>1)
1477 0 : ftm_p[model]->reset();
1478 : }
1479 :
1480 0 : Int ntaylors=sm_->numberOfTaylorTerms(),
1481 0 : nfields = sm_->numberOfModels()/ntaylors;
1482 :
1483 0 : for(Int field=0; field < nfields; ++field)
1484 : {
1485 : //the different apply...jones user ft_ and ift_
1486 0 : ft_=&(*ftm_p[field]);
1487 0 : ift_=&(*iftm_p[field]);
1488 :
1489 0 : Block<CountedPtr<ImageInterface<Complex> > > imGetSliceVec(sm_->numberOfTaylorTerms());
1490 0 : PtrBlock<SubImage<Float> *> modelSliceVec(sm_->numberOfTaylorTerms());
1491 0 : PtrBlock<SubImage<Float> *> weightSliceVec(sm_->numberOfTaylorTerms());
1492 0 : PtrBlock<SubImage<Float> *> fluxScaleVec(sm_->numberOfTaylorTerms());
1493 0 : Block<Matrix<Float> > weightVec(sm_->numberOfTaylorTerms()); // this is by value
1494 :
1495 0 : for(Int taylor=0; taylor < sm_->numberOfTaylorTerms(); ++taylor)
1496 : {
1497 0 : Int model = sm_->getModelIndex(field,taylor);
1498 : // NEW : Do the applySkyJones slice-by-slice -- to make it go into initializeToVis :(
1499 : ///cerr << "Taylor, Model, Field: " << taylor << " " << model << " " << field << endl;
1500 0 : if(incremental)
1501 0 : sliceCube(modelSliceVec[taylor], sm_->deltaImage(model), cubeSlice, nCubeSlice);
1502 : else
1503 0 : sliceCube(modelSliceVec[taylor], sm_->image(model), cubeSlice, nCubeSlice);
1504 :
1505 0 : sliceCube(fluxScaleVec[taylor], sm_->fluxScale(model), cubeSlice, nCubeSlice);
1506 0 : sliceCube(weightSliceVec[taylor], sm_->ggS(model), cubeSlice, nCubeSlice);
1507 0 : sliceCube(imGetSlice_p[model], model, cubeSlice, nCubeSlice, 1);
1508 0 : imGetSliceVec[taylor] = imGetSlice_p[model];
1509 0 : weightVec[taylor] = weightSlice_p[model];
1510 : }// end of taylor
1511 :
1512 : //U// cout << "CubeSkyEquation :: Calling new initializeToVis with " << modelSliceVec.nelements() << " models and " << imGetSliceVec.nelements() << " complex grids " << endl;
1513 : //U// LatticeExprNode LEN = max( *(modelSliceVec[0] ) );
1514 : //U// cout << "CubeSkyEq : Peak in image to be predicted : " << LEN.getFloat() << endl;
1515 :
1516 0 : ftm_p[field]->initializeToVis(imGetSliceVec, modelSliceVec, weightSliceVec,
1517 : fluxScaleVec, weightVec,vb);
1518 :
1519 0 : for (Int taylor=0; taylor < sm_->numberOfTaylorTerms(); ++taylor)
1520 : {
1521 : // Int model = sm_->getModelIndex(field,taylor);
1522 : // weightSlice_p[model] = weightSliceVec[taylor]; // because this is by value...
1523 0 : delete modelSliceVec[taylor];
1524 0 : delete weightSliceVec[taylor];
1525 0 : delete fluxScaleVec[taylor];
1526 : }
1527 0 : }//end of field
1528 0 : ft_=&(*ftm_p[0]);
1529 0 : ift_=&(*iftm_p[0]);
1530 0 : }
1531 :
1532 12 : void CubeSkyEquation::oldInitializeGetSlice(const VisBuffer& vb,
1533 : Int row,
1534 : Bool incremental, Int cubeSlice,
1535 : Int nCubeSlice){
1536 12 : sm_->setImageNormalization(false);
1537 :
1538 12 : imGetSlice_p.resize(sm_->numberOfModels(), true, false);
1539 24 : for(Int model=0; model < sm_->numberOfModels(); ++model){
1540 12 : if(nCubeSlice>1){
1541 0 : ftm_p[model]->reset();
1542 : }
1543 : //the different apply...jones user ft_ and ift_
1544 12 : ft_=&(*ftm_p[model]);
1545 12 : ift_=&(*iftm_p[model]);
1546 12 : if(cubeSlice==0){
1547 12 : if(incremental) {
1548 1 : applySkyJones(vb, row, sm_->deltaImage(model), sm_->cImage(model));
1549 : }
1550 : else {
1551 11 : applySkyJones(vb, row, sm_->image(model), sm_->cImage(model));
1552 : }
1553 : }
1554 12 : sliceCube(imGetSlice_p[model], model, cubeSlice, nCubeSlice, 1);
1555 12 : ftm_p[model]->initializeToVis(*(imGetSlice_p[model]), vb);
1556 : }
1557 12 : ft_=&(*ftm_p[0]);
1558 12 : ift_=&(*iftm_p[0]);
1559 :
1560 :
1561 12 : }
1562 :
1563 26 : void CubeSkyEquation::sliceCube(CountedPtr<ImageInterface<Complex> >& slice,Int model, Int cubeSlice,
1564 : Int nCubeSlice, Int typeOfSlice){
1565 :
1566 26 : IPosition blc(4,0,0,0,0);
1567 52 : IPosition trc(4,sm_->cImage(model).shape()(0)-1,
1568 78 : sm_->cImage(model).shape()(1)-1,sm_->cImage(model).shape()(2)-1,
1569 78 : 0);
1570 26 : Int beginChannel=cubeSlice*nchanPerSlice_p;
1571 26 : Int endChannel=beginChannel+nchanPerSlice_p-1;
1572 26 : if(cubeSlice==(nCubeSlice-1))
1573 26 : endChannel=sm_->image(model).shape()(3)-1;
1574 26 : blc(3)=beginChannel;
1575 26 : trc(3)=endChannel;
1576 26 : sl_p=Slicer (blc, trc, Slicer::endIsLast);
1577 26 : SubImage<Complex>* sliceIm= new SubImage<Complex>(sm_->cImage(model), sl_p, true); /// UUU changes to true
1578 : // cerr << "SliceCube: " << beginChannel << " " << endChannel << endl;
1579 26 : if(typeOfSlice==0){
1580 :
1581 14 : Double memoryMB=HostInfo::memoryFree()/1024.0/(5.0*(sm_->numberOfModels()));
1582 28 : slice=new TempImage<Complex> (TiledShape(sliceIm->shape(),
1583 42 : IPosition(4, min(sliceIm->shape()(0)/4, 1000), min(sliceIm->shape()(1)/4, 1000),sliceIm->shape()(2) , 1)), sliceIm->coordinates(), sm_->getMemoryUse() ? memoryMB: 0);
1584 :
1585 :
1586 : //slice->setMaximumCacheSize((sliceIm->shape()[0])*(sliceIm->shape()[1])/4);
1587 14 : slice->setMaximumCacheSize(sliceIm->shape().product());
1588 : //slice.copyData(sliceIm);
1589 14 : slice->set(Complex(0.0, 0.0));
1590 : //slice->setCoordinateInfo(sm_->image(model).coordinates());
1591 14 : delete sliceIm;
1592 : }
1593 : else{
1594 12 : slice=sliceIm;
1595 : }
1596 :
1597 26 : }
1598 :
1599 42 : void CubeSkyEquation::sliceCube(SubImage<Float>*& slice,
1600 : ImageInterface<Float>& image, Int cubeSlice,
1601 : Int nCubeSlice){
1602 :
1603 42 : IPosition blc(4,0,0,0,0);
1604 84 : IPosition trc(4,image.shape()(0)-1,
1605 126 : image.shape()(1)-1,image.shape()(2)-1,
1606 126 : 0);
1607 42 : Int beginChannel=cubeSlice*nchanPerSlice_p;
1608 42 : Int endChannel=beginChannel+nchanPerSlice_p-1;
1609 42 : if(cubeSlice==(nCubeSlice-1))
1610 42 : endChannel=image.shape()(3)-1;
1611 42 : blc(3)=beginChannel;
1612 42 : trc(3)=endChannel;
1613 42 : sl_p=Slicer(blc, trc, Slicer::endIsLast);
1614 : //writeable if possible
1615 42 : slice= new SubImage<Float> (image, sl_p, true);
1616 42 : }
1617 :
1618 268 : VisBuffer& CubeSkyEquation::getSlice(VisBuffer& result,
1619 : Bool incremental,
1620 : Int cubeSlice, Int nCubeSlice) {
1621 :
1622 268 : Int nRow=result.nRow();
1623 :
1624 268 : result.modelVisCube(); // get the visibility so vb will have it
1625 268 : VisBuffer vb(result); // method only called using writable VI so no ROVIA
1626 :
1627 268 : Int nmodels=sm_->numberOfModels()/( isNewFTM()? sm_->numberOfTaylorTerms() : 1 );
1628 268 : Bool FTChanged=ftm_p[0]->changed(vb);
1629 :
1630 : // we might need to recompute the "sky" for every single row, but we
1631 : // avoid this if possible.
1632 268 : internalChangesGet_p=false; // Does this VB change inside itself?
1633 268 : firstOneChangesGet_p=false; // Has this VB changed from the previous one?
1634 1072 : if((ftm_p[0]->name() != "MosaicFT") && (ftm_p[0]->name() != "PBWProjectFT") &&
1635 1072 : (ftm_p[0]->name() != "AWProjectFT") && (ftm_p[0]->name() != "AWProjectWBFT")) {
1636 268 : changedSkyJonesLogic(result, firstOneChangesGet_p, internalChangesGet_p);
1637 : }
1638 :
1639 268 : if(internalChangesGet_p || internalChangesPut_p) {
1640 0 : if(internalChangesPut_p)
1641 0 : internalChangesPut_p=false;
1642 : // Yes there are changes within this buffer: go row by row.
1643 : // This will automatically catch a change in the FTMachine so
1644 : // we don't have to check for that.
1645 :
1646 0 : Matrix<Complex> refres;
1647 0 : Matrix<Complex> refvb;
1648 0 : for (Int row=0; row<nRow; row++) {
1649 0 : finalizeGetSlice();
1650 0 : initializeGetSlice(result, row, false, cubeSlice,
1651 : nCubeSlice);
1652 0 : if(incremental || (nmodels > 1) ){
1653 0 : for (Int model=0; model < nmodels; ++model){
1654 0 : ftm_p[model]->get(vb,row);
1655 0 : refvb.reference(vb.modelVisCube().xyPlane(row));
1656 0 : refres.reference(result.modelVisCube().xyPlane(row));
1657 0 : refres += refvb;
1658 : }
1659 0 : }
1660 : else
1661 0 : ftm_p[0]->get(result, row);
1662 : }
1663 0 : }
1664 268 : else if (FTChanged || firstOneChangesGet_p || firstOneChangesPut_p) {
1665 0 : if(firstOneChangesPut_p)
1666 0 : firstOneChangesPut_p=false;
1667 : // This buffer has changed wrt the previous buffer, but
1668 : // this buffer has no changes within it. Again we don't need to
1669 : // check for the FTMachine changing.
1670 :
1671 0 : finalizeGetSlice();
1672 0 : initializeGetSlice(result, 0, false, cubeSlice, nCubeSlice);
1673 0 : if(incremental || (nmodels > 1) ){
1674 0 : for (Int model=0; model < nmodels; ++model){
1675 0 : ftm_p[model]->get(vb);
1676 0 : result.modelVisCube()+=vb.modelVisCube();
1677 : }
1678 0 : }
1679 : else
1680 0 : ftm_p[0]->get(result);
1681 0 : }
1682 : else {
1683 268 : if(incremental || (nmodels >1) ){
1684 36 : for (Int model=0; model < nmodels; ++model){
1685 18 : ftm_p[model]->get(vb);
1686 18 : result.modelVisCube()+=vb.modelVisCube();
1687 : }
1688 18 : }
1689 : else
1690 250 : ftm_p[0]->get(result);
1691 : }
1692 268 : return result;
1693 :
1694 268 : }
1695 :
1696 : void
1697 12 : CubeSkyEquation::finalizeGetSlice(){
1698 : //// place-holders.... there is nothing to do after degridding
1699 :
1700 12 : if( ftm_p[0]->name() != "NewMultiTermFT" ) // Added this filter to prevent seg-fault in widebandmosaic regression that started after r30978.
1701 : {
1702 24 : for (Int model=0; model < sm_->numberOfModels(); ++model)
1703 12 : ftm_p[model]->finalizeToVis();
1704 : }
1705 12 : }
1706 :
1707 : Bool
1708 26 : CubeSkyEquation::getFreqRange(ROVisibilityIterator& vi,
1709 : const CoordinateSystem& coords,
1710 : Int slice, Int nslice){
1711 : //bypass this for now
1712 26 : return false;
1713 : // Enforce that all SPWs are in the same frequency frame.
1714 : //
1715 : // If all the SPWs in the MS are in LSRK frame, we can do data
1716 : // selection (since image is always in LSRK).
1717 : //
1718 : // If not all SPWs in the MS are in the same frequency frame and
1719 : // in LSRK frame, for now, disable data selection since the
1720 : // mapping between image (in LSRK) and MS channels will be time
1721 : // variable.
1722 : VisBufferAutoPtr vb (vi);
1723 : ROScalarMeasColumn<MFrequency> freqFrame=vb->msColumns().spectralWindow().refFrequencyMeas();
1724 : uInt nrows=vb->msColumns().spectralWindow().nrow();
1725 : String firstString = freqFrame(0).getRefString();
1726 : Bool allFramesSame=true;
1727 : for (uInt i=0;i<nrows;i++)
1728 : if (freqFrame(i).getRefString() != firstString)
1729 : {allFramesSame = false;break;}
1730 :
1731 : if (!allFramesSame || (firstString!="LSRK"))
1732 : return false;
1733 :
1734 : // Only one slice lets keep what the user selected
1735 : if(nslice==1)
1736 : return false;
1737 :
1738 : Double start=0.0;
1739 : Double end=0.0;
1740 : Double chanwidth=1.0;
1741 : Int specIndex=coords.findCoordinate(Coordinate::SPECTRAL);
1742 : SpectralCoordinate specCoord=coords.spectralCoordinate(specIndex);
1743 : Vector<Int>spectralPixelAxis=coords.pixelAxes(specIndex);
1744 : if(nchanPerSlice_p>0){
1745 : specCoord.toWorld(start,Double(slice*nchanPerSlice_p)-0.5);
1746 : specCoord.toWorld(end, Double(nchanPerSlice_p*(slice+1))+0.5);
1747 : chanwidth=fabs(end-start)/Double(nchanPerSlice_p);
1748 : }
1749 : if(end < start){
1750 : Double tempoo=start;
1751 : start=end;
1752 : end=tempoo;
1753 : }
1754 :
1755 : Block<Vector<Int> > spwb;
1756 : Block<Vector<Int> > startb;
1757 : Block<Vector<Int> > nchanb;
1758 : Block<Vector<Int> > incrb=blockChanInc_p;
1759 : vi.getSpwInFreqRange(spwb, startb, nchanb, start, end, chanwidth);
1760 : if(spwb.nelements()==0)
1761 : return false;
1762 :
1763 : //cerr << "Original is " << blockChanStart_p[0] << " " << blockChanWidth_p[0] << " " << blockChanInc_p[0] << " "
1764 : // << blockSpw_p[0] << endl;
1765 : //vi.selectChannel(1, startb[0][0], nchanb[0][0], 1, spwb[0][0]);
1766 : vi.selectChannel(blockNumChanGroup_p, startb, nchanb, incrb, spwb);
1767 :
1768 : return true;
1769 :
1770 : }
1771 :
1772 14 : void CubeSkyEquation::fixImageScale()
1773 : {
1774 28 : LogIO os(LogOrigin("CubeSkyEquation", "fixImageScale"));
1775 :
1776 : // make a minimum value to ggS
1777 : // This has the same effect as Sault Weighting, but
1778 : // is implemented somewhat differently.
1779 : // We also keep the fluxScale(mod) images around to
1780 : // undo the weighting.
1781 :
1782 28 : for (Int model=0;model<sm_->numberOfModels()/( isNewFTM()? sm_->numberOfTaylorTerms() : 1 );model++) {
1783 :
1784 14 : Float ggSMax=0.0;
1785 : {
1786 :
1787 14 : LatticeExprNode LEN = max( sm_->ggS(model) );
1788 14 : ggSMax = max(ggSMax,LEN.getFloat());
1789 14 : }
1790 : Float ggSMin1;
1791 : Float ggSMin2;
1792 :
1793 :
1794 14 : ggSMin1 = ggSMax * constPB_p * constPB_p;
1795 14 : ggSMin2 = ggSMax * minPB_p * minPB_p;
1796 :
1797 :
1798 14 : if(ej_ || (ftm_p[model]->name() == "MosaicFT") ) {
1799 :
1800 :
1801 :
1802 : /*Don't print this for now
1803 : if (scaleType_p == "SAULT") {
1804 : os << "Using SAULT image plane weighting" << LogIO::POST;
1805 : }
1806 : else {
1807 : os << "Using No image plane weighting" << LogIO::POST;
1808 : }
1809 : */
1810 0 : sm_->fluxScale(model).removeRegion ("mask0", RegionHandler::Any, false);
1811 0 : if ((ftm_p[model]->name()!="MosaicFT")) {
1812 0 : if(scaleType_p=="SAULT"){
1813 :
1814 : // Adjust flux scale to account for ggS being truncated at ggSMin1
1815 : // Below ggSMin2, set flux scale to 0.0
1816 : // FluxScale * image => true brightness distribution, but
1817 : // noise increases at edge.
1818 : // if ggS < ggSMin2, set to Zero;
1819 : // if ggS > ggSMin2 && < ggSMin1, set to ggSMin1/ggS
1820 : // if ggS > ggSMin1, set to 1.0
1821 :
1822 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1823 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1824 0 : sqrt((sm_->ggS(model))/ggSMin1) )) );
1825 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1826 0 : (iif(sm_->ggS(model) > (ggSMin1), 1.0,
1827 0 : (sm_->fluxScale(model)) )) );
1828 : // truncate ggS at ggSMin1
1829 0 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1830 0 : (iif(sm_->ggS(model) < (ggSMin1), ggSMin1*(sm_->fluxScale(model)),
1831 0 : sm_->ggS(model)) )
1832 : );
1833 :
1834 : }
1835 :
1836 : else{
1837 :
1838 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1839 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1840 0 : sqrt((sm_->ggS(model))/ggSMax) )) );
1841 0 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1842 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1843 0 : sqrt((sm_->ggS(model))*ggSMax) )) );
1844 :
1845 : }
1846 :
1847 : } else {
1848 :
1849 0 : Int nXX=sm_->ggS(model).shape()(0);
1850 0 : Int nYY=sm_->ggS(model).shape()(1);
1851 0 : Int npola= sm_->ggS(model).shape()(2);
1852 0 : Int nchana= sm_->ggS(model).shape()(3);
1853 0 : IPosition blc(4,nXX, nYY, npola, nchana);
1854 0 : IPosition trc(4, nXX, nYY, npola, nchana);
1855 0 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
1856 :
1857 0 : Vector<Float> planeWeightMax(npola*nchana, 0.0);
1858 : //Those damn weights per plane can be wildly different so
1859 : //deal with it properly here
1860 0 : for (Int j=0; j < npola; ++j){
1861 0 : for (Int k=0; k < nchana ; ++k){
1862 :
1863 0 : blc(2)=j; trc(2)=j;
1864 0 : blc(3)=k; trc(3)=k;
1865 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1866 0 : SubImage<Float> fscalesub(sm_->fluxScale(model), sl, true);
1867 0 : SubImage<Float> ggSSub(sm_->ggS(model), sl, true);
1868 : Float planeMax;
1869 0 : LatticeExprNode LEN = max( ggSSub );
1870 0 : planeMax = LEN.getFloat();
1871 0 : planeWeightMax(j*nchana+k)=planeMax;
1872 : ///////////
1873 0 : LatticeExprNode LEN1 = min( ggSSub );
1874 : os << LogIO::DEBUG1
1875 0 : << "Max " << planeMax << " min " << LEN1.getFloat() << LogIO::POST;
1876 : ////////////
1877 : ///As we chop the image later...the weight can vary per channel
1878 : ///lets be conservative and go to 1% of ggsMin2
1879 0 : if(planeMax !=0){
1880 : //cerr << "DOFLAT " << doflat_p << endl;
1881 0 : if(doflat_p){
1882 0 : fscalesub.copyData( (LatticeExpr<Float>)
1883 0 : (iif(ggSSub < (ggSMin2/100.0),
1884 0 : 0.0, sqrt(ggSSub/planeMax))));
1885 0 : ggSSub.copyData( (LatticeExpr<Float>)
1886 0 : (iif(ggSSub < (ggSMin2/100.0), 0.0,
1887 0 : sqrt(planeMax*ggSSub))));
1888 : }
1889 : else{
1890 0 : fscalesub.copyData( (LatticeExpr<Float>)
1891 0 : (iif(ggSSub < (ggSMin2/100.0),
1892 0 : 0.0, (ggSSub/planeMax))));
1893 0 : ggSSub.copyData( (LatticeExpr<Float>)
1894 0 : (iif(ggSSub < (ggSMin2/100.0), 0.0,
1895 : (planeMax))));
1896 : }
1897 :
1898 :
1899 :
1900 : //ggSSub.copyData( (LatticeExpr<Float>)
1901 : // (iif(ggSSub < (ggSMin2/100.0), 0.0,
1902 : // planeMax)));
1903 :
1904 :
1905 : }
1906 0 : }
1907 :
1908 : }
1909 :
1910 0 : TableRecord info=(sm_->fluxScale(model)).miscInfo();
1911 0 : info.define("weightMax", planeWeightMax);
1912 0 : sm_->fluxScale(model).setMiscInfo(info);
1913 : /*
1914 :
1915 : ftm_p[model]->getFluxImage(sm_->fluxScale(model));
1916 :
1917 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1918 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1919 : (sm_->ggS(model)/ggSMax) )) );
1920 :
1921 : */
1922 : //}
1923 0 : }
1924 :
1925 : //because for usual ft machines a applySJoneInv is done on the gS
1926 : //in the finalizepu tstage...need to understand if its necessary
1927 : /*need to understand that square business
1928 : if( (ft_->name() != "MosaicFT") && (!isPSFWork_p)){
1929 : sm_->gS(model).copyData( (LatticeExpr<Float>)
1930 : (iif(sm_->fluxScale(model) > 0.0,
1931 : ((sm_->gS(model))/(sm_->fluxScale(model))), 0.0 )) );
1932 :
1933 : }
1934 : */
1935 : ///
1936 0 : sm_->fluxScale(model).clearCache();
1937 0 : sm_->ggS(model).clearCache();
1938 : }
1939 :
1940 : }
1941 14 : }
1942 :
1943 0 : void CubeSkyEquation::tmpWBNormalizeImage(Bool& dopsf, const Float& pbLimit)
1944 : {
1945 0 : LogIO os(LogOrigin("CubeSkyEquation", "tmpNormalizeImage"));
1946 :
1947 0 : if (dopsf) return;
1948 :
1949 :
1950 : Int nCubeSlice;
1951 : // Number of Taylor terms per field
1952 0 : Int ntaylors = sm_->numberOfTaylorTerms();
1953 0 : isLargeCube(sm_->cImage(0), nCubeSlice);
1954 :
1955 : // PSFs are normalized in makeApproxPSF()
1956 0 : if(dopsf) ntaylors = 2 * sm_->numberOfTaylorTerms() - 1;
1957 :
1958 0 : Int nfields = sm_->numberOfModels()/ntaylors;
1959 :
1960 0 : for (Int cubeSlice=0; cubeSlice<nCubeSlice;cubeSlice++)
1961 : {
1962 0 : for (Int field=0; field<nfields; field++)
1963 : {
1964 0 : Int baseindex = sm_->getModelIndex(field,0); // field,taylorterm
1965 :
1966 : SubImage<Float> *ggSSliceVec;
1967 0 : sliceCube(ggSSliceVec, sm_->ggS(baseindex), cubeSlice, nCubeSlice);
1968 :
1969 0 : for (Int taylor=0; taylor < ntaylors; ++taylor)
1970 : {
1971 0 : Int index = sm_->getModelIndex(field, taylor);
1972 :
1973 : SubImage<Float> *gSSliceVec;
1974 0 : sliceCube(gSSliceVec, sm_->gS(index), cubeSlice, nCubeSlice);
1975 :
1976 : //
1977 : // If the FTM is NewMultiTermFT and is configure to not
1978 : // apply PB corrections, don't apply the PB correction
1979 : // here either.
1980 : //
1981 0 : LatticeExpr<Float> le;
1982 0 : if ((ft_->name()=="NewMultiTermFT"))
1983 : {
1984 0 : if (((NewMultiTermFT *)ft_)->getDOPBCorrection())
1985 : {
1986 : ////// PBSQWeight
1987 0 : le=LatticeExpr<Float>(iif((*ggSSliceVec)>(pbLimit), (*gSSliceVec)/(sqrt(*ggSSliceVec)), 0.0)); // The negative sign is in FTM::normalizeImage()
1988 : ///// PBWeight
1989 : //le=LatticeExpr<Float>(iif((*ggSSliceVec)>(pbLimit), (*gSSliceVec)/((*ggSSliceVec)), 0.0)); // The negative sign is in FTM::normalizeImage()
1990 0 : gSSliceVec->copyData(le);
1991 : }
1992 : }
1993 : else
1994 : {
1995 : ////// PBSQWeight
1996 0 : le=LatticeExpr<Float>(iif((*ggSSliceVec)>(pbLimit), (*gSSliceVec)/(sqrt(*ggSSliceVec)), 0.0)); // The negative sign is in FTM::normalizeImage()
1997 : ////// PBWeight
1998 : //le=LatticeExpr<Float>(iif((*ggSSliceVec)>(pbLimit), (*gSSliceVec)/((*ggSSliceVec)), 0.0)); // The negative sign is in FTM::normalizeImage()
1999 0 : gSSliceVec->copyData(le);
2000 : }
2001 :
2002 : // if (dopsf)
2003 : // {
2004 : // storeImg(String("thePSF.im"), *gSSliceVec);
2005 : // storeImg(String("thePB.im"), *ggSSliceVec);
2006 : // }
2007 0 : delete gSSliceVec;
2008 :
2009 0 : }
2010 0 : delete ggSSliceVec;
2011 :
2012 : }
2013 : }
2014 :
2015 0 : }
2016 : } //# NAMESPACE CASA - END
|
