Line data Source code
1 : //# CalibratingVi2.cc: Implementation of the CalibratingVi2 class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <synthesis/MeasurementComponents/CalibratingVi2.h>
24 : #include <synthesis/MeasurementComponents/Calibrater.h>
25 : #include <synthesis/CalLibrary/CalLibraryTools.h>
26 : #include <casacore/casa/Arrays/ArrayPartMath.h>
27 : #include <casacore/casa/Arrays/MaskArrMath.h>
28 :
29 : #ifdef _OPENMP
30 : #include <omp.h>
31 : #endif
32 :
33 : using namespace casacore;
34 : namespace casa { //# NAMESPACE CASA - BEGIN
35 :
36 : namespace vi { //# NAMESPACE VI - BEGIN
37 :
38 :
39 :
40 :
41 : // -----------------------------------------------------------------------
42 : //
43 :
44 0 : CalibratingParameters::CalibratingParameters() :
45 : byCalLib_p(false),
46 : calLibRecord_p(Record()),
47 0 : corrFactor_p(1.0) // temporary, for initial testing (default is a non-trivial factor)
48 0 : {}
49 :
50 0 : CalibratingParameters::CalibratingParameters(const Record& calLibRecord) :
51 : byCalLib_p(false),
52 : calLibRecord_p(calLibRecord),
53 0 : corrFactor_p(1.0)
54 : {
55 :
56 0 : if (calLibRecord_p.isDefined("calfactor")) {
57 : //cout << "CalibratingParameters::ctor: found calfactor." << endl;
58 : // Detect calfactor in the specified Record
59 0 : corrFactor_p=calLibRecord_p.asFloat("calfactor");
60 0 : byCalLib_p = false; // signal not using a real callib
61 : }
62 0 : else if (calLibRecord_p.nfields()>0) {
63 : //cout << "CalibratingParameters::ctor: found non-trivial callib." << endl;
64 : // Apparently this will be a real callib
65 0 : byCalLib_p = true; // signal using a real callib
66 : }
67 : else
68 0 : throw(AipsError("Invalid use of callib Record"));
69 0 : validate();
70 0 : }
71 :
72 : // Construct using callib parser
73 0 : CalibratingParameters::CalibratingParameters(const String& callib) :
74 : byCalLib_p(true),
75 : calLibRecord_p(callibSetParams(callib)),
76 0 : corrFactor_p(1.0)
77 : {
78 0 : validate();
79 0 : }
80 :
81 0 : CalibratingParameters::CalibratingParameters(Float corrFactor) :
82 : byCalLib_p(false),
83 : calLibRecord_p(Record()),
84 0 : corrFactor_p(corrFactor) // temporary, for initial testing
85 : {
86 0 : validate();
87 0 : }
88 :
89 0 : CalibratingParameters::CalibratingParameters(const CalibratingParameters& other)
90 : {
91 0 : *this = other;
92 0 : }
93 :
94 0 : CalibratingParameters& CalibratingParameters::operator=(const CalibratingParameters& other)
95 : {
96 0 : if (this != &other) {
97 0 : byCalLib_p = other.byCalLib_p;
98 0 : calLibRecord_p = other.calLibRecord_p;
99 0 : corrFactor_p = other.corrFactor_p;
100 0 : validate();
101 : }
102 0 : return *this;
103 : }
104 :
105 0 : Bool CalibratingParameters::byCalLib() const
106 : {
107 0 : return byCalLib_p;
108 : }
109 :
110 0 : const Record& CalibratingParameters::getCalLibRecord() const
111 : {
112 0 : return calLibRecord_p;
113 : }
114 :
115 : // temporary, for initial testing
116 0 : Float CalibratingParameters::getCorrFactor() const
117 : {
118 0 : return corrFactor_p;
119 : }
120 :
121 :
122 0 : void CalibratingParameters::setCalLibRecord(const Record& calLibRecord)
123 : {
124 0 : calLibRecord_p = calLibRecord;
125 0 : }
126 :
127 : // temporary, for initial testing
128 0 : void CalibratingParameters::setCorrFactor(Float corrFactor)
129 : {
130 0 : corrFactor_p = corrFactor;
131 0 : }
132 :
133 0 : void CalibratingParameters::validate() const
134 : {
135 : // nothing meaningful to do yet
136 0 : }
137 :
138 :
139 :
140 : // -----------------------------------------------------------------------
141 0 : CalibratingVi2::CalibratingVi2( vi::ViImplementation2 * inputVii,
142 0 : const CalibratingParameters& calpar) :
143 : TransformingVi2 (inputVii),
144 : cb_p(0),
145 : ve_p(0),
146 0 : corrFactor_p(calpar.getCorrFactor()), // temporary
147 0 : visCalibrationOK_p(False)
148 : {
149 :
150 : // Initialize underlying ViImpl2
151 0 : getVii()->originChunks();
152 0 : getVii()->origin();
153 :
154 : // Make the internal VisBuffer2 for CalibratingVi2 clients
155 0 : setVisBuffer(createAttachedVisBuffer (VbRekeyable));
156 :
157 0 : }
158 :
159 : // -----------------------------------------------------------------------
160 0 : CalibratingVi2::CalibratingVi2( vi::ViImplementation2 * inputVii,
161 : const CalibratingParameters& calpar,
162 0 : String msname) :
163 : TransformingVi2 (inputVii),
164 : //cb_p(new OldCalibrater(msname)),
165 0 : cb_p(Calibrater::factory(msname)),
166 : ve_p(0),
167 : corrFactor_p(1.0),
168 0 : visCalibrationOK_p(False)
169 : {
170 :
171 0 : if (calpar.byCalLib()) {
172 : // Arrange calibration
173 0 : cb_p->validatecallib(calpar.getCalLibRecord());
174 0 : cb_p->setcallib2(calpar.getCalLibRecord(),&(inputVii->ms())); // Use underlying MS!
175 0 : cb_p->applystate();
176 : // Point to VisEquation
177 0 : ve_p = cb_p->ve();
178 : }
179 : else {
180 : // Simple mode using only the calfactor (good for tests)
181 0 : corrFactor_p=calpar.getCorrFactor();
182 : }
183 :
184 : // Initialize underlying ViImpl2
185 0 : getVii()->originChunks();
186 0 : getVii()->origin();
187 :
188 : // Make a VisBuffer for CalibratingVi2 clients
189 0 : setVisBuffer(createAttachedVisBuffer (VbRekeyable));
190 :
191 0 : }
192 : // -----------------------------------------------------------------------
193 0 : CalibratingVi2::CalibratingVi2( vi::ViImplementation2 * inputVii,
194 0 : VisEquation *ve) :
195 : TransformingVi2 (inputVii),
196 : cb_p(0),
197 : ve_p(ve),
198 : corrFactor_p(1.0),
199 0 : visCalibrationOK_p(False)
200 : {
201 :
202 : // Initialize underlying ViImpl2
203 0 : getVii()->originChunks();
204 0 : getVii()->origin();
205 :
206 : // Make a VisBuffer for CalibratingVi2 clients
207 0 : setVisBuffer(createAttachedVisBuffer (VbRekeyable));
208 :
209 0 : }
210 :
211 : // -----------------------------------------------------------------------
212 : //
213 : // -----------------------------------------------------------------------
214 0 : CalibratingVi2::~CalibratingVi2()
215 : {
216 : // cout << " ~CalVi2: " << this << endl;
217 : // ve_p is a borrowed pointer, so need not delete here
218 0 : ve_p=0;
219 :
220 : // Delete Calibrater object if present
221 0 : if (cb_p) delete cb_p;
222 0 : cb_p=0;
223 :
224 0 : }
225 :
226 :
227 :
228 : // -----------------------------------------------------------------------
229 : //
230 : // -----------------------------------------------------------------------
231 : void
232 0 : CalibratingVi2::origin()
233 : {
234 :
235 : // Drive underlying VII2
236 0 : getVii()->origin();
237 :
238 : // Keep my VB2 happily synchronized
239 : // (this comes from TransformingVi2)
240 0 : configureNewSubchunk();
241 :
242 : // Data/wts not yet corrected
243 0 : visCalibrationOK_p=False;
244 0 : }
245 :
246 : // -----------------------------------------------------------------------
247 : //
248 : // -----------------------------------------------------------------------
249 0 : void CalibratingVi2::next()
250 : {
251 :
252 : // Drive underlying VII2
253 0 : getVii()->next();
254 :
255 : // Keep my VB2 happily synchronized
256 : // (this comes from TransformingVi2)
257 0 : configureNewSubchunk();
258 :
259 : // Data/wts not yet corrected
260 0 : visCalibrationOK_p=False;
261 :
262 0 : }
263 :
264 : // -----------------------------------------------------------------------
265 : //
266 : // -----------------------------------------------------------------------
267 0 : void CalibratingVi2::flag(Cube<Bool>& flagC) const
268 : {
269 : // cout << "CVI2::flag(Cube)...";
270 :
271 : // Call for correction, which might set some flags
272 0 : calibrateCurrentVB();
273 :
274 : // copy result to caller's Cube<Bool>
275 0 : flagC.assign(getVii()->getVisBuffer()->flagCube());
276 :
277 0 : }
278 :
279 :
280 : // -----------------------------------------------------------------------
281 : //
282 : // -----------------------------------------------------------------------
283 0 : void CalibratingVi2::weight(Matrix<Float>& wt) const
284 : {
285 :
286 : // cout << "CVI2::weight...";
287 :
288 : // Call for correction
289 : // TBD: optimize w.r.t. calibrating only the weights?
290 0 : calibrateCurrentVB();
291 :
292 : // copy result to caller's Matrix<Float>
293 0 : wt.assign(getVii()->getVisBuffer()->weight());
294 :
295 0 : }
296 :
297 :
298 : // -----------------------------------------------------------------------
299 : //
300 : // -----------------------------------------------------------------------
301 0 : void CalibratingVi2::weightSpectrum(Cube<Float>& wtsp) const
302 : {
303 :
304 : // cout << "CVI2::weightSpectrum...";
305 :
306 0 : if (this->weightSpectrumExists()) {
307 :
308 : // Call for correction
309 : // TBD: optimize w.r.t. calibrating only the weights?
310 0 : calibrateCurrentVB();
311 :
312 : // copy result to caller's Matrix<Float>
313 0 : wtsp.assign(getVii()->getVisBuffer()->weightSpectrum());
314 :
315 : }
316 : else {
317 : // same as ordinary VisibilityIteratorImpl2
318 0 : wtsp.resize(0,0,0);
319 : }
320 0 : }
321 :
322 :
323 :
324 :
325 : // -----------------------------------------------------------------------
326 : //
327 : // -----------------------------------------------------------------------
328 0 : void CalibratingVi2::visibilityCorrected(Cube<Complex>& vis) const
329 : {
330 :
331 : // cout << "CVI2::visibilityCorrected...";
332 :
333 : // TBD:
334 : // o consider if underlying VisBuffer should be maintained const?
335 : // (and this layer's VB set and adjusted)
336 : // (will this break on-demand VB stuff?)
337 : // o make corresponding visibilityModel method... (solve context)
338 :
339 : // Call the actual correction method
340 : // (only does the actual work, if needed)
341 0 : calibrateCurrentVB();
342 :
343 : // copy result to caller's Cube<Complex>
344 0 : vis.assign(getVii()->getVisBuffer()->visCubeCorrected());
345 :
346 0 : }
347 :
348 :
349 :
350 : // -----------------------------------------------------------------------
351 : //
352 : // -----------------------------------------------------------------------
353 0 : Bool CalibratingVi2::existsColumn(VisBufferComponent2 id) const
354 : {
355 : Bool result;
356 0 : switch (id){
357 :
358 0 : case VisBufferComponent2::VisibilityCorrected:
359 : case VisBufferComponent2::VisibilityCubeCorrected:
360 : case VisBufferComponent2::VisibilityCubeModel:
361 :
362 0 : result = true;
363 0 : break;
364 :
365 0 : default:
366 0 : result = getVii()->existsColumn(id);
367 0 : break;
368 : }
369 :
370 0 : return result;
371 :
372 : }
373 :
374 : // -----------------------------------------------------------------------
375 : //
376 : // -----------------------------------------------------------------------
377 0 : void CalibratingVi2::calibrateCurrentVB() const
378 : {
379 : // This method must call NO ordinary VB2 accessors that require ViImpl
380 : // methods that are defined in this class, _even_implicitly_, because
381 : // the VB2 uses the VI2 that has its ViImpl overridden by these local
382 : // methods. This causes infinite loops!!!
383 : // One way to avoid this is to ensure that every method in this class
384 : // initializes the VB2 fields via getVii() methods.....
385 :
386 : // cout << " calibrateCurrentVB(): " << boolalpha << visCalibrationOK_p;
387 :
388 : // Do the correction, if not done yet
389 0 : if (!visCalibrationOK_p) {
390 :
391 : // We will operate on the underlying VB
392 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
393 :
394 : // sense if WEIGHT_SPECTRUM exists
395 0 : Bool doWtSp = getVii()->weightSpectrumExists();
396 :
397 : // Fill current flags, raw weights, and raw vis
398 0 : vb->flagCube();
399 0 : vb->resetWeightsUsingSigma(); // this is smart re spec weights or not
400 :
401 : // Initialize "corrected" data from "raw" data
402 0 : initCorrected(vb);
403 :
404 : // If the VisEquation is set, use it, otherwise use the corrFactor_p
405 0 : if (ve_p) {
406 : // Apply calibration via the VisEquation
407 0 : ve_p->correct2(*vb,False,doWtSp);
408 :
409 : // Set unchan'd weights, in case they are requested
410 0 : if (doWtSp)
411 0 : vb->setWeight(partialMedians(vb->weightSpectrum(),IPosition(1,1)));
412 :
413 : }
414 : else {
415 :
416 : // Use supplied corrFactor_p to make this corrected data look changed
417 0 : Cube<Complex> vCC(vb->visCubeCorrected());
418 0 : vCC*=corrFactor_p;
419 0 : vb->setVisCubeCorrected(vCC);
420 :
421 0 : if (doWtSp) {
422 : // Calibrate the WS
423 0 : Cube<Float> wS(vb->weightSpectrum()); // Was set above
424 0 : wS/=(corrFactor_p*corrFactor_p);
425 0 : vb->setWeightSpectrum(wS);
426 : // Set W via median on chan axis
427 0 : vb->setWeight(partialMedians(wS,IPosition(1,1)));
428 : }
429 : else {
430 : // Just calibrate the W
431 0 : Matrix<Float> w(vb->weight()); // Was set above
432 0 : w/=(corrFactor_p*corrFactor_p);
433 0 : vb->setWeight(w);
434 : }
435 : }
436 :
437 : // Signal that we have applied the correction, to avoid unnecessary redundancy
438 0 : visCalibrationOK_p=True;
439 :
440 : // cout << "-->" << visCalibrationOK_p;
441 :
442 : }
443 : // cout << endl;
444 0 : }
445 :
446 0 : void CalibratingVi2::initCorrected(casa::vi::VisBuffer2* vb) const {
447 :
448 0 : if (getVii()->existsColumn(VisBufferComponent2::VisibilityCubeFloat)) {
449 : // Convert FLOAT_DATA to Complex, and assign
450 0 : Cube<Float> f(vb->visCubeFloat());
451 0 : Cube<Complex> c;
452 0 : c.resize(f.shape());
453 0 : convertArray(c,f);
454 0 : vb->setVisCubeCorrected(c);
455 : }
456 : else
457 : // Just copy the (already-Complex DATA column)
458 0 : vb->setVisCubeCorrected(vb->visCube());
459 :
460 : // cout << "CalVi2::setCorrected()" << endl;
461 0 : }
462 :
463 0 : CalVi2LayerFactory::CalVi2LayerFactory(const CalibratingParameters& pars)
464 : : ViiLayerFactory(),
465 0 : calpars_(pars)
466 0 : {}
467 :
468 :
469 : // CalVi2-specific layer-creater
470 0 : ViImplementation2 * CalVi2LayerFactory::createInstance (ViImplementation2* vii0) const {
471 : // Make the CalibratingVi2, using supplied ViImplementation2, and return it
472 0 : ViImplementation2 *vii = new CalibratingVi2(vii0,calpars_);
473 0 : return vii;
474 : }
475 :
476 :
477 : /********************************/
478 :
479 :
480 : // -----------------------------------------------------------------------
481 0 : CalSolvingVi2::CalSolvingVi2( vi::ViImplementation2 * inputVii,
482 0 : const CalibratingParameters& calpar) :
483 : CalibratingVi2 (inputVii,calpar),
484 0 : corrDepFlags_(false)
485 : {
486 0 : corrFactor_p=sqrt(corrFactor_p);
487 0 : }
488 :
489 : // -----------------------------------------------------------------------
490 0 : CalSolvingVi2::CalSolvingVi2( vi::ViImplementation2 * inputVii,
491 : const CalibratingParameters& calpar,
492 0 : String msname) :
493 : CalibratingVi2 (inputVii,calpar,msname),
494 0 : corrDepFlags_(false)
495 : {
496 : // Nothing specialized to do here (except ctor parent, above)
497 0 : }
498 :
499 : // -----------------------------------------------------------------------
500 0 : CalSolvingVi2::CalSolvingVi2( vi::ViImplementation2 * inputVii,
501 : VisEquation *ve,
502 0 : const Bool& corrDepFlags) :
503 : CalibratingVi2 (inputVii,ve),
504 0 : corrDepFlags_(corrDepFlags)
505 : {
506 : // Nothing specialized to do here (except ctor parent, above)
507 0 : }
508 :
509 : //#define REPORTTIMING
510 :
511 : // -----------------------------------------------------------------------
512 : //
513 : // -----------------------------------------------------------------------
514 0 : CalSolvingVi2::~CalSolvingVi2()
515 : {
516 : // cout << " ~CalSolVi2: " << this << endl;
517 : // cout << "CSVi2::calibrateCurrentVB: Is WtSpec automatic?" << endl;
518 : // Nothing specialized to do here (except ctor parent, above)
519 :
520 : #ifdef REPORTTIMEING
521 : cout << "CalSolvingVi2 stats: " << endl;
522 : cout << " nVB_=" << nVB_
523 : << " nVB0_=" << nVB0_ << endl;
524 :
525 : cout << " Tcal_="
526 : << Tcalws_
527 : << "+" << Tcalfl_
528 : << "+" << Tcal2_
529 : << "=" << Tcalws_+Tcalfl_+Tcal2_
530 : << endl;
531 :
532 : cout << " Tio_=" << Tio_
533 : << endl;
534 : cout << " Samples="
535 : << " total=" << ntotal_
536 : << " flagged=" << nflagged_ << " (" << Double(nflagged_)/Double(ntotal_) << " of total)"
537 : << " skippable=" << nskipped_ << " (" << Double(nskipped_)/Double(nflagged_) << " of flagged)"
538 : << endl;
539 : cout << " spurious fraction=" << Double(nflagged_-nskipped_)/Double(ntotal_-nskipped_) << endl;
540 : cout << " good fraction=" << Double(ntotal_-nflagged_)/Double(ntotal_-nskipped_) << endl;
541 : #endif
542 :
543 :
544 0 : }
545 :
546 : // -----------------------------------------------------------------------
547 : //
548 : // -----------------------------------------------------------------------
549 : void
550 0 : CalSolvingVi2::originChunks(casacore::Bool forceRewind)
551 : {
552 :
553 : // Zero counters/timers
554 0 : ntotal_=nflagged_=nskipped_=nVB_=nVB0_=0;
555 0 : Tio_=Tcalws_=Tcalfl_=Tcal2_=0.0;
556 :
557 : // Call next layer
558 0 : getVii()->originChunks(forceRewind);
559 0 : }
560 :
561 : // -----------------------------------------------------------------------
562 : //
563 : // -----------------------------------------------------------------------
564 : void
565 0 : CalSolvingVi2::origin()
566 : {
567 :
568 : // Drive underlying VII2
569 0 : getVii()->origin();
570 :
571 : //ntotal_+=getVii()->getVisBuffer()->flagCube().nelements();
572 : //nflagged_+=ntrue(getVii()->getVisBuffer()->flagCube());
573 :
574 : // if (nfalse(getVii()->getVisBuffer()->flagCube())==0) {
575 : // cout << "Using (origin)..." << endl;
576 : // }
577 :
578 : /*
579 :
580 : // Step over completely flagged VB2s in the current chunk
581 : // NB: last one will be used in any case
582 : if (getVii()->more()) { // not already the last VB2 in this chunk
583 : while (getVii()->more()) {
584 : if (nfalse(getVii()->getVisBuffer()->flagCube())==0) {
585 : // This VB2 is entirely flagged, step to next one
586 : getVii()->next();
587 : cout << "skipping (origin)..." << endl;
588 : }
589 : else
590 : // This VB2 has some unflagged data, so use it
591 : break;
592 : }
593 : }
594 : */
595 : /*
596 : cout << "Samples (origin) ="
597 : << " total=" << ntotal_
598 : << " flagged=" << nflagged_
599 : << " skipped=" << nskipped_
600 : << endl;
601 : */
602 : // Keep my VB2 happily synchronized
603 : // (this comes from TransformingVi2)
604 0 : configureNewSubchunk();
605 :
606 : // Data/wts not yet corrected
607 0 : visCalibrationOK_p=False;
608 0 : }
609 :
610 : // -----------------------------------------------------------------------
611 : //
612 : // -----------------------------------------------------------------------
613 0 : void CalSolvingVi2::next()
614 : {
615 :
616 : // Drive underlying VII2
617 0 : getVii()->next();
618 :
619 : /*
620 :
621 : // Step over completely flagged VB2s in the current chunk
622 : // NB: last one will be used in any case
623 : if (getVii()->more()) { // not already the last VB2 in this chunk
624 : while (getVii()->more()) {
625 : if (nfalse(getVii()->getVisBuffer()->flagCube())==0) {
626 : // This VB2 is entirely flagged, step to next one
627 : ntotal_+=getVii()->getVisBuffer()->flagCube().nelements();
628 : nflagged_+=ntrue(getVii()->getVisBuffer()->flagCube());
629 : nskipped_+=ntrue(getVii()->getVisBuffer()->flagCube());
630 : getVii()->next();
631 : }
632 : else
633 : // This VB2 has some unflagged data, so use it
634 : break;
635 : }
636 : }
637 : */
638 :
639 : // Keep my VB2 happily synchronized
640 : // (this comes from TransformingVi2)
641 0 : configureNewSubchunk();
642 :
643 : // Data/wts not yet corrected
644 0 : visCalibrationOK_p=False;
645 :
646 : /*
647 : cout << "Samples (next) ="
648 : << " total=" << ntotal_
649 : << " flagged=" << nflagged_
650 : << " skipped=" << nskipped_
651 : << endl;
652 : */
653 :
654 0 : }
655 :
656 :
657 :
658 : // -----------------------------------------------------------------------
659 0 : void CalSolvingVi2::visibilityModel(Cube<Complex>& mod) const
660 : {
661 :
662 : // cout << "CSVI2::visibilityModel...";
663 :
664 : // Call the actual correction method
665 : // (only does the actual work, if needed)
666 0 : calibrateCurrentVB();
667 :
668 : // copy result to caller's Cube<Complex>
669 : // (this is no-op when mod has same data address as visCubeModel())
670 0 : mod.assign(getVii()->getVisBuffer()->visCubeModel());
671 :
672 0 : }
673 :
674 : // -----------------------------------------------------------------------
675 0 : void CalSolvingVi2::calibrateCurrentVB() const
676 : {
677 :
678 : //cout << " CalSolvingVi2::calibrateCurrentVB(): " << boolalpha << visCalibrationOK_p;
679 :
680 :
681 : // Do the correction, if not done yet
682 0 : if (!visCalibrationOK_p) {
683 :
684 : #ifdef _OPENMP
685 0 : Double time0=omp_get_wtime();
686 : #endif
687 :
688 : // Count the VBs that we are processing
689 0 : ++nVB_;
690 :
691 : // Get the underlying ViImpl2's VisBuffer, to munge it
692 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
693 :
694 : // sense if WEIGHT_SPECTRUM exists
695 : // (should be true in CalSolvingVi2!)
696 0 : Bool doWtSp = this->weightSpectrumExists();
697 :
698 : // Fill flags
699 0 : Cube<Bool> fl(vb->flagCube());
700 :
701 : // Fill model
702 : // NB: model-filling handled below, in ve_p (or on-demand)
703 0 : vb->visCubeModel();
704 :
705 : // Initialize the to-be-calibrated weights
706 : // (this fills WS, if available from below)
707 : // TBD: better semantics: vb->setCorrDataWtSpec(vb->dataWtSpec());
708 0 : vb->resetWeightsUsingSigma();
709 :
710 : // Initialize corrected data w/ data
711 0 : initCorrected(vb);
712 :
713 : // Counting
714 0 : Int64 nsamp=fl.nelements();
715 0 : if (nfalse(fl)==0) {
716 0 : nskipped_+=nsamp;
717 0 : ++nVB0_;
718 : }
719 0 : ntotal_+=nsamp;
720 0 : nflagged_+=ntrue(fl);
721 :
722 : #ifdef _OPENMP
723 0 : Double time1a=omp_get_wtime();
724 : #endif
725 :
726 : // Ensure we got weightSpectrum (fill it, if not)
727 0 : verifyWeightSpectrum(vb);
728 :
729 : #ifdef _OPENMP
730 0 : Double time1b=omp_get_wtime();
731 : #endif
732 :
733 : // Set old-style flags, FOR NOW
734 0 : if (!corrDepFlags_)
735 0 : corrIndepFlags(vb);
736 :
737 : #ifdef _OPENMP
738 0 : Double time2=omp_get_wtime();
739 : #endif
740 :
741 : // If the VisEquation is set, use it, otherwise use the corrFactor_p
742 0 : if (ve_p) {
743 :
744 : // Apply calibration via the VisEquation
745 0 : ve_p->collapse2(*vb); // ,False,doWtSp);
746 :
747 : // Set unchan'd weights, in case they are requested
748 : //if (doWtSp)
749 : // vb->setWeight(partialMedians(vb->weightSpectrum(),IPosition(1,1)));
750 :
751 : }
752 : else {
753 :
754 : // Use supplied corrFactor_p to make this corrected data look changed
755 : // In leiu of working VisEquation that is TBD
756 :
757 : // Correct data
758 0 : Cube<Complex> vCC(vb->visCubeCorrected());
759 0 : vCC*=corrFactor_p;
760 0 : vb->setVisCubeCorrected(vCC);
761 :
762 : // Corrupt the model
763 0 : Cube<Complex> vCM(vb->visCubeModel());
764 0 : vCM/=corrFactor_p;
765 0 : vb->setVisCubeModel(vCM);
766 :
767 0 : Cube<Float> modA2(square(amplitude(vCM)));
768 0 : Cube<Bool> modmask(modA2>0.0f);
769 :
770 0 : MaskedArray<Complex> vCCm=vCC(modmask);
771 0 : MaskedArray<Complex> vCMm=vCM(modmask);
772 0 : MaskedArray<Float> modA2m(modA2(modmask));
773 :
774 : // Divide corr data by mod data (where mod non-zero)
775 0 : vCCm=vCCm/vCMm;
776 0 : vCMm=Complex(1.0); // model now unity
777 : //vCMm=vCMm/vCMm; // less efficient?
778 :
779 0 : if (doWtSp) {
780 : // Calibrate the WS
781 0 : Cube<Float> wS(vb->weightSpectrum()); // Was set above
782 0 : wS/=(corrFactor_p*corrFactor_p);
783 :
784 : // adjust for model division
785 0 : MaskedArray<Float> wSm(wS(modmask));
786 0 : wSm=wSm*modA2m;
787 :
788 0 : vb->setWeightSpectrum(wS);
789 : // Set W via median on chan axis
790 0 : vb->setWeight(partialMedians(wS,IPosition(1,1)));
791 : }
792 : else {
793 : // Just calibrate the W
794 0 : Matrix<Float> w(vb->weight()); // Was set above
795 0 : w/=(corrFactor_p*corrFactor_p);
796 0 : vb->setWeight(w);
797 : }
798 : }
799 :
800 : // Signal that we have applied the correction, to avoid unnecessary redundancy
801 0 : visCalibrationOK_p=True;
802 :
803 : #ifdef _OPENMP
804 0 : Double time3=omp_get_wtime();
805 :
806 0 : Double dTio=time1a-time0;
807 0 : Double dTcalws=time1b-time1a;
808 0 : Double dTcalfl=time2-time1b;
809 0 : Double dTcal2=time3-time2;
810 0 : Tio_+=dTio;
811 0 : Tcalws_+=dTcalws;
812 0 : Tcalfl_+=dTcalfl;
813 0 : Tcal2_+=dTcal2;
814 : #endif
815 :
816 : /*
817 : cout << "nVB_=" << nVB_
818 : << " dTio = " << dTio
819 : << " dTcal1 = " << dTcal1
820 : << " dTcal2 = " << dTcal2
821 : << " Tio_ = " << Tio_
822 : << " Tcal = " << Tcal1
823 : << "+" << Tcal2_
824 : << "=" << Tcal1_+Tcal2_
825 : << endl;
826 : */
827 :
828 : }
829 :
830 : // cout << endl;
831 0 : }
832 :
833 0 : void CalSolvingVi2::corrIndepFlags(casa::vi::VisBuffer2* vb) const {
834 :
835 0 : if (vb->nCorrelations()==1)
836 : // Nothing to do if only one correlation
837 0 : return;
838 :
839 : // If more than one correlation, if any is flagged, all are (per row, chan)
840 0 : Cube<Bool> flc(vb->flagCube());
841 :
842 0 : VectorIterator<Bool> fvi(flc,0);
843 0 : Vector<Bool>& fviv(fvi.vector()); // stays sync'd!
844 0 : Int nCor=fviv.nelements();
845 0 : while (!fvi.pastEnd()) {
846 0 : for (Int icor=0;icor<nCor;++icor) {
847 0 : if (fviv(icor)) {
848 0 : fviv.set(true);
849 0 : continue; // jump out of icor loop if flag=T found
850 : }
851 : }
852 0 : fvi.next();
853 : }
854 :
855 0 : return;
856 : }
857 :
858 0 : void CalSolvingVi2::verifyWeightSpectrum(casa::vi::VisBuffer2* vb) const {
859 :
860 : // If we didn't get WS from below, populate it in the specified vb
861 0 : if (!getVii()->weightSpectrumExists()) {
862 0 : IPosition sh=vb->getShape();
863 0 : Cube<Float> wtsp(sh,0.0f);
864 :
865 : // Unchan'd weight as Cube w/ degenerate chan axis
866 0 : Cube<Float> wtsp0(vb->weight().reform(IPosition(3,sh(0),1,sh(2))));
867 :
868 0 : VectorIterator<Float> ivi(wtsp0,1);
869 0 : Vector<Float>& iviv(ivi.vector()); // stays sync'd!
870 0 : VectorIterator<Float> ovi(wtsp,1);
871 0 : Vector<Float>& oviv(ovi.vector()); // stays sync'd!
872 :
873 0 : while (!ivi.pastEnd()) {
874 0 : oviv.set(iviv(0));
875 0 : ivi.next();
876 0 : ovi.next();
877 : }
878 :
879 : // Set it in the vb
880 0 : vb->setWeightSpectrum(wtsp);
881 : }
882 :
883 : // In all cases set flagged cells' weights to zero
884 : // (This is ok in solving context.)
885 0 : Cube<Float> wtsp(vb->weightSpectrum());
886 0 : wtsp(vb->flagCube())=0.0;
887 :
888 0 : }
889 :
890 0 : CalSolvingVi2LayerFactory::CalSolvingVi2LayerFactory(const CalibratingParameters& pars)
891 0 : : CalVi2LayerFactory(pars)
892 0 : {}
893 :
894 :
895 : // CalSolvingVi2-specific layer-creater
896 0 : ViImplementation2 * CalSolvingVi2LayerFactory::createInstance (ViImplementation2* vii0) const {
897 : // Make the CalibratingVi2, using supplied ViImplementation2, and return it
898 0 : ViImplementation2 *vii = new CalSolvingVi2(vii0,calpars_);
899 0 : return vii;
900 : }
901 :
902 :
903 :
904 0 : CalSolvingVi2LayerFactoryByVE::CalSolvingVi2LayerFactoryByVE(VisEquation *ve,const casacore::Bool& corrDepFlags)
905 : : ViiLayerFactory(),
906 : ve_p(ve),
907 0 : corrDepFlags_(corrDepFlags)
908 : {
909 : // ve_p->state();
910 :
911 0 : }
912 :
913 : // CalSolvingVi2-specific layer-creater
914 0 : ViImplementation2 * CalSolvingVi2LayerFactoryByVE::createInstance (ViImplementation2* vii0) const {
915 : // Make the CalibratingVi2, using supplied ViImplementation2, and return it
916 0 : ViImplementation2 *vii = new CalSolvingVi2(vii0,ve_p,corrDepFlags_);
917 0 : return vii;
918 : }
919 :
920 :
921 :
922 :
923 :
924 :
925 :
926 :
927 : } //# NAMESPACE VI - END
928 : } //# NAMESPACE CASA - END
929 :
930 :
|
