Line data Source code
1 : #include <casacore/casa/Arrays/ArrayMath.h>
2 : #include <casacore/casa/Arrays/ArrayPartMath.h>
3 : #include <casacore/casa/BasicMath/Functors.h>
4 : #include <stdcasa/UtilJ.h>
5 : #include <msvis/MSVis/AveragingTvi2.h>
6 : #include <msvis/MSVis/AveragingVi2Factory.h>
7 : #include <msvis/MSVis/MsRows.h>
8 : #include <msvis/MSVis/VisBufferComponents2.h>
9 : #include <msvis/MSVis/VisBuffer2.h>
10 : #include <msvis/MSVis/VisBufferImpl2.h>
11 : #include <msvis/MSVis/Vbi2MsRow.h>
12 : #include <msvis/MSVis/VisibilityIterator2.h>
13 : #include <tuple>
14 : #include <set>
15 :
16 : using std::set;
17 :
18 : using namespace casacore;
19 : namespace casa {
20 :
21 : namespace vi {
22 :
23 : namespace avg {
24 :
25 : using casa::ms::MsRow;
26 :
27 :
28 : ///////////////////////////////////////////////////////////
29 : //
30 : // VbAvg: Averaging VisBuffer
31 : //
32 : /*
33 :
34 : Definition: A baseline sample (baseline for short) is a table row
35 : with a particular (antenna1, antenna2) pair at a particular time.
36 :
37 :
38 : Averaging does not cross chunk boundaries of the input VI so the
39 : input VI determines what values are averaged together. For example,
40 : if the input VI is allows data from multiple scans into the same chunk
41 : then some averaging across scans can occur; in this case the scan number
42 : of the result will be the scan number of the first baseline instance
43 : seen.
44 :
45 : Row-level value treatment
46 : =========================
47 :
48 : The average is approximately computed on a per baseline basis:
49 :
50 : averaged_baseline (antenna1, antenna2) =
51 : sumOverAveragingInterval (baseline (antenna1, antenna2)) /
52 : nBaselinesFoundInInterval
53 :
54 : How row-level values are computed
55 : ---------------------------------
56 :
57 : Time - Set to time of first baseline making up the average plus
58 : half of the averaging interval.
59 : Antenna1 - Copied from instance of baseline
60 : Antenna2 - Copied from instance of baseline
61 : Feed1 - Copied from instance of baseline
62 : Feed2 - Copied from instance of baseline
63 : Flag_Row - The flag row is the logical "and" of the flag rows
64 : averaged together.
65 : Data_Desc_Id - Copied from instance of baseline
66 : Processor_Id - Copied from instance of baseline
67 : Field_Id - Copied from instance of baseline
68 : Interval - Set to averaging interval
69 : Exposure - Minimum of the interval and the sum of the exposures
70 : from unflagged rows
71 : Time_Centroid - sum (timeCentroid[i] * exposure[i]) / sum (exposure[i])
72 : Scan_Number - Copied from instance of baseline
73 : Sigma - ???
74 : Array_Id - Copied from instance of baseline
75 : Observation_Id - Copied from instance of baseline
76 : State_Id - Copied from instance of baseline
77 : Uvw - Weighted average of the UVW values for the baseline
78 : Weight - ???
79 :
80 : Cube-level value treatment
81 : --------------------------
82 :
83 : For each baseline (i.e., antenna1, antenna2 pair) the average is
84 : computed for each correlation (co) and channel (ch) of the data cube.
85 :
86 : Data - sum (f(weightSpectrum (co, ch)) * data (co, ch)) /
87 : sum (f(weightSpectrum (co, ch)))
88 : f :== optional function applied to the weights; default is unity function.
89 : Corrected_Data - Same was for Data however, VI setup can disable producing
90 : averaged Corrected_Data
91 : Model_Data - Same was for Data however, VI setup can disable producing
92 : averaged Model_Data
93 : Weight_Spectrum - sum (weightSpectrum (co, ch))
94 : Flag - Each averaged flag (correlation, channel) is the logical "and"
95 : of the flag (correlation, channel) making up the average.
96 :
97 : */
98 :
99 :
100 : class BaselineIndex {
101 :
102 : public:
103 :
104 : // make noncopyable...
105 : BaselineIndex( const BaselineIndex& ) = delete;
106 : BaselineIndex& operator=( const BaselineIndex& ) = delete;
107 :
108 : BaselineIndex ();
109 : ~BaselineIndex ();
110 :
111 : Int operator () (Int antenna1, Int antenna2, Int spectralWindow);
112 : void configure (Int nAntennas, Int nSpw, const VisBuffer2 * vb);
113 :
114 :
115 : private:
116 :
117 : enum {Empty = -1};
118 :
119 : class SpwIndex : public Matrix<Int>{
120 :
121 : public:
122 :
123 0 : SpwIndex (Int n) : Matrix<Int> (n, n, Empty) {}
124 :
125 : Int
126 0 : getBaselineNumber (Int antenna1, Int antenna2, Int & nBaselines)
127 : {
128 0 : Int i = (* this) (antenna1, antenna2);
129 :
130 0 : if (i == Empty){
131 :
132 0 : i = nBaselines ++;
133 0 : (* this) (antenna1, antenna2) = i;
134 : }
135 :
136 0 : return i;
137 : }
138 :
139 : private:
140 :
141 : };
142 :
143 : typedef vector<SpwIndex *> Index;
144 :
145 : SpwIndex * addSpwIndex (Int spw);
146 : Matrix<Int> * createMatrix ();
147 : void destroy();
148 : SpwIndex * getSpwIndex (Int spw);
149 :
150 : Index index_p;
151 : Int nAntennas_p;
152 : Int nBaselines_p;
153 : Int nSpw_p;
154 :
155 : };
156 :
157 0 : BaselineIndex::BaselineIndex ()
158 : : nAntennas_p (0),
159 : nBaselines_p (0),
160 0 : nSpw_p (0)
161 0 : {}
162 :
163 0 : BaselineIndex::~BaselineIndex ()
164 : {
165 0 : destroy();
166 0 : }
167 :
168 : Int
169 0 : BaselineIndex::operator () (Int antenna1, Int antenna2, Int spectralWindow)
170 : {
171 0 : SpwIndex * spwIndex = getSpwIndex (spectralWindow);
172 0 : if (spwIndex == 0){
173 0 : addSpwIndex (spectralWindow);
174 : }
175 :
176 0 : Int i = spwIndex->getBaselineNumber (antenna1, antenna2, nBaselines_p);
177 :
178 0 : return i;
179 : }
180 :
181 :
182 :
183 : BaselineIndex::SpwIndex *
184 0 : BaselineIndex::addSpwIndex (Int i)
185 : {
186 : // Delete an existing SpwIndex so that we start fresh
187 :
188 0 : delete index_p [i];
189 :
190 : // Create a new SpwIndex and insert it into the main index.
191 :
192 0 : index_p [i] = new SpwIndex (nAntennas_p);
193 :
194 0 : return index_p [i];
195 : }
196 :
197 : void
198 0 : BaselineIndex::configure (Int nAntennas, Int nSpw, const VisBuffer2 * vb)
199 : {
200 : // Capture the shape parameters
201 :
202 0 : nAntennas_p = nAntennas;
203 0 : nSpw_p = nSpw;
204 0 : nBaselines_p = 0;
205 :
206 : // Get rid of the existing index
207 :
208 0 : destroy ();
209 0 : index_p = Index (nSpw_p, (SpwIndex *) 0);
210 :
211 : // Fill out the index based on the contents of the first VB.
212 : // Usually this will determine the pattern for all of the VBs to be
213 : // averaged together so that is the ordering the index should
214 : // capture.
215 :
216 0 : for (rownr_t i = 0; i < vb->nRows(); i++){
217 :
218 : // Eagerly flesh out the Spw's index
219 :
220 0 : Int spw = vb->spectralWindows() (i);
221 0 : Int antenna1 = vb->antenna1()(i);
222 0 : Int antenna2 = vb->antenna2()(i);
223 :
224 0 : (* this) (antenna1, antenna2, spw);
225 : }
226 0 : }
227 :
228 : void
229 0 : BaselineIndex::destroy ()
230 : {
231 : // Delete all the dynamically allocated spectral window indices.
232 : // The vector destructor will take care of index_p itself.
233 :
234 0 : for (Index::iterator i = index_p.begin();
235 0 : i != index_p.end();
236 0 : i++){
237 :
238 0 : delete (* i);
239 : }
240 0 : }
241 :
242 : BaselineIndex::SpwIndex *
243 0 : BaselineIndex::getSpwIndex (Int spw)
244 : {
245 0 : AssertCc (spw < (int) index_p.size());
246 :
247 0 : SpwIndex * spwIndex = index_p [spw];
248 :
249 0 : if (spwIndex == 0){
250 0 : spwIndex = addSpwIndex (spw);
251 : }
252 :
253 0 : return spwIndex;
254 : }
255 :
256 : template <typename T>
257 : class PrefilledMatrix {
258 :
259 : public:
260 :
261 0 : PrefilledMatrix () : matrix_p (0, 0, 0), nChannels_p (0), nCorrelations_p (0) {}
262 :
263 : const Matrix<T> &
264 0 : getMatrix (Int nCorrelations, Int nChannels, const T & value)
265 : {
266 0 : if (nCorrelations != nCorrelations_p || nChannels != nChannels_p ||
267 0 : value != value_p){
268 :
269 0 : nCorrelations_p = nCorrelations;
270 0 : nChannels_p = nChannels;
271 0 : value_p = value;
272 :
273 0 : matrix_p.assign (Matrix<T> (nCorrelations_p, nChannels_p, value_p));
274 : }
275 :
276 0 : return matrix_p;
277 : }
278 :
279 : private:
280 :
281 : Matrix<T> matrix_p;
282 : Int nChannels_p;
283 : Int nCorrelations_p;
284 : T value_p;
285 :
286 : };
287 :
288 : template <typename T>
289 : class CachedPlaneAvg : public ms::CachedArrayBase {
290 :
291 : public:
292 :
293 : typedef const Cube<T> & (casa::vi::avg::VbAvg::* Accessor) () const;
294 :
295 0 : CachedPlaneAvg (Accessor accessor) : accessor_p (accessor) {}
296 :
297 : Matrix<T> &
298 0 : getCachedPlane (casa::vi::avg::VbAvg * vb, Int row)
299 : {
300 0 : if (! isCached()){
301 :
302 : //cache_p.reference ((vb ->* accessor_p)().xyPlane (row)); // replace with something more efficient
303 0 : referenceMatrix (cache_p, (vb ->* accessor_p)(), row);
304 0 : setCached ();
305 : }
306 :
307 0 : return cache_p;
308 : }
309 :
310 : private:
311 :
312 : static void
313 0 : referenceMatrix (Matrix<T> & cache, const Cube<T> & src, Int row)
314 : {
315 0 : IPosition shape = src.shape ();
316 0 : shape.resize (2);
317 :
318 : // This is a bit sleazy but it seems to be helpful to performance.
319 : // Assumes contiguously stored cube.
320 :
321 0 : T * storage = const_cast <T *> (& src (IPosition (3, 0, 0, row)));
322 :
323 0 : cache.takeStorage (shape, storage, casacore::SHARE);
324 0 : }
325 :
326 : Accessor accessor_p;
327 : Matrix<T> cache_p;
328 : };
329 :
330 :
331 : class VbAvg;
332 :
333 : /**
334 : * Holds multiple rows in a buffer. changeRow() alternates between rows. The details
335 : * of how the rows are handled need to be traced to its parent class, Vbi2MsRow, and its
336 : * parent MsRow.
337 : */
338 : class MsRowAvg : public ms::Vbi2MsRow {
339 :
340 : public:
341 :
342 : MsRowAvg (rownr_t row, const VbAvg * vb);
343 :
344 : // Constructor for read/write access
345 :
346 : MsRowAvg (rownr_t row, VbAvg * vb);
347 :
348 0 : virtual ~MsRowAvg () {}
349 :
350 : Bool baselinePresent () const;
351 : /// For how many of the rows reachable by changeRow() does baselinePresent(() hold?
352 : /// That's equivalent to asking how many baselines are being
353 : Int nBaselinesPresent () const;
354 :
355 : Vector<Bool> correlationFlagsMutable ();
356 : const Matrix<Int> & counts () const;
357 : Int countsBaseline () const;
358 0 : Matrix<Bool> & flagsMutable () { return Vbi2MsRow::flagsMutable();}
359 : Double intervalLast () const;
360 : Double timeFirst () const;
361 : Double timeLast () const;
362 : Vector<Double> uvwFirst ();
363 :
364 : void setBaselinePresent (Bool isPresent);
365 : void setCounts (const Matrix<Int> &);
366 : void setCountsBaseline (Int);
367 : void setIntervalLast (Double);
368 : void setTimeFirst (Double);
369 : void setTimeLast (Double);
370 :
371 : Double getNormalizationFactor();
372 : void setNormalizationFactor(Double normalizationFactor);
373 : void accumulateNormalizationFactor(Double normalizationFactor);
374 :
375 : /**
376 : * For bookkeeping. Input row indices that have been added so far in the current row
377 : * (set with changeRow). This is a list of input rows attached to this averaged row,
378 : * needed to build the map of input->output row indices when this row is transferred to
379 : * the output/averaged buffer.
380 : */
381 : std::vector<Int> inputRowIdxs() { return inputRowIdxs_p[row()]; }
382 : /**
383 : * Adds into the 'inputRowIdxs' list the index of an input row from the buffer being
384 : * averaged.
385 : * @param idx the index of the input row in the input buffer
386 : */
387 0 : void addInputRowIdx(Int idx) {
388 0 : inputRowIdxs_p[row()].push_back(idx);
389 0 : }
390 : /**
391 : * Clear the list of input rows attached to this row. This is used once the row is
392 : * transferred to the output/averaged buffer (typically after every average interval).
393 : */
394 : void clearRowIdxs() { inputRowIdxs_p[row()].clear(); }
395 :
396 : private:
397 :
398 : void configureCountsCache ();
399 :
400 : mutable CachedPlaneAvg<Int> countsCache_p;
401 : Vector<Double> normalizationFactor_p;
402 : VbAvg * vbAvg_p; // [use]
403 : // map: input buffer row index -> output buffer row index
404 : std::map<Int, std::vector<Int>> inputRowIdxs_p;
405 : };
406 :
407 : /**
408 : * It looks like the intended usage of this VbAvg class (from AveragingTvi2::
409 : * produceSubchunk()) is as follows:
410 : *
411 : * // Use a "VbAvg vbAvg":
412 : * VisBufferImpl2 * vbToFill = // get/create output (averaged) buffer
413 : * vbToFill->setFillable(true);
414 : * vbAvg.setBufferToFill(vbToFill);
415 : * // we have the input buffer (to be averaged) in "VisibilityIteratorImpl2 vii;
416 : * while (vii->more()) {
417 : * ...
418 : * vbAvg.accumulate(vb, subhunk);
419 : * }
420 : * vbAvg.finalizeAverages();
421 : * vbAvg.finalizeBufferFillnig();
422 : */
423 : class VbAvg : public VisBufferImpl2 {
424 :
425 : public:
426 :
427 : friend class MsRowAvg;
428 :
429 : VbAvg (const AveragingParameters & averagingParameters, ViImplementation2 * vi);
430 :
431 : void accumulate (const VisBuffer2 * vb, const Subchunk & subchunk);
432 : const Cube<Int> & counts () const;
433 : Bool empty () const;
434 : void finalizeBufferFilling ();
435 : void finalizeAverages ();
436 : MsRowAvg * getRow (Int row) const;
437 : MsRowAvg * getRowMutable (Int row);
438 : Bool isComplete () const;
439 : Bool isUsingUvwDistance () const;
440 : void markEmpty ();
441 : int nSpectralWindowsInBuffer () const;
442 : void setBufferToFill (VisBufferImpl2 *);
443 : void startChunk (ViImplementation2 *);
444 : Int getBaselineIndex (Int antenna1, Int antenna2, Int spw) const;
445 : // Vector with row idx in the averaged/ooutput buffers
446 0 : const std::vector<size_t> & row2AvgRow() const { return row2AvgRow_p; };
447 :
448 : protected:
449 :
450 : class Doing {
451 : public:
452 :
453 0 : Doing ()
454 0 : : correctedData_p (false),
455 : modelData_p (false),
456 : observedData_p (false),
457 : floatData_p(false),
458 : onlymetadata_p(true),
459 : weightSpectrumIn_p (false),
460 : sigmaSpectrumIn_p (false),
461 : weightSpectrumOut_p (false),
462 0 : sigmaSpectrumOut_p (false)
463 0 : {}
464 :
465 : Bool correctedData_p;
466 : Bool modelData_p;
467 : Bool observedData_p;
468 : Bool floatData_p;
469 : Bool onlymetadata_p;
470 : Bool weightSpectrumIn_p;
471 : Bool sigmaSpectrumIn_p;
472 : Bool weightSpectrumOut_p;
473 : Bool sigmaSpectrumOut_p;
474 : };
475 :
476 : class AccumulationParameters {
477 :
478 : public:
479 :
480 0 : AccumulationParameters (MsRow * rowInput, MsRowAvg * rowAveraged, const Doing & doing)
481 0 : : correctedIn_p (doing.correctedData_p ? rowInput->corrected().data() : 0),
482 0 : correctedOut_p (doing.correctedData_p ? rowAveraged->correctedMutable ().data(): 0),
483 0 : flagCubeIn_p (rowInput->flags().data()),
484 0 : flagCubeOut_p (rowAveraged->flagsMutable().data()),
485 0 : floatDataIn_p (doing.floatData_p ? rowInput->singleDishData().data() : 0),
486 0 : floatDataOut_p (doing.floatData_p ? rowAveraged->singleDishDataMutable().data() : 0),
487 0 : modelIn_p (doing.modelData_p ? rowInput->model().data(): 0),
488 0 : modelOut_p (doing.modelData_p ? rowAveraged->modelMutable().data() : 0),
489 0 : observedIn_p (doing.observedData_p ? rowInput->observed().data() : 0),
490 0 : observedOut_p (doing.observedData_p ? rowAveraged->observedMutable().data() : 0),
491 0 : sigmaIn_p (& rowInput->sigma()),
492 0 : sigmaOut_p (& rowAveraged->sigmaMutable()),
493 0 : sigmaSpectrumIn_p (doing.sigmaSpectrumIn_p ? rowInput->sigmaSpectrum().data() : 0),
494 0 : sigmaSpectrumOut_p (doing.sigmaSpectrumOut_p ? rowAveraged->sigmaSpectrumMutable().data() : 0),
495 0 : weightIn_p (& rowInput->weight()),
496 0 : weightOut_p (& rowAveraged->weightMutable()),
497 0 : weightSpectrumIn_p (doing.weightSpectrumIn_p ? rowInput->weightSpectrum().data() : 0),
498 0 : weightSpectrumOut_p (doing.weightSpectrumOut_p ? rowAveraged->weightSpectrumMutable().data() : 0)
499 0 : {}
500 :
501 0 : void incrementCubePointers()
502 : {
503 : // For improved performance this class is designed to sweep the cube data elements in this row
504 : // in the order (correlation0, channel0), (correlation1, channel1), etc.
505 :
506 0 : correctedIn_p && correctedIn_p ++;
507 0 : correctedOut_p && correctedOut_p ++;
508 0 : flagCubeIn_p && flagCubeIn_p ++;
509 0 : flagCubeOut_p && flagCubeOut_p ++;
510 0 : floatDataIn_p && floatDataIn_p ++;
511 0 : floatDataOut_p && floatDataOut_p ++;
512 0 : modelIn_p && modelIn_p ++;
513 0 : modelOut_p && modelOut_p ++;
514 0 : observedIn_p && observedIn_p ++;
515 0 : observedOut_p && observedOut_p ++;
516 0 : sigmaSpectrumIn_p && sigmaSpectrumIn_p ++;
517 0 : sigmaSpectrumOut_p && sigmaSpectrumOut_p ++;
518 0 : weightSpectrumIn_p && weightSpectrumIn_p ++;
519 0 : weightSpectrumOut_p && weightSpectrumOut_p ++;
520 0 : }
521 :
522 : inline const Complex *
523 0 : correctedIn ()
524 : {
525 0 : assert (correctedIn_p != 0);
526 0 : return correctedIn_p;
527 : }
528 :
529 : inline Complex *
530 0 : correctedOut ()
531 : {
532 0 : assert (correctedOut_p != 0);
533 0 : return correctedOut_p;
534 : }
535 :
536 : inline const Float *
537 0 : floatDataIn ()
538 : {
539 0 : return floatDataIn_p;
540 : }
541 :
542 : inline Float *
543 0 : floatDataOut ()
544 : {
545 0 : return floatDataOut_p;
546 : }
547 :
548 : inline const Complex *
549 0 : modelIn ()
550 : {
551 0 : assert (modelIn_p != 0);
552 0 : return modelIn_p;
553 : }
554 :
555 : inline Complex *
556 0 : modelOut ()
557 : {
558 0 : assert (modelOut_p != 0);
559 0 : return modelOut_p;
560 : }
561 :
562 : inline const Complex *
563 0 : observedIn ()
564 : {
565 0 : assert (observedIn_p != 0);
566 0 : return observedIn_p;
567 : }
568 :
569 : inline Complex *
570 0 : observedOut ()
571 : {
572 0 : assert (observedOut_p != 0);
573 0 : return observedOut_p;
574 : }
575 :
576 : inline const Bool *
577 : flagCubeIn ()
578 : {
579 : assert (flagCubeIn_p != 0);
580 : return flagCubeIn_p;
581 : }
582 :
583 : inline Bool *
584 : flagCubeOut ()
585 : {
586 : assert (flagCubeOut_p != 0);
587 : return flagCubeOut_p;
588 : }
589 :
590 : inline const Float *
591 0 : sigmaSpectrumIn ()
592 : {
593 0 : return sigmaSpectrumIn_p;
594 : }
595 :
596 : inline Float *
597 0 : sigmaSpectrumOut ()
598 : {
599 0 : assert (sigmaSpectrumOut_p != 0);
600 0 : return sigmaSpectrumOut_p;
601 : }
602 :
603 : inline const Float *
604 0 : weightSpectrumIn ()
605 : {
606 0 : assert (weightSpectrumIn_p != 0);
607 0 : return weightSpectrumIn_p;
608 : }
609 :
610 : inline Float *
611 0 : weightSpectrumOut ()
612 : {
613 0 : assert (weightSpectrumOut_p != 0);
614 0 : return weightSpectrumOut_p;
615 : }
616 :
617 : inline const Vector<Float> &
618 : weightIn ()
619 : {
620 : assert (weightIn_p != 0);
621 : return *weightIn_p;
622 : }
623 :
624 : inline Vector<Float> &
625 : weightOut ()
626 : {
627 : assert (weightOut_p != 0);
628 : return *weightOut_p;
629 : }
630 :
631 : inline const Vector<Float> &
632 : sigmaIn ()
633 : {
634 : assert (sigmaIn_p != 0);
635 : return *sigmaIn_p;
636 : }
637 :
638 : inline Vector<Float> &
639 : sigmaOut ()
640 : {
641 : assert (sigmaOut_p != 0);
642 : return *sigmaOut_p;
643 : }
644 :
645 :
646 :
647 : private:
648 :
649 : const Complex * correctedIn_p;
650 : Complex * correctedOut_p;
651 : const Bool * flagCubeIn_p;
652 : Bool * flagCubeOut_p;
653 : const Float * floatDataIn_p;
654 : Float * floatDataOut_p;
655 : const Complex * modelIn_p;
656 : Complex * modelOut_p;
657 : const Complex * observedIn_p;
658 : Complex * observedOut_p;
659 : const Vector<Float> * sigmaIn_p;
660 : Vector<Float> * sigmaOut_p;
661 : const Float * sigmaSpectrumIn_p;
662 : Float * sigmaSpectrumOut_p;
663 : const Vector<Float> * weightIn_p;
664 : Vector<Float> * weightOut_p;
665 : const Float * weightSpectrumIn_p;
666 : Float * weightSpectrumOut_p;
667 :
668 :
669 : };
670 :
671 : std::pair<Bool, Vector<Double> > accumulateCubeData (MsRow * rowInput, MsRowAvg * rowAveraged);
672 : void accumulateElementForCubes (AccumulationParameters & accumulationParameters,
673 : Bool zeroAccumulation);
674 : template<typename T>
675 : void
676 : accumulateElementForCube (const T * unweightedValue,
677 : Float weight,
678 : Bool zeroAccumulation,
679 : T * accumulator);
680 :
681 : template <typename T>
682 : T accumulateRowDatum (const T & averagedValue, const T & inputValue,
683 : Bool resetAverage);
684 :
685 :
686 :
687 : void accumulateExposure (const VisBuffer2 *);
688 : /*
689 : * Called once per row in the input buffer
690 : * @param rowInput row from the input buffer being averaged
691 : * @param rowAveraged changing "accumulator" row for the output buffer
692 : * @param subchunk - hard to understand
693 : * @param iidx index of the input row in the input buffer
694 : */
695 : void accumulateOneRow (MsRow * rowInput, MsRowAvg * rowAveraged,
696 : const Subchunk & subchunk, Int iidx);
697 : void accumulateRowData (MsRow * rowInput, MsRowAvg * rowAveraged, Double adjustedWeight,
698 : Bool rowFlagged);
699 : void accumulateTimeCentroid (const VisBuffer2 * input);
700 : void captureIterationInfo (VisBufferImpl2 * dstVb, const VisBuffer2 * srcVb,
701 : const Subchunk & subchunk);
702 : void copyBaseline (Int sourceBaseline, Int targetBaseline);
703 : template<typename T>
704 : void copyBaselineScalar (Int sourceBaseline, Int targetBaseline,
705 : Vector<T> & columnVector);
706 : template<typename T>
707 : void copyCubePlaneIf (Bool condition, Int sourceBaseline,
708 : Int targetBaseline, Cube<T> & cube);
709 : void copyIdValues (MsRow * rowInput, MsRowAvg * rowAveraged);
710 : void copyIdValue (Int inputId, Int & averagedId);
711 : void finalizeBaseline (MsRowAvg *);
712 : void finalizeBaselineIfNeeded (MsRow * rowInput, MsRowAvg * rowAveraged,
713 : const Subchunk & subchunk);
714 : void finalizeCubeData (MsRowAvg *);
715 : void finalizeRowData (MsRowAvg *);
716 : AccumulationParameters * getAccumulationParameters (MsRow * rowInput,
717 : MsRowAvg * rowAveraged);
718 : Int getBaselineIndex (const MsRow *) const;
719 : void initializeBaseline (MsRow * rowInput, MsRowAvg * rowAveraged,
720 : const Subchunk & subchunk);
721 : Int nBaselines () const;
722 : void prepareIds (const VisBuffer2 * vb);
723 : void removeMissingBaselines ();
724 :
725 : private:
726 :
727 : void setupVbAvg (const VisBuffer2 *);
728 : void setupArrays (Int nCorrelations, Int nChannels, Int nBaselines);
729 : void setupBaselineIndices (Int nAntennas, const VisBuffer2 * vb);
730 :
731 : /// A baseline being averaged into a MSRowAvg is put into the output/averaged buffer and
732 : /// set as not present
733 : void transferBaseline (MsRowAvg *);
734 :
735 : template <typename T>
736 : static T
737 0 : distanceSquared (const Vector<T> & p1, const Vector<T> & p2)
738 : {
739 0 : assert (p1.size() == 3 && p2.size() == 3);
740 :
741 0 : T distanceSquared = 0;
742 :
743 0 : for (Int i = 0; i < 3; i++){
744 0 : T delta = p1[i] - p2[i];
745 0 : distanceSquared += delta * delta;
746 : }
747 :
748 0 : return distanceSquared;
749 : }
750 :
751 : Double averagingInterval_p;
752 : AveragingOptions averagingOptions_p;
753 : const ViImplementation2 * averagingVii_p;
754 : mutable BaselineIndex baselineIndex_p; // map of antenna1,antenna2 to row number in this VB.
755 : Vector<Bool> baselinePresent_p; // indicates whether baseline has any data
756 : Vector<Double> normalizationFactor_p; // indicates whether baseline has any data
757 : VisBufferImpl2 * bufferToFill_p;
758 : Bool complete_p; // average is complete
759 : Matrix<Bool> correlationFlags_p; // used for weight accumulation
760 : Cube<Int> counts_p; // number of items summed together for each correlation, channel, baseline item.
761 : Vector<Int> countsBaseline_p; // number of items summed together for each baseline.
762 : Doing doing_p;
763 : Bool empty_p; // true when buffer hasn't seen any data
764 : Vector<Double> intervalLast_p;
765 : Double maxTimeDistance_p;
766 : Double maxUvwDistanceSquared_p;
767 : Bool needIterationInfo_p;
768 : VisBufferComponents2 optionalComponentsToCopy_p;
769 : Int rowIdGenerator_p;
770 : Double sampleInterval_p;
771 : Double startTime_p; // time of the first sample in average
772 : Vector<Double> timeFirst_p;
773 : Vector<Double> timeLast_p;
774 : mutable PrefilledMatrix<Bool> trueBool_p;
775 : Matrix<Double> uvwFirst_p;
776 : Bool usingUvwDistance_p;
777 : mutable PrefilledMatrix<Int> zeroInt_p;
778 : mutable PrefilledMatrix<Float> zeroFloat_p;
779 :
780 : std::vector<size_t> row2AvgRow_p;
781 :
782 : LogIO logger_p;
783 : };
784 :
785 : ///////////////////////////////////////////////////////////
786 : //
787 : // Set of Averaging VisBuffers, one per active DD ID
788 :
789 :
790 : //class VbSet {
791 : //
792 : //public:
793 : //
794 : // VbSet (const AveragingParameters & averagingParameters);
795 : // ~VbSet ();
796 : //
797 : // void accumulate (const VisBuffer2 *, const Subchunk & subchunk);
798 : // void finalizeBufferFilling ();
799 : // void setBufferToFill (VisBuffer2 *);
800 : //
801 : //
802 : // VbAvg * add (Int ddId);
803 : // Bool anyAveragesReady(Int ddid = -1) const;
804 : // Bool contains (Int ddId) const;
805 : //// void finalizeAverage (Int ddId);
806 : // void finalizeAverages ();
807 : // void finalizeRowIfNeeded (ms::MsRow * rowInput, avg::MsRowAvg * rowAveraged, const Subchunk & subchunk);
808 : // void flush (Bool okIfNonempty = false, ViImplementation2 * vi = 0); // delete all averagers
809 : // Int getFirstReadyDdid () const;
810 : // void transferAverage (Int ddId, VisBuffer2 * vb);
811 : // Bool vbPastAveragingInterval (const VisBuffer2 * vb) const;
812 : // void zero ();
813 : //
814 : //protected:
815 : //
816 : // void seeIfCubeColumnsExist (ViImplementation2 * vi);
817 : //
818 : //private:
819 : //
820 : // typedef map<Int, VbAvg *> Averagers;
821 : //
822 : // const Double averagingInterval_p;
823 : // AveragingOptions averagingOptions_p;
824 : // const AveragingParameters averagingParameters_p;
825 : // Bool doingCorrectedData_p;
826 : // Bool doingModelData_p;
827 : // Bool doingObservedData_p;
828 : // Bool doingWeightSpectrum_p;
829 : // Bool doingsigmaSpectrum_p;
830 : // Averagers vbAveragers_p;
831 : //};
832 :
833 0 : MsRowAvg::MsRowAvg (rownr_t row, const VbAvg * vb)
834 : : Vbi2MsRow (row, vb),
835 : countsCache_p (& VbAvg::counts),
836 : normalizationFactor_p(0.0),
837 0 : vbAvg_p (const_cast<VbAvg *> (vb))
838 : {
839 0 : configureCountsCache();
840 0 : }
841 :
842 :
843 : // Constructor for read/write access
844 :
845 0 : MsRowAvg::MsRowAvg (rownr_t row, VbAvg * vb)
846 : : Vbi2MsRow (row, vb),
847 : countsCache_p (& VbAvg::counts),
848 : normalizationFactor_p(0.0),
849 0 : vbAvg_p (vb)
850 : {
851 0 : configureCountsCache();
852 0 : }
853 :
854 : Bool
855 0 : MsRowAvg::baselinePresent () const
856 : {
857 0 : return vbAvg_p->baselinePresent_p (row ());
858 : }
859 :
860 : Int
861 0 : MsRowAvg::nBaselinesPresent () const
862 : {
863 0 : return std::count(vbAvg_p->baselinePresent_p.begin(), vbAvg_p->baselinePresent_p.end(),
864 0 : true);
865 : }
866 :
867 : void
868 0 : MsRowAvg::configureCountsCache ()
869 : {
870 0 : addToCachedArrays (countsCache_p);
871 0 : }
872 :
873 : const Matrix<Int> &
874 0 : MsRowAvg::counts () const
875 : {
876 0 : return countsCache_p.getCachedPlane (dynamic_cast<VbAvg *> (getVbi()), row());
877 : }
878 :
879 : Vector<Bool>
880 0 : MsRowAvg::correlationFlagsMutable ()
881 : {
882 0 : return vbAvg_p->correlationFlags_p.column (row());
883 : }
884 :
885 : Int
886 0 : MsRowAvg::countsBaseline () const
887 : {
888 0 : return vbAvg_p->countsBaseline_p (row ());
889 : }
890 :
891 : void
892 0 : MsRowAvg::setBaselinePresent (Bool value)
893 : {
894 0 : vbAvg_p->baselinePresent_p (row ()) = value;
895 0 : }
896 :
897 :
898 : void
899 0 : MsRowAvg::setCountsBaseline (Int value)
900 : {
901 0 : vbAvg_p->countsBaseline_p (row ()) = value;
902 0 : }
903 :
904 : Double
905 0 : MsRowAvg::intervalLast () const
906 : {
907 0 : return vbAvg_p->intervalLast_p (row ());
908 : }
909 :
910 :
911 : Double
912 0 : MsRowAvg::timeFirst () const
913 : {
914 0 : return vbAvg_p->timeFirst_p (row ());
915 : }
916 :
917 : Double
918 0 : MsRowAvg::timeLast () const
919 : {
920 0 : return vbAvg_p->timeLast_p (row ());
921 : }
922 :
923 : Vector<Double>
924 0 : MsRowAvg::uvwFirst ()
925 : {
926 0 : return vbAvg_p->uvwFirst_p.column (row());
927 : }
928 :
929 :
930 : void
931 0 : MsRowAvg::setCounts (const Matrix<Int> & value)
932 : {
933 0 : Matrix<Int> & theCounts = countsCache_p.getCachedPlane (dynamic_cast<VbAvg *> (getVbi()), row());
934 0 : theCounts = value;
935 0 : }
936 :
937 : void
938 0 : MsRowAvg::setIntervalLast (Double value)
939 : {
940 0 : vbAvg_p->intervalLast_p (row ()) = value;
941 0 : }
942 :
943 :
944 : void
945 0 : MsRowAvg::setTimeFirst (Double value)
946 : {
947 0 : vbAvg_p->timeFirst_p (row ()) = value;
948 0 : }
949 :
950 : void
951 0 : MsRowAvg::setTimeLast (Double value)
952 : {
953 0 : vbAvg_p->timeLast_p (row ()) = value;
954 0 : }
955 :
956 0 : Double MsRowAvg::getNormalizationFactor()
957 : {
958 0 : return vbAvg_p->normalizationFactor_p (row ());
959 : }
960 :
961 0 : void MsRowAvg::setNormalizationFactor(Double normalizationFactor)
962 : {
963 0 : vbAvg_p->normalizationFactor_p (row ()) = normalizationFactor;
964 0 : }
965 :
966 0 : void MsRowAvg::accumulateNormalizationFactor(Double normalizationFactor)
967 : {
968 0 : vbAvg_p->normalizationFactor_p (row ()) += normalizationFactor;
969 0 : }
970 :
971 0 : VbAvg::VbAvg (const AveragingParameters & averagingParameters, ViImplementation2 * vii)
972 : : VisBufferImpl2 (vii, VbRekeyable),
973 0 : averagingInterval_p (averagingParameters.getAveragingInterval ()),
974 0 : averagingOptions_p (averagingParameters.getOptions()),
975 : averagingVii_p (vii),
976 : bufferToFill_p (0),
977 : complete_p (false),
978 : doing_p (), // all false until determined later on
979 : empty_p (true),
980 0 : maxTimeDistance_p (averagingParameters.getAveragingInterval() * (0.999)),
981 : // Shrink it just a bit for roundoff
982 0 : maxUvwDistanceSquared_p (pow(averagingParameters.getMaxUvwDistance(),2)),
983 : needIterationInfo_p (true),
984 : rowIdGenerator_p (0),
985 : sampleInterval_p (0),
986 : startTime_p (0),
987 0 : usingUvwDistance_p (averagingParameters.getOptions().contains (AveragingOptions::BaselineDependentAveraging))
988 0 : {}
989 :
990 : /**
991 : * Calculates the row index in the output buffer, given an averaged row, a baseline index
992 : * corresponding to this averaged row, and the magic "rowIdGenerator" of the VbAvg.
993 : *
994 : * @param nrows number of rows in the input buffer being averaged
995 : * @param baselineIndex index of the baseline being averaged into the rowAveraged
996 : * @param rowAveraged "accumulator" row being produced for the output buffer
997 : * @param rowIdGen the rowIdGenerator of the VbAvg, which increases (in a not so clean way)
998 : * for every new baseline inside every chunk
999 : */
1000 : Int
1001 0 : calcOutRowIdx(Int nrows, Int baselineIndex, const MsRowAvg* rowAveraged, Int rowIdGen)
1002 : {
1003 0 : auto nBasePresent = rowAveraged->nBaselinesPresent();
1004 : // check whether multiple time steps are being averaged into the output buffer
1005 : // (row blocking or similar feature is enabled)
1006 0 : const bool multitime = nrows > nBasePresent;
1007 :
1008 0 : if (!multitime) {
1009 : // There is only one time step, so the index must be simply the baseline index.
1010 : // with row blocking disabled, it doesn't seem to be possible to make sense out of
1011 : // rowIdgenerator_p for the purpose of this calculation -> skip the more general
1012 : // calculation from below and just use baseline index.
1013 0 : return baselineIndex;
1014 : }
1015 :
1016 : // the rowIdgenerator_p that we get in rowIdGen increases +1 for every new baseline.
1017 : // It is not really a proper (input) row id. After all baselines have been seen for a
1018 : // time step, the rows of the next time steps will get the same id.
1019 : // And to turn it into a valid output id we need the following
1020 : // "divide_with_roundup + multiply + baseline_index"
1021 0 : Int rowIdG_div_baselines_roundup = 0;
1022 0 : if (nBasePresent > 0)
1023 0 : rowIdG_div_baselines_roundup = (rowIdGen + nBasePresent - 1)/ nBasePresent;
1024 0 : const Int id = rowIdG_div_baselines_roundup * nBasePresent + baselineIndex;
1025 :
1026 0 : return id;
1027 : }
1028 :
1029 : void
1030 0 : VbAvg::accumulate (const VisBuffer2 * vb, const Subchunk & subchunk)
1031 : {
1032 0 : if (empty_p){
1033 0 : setupVbAvg (vb);
1034 : }
1035 :
1036 0 : if (needIterationInfo_p){
1037 0 : captureIterationInfo (bufferToFill_p, vb, subchunk);
1038 0 : needIterationInfo_p = false;
1039 : }
1040 :
1041 0 : assert (bufferToFill_p != 0);
1042 :
1043 0 : MsRowAvg * rowAveraged = getRowMutable (0);
1044 0 : MsRow * rowInput = vb->getRow (0);
1045 :
1046 0 : auto nrows = vb->nRows();
1047 0 : row2AvgRow_p.resize(nrows);
1048 0 : for (rownr_t row = 0; row < nrows; ++row){
1049 :
1050 0 : rowInput->changeRow (row);
1051 0 : Int baselineIndex = getBaselineIndex (rowInput);
1052 :
1053 0 : rowAveraged->changeRow (baselineIndex);
1054 :
1055 0 : accumulateOneRow (rowInput, rowAveraged, subchunk, row);
1056 :
1057 0 : row2AvgRow_p[row] = calcOutRowIdx(nrows, baselineIndex, rowAveraged,
1058 : rowIdGenerator_p);
1059 : }
1060 :
1061 0 : delete rowAveraged;
1062 0 : delete rowInput;
1063 :
1064 0 : }
1065 :
1066 : void
1067 0 : VbAvg::accumulateOneRow (MsRow * rowInput, MsRowAvg * rowAveraged, const Subchunk & subchunk,
1068 : Int iidx)
1069 : {
1070 0 : finalizeBaselineIfNeeded (rowInput, rowAveraged, subchunk);
1071 :
1072 0 : if (! rowAveraged->baselinePresent())
1073 : {
1074 :
1075 0 : initializeBaseline (rowInput, rowAveraged, subchunk);
1076 : }
1077 :
1078 : // bookkeeping - save for later that this input row is being averaged into the output row
1079 0 : rowAveraged->addInputRowIdx(iidx);
1080 :
1081 : // Accumulate data that is matrix-valued (e.g., vis, visModel, etc.).
1082 : // The computation of time centroid requires the use of the weight column
1083 : // adjusted for the flag cube. Get the before and after weight accumulation
1084 : // and the difference is the adjusted weight for this row.
1085 :
1086 0 : Vector<Double> adjustedWeights;
1087 0 : Bool rowFlagged = false;
1088 :
1089 0 : std::tie (rowFlagged, adjustedWeights) = accumulateCubeData (rowInput, rowAveraged);
1090 :
1091 0 : Double adjustedWeight = 0;
1092 0 : for (Int c = 0; c < nCorrelations(); c++){
1093 :
1094 0 : adjustedWeight += adjustedWeights (c);
1095 : }
1096 :
1097 : // Accumulate the non matrix-valued data
1098 0 : accumulateRowData (rowInput, rowAveraged, adjustedWeight, rowFlagged);
1099 0 : }
1100 :
1101 : //void
1102 : //VbAvg::accumulate (const VisBuffer2 * vb, const Subchunk & subchunk)
1103 : //{
1104 : // // "Add" the contents of this buffer to the accumulation.
1105 : //
1106 : // if (empty_p){
1107 : //
1108 : // // Initialize the buffer if this is the first time bit of data that it is
1109 : // // being used after either creation or clearing.
1110 : //
1111 : // prepareForFirstData (vb, subchunk);
1112 : //
1113 : // empty_p = false;
1114 : // }
1115 : //
1116 : // // Averaging can be computed as flagged or unflagged. If all the inputs to a
1117 : // // datum are flagged, then the averaged datum (e.g., a visibility point)
1118 : // // will also be flagged. For an unflagged averaged datum, it will represent
1119 : // // the average of all of the unflagged points for that datum. This is done
1120 : // // by assuming the accumulation is flagged and continuing to accumulate
1121 : // // flagged points until the first unflagged point for a datum is encountered;
1122 : // // when this happens the count is zeroed and the averaged datum's flag is cleared.
1123 : //
1124 : // // Loop over all of the rows in the VB. Map each one to a baseline and then accumulate
1125 : // // the values for each correlation,channel cell. Each row in the accumulating VB corresponds
1126 : // // to one baseline (i.e., pair of (antenna1, antenna2) where antenna1 <= antenna2).
1127 : //
1128 : // ThrowIf (vb->nRows() > nBaselines(),
1129 : // String::format ("Expected %d baselines in VisBuffer but it contained %d rows",
1130 : // nBaselines(), nRows()));
1131 : //
1132 : // for (Int row = 0; row < vb->nRows(); row ++){
1133 : //
1134 : // // Accumulate data for fields that are scalars (and uvw) in each row
1135 : //
1136 : // accumulateRowData (vb, row);
1137 : //
1138 : // // Accumulate data that is matrix-valued (e.g., vis, visModel, etc.)
1139 : //
1140 : // accumulateCubeData (vb, row);
1141 : // }
1142 : //}
1143 :
1144 : void
1145 0 : VbAvg::finalizeBufferFilling ()
1146 : {
1147 0 : bufferToFill_p->appendRowsComplete();
1148 0 : bufferToFill_p = 0; // decouple
1149 0 : }
1150 :
1151 : template<typename T>
1152 : inline void
1153 0 : VbAvg::accumulateElementForCube (const T * unweightedValue,
1154 : Float weight,
1155 : Bool zeroAccumulation,
1156 : T * accumulator)
1157 : {
1158 : // Update the sum for this model visibility cube element.
1159 :
1160 0 : if (zeroAccumulation){
1161 0 : * accumulator = (* unweightedValue) * weight;
1162 : }
1163 : else{
1164 0 : * accumulator += (* unweightedValue) * weight;
1165 : }
1166 0 : }
1167 :
1168 :
1169 : std::pair<Bool, Vector<Double> >
1170 0 : VbAvg::accumulateCubeData (MsRow * rowInput, MsRowAvg * rowAveraged)
1171 : {
1172 : // Accumulate the sums needed for averaging of cube data (e.g., visibility).
1173 :
1174 0 : const Matrix<Bool> & inputFlags = rowInput->flags ();
1175 0 : Matrix<Bool> & averagedFlags = rowAveraged->flagsMutable ();
1176 0 : Matrix<Int> counts = rowAveraged->counts ();
1177 0 : Vector<Bool> correlationFlagged = rowAveraged->correlationFlagsMutable ();
1178 :
1179 0 : AccumulationParameters accumulationParameters (rowInput, rowAveraged, doing_p);
1180 : // is a member variable to reduce memory allocations (jhj)
1181 :
1182 0 : IPosition shape = inputFlags.shape();
1183 0 : const Int nChannels = shape (1);
1184 0 : const Int nCorrelations = shape (0);
1185 :
1186 0 : Bool rowFlagged = true; // true if all correlations and all channels flagged
1187 :
1188 0 : for (Int channel = 0; channel < nChannels; channel ++){
1189 :
1190 0 : for (Int correlation = 0; correlation < nCorrelations; correlation ++){
1191 :
1192 : // Based on the current flag state of the accumulation and the current flag
1193 : // state of the correlation,channel, accumulate the data (or not). Accumulate
1194 : // flagged data until the first unflagged datum appears. Then restart the
1195 : // accumulation with that datum.
1196 :
1197 0 : Bool inputFlagged = inputFlags (correlation, channel);
1198 0 : if (rowFlagged && ! inputFlagged){
1199 0 : rowFlagged = false;
1200 : }
1201 : //rowFlagged = rowFlagged && inputFlagged;
1202 0 : Bool accumulatorFlagged = averagedFlags (correlation, channel);
1203 :
1204 0 : if (! accumulatorFlagged && inputFlagged){
1205 0 : accumulationParameters.incrementCubePointers();
1206 0 : continue;// good accumulation, bad data so toss it.
1207 : }
1208 :
1209 : // If changing from flagged to unflagged for this cube element, reset the
1210 : // accumulation count to 1; otherwise increment the count.
1211 :
1212 0 : Bool flagChange = (accumulatorFlagged && ! inputFlagged);
1213 0 : Bool zeroAccumulation = flagChange || counts (correlation, channel) == 0;
1214 :
1215 0 : if (flagChange){
1216 0 : averagedFlags (correlation, channel) = false;
1217 : }
1218 :
1219 0 : if (zeroAccumulation){
1220 0 : counts (correlation, channel) = 1;
1221 : }
1222 : else{
1223 0 : counts (correlation, channel) += 1;
1224 : }
1225 :
1226 : // Accumulate the sum for each cube element
1227 :
1228 0 : accumulateElementForCubes (accumulationParameters,
1229 : zeroAccumulation); // zeroes out accumulation
1230 :
1231 0 : accumulationParameters.incrementCubePointers();
1232 :
1233 : // Update correlation Flag
1234 :
1235 0 : if (correlationFlagged (correlation) && ! inputFlagged){
1236 0 : correlationFlagged (correlation) = false;
1237 : }
1238 : }
1239 : }
1240 :
1241 0 : Vector<Double> adjustedWeight = Vector<Double> (nCorrelations, 1);
1242 0 : if (doing_p.correctedData_p)
1243 : {
1244 0 : for (Int correlation = 0; correlation < nCorrelations; correlation ++)
1245 : {
1246 0 : adjustedWeight(correlation) = rowInput->weight(correlation);
1247 : }
1248 : }
1249 0 : else if (doing_p.observedData_p || doing_p.floatData_p)
1250 : {
1251 0 : for (Int correlation = 0; correlation < nCorrelations; correlation ++)
1252 : {
1253 0 : adjustedWeight(correlation) = AveragingTvi2::sigmaToWeight(rowInput->sigma (correlation));
1254 : }
1255 : }
1256 :
1257 0 : return std::make_pair (rowFlagged, adjustedWeight);
1258 : }
1259 :
1260 :
1261 : inline void
1262 0 : VbAvg::accumulateElementForCubes (AccumulationParameters & accumulationParameters,
1263 : Bool zeroAccumulation)
1264 : {
1265 :
1266 : // NOTE: The channelized flag check comes from the calling context (continue statement)
1267 0 : float weightCorrected = 1.0f;
1268 0 : float weightObserved = 1.0f;
1269 0 : const float One = 1.0f;
1270 :
1271 0 : if (doing_p.correctedData_p)
1272 : {
1273 : // The weight corresponding to CORRECTED_DATA is that stored in WEIGHT
1274 0 : weightCorrected = * accumulationParameters.weightSpectrumIn ();
1275 :
1276 :
1277 : // Accumulate weighted average contribution (normalization will come at the end)
1278 0 : accumulateElementForCube ( accumulationParameters.correctedIn (),
1279 : weightCorrected, zeroAccumulation,
1280 : accumulationParameters.correctedOut ());
1281 :
1282 : // The weight resulting from weighted average is the sum of the weights
1283 0 : accumulateElementForCube ( & weightCorrected,
1284 : One, zeroAccumulation,
1285 : accumulationParameters.weightSpectrumOut ());
1286 : }
1287 :
1288 0 : if (doing_p.observedData_p)
1289 : {
1290 : // The weight corresponding to DATA is that derived from the rms stored in SIGMA
1291 : // This has to
1292 0 : weightObserved = AveragingTvi2::sigmaToWeight(* accumulationParameters.sigmaSpectrumIn ());
1293 :
1294 : // Accumulate weighted average contribution (normalization will come at the end)
1295 :
1296 0 : accumulateElementForCube ( accumulationParameters.observedIn (),
1297 : weightObserved, zeroAccumulation,
1298 : accumulationParameters.observedOut ());
1299 :
1300 0 : if (not doing_p.correctedData_p)
1301 : {
1302 : // The weight resulting from weighted average is the sum of the weights
1303 0 : accumulateElementForCube ( & weightObserved,
1304 : One, zeroAccumulation,
1305 : accumulationParameters.weightSpectrumOut ());
1306 : }
1307 :
1308 : // This will always create a sigma spectrum column which is not empty.
1309 : // This is useful in particular if not doing_p.correctedData_p but doing_p.modelData_p,
1310 : // so that modelData can be properly divided by sigmaSpectrumOut in finalizeCubeData
1311 : // We store the accumulated weight in sigmaSpectrumOut pending of
1312 : // - normalization
1313 : // - SIGMA = 1/sqrt(WEIGHT) in-place transformation
1314 0 : accumulateElementForCube (&weightObserved,
1315 : One, zeroAccumulation,
1316 : accumulationParameters.sigmaSpectrumOut ());
1317 :
1318 : }
1319 :
1320 : // For model data is less clear what to do, what in order to convert to
1321 : // split we use WEIGHT if averaging CORRECTED_DATA and SIGMA if avg. DATA.
1322 : // Finally we use WEIGHT by default when averaging MODEL_DATA only
1323 0 : if (doing_p.modelData_p)
1324 : {
1325 0 : if (doing_p.correctedData_p)
1326 : {
1327 0 : accumulateElementForCube ( accumulationParameters.modelIn (),
1328 : weightCorrected, zeroAccumulation,
1329 : accumulationParameters.modelOut ());
1330 : }
1331 0 : else if (doing_p.observedData_p)
1332 : {
1333 0 : accumulateElementForCube ( accumulationParameters.modelIn (),
1334 : weightObserved, zeroAccumulation,
1335 : accumulationParameters.modelOut ());
1336 : }
1337 : else
1338 : {
1339 0 : accumulateElementForCube ( accumulationParameters.modelIn (),
1340 : One, zeroAccumulation,
1341 : accumulationParameters.modelOut ());
1342 :
1343 : // When doing MODEL_DATA only the accumulated weight spectrum should just represent counts
1344 0 : accumulateElementForCube ( & One,
1345 : 1.0f, zeroAccumulation,
1346 : accumulationParameters.weightSpectrumOut ());
1347 : }
1348 : }
1349 :
1350 0 : if (doing_p.floatData_p)
1351 : {
1352 :
1353 : // The weight corresponding to FLOAT_DATA is that derived from the rms stored in SIGMA
1354 0 : weightObserved = AveragingTvi2::sigmaToWeight(* accumulationParameters.sigmaSpectrumIn ());
1355 :
1356 : // Accumulate weighted average contribution (normalization will come at the end)
1357 0 : accumulateElementForCube ( accumulationParameters.floatDataIn (),
1358 : weightObserved, zeroAccumulation,
1359 : accumulationParameters.floatDataOut ());
1360 :
1361 : // The weight resulting from weighted average is the sum of the weights
1362 0 : accumulateElementForCube ( & weightObserved,
1363 : 1.0f, zeroAccumulation,
1364 : accumulationParameters.weightSpectrumOut ());
1365 : }
1366 :
1367 0 : return;
1368 : }
1369 :
1370 :
1371 :
1372 : template <typename T>
1373 : T
1374 0 : VbAvg::accumulateRowDatum (const T & averagedValue, const T & inputValue, Bool resetAverage)
1375 : {
1376 0 : if (resetAverage){
1377 0 : return inputValue;
1378 : }
1379 : else{
1380 0 : return inputValue + averagedValue;
1381 : }
1382 : }
1383 :
1384 : void
1385 0 : VbAvg::accumulateRowData (MsRow * rowInput, MsRowAvg * rowAveraged,
1386 : Double adjustedWeight, Bool rowFlagged)
1387 : {
1388 :
1389 : // Grab working copies of the values to be accumulated.
1390 :
1391 0 : Bool accumulatorRowFlagged = rowAveraged->isRowFlagged();
1392 0 : Bool flagChange = accumulatorRowFlagged != rowFlagged || // actual change
1393 0 : rowAveraged->countsBaseline() == 0; // first time
1394 :
1395 0 : if (! accumulatorRowFlagged && rowFlagged){
1396 : // good accumulation, bad data --> skip it
1397 : }
1398 : else{
1399 :
1400 : // Update the row's accumulations; if the flagChanged then zero out the
1401 : // previous (flagged) accumulation first.
1402 :
1403 0 : rowAveraged->setCountsBaseline (accumulateRowDatum (rowAveraged->countsBaseline(),
1404 0 : 1,
1405 : flagChange));
1406 :
1407 : // The WEIGHT column is handled under accumulateCubeData because of the
1408 : // interrelationship between weight and weightSpectrum. The SIGMA column is
1409 : // handled in finalizeBaseline for similar reasons.
1410 :
1411 0 : accumulatorRowFlagged = accumulatorRowFlagged && rowFlagged;
1412 0 : rowAveraged->setRowFlag (accumulatorRowFlagged);
1413 :
1414 0 : rowAveraged->setExposure (accumulateRowDatum (rowAveraged->exposure(),
1415 0 : rowInput->exposure (),
1416 0 : flagChange));
1417 :
1418 : // While accumulating flagged values, the weights will be zero, so accumulate
1419 : // an arithmetic average until the accumulator becomes unflagged.
1420 :
1421 : Double weightToUse;
1422 :
1423 0 : if (accumulatorRowFlagged){
1424 0 : weightToUse = 1;
1425 : }
1426 : else{
1427 0 : weightToUse = adjustedWeight;
1428 : }
1429 :
1430 0 : if (flagChange){
1431 0 : rowAveraged->setNormalizationFactor(0.0);
1432 : }
1433 :
1434 0 : rowAveraged->accumulateNormalizationFactor(weightToUse);
1435 :
1436 0 : Double weightedTC = (rowInput->timeCentroid() - rowAveraged->timeFirst()) * weightToUse;
1437 0 : rowAveraged->setTimeCentroid (accumulateRowDatum (rowAveraged->timeCentroid(),
1438 : weightedTC,
1439 0 : flagChange));
1440 :
1441 0 : Vector<Double> weightedUvw = rowInput->uvw() * weightToUse;
1442 :
1443 0 : rowAveraged->setUvw (accumulateRowDatum (rowAveraged->uvw (),
1444 : weightedUvw,
1445 0 : flagChange));
1446 :
1447 : // Capture a couple pieces of data
1448 :
1449 0 : rowAveraged->setTimeLast (rowInput->time());
1450 :
1451 0 : rowAveraged->setIntervalLast (rowInput->interval());
1452 : }
1453 0 : }
1454 :
1455 : //Vector<Float>
1456 : //VbAvg::adjustWeightForFlags (MsRow * rowInput)
1457 : //{
1458 : // Matrix<Bool> flags = rowInput->flags();
1459 : // Vector<Float> adjustedWeight = rowInput->weight();
1460 : //
1461 : // for (Int correlation = 0; correlation < nCorrelations(); correlation++){
1462 : //
1463 : // // Sum up the number of good channels in this correlation
1464 : //
1465 : // Int sum = 0;
1466 : //
1467 : // for (Int channel = 0; channel < nChannels(); channel ++){
1468 : //
1469 : // if (! flags (correlation, channel)){
1470 : //
1471 : // sum ++;
1472 : // }
1473 : // }
1474 : //
1475 : // // Adjust the weight by multiplying by the fraction of good channels.
1476 : //
1477 : // Float factor = ((float) sum) / nChannels();
1478 : // adjustedWeight [correlation] *= factor;
1479 : // }
1480 : //
1481 : // return adjustedWeight;
1482 : //}
1483 :
1484 : const Cube<Int> &
1485 0 : VbAvg::counts () const
1486 : {
1487 0 : return counts_p;
1488 : }
1489 :
1490 :
1491 : void
1492 0 : VbAvg::copyIdValues (MsRow * rowInput, MsRowAvg * rowAveraged)
1493 : {
1494 0 : rowAveraged->setAntenna1 (rowInput->antenna1());
1495 0 : rowAveraged->setAntenna2 (rowInput->antenna2());
1496 0 : rowAveraged->setArrayId (rowInput->arrayId());
1497 0 : rowAveraged->setDataDescriptionId (rowInput->dataDescriptionId());
1498 0 : rowAveraged->setFeed1 (rowInput->feed1());
1499 0 : rowAveraged->setFeed2 (rowInput->feed2());
1500 0 : rowAveraged->setFieldId (rowInput->fieldId());
1501 0 : rowAveraged->setProcessorId (rowInput->processorId());
1502 0 : rowAveraged->setScanNumber (rowInput->scanNumber());
1503 0 : rowAveraged->setSpectralWindow (rowInput->spectralWindow());
1504 0 : rowAveraged->setObservationId (rowInput->observationId());
1505 0 : rowAveraged->setStateId (rowInput->stateId());
1506 0 : }
1507 :
1508 : void
1509 0 : VbAvg::copyIdValue (Int inputId, Int & averagedId)
1510 : {
1511 0 : if (averagedId < 0){
1512 0 : averagedId = inputId;
1513 : }
1514 0 : }
1515 :
1516 : Bool
1517 0 : VbAvg::empty () const
1518 : {
1519 0 : return empty_p;
1520 : }
1521 :
1522 : Int
1523 0 : VbAvg::getBaselineIndex (const MsRow * msRow) const
1524 : {
1525 : // Lookup the baseline index using the prebuilt lookup table.
1526 : //
1527 : // The baseline index is the index in the sequence
1528 : // {(0,0),(1,0),(1,1),(2,0),(2,1),(2,2), ...} (i.e., the index in a
1529 : // 1-d representation of the lower right half of the square matrix of size
1530 : // nAntennas).
1531 : //
1532 : // This handles the case where the baseline ordering in the input VB might
1533 : // shift from VB to VB.
1534 :
1535 0 : const Int antenna1 = msRow->antenna1 ();
1536 0 : const Int antenna2 = msRow->antenna2 ();
1537 0 : const Int spw = msRow->spectralWindow ();
1538 :
1539 0 : const Int index = baselineIndex_p (antenna1, antenna2, spw);
1540 :
1541 0 : return index;
1542 : }
1543 :
1544 : Int
1545 0 : VbAvg::getBaselineIndex (Int antenna1, Int antenna2, Int spw) const
1546 : {
1547 0 : const Int index = baselineIndex_p (antenna1, antenna2, spw);
1548 :
1549 0 : return index;
1550 : }
1551 :
1552 : void
1553 0 : VbAvg::finalizeAverages ()
1554 : {
1555 0 : if (empty_p){
1556 0 : return; // nothing to finalize
1557 : }
1558 :
1559 0 : MsRowAvg * msRowAvg = getRowMutable (0);
1560 :
1561 0 : for (Int baseline = 0; baseline < nBaselines(); baseline ++){
1562 :
1563 0 : msRowAvg->changeRow (baseline);
1564 :
1565 0 : if (msRowAvg->baselinePresent()){
1566 0 : finalizeBaseline (msRowAvg);
1567 : }
1568 :
1569 : }
1570 :
1571 0 : delete msRowAvg;
1572 :
1573 0 : empty_p = true;
1574 : }
1575 :
1576 : void
1577 0 : VbAvg::finalizeBaseline (MsRowAvg * msRowAvg)
1578 : {
1579 : // Software is no longer supposed to rely on the row flag.
1580 : // Setting it to false insures that legacy software will
1581 : // have to look at the flag cubes.
1582 :
1583 0 : msRowAvg->setRowFlag(false);
1584 :
1585 0 : finalizeCubeData (msRowAvg);
1586 :
1587 0 : finalizeRowData (msRowAvg);
1588 :
1589 0 : transferBaseline (msRowAvg);
1590 0 : }
1591 :
1592 : // Functor to divide variables of possibly different types.
1593 : // This is unlike std::divides which requires equal types.
1594 : template <typename L, typename R=L, typename RES=L>
1595 : struct DividesNonZero : public std::binary_function<L,R,RES>
1596 : {
1597 0 : RES operator() (const L& x, const R& y) const
1598 : {
1599 0 : { return y > 0? RES(x)/y : RES(x); }
1600 : }
1601 : };
1602 :
1603 :
1604 : void
1605 0 : VbAvg::finalizeCubeData (MsRowAvg * msRowAvg)
1606 : {
1607 : // Divide each of the data cubes in use by the sum of the appropriate weights.
1608 :
1609 : typedef DividesNonZero <Complex, Float, Complex> DivideOp;
1610 : DivideOp op;
1611 :
1612 0 : if (doing_p.correctedData_p)
1613 : {
1614 0 : Matrix<Complex> corrected = msRowAvg->correctedMutable();
1615 0 : arrayTransformInPlace<Complex, Float, DivideOp > (corrected,msRowAvg->weightSpectrum (), op);
1616 : }
1617 :
1618 0 : if (doing_p.observedData_p)
1619 : {
1620 0 : Matrix<Complex> observed = msRowAvg->observedMutable();
1621 0 : if (not doing_p.correctedData_p)
1622 0 : arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->weightSpectrum (), op);
1623 : else
1624 0 : arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->sigmaSpectrum (), op);
1625 : }
1626 :
1627 0 : if (doing_p.modelData_p)
1628 : {
1629 0 : Matrix<Complex> model = msRowAvg->modelMutable();
1630 :
1631 0 : if (doing_p.correctedData_p)
1632 0 : arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->weightSpectrum (), op);
1633 0 : else if (doing_p.observedData_p)
1634 0 : arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->sigmaSpectrum (), op);
1635 : else
1636 0 : arrayTransformInPlace<Complex, Int, DivideOp > (model,msRowAvg->counts (), op);
1637 : }
1638 :
1639 0 : if (doing_p.floatData_p)
1640 : {
1641 : typedef Divides <Float, Float, Float> DivideOpFloat;
1642 : DivideOpFloat opFloat;
1643 :
1644 0 : Matrix<Float> visCubeFloat = msRowAvg->singleDishDataMutable();
1645 0 : arrayTransformInPlace<Float, Float, DivideOpFloat > (visCubeFloat,msRowAvg->weightSpectrum (), opFloat);
1646 : }
1647 :
1648 :
1649 0 : return;
1650 : }
1651 :
1652 : void
1653 0 : VbAvg::finalizeRowData (MsRowAvg * msRow)
1654 : {
1655 0 : Int n = msRow->countsBaseline ();
1656 :
1657 : // Obtain row-level WEIGHT by calculating the mean of WEIGHT_SPECTRUM
1658 : // msRow->setWeight(partialMedians(msRow->weightSpectrum(),IPosition(1,1),true));
1659 0 : msRow->setWeight(AveragingTvi2::average(msRow->weightSpectrum(),msRow->flags()));
1660 :
1661 : // If doing both DATA and CORRECTED_DATA then SIGMA_SPECTRUM contains the weight
1662 : // (not sigma) accumulation for DATA, and we have to derive SIGMA from it
1663 0 : if (doing_p.correctedData_p and doing_p.observedData_p)
1664 : {
1665 : // jagonzal: SIGMA is not derived from the mean of SIGMA_SPECTRUM but from the mean
1666 : // of the WEIGHT format of SIGMA_SPECTRUM turned into SIGMA by using 1/pow(weight,2)
1667 : // Vector<Float> weight = partialMedians(msRow->sigmaSpectrum(),IPosition(1,1),true);
1668 0 : Vector<Float> weight = AveragingTvi2::average(msRow->sigmaSpectrum(),msRow->flags());
1669 0 : arrayTransformInPlace (weight, AveragingTvi2::weightToSigma);
1670 0 : msRow->setSigma (weight);
1671 :
1672 : // Now convert the DATA weight accumulation stored in sigmaSpectrum into the real SIGMA_SPECTRUM
1673 : // TODO: This should happen only if we are writing out SIGMA_SPECTRUM but
1674 : // multiple column operation is rare and might be forbidden in the future
1675 0 : Matrix<Float> sigmaSpectrun = msRow->sigmaSpectrum(); // Reference copy
1676 0 : arrayTransformInPlace (sigmaSpectrun, AveragingTvi2::weightToSigma);
1677 : }
1678 : // Otherwise (doing only DATA/FLOAT_DATA or CORRECTED_DATA) we can derive SIGMA from WEIGHT directly
1679 0 : else if ( not doing_p.onlymetadata_p)
1680 : {
1681 : // jagonzal: SIGMA is not derived from the mean of SIGMA_SPECTRUM
1682 : // but from WEIGHT turned into SIGMA by using 1/pow(weight,2)
1683 0 : Vector<Float> sigma = msRow->sigma(); // Reference copy
1684 0 : sigma = msRow->weight(); // Normal copy (transfer Weight values to Sigma)
1685 0 : arrayTransformInPlace (sigma, AveragingTvi2::weightToSigma);
1686 :
1687 : // Derive SIGMA_SPECTRUM from computed WEIGHT_SPECTRUM
1688 0 : Matrix<Float> sigmaSpectrun = msRow->sigmaSpectrum(); // Reference copy
1689 0 : sigmaSpectrun = msRow->weightSpectrum(); // Normal copy (transfer WeightSpectrum values to SigmaSpectrum)
1690 0 : arrayTransformInPlace (sigmaSpectrun, AveragingTvi2::weightToSigma);
1691 : }
1692 :
1693 : // Get the normalization factor for this baseline, containing
1694 : // the accumulation of row-level) weight contributions
1695 0 : Double weight = msRow->getNormalizationFactor();
1696 :
1697 0 : if (n != 0){
1698 :
1699 0 : if (weight == 0){
1700 :
1701 : // The weights are all zero so compute an arithmetic average
1702 : // so that a somewhat value can go into these columns (i.e. rather than NaN).
1703 :
1704 0 : weight = msRow->countsBaseline();
1705 : }
1706 :
1707 0 : msRow->setTimeCentroid (msRow->timeCentroid() / weight + msRow->timeFirst());
1708 :
1709 0 : msRow->setUvw (msRow->uvw() / weight);
1710 :
1711 : // Exposure is a simple sum, not an average so it is already
1712 : // done at this point.
1713 : }
1714 :
1715 : // Fill in the time and the interval
1716 : //
1717 : // The time of a sample is centered over the integration time period.
1718 : // Thus find the center of the averaged interval it is necessary to
1719 : // slide it back by 1/2 an interval.
1720 :
1721 0 : Double dT = msRow->timeLast () - msRow->timeFirst();
1722 :
1723 0 : Double centerOfInterval = msRow->timeFirst () + dT / 2;
1724 :
1725 0 : msRow->setTime (centerOfInterval);
1726 :
1727 0 : if (dT != 0){
1728 :
1729 : // The interval is the center-to-center time + half of the intervals of
1730 : // the first and the last sample (if dT == 0 then the interval is
1731 : // already the interval of the first sample and is correct)
1732 :
1733 0 : Double interval = dT + msRow->interval() / 2 + msRow->intervalLast() / 2;
1734 0 : msRow->setInterval (interval);
1735 : }
1736 0 : }
1737 :
1738 : void
1739 0 : VbAvg::finalizeBaselineIfNeeded (MsRow * rowInput, MsRowAvg * rowAveraged, const Subchunk & /*subchunk*/)
1740 : {
1741 0 : if (! rowAveraged->baselinePresent()){
1742 0 : return;
1743 : }
1744 :
1745 : // Finalization is needed if either the uvw distance or the time distance between the input
1746 : // baseline and the averaged baseline is above the maximum
1747 :
1748 0 : Bool needed = usingUvwDistance_p;
1749 :
1750 0 : if (needed) {
1751 0 : Double deltaUvw = distanceSquared (rowInput->uvw(), rowAveraged->uvwFirst ());
1752 0 : needed = deltaUvw > maxUvwDistanceSquared_p;
1753 : }
1754 :
1755 0 : needed = needed || (rowInput->time() - rowAveraged->timeFirst()) > maxTimeDistance_p;
1756 :
1757 0 : if (needed){
1758 :
1759 : // Finalize the data so that the final averaging products and then move them to
1760 : // output buffer.
1761 :
1762 0 : finalizeBaseline (rowAveraged);
1763 : }
1764 : }
1765 :
1766 : MsRowAvg *
1767 0 : VbAvg::getRow (Int row) const
1768 : {
1769 0 : return new MsRowAvg (row, this);
1770 : }
1771 :
1772 : MsRowAvg *
1773 0 : VbAvg::getRowMutable (Int row)
1774 : {
1775 0 : return new MsRowAvg (row, this);
1776 : }
1777 :
1778 : void
1779 0 : VbAvg::initializeBaseline (MsRow * rowInput, MsRowAvg * rowAveraged,
1780 : const Subchunk & /*subchunk*/)
1781 : {
1782 0 : copyIdValues (rowInput, rowAveraged);
1783 :
1784 : // Size and zero out the counters
1785 :
1786 0 : rowAveraged->setInterval (rowInput->interval()); // capture first one
1787 0 : rowAveraged->setTimeFirst (rowInput->time());
1788 0 : rowAveraged->setTimeLast (rowInput->time());
1789 0 : rowAveraged->uvwFirst () = rowInput->uvw ();
1790 :
1791 0 : rowAveraged->setCountsBaseline (0);
1792 :
1793 0 : IPosition shape = rowInput->flags().shape();
1794 0 : Int nCorrelations = shape (0);
1795 0 : Int nChannels = shape (1);
1796 :
1797 0 : rowAveraged->setCounts (zeroInt_p.getMatrix (nCorrelations, nChannels, 0));
1798 0 : rowAveraged->setWeight (Vector<Float> (nCorrelations, 0));
1799 0 : rowAveraged->setTimeCentroid (0.0);
1800 :
1801 0 : if (doing_p.weightSpectrumOut_p){
1802 0 : rowAveraged->setWeightSpectrum (zeroFloat_p.getMatrix (nCorrelations, nChannels, 0));
1803 : }
1804 :
1805 0 : if (doing_p.sigmaSpectrumOut_p){
1806 0 : rowAveraged->setSigmaSpectrum (zeroFloat_p.getMatrix (nCorrelations, nChannels, 0));
1807 : }
1808 :
1809 : // VisBufferComponents2 exclusions =
1810 : // VisBufferComponents2::these(VisibilityObserved, VisibilityCorrected,
1811 : // VisibilityModel, CorrType, JonesC, Unknown);
1812 : //
1813 : // cacheResizeAndZero(exclusions);
1814 :
1815 : // Flag everything to start with
1816 :
1817 0 : rowAveraged->setRowFlag (true); // only for use during row-value accumulation
1818 0 : rowAveraged->setFlags(trueBool_p.getMatrix (nCorrelations, nChannels, true));
1819 0 : rowAveraged->correlationFlagsMutable() = Vector<Bool> (nCorrelations, true);
1820 :
1821 0 : rowAveraged->setBaselinePresent(true);
1822 :
1823 0 : rowAveraged->setNormalizationFactor(0.0);
1824 0 : }
1825 :
1826 :
1827 : Bool
1828 0 : VbAvg::isComplete () const
1829 : {
1830 0 : return complete_p;
1831 : }
1832 :
1833 : Bool
1834 0 : VbAvg::isUsingUvwDistance () const
1835 : {
1836 0 : return usingUvwDistance_p;
1837 : }
1838 :
1839 :
1840 : //void
1841 : //VbAvg::markEmpty ()
1842 : //{
1843 : // empty_p = true;
1844 : // complete_p = false;
1845 : //}
1846 :
1847 : Int
1848 0 : VbAvg::nBaselines () const
1849 : {
1850 0 : return countsBaseline_p.nelements();
1851 : }
1852 :
1853 : Int
1854 0 : VbAvg::nSpectralWindowsInBuffer () const
1855 : {
1856 0 : const Vector<Int> & windows = spectralWindows();
1857 0 : set <Int> s;
1858 :
1859 0 : for (uInt i = 0; i < windows.nelements(); i ++){
1860 0 : s.insert (windows(i));
1861 : }
1862 :
1863 0 : return (Int) s.size();
1864 :
1865 : }
1866 :
1867 :
1868 : void
1869 0 : VbAvg::captureIterationInfo (VisBufferImpl2 * dstVb, const VisBuffer2 * srcVb,
1870 : const Subchunk & subchunk)
1871 : {
1872 0 : dstVb->setIterationInfo (srcVb->msId(),
1873 0 : srcVb->msName(),
1874 0 : srcVb->isNewMs(),
1875 0 : srcVb->isNewArrayId(),
1876 0 : srcVb->isNewFieldId(),
1877 0 : srcVb->isNewSpectralWindow(),
1878 : subchunk,
1879 0 : srcVb->getCorrelationTypes (),
1880 0 : srcVb->getCorrelationTypesDefined(),
1881 0 : srcVb->getCorrelationTypesSelected(),
1882 0 : CountedPtr <WeightScaling> ( ));
1883 :
1884 : // Request info from the VB which will cause it to be filled
1885 : // into cache from the input VII at its current position.
1886 :
1887 0 : dstVb->setRekeyable(true);
1888 0 : dstVb->setShape(srcVb->nCorrelations(), srcVb->nChannels(), nBaselines(), false);
1889 : // Do not clear the cache since we're resuing the storage
1890 :
1891 0 : dstVb->phaseCenter();
1892 0 : dstVb->nAntennas();
1893 0 : dstVb->correlationTypes();
1894 0 : dstVb->polarizationFrame();
1895 0 : dstVb->polarizationId();
1896 0 : }
1897 :
1898 : //void
1899 : //VbAvg::prepareForFirstData (const VisBuffer2 * vb, const Subchunk & subchunk)
1900 : //{
1901 : // startTime_p = vb->time() (0);
1902 : // sampleInterval_p = vb->timeInterval() (0);
1903 : //
1904 : // Int nAntennas = vb->nAntennas();
1905 : // Int nSpw = vb->getVi()->nSpectralWindows();
1906 : // Int nBaselines = ((nAntennas * (nAntennas + 1)) / 2) * nSpw;
1907 : //
1908 : // // Size and zero out the counters
1909 : //
1910 : // timeFirst_p = Vector<Double> (nBaselines, vb->time() (0));
1911 : // timeLast_p = Vector<Double> (nBaselines, vb->time() (0));
1912 : // uvwFirst_p = Vector<Double> (nBaselines, vb->uvw().column(0));
1913 : //
1914 : // countsBaseline_p = Vector<Int> (nBaselines, 0);
1915 : // counts_p = Cube<Int> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1916 : // weightSum_p = Cube<Float> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1917 : // if (doing_p.sigmaSpectrum_p){
1918 : // weightCorrectedSum_p = Cube<Float> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1919 : // }
1920 : //
1921 : // baselineIndex_p.configure (nAntennas, nSpw, vb);
1922 : //
1923 : // // Reshape the inherited members from VisBuffer2
1924 : //
1925 : // captureIterationInfo (vb, subchunk);
1926 : //
1927 : // setShape (vb->nCorrelations(), vb->nChannels(), nBaselines, false);
1928 : //
1929 : // VisBufferComponents2 exclusions =
1930 : // VisBufferComponents2::these(VisibilityObserved, VisibilityCorrected,
1931 : // VisibilityModel, CorrType, JonesC, Unknown);
1932 : // cacheResizeAndZero(exclusions);
1933 : //
1934 : // prepareIds (vb);
1935 : //
1936 : // // Flag everything to start with
1937 : //
1938 : // setFlagCube (Cube<Bool> (vb->nCorrelations(), vb->nChannels(), nBaselines, true));
1939 : // setFlagRow (Vector<Bool> (nBaselines, true));
1940 : //
1941 : // complete_p = false;
1942 : //}
1943 :
1944 : void
1945 0 : VbAvg::prepareIds (const VisBuffer2 * vb)
1946 : {
1947 : // Set these row ID values to indicate they are unknown
1948 :
1949 0 : Vector<Int> minusOne (nBaselines(), -1);
1950 :
1951 0 : setAntenna1 (minusOne);
1952 0 : setAntenna2 (minusOne);
1953 0 : setDataDescriptionIds (minusOne);
1954 0 : setFeed1 (minusOne);
1955 0 : setFeed2 (minusOne);
1956 0 : setProcessorId (minusOne);
1957 0 : setScan (minusOne);
1958 0 : setObservationId (minusOne);
1959 0 : setSpectralWindows (minusOne);
1960 0 : setStateId (minusOne);
1961 :
1962 : // Copy the value from the input VB
1963 :
1964 0 : Vector<Int> tmp (nBaselines(), vb->arrayId()(0));
1965 :
1966 0 : setArrayId (tmp);
1967 :
1968 0 : tmp = vb->dataDescriptionIds()(0);
1969 0 : setDataDescriptionIds (tmp);
1970 :
1971 0 : tmp = vb->fieldId()(0);
1972 0 : setFieldId (tmp);
1973 0 : }
1974 :
1975 : //void
1976 : //VbAvg::removeMissingBaselines ()
1977 : //{
1978 : // // Some baselines may not be present in the portion of the input data
1979 : // // that made up this average. However, this is not known until after
1980 : // // all of the data is seen. Thus at finalization time these bogus
1981 : // // baselines should be removed from the average so as not to pass
1982 : // // flagged but zero-exposure baselines into the output.
1983 : //
1984 : //
1985 : // Vector<Int> rowsToDelete (nBaselines());
1986 : //
1987 : // Int nBaselinesDeleted = 0;
1988 : //
1989 : // for (Int baseline = 0; baseline < nBaselines(); baseline ++){
1990 : //
1991 : // if (countsBaseline_p (baseline) == 0){
1992 : // rowsToDelete (nBaselinesDeleted) = baseline;
1993 : // nBaselinesDeleted ++;
1994 : // }
1995 : // }
1996 : //
1997 : // rowsToDelete.resize (nBaselinesDeleted, true);
1998 : //
1999 : // deleteRows (rowsToDelete);
2000 : //}
2001 :
2002 : void
2003 0 : VbAvg::setBufferToFill(VisBufferImpl2 * vb)
2004 : {
2005 0 : bufferToFill_p = vb;
2006 :
2007 : // Set flag so that iteration information will be captured into
2008 : // this VB from the first input VB.
2009 :
2010 0 : needIterationInfo_p = true;
2011 0 : }
2012 :
2013 : void
2014 0 : VbAvg::setupVbAvg (const VisBuffer2 * vb)
2015 : {
2016 : // Configure the index
2017 :
2018 0 : Int nAntennas = vb->nAntennas();
2019 :
2020 : // This is a kluge to allow multiple spectral windows (of the same shape)
2021 : // to be combined into a single VB. This really shouldn't be allowed!!!
2022 :
2023 0 : set<uInt> spwInVb;
2024 :
2025 0 : for (rownr_t i = 0; i < vb->nRows(); i++){
2026 0 : spwInVb.insert (vb->dataDescriptionIds()(i));
2027 : }
2028 :
2029 0 : uInt nSpwInVb = spwInVb.size();
2030 :
2031 0 : Int nSpw = averagingVii_p->nSpectralWindows();
2032 :
2033 0 : baselineIndex_p.configure (nAntennas, nSpw, vb);
2034 :
2035 0 : Int nBaselines = ((nAntennas * (nAntennas + 1)) / 2)* nSpwInVb;
2036 :
2037 0 : setShape (vb->nCorrelations(), vb->nChannels(), nBaselines);
2038 :
2039 0 : setupArrays (vb->nCorrelations(), vb->nChannels(), nBaselines);
2040 :
2041 0 : baselinePresent_p.resize(nBaselines);
2042 0 : baselinePresent_p = False;
2043 :
2044 0 : normalizationFactor_p.resize(nBaselines);
2045 0 : normalizationFactor_p = 0.0;
2046 :
2047 0 : empty_p = false;
2048 0 : }
2049 :
2050 : void
2051 0 : VbAvg::setupArrays (Int nCorrelations, Int nChannels, Int nBaselines)
2052 : {
2053 :
2054 0 : setShape (nCorrelations, nChannels, nBaselines);
2055 :
2056 : // Start out with all of the array-valued components except the
2057 : // optional ones.
2058 :
2059 : VisBufferComponents2 including =
2060 : VisBufferComponents2::these ({VisBufferComponent2::Antenna1,
2061 : VisBufferComponent2::Antenna2,
2062 : VisBufferComponent2::ArrayId,
2063 : VisBufferComponent2::CorrType,
2064 : VisBufferComponent2::DataDescriptionIds,
2065 : VisBufferComponent2::Exposure,
2066 : VisBufferComponent2::Feed1,
2067 : VisBufferComponent2::Feed2,
2068 : VisBufferComponent2::FieldId,
2069 : VisBufferComponent2::FlagCategory,
2070 : VisBufferComponent2::FlagCube,
2071 : VisBufferComponent2::FlagRow,
2072 : VisBufferComponent2::ObservationId,
2073 : VisBufferComponent2::ProcessorId,
2074 : VisBufferComponent2::RowIds,
2075 : VisBufferComponent2::Scan,
2076 : VisBufferComponent2::Sigma,
2077 : VisBufferComponent2::SpectralWindows,
2078 : VisBufferComponent2::StateId,
2079 : VisBufferComponent2::Time,
2080 : VisBufferComponent2::TimeCentroid,
2081 : VisBufferComponent2::TimeInterval,
2082 : VisBufferComponent2::Weight,
2083 0 : VisBufferComponent2::Uvw});
2084 :
2085 0 : if (doing_p.modelData_p){
2086 0 : including += VisBufferComponent2::VisibilityCubeModel;
2087 : }
2088 :
2089 0 : if (doing_p.correctedData_p){
2090 0 : including += VisBufferComponent2::VisibilityCubeCorrected;
2091 : }
2092 :
2093 0 : if (doing_p.observedData_p){
2094 0 : including += VisBufferComponent2::VisibilityCubeObserved;
2095 : }
2096 :
2097 0 : if (doing_p.floatData_p){
2098 0 : including += VisBufferComponent2::VisibilityCubeFloat;
2099 : }
2100 :
2101 0 : if (doing_p.weightSpectrumOut_p){
2102 0 : including += VisBufferComponent2::WeightSpectrum;
2103 : }
2104 :
2105 0 : if (doing_p.sigmaSpectrumOut_p){
2106 0 : including += VisBufferComponent2::SigmaSpectrum;
2107 : }
2108 :
2109 0 : cacheResizeAndZero (including);
2110 :
2111 0 : correlationFlags_p.reference (Matrix<Bool> (IPosition (2, nCorrelations, nBaselines), true));
2112 0 : counts_p.reference (Cube<Int> (IPosition (3, nCorrelations, nChannels, nBaselines), 0));
2113 0 : countsBaseline_p.reference (Vector<Int> (nBaselines, 0)); // number of items summed together for each baseline.
2114 0 : intervalLast_p.reference (Vector<Double> (nBaselines, 0));
2115 0 : timeFirst_p.reference (Vector<Double> (nBaselines, 0));
2116 0 : timeLast_p.reference (Vector<Double> (nBaselines, 0));
2117 0 : uvwFirst_p.reference (Matrix<Double> (IPosition (2, 3, nBaselines), 0.0));
2118 0 : }
2119 :
2120 : void
2121 0 : VbAvg::startChunk (ViImplementation2 * vi)
2122 : {
2123 0 : empty_p = true;
2124 :
2125 0 : rowIdGenerator_p = 0;
2126 0 : row2AvgRow_p.resize(0);
2127 :
2128 : // See if the new MS has corrected and/or model data columns
2129 :
2130 0 : doing_p.observedData_p = averagingOptions_p.contains (AveragingOptions::AverageObserved);
2131 0 : doing_p.correctedData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeCorrected) &&
2132 0 : averagingOptions_p.contains (AveragingOptions::AverageCorrected);
2133 0 : doing_p.modelData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeModel) &&
2134 0 : averagingOptions_p.contains (AveragingOptions::AverageModel);
2135 0 : doing_p.floatData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeFloat) &&
2136 0 : averagingOptions_p.contains (AveragingOptions::AverageFloat);
2137 :
2138 0 : doing_p.weightSpectrumIn_p = doing_p.correctedData_p;
2139 0 : doing_p.sigmaSpectrumIn_p = doing_p.observedData_p || doing_p.floatData_p;
2140 0 : doing_p.weightSpectrumOut_p = true; // We always use the output WeightSpectrum
2141 0 : doing_p.sigmaSpectrumOut_p = true; // We always use the output SigmaSpectrum
2142 :
2143 0 : if (doing_p.observedData_p or doing_p.correctedData_p or doing_p.modelData_p or doing_p.floatData_p)
2144 : {
2145 0 : doing_p.onlymetadata_p = false;
2146 : }
2147 :
2148 : // Set up the flags for row copying
2149 :
2150 0 : optionalComponentsToCopy_p = VisBufferComponents2::none();
2151 :
2152 0 : if (doing_p.observedData_p){
2153 0 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeObserved;
2154 : }
2155 :
2156 0 : if (doing_p.correctedData_p){
2157 0 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeCorrected;
2158 : }
2159 :
2160 0 : if (doing_p.modelData_p){
2161 0 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeModel;
2162 : }
2163 :
2164 0 : if (doing_p.floatData_p){
2165 0 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeFloat;
2166 : }
2167 :
2168 0 : if (doing_p.weightSpectrumOut_p){
2169 0 : optionalComponentsToCopy_p += VisBufferComponent2::WeightSpectrum;
2170 : }
2171 :
2172 0 : if (doing_p.sigmaSpectrumOut_p){
2173 0 : optionalComponentsToCopy_p += VisBufferComponent2::SigmaSpectrum;
2174 : }
2175 0 : }
2176 :
2177 : void
2178 0 : VbAvg::transferBaseline (MsRowAvg * rowAveraged)
2179 : {
2180 0 : rowAveraged->setRowId (rowIdGenerator_p ++);
2181 0 : bufferToFill_p->appendRow (rowAveraged, nBaselines(), optionalComponentsToCopy_p);
2182 :
2183 0 : rowAveraged->setBaselinePresent(false);
2184 0 : }
2185 :
2186 :
2187 : //VbSet::VbSet (const AveragingParameters & averagingParameters)
2188 : //: averagingInterval_p (averagingParameters.getAveragingInterval ()),
2189 : // averagingOptions_p (averagingParameters.getOptions()),
2190 : // averagingParameters_p (averagingParameters),
2191 : // doingCorrectedData_p (false),
2192 : // doingModelData_p (false),
2193 : // doingObservedData_p (false),
2194 : // doingWeightSpectrum_p (false),
2195 : // doingsigmaSpectrum_p (false),
2196 : // vbAveragers_p ()
2197 : //{}
2198 :
2199 : //VbSet::~VbSet ()
2200 : //{
2201 : // flush (true); // allow killing nonempty buffers
2202 : //}
2203 : //
2204 : //void
2205 : //VbSet::accumulate (const VisBuffer2 * input, const Subchunk & subchunk)
2206 : //{
2207 : // Int ddId = input->dataDescriptionIds()(0);
2208 : //
2209 : // if (! utilj::containsKey (ddId, vbAveragers_p)){ // Need a new averager
2210 : // add (ddId);
2211 : // }
2212 : //
2213 : // VbAvg * vba = vbAveragers_p [ddId];
2214 : //
2215 : // vba->accumulate (input, subchunk);
2216 : //}
2217 : //
2218 : //VbAvg *
2219 : //VbSet::add (Int ddId)
2220 : //{
2221 : // VbAvg::Doing doing;
2222 : //
2223 : // doing.correctedData_p = doingCorrectedData_p;
2224 : // doing.modelData_p = doingModelData_p;
2225 : // doing.observedData_p = doingObservedData_p;
2226 : // doing.weightSpectrum_p = doingWeightSpectrum_p;
2227 : // doing.sigmaSpectrum_p = doingsigmaSpectrum_p;
2228 : //
2229 : // VbAvg * newAverager = new VbAvg (doing, averagingParameters_p);
2230 : //
2231 : // vbAveragers_p [ddId] = newAverager;
2232 : //
2233 : // return newAverager;
2234 : //}
2235 : //
2236 : //Bool
2237 : //VbSet::anyAveragesReady(Int ddid) const
2238 : //{
2239 : // Bool any = false;
2240 : //
2241 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2242 : // a != vbAveragers_p.end();
2243 : // a++){
2244 : //
2245 : // if (a->second->isComplete() &&
2246 : // (ddid < 0 || ddid == a->second->dataDescriptionIds()(0))){
2247 : //
2248 : // any = true;
2249 : // break;
2250 : // }
2251 : // }
2252 : //
2253 : // return any;
2254 : //}
2255 : //
2256 : //Bool
2257 : //VbSet::contains (Int ddId) const
2258 : //{
2259 : // return vbAveragers_p.find (ddId) != vbAveragers_p.end();
2260 : //}
2261 : //
2262 : ////void
2263 : ////VbSet::finalizeAverage (Int ddId)
2264 : ////{
2265 : //// vbAveragers_p [ddId]->finalizeAverage();
2266 : ////}
2267 : //
2268 : //void
2269 : //VbSet::finalizeAverages ()
2270 : //{
2271 : //// for (Averagers::iterator a = vbAveragers_p.begin();
2272 : //// a != vbAveragers_p.end();
2273 : //// a ++){
2274 : ////
2275 : //// finalizeAverage (a->first);
2276 : //// }
2277 : //}
2278 : //
2279 : //void
2280 : //VbSet::flush (Bool okIfNonempty, ViImplementation2 * vi)
2281 : //{
2282 : // // Get rid of all of the averagers. This is done at
2283 : // // destruction time and when a sweeping into a new MS.
2284 : //
2285 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2286 : // a != vbAveragers_p.end();
2287 : // a ++){
2288 : //
2289 : // Assert (okIfNonempty || (a->second)->empty());
2290 : // // should have been emptied before calling this
2291 : //
2292 : // delete a->second;
2293 : // }
2294 : //
2295 : // vbAveragers_p.clear();
2296 : //
2297 : // seeIfCubeColumnsExist (vi);
2298 : //}
2299 : //
2300 : //void
2301 : //VbSet::seeIfCubeColumnsExist (ViImplementation2 * vi)
2302 : //{
2303 : // if (vi != 0){
2304 : //
2305 : // // See if the new MS has corrected and/or model data columns
2306 : //
2307 : // doingObservedData_p = averagingOptions_p.contains (AveragingOptions::AverageObserved);
2308 : // doingCorrectedData_p = vi->existsColumn (VisibilityCubeCorrected) &&
2309 : // averagingOptions_p.contains (AveragingOptions::AverageCorrected);
2310 : // doingModelData_p = vi->existsColumn (VisibilityCubeModel) &&
2311 : // averagingOptions_p.contains (AveragingOptions::AverageModel);
2312 : // doingWeightSpectrum_p = vi->existsColumn (WeightSpectrum);
2313 : //
2314 : // // If the use of corrected weights were specified for one of the averages, it's an
2315 : // // error if the column does not exist. Also set the doing flag for corrected weights
2316 : // // if it's being used in some way.
2317 : //
2318 : // Bool needCorrectedWeights =
2319 : // averagingOptions_p.contains (AveragingOptions::ModelUseCorrectedWeights) ||
2320 : // averagingOptions_p.contains (AveragingOptions::CorrectedUseCorrectedWeights);
2321 : //
2322 : // Bool correctedWeightsMissing = needCorrectedWeights &&
2323 : // ! vi->existsColumn (sigmaSpectrum);
2324 : //
2325 : // ThrowIf (correctedWeightsMissing,
2326 : // "Corrected_weight_spectrum not present but required by provided averaging options");
2327 : //
2328 : // doingsigmaSpectrum_p = needCorrectedWeights;
2329 : // }
2330 : //}
2331 : //
2332 : //Int
2333 : //VbSet::getFirstReadyDdid () const
2334 : //{
2335 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2336 : // a != vbAveragers_p.end();
2337 : // a ++){
2338 : //
2339 : // if (a->second->isComplete()){
2340 : // return a->first;
2341 : // }
2342 : // }
2343 : //
2344 : // Assert (false); // ought to have been one that's ready
2345 : //
2346 : // return -1; // shouldn't be called
2347 : //}
2348 : //
2349 : ////void
2350 : ////VbSet::transferAverage (Int ddId, VisBuffer2 * vb)
2351 : ////{
2352 : //// Assert (utilj::containsKey (ddId, vbAveragers_p));
2353 : ////
2354 : //// VbAvg * vba = vbAveragers_p [ddId];
2355 : ////
2356 : //// Assert (vba != 0 && ! vba->empty ());
2357 : ////
2358 : //// // Copy< the completed average into the provided VisBuffer, but
2359 : //// // first reshape the VB if it's shape is different.
2360 : ////
2361 : //// vba->transferAverage (vb);
2362 : ////
2363 : ////}
2364 : //
2365 : //void
2366 : //VbSet::zero ()
2367 : //{
2368 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2369 : // a != vbAveragers_p.end();
2370 : // a ++){
2371 : //
2372 : // a->second->markEmpty();
2373 : // }
2374 : //}
2375 :
2376 : ///////////////////////
2377 : // //
2378 : // End Namespace avg //
2379 : // //
2380 : ///////////////////////
2381 :
2382 : } // end avg
2383 :
2384 : using namespace avg;
2385 :
2386 0 : AveragingTvi2::AveragingTvi2 (ViImplementation2 * inputVi,
2387 0 : const AveragingParameters & averagingParameters)
2388 : : TransformingVi2 (inputVi),
2389 0 : averagingInterval_p (averagingParameters.getAveragingInterval()),
2390 0 : averagingOptions_p (averagingParameters.getOptions()),
2391 : averagingParameters_p (averagingParameters),
2392 : ddidLastUsed_p (-1),
2393 : inputViiAdvanced_p (false),
2394 0 : vbAvg_p (new VbAvg (averagingParameters, this))
2395 : {
2396 0 : validateInputVi (inputVi);
2397 :
2398 : // Position input Vi to the first subchunk
2399 :
2400 0 : getVii()->originChunks();
2401 0 : getVii()->origin ();
2402 :
2403 0 : setVisBuffer (createAttachedVisBuffer (VbNoOptions));
2404 0 : }
2405 :
2406 0 : AveragingTvi2::~AveragingTvi2 ()
2407 : {
2408 0 : delete vbAvg_p;
2409 0 : }
2410 :
2411 : //void
2412 : //AveragingTvi2::advanceInputVii ()
2413 : //{
2414 : // Assert (false);
2415 : //
2416 : // // Attempt to advance to the next subchunk
2417 : //
2418 : // getVii()->next ();
2419 : //
2420 : // if (! getVii()->more()){
2421 : //
2422 : // // No more subchunks so advance to the next chunk
2423 : //
2424 : // getVii()->nextChunk();
2425 : //
2426 : // if (! getVii()->moreChunks()){
2427 : // return; // no more chunks
2428 : // }
2429 : //
2430 : // // Position to the first subchunk
2431 : //
2432 : // getVii()->origin();
2433 : // }
2434 : //}
2435 :
2436 : //Int
2437 : //AveragingTvi2::determineDdidToUse() const
2438 : //{
2439 : // if (ddidLastUsed_p >= 0 && vbSet_p->anyAveragesReady (ddidLastUsed_p)){
2440 : // return ddidLastUsed_p;
2441 : // }
2442 : //
2443 : // return vbSet_p->getFirstReadyDdid();
2444 : //}
2445 :
2446 : Bool
2447 0 : AveragingTvi2::more () const
2448 : {
2449 0 : return more_p;
2450 : }
2451 :
2452 : Bool
2453 0 : AveragingTvi2::moreChunks () const
2454 : {
2455 0 : return getVii()->moreChunks();
2456 : }
2457 :
2458 : void
2459 0 : AveragingTvi2::next ()
2460 : {
2461 0 : subchunkExists_p = false;
2462 :
2463 0 : startBuffer_p = endBuffer_p + 1;
2464 0 : endBuffer_p = startBuffer_p - 1;
2465 :
2466 0 : if (getVii()->more()){
2467 0 : getVii()->next();
2468 0 : endBuffer_p++;
2469 : }
2470 :
2471 0 : produceSubchunk ();
2472 :
2473 0 : subchunk_p.incrementSubChunk();
2474 0 : }
2475 :
2476 : void
2477 0 : AveragingTvi2::nextChunk ()
2478 : {
2479 : // New chunk, so toss all of the accumulators
2480 :
2481 0 : vbAvg_p->startChunk (getVii());
2482 :
2483 : // Advance the input to the next chunk as well.
2484 :
2485 0 : getVii()->nextChunk ();
2486 :
2487 0 : subchunk_p.incrementChunk();
2488 :
2489 0 : more_p = false;
2490 0 : }
2491 :
2492 : void
2493 0 : AveragingTvi2::origin ()
2494 : {
2495 : // Position input VI to the start of the chunk
2496 :
2497 0 : subchunk_p.resetSubChunk();
2498 :
2499 0 : getVii()->origin();
2500 :
2501 0 : startBuffer_p = 0;
2502 0 : endBuffer_p = -1;
2503 :
2504 : // Get the first subchunk ready.
2505 :
2506 0 : produceSubchunk ();
2507 0 : }
2508 :
2509 : void
2510 0 : AveragingTvi2::originChunks (Bool forceRewind)
2511 : {
2512 : // Ensure that the underlying VI is in a state where some metadata
2513 : // can be retrieved
2514 0 : getVii()->originChunks(forceRewind);
2515 :
2516 : // Initialize the chunk
2517 0 : vbAvg_p->startChunk (getVii());
2518 :
2519 0 : more_p = false;
2520 :
2521 0 : subchunk_p.resetToOrigin();
2522 0 : }
2523 :
2524 : /**
2525 : * Configure a VisBuffer with given averagingOptions related to phase shifting
2526 : *
2527 : * @param vb a VisBuffer to set up in terms of phase shifting
2528 : * @param averagingOptions averaging options enabled
2529 : * @param avgPars AveragingParmeters object to set into the buffer
2530 : */
2531 : void
2532 0 : setPhaseShiftingOptions(VisBuffer2 * vb, const AveragingOptions &averagingOptions,
2533 : const AveragingParameters &avgPars)
2534 : {
2535 0 : if (averagingOptions.contains (AveragingOptions::phaseShifting))
2536 : {
2537 0 : if (averagingOptions.contains (AveragingOptions::AverageObserved))
2538 : {
2539 0 : vb->visCube();
2540 : }
2541 :
2542 0 : if (averagingOptions.contains (AveragingOptions::AverageCorrected))
2543 : {
2544 0 : vb->visCubeCorrected();
2545 : }
2546 :
2547 0 : if (averagingOptions.contains (AveragingOptions::AverageModel))
2548 : {
2549 0 : vb->visCubeModel();
2550 : }
2551 :
2552 0 : vb->phaseCenterShift(avgPars.getXpcOffset(),
2553 0 : avgPars.getYpcOffset());
2554 : }
2555 0 : }
2556 :
2557 : /**
2558 : * The iteration of this method is different depending on whether "row blocking" is used or
2559 : * not. The reason is that the while loop already had enough complexity embedded when fixes
2560 : * were done to get flagdata+time-average+row-blocking working (CAS-11910). Hopefully in the
2561 : * near future we can get rid of the hacky "row blocking" feature. For the time being, it is
2562 : * not clear how it could possibly work together with the "uvwdistance" feature. So better
2563 : * to keep those features separate.
2564 : *
2565 : * So the "if (block > 0)" separates iteration when using row blocking. That implies that
2566 : * row blocking takes precedence over (and disables) other features like
2567 : * "isUsingUvwDistance()".
2568 : * An alternative would be to add comparisons between block and vbToFill->appendSize() in
2569 : * the ifs below. Something like:
2570 : * if (! vbAvg_p->isUsingUvwDistance()
2571 : * && (block == 0 && vbToFill->appendSize() > 0
2572 : * || (block > 0 && vbToFill->appendSize() >= block)
2573 : * )
2574 : * ){
2575 : * ...
2576 : * else if ((block > 0 && vbToFill->appendSize() < block) ||
2577 : * vbToFill->appendSize() < nBaselines * nWindows){
2578 : * ...
2579 : * } else {
2580 : *
2581 : * But I prefer not adding more complexity to those ifs. The potential combinations would
2582 : * be too many to handle in a handful of if branches, and they are not well understood let
2583 : * alone well tested.
2584 : */
2585 : void
2586 0 : AveragingTvi2::produceSubchunk ()
2587 : {
2588 0 : VisBufferImpl2 * vbToFill = dynamic_cast<VisBufferImpl2 *> (getVisBuffer());
2589 0 : assert (vbToFill != 0);
2590 :
2591 0 : vbToFill->setFillable (true); // New subchunk, so it's fillable
2592 :
2593 0 : vbAvg_p->setBufferToFill (vbToFill);
2594 :
2595 0 : Int nBaselines = nAntennas() * (nAntennas() -1) / 2;
2596 : // This is just a heuristic to keep output VBs from being too small
2597 :
2598 : // jagonzal: Handle nBaselines for SD case
2599 0 : if (nBaselines == 0) nBaselines = 1;
2600 :
2601 0 : auto block = getVii()->getRowBlocking();
2602 0 : while (getVii()->more()){
2603 :
2604 0 : VisBuffer2 * vb = getVii()->getVisBuffer();
2605 :
2606 0 : setPhaseShiftingOptions(vb, averagingOptions_p, averagingParameters_p);
2607 :
2608 0 : bool rowBlocking = block > 0;
2609 0 : vbAvg_p->accumulate (vb, subchunk_p);
2610 :
2611 0 : if (rowBlocking) {
2612 0 : auto app = vbToFill->appendSize();
2613 0 : if (app <= block) {
2614 0 : getVii()->next();
2615 0 : endBuffer_p++;
2616 : } else {
2617 0 : break;
2618 : }
2619 : } else {
2620 0 : Int nWindows = vbAvg_p->nSpectralWindowsInBuffer ();
2621 0 : if (! vbAvg_p->isUsingUvwDistance() && vbToFill->appendSize() > 0){
2622 : // Doing straight average and some data has been produced so
2623 : // output it to the user
2624 0 : break;
2625 : }
2626 0 : else if (vbToFill->appendSize() < nBaselines * nWindows) {
2627 0 : getVii()->next();
2628 0 : endBuffer_p++;
2629 : }
2630 : else {
2631 0 : break;
2632 : }
2633 : }
2634 : };
2635 :
2636 0 : if (! getVii()->more()){
2637 0 : vbAvg_p->finalizeAverages ();
2638 : }
2639 :
2640 0 : vbAvg_p->finalizeBufferFilling ();
2641 :
2642 0 : more_p = getVii()->more() || // more to read
2643 0 : vbToFill->nRows() > 0; // some to process
2644 0 : }
2645 :
2646 : // Bool
2647 : // AveragingTvi2::reachedAveragingBoundary()
2648 : // {
2649 : // // An average can be terminated for a variety of reasons:
2650 : // // o the time interval has elapsed
2651 : // // o the current MS is completed
2652 : // // o no more input data
2653 : // // o other (e.g, change of scan, etc.)
2654 :
2655 : // Bool reachedIt = false;
2656 : // VisBuffer2 * vb = getVii()->getVisBuffer();
2657 :
2658 : // if (! getVii()->more() && ! getVii ()->moreChunks()){
2659 :
2660 : // reachedIt = true; // no more input data
2661 : // }
2662 : // else if (vb->isNewMs()){
2663 : // reachedIt = true; // new MS
2664 : // }
2665 :
2666 : // return reachedIt;
2667 : // }
2668 :
2669 : //Bool
2670 : //AveragingTvi2::subchunksReady() const
2671 : //{
2672 : // Bool ready = vbSet_p->anyAveragesReady();
2673 : //
2674 : // return ready;
2675 : //}
2676 :
2677 : void
2678 0 : AveragingTvi2::validateInputVi (ViImplementation2 *)
2679 : {
2680 : // Validate that the input VI is compatible with this VI.
2681 :
2682 : // No implmented right now
2683 0 : }
2684 :
2685 0 : Float AveragingTvi2::weightToSigma (Float weight)
2686 : {
2687 0 : return weight > FLT_MIN ? 1.0 / std::sqrt (weight) : -1.0; // bogosity indicator
2688 : }
2689 :
2690 : Vector<Float>
2691 0 : AveragingTvi2::average (const Matrix<Float> &data, const Matrix<Bool> &flags)
2692 : {
2693 0 : IPosition shape = data.shape();
2694 0 : Vector<Float> result(shape(0),0);
2695 0 : Vector<uInt> samples(shape(0),0);
2696 0 : uInt nCorrelations = shape (0);
2697 0 : uInt nChannels = shape (1);
2698 :
2699 0 : for (uInt correlation = 0; correlation < nCorrelations; correlation++)
2700 : {
2701 0 : int nSamples = 0;
2702 0 : float sum = 0;
2703 0 : bool accumulatorFlag = true;
2704 :
2705 0 : for (uInt channel=0; channel< nChannels; channel++)
2706 : {
2707 0 : Bool inputFlag = flags(correlation,channel);
2708 : // true/true or false/false
2709 0 : if (accumulatorFlag == inputFlag)
2710 : {
2711 0 : nSamples ++;
2712 0 : sum += data (correlation, channel);
2713 : }
2714 : // true/false: Reset accumulation when accumulator switches from flagged to unflagged
2715 0 : else if ( accumulatorFlag and ! inputFlag )
2716 : {
2717 0 : accumulatorFlag = false;
2718 0 : nSamples = 1;
2719 0 : sum = data (correlation, channel);
2720 : }
2721 : // else ignore case where accumulator is valid and data is not
2722 :
2723 : }
2724 :
2725 0 : result (correlation) = sum / (nSamples != 0 ? nSamples : 1);
2726 : }
2727 :
2728 0 : return result;
2729 : }
2730 :
2731 : Matrix<Float>
2732 0 : AveragingTvi2::average (const Cube<Float> &data, const Cube<Bool> &flags)
2733 : {
2734 0 : IPosition shape = data.shape();
2735 0 : uInt nRows = shape(2);
2736 0 : uInt nCorrelations = shape (0);
2737 :
2738 0 : Matrix<Float> result(nCorrelations, nRows, 0);
2739 :
2740 0 : for (uInt row=0; row < nRows; row++)
2741 : {
2742 0 : result.column(row) = AveragingTvi2::average (data.xyPlane(row), flags.xyPlane(row));
2743 : }
2744 :
2745 0 : return result;
2746 : }
2747 :
2748 : /**
2749 : * Strategy to support different ways of propagating flags from the 'transformed' cube to
2750 : * the original ('propagated') cube. Iterates through rows, channels, correlations.
2751 : *
2752 : * This is meant to be used from propagateChanAvgFlags with at least two alternative
2753 : * functors. One to propagate flags as required by flagdata (preserve pre-existing flags
2754 : * in the original data cube), and a second one to propagate flags as required by plotms.
2755 : * CAS-12737, CAS-9985, CAS-12205.
2756 : *
2757 : * @param flagrow per row FLAG_ROW value
2758 : * @param flagMapped propagated FLAG_ROW
2759 : * @param flagCubeMapped Cube of flags in which flags are to be propagated
2760 : * @param row2AvgRow map of input_buffer_row_index->output_buffer_row_index (this is pre-
2761 : * calculated by the "VbAvg" when averaging rows and is needed here).
2762 : */
2763 : template <typename Functor>
2764 0 : void cubePropagateFlags(const Vector<Bool> &flagRow,
2765 : Vector<Bool> &flagMapped,
2766 : Cube<Bool> &flagCubeMapped,
2767 : std::vector<size_t> row2AvgRow,
2768 : Functor propagate)
2769 : {
2770 0 : Int nRowsMapped = flagCubeMapped.shape()[2];
2771 0 : for(Int rowMapped=0; rowMapped < nRowsMapped; ++rowMapped) {
2772 0 : size_t index = row2AvgRow[rowMapped];
2773 0 : flagMapped(rowMapped) = flagRow(index);
2774 :
2775 0 : for (Int chan_i=0; chan_i < flagCubeMapped.shape()[1]; ++chan_i) {
2776 0 : for (Int corr_i=0; corr_i<flagCubeMapped.shape()[0]; ++corr_i) {
2777 0 : propagate(rowMapped, chan_i, corr_i, index);
2778 : }
2779 : }
2780 : }
2781 0 : }
2782 :
2783 : // -----------------------------------------------------------------------
2784 : //
2785 : // -----------------------------------------------------------------------
2786 0 : void AveragingTvi2::writeFlag (const Cube<Bool> & flag)
2787 : {
2788 : // Calculate FLAG_ROW from flag
2789 0 : Vector<Bool> flagRow;
2790 0 : TransformingVi2::calculateFlagRowFromFlagCube(flag, flagRow);
2791 :
2792 : const auto flagdataFlagPropagation =
2793 0 : averagingOptions_p.contains(AveragingOptions::flagdataFlagPropagation);
2794 :
2795 : // Propagate transformed flags
2796 0 : getVii()->origin();
2797 0 : Int currentBuffer = 0;
2798 0 : while (getVii()->more())
2799 : {
2800 0 : if ((currentBuffer >= startBuffer_p) and (currentBuffer <= endBuffer_p))
2801 : {
2802 : // Allocated propagated flag vector/cube
2803 0 : uInt nOriginalRows = getVii()->getVisBuffer()->nRows();
2804 0 : Vector<Bool> flagMapped(nOriginalRows,false);
2805 0 : Cube<Bool> flagCubeMapped;
2806 0 : flagCubeMapped = getVii()->getVisBuffer()->flagCube();
2807 :
2808 : // Keeping two separate blocks for 'flagdataFlagPropagation' (CAS-12737,
2809 : // CAS-12205, CAS-9985) until this issue is better settled.
2810 : // Fill propagated flag vector/cube
2811 0 : if (flagdataFlagPropagation)
2812 : {
2813 0 : cubePropagateFlags(flagRow, flagMapped, flagCubeMapped, vbAvg_p->row2AvgRow(),
2814 : [&flag, &flagCubeMapped]
2815 0 : (uInt rowMapped, uInt chan_i, uInt corr_i, uInt index) {
2816 0 : if (flag(corr_i,chan_i,index))
2817 0 : flagCubeMapped(corr_i,chan_i,rowMapped) = true;
2818 0 : });
2819 : }
2820 : else
2821 : {
2822 0 : cubePropagateFlags(flagRow, flagMapped, flagCubeMapped, vbAvg_p->row2AvgRow(),
2823 : [&flag, &flagCubeMapped]
2824 0 : (uInt rowMapped, uInt chan_i, uInt corr_i, uInt index) {
2825 0 : flagCubeMapped(corr_i, chan_i, rowMapped) =
2826 0 : flag(corr_i, chan_i, index);
2827 0 : });
2828 : }
2829 :
2830 : // Write propagated flag vector/cube
2831 0 : getVii()->writeFlag(flagCubeMapped);
2832 0 : getVii()->writeFlagRow(flagMapped);
2833 : }
2834 :
2835 0 : currentBuffer++;
2836 0 : getVii()->next();
2837 0 : if (currentBuffer > endBuffer_p) break;
2838 : }
2839 0 : }
2840 :
2841 : // -----------------------------------------------------------------------
2842 : //
2843 : // -----------------------------------------------------------------------
2844 0 : void AveragingTvi2::writeFlagRow (const Vector<Bool> & flagRow)
2845 : {
2846 : // Create index map for averaged data
2847 0 : VisBuffer2 *avgVB = getVisBuffer();
2848 0 : Vector<Int> avgAnt1 = avgVB->antenna1();
2849 0 : Vector<Int> avgAnt2 = avgVB->antenna2();
2850 0 : Vector<Int> avgSPW = avgVB->spectralWindows();
2851 :
2852 0 : std::map< Int, std::map <Int, std::map< Int, uInt> > > spwAnt1Ant2IndexMap;
2853 0 : for (uInt avgRow=0;avgRow<avgAnt1.size();avgRow++)
2854 : {
2855 0 : spwAnt1Ant2IndexMap[avgSPW(avgRow)][avgAnt1(avgRow)][avgAnt2(avgRow)] = avgRow;
2856 : }
2857 :
2858 : // Propagate transformed flags
2859 0 : getVii()->origin();
2860 0 : Int currentBuffer = 0;
2861 0 : while (getVii()->more())
2862 : {
2863 0 : if ((currentBuffer >= startBuffer_p) and (currentBuffer <= endBuffer_p))
2864 : {
2865 : // Allocated propagated flag vector/cube
2866 0 : uInt nOriginalRows = getVii()->getVisBuffer()->nRows();
2867 0 : Vector<Bool> flagMapped(nOriginalRows,false);
2868 :
2869 : // Get original ant1/ant2/spw cols. to determine the mapped index
2870 0 : Vector<Int> orgAnt1 = getVii()->getVisBuffer()->antenna1();
2871 0 : Vector<Int> orgAnt2 = getVii()->getVisBuffer()->antenna2();
2872 0 : Vector<Int> orgSPW = getVii()->getVisBuffer()->spectralWindows();
2873 :
2874 : // Fill propagated flag vector/cube
2875 0 : for (uInt row=0;row<nOriginalRows;row++)
2876 : {
2877 0 : uInt index = spwAnt1Ant2IndexMap[orgSPW(row)][orgAnt1(row)][orgAnt2(row)];
2878 0 : flagMapped(row) = flagRow(index);
2879 : }
2880 :
2881 : // Write propagated flag vector/cube
2882 0 : getVii()->writeFlagRow(flagMapped);
2883 :
2884 : }
2885 :
2886 0 : currentBuffer += 1;
2887 0 : getVii()->next();
2888 0 : if (currentBuffer > endBuffer_p) break;
2889 : }
2890 :
2891 0 : return;
2892 : }
2893 :
2894 0 : void AveragingTvi2::visibilityObserved(casacore::Cube<casacore::Complex>& vis) const
2895 : {
2896 0 : if(!averagingOptions_p.contains(AveragingOptions::AverageObserved))
2897 0 : throw AipsError("Requesting visibilityCorrected but column has not been averaged");
2898 0 : VisBuffer2* vb = getVisBuffer();
2899 0 : vis = vb->visCube();
2900 0 : return;
2901 : }
2902 :
2903 0 : void AveragingTvi2::visibilityCorrected(casacore::Cube<casacore::Complex>& vis) const
2904 : {
2905 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeCorrected) ||
2906 0 : !averagingOptions_p.contains(AveragingOptions::AverageCorrected))
2907 : throw AipsError("Requesting visibilityCorrected but column not "
2908 0 : "provided by underlying VI or column has not been averaged");
2909 0 : VisBuffer2* vb = getVisBuffer();
2910 0 : vis = vb->visCubeCorrected();
2911 0 : return;
2912 : }
2913 :
2914 0 : void AveragingTvi2::visibilityModel(casacore::Cube<casacore::Complex>& vis) const
2915 : {
2916 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeModel) ||
2917 0 : !averagingOptions_p.contains(AveragingOptions::AverageModel))
2918 : throw AipsError("Requesting visibilityModel but column not "
2919 0 : "provided by underlying VI or column has not been averaged");
2920 0 : VisBuffer2* vb = getVisBuffer();
2921 0 : vis = vb->visCubeModel();
2922 0 : return;
2923 : }
2924 :
2925 0 : void AveragingTvi2::floatData(casacore::Cube<casacore::Float>& fcube) const
2926 : {
2927 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeFloat) ||
2928 0 : !averagingOptions_p.contains(AveragingOptions::AverageFloat))
2929 : throw AipsError("Requesting floatData but column not "
2930 0 : "provided by underlying VI or column has not been averaged");
2931 0 : VisBuffer2* vb = getVisBuffer();
2932 0 : fcube = vb->visCubeFloat();
2933 0 : return;
2934 : }
2935 :
2936 0 : void AveragingTvi2::flag(casacore::Cube<casacore::Bool>& flags) const
2937 : {
2938 0 : VisBuffer2* vb = getVisBuffer();
2939 0 : flags = vb->flagCube();
2940 0 : return;
2941 : }
2942 :
2943 0 : void AveragingTvi2::flagRow(casacore::Vector<casacore::Bool>& rowflags) const
2944 : {
2945 0 : VisBuffer2* vb = getVisBuffer();
2946 0 : rowflags = vb->flagRow();
2947 0 : return;
2948 : }
2949 :
2950 0 : void AveragingTvi2::sigma(casacore::Matrix<casacore::Float> & sigmat) const
2951 : {
2952 0 : VisBuffer2* vb = getVisBuffer();
2953 0 : sigmat = vb->sigma();
2954 0 : return;
2955 : }
2956 :
2957 0 : void AveragingTvi2::weight (casacore::Matrix<casacore::Float> & wtmat) const
2958 : {
2959 0 : VisBuffer2* vb = getVisBuffer();
2960 0 : wtmat = vb->weight();
2961 0 : return;
2962 : }
2963 :
2964 0 : void AveragingTvi2::weightSpectrum (casacore::Cube<casacore::Float> & wtsp) const
2965 : {
2966 0 : VisBuffer2* vb = getVisBuffer();
2967 0 : wtsp = vb->weightSpectrum();
2968 0 : return;
2969 : }
2970 :
2971 0 : void AveragingTvi2::sigmaSpectrum (casacore::Cube<casacore::Float> & sigsp) const
2972 : {
2973 0 : VisBuffer2* vb = getVisBuffer();
2974 0 : sigsp = vb->sigmaSpectrum();
2975 0 : return;
2976 : }
2977 :
2978 0 : void AveragingTvi2::exposure (casacore::Vector<double> & expo) const
2979 : {
2980 0 : VisBuffer2* vb = getVisBuffer();
2981 0 : expo = vb->exposure();
2982 0 : return;
2983 : }
2984 :
2985 0 : void AveragingTvi2::getRowIds (Vector<rownr_t> & rowids) const
2986 : {
2987 0 : VisBuffer2* vb = getVisBuffer();
2988 0 : rowids = vb->rowIds();
2989 0 : return;
2990 : }
2991 :
2992 0 : void AveragingTvi2::time (casacore::Vector<double> & t) const
2993 : {
2994 0 : VisBuffer2* vb = getVisBuffer();
2995 0 : t = vb->time();
2996 0 : return;
2997 : }
2998 :
2999 0 : void AveragingTvi2::timeInterval (casacore::Vector<double> & ti) const
3000 : {
3001 0 : VisBuffer2* vb = getVisBuffer();
3002 0 : ti = vb->timeInterval();
3003 0 : return;
3004 : }
3005 :
3006 0 : void AveragingTvi2::timeCentroid (casacore::Vector<double> & t) const
3007 : {
3008 0 : VisBuffer2* vb = getVisBuffer();
3009 0 : t = vb->timeCentroid();
3010 0 : return;
3011 : }
3012 :
3013 0 : void AveragingTvi2::uvw (casacore::Matrix<double> & uvwmat) const
3014 : {
3015 0 : VisBuffer2* vb = getVisBuffer();
3016 0 : uvwmat = vb->uvw();
3017 0 : return;
3018 : }
3019 :
3020 0 : void AveragingTvi2::antenna1 (casacore::Vector<casacore::Int> & ant1) const
3021 : {
3022 0 : VisBuffer2* vb = getVisBuffer();
3023 0 : ant1 = vb->antenna1();
3024 0 : return;
3025 : }
3026 :
3027 0 : void AveragingTvi2::antenna2 (casacore::Vector<casacore::Int> & ant2) const
3028 : {
3029 0 : VisBuffer2* vb = getVisBuffer();
3030 0 : ant2 = vb->antenna2();
3031 0 : return;
3032 : }
3033 :
3034 0 : casacore::Bool AveragingTvi2::weightSpectrumExists () const
3035 : {
3036 : //According to VbAvg::startChunk code comments,
3037 : //there is always an output weightSpectrum. See also CAS-11559.
3038 0 : return true;
3039 : }
3040 :
3041 0 : casacore::Bool AveragingTvi2::sigmaSpectrumExists () const
3042 : {
3043 : //According to VbAvg::startChunk code comments,
3044 : //there is always an output sigmaSpectrum. See also CAS-11559.
3045 0 : return true;
3046 : }
3047 :
3048 : } // end namespace vi
3049 :
3050 : using namespace casacore;
3051 : } // end namespace casa
|
