Line data Source code
1 : //# VisEquation: Vis Equation
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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: VisEquation.cc,v 19.18 2006/01/13 01:06:55 gmoellen Exp $
27 :
28 :
29 : #include <synthesis/MeasurementEquations/VisEquation.h>
30 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
31 : #include <casacore/casa/Arrays/ArrayMath.h>
32 : #include <casacore/casa/Arrays/ArrayPartMath.h>
33 : #include <casacore/casa/Utilities/Assert.h>
34 : #include <casacore/casa/BasicSL/String.h>
35 :
36 : #include <casacore/casa/Exceptions/Error.h>
37 : #include <msvis/MSVis/VisBuffer.h>
38 : //#include <casa/Quanta/MVTime.h>
39 : #include <iostream>
40 :
41 : #include <casacore/casa/OS/Timer.h>
42 :
43 : #define VISEQPRTLEV 0
44 :
45 : using namespace casacore;
46 : namespace casa { //# NAMESPACE CASA - BEGIN
47 :
48 : // ***************************************************************************
49 : // ******************** Start of public member functions ********************
50 : // ***************************************************************************
51 :
52 : //----------------------------------------------------------------------
53 1 : VisEquation::VisEquation() :
54 : vcpb_(NULL),
55 : napp_(0),
56 : lfd_(-1),
57 : rfd_(9999),
58 : freqAveOK_(false),
59 : svc_(NULL),
60 : pivot_(VisCal::ALL), // at the sky
61 : // spwOK_(),
62 : useInternalModel_(false),
63 : nVisTotal_(0),
64 1 : prtlev_(VISEQPRTLEV)
65 : {
66 1 : if (prtlev()>0) cout << "VE::VE()" << endl;
67 1 : };
68 :
69 : //----------------------------------------------------------------------
70 1 : VisEquation::~VisEquation() {};
71 :
72 : //----------------------------------------------------------------------
73 0 : VisEquation& VisEquation::operator=(const VisEquation& other)
74 : {
75 0 : if(this!=&other) {
76 0 : vcpb_=other.vcpb_;
77 0 : napp_=other.napp_;
78 0 : svc_=other.svc_;
79 : };
80 0 : return *this;
81 : };
82 :
83 : //----------------------------------------------------------------------
84 0 : VisEquation::VisEquation(const VisEquation& other)
85 : {
86 0 : operator=(other);
87 0 : }
88 :
89 :
90 : //----------------------------------------------------------------------
91 1 : void VisEquation::setapply(PtrBlock<VisCal*>& vcin) {
92 :
93 1 : if (prtlev()>0) cout << "VE::setapply()" << endl;
94 :
95 : // Be sure internal pointer points to same PB
96 1 : vcpb_ = &vcin;
97 :
98 : // How many VisCals in list?
99 1 : napp_=vc().nelements();
100 :
101 : // only if at least one VisCal in list
102 1 : if (napp_>0) {
103 :
104 : // only sort if a non-trivial list
105 0 : if (napp_>1) {
106 :
107 : // A temporary local copy for sorting:
108 0 : PtrBlock<VisCal*> lvc;
109 0 : lvc.resize(napp_,true,true);
110 0 : lvc=vc();
111 :
112 : // Sorted index will go here
113 0 : Vector<uInt> order(napp_,0);
114 :
115 : // Fill in the sort key
116 0 : Vector<Int> key(napp_,0);
117 0 : for (Int i=0;i<napp_;i++)
118 0 : key(i)=Int(vc()[i]->type());
119 :
120 : // Do the sort
121 : {
122 0 : Sort sort;
123 0 : sort.sortKey(&key(0),TpInt);
124 0 : sort.sort(order,uInt(napp_));
125 : }
126 :
127 : // Assign VisCals in sorted order
128 0 : if (prtlev()>2) cout << "Sorted VisCals:" << endl;
129 0 : vc().set(NULL);
130 0 : for (Int i=0;i<napp_;i++) {
131 0 : vc()[i]=lvc[order(i)];
132 :
133 0 : if (prtlev()>2) cout << vc()[i]->typeName() << " (" << vc()[i]->type() << ")" << endl;
134 0 : vc()[i]->initCalFlagCount();
135 : }
136 :
137 : }
138 : }
139 :
140 : // Set up freq-dependence of the Vis Equation
141 : // (for currently specified solve/apply types)
142 : // TBD: only needed in setsolve?
143 : // setFreqDep();
144 :
145 : // Initialize visibility counter
146 1 : nVisTotal_=0;
147 :
148 1 : }
149 :
150 : //----------------------------------------------------------------------
151 0 : void VisEquation::setsolve(SolvableVisCal& svc) {
152 :
153 0 : if (prtlev()>0) cout << "VE::setsolve()" << endl;
154 :
155 : // VE's own pointer to solvable component
156 0 : svc_=&svc;
157 :
158 : // Set up freq-dependence of the Vis Equation
159 : // (for currently specified solve/apply components)
160 0 : setFreqDep();
161 :
162 0 : }
163 :
164 :
165 : //----------------------------------------------------------------------
166 0 : void VisEquation::setPivot(VisCal::Type pivot) {
167 :
168 0 : if (prtlev()>0) cout << "VE::setPivot()" << endl;
169 :
170 0 : pivot_ = pivot;
171 :
172 0 : }
173 :
174 :
175 : //----------------------------------------------------------------------
176 0 : void VisEquation::setModel(const Vector<Float>& stokes) {
177 :
178 0 : if (prtlev()>0) cout << "VE::setModel()" << endl;
179 :
180 : // Set the internal point source model
181 0 : useInternalModel_=true;
182 0 : stokesModel_.resize(4);
183 0 : stokesModel_.set(0.0);
184 0 : stokesModel_(Slice(0,stokes.nelements(),1))=stokes;
185 :
186 0 : }
187 :
188 : //----------------------------------------------------------------------
189 : // List the terms
190 0 : Vector<VisCal::Type> VisEquation::listTypes() const {
191 :
192 0 : Vector<VisCal::Type> typelist(nTerms());
193 0 : for (Int i=0;i<napp_;++i)
194 0 : typelist(i)=(*vcpb_)[i]->type();
195 :
196 0 : return typelist;
197 :
198 : }
199 :
200 :
201 : //----------------------------------------------------------------------
202 : // Get spwOK for a single spw by aggregating from the vc's
203 0 : Bool VisEquation::spwOK(const Int& spw) {
204 0 : if (napp_<1) return true; // if there are none, all is ok
205 : // otherwise, accumulate from the VisCals
206 0 : Bool spwok=vc()[0]->spwOK(spw);
207 0 : for (Int i=1;i<napp_;++i)
208 0 : spwok = spwok && vc()[i]->spwOK(spw);
209 0 : return spwok;
210 : }
211 :
212 :
213 : //----------------------------------------------------------------------
214 : // Report if calibration is collectively calibrateable by all VCs
215 : // (including possible agnosticism by somein CalLibrary context;
216 : // see SolvableVisCal::VBOKforCalApply)
217 0 : Bool VisEquation::VBOKforCalApply(vi::VisBuffer2& vb) {
218 :
219 0 : Bool okForCal(True); // nominal
220 0 : for (Int i=0;i<napp_;++i)
221 0 : okForCal = okForCal && vc()[i]->VBOKforCalApply(vb);
222 :
223 0 : return okForCal;
224 :
225 : }
226 :
227 :
228 :
229 : //----------------------------------------------------------------------
230 : // Correct in place the OBSERVED visibilities in a VisBuffer
231 0 : void VisEquation::correct(VisBuffer& vb, Bool trial) {
232 :
233 0 : if (prtlev()>0) cout << "VE::correct()" << endl;
234 :
235 0 : AlwaysAssert(ok(),AipsError);
236 :
237 0 : if (napp_==0) throw(AipsError("Nothing to Apply"));
238 :
239 : // Ensure correlations in canonical order
240 : // (this is a no-op if no sort necessary)
241 0 : vb.sortCorr();
242 :
243 : // Accumulate visibility count
244 0 : nVisTotal_+=(vb.nCorr()*vb.nChannel()*vb.nRow());
245 :
246 : // Apply each VisCal in left-to-right order
247 0 : for (Int iapp=0;iapp<napp_;iapp++)
248 0 : vc()[iapp]->correct(vb,trial);
249 :
250 : // Ensure correlations restored to original order
251 : // (this is a no-op if no sort necessary)
252 0 : vb.unSortCorr();
253 :
254 0 : }
255 : //----------------------------------------------------------------------
256 : // Correct in place the OBSERVED visibilities in a VisBuffer
257 0 : void VisEquation::correct2(vi::VisBuffer2& vb, Bool trial, Bool doWtSp) {
258 :
259 0 : if (prtlev()>0) cout << "VE::correct2()" << endl;
260 :
261 0 : AlwaysAssert(ok(),AipsError);
262 :
263 0 : if (napp_==0) throw(AipsError("Nothing to Apply"));
264 :
265 : // Accumulate visibility count
266 0 : nVisTotal_+=(vb.nCorrelations()*vb.nChannels()*vb.nRows());
267 :
268 : // Apply each VisCal in left-to-right order
269 0 : for (Int iapp=0;iapp<napp_;iapp++)
270 0 : vc()[iapp]->correct2(vb,trial,doWtSp);
271 :
272 0 : }
273 :
274 : //----------------------------------------------------------------------
275 : // Report action record info (derived from consituent VisCals
276 0 : Record VisEquation::actionRec() {
277 0 : Record cf;
278 : // Add in each VisCal's record
279 0 : for (Int iapp=0;iapp<napp_;iapp++)
280 0 : cf.defineRecord(iapp,vc()[iapp]->actionRec());
281 :
282 : // The TOTAL number of visibilities that passed through the VisEquationp
283 0 : cf.define("nVisTotal",nVisTotal_);
284 :
285 0 : return cf;
286 : }
287 :
288 : //----------------------------------------------------------------------
289 : // Corrupt in place the MODEL visibilities in a VisBuffer
290 0 : void VisEquation::corrupt(VisBuffer& vb) {
291 :
292 0 : if (prtlev()>0) cout << "VE::corrupt()" << endl;
293 :
294 0 : AlwaysAssert(ok(),AipsError);
295 :
296 0 : if (napp_==0) throw(AipsError("Nothing to Apply"));
297 :
298 : // Ensure correlations in canonical order
299 : // (this is a no-op if no sort necessary)
300 0 : vb.sortCorr();
301 :
302 : // Apply each VisCal in right-to-left order
303 0 : for (Int iapp=(napp_-1);iapp>-1;iapp--)
304 0 : vc()[iapp]->corrupt(vb);
305 :
306 : // Ensure correlations restored to original order
307 : // (this is a no-op if no sort necessary)
308 0 : vb.unSortCorr();
309 :
310 0 : }
311 :
312 : //----------------------------------------------------------------------
313 : // Corrupt in place the MODEL visibilities in a VisBuffer
314 0 : void VisEquation::corrupt2(vi::VisBuffer2& vb) {
315 :
316 0 : if (prtlev()>0) cout << "VE::corrupt2(VB2)" << endl;
317 :
318 0 : AlwaysAssert(ok(),AipsError);
319 :
320 0 : if (napp_==0) throw(AipsError("Nothing to Apply"));
321 :
322 : // Apply each VisCal in right-to-left order
323 0 : for (Int iapp=(napp_-1);iapp>-1;iapp--)
324 0 : vc()[iapp]->corrupt2(vb);
325 :
326 0 : }
327 :
328 : //----------------------------------------------------------------------
329 0 : void VisEquation::collapse(VisBuffer& vb) {
330 :
331 0 : if (prtlev()>0) cout << "VE::collapse()" << endl;
332 :
333 : // Handle origin of model data here:
334 0 : if (useInternalModel_)
335 0 : vb.setModelVisCube(stokesModel_);
336 : else
337 0 : vb.modelVisCube();
338 :
339 : // Ensure required columns are present!
340 0 : vb.visCube();
341 0 : vb.weightMat();
342 :
343 : // Re-calculate weights from sigma column
344 : // TBD: somehow avoid if not necessary?
345 0 : vb.resetWeightMat();
346 :
347 : // Ensure correlations in canonical order
348 : // (this is a no-op if no sort necessary)
349 : // TBD: optimize in combo with model origination?
350 0 : vb.sortCorr();
351 :
352 : // If we are solving for the polarization:
353 : // 1. Normalize data and model by I model
354 : // 2. Set cross-hands to (1,0) so P factors multiplying them
355 : // are propogated, and we can solve for pol factors
356 0 : if (svc().solvePol())
357 0 : svc().setUpForPolSolve(vb);
358 :
359 : // TBD: When warranted, do freqAve before setUpForPolSolve?
360 :
361 : // initialize LHS/RHS indices
362 0 : Int lidx=0;
363 0 : Int ridx=napp_-1;
364 :
365 : // If solve NOT freqDep, and data is, we want
366 : // to freqAve as soon as possible before solve;
367 : // apply any freqDep cal first
368 0 : if ( freqAveOK_ && !svc().freqDepMat() && vb.nChannel()>1 ) {
369 :
370 : // Correct OBSERVED data up to last freqDep term on LHS
371 : // (type(lfd_) guaranteed < type(svc))
372 0 : while (lidx<napp_ && lidx <= lfd_) {
373 0 : vc()[lidx]->correct(vb);
374 0 : lidx++;
375 : }
376 :
377 : // Corrupt MODEL data down to last freqDep term on RHS
378 : // (type(rfd_) guaranteed >= type(svc))
379 0 : while (ridx>-1 && ridx >= rfd_) {
380 0 : vc()[ridx]->corrupt(vb);
381 0 : ridx--;
382 : }
383 :
384 : // All freq-dep cal has been applied, now freqAve
385 0 : vb.freqAveCubes(); // BOTH sides (if present)!
386 :
387 : }
388 :
389 : // Correct DATA up to solved-for term
390 0 : while (lidx<napp_ && vc()[lidx]->type() < svc().type()) {
391 0 : vc()[lidx]->correct(vb);
392 0 : lidx++;
393 : }
394 :
395 : // Corrupt MODEL down to solved-for term (incl. same type as solved-for term)
396 0 : while (ridx>-1 && vc()[ridx]->type() >= svc().type()) {
397 0 : vc()[ridx]->corrupt(vb);
398 0 : ridx--;
399 : }
400 :
401 0 : }
402 :
403 : //----------------------------------------------------------------------
404 0 : void VisEquation::diffModelStokes(vi::VisBuffer2& vb, std::map<String,Cube<Complex> > dMdS) {
405 :
406 0 : Int nCorr(vb.nCorrelations());
407 0 : if (nCorr<4)
408 0 : throw(AipsError("Cannot differentiate w.r.t. Model Stokes unless data has four correlations"));
409 :
410 0 : Int nChan(vb.nChannels());
411 0 : Int nRow(vb.nRows());
412 :
413 : // Incoming map should be empty
414 0 : dMdS.clear();
415 :
416 : // Basis-dependent indexing
417 0 : Slice Isl(0,2,3); // not basis-specific
418 0 : Slice Qsl0,Qsl1,Usl0,Usl1,Vsl0,Vsl1;
419 0 : Bool doCirc(true);
420 0 : if (vb.correlationTypes()(0)==5) {
421 : // Circular
422 0 : doCirc=true;
423 0 : Qsl0=Slice(1); // cross-hand
424 0 : Qsl1=Slice(2);
425 0 : Usl0=Slice(1); // cross-hand
426 0 : Usl1=Slice(2);
427 0 : Vsl0=Slice(0); // parallel-hand
428 0 : Vsl1=Slice(3);
429 : }
430 0 : else if (vb.correlationTypes()(0)==9) {
431 : // Linear
432 0 : doCirc=false;
433 0 : Qsl0=Slice(0); // parallel-hand
434 0 : Qsl1=Slice(3);
435 0 : Usl0=Slice(1); // cross-hand
436 0 : Usl1=Slice(2);
437 0 : Vsl0=Slice(1); // cross-hand
438 0 : Vsl1=Slice(2);
439 : }
440 : else
441 0 : throw(AipsError("Cannot differentiate w.r.t. Model Stokes for unrecognized polarization basis"));
442 :
443 0 : Complex cOne(1,0), cIm(0,1);
444 :
445 0 : String I("I");
446 0 : dMdS[I]=Cube<Complex>(nCorr,nChan,nRow,0.0);
447 0 : Cube<Complex>& dMdI(dMdS[I]);
448 0 : dMdI(Isl,Slice(),Slice()).set(cOne); // not basis-specific
449 :
450 0 : String Q("Q");
451 0 : dMdS[Q]=Cube<Complex>(nCorr,nChan,nRow,0.0);
452 0 : Cube<Complex>& dMdQ(dMdS[Q]);
453 0 : dMdQ(Qsl0,Slice(),Slice()).set(cOne);
454 0 : dMdQ(Qsl1,Slice(),Slice()).set( (doCirc ? cOne : -cOne) );
455 :
456 0 : String U("U");
457 0 : dMdS[U]=Cube<Complex>(nCorr,nChan,nRow,0.0);
458 0 : Cube<Complex>& dMdU(dMdS[U]);
459 0 : dMdU(Usl0,Slice(),Slice()).set( (doCirc ? cIm : cOne) );
460 0 : dMdU(Usl1,Slice(),Slice()).set( (doCirc ? -cIm : cOne) );
461 :
462 0 : String V("V");
463 0 : dMdS[V]=Cube<Complex>(nCorr,nChan,nRow,0.0);
464 0 : Cube<Complex>& dMdV(dMdS[V]);
465 0 : dMdV(Vsl0,Slice(),Slice()).set( (doCirc ? cOne : cIm) );
466 0 : dMdV(Vsl1,Slice(),Slice()).set( (doCirc ? -cOne : -cIm) );
467 :
468 0 : Int ridx=napp_-1;
469 :
470 : // Corrupt MODEL down to solved-for term (incl. same type as solved-for term)
471 0 : while (ridx>-1 && vc()[ridx]->type() >= svc().type()) {
472 0 : vc()[ridx]->corrupt2(vb,dMdI);
473 0 : vc()[ridx]->corrupt2(vb,dMdQ);
474 0 : vc()[ridx]->corrupt2(vb,dMdU);
475 0 : vc()[ridx]->corrupt2(vb,dMdV);
476 0 : ridx--;
477 : }
478 :
479 0 : }
480 :
481 :
482 :
483 :
484 :
485 : //----------------------------------------------------------------------
486 0 : void VisEquation::collapse2(vi::VisBuffer2& vb) {
487 :
488 0 : if (prtlev()>0) cout << "VE::collapse2(VB2)" << endl;
489 :
490 : // Trap case of unavailable calibration in any vc we intend to apply below
491 : // In the solve context, if we can't pre-cal, we flag it
492 : // NB: this assumes only one spw in the VB2!
493 : //if (!this->spwOK(vb.spectralWindows()(0))) {
494 : // Use new VBOKforCalApply, which is f(obs,fld,intent,spw) (not just f(spw))
495 0 : if (!this->VBOKforCalApply(vb)) {
496 :
497 : //cout << "UNCALIBRATEABLE VB2 in VE::collapse2" << endl;
498 :
499 0 : Cube<Bool> fl(vb.flagCube()); fl.set(true);
500 0 : Cube<Float> wtsp(vb.weightSpectrum()); wtsp.set(0.0f);
501 0 : Matrix<Float> wt(vb.weight()); wt.set(0.0f);
502 0 : return;
503 : }
504 :
505 : // Handle origin of model data here:
506 0 : if (useInternalModel_)
507 : // Use specified (point-source) stokes model
508 0 : vb.setVisCubeModel(stokesModel_);
509 : else
510 : // from MS
511 0 : vb.visCubeModel();
512 :
513 : // If we are solving for the polarization:
514 : // 1. Normalize data and model by I model
515 : // 2. Set cross-hands to (1,0) so P factors multiplying them
516 : // are propogated, and we can solve for pol factors
517 0 : if (svc().solvePol())
518 0 : svc().setUpForPolSolve(vb);
519 :
520 : // initialize LHS/RHS indices
521 0 : Int lidx=0;
522 0 : Int ridx=napp_-1;
523 :
524 : // Correct DATA up to solved-for term
525 0 : while (lidx<napp_ && vc()[lidx]->type() < svc().type()) {
526 0 : vc()[lidx]->correct2(vb,False,True);
527 0 : lidx++;
528 : }
529 :
530 : // Corrupt MODEL down to solved-for term (incl. same type as solved-for term)
531 0 : while (ridx>-1 && vc()[ridx]->type() >= svc().type()) {
532 0 : vc()[ridx]->corrupt2(vb);
533 0 : ridx--;
534 : }
535 :
536 0 : if (svc().normalizable())
537 0 : divideCorrByModel(vb);
538 :
539 :
540 : // Make fractional visibilities, if appropriate
541 : // (e.g., for some polarization calibration solves on point-like calibrators)
542 0 : if (svc().divideByStokesIModelForSolve())
543 0 : divideByStokesIModel(vb);
544 :
545 :
546 : }
547 :
548 0 : void VisEquation::divideCorrByModel(vi::VisBuffer2& vb) {
549 :
550 : // This divides corrected data by model
551 : // ... and updates weightspec accordingly
552 :
553 0 : Cube<Complex> c(vb.visCubeCorrected());
554 0 : Cube<Complex> m(vb.visCubeModel());
555 0 : Cube<Bool> fl(vb.flagCube());
556 0 : Cube<Float> w(vb.weightSpectrum());
557 :
558 0 : Complex cOne(1.0);
559 :
560 0 : for (rownr_t irow=0;irow<vb.nRows();++irow) {
561 0 : if (vb.flagRow()(irow)) {
562 : // Row flagged, make sure cube also flagged, weight/data zeroed
563 0 : c(Slice(),Slice(),Slice(irow,1,1))=0.0f;
564 0 : w(Slice(),Slice(),Slice(irow,1,1))=0.0f;
565 0 : fl(Slice(),Slice(),Slice(irow,1,1))=True;
566 : }
567 : else {
568 0 : Bool *flp=&fl(0,0,irow);
569 0 : Float *wtp=&w(0,0,irow);
570 0 : Complex *cvp=&c(0,0,irow);
571 0 : Complex *mvp=&m(0,0,irow);
572 :
573 0 : for (Int ichan=0;ichan<vb.nChannels();++ichan) {
574 0 : for (Int icorr=0;icorr<vb.nCorrelations();++icorr) {
575 0 : Bool& Fl(*flp);
576 0 : Float& W(*wtp);
577 : //Bool& Fl(fl(icorr,ichan,irow));
578 : //Float& W(w(icorr,ichan,irow));
579 0 : if (!Fl) {
580 : // Not flagged...
581 0 : Float A=abs(*mvp);
582 : //Float A=abs(m(icorr,ichan,irow));
583 0 : if (A >0.0f) {
584 : // ...and model non-zero
585 0 : Complex& C(*cvp);
586 0 : Complex& M(*mvp);
587 : //Complex& C(c(icorr,ichan,irow));
588 : //Complex& M(m(icorr,ichan,irow));
589 :
590 : // divide corr data by model
591 : // NB: Use of DComplex here increased cost of this calculation by ~33%
592 0 : C=Complex(DComplex(C)/DComplex(M));
593 : //C=C/M;
594 :
595 0 : W*=square(A); // multiply weight by model**2
596 0 : M=cOne; // divide model by itself
597 :
598 : }
599 : }
600 : else {
601 : // be sure it is flagged and weightless
602 0 : Fl=True;
603 0 : W=0.0f;
604 : }
605 0 : ++cvp;
606 0 : ++mvp;
607 0 : ++flp;
608 0 : ++wtp;
609 : } // icorr
610 : } // ichan
611 : } // !flagRow
612 : } // irow
613 :
614 : // Set unchan'd weight, in case someone wants it
615 : // NB: Use of median increases cost by ~100%
616 : // NB: use of mean increases cost by ~50%
617 : // ...but both are inaccurate if some channels flagged,
618 : // and it should not be necessary to do this here
619 0 : vb.setWeight(partialMedians(vb.weightSpectrum(),IPosition(1,1),True));
620 : //vb.setWeight(partialMeans(vb.weightSpectrum(),IPosition(1,1)));
621 :
622 0 : }
623 :
624 0 : void VisEquation::divideByStokesIModel(vi::VisBuffer2& vb) {
625 :
626 0 : Int nCorr(vb.nCorrelations());
627 :
628 : // This divides corrected data and model by the Stokes I model
629 : // ... and updates weightspec accordingly
630 :
631 0 : Cube<Complex> c(vb.visCubeCorrected());
632 0 : Cube<Complex> m(vb.visCubeModel());
633 0 : Cube<Bool> fl(vb.flagCube());
634 0 : Cube<Float> w(vb.weightSpectrum());
635 :
636 : //Complex cOne(1.0);
637 :
638 0 : for (rownr_t irow=0;irow<vb.nRows();++irow) {
639 0 : if (vb.flagRow()(irow)) {
640 : // Row flagged, make sure cube also flagged, weight/data zeroed
641 0 : c(Slice(),Slice(),Slice(irow,1,1))=0.0f;
642 0 : w(Slice(),Slice(),Slice(irow,1,1))=0.0f;
643 0 : fl(Slice(),Slice(),Slice(irow,1,1))=True;
644 : }
645 : else {
646 : // Bool *flp=&fl(0,0,irow);
647 0 : Float *wtp=&w(0,0,irow);
648 0 : Complex *cvp=&c(0,0,irow);
649 0 : Complex *mvp=&m(0,0,irow);
650 :
651 0 : for (Int ichan=0;ichan<vb.nChannels();++ichan) {
652 0 : Complex Imod(0.0);
653 0 : Float Iamp2(1.0);
654 0 : if (!fl(0,ichan,irow) && !fl(nCorr-1,ichan,irow)) {
655 0 : Imod=(m(0,ichan,irow) + m(nCorr-1,ichan,irow))/2.0f;
656 0 : Iamp2=real(Imod*conj(Imod)); // squared model amp (for weight adjust)
657 0 : if (Iamp2>0.0f) {
658 0 : for (Int icorr=0;icorr<nCorr;++icorr) {
659 0 : Float& W(*wtp);
660 0 : Complex& C(*cvp);
661 0 : Complex& M(*mvp);
662 0 : C/=Imod;
663 0 : M/=Imod;
664 0 : W*=Iamp2;
665 0 : ++cvp;
666 0 : ++mvp;
667 0 : ++wtp;
668 : } // icorr
669 : } // non-zero Imod
670 : } // parallel hands not flagged
671 : } // ichan
672 : } // !flagRow
673 : } // irow
674 :
675 : // Set unchan'd weight, in case someone wants it
676 : // NB: Use of median increases cost by ~100%
677 : // NB: use of mean increases cost by ~50%
678 : // ...but both are inaccurate if some channels flagged,
679 : // and it should not be necessary to do this here
680 0 : vb.setWeight(partialMedians(vb.weightSpectrum(),IPosition(1,1),True));
681 : //vb.setWeight(partialMeans(vb.weightSpectrum(),IPosition(1,1)));
682 :
683 0 : }
684 :
685 :
686 :
687 :
688 : //----------------------------------------------------------------------
689 : // void VisEquation::collapseForSim(VisBuffer& vb) {
690 0 : void VisEquation::collapseForSim(VisBuffer& vb) {
691 :
692 0 : if (prtlev()>0) cout << "VE::collapseforSim()" << endl;
693 :
694 : // Handle origin of model data here (?):
695 :
696 : // Ensure correlations in canonical order
697 : // (this is a no-op if no sort necessary)
698 : // TBD: optimize in combo with model origination?
699 0 : vb.sortCorr();
700 :
701 : // initialize LHS/RHS indices
702 0 : Int lidx=0;
703 0 : Int ridx=napp_-1;
704 :
705 : #ifdef RI_DEBUG
706 : cout << "vb.visCube original: " << vb.visCube()(0,0,500) << vb.visCube()(0,0,1216) << vb.visCube()(0,0,1224) << endl;
707 : cout << "vb.model original: " << vb.modelVisCube()(0,0,500) << vb.modelVisCube()(0,0,1216) << vb.modelVisCube()(0,0,1224) << endl;
708 : #endif
709 :
710 : // copy data to model, to be corrupted in place there.
711 : // 20091030 RI changed skyequation to use Observed. the below
712 : // should not require scratch columns
713 0 : vb.setModelVisCube(vb.visCube());
714 :
715 : // zero the data. correct will operate in place on data, so
716 : // if we don't have an AMueller we don't get anything from this.
717 0 : vb.setVisCube(0.0);
718 :
719 : #ifdef RI_DEBUG
720 : cout << "vb.visCube before crct: " << vb.visCube()(0,0,500) << vb.visCube()(0,0,1216) << vb.visCube()(0,0,1224) << endl;
721 : #endif
722 :
723 : // Correct DATA up to pivot
724 0 : while (lidx<napp_ && vc()[lidx]->type() < pivot_) {
725 0 : if (prtlev()>2) cout << vc()[lidx]->typeName();
726 0 : if (vc()[ridx]->extraTag()!="NoiseScale" or vc()[lidx]->type()!=VisCal::T) {
727 0 : vc()[lidx]->correct(vb);
728 0 : if (prtlev()>2) cout << " -> correct";
729 : }
730 0 : if (prtlev()>2) cout << endl;
731 0 : lidx++;
732 : }
733 :
734 : #ifdef RI_DEBUG
735 : cout << "vb.visCube after crct: " << vb.visCube()(0,0,500) << vb.visCube()(0,0,1216) << vb.visCube()(0,0,1224) << endl;
736 : cout << "vb.model before crpt: " << vb.modelVisCube()(0,0,500) << vb.modelVisCube()(0,0,1216) << vb.modelVisCube()(0,0,1224) << endl;
737 : #endif
738 :
739 : // Corrupt Model down to (and including) the pivot
740 0 : while (ridx>-1 && vc()[ridx]->type() >= pivot_) {
741 0 : if (prtlev()>2) cout << vc()[lidx]->typeName();
742 : // manually pick off a T intended to be noise scaling T:
743 0 : if (pivot_ <= VisCal::T and vc()[ridx]->type()==VisCal::T) {
744 0 : if (vc()[ridx]->extraTag()=="NoiseScale") {
745 0 : vc()[ridx]->correct(vb); // correct DATA
746 0 : if (prtlev()>2) cout << " -> correct";
747 : } else {
748 0 : vc()[ridx]->corrupt(vb);
749 0 : if (prtlev()>2) cout << " -> corrupt";
750 : }
751 : } else {
752 0 : vc()[ridx]->corrupt(vb);
753 0 : if (prtlev()>2) cout << " -> corrupt";
754 : }
755 0 : if (prtlev()>2) cout << endl;
756 0 : ridx--;
757 : }
758 :
759 : #ifdef RI_DEBUG
760 : cout << "vb.model after crpt: " << vb.modelVisCube()(0,0,500) << vb.modelVisCube()(0,0,1216) << vb.modelVisCube()(0,0,1224) << endl;
761 : cout << "vb.visCube after crpt: " << vb.visCube()(0,0,500) << vb.visCube()(0,0,1216) << vb.visCube()(0,0,1224) << endl << endl;
762 : #endif
763 :
764 : // add corrected/scaled data (e.g. noise) to corrupted model
765 : // vb.modelVisCube()+=vb.visCube();
766 :
767 : // add corrupted Model to corrected/scaled data (i.e.. noise)
768 0 : vb.visCube()+=vb.modelVisCube();
769 :
770 : // Put correlations back in original order
771 : // (~no-op if already canonical)
772 0 : vb.unSortCorr();
773 :
774 0 : }
775 :
776 0 : void VisEquation::state() {
777 :
778 0 : if (prtlev()>0) cout << "VE::state()" << endl;
779 :
780 0 : cout << boolalpha;
781 :
782 : // Order in which DATA is corrected
783 0 : cout << "Correct order:" << endl;
784 0 : if (napp_>0) {
785 0 : for (Int iapp=0;iapp<napp_;iapp++)
786 0 : cout << " " << iapp << " "
787 0 : << vc()[iapp]->typeName() << " ("
788 0 : << vc()[iapp]->type() << ")" << endl;
789 : }
790 : else
791 0 : cout << " <none>" << endl;
792 :
793 0 : cout << endl;
794 :
795 0 : cout << "Corrupt order:" << endl;
796 0 : if (napp_>0) {
797 0 : for (Int iapp=(napp_-1);iapp>-1;iapp--)
798 0 : cout << " " << iapp << " "
799 0 : << vc()[iapp]->typeName() << " ("
800 0 : << vc()[iapp]->type() << ")" << endl;
801 : }
802 : else
803 0 : cout << " <none>" << endl;
804 :
805 0 : cout << endl;
806 :
807 0 : cout << "Source model:" << endl;
808 0 : cout << " useInternalModel_ = " << useInternalModel_ << endl;
809 0 : if (useInternalModel_)
810 0 : cout << " Stokes = " << stokesModel_ << endl;
811 : else
812 0 : cout << " <using MODEL_DATA column>" << endl;
813 0 : cout << endl;
814 :
815 0 : cout << "Collapse order:" << endl;
816 0 : cout << " freqAveOK_ = " << freqAveOK_ << endl;
817 :
818 0 : if (svc_) {
819 0 : Int lidx=0;
820 0 : Int ridx=napp_-1;
821 :
822 0 : if ( freqAveOK_ && !svc().freqDepMat() ) {
823 : // Correct DATA up to last freqDep term on LHS
824 0 : cout << " LHS (pre-freqAve):" << endl;
825 0 : if (lidx <= lfd_)
826 0 : while (lidx <= lfd_) { // up to and including the lfd_th term
827 0 : cout << " " << lidx << " "
828 0 : << vc()[lidx]->typeName() << " ("
829 0 : << vc()[lidx]->type() << ")"
830 0 : << " (freqDep correct)"
831 0 : << endl;
832 0 : lidx++;
833 : }
834 : else
835 0 : cout << " <none>" << endl;
836 :
837 : // Corrupt MODEL down to last freqDep term on RHS
838 0 : cout << " RHS (pre-freqAve):" << endl;
839 0 : if (ridx >= rfd_)
840 0 : while (ridx >= rfd_) { // down to and including the rfd_th term
841 0 : cout << " " << (ridx) << " "
842 0 : << vc()[ridx]->typeName() << " ("
843 0 : << vc()[ridx]->type() << ")"
844 0 : << " (freqDep corrupt)"
845 0 : << endl;
846 0 : ridx--;
847 : }
848 : else
849 0 : cout << " <none>" << endl;
850 :
851 0 : cout << " Frequency average both sides" << endl;
852 : }
853 :
854 0 : cout << " LHS:" << endl;
855 0 : if (lidx<napp_ && vc()[lidx]->type() < svc().type())
856 0 : while (lidx<napp_ && vc()[lidx]->type() < svc().type()) {
857 0 : cout << " " << (lidx) << " "
858 0 : << vc()[lidx]->typeName() << " ("
859 0 : << vc()[lidx]->type() << ")"
860 0 : << endl;
861 0 : lidx++;
862 : }
863 : else
864 0 : cout << " <none>" << endl;
865 :
866 0 : cout << " RHS:" << endl;
867 0 : if (ridx>-1 && vc()[ridx]->type() >= svc().type())
868 0 : while (ridx>-1 && vc()[ridx]->type() >= svc().type()) {
869 0 : cout << " " << (ridx) << " "
870 0 : << vc()[ridx]->typeName() << " ("
871 0 : << vc()[ridx]->type() << ")"
872 0 : << endl;
873 0 : ridx--;
874 : }
875 : else
876 0 : cout << " <none>" << endl;
877 :
878 :
879 0 : cout << "Ready to solve for " << svc().typeName() << " (" << svc().type() << ")" << endl;
880 :
881 : }
882 : else
883 0 : cout << " <Nothing to solve for>" << endl;
884 :
885 0 : }
886 :
887 :
888 :
889 : //----------------------------------------------------------------------
890 : // Determine residuals of VisBuffer data w.r.t. svc_
891 0 : void VisEquation::residuals(VisBuffer& vb,
892 : Cube<Complex>& R,
893 : const Int chan) {
894 :
895 :
896 0 : if (prtlev()>3) cout << "VE::residuals()" << endl;
897 0 : if (prtlev()>13) // Here only to shush a
898 0 : cout << "vb.nRow(): " << vb.nRow() // compiler warning about
899 0 : << "\nR.shape(): " << R.shape() // unused variables.
900 0 : << "\nchan: " << chan
901 0 : << endl;
902 :
903 : // Trap unspecified solvable term
904 0 : if (!svc_)
905 0 : throw(AipsError("No calibration term to differentiate."));
906 :
907 : // TBD: support for public non-in-place corruption and
908 : // corrupt for specific channel
909 :
910 0 : throw(AipsError("Simple residual calculation NYI."));
911 :
912 : }
913 :
914 :
915 : // Calculate residuals and differentiated residuals
916 0 : void VisEquation::diffResiduals(CalVisBuffer& cvb) {
917 :
918 0 : if (prtlev()>3) cout << "VE::diffResiduals(CVB)" << endl;
919 :
920 : // Trap unspecified solvable term
921 0 : if (!svc_)
922 0 : throw(AipsError("No calibration term to differentiate."));
923 :
924 : // Delegate to the SVC
925 0 : svc().differentiate(cvb);
926 :
927 0 : }
928 :
929 : // Calculate residuals and differentiated residuals
930 0 : void VisEquation::differentiate(SolveDataBuffer& sdb) {
931 :
932 0 : if (prtlev()>3) cout << "VE::differentiate(SDB)" << endl;
933 :
934 : // Trap unspecified solvable term
935 0 : if (!svc_)
936 0 : throw(AipsError("No calibration term to differentiate."));
937 :
938 : // Delegate to the SVC
939 0 : svc().differentiate(sdb);
940 :
941 0 : }
942 :
943 :
944 : // Calculate residuals and differentiated residuals
945 0 : void VisEquation::diffResiduals(VisBuffer& vb,
946 : Cube<Complex>& R,
947 : Array<Complex>& dR,
948 : Matrix<Bool>& Rflg) {
949 :
950 0 : if (prtlev()>3) cout << "VE::diffResiduals()" << endl;
951 :
952 : // Trap unspecified solvable term
953 0 : if (!svc_)
954 0 : throw(AipsError("No calibration term to differentiate."));
955 :
956 : // Get trial model corrupt and differentiation from the SVC
957 : // (R, dR, and Rflg will be sized by svc_)
958 0 : svc().differentiate(vb,R,dR,Rflg);
959 :
960 : // Now, subtract obs'd data from trial model corruption
961 :
962 : // Shape info
963 0 : Int nCorr(vb.corrType().nelements());
964 0 : Int& nChan(vb.nChannel());
965 0 : Int& nRow(vb.nRow());
966 0 : IPosition blc(3,0,0,0);
967 0 : IPosition trc(3,nCorr-1,nChan-1,nRow-1);
968 :
969 : // Slice if specific channel requested
970 0 : if (svc().freqDepPar())
971 0 : blc(1)=trc(1)=svc().focusChan();
972 :
973 : // If shapes match, subtract, else throw
974 : // TBD: avoid subtraction in flagged channels?
975 0 : Cube<Complex> Vo(vb.visCube()(blc,trc));
976 0 : if (R.shape()==Vo.shape())
977 0 : R-=Vo;
978 : else
979 0 : throw(AipsError("Shape mismatch in residual calculation"));
980 :
981 0 : }
982 :
983 0 : void VisEquation::diffResiduals(VisBuffer& R,
984 : VisBuffer& dR0,
985 : VisBuffer& dR1,
986 : Matrix<Bool>& Rflg)
987 : {
988 0 : if (prtlev()>3) cout << "VE::diffResiduals()" << endl;
989 :
990 : // Trap unspecified solvable term
991 0 : if (!svc_)
992 0 : throw(AipsError("No calibration term to differentiate."));
993 :
994 : // Get trial model corrupt and differentiation from the SVC
995 : // (R, dR, and Rflg will be sized by svc_)
996 0 : svc().differentiate(R,dR0,dR1, Rflg);
997 :
998 : // Now, subtract obs'd data from trial model corruption
999 :
1000 : // Shape info
1001 0 : Int nCorr(R.corrType().nelements());
1002 0 : Int& nChan(R.nChannel());
1003 0 : Int& nRow(R.nRow());
1004 0 : IPosition blc(3,0,0,0);
1005 0 : IPosition trc(3,nCorr-1,nChan-1,nRow-1);
1006 :
1007 : // Slice if specific channel requested
1008 0 : if (svc().freqDepPar())
1009 0 : blc(1)=trc(1)=svc().focusChan();
1010 : // If shapes match, subtract, else throw
1011 : // TBD: avoid subtraction in flagged channels?
1012 : // Cube<Complex> Vo(R.visCube()(blc,trc));
1013 : // cout << R.correctedVisCube().shape() << " "
1014 : // << R.visCube().shape() << " "
1015 : // << blc << " " << trc << endl;
1016 : // Cube<Complex> Vo(R.correctedVisCube()(blc,trc));
1017 0 : Cube<Complex> Vo(R.visCube()(blc,trc));
1018 0 : Cube<Complex> Res;Res.reference(R.modelVisCube());
1019 : /*
1020 : for(Int i=0;i<Res.shape()(2);i++)
1021 : {
1022 : cout << i
1023 : << " " << Res(0,0,i)
1024 : << " " << Res(3,0,i)
1025 : << " " << R.visCube()(0,0,i)
1026 : << " " << R.visCube()(3,0,i)
1027 : << " " << R.flag()(0,i)
1028 : << " " << R.flag()(3,i)
1029 : << " " << R.antenna1()(i)
1030 : << " " << R.antenna2()(i)
1031 : << " " << R.flagRow()(i)
1032 : << endl;
1033 : }
1034 : exit(0);
1035 : */
1036 0 : if (Res.shape()==Vo.shape())
1037 0 : Res-=Vo;
1038 : else
1039 0 : throw(AipsError("Shape mismatch in residual calculation"));
1040 0 : }
1041 :
1042 :
1043 : //----------------------------------------------------------------------
1044 : // ***************************************************************************
1045 : // ******************** Start of protected member functions *****************
1046 : // ***************************************************************************
1047 :
1048 :
1049 0 : void VisEquation::setFreqDep() {
1050 :
1051 0 : if (prtlev()>2) cout << "VE::setFreqDep()" << endl;
1052 :
1053 : // Nominal freq-dep is none on either side
1054 0 : lfd_=-1; // right-most freq-dep term on LHS
1055 0 : rfd_=napp_; // left-most freq-dep term on RHS
1056 :
1057 : // Nominally averaging in frequency before normalization is NOT OK
1058 : // (we will revise this when we can assert constant MODEL_DATA)
1059 0 : freqAveOK_=false;
1060 :
1061 : // Only if there are both apply-able and solve-able terms
1062 0 : if (svc_ && napp_>0) {
1063 :
1064 : // freqdep to LEFT of solvable type
1065 0 : for (Int idx=0; (idx<napp_ && vc()[idx]->type()<svc().type()); idx++)
1066 0 : if (vc()[idx]->freqDepMat()) lfd_=idx;
1067 :
1068 : // freqdep to RIGHT of solvable type
1069 0 : for (Int idx=(napp_-1); (idx>-1 && vc()[idx]->type()>=svc().type()); idx--)
1070 0 : if (vc()[idx]->freqDepMat()) {
1071 0 : rfd_=idx;
1072 : // If we will corrupt the model with something freqdep, we can't
1073 : // frequency average in collapse
1074 0 : freqAveOK_=false;
1075 : }
1076 :
1077 : }
1078 :
1079 :
1080 0 : }
1081 :
1082 0 : Bool VisEquation::ok() {
1083 :
1084 0 : if (napp_ != Int(vc().nelements()))
1085 0 : return false;
1086 :
1087 0 : return(true);
1088 : }
1089 :
1090 :
1091 : } //# NAMESPACE CASA - END
1092 :
|
