Line data Source code
1 : //# SkyEquation.cc: Implementation of Sky Equation classes
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <casacore/casa/BasicSL/Complex.h>
28 : #include <casacore/casa/Arrays/Matrix.h>
29 : #include <casacore/measures/Measures/MeasConvert.h>
30 : #include <synthesis/MeasurementEquations/SkyEquation.h>
31 : #include <casacore/images/Images/ImageInterface.h>
32 : #include <casacore/images/Images/SubImage.h>
33 : #include <casacore/images/Regions/ImageRegion.h>
34 : #include <synthesis/TransformMachines/SkyJones.h>
35 : #include <synthesis/TransformMachines/FTMachine.h>
36 :
37 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
38 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
39 :
40 : #include <components/ComponentModels/Flux.h>
41 : #include <components/ComponentModels/PointShape.h>
42 : #include <components/ComponentModels/ConstantSpectrum.h>
43 :
44 : #include <synthesis/TransformMachines/ComponentFTMachine.h>
45 : #include <synthesis/MeasurementComponents/SkyModel.h>
46 :
47 : #include <msvis/MSVis/VisSet.h>
48 : #include <synthesis/TransformMachines/StokesImageUtil.h>
49 : #include <msvis/MSVis/StokesVector.h>
50 : #include <msvis/MSVis/VisBufferUtil.h>
51 :
52 : #include <casacore/casa/Arrays/ArrayMath.h>
53 : #include <casacore/casa/Arrays/MatrixMath.h>
54 : #include <casacore/casa/Utilities/Assert.h>
55 : #include <casacore/casa/BasicSL/String.h>
56 : #include <casacore/lattices/Lattices/Lattice.h>
57 : #include <casacore/measures/Measures/UVWMachine.h>
58 : #include <casacore/lattices/Lattices/ArrayLattice.h>
59 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
60 : #include <casacore/lattices/LEL/LatticeExpr.h>
61 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
62 : #include <casacore/lattices/Lattices/LatticeIterator.h>
63 : #include <casacore/lattices/Lattices/LatticeStepper.h>
64 : #include <casacore/lattices/LRegions/LCRegion.h>
65 : #include <casacore/casa/Containers/Block.h>
66 :
67 : #include <casacore/casa/Exceptions/Error.h>
68 : #include <msvis/MSVis/VisibilityIterator.h>
69 : #include <msvis/MSVis/VisBuffer.h>
70 : ///////////////#include <msvis/MSVis/VisBufferAsync.h>
71 : #include <iostream>
72 :
73 : #include <casacore/casa/System/ProgressMeter.h>
74 :
75 : #include <memory>
76 :
77 : using namespace casacore;
78 : namespace casa { //# NAMESPACE CASA - BEGIN
79 :
80 : // ***************************************************************************
81 : // ******************** Start of public member functions ********************
82 : // ***************************************************************************
83 :
84 :
85 : //----------------------------------------------------------------------
86 0 : SkyEquation::SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft,
87 0 : FTMachine& ift)
88 0 : : sm_(&sm), vs_(&vs), ft_(&ft), ift_(&ift), cft_(0), ej_(0), dj_(0),
89 0 : tj_(0), fj_(0), iDebug_p(0), isPSFWork_p(false), noModelCol_p(false), isBeginingOfSkyJonesCache_p(true), doflat_p(false)
90 : {
91 0 : rvi_p=&(vs_->iter());
92 0 : wvi_p=&(vs_->iter());
93 :
94 0 : vb_p.set (rvi_p);
95 0 : };
96 :
97 : //----------------------------------------------------------------------
98 0 : SkyEquation::SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft)
99 0 : : sm_(&sm), vs_(&vs), ft_(&ft), ift_(&ft), cft_(0), ej_(0), dj_(0),
100 0 : tj_(0), fj_(0), iDebug_p(0), isPSFWork_p(false),noModelCol_p(false), isBeginingOfSkyJonesCache_p(true), doflat_p(false)
101 : {
102 0 : rvi_p=&(vs_->iter());
103 0 : wvi_p=&(vs_->iter());
104 :
105 :
106 0 : vb_p.set (rvi_p);
107 0 : };
108 :
109 : //----------------------------------------------------------------------
110 0 : SkyEquation::SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft,
111 0 : FTMachine& ift, ComponentFTMachine& cft)
112 0 : : sm_(&sm), vs_(&vs), ft_(&ft), ift_(&ift), cft_(&cft), ej_(0),
113 0 : dj_(0), tj_(0), fj_(0), iDebug_p(0), isPSFWork_p(false),noModelCol_p(false),isBeginingOfSkyJonesCache_p(true), doflat_p(false)
114 : {
115 0 : rvi_p=&(vs_->iter());
116 0 : wvi_p=&(vs_->iter());
117 :
118 :
119 0 : vb_p.set (rvi_p);
120 0 : };
121 :
122 : //----------------------------------------------------------------------
123 0 : SkyEquation::SkyEquation(SkyModel& sm, VisSet& vs, FTMachine& ft,
124 0 : ComponentFTMachine& cft, Bool noModelCol)
125 0 : : sm_(&sm), vs_(&vs), ft_(&ft), ift_(&ft), cft_(&cft), ej_(0),
126 0 : dj_(0), tj_(0), fj_(0), iDebug_p(0), isPSFWork_p(false),
127 0 : noModelCol_p(noModelCol),isBeginingOfSkyJonesCache_p(true),doflat_p(false)
128 : {
129 :
130 0 : rvi_p=&(vs_->iter());
131 :
132 0 : wvi_p=&(vs_->iter());
133 0 : vb_p.set (rvi_p);
134 0 : };
135 :
136 37 : SkyEquation::SkyEquation(SkyModel& sm, ROVisibilityIterator& vi, FTMachine& ft,
137 37 : ComponentFTMachine& cft, Bool noModelCol)
138 37 : : sm_(&sm), ft_(&ft), ift_(&ft), cft_(&cft), ej_(0),
139 111 : dj_(0), tj_(0), fj_(0), iDebug_p(0), isPSFWork_p(false),
140 37 : noModelCol_p(noModelCol),isBeginingOfSkyJonesCache_p(true), doflat_p(false){
141 :
142 :
143 : //visiter is read only
144 37 : rvi_p=&vi;
145 37 : if(rvi_p->ms().isWritable())
146 37 : wvi_p=static_cast<VisibilityIterator *>(&vi);
147 : else
148 0 : wvi_p=NULL;
149 :
150 37 : vb_p.set (rvi_p);
151 37 : }
152 :
153 : //----------------------------------------------------------------------
154 37 : SkyEquation::~SkyEquation() {
155 37 : }
156 :
157 : //----------------------------------------------------------------------
158 0 : SkyEquation& SkyEquation::operator=(const SkyEquation& other)
159 : {
160 0 : if(this!=&other) {
161 0 : sm_=other.sm_;
162 0 : vs_=other.vs_;
163 0 : ft_=other.ft_;
164 0 : ift_=other.ift_;
165 0 : cft_=other.cft_;
166 0 : ej_=other.ej_;
167 0 : dj_=other.dj_;
168 0 : tj_=other.tj_;
169 0 : fj_=other.fj_;
170 0 : rvi_p=other.rvi_p;
171 0 : wvi_p=other.wvi_p;
172 0 : doflat_p=other.doflat_p;
173 0 : vb_p = other.vb_p;
174 : };
175 0 : return *this;
176 : };
177 :
178 : //----------------------------------------------------------------------
179 0 : SkyEquation::SkyEquation(const SkyEquation& other)
180 : {
181 0 : operator=(other);
182 0 : }
183 :
184 : //----------------------------------------------------------------------
185 0 : void SkyEquation::setSkyJones(SkyJones& j) {
186 0 : SkyJones::Type t=j.type();
187 0 : switch(t) {
188 0 : case SkyJones::E:
189 0 : ej_=&j;
190 0 : break;
191 0 : case SkyJones::D:
192 0 : dj_=&j;
193 0 : break;
194 0 : case SkyJones::T:
195 0 : tj_=&j;
196 0 : break;
197 0 : case SkyJones::F:
198 0 : fj_=&j;
199 0 : break;
200 : }
201 0 : };
202 :
203 : //----------------------------------------------------------------------
204 : // Predict the Sky coherence
205 : //void SkyEquation::predictComponents(Bool& incremental, Bool& initialized){
206 : // predictComponents(incrementa,initialized,MS::MODEL_DATA);
207 : //}
208 :
209 37 : void SkyEquation::predictComponents(Bool& incremental, Bool& initialized, MS::PredefinedColumns Type){
210 : // Initialize
211 :
212 : // Do the component model only if this is not an incremental update;
213 37 : if(sm_->hasComponentList() && !incremental ) {
214 : //if(noModelCol_p)
215 : // throw(AipsError("Cannot deal with componentlists without using scratch columns yet"));
216 37 : if(wvi_p==NULL)
217 0 : throw(AipsError("Cannot save model in non-writable ms"));
218 37 : VisIter& vi=*wvi_p;
219 37 : checkVisIterNumRows(vi);
220 37 : VisBufferAutoPtr vb(vi);
221 :
222 : // Reset the various SkyJones
223 37 : resetSkyJones();
224 :
225 : // Loop over all visibilities
226 :
227 : // This code is currently commented out because numberCoh
228 : // called on the VI is weird relative to the VI Async semantics,
229 : // at least at first glance
230 : //
231 : // Int cohDone=0;
232 : // ProgressMeter pm(1.0, Double(vi.numberCoh()),
233 : // "Predicting component coherences",
234 : // "", "", "", true);
235 37 : Int oldmsid=-1;
236 :
237 74 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
238 5825 : for (vi.origin(); vi.more(); vi++) {
239 5788 : if(!noModelCol_p){
240 2888 : if(!incremental&&!initialized) {
241 2888 : vb->setModelVisCube(Complex(0.0,0.0));
242 : // vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
243 : }
244 :
245 :
246 : // get always fills Model
247 2888 : get(* vb, sm_->componentList() );
248 :
249 : // and write it to the model MS
250 2888 : switch(Type) {
251 2888 : case MS::MODEL_DATA:
252 2888 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
253 2888 : break;
254 0 : case MS::DATA:
255 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Observed);
256 0 : break;
257 0 : case MS::CORRECTED_DATA:
258 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Corrected);
259 0 : break;
260 0 : default:
261 0 : throw (AipsError("Programmer made a wrong call"));
262 : }
263 : }
264 : else{
265 2900 : if(vb->msId() != oldmsid){
266 19 : oldmsid=vb->msId();
267 19 : String err;
268 19 : Record clrec;
269 : //cerr << "in componentlist saving" << endl;
270 19 : if(anyGT(Int(Stokes::RR), vb->corrType())){
271 0 : throw(AipsError("Invalid MS: correlation type cannot be decomposed into feed product; DataDescId="+String::toString(vb->dataDescriptionId()) ));
272 : }
273 19 : if(!sm_->componentList().toRecord(err, clrec))
274 0 : throw(AipsError("Error saving componentlist: "+err));
275 19 : vi.putModel(clrec, true, true);
276 19 : }
277 :
278 : }
279 : // See above commented-out code
280 : // cohDone+=vb->nRow();
281 : // pm.update(Double(cohDone));
282 : }
283 : }
284 37 : if(!incremental&&!initialized) initialized=true;
285 37 : }
286 :
287 :
288 37 : }
289 :
290 0 : void SkyEquation::predict(Bool incremental, MS::PredefinedColumns Type) {
291 :
292 :
293 :
294 : // if(noModelCol_p)
295 : // throw(AipsError("Cannot predict model vis without using scratch columns yet"));
296 0 : AlwaysAssert(cft_, AipsError);
297 0 : AlwaysAssert(sm_, AipsError);
298 : //AlwaysAssert(vs_, AipsError);
299 0 : if(sm_->numberOfModels()!= 0) AlwaysAssert(ok(),AipsError);
300 : // Initialize
301 0 : VisIter& vi= *wvi_p;
302 0 : checkVisIterNumRows(vi);
303 :
304 0 : VisBufferAutoPtr vb(vi);
305 :
306 : // Reset the visibilities only if this is not an incremental
307 : // change to the model
308 0 : Bool initialized=false;
309 : // **** predictcomponents doesn't do anything if incremental!
310 : // it gets components into vb->model, and writes to vi::(Type)
311 0 : predictComponents(incremental, initialized, Type);
312 :
313 : // Now do the images
314 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
315 :
316 0 : if( (sm_->isEmpty(model)) && (model==0) && !initialized && !incremental){
317 : // We are at the begining with an empty model start
318 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
319 0 : for (vi.origin(); vi.more(); vi++) {
320 0 : switch(Type) {
321 0 : case MS::MODEL_DATA:
322 0 : vb->setModelVisCube(Complex(0.0,0.0));
323 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
324 0 : break;
325 0 : case MS::DATA:
326 0 : vb->setModelVisCube(Complex(0.0,0.0));
327 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Observed);
328 0 : break;
329 0 : case MS::CORRECTED_DATA:
330 0 : vb->setModelVisCube(Complex(0.0,0.0));
331 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Corrected);
332 0 : break;
333 0 : default:
334 0 : throw (AipsError("Programmer made a wrong call"));
335 : }
336 : }
337 : }
338 : //if (!incremental&&!initialized) initialized=true;
339 :
340 :
341 : }
342 : // Don't bother with empty images
343 0 : if(!sm_->isEmpty(model)) {
344 :
345 : // Change the model polarization frame
346 0 : if(vb->polFrame()==MSIter::Linear) {
347 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
348 : StokesImageUtil::LINEAR);
349 : }
350 : else {
351 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
352 : StokesImageUtil::CIRCULAR);
353 : }
354 :
355 0 : scaleImage(model, incremental);
356 :
357 : // Reset the various SkyJones
358 0 : resetSkyJones();
359 :
360 :
361 : // Initialize get (i.e. Transform from Sky)
362 0 : vi.originChunks();
363 0 : vi.origin();
364 0 : initializeGet(* vb, 0, model, incremental);
365 0 : Int cohDone=0;
366 :
367 0 : ostringstream modelName;modelName<<"Model "<<model+1
368 0 : <<" : predicting coherences";
369 0 : ProgressMeter pm(1.0, Double(vi.numberCoh()),
370 0 : modelName, "", "", "", true);
371 : // Loop over all visibilities
372 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
373 0 : for (vi.origin(); vi.more(); vi++) {
374 0 : if(!incremental&&!initialized) {
375 0 : vb->setModelVisCube(Complex(0.0,0.0));
376 : // vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
377 0 : switch(Type) {
378 0 : case MS::MODEL_DATA:
379 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
380 0 : break;
381 0 : case MS::DATA:
382 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Observed);
383 0 : break;
384 0 : case MS::CORRECTED_DATA:
385 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Corrected);
386 0 : break;
387 0 : default:
388 0 : throw (AipsError("Programmer made a wrong call"));
389 : }
390 : }
391 : // get the model visibility (adds to vb->model)
392 : //get(vb,model,incremental);
393 : // this version takes Type and reads existing Type col instead
394 : // of hardcoded existing MODEL column
395 0 : get(* vb,model,incremental,Type);
396 : // and write it to VI
397 0 : switch(Type) {
398 0 : case MS::MODEL_DATA:
399 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Model);
400 0 : break;
401 0 : case MS::DATA:
402 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Observed);
403 0 : break;
404 0 : case MS::CORRECTED_DATA:
405 0 : vi.setVis(vb->modelVisCube(),VisibilityIterator::Corrected);
406 0 : break;
407 0 : default:
408 0 : throw (AipsError("Programmer made a wrong call"));
409 : }
410 0 : cohDone+=vb->nRow();
411 0 : pm.update(Double(cohDone));
412 : }
413 : }
414 0 : finalizeGet();
415 0 : unScaleImage(model, incremental);
416 0 : if(!incremental&&!initialized) initialized=true;
417 :
418 0 : } else {
419 0 : cout << "WARN: model is empty." << endl;
420 : }
421 : }
422 0 : }
423 :
424 : //----------------------------------------------------------------------
425 0 : void SkyEquation::gradientsChiSquared(const Matrix<Bool>&,
426 : SkyJones&) {
427 0 : throw(AipsError("SkyEquation:: solution for SkyJones not yet implemented"));
428 : }
429 :
430 : //----------------------------------------------------------------------
431 0 : void SkyEquation::gradientsChiSquared(Bool incremental, Bool commitModel) {
432 0 : AlwaysAssert(ok(),AipsError);
433 :
434 0 : if ((ft_->name() == "PBWProjectFT"))
435 : {
436 0 : ft_->setNoPadding(false);
437 0 : fullGradientsChiSquared(incremental);
438 : }
439 : else
440 : {
441 0 : Bool forceFull=true;
442 : // for these 2 gridders force incremental
443 0 : if((ft_->name() == "MosaicFT") || (ft_->name() == "WProjectFT"))
444 0 : forceFull=true;
445 :
446 0 : if( (sm_->numberOfModels() != 1) || !ft_->isFourier() || !incremental
447 0 : || forceFull){
448 0 : if(commitModel || !noModelCol_p){
449 0 : ft_->setNoPadding(false);
450 0 : fullGradientsChiSquared(incremental);
451 : }
452 : else{
453 : // For now use corrected_data...
454 0 : ft_->setNoPadding(true);
455 0 : cout << "This mode is gone...we should not be coming here" << endl;
456 : }
457 : }
458 : else {
459 0 : incrementGradientsChiSquared();
460 : }
461 : }
462 0 : }
463 : //----------------------------------------------------------------------
464 0 : void SkyEquation::fullGradientsChiSquared(Bool incremental) {
465 :
466 0 : AlwaysAssert(ok(),AipsError);
467 :
468 : // Predict the visibilities
469 0 : predict(incremental);
470 :
471 0 : sumwt = 0.0;
472 0 : chisq = 0.0;
473 :
474 : // Initialize the gradients
475 0 : sm_->initializeGradients();
476 :
477 0 : ROVisIter& vi(*rvi_p);
478 :
479 : // Loop over all models in SkyModel
480 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
481 :
482 : // Reset the various SkyJones
483 0 : resetSkyJones();
484 :
485 : // Loop over all visibilities and pixels
486 0 : checkVisIterNumRows(vi);
487 0 : VisBuffer vb(vi);
488 :
489 0 : if(sm_->isSolveable(model)) {
490 :
491 : // Initialize
492 0 : scaleImage(model, incremental);
493 :
494 0 : vi.originChunks();
495 0 : vi.origin();
496 :
497 : // Change the model polarization frame
498 0 : if(vb.polFrame()==MSIter::Linear) {
499 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
500 : StokesImageUtil::LINEAR);
501 : }
502 : else {
503 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
504 : StokesImageUtil::CIRCULAR);
505 : }
506 :
507 0 : initializePut(vb, model);
508 0 : Int cohDone=0;
509 :
510 0 : ostringstream modelName;modelName<<"Model "<<model+1
511 0 : <<" : transforming residuals";
512 0 : ProgressMeter pm(1.0, Double(vi.numberCoh()),
513 0 : modelName, "", "", "", true);
514 : // Loop over the visibilities, putting VisBuffers
515 :
516 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
517 0 : for (vi.origin(); vi.more(); vi++) {
518 0 : vb.modelVisCube()-=vb.correctedVisCube();
519 : // vb.setVisCube(vb.modelVisCube());
520 0 : put(vb, model, false, FTMachine::MODEL);
521 0 : cohDone+=vb.nRow();
522 0 : pm.update(Double(cohDone));
523 : }
524 : }
525 : // Do the transform, apply the SkyJones transformation
526 : // and sum the statistics for this model
527 0 : finalizePut(* vb_p, model);
528 0 : unScaleImage(model, incremental);
529 0 : }
530 0 : }
531 :
532 0 : fixImageScale();
533 :
534 : // Finish off any calculations needed internal to SkyModel
535 0 : sm_->finalizeGradients();
536 :
537 0 : }
538 :
539 : //----------------------------------------------------------------------
540 :
541 :
542 : //----------------------------------------------------------------------
543 :
544 :
545 0 : void SkyEquation::incrementGradientsChiSquared() {
546 :
547 0 : AlwaysAssert(ok(),AipsError);
548 :
549 : // Check to see if we need to make the XFRs
550 0 : if(!sm_->hasXFR(0)) makeComplexXFRs();
551 :
552 0 : ROVisIter& vi(*rvi_p);
553 :
554 : // Reset the various SkyJones
555 0 : resetSkyJones();
556 :
557 : // Loop over all models in SkyModel
558 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
559 :
560 0 : if(sm_->isSolveable(model)) {
561 :
562 0 : iDebug_p++;
563 :
564 0 : scaleDeltaImage(model);
565 0 : VisBuffer vb(vi);
566 0 : vi.origin();
567 :
568 : // Change the model polarization frame
569 0 : if(vb.polFrame()==MSIter::Linear) {
570 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
571 : StokesImageUtil::LINEAR);
572 : }
573 : else {
574 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
575 : StokesImageUtil::CIRCULAR);
576 : }
577 :
578 0 : Int numXFR=0;
579 0 : vi.originChunks();
580 0 : vi.origin();
581 0 : initializePutConvolve(vb, model, numXFR);
582 : // Iterate
583 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
584 0 : for (vi.origin(); vi.more(); vi++) {
585 0 : putConvolve(vb, model, numXFR);
586 : }
587 : }
588 : // Finish off
589 0 : finalizePutConvolve(* vb_p, model, numXFR);
590 0 : unScaleDeltaImage(model);
591 0 : }
592 : }
593 :
594 :
595 : // Finish off any calculations needed internal to SkyModel
596 0 : sm_->finalizeGradients();
597 0 : };
598 :
599 : //----------------------------------------------------------------------
600 0 : void SkyEquation::makeComplexXFRs()
601 : {
602 :
603 0 : AlwaysAssert(ok(),AipsError);
604 :
605 0 : ROVisIter& vi(*rvi_p);
606 :
607 : // Loop over all models in SkyModel
608 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
609 :
610 0 : if(sm_->isSolveable(model)) {
611 :
612 : // Loop over all visibilities and pixels
613 0 : VisBuffer vb(vi);
614 :
615 : // Change the model polarization frame
616 0 : if(vb.polFrame()==MSIter::Linear) {
617 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
618 : StokesImageUtil::LINEAR);
619 : }
620 : else {
621 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
622 : StokesImageUtil::CIRCULAR);
623 : }
624 :
625 : // Initialize put (i.e. transform to Sky) for this model
626 : // and XFR
627 0 : vi.origin();
628 0 : Int numXFR=0;
629 0 : vi.originChunks();
630 0 : vi.origin();
631 0 : initializePutXFR(vb, model, numXFR);
632 0 : Int cohDone=0;
633 :
634 0 : ostringstream modelName;modelName<<"Model "<<model+1
635 0 : <<" : transforming to PSF";
636 0 : ProgressMeter pm(1.0, Double(vi.numberCoh()),
637 0 : modelName, "", "", "", true);
638 : // Loop over the visibilities, putting VisBuffers
639 :
640 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
641 0 : for (vi.origin(); vi.more(); vi++) {
642 :
643 0 : vb.setVisCube(Complex(1.0,0.0));
644 0 : putXFR(vb, model, numXFR);
645 0 : cohDone+=vb.nRow();
646 0 : pm.update(Double(cohDone));
647 : }
648 : }
649 : // Do the transform, apply the SkyJones transformation
650 : // and sum the statistics for this model
651 0 : finalizePutXFR(* vb_p, model, numXFR);
652 0 : }
653 : }
654 0 : }
655 :
656 : //----------------------------------------------------------------------
657 : // Solve for a SkyJones
658 0 : Bool SkyEquation::solveSkyJones(SkyJones& sj) {
659 0 : setSkyJones(sj);
660 0 : return sj.solve(*this);
661 : }
662 :
663 : //----------------------------------------------------------------------
664 : // We make an approximate PSF for each plane.We only do this per model
665 : // since we may not need all PSFs.
666 : // ************* Note that this overwrites the model! ******************
667 0 : void SkyEquation::makeApproxPSF(Int model, ImageInterface<Float>& psf) {
668 :
669 0 : LogIO os(LogOrigin("SkyEquation", "makeApproxPSF"));
670 :
671 0 : AlwaysAssert(ok(), AipsError);
672 0 : AlwaysAssert(cft_, AipsError);
673 0 : AlwaysAssert(sm_, AipsError);
674 : //AlwaysAssert(vs_, AipsError);
675 :
676 0 : ft_->setNoPadding(noModelCol_p);
677 :
678 0 : isPSFWork_p= true; // avoid PB correction etc for PSF estimation
679 0 : Bool doPSF=true;
680 0 : if(ft_->name() == "MosaicFT") {
681 : // Reset the various SkyJones
682 0 : doPSF=false;
683 0 : resetSkyJones();
684 :
685 0 : VisIter& vi(*wvi_p);
686 0 : checkVisIterNumRows(vi);
687 : // Loop over all visibilities and pixels
688 0 : VisBuffer vb(vi);
689 :
690 0 : vi.originChunks();
691 0 : vi.origin();
692 :
693 : // Change the model polarization frame
694 0 : if(vb.polFrame()==MSIter::Linear) {
695 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
696 : StokesImageUtil::LINEAR);
697 : }
698 : else {
699 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
700 : StokesImageUtil::CIRCULAR);
701 : }
702 :
703 0 : IPosition start(4, sm_->image(model).shape()(0)/2,
704 0 : sm_->image(model).shape()(1)/2, 0, 0);
705 0 : IPosition shape(4, 1, 1, sm_->image(model).shape()(2), sm_->image(model).shape()(3));
706 0 : Array<Float> line(shape);
707 0 : TempImage<Float> savedModel(sm_->image(model).shape(),
708 0 : sm_->image(model).coordinates());
709 0 : savedModel.copyData(sm_->image(model));
710 0 : sm_->image(model).set(0.0);
711 0 : line=1.0;
712 0 : sm_->image(model).putSlice(line, start);
713 : //initializeGet(vb, -1, model, false);
714 0 : StokesImageUtil::From(sm_->cImage(model), static_cast <const ImageInterface<Float>& >(sm_->image(model)));
715 0 : ft_->initializeToVis(sm_->cImage(model),vb);
716 : // Loop over all visibilities
717 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
718 0 : for (vi.origin(); vi.more(); vi++) {
719 0 : vb.setModelVisCube(Complex(0.0,0.0));
720 : // get(vb, model, false);
721 0 : ft_->get(vb);
722 0 : vi.setVis(vb.modelVisCube(),VisibilityIterator::Model);
723 : }
724 : }
725 0 : finalizeGet();
726 :
727 0 : sm_->image(model).copyData(savedModel);
728 0 : }
729 :
730 : // Initialize the gradients
731 0 : sm_->initializeGradients();
732 :
733 :
734 0 : ROVisIter& vi(*rvi_p);
735 :
736 : // Reset the various SkyJones
737 0 : resetSkyJones();
738 :
739 0 : checkVisIterNumRows(vi);
740 : // Loop over all visibilities and pixels
741 0 : VisBuffer vb(vi);
742 :
743 :
744 0 : vi.originChunks();
745 0 : vi.origin();
746 :
747 : // Change the model polarization frame
748 0 : if(vb.polFrame()==MSIter::Linear) {
749 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
750 : StokesImageUtil::LINEAR);
751 : }
752 : else {
753 0 : StokesImageUtil::changeCStokesRep(sm_->cImage(model),
754 : StokesImageUtil::CIRCULAR);
755 : }
756 :
757 :
758 0 : initializePut(vb, model);
759 :
760 : // Loop over the visibilities, putting VisBuffers
761 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
762 0 : for (vi.origin(); vi.more(); vi++) {
763 : // vb.setVisCube(vb.modelVisCube());
764 0 : put(vb, model, doPSF, FTMachine::MODEL);
765 : }
766 : }
767 :
768 : // Do the transform, apply the SkyJones transformation
769 0 : finalizePut(* vb_p, model);
770 0 : sm_->finalizeGradients();
771 0 : fixImageScale();
772 0 : LatticeExpr<Float> le(iif(sm_->ggS(model)>(0.0),
773 0 : (sm_->gS(model)/sm_->ggS(model)), 0.0));
774 0 : psf.copyData(le);
775 0 : LatticeExprNode maxPSF=max(psf);
776 0 : Float maxpsf=maxPSF.getFloat();
777 0 : if(abs(maxpsf-1.0) > 1e-3) {
778 : os << "Maximum of approximate PSF for field " << model << " = "
779 0 : << maxpsf << " : renormalizing to unity" << LogIO::POST;
780 : }
781 0 : if(maxpsf > 0.0 ){
782 0 : LatticeExpr<Float> len(psf/maxpsf);
783 0 : psf.copyData(len);
784 :
785 0 : }
786 : else{
787 0 : throw(AipsError("SkyEquation:: PSF calculation resulted in a PSF lesser than 0 !"));
788 :
789 : }
790 :
791 :
792 :
793 0 : isPSFWork_p=false; // resseting this flag so that subsequent calculation uses
794 : // the right SkyJones correction;
795 :
796 :
797 0 : }
798 :
799 0 : void SkyEquation::makeApproxPSF(PtrBlock<ImageInterface<Float> *>& PSFs) {
800 :
801 0 : Int nmodels=PSFs.nelements();
802 0 : for (Int model=0; model < nmodels; ++model){
803 0 : makeApproxPSF(model, *(PSFs[model]));
804 :
805 : }
806 :
807 :
808 0 : }
809 : //----------------------------------------------------------------------
810 : // Solve for a SkyModel
811 0 : Bool SkyEquation::solveSkyModel() {
812 0 : return sm_->solve(*this);
813 : }
814 :
815 : // ***************************************************************************
816 : // ******************** Start of private member functions *******************
817 : // ***************************************************************************
818 :
819 : // Initialize
820 0 : void SkyEquation::initializeGet(const VisBuffer& vb, Int row, Int model,
821 : Bool incremental) {
822 :
823 0 : AlwaysAssert(ok(),AipsError);
824 0 : if(incremental) {
825 0 : applySkyJones(vb, row, sm_->deltaImage(model), sm_->cImage(model));
826 : }
827 : else {
828 0 : applySkyJones(vb, row, sm_->image(model), sm_->cImage(model));
829 : }
830 0 : ft_->initializeToVis(sm_->cImage(model),vb);
831 0 : }
832 : //
833 :
834 : // Add the sky visibility for this coherence sample
835 0 : VisBuffer& SkyEquation::get(VisBuffer& result, Int model,
836 : Bool incremental,
837 : MS::PredefinedColumns Type) {
838 :
839 0 : AlwaysAssert(ok(),AipsError);
840 0 : Int nRow=result.nRow();
841 :
842 : // we'll always return MODEL, but maybe expect existing data in another col
843 : // yes this is a bit convoluted - probably should change whatever calls
844 : // get to expect different columns.
845 0 : switch(Type) {
846 0 : case MS::MODEL_DATA:
847 0 : result.modelVisCube();
848 0 : break;
849 0 : case MS::DATA:
850 0 : result.visCube();
851 0 : result.setModelVisCube(result.visCube());
852 0 : break;
853 0 : case MS::CORRECTED_DATA:
854 0 : result.correctedVisCube();
855 0 : result.setModelVisCube(result.correctedVisCube());
856 0 : break;
857 0 : default:
858 0 : throw (AipsError("Programmer made a wrong call"));
859 : }
860 :
861 : //result.modelVisCube(); // get the visibility so vb will have it
862 : //VisBufferAutoPtr vb (result);
863 : // We need msId and mscolumns from result in the vb passed to FT below.
864 : // I can't find a public method to copy that ms info over, so instead,
865 : // save the datacube and add it back at the end of this method.
866 : VisBuffer *oldvb,*vb;
867 0 : vb = &result; // this will get replaced by FT->get
868 0 : oldvb = new VisBuffer(result);
869 0 : oldvb->setModelVisCube(result.modelVisCube());
870 :
871 0 : Bool FTChanged=changedFTMachine(* vb);
872 :
873 : // we might need to recompute the "sky" for every single row, but we
874 : // avoid this if possible.
875 0 : Bool internalChanges=false; // Does this VB change inside itself?
876 0 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
877 0 : if(ft_->name()!="MosaicFT"){
878 0 : changedSkyJonesLogic(result, firstOneChanges, internalChanges);
879 : }
880 0 : if(internalChanges) {
881 : // Yes there are changes within this buffer: go row by row.
882 : // This will automatically catch a change in the FTMachine so
883 : // we don't have to check for that.
884 0 : for (Int row=0; row<nRow; row++) {
885 0 : finalizeGet();
886 0 : initializeGet(result, row, model, incremental);
887 0 : ft_->get(* vb,row);
888 : }
889 : }
890 0 : else if (FTChanged||firstOneChanges) {
891 : // This buffer has changed wrt the previous buffer, but
892 : // this buffer has no changes within it. Again we don't need to
893 : // check for the FTMachine changing.
894 0 : finalizeGet();
895 0 : initializeGet(result, 0, model, incremental);
896 0 : ft_->get(* vb);
897 : }
898 : else {
899 0 : ft_->get(* vb);
900 : }
901 0 : result.modelVisCube()+=oldvb->modelVisCube();
902 0 : return result;
903 : }
904 :
905 0 : void SkyEquation::finalizeGet()
906 : {
907 : // Do nothing
908 : // MFSkyEquation has something to do, we just
909 : // need to have "get()" calling finalizeGet()
910 0 : };
911 :
912 :
913 : // Add the sky visibility for this component
914 0 : VisBuffer& SkyEquation::get(VisBuffer& result,
915 : const SkyComponent& component)
916 : {
917 :
918 0 : AlwaysAssert(sm_,AipsError);
919 0 : AlwaysAssert(cft_,AipsError);
920 :
921 0 : Int nRow=result.nRow();
922 :
923 0 : result.modelVisCube(); // get the visibility so vb will have it
924 0 : VisBufferAutoPtr vb (result);
925 0 : SkyComponent corruptedComponent = component.copy();
926 0 : if(vb->polFrame()==MSIter::Linear) {
927 0 : if(corruptedComponent.flux().pol()==ComponentType::STOKES) {
928 0 : corruptedComponent.flux().convertPol(ComponentType::LINEAR);
929 : }
930 : }
931 : else {
932 0 : if(corruptedComponent.flux().pol()==ComponentType::STOKES) {
933 0 : corruptedComponent.flux().convertPol(ComponentType::CIRCULAR);
934 : }
935 : }
936 : // we might need to recompute the "sky" for every single row, but
937 : // we avoid this, if possible
938 0 : Bool internalChanges=false; // Does this VB change inside itself
939 0 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
940 0 : changedSkyJonesLogic(result,firstOneChanges,internalChanges);
941 0 : if (internalChanges) // yes, we have to go row by row
942 0 : for (Int row=0;row<nRow;++row) {
943 0 : SkyComponent tempComponent=corruptedComponent.copy();
944 0 : applySkyJones(tempComponent,result,row);
945 0 : cft_->get(* vb,tempComponent,row);
946 0 : }
947 : else { // we don't have a cache of corruptedComponent, therefore
948 : // firstOneChanges is equivalent to the default case
949 0 : applySkyJones(corruptedComponent, result, 0);
950 0 : cft_->get(* vb, corruptedComponent);
951 : }
952 0 : result.modelVisCube()+=vb->modelVisCube();
953 0 : return result;
954 0 : }
955 :
956 :
957 :
958 : // Add the sky visibility for this component
959 2888 : VisBuffer& SkyEquation::get(VisBuffer& result,
960 : const ComponentList& compList)
961 : // MS::PredefinedColumns Type) {
962 : {
963 :
964 2888 : AlwaysAssert(sm_,AipsError);
965 2888 : AlwaysAssert(cft_,AipsError);
966 :
967 2888 : Int nRow=result.nRow();
968 :
969 2888 : result.modelVisCube(); // get the visibility so vb will have it
970 : //VisBufferAutoPtr vb(result);
971 2888 : VisBuffer vb(result); // method only called using writable VI so no ROVIA
972 : //Above state copy is incomplete like msId is not copied...make sure to use
973 : //result when metadata is important.
974 2888 : result.setModelVisCube(Complex(0.0,0.0));
975 :
976 : // CURRENTLY we do not have the applySkyJones code which
977 : // works on ComponentLists; so we need to get out the
978 : // individual SkyComponents here
979 :
980 : // we might need to recompute the "sky" for every single row, but
981 : // we avoid this, if possible
982 2888 : Bool internalChanges=false; // Does this VB change inside itself
983 2888 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
984 2888 : changedSkyJonesLogic(result,firstOneChanges,internalChanges);
985 :
986 2888 : uInt ncomponents=compList.nelements();
987 5776 : for (uInt icomp=0;icomp<ncomponents;icomp++) {
988 2888 : SkyComponent component=compList.component(icomp).copy();
989 2888 : if(vb.polFrame()==MSIter::Linear) {
990 0 : if(component.flux().pol()==ComponentType::STOKES) {
991 0 : component.flux().convertPol(ComponentType::LINEAR);
992 : }
993 : }
994 : else {
995 2888 : if(component.flux().pol()==ComponentType::STOKES) {
996 2888 : component.flux().convertPol(ComponentType::CIRCULAR);
997 : }
998 : }
999 2888 : if (internalChanges) // yes, we have to go row by row
1000 0 : for (Int row=0;row<nRow;++row) {
1001 0 : SkyComponent tempComponent=component.copy();
1002 0 : applySkyJones(tempComponent,result,row);
1003 0 : cft_->get(vb,tempComponent,row);
1004 0 : }
1005 : else { // we don't have a cache for component, therefore
1006 : // firstOneChanges is equivalent to the default case
1007 2888 : applySkyJones(component, result, 0);
1008 2888 : cft_->get(vb, component);
1009 : }
1010 2888 : result.modelVisCube()+=vb.modelVisCube();
1011 2888 : }
1012 : // Now add into the existing model visibility cube
1013 2888 : return result;
1014 2888 : }
1015 :
1016 :
1017 : // Corrupt a SkyComponent
1018 2888 : SkyComponent& SkyEquation::applySkyJones(SkyComponent& corruptedComponent,
1019 : const VisBuffer& vb,
1020 : Int row)
1021 : {
1022 2888 : if(!isPSFWork_p){
1023 : //The last of the bolean in the following line forces a full spectral corruption
1024 : // May be should do it on detection of fractional bandwidth
1025 2888 : if(ej_) ej_->apply(corruptedComponent,corruptedComponent,vb,row, True, True);
1026 2888 : if(dj_) dj_->apply(corruptedComponent,corruptedComponent,vb,row);
1027 2888 : if(tj_) tj_->apply(corruptedComponent,corruptedComponent,vb,row);
1028 2888 : if(fj_) fj_->apply(corruptedComponent,corruptedComponent,vb,row);
1029 : }
1030 2888 : return corruptedComponent;
1031 : }
1032 :
1033 0 : void SkyEquation::initializePut(const VisBuffer& vb, Int model) {
1034 0 : AlwaysAssert(ok(),AipsError);
1035 0 : ift_->initializeToSky(sm_->cImage(model),sm_->weight(model),vb);
1036 0 : assertSkyJones(vb, -1);
1037 0 : vb_p->assign(vb, false);
1038 0 : vb_p->updateCoordInfo();
1039 0 : }
1040 :
1041 :
1042 0 : void SkyEquation::put(const VisBuffer & vb, Int model, Bool dopsf, FTMachine::Type col) {
1043 :
1044 0 : AlwaysAssert(ok(),AipsError);
1045 :
1046 0 : Bool IFTChanged=changedIFTMachine(vb);
1047 :
1048 : // we might need to recompute the "sky" for every single row, but we
1049 : // avoid this if possible.
1050 :
1051 :
1052 0 : Int nRow=vb.nRow();
1053 0 : Bool internalChanges=false; // Does this VB change inside itself?
1054 0 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
1055 0 : if(ft_->name() != "MosaicFT"){
1056 0 : changedSkyJonesLogic(vb, firstOneChanges, internalChanges);
1057 : }
1058 0 : if(internalChanges) {
1059 : // Yes there are changes: go row by row.
1060 0 : for (Int row=0; row<nRow; row++) {
1061 0 : if(IFTChanged||changedSkyJones(vb,row)) {
1062 : // Need to apply the SkyJones from the previous row
1063 : // and finish off before starting with this row
1064 0 : finalizePut(* vb_p, model);
1065 0 : initializePut(vb, model);
1066 : }
1067 0 : ift_->put(vb, row, dopsf, col);
1068 : }
1069 : }
1070 0 : else if (IFTChanged||firstOneChanges) {
1071 :
1072 0 : if(!isBeginingOfSkyJonesCache_p){
1073 0 : finalizePut(* vb_p, model);
1074 : }
1075 0 : initializePut(vb, model);
1076 0 : isBeginingOfSkyJonesCache_p=false;
1077 0 : ift_->put(vb, -1, dopsf, col);
1078 : }
1079 : else {
1080 0 : ift_->put(vb, -1, dopsf, col);
1081 : }
1082 :
1083 0 : isBeginingOfSkyJonesCache_p=false;
1084 0 : }
1085 :
1086 :
1087 :
1088 :
1089 0 : void SkyEquation::finalizePut(const VisBuffer& vb, Int model) {
1090 :
1091 : // Actually do the transform. Update weights as we do so.
1092 0 : ift_->finalizeToSky();
1093 : // 1. Now get the (unnormalized) image and add the
1094 : // weight to the summed weights
1095 0 : Matrix<Float> delta;
1096 0 : sm_->cImage(model).copyData(ift_->getImage(delta, false));
1097 0 : sm_->weight(model)+=delta;
1098 : // 2. Apply the SkyJones and add to grad chisquared
1099 0 : applySkyJonesInv(vb, -1, sm_->cImage(model), sm_->work(model),
1100 0 : sm_->gS(model));
1101 :
1102 : // 3. Apply the square of the SkyJones and add this to gradgrad chisquared
1103 0 : applySkyJonesSquare(vb, -1, sm_->weight(model), sm_->work(model),
1104 0 : sm_->ggS(model));
1105 :
1106 :
1107 : // 4. Finally, we add the statistics
1108 0 : sm_->addStatistics(sumwt, chisq);
1109 0 : }
1110 :
1111 :
1112 :
1113 :
1114 0 : void SkyEquation::initializePutXFR(const VisBuffer& vb, Int model,
1115 : Int numXFR) {
1116 0 : AlwaysAssert(ok(),AipsError);
1117 0 : Matrix<Float> weight;
1118 0 : ift_->initializeToSky(sm_->XFR(model, numXFR), weight, vb);
1119 0 : assertSkyJones(vb, -1);
1120 0 : vb_p->assign(vb, false);
1121 0 : vb_p->updateCoordInfo();
1122 0 : }
1123 :
1124 0 : void SkyEquation::putXFR(const VisBuffer & vb, Int model, Int& numXFR) {
1125 :
1126 0 : AlwaysAssert(ok(),AipsError);
1127 :
1128 0 : Bool IFTChanged=changedIFTMachine(vb);
1129 :
1130 0 : Bool internalChanges=false; // Does this VB change inside itself?
1131 0 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
1132 0 : changedSkyJonesLogic(vb, firstOneChanges, internalChanges);
1133 0 : if(internalChanges) {
1134 : // Yes there are changes within this buffer: go row by row.
1135 0 : Int nRow=vb.nRow();
1136 0 : for (Int row=0; row<nRow; row++) {
1137 0 : if(IFTChanged||changedSkyJones(vb,row)) {
1138 : // Need to apply the SkyJones from the previous row
1139 : // and finish off before starting with this row
1140 0 : finalizePutXFR(* vb_p, model, numXFR); // also, this needs to know about
1141 : // the vb row number
1142 0 : numXFR++;
1143 0 : initializePutXFR(vb, model, numXFR);
1144 : }
1145 0 : ift_->put(vb, row, true);
1146 : }
1147 : }
1148 0 : else if (IFTChanged||firstOneChanges) {
1149 0 : finalizePutXFR(* vb_p, model, numXFR);
1150 0 : numXFR++;
1151 0 : initializePutXFR(vb, model, numXFR);
1152 0 : ift_->put(vb, -1, true);
1153 : } else {
1154 0 : ift_->put(vb, -1, true);
1155 : }
1156 0 : }
1157 :
1158 0 : void SkyEquation::finalizePutXFR(const VisBuffer& vb, Int model, Int numXFR)
1159 : {
1160 : // Actually do the transform. FFT back to the visibility plane
1161 0 : ift_->finalizeToSky();
1162 0 : Matrix<Float> weights;
1163 0 : sm_->XFR(model, numXFR).copyData(ift_->getImage(weights, false));
1164 0 : LatticeFFT::cfft2d(sm_->XFR(model, numXFR));
1165 0 : assertSkyJones(vb, -1);
1166 : // TempImage<Complex> *tip = ( TempImage<Complex> *) ( &(sm_->XFR(model, numXFR)) );
1167 : // tip->tempClose();
1168 0 : }
1169 :
1170 : // Here we do the whole thing: apply SkyJones, FFT, cross-multiply
1171 : // by the XFR, inverseFFT and then apply SkyJones again
1172 0 : void SkyEquation::initializePutConvolve(const VisBuffer& vb, Int model,
1173 : Int numXFR)
1174 : {
1175 0 : LogIO os(LogOrigin("SkyEquation", "initializePutConvolve"));
1176 0 : AlwaysAssert(ok(),AipsError);
1177 0 : AlwaysAssert(model>-1, AipsError);
1178 0 : AlwaysAssert(numXFR>-1, AipsError);
1179 0 : assertSkyJones(vb, -1);
1180 0 : vb_p->assign(vb, false);
1181 0 : vb_p->updateCoordInfo();
1182 0 : }
1183 :
1184 : // Here we step through the visbuffer and do the convolution
1185 : // if something changes
1186 0 : void SkyEquation::putConvolve(const VisBuffer & vb, Int model, Int& numXFR) {
1187 :
1188 0 : AlwaysAssert(ok(),AipsError);
1189 :
1190 0 : Int nRow=vb.nRow();
1191 0 : Bool internalChanges=false; // Does this VB change inside itself?
1192 0 : Bool firstOneChanges=false; // Has this VB changed from the previous one?
1193 0 : changedSkyJonesLogic(vb, firstOneChanges, internalChanges);
1194 0 : if(internalChanges) {
1195 : // Yes there are changes within this buffer: go row by row.
1196 0 : for (Int row=0; row<nRow; row++) {
1197 0 : if(changedSkyJones(vb,row)) {
1198 : // Need to apply the SkyJones from the previous row
1199 : // and finish off before starting with this row
1200 0 : finalizePutConvolve(* vb_p, model, numXFR);
1201 0 : numXFR++;
1202 0 : initializePutConvolve(vb, model, numXFR);
1203 : }
1204 : }
1205 : }
1206 0 : else if (firstOneChanges) {
1207 0 : finalizePutConvolve(* vb_p, model, numXFR);
1208 0 : numXFR++;
1209 0 : initializePutConvolve(vb, model, numXFR);
1210 : }
1211 0 : };
1212 :
1213 : // Here we do the convolution and transform back
1214 0 : void SkyEquation::finalizePutConvolve(const VisBuffer& vb, Int model,
1215 : Int numXFR)
1216 : {
1217 0 : LogIO os(LogOrigin("SkyEquation", "finalizePutConvolve"));
1218 :
1219 0 : applySkyJones(vb, -1, sm_->deltaImage(model), sm_->cImage(model));
1220 0 : LatticeFFT::cfft2d(sm_->cImage(model));
1221 0 : LatticeExpr<Complex> latticeExpr(conj(sm_->XFR(model, numXFR))*sm_->cImage(model));
1222 0 : sm_->cImage(model).copyData(latticeExpr);
1223 0 : LatticeFFT::cfft2d(sm_->cImage(model), false);
1224 0 : applySkyJonesInv(vb, -1, sm_->cImage(model), sm_->work(model),
1225 0 : sm_->gS(model));
1226 0 : }
1227 :
1228 2888 : void SkyEquation::changedSkyJonesLogic(const VisBuffer& vb,
1229 : Bool& firstOneChanges,
1230 : Bool& internalChanges)
1231 : {
1232 2888 : if(changedSkyJones(vb,0)) {
1233 0 : firstOneChanges=true;
1234 : }
1235 2888 : Int lastrow = -1;
1236 2888 : if(changedSkyJonesBuffer(vb, 0, lastrow)) {
1237 0 : internalChanges=true;
1238 : }
1239 5776 : return;
1240 : };
1241 :
1242 :
1243 : // Has the total SkyJones changed since the last application of the SkyJones?
1244 2888 : Bool SkyEquation::changedSkyJones(const VisBuffer& vb, Int row) {
1245 2888 : if(ej_) if(ej_->changed(vb,row)) return true; // Electric field pattern
1246 2888 : if(dj_) if(dj_->changed(vb,row)) return true; // Polarization field pattern
1247 2888 : if(tj_) if(tj_->changed(vb,row)) return true; // Atmospheric gain
1248 2888 : if(fj_) if(fj_->changed(vb,row)) return true; // Faraday rotation
1249 2888 : return false;
1250 : };
1251 :
1252 :
1253 : // Has the FT Machine changed since the last application ?
1254 0 : Bool SkyEquation::changedFTMachine(const VisBuffer& vb) {
1255 0 : return ft_->changed(vb);
1256 : };
1257 :
1258 : // Has the Inverse FT Machine changed since the last application ?
1259 0 : Bool SkyEquation::changedIFTMachine(const VisBuffer& vb) {
1260 0 : return ift_->changed(vb);
1261 : };
1262 :
1263 :
1264 : // Does the SkyJones change in this buffer starting at row1?
1265 2888 : Bool SkyEquation::changedSkyJonesBuffer
1266 : (const VisBuffer& vb, Int row1, Int& row2) {
1267 2888 : Int row2temp = 0;
1268 2888 : row2 = vb.nRow() - 1;
1269 2888 : Bool didChange = false;
1270 2888 : if(ej_) { // Electric field pattern
1271 0 : if(ej_->changedBuffer(vb,row1, row2temp)) {
1272 0 : didChange = true;
1273 0 : row2 = min (row2, row2temp);
1274 : }
1275 : }
1276 2888 : if(dj_) { // Polarization field pattern
1277 0 : if(dj_->changedBuffer(vb,row1, row2temp)) {
1278 0 : didChange = true;
1279 0 : row2 = min (row2, row2temp);
1280 : }
1281 : }
1282 2888 : if(tj_) { // Atmospheric gain
1283 0 : if(tj_->changedBuffer(vb,row1, row2temp)) {
1284 0 : didChange = true;
1285 0 : row2 = min (row2, row2temp);
1286 : }
1287 : }
1288 2888 : if(fj_) { // Faraday rotation
1289 0 : if(fj_->changedBuffer(vb,row1, row2temp)) {
1290 0 : didChange = true;
1291 0 : row2 = min (row2, row2temp);
1292 : }
1293 : }
1294 2888 : return didChange;
1295 : };
1296 :
1297 :
1298 37 : void SkyEquation::setPhaseCenterTime(const Double time){
1299 :
1300 37 : if(cft_)
1301 37 : cft_->setPhaseCenterTime(time);
1302 37 : if(ft_)
1303 37 : ft_->setPhaseCenterTime(time);
1304 37 : if(ift_)
1305 37 : ift_->setPhaseCenterTime(time);
1306 37 : }
1307 0 : Double SkyEquation::getPhaseCenterTime(){
1308 0 : if(cft_)
1309 0 : return cft_->getPhaseCenterTime();
1310 0 : if(ft_)
1311 0 : return ft_->getPhaseCenterTime();
1312 0 : return -1.0;
1313 : }
1314 :
1315 : // Reset all of the SkyJones to initial state
1316 37 : void SkyEquation::resetSkyJones()
1317 : {
1318 37 : if(ej_) ej_->reset(); // Electric field pattern
1319 37 : if(dj_) dj_->reset(); // Polarization field pattern
1320 37 : if(tj_) tj_->reset(); // Atmospheric gain
1321 37 : if(fj_) fj_->reset(); // Faraday rotation
1322 :
1323 37 : isBeginingOfSkyJonesCache_p=true;
1324 :
1325 37 : };
1326 :
1327 : // Assure that we've taken care of the SkyJones terms
1328 : // (needed for e.g. putXFR)
1329 0 : void SkyEquation::assertSkyJones(const VisBuffer& vb, Int row) {
1330 0 : if(ej_) ej_->assure(vb,row); // Electric field pattern
1331 0 : if(dj_) dj_->assure(vb,row); // Polarization field pattern
1332 0 : if(tj_) tj_->assure(vb,row); // Atmospheric gain
1333 0 : if(fj_) fj_->assure(vb,row); // Faraday rotation
1334 0 : };
1335 :
1336 0 : ImageInterface<Complex>& SkyEquation::applySkyJones(const VisBuffer& vb,
1337 : Int row,
1338 : ImageInterface<Float>& in,
1339 : ImageInterface<Complex>& out) {
1340 :
1341 : //Pol axis need not be same
1342 0 : AlwaysAssert(in.shape()[0]==out.shape()[0], AipsError);
1343 0 : AlwaysAssert(in.shape()[1]==out.shape()[1], AipsError);
1344 0 : AlwaysAssert(in.shape()[3]==out.shape()[3], AipsError);
1345 : // Convert from Stokes to Complex
1346 0 : StokesImageUtil::From(out, in);
1347 :
1348 : // Now apply the SkyJones as needed
1349 0 : if(!isPSFWork_p && (!(ift_->name().contains("MosaicFT")))){
1350 0 : if(ej_) ej_->apply(out,out,vb,row,true);
1351 0 : if(dj_) dj_->apply(out,out,vb,row,true);
1352 0 : if(tj_) tj_->apply(out,out,vb,row,true);
1353 0 : if(fj_) fj_->apply(out,out,vb,row,true);
1354 : }
1355 0 : return out;
1356 : };
1357 :
1358 :
1359 : // Calculate gradChisq. In the SkyModel, this is used to update
1360 : // the estimated image.
1361 0 : void SkyEquation::applySkyJonesInv(const VisBuffer& vb, Int row,
1362 : ImageInterface<Complex>& in,
1363 : ImageInterface<Float>& work,
1364 : ImageInterface<Float>& gS) {
1365 :
1366 0 : AlwaysAssert(in.shape()[0]==work.shape()[0], AipsError);
1367 0 : AlwaysAssert(in.shape()[1]==work.shape()[1], AipsError);
1368 0 : AlwaysAssert(in.shape()[3]==work.shape()[3], AipsError);
1369 0 : AlwaysAssert(gS.shape()==work.shape(), AipsError);
1370 :
1371 : // Apply the various SkyJones to the current image
1372 : // remembering to apply the Jones in the backward
1373 : // direction
1374 0 : if(!isPSFWork_p && (ift_->name() != "MosaicFT")){
1375 :
1376 0 : if(ej_) ej_->apply(in,in,vb,row,false);
1377 0 : if(dj_) dj_->apply(in,in,vb,row,false);
1378 0 : if(tj_) tj_->apply(in,in,vb,row,false);
1379 0 : if(fj_) fj_->apply(in,in,vb,row,false);
1380 : }
1381 : // Convert to IQUV format
1382 0 : if(isPSFWork_p)
1383 : {
1384 : // For the PSF, choose only those stokes planes that have a valid PSF
1385 0 : StokesImageUtil::ToStokesPSF(work,in);
1386 : }
1387 : else
1388 : {
1389 0 : StokesImageUtil::To(work,in);
1390 : }
1391 :
1392 : // Now add to the existing gradChisq image
1393 0 : LatticeExpr<Float> le(gS+work);
1394 0 : gS.copyData(le);
1395 :
1396 0 : };
1397 :
1398 :
1399 :
1400 : // Calculate gradgradChisq. In the SkyModel, this is
1401 : // used to update the estimated image.
1402 0 : void SkyEquation::applySkyJonesSquare(const VisBuffer& vb, Int row,
1403 : Matrix<Float>& weights,
1404 : ImageInterface<Float>& work,
1405 : ImageInterface<Float>& ggS) {
1406 :
1407 0 : AlwaysAssert(work.shape()==ggS.shape(), AipsError);
1408 :
1409 : // First fill the work image with the appropriate value
1410 : // of the weight.
1411 0 : ift_->getWeightImage(work, weights);
1412 : // Apply SkyJones as needed
1413 0 : if((ft_->name() != "MosaicFT") ){
1414 0 : if(ej_) ej_->applySquare(work,work,vb,row);
1415 0 : if(dj_) dj_->applySquare(work,work,vb,row);
1416 0 : if(tj_) tj_->applySquare(work,work,vb,row);
1417 0 : if(fj_) fj_->applySquare(work,work,vb,row);
1418 : }
1419 : // Now add to the existing gradgradChisq image
1420 0 : if((ft_->name() != "MosaicFT")){
1421 0 : LatticeExpr<Float> le(ggS+work);
1422 0 : ggS.copyData(le);
1423 0 : }
1424 : else{
1425 0 : ggS.copyData(work);
1426 : }
1427 :
1428 0 : };
1429 :
1430 :
1431 0 : Bool SkyEquation::ok() {
1432 :
1433 0 : AlwaysAssert(sm_,AipsError);
1434 : //AlwaysAssert(vs_,AipsError);
1435 0 : AlwaysAssert(ft_,AipsError);
1436 0 : AlwaysAssert(ift_,AipsError);
1437 0 : AlwaysAssert(rvi_p, AipsError);
1438 :
1439 0 : return(true);
1440 : }
1441 :
1442 :
1443 0 : void SkyEquation::scaleImage(Int model)
1444 : {
1445 :
1446 0 : if (sm_->doFluxScale(model)){
1447 :
1448 0 : LatticeExpr<Float> latticeExpr( iif(sm_->fluxScale(model) <= (0.0), 0.0, (sm_->image(model))/(sm_->fluxScale(model))) );
1449 0 : sm_->image(model).copyData(latticeExpr);
1450 0 : sm_->fluxScale(model).clearCache();
1451 0 : sm_->image(model).clearCache();
1452 0 : }
1453 0 : };
1454 :
1455 0 : void SkyEquation::unScaleImage(Int model)
1456 : {
1457 :
1458 0 : if ( sm_->doFluxScale(model)){
1459 :
1460 0 : LatticeExpr<Float> latticeExpr( sm_->image(model) * (sm_->fluxScale(model)) );
1461 0 : sm_->image(model).copyData(latticeExpr);
1462 0 : sm_->fluxScale(model).clearCache();
1463 0 : sm_->image(model).clearCache();
1464 0 : }
1465 0 : };
1466 :
1467 0 : void SkyEquation::scaleImage(Int model, Bool incremental)
1468 : {
1469 0 : if (incremental) {
1470 0 : scaleDeltaImage(model);
1471 : } else {
1472 0 : scaleImage(model);
1473 : }
1474 0 : };
1475 :
1476 0 : void SkyEquation::unScaleImage(Int model, Bool incremental)
1477 : {
1478 0 : if (incremental) {
1479 0 : unScaleDeltaImage(model);
1480 : } else {
1481 0 : unScaleImage(model);
1482 : }
1483 0 : };
1484 :
1485 0 : void SkyEquation::scaleDeltaImage(Int model)
1486 : {
1487 0 : if ((sm_->doFluxScale(model))){
1488 0 : sm_->deltaImage(model).copyData( (LatticeExpr<Float>)
1489 0 : (iif(sm_->fluxScale(model) <= (0.0), 0.0,
1490 0 : ((sm_->deltaImage(model))/(sm_->fluxScale(model))) )) );
1491 0 : sm_->fluxScale(model).clearCache();
1492 0 : sm_->deltaImage(model).clearCache();
1493 : }
1494 :
1495 0 : };
1496 :
1497 0 : void SkyEquation::getCoverageImage(Int model, ImageInterface<Float>& im){
1498 0 : if ((sm_->doFluxScale(model))){
1499 0 : ift_->getFluxImage(im);
1500 : }
1501 :
1502 0 : }
1503 :
1504 :
1505 0 : void SkyEquation::unScaleDeltaImage(Int model)
1506 : {
1507 0 : if ( (sm_->doFluxScale(model))){
1508 0 : LatticeExpr<Float> latticeExpr( sm_->deltaImage(model) * (sm_->fluxScale(model)) );
1509 0 : sm_->deltaImage(model).copyData(latticeExpr);
1510 0 : sm_->fluxScale(model).clearCache();
1511 0 : sm_->deltaImage(model).clearCache();
1512 0 : }
1513 0 : };
1514 :
1515 0 : void SkyEquation::fixImageScale()
1516 : {
1517 0 : LogIO os(LogOrigin("SkyEquation", "fixImageScale"));
1518 :
1519 : // make a minimum value to ggS
1520 : // This has the same effect as Sault Weighting, but
1521 : // is implemented somewhat differently.
1522 : // We also keep the fluxScale(mod) images around to
1523 : // undo the weighting.
1524 :
1525 0 : if(ej_ || (ft_->name() == "MosaicFT") ) {
1526 0 : Float ggSMax=0.0;
1527 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
1528 :
1529 0 : LatticeExprNode LEN = max( sm_->ggS(model) );
1530 0 : ggSMax = max(ggSMax,LEN.getFloat());
1531 0 : }
1532 :
1533 0 : ggSMax_p=ggSMax;
1534 : Float ggSMin1;
1535 : Float ggSMin2;
1536 :
1537 0 : ggSMin1 = ggSMax * constPB_p * constPB_p;
1538 0 : ggSMin2 = ggSMax * minPB_p * minPB_p;
1539 :
1540 :
1541 : /*Don't print this for now
1542 : if (scaleType_p == "SAULT") {
1543 : os << "Using SAULT image plane weighting" << LogIO::POST;
1544 : }
1545 : else {
1546 : os << "Using No image plane weighting" << LogIO::POST;
1547 : }
1548 : */
1549 :
1550 0 : for (Int model=0;model<sm_->numberOfModels();model++) {
1551 0 : sm_->fluxScale(model).removeRegion ("mask0", RegionHandler::Any, false);
1552 0 : if ((ft_->name()!="MosaicFT")) {
1553 0 : if(scaleType_p=="SAULT"){
1554 :
1555 : // Adjust flux scale to account for ggS being truncated at ggSMin1
1556 : // Below ggSMin2, set flux scale to 0.0
1557 : // FluxScale * image => true brightness distribution, but
1558 : // noise increases at edge.
1559 : // if ggS < ggSMin2, set to Zero;
1560 : // if ggS > ggSMin2 && < ggSMin1, set to ggSMin1/ggS
1561 : // if ggS > ggSMin1, set to 1.0
1562 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1563 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1564 0 : sqrt((sm_->ggS(model))/ggSMin1) )) );
1565 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1566 0 : (iif(sm_->ggS(model) > (ggSMin1), 1.0,
1567 0 : (sm_->fluxScale(model)) )) );
1568 : // truncate ggS at ggSMin1
1569 0 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1570 0 : (iif(sm_->ggS(model) < (ggSMin1), ggSMin1*(sm_->fluxScale(model)),
1571 0 : sm_->ggS(model)) )
1572 : );
1573 : }
1574 :
1575 : else{
1576 :
1577 0 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>)
1578 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1579 0 : sqrt((sm_->ggS(model))/ggSMax) )) );
1580 0 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1581 0 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1582 0 : sqrt((sm_->ggS(model))*ggSMax) )) );
1583 :
1584 : }
1585 :
1586 : } else {
1587 : /*
1588 : if(ft_->name() != "MosaicFT"){
1589 : sm_->fluxScale(model).copyData( (LatticeExpr<Float>) 1.0 );
1590 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1591 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1592 : sm_->ggS(model)) ));
1593 :
1594 :
1595 : }
1596 : else{
1597 :
1598 : */
1599 :
1600 :
1601 0 : Int nXX=sm_->ggS(model).shape()(0);
1602 0 : Int nYY=sm_->ggS(model).shape()(1);
1603 0 : Int npola= sm_->ggS(model).shape()(2);
1604 0 : Int nchana= sm_->ggS(model).shape()(3);
1605 0 : IPosition blc(4,nXX, nYY, npola, nchana);
1606 0 : IPosition trc(4, nXX, nYY, npola, nchana);
1607 0 : blc(0)=0; blc(1)=0; trc(0)=nXX-1; trc(1)=nYY-1;
1608 :
1609 :
1610 : //Those damn weights per plane can be wildly different so
1611 : //deal with it properly here
1612 0 : for (Int j=0; j < npola; ++j){
1613 0 : for (Int k=0; k < nchana ; ++k){
1614 :
1615 0 : blc(2)=j; trc(2)=j;
1616 0 : blc(3)=k; trc(3)=k;
1617 0 : Slicer sl(blc, trc, Slicer::endIsLast);
1618 0 : SubImage<Float> fscalesub(sm_->fluxScale(model), sl, true);
1619 0 : SubImage<Float> ggSSub(sm_->ggS(model), sl, true);
1620 : Float planeMax;
1621 0 : LatticeExprNode LEN = max( ggSSub );
1622 0 : planeMax = LEN.getFloat();
1623 0 : if(planeMax !=0){
1624 0 : fscalesub.copyData( (LatticeExpr<Float>)
1625 0 : (iif(ggSSub < (ggSMin2),
1626 0 : 0.0, (ggSSub/planeMax))));
1627 0 : ggSSub.copyData( (LatticeExpr<Float>)
1628 0 : (iif(ggSSub < (ggSMin2), 0.0,
1629 : planeMax) ));
1630 :
1631 :
1632 :
1633 : }
1634 0 : }
1635 :
1636 : }
1637 :
1638 : /*
1639 : ft_->getFluxImage(sm_->fluxScale(model));
1640 :
1641 : sm_->ggS(model).copyData( (LatticeExpr<Float>)
1642 : (iif(sm_->ggS(model) < (ggSMin2), 0.0,
1643 : (sm_->ggS(model)) )) );
1644 : */
1645 : //}
1646 0 : }
1647 :
1648 : //because for usual ft machines a applySJoneInv is done on the gS
1649 : //in the finalizeput stage...need to understand if its necessary
1650 : /*need to understand that square business
1651 : if( (ft_->name() != "MosaicFT") && (!isPSFWork_p)){
1652 : sm_->gS(model).copyData( (LatticeExpr<Float>)
1653 : (iif(sm_->fluxScale(model) > 0.0,
1654 : ((sm_->gS(model))/(sm_->fluxScale(model))), 0.0 )) );
1655 :
1656 :
1657 : }
1658 : */
1659 : ///
1660 0 : sm_->ggS(model).clearCache();
1661 0 : sm_->fluxScale(model).clearCache();
1662 : }
1663 :
1664 : }
1665 0 : }
1666 :
1667 74 : void SkyEquation::checkVisIterNumRows(ROVisibilityIterator& vi){
1668 :
1669 74 : VisBuffer * vb = vi.getVisBuffer();
1670 74 : VisBufferAutoPtr vbap;
1671 74 : if (vb == NULL){
1672 37 : VisBufferAutoPtr tmp (vi);
1673 37 : vbap = tmp;
1674 37 : vb = vbap.get();
1675 37 : }
1676 :
1677 74 : vi.originChunks();
1678 74 : vi.origin();
1679 74 : Int nAnt=vb->numberAnt();
1680 74 : if(nAnt >1){
1681 74 : if (vb->nRow() < (nAnt*(nAnt-1)/4)){
1682 0 : vi.setRowBlocking( nAnt*(nAnt-1)/2+nAnt);
1683 0 : vi.originChunks();
1684 0 : vi.origin();
1685 : }
1686 : }
1687 74 : }
1688 :
1689 0 : String SkyEquation::associatedMSName(){
1690 0 : return String("");
1691 : };
1692 :
1693 0 : void SkyEquation::lock(){
1694 :
1695 : //Do nothing for now
1696 :
1697 0 : }
1698 :
1699 0 : void SkyEquation::unlock(){
1700 : // Do nothing for now
1701 :
1702 0 : }
1703 :
1704 : } //# NAMESPACE CASA - END
1705 :
|
