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