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 14733 : SpwIndex (Int n) : Matrix<Int> (n, n, Empty) {}
125 :
126 : Int
127 11256450 : getBaselineNumber (Int antenna1, Int antenna2, Int & nBaselines)
128 : {
129 11256450 : Int i = (* this) (antenna1, antenna2);
130 :
131 11256450 : if (i == Empty){
132 :
133 278650 : i = nBaselines ++;
134 278650 : (* this) (antenna1, antenna2) = i;
135 : }
136 :
137 11256450 : 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 363 : BaselineIndex::BaselineIndex ()
159 363 : : nAntennas_p (0),
160 363 : nBaselines_p (0),
161 363 : nSpw_p (0)
162 363 : {}
163 :
164 363 : BaselineIndex::~BaselineIndex ()
165 : {
166 363 : destroy();
167 363 : }
168 :
169 : Int
170 11256450 : BaselineIndex::operator () (Int antenna1, Int antenna2, Int spectralWindow)
171 : {
172 11256450 : SpwIndex * spwIndex = getSpwIndex (spectralWindow);
173 11256450 : if (spwIndex == 0){
174 0 : addSpwIndex (spectralWindow);
175 : }
176 :
177 11256450 : Int i = spwIndex->getBaselineNumber (antenna1, antenna2, nBaselines_p);
178 :
179 11256450 : return i;
180 : }
181 :
182 :
183 :
184 : BaselineIndex::SpwIndex *
185 14733 : BaselineIndex::addSpwIndex (Int i)
186 : {
187 : // Delete an existing SpwIndex so that we start fresh
188 :
189 14733 : delete index_p [i];
190 :
191 : // Create a new SpwIndex and insert it into the main index.
192 :
193 14733 : index_p [i] = new SpwIndex (nAntennas_p);
194 :
195 14733 : return index_p [i];
196 : }
197 :
198 : void
199 14691 : BaselineIndex::configure (Int nAntennas, Int nSpw, const VisBuffer2 * vb)
200 : {
201 : // Capture the shape parameters
202 :
203 14691 : nAntennas_p = nAntennas;
204 14691 : nSpw_p = nSpw;
205 14691 : nBaselines_p = 0;
206 :
207 : // Get rid of the existing index
208 :
209 14691 : destroy ();
210 14691 : 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 313615 : for (rownr_t i = 0; i < vb->nRows(); i++){
218 :
219 : // Eagerly flesh out the Spw's index
220 :
221 298924 : Int spw = vb->spectralWindows() (i);
222 298924 : Int antenna1 = vb->antenna1()(i);
223 298924 : Int antenna2 = vb->antenna2()(i);
224 :
225 298924 : (* this) (antenna1, antenna2, spw);
226 : }
227 14691 : }
228 :
229 : void
230 15054 : 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 15054 : for (Index::iterator i = index_p.begin();
236 87635 : i != index_p.end();
237 72581 : i++){
238 :
239 72581 : delete (* i);
240 : }
241 15054 : }
242 :
243 : BaselineIndex::SpwIndex *
244 11256450 : BaselineIndex::getSpwIndex (Int spw)
245 : {
246 11256450 : AssertCc (spw < (int) index_p.size());
247 :
248 11256450 : SpwIndex * spwIndex = index_p [spw];
249 :
250 11256450 : if (spwIndex == 0){
251 14733 : spwIndex = addSpwIndex (spw);
252 : }
253 :
254 11256450 : return spwIndex;
255 : }
256 :
257 : template <typename T>
258 : class PrefilledMatrix {
259 :
260 : public:
261 :
262 1089 : PrefilledMatrix () : matrix_p (0, 0, 0), nChannels_p (0), nCorrelations_p (0) {}
263 :
264 : const Matrix<T> &
265 1960876 : getMatrix (Int nCorrelations, Int nChannels, const T & value)
266 : {
267 1960876 : if (nCorrelations != nCorrelations_p || nChannels != nChannels_p ||
268 1959706 : value != value_p){
269 :
270 1170 : nCorrelations_p = nCorrelations;
271 1170 : nChannels_p = nChannels;
272 1170 : value_p = value;
273 :
274 1170 : matrix_p.assign (Matrix<T> (nCorrelations_p, nChannels_p, value_p));
275 : }
276 :
277 1960876 : 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 607269 : CachedPlaneAvg (Accessor accessor) : accessor_p (accessor) {}
297 :
298 : Matrix<T> &
299 11476550 : getCachedPlane (casa::vi::avg::VbAvg * vb, Int row)
300 : {
301 11476550 : if (! isCached()){
302 :
303 : //cache_p.reference ((vb ->* accessor_p)().xyPlane (row)); // replace with something more efficient
304 10960983 : referenceMatrix (cache_p, (vb ->* accessor_p)(), row);
305 10960983 : setCached ();
306 : }
307 :
308 11476550 : return cache_p;
309 : }
310 :
311 : private:
312 :
313 : static void
314 10960983 : referenceMatrix (Matrix<T> & cache, const Cube<T> & src, Int row)
315 : {
316 10960983 : IPosition shape = src.shape ();
317 10960983 : 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 10960983 : T * storage = const_cast <T *> (& src (IPosition (3, 0, 0, row)));
323 :
324 10960983 : cache.takeStorage (shape, storage, casacore::SHARE);
325 10960983 : }
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 1214538 : 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 21915052 : 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 10957526 : void addInputRowIdx(Int idx) {
389 10957526 : inputRowIdxs_p[row()].push_back(idx);
390 10957526 : }
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 363 : Doing ()
455 363 : : correctedData_p (false),
456 363 : modelData_p (false),
457 363 : observedData_p (false),
458 363 : floatData_p(false),
459 363 : onlymetadata_p(true),
460 363 : weightSpectrumIn_p (false),
461 363 : sigmaSpectrumIn_p (false),
462 363 : weightSpectrumOut_p (false),
463 363 : sigmaSpectrumOut_p (false)
464 363 : {}
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 10957526 : AccumulationParameters (MsRow * rowInput, MsRowAvg * rowAveraged, const Doing & doing)
482 10957526 : : correctedIn_p (doing.correctedData_p ? rowInput->corrected().data() : 0),
483 10957526 : correctedOut_p (doing.correctedData_p ? rowAveraged->correctedMutable ().data(): 0),
484 10957526 : flagCubeIn_p (rowInput->flags().data()),
485 10957526 : flagCubeOut_p (rowAveraged->flagsMutable().data()),
486 10957526 : floatDataIn_p (doing.floatData_p ? rowInput->singleDishData().data() : 0),
487 10957526 : floatDataOut_p (doing.floatData_p ? rowAveraged->singleDishDataMutable().data() : 0),
488 10957526 : modelIn_p (doing.modelData_p ? rowInput->model().data(): 0),
489 10957526 : modelOut_p (doing.modelData_p ? rowAveraged->modelMutable().data() : 0),
490 10957526 : observedIn_p (doing.observedData_p ? rowInput->observed().data() : 0),
491 10957526 : observedOut_p (doing.observedData_p ? rowAveraged->observedMutable().data() : 0),
492 10957526 : sigmaIn_p (& rowInput->sigma()),
493 10957526 : sigmaOut_p (& rowAveraged->sigmaMutable()),
494 10957526 : sigmaSpectrumIn_p (doing.sigmaSpectrumIn_p ? rowInput->sigmaSpectrum().data() : 0),
495 10957526 : sigmaSpectrumOut_p (doing.sigmaSpectrumOut_p ? rowAveraged->sigmaSpectrumMutable().data() : 0),
496 10957526 : weightIn_p (& rowInput->weight()),
497 10957526 : weightOut_p (& rowAveraged->weightMutable()),
498 10957526 : weightSpectrumIn_p (doing.weightSpectrumIn_p ? rowInput->weightSpectrum().data() : 0),
499 10957526 : weightSpectrumOut_p (doing.weightSpectrumOut_p ? rowAveraged->weightSpectrumMutable().data() : 0)
500 10957526 : {}
501 :
502 779326154 : 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 779326154 : correctedIn_p && correctedIn_p ++;
508 779326154 : correctedOut_p && correctedOut_p ++;
509 779326154 : flagCubeIn_p && flagCubeIn_p ++;
510 779326154 : flagCubeOut_p && flagCubeOut_p ++;
511 779326154 : floatDataIn_p && floatDataIn_p ++;
512 779326154 : floatDataOut_p && floatDataOut_p ++;
513 779326154 : modelIn_p && modelIn_p ++;
514 779326154 : modelOut_p && modelOut_p ++;
515 779326154 : observedIn_p && observedIn_p ++;
516 779326154 : observedOut_p && observedOut_p ++;
517 779326154 : sigmaSpectrumIn_p && sigmaSpectrumIn_p ++;
518 779326154 : sigmaSpectrumOut_p && sigmaSpectrumOut_p ++;
519 779326154 : weightSpectrumIn_p && weightSpectrumIn_p ++;
520 779326154 : weightSpectrumOut_p && weightSpectrumOut_p ++;
521 779326154 : }
522 :
523 : inline const Complex *
524 239440858 : correctedIn ()
525 : {
526 239440858 : assert (correctedIn_p != 0);
527 239440858 : return correctedIn_p;
528 : }
529 :
530 : inline Complex *
531 239440858 : correctedOut ()
532 : {
533 239440858 : assert (correctedOut_p != 0);
534 239440858 : return correctedOut_p;
535 : }
536 :
537 : inline const Float *
538 192072192 : floatDataIn ()
539 : {
540 192072192 : return floatDataIn_p;
541 : }
542 :
543 : inline Float *
544 192072192 : floatDataOut ()
545 : {
546 192072192 : return floatDataOut_p;
547 : }
548 :
549 : inline const Complex *
550 234538916 : modelIn ()
551 : {
552 234538916 : assert (modelIn_p != 0);
553 234538916 : return modelIn_p;
554 : }
555 :
556 : inline Complex *
557 234538916 : modelOut ()
558 : {
559 234538916 : assert (modelOut_p != 0);
560 234538916 : return modelOut_p;
561 : }
562 :
563 : inline const Complex *
564 335197513 : observedIn ()
565 : {
566 335197513 : assert (observedIn_p != 0);
567 335197513 : return observedIn_p;
568 : }
569 :
570 : inline Complex *
571 335197513 : observedOut ()
572 : {
573 335197513 : assert (observedOut_p != 0);
574 335197513 : 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 527269705 : sigmaSpectrumIn ()
593 : {
594 527269705 : return sigmaSpectrumIn_p;
595 : }
596 :
597 : inline Float *
598 335197513 : sigmaSpectrumOut ()
599 : {
600 335197513 : assert (sigmaSpectrumOut_p != 0);
601 335197513 : return sigmaSpectrumOut_p;
602 : }
603 :
604 : inline const Float *
605 239440858 : weightSpectrumIn ()
606 : {
607 239440858 : assert (weightSpectrumIn_p != 0);
608 239440858 : return weightSpectrumIn_p;
609 : }
610 :
611 : inline Float *
612 776865919 : weightSpectrumOut ()
613 : {
614 776865919 : assert (weightSpectrumOut_p != 0);
615 776865919 : 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 16992 : distanceSquared (const Vector<T> & p1, const Vector<T> & p2)
739 : {
740 16992 : assert (p1.size() == 3 && p2.size() == 3);
741 :
742 16992 : T distanceSquared = 0;
743 :
744 67968 : for (Int i = 0; i < 3; i++){
745 50976 : T delta = p1[i] - p2[i];
746 50976 : distanceSquared += delta * delta;
747 : }
748 :
749 16992 : 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 607269 : MsRowAvg::MsRowAvg (rownr_t row, VbAvg * vb)
847 : : Vbi2MsRow (row, vb),
848 607269 : countsCache_p (& VbAvg::counts),
849 607269 : normalizationFactor_p(0.0),
850 1214538 : vbAvg_p (vb)
851 : {
852 607269 : configureCountsCache();
853 607269 : }
854 :
855 : Bool
856 22271554 : MsRowAvg::baselinePresent () const
857 : {
858 22271554 : return vbAvg_p->baselinePresent_p (row ());
859 : }
860 :
861 : Int
862 10957526 : MsRowAvg::nBaselinesPresent () const
863 : {
864 21915052 : return std::count(vbAvg_p->baselinePresent_p.begin(), vbAvg_p->baselinePresent_p.end(),
865 32872578 : true);
866 : }
867 :
868 : void
869 607269 : MsRowAvg::configureCountsCache ()
870 : {
871 607269 : addToCachedArrays (countsCache_p);
872 607269 : }
873 :
874 : const Matrix<Int> &
875 10986331 : MsRowAvg::counts () const
876 : {
877 10986331 : return countsCache_p.getCachedPlane (dynamic_cast<VbAvg *> (getVbi()), row());
878 : }
879 :
880 : Vector<Bool>
881 11447745 : MsRowAvg::correlationFlagsMutable ()
882 : {
883 11447745 : return vbAvg_p->correlationFlags_p.column (row());
884 : }
885 :
886 : Int
887 21944253 : MsRowAvg::countsBaseline () const
888 : {
889 21944253 : return vbAvg_p->countsBaseline_p (row ());
890 : }
891 :
892 : void
893 980438 : MsRowAvg::setBaselinePresent (Bool value)
894 : {
895 980438 : vbAvg_p->baselinePresent_p (row ()) = value;
896 980438 : }
897 :
898 :
899 : void
900 11447745 : MsRowAvg::setCountsBaseline (Int value)
901 : {
902 11447745 : vbAvg_p->countsBaseline_p (row ()) = value;
903 11447745 : }
904 :
905 : Double
906 440710 : MsRowAvg::intervalLast () const
907 : {
908 440710 : return vbAvg_p->intervalLast_p (row ());
909 : }
910 :
911 :
912 : Double
913 23107059 : MsRowAvg::timeFirst () const
914 : {
915 23107059 : return vbAvg_p->timeFirst_p (row ());
916 : }
917 :
918 : Double
919 490219 : MsRowAvg::timeLast () const
920 : {
921 490219 : return vbAvg_p->timeLast_p (row ());
922 : }
923 :
924 : Vector<Double>
925 507211 : MsRowAvg::uvwFirst ()
926 : {
927 507211 : return vbAvg_p->uvwFirst_p.column (row());
928 : }
929 :
930 :
931 : void
932 490219 : MsRowAvg::setCounts (const Matrix<Int> & value)
933 : {
934 490219 : Matrix<Int> & theCounts = countsCache_p.getCachedPlane (dynamic_cast<VbAvg *> (getVbi()), row());
935 490219 : theCounts = value;
936 490219 : }
937 :
938 : void
939 10957526 : MsRowAvg::setIntervalLast (Double value)
940 : {
941 10957526 : vbAvg_p->intervalLast_p (row ()) = value;
942 10957526 : }
943 :
944 :
945 : void
946 490219 : MsRowAvg::setTimeFirst (Double value)
947 : {
948 490219 : vbAvg_p->timeFirst_p (row ()) = value;
949 490219 : }
950 :
951 : void
952 11447745 : MsRowAvg::setTimeLast (Double value)
953 : {
954 11447745 : vbAvg_p->timeLast_p (row ()) = value;
955 11447745 : }
956 :
957 490219 : Double MsRowAvg::getNormalizationFactor()
958 : {
959 490219 : return vbAvg_p->normalizationFactor_p (row ());
960 : }
961 :
962 980538 : void MsRowAvg::setNormalizationFactor(Double normalizationFactor)
963 : {
964 980538 : vbAvg_p->normalizationFactor_p (row ()) = normalizationFactor;
965 980538 : }
966 :
967 10957526 : void MsRowAvg::accumulateNormalizationFactor(Double normalizationFactor)
968 : {
969 10957526 : vbAvg_p->normalizationFactor_p (row ()) += normalizationFactor;
970 10957526 : }
971 :
972 363 : VbAvg::VbAvg (const AveragingParameters & averagingParameters, ViImplementation2 * vii)
973 : : VisBufferImpl2 (vii, VbRekeyable),
974 726 : averagingInterval_p (averagingParameters.getAveragingInterval ()),
975 363 : averagingOptions_p (averagingParameters.getOptions()),
976 363 : averagingVii_p (vii),
977 363 : bufferToFill_p (0),
978 363 : complete_p (false),
979 363 : doing_p (), // all false until determined later on
980 363 : empty_p (true),
981 363 : maxTimeDistance_p (averagingParameters.getAveragingInterval() * (0.999)),
982 : // Shrink it just a bit for roundoff
983 363 : maxUvwDistanceSquared_p (pow(averagingParameters.getMaxUvwDistance(),2)),
984 363 : needIterationInfo_p (true),
985 363 : rowIdGenerator_p (0),
986 363 : sampleInterval_p (0),
987 363 : startTime_p (0),
988 1452 : usingUvwDistance_p (averagingParameters.getOptions().contains (AveragingOptions::BaselineDependentAveraging))
989 363 : {}
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 10957526 : calcOutRowIdx(Int nrows, Int baselineIndex, const MsRowAvg* rowAveraged, Int rowIdGen)
1003 : {
1004 10957526 : 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 10957526 : const bool multitime = nrows > nBasePresent;
1008 :
1009 10957526 : 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 10669943 : 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 287583 : Int rowIdG_div_baselines_roundup = 0;
1023 287583 : if (nBasePresent > 0)
1024 287583 : rowIdG_div_baselines_roundup = (rowIdGen + nBasePresent - 1)/ nBasePresent;
1025 287583 : const Int id = rowIdG_div_baselines_roundup * nBasePresent + baselineIndex;
1026 :
1027 287583 : return id;
1028 : }
1029 :
1030 : void
1031 592578 : VbAvg::accumulate (const VisBuffer2 * vb, const Subchunk & subchunk)
1032 : {
1033 592578 : if (empty_p){
1034 14691 : setupVbAvg (vb);
1035 : }
1036 :
1037 592578 : if (needIterationInfo_p){
1038 39522 : captureIterationInfo (bufferToFill_p, vb, subchunk);
1039 39522 : needIterationInfo_p = false;
1040 : }
1041 :
1042 592578 : assert (bufferToFill_p != 0);
1043 :
1044 592578 : MsRowAvg * rowAveraged = getRowMutable (0);
1045 592578 : MsRow * rowInput = vb->getRow (0);
1046 :
1047 592578 : auto nrows = vb->nRows();
1048 592578 : row2AvgRow_p.resize(nrows);
1049 11550104 : for (rownr_t row = 0; row < nrows; ++row){
1050 :
1051 10957526 : rowInput->changeRow (row);
1052 10957526 : Int baselineIndex = getBaselineIndex (rowInput);
1053 :
1054 10957526 : rowAveraged->changeRow (baselineIndex);
1055 :
1056 10957526 : accumulateOneRow (rowInput, rowAveraged, subchunk, row);
1057 :
1058 10957526 : row2AvgRow_p[row] = calcOutRowIdx(nrows, baselineIndex, rowAveraged,
1059 : rowIdGenerator_p);
1060 : }
1061 :
1062 592578 : delete rowAveraged;
1063 592578 : delete rowInput;
1064 :
1065 592578 : }
1066 :
1067 : void
1068 10957526 : VbAvg::accumulateOneRow (MsRow * rowInput, MsRowAvg * rowAveraged, const Subchunk & subchunk,
1069 : Int iidx)
1070 : {
1071 10957526 : finalizeBaselineIfNeeded (rowInput, rowAveraged, subchunk);
1072 :
1073 10957526 : if (! rowAveraged->baselinePresent())
1074 : {
1075 :
1076 490219 : initializeBaseline (rowInput, rowAveraged, subchunk);
1077 : }
1078 :
1079 : // bookkeeping - save for later that this input row is being averaged into the output row
1080 10957526 : 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 10957526 : Vector<Double> adjustedWeights;
1088 10957526 : Bool rowFlagged = false;
1089 :
1090 10957526 : std::tie (rowFlagged, adjustedWeights) = accumulateCubeData (rowInput, rowAveraged);
1091 :
1092 10957526 : Double adjustedWeight = 0;
1093 53591530 : for (Int c = 0; c < nCorrelations(); c++){
1094 :
1095 42634004 : adjustedWeight += adjustedWeights (c);
1096 : }
1097 :
1098 : // Accumulate the non matrix-valued data
1099 10957526 : accumulateRowData (rowInput, rowAveraged, adjustedWeight, rowFlagged);
1100 10957526 : }
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 52351 : VbAvg::finalizeBufferFilling ()
1147 : {
1148 52351 : bufferToFill_p->appendRowsComplete();
1149 52351 : bufferToFill_p = 0; // decouple
1150 52351 : }
1151 :
1152 : template<typename T>
1153 : inline void
1154 2113312911 : 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 2113312911 : if (zeroAccumulation){
1162 466163225 : * accumulator = (* unweightedValue) * weight;
1163 : }
1164 : else{
1165 1647149686 : * accumulator += (* unweightedValue) * weight;
1166 : }
1167 2113312911 : }
1168 :
1169 :
1170 : std::pair<Bool, Vector<Double> >
1171 10957526 : VbAvg::accumulateCubeData (MsRow * rowInput, MsRowAvg * rowAveraged)
1172 : {
1173 : // Accumulate the sums needed for averaging of cube data (e.g., visibility).
1174 :
1175 10957526 : const Matrix<Bool> & inputFlags = rowInput->flags ();
1176 10957526 : Matrix<Bool> & averagedFlags = rowAveraged->flagsMutable ();
1177 10957526 : Matrix<Int> counts = rowAveraged->counts ();
1178 10957526 : Vector<Bool> correlationFlagged = rowAveraged->correlationFlagsMutable ();
1179 :
1180 10957526 : AccumulationParameters accumulationParameters (rowInput, rowAveraged, doing_p);
1181 : // is a member variable to reduce memory allocations (jhj)
1182 :
1183 10957526 : IPosition shape = inputFlags.shape();
1184 10957526 : const Int nChannels = shape (1);
1185 10957526 : const Int nCorrelations = shape (0);
1186 :
1187 10957526 : Bool rowFlagged = true; // true if all correlations and all channels flagged
1188 :
1189 290129283 : for (Int channel = 0; channel < nChannels; channel ++){
1190 :
1191 1058497911 : 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 779326154 : Bool inputFlagged = inputFlags (correlation, channel);
1199 779326154 : if (rowFlagged && ! inputFlagged){
1200 10295784 : rowFlagged = false;
1201 : }
1202 : //rowFlagged = rowFlagged && inputFlagged;
1203 779326154 : Bool accumulatorFlagged = averagedFlags (correlation, channel);
1204 :
1205 779326154 : if (! accumulatorFlagged && inputFlagged){
1206 2460235 : accumulationParameters.incrementCubePointers();
1207 2460235 : 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 776865919 : Bool flagChange = (accumulatorFlagged && ! inputFlagged);
1214 776865919 : Bool zeroAccumulation = flagChange || counts (correlation, channel) == 0;
1215 :
1216 776865919 : if (flagChange){
1217 149937959 : averagedFlags (correlation, channel) = false;
1218 : }
1219 :
1220 776865919 : if (zeroAccumulation){
1221 156684479 : counts (correlation, channel) = 1;
1222 : }
1223 : else{
1224 620181440 : counts (correlation, channel) += 1;
1225 : }
1226 :
1227 : // Accumulate the sum for each cube element
1228 :
1229 776865919 : accumulateElementForCubes (accumulationParameters,
1230 : zeroAccumulation); // zeroes out accumulation
1231 :
1232 776865919 : accumulationParameters.incrementCubePointers();
1233 :
1234 : // Update correlation Flag
1235 :
1236 776865919 : if (correlationFlagged (correlation) && ! inputFlagged){
1237 1575714 : correlationFlagged (correlation) = false;
1238 : }
1239 : }
1240 : }
1241 :
1242 10957526 : Vector<Double> adjustedWeight = Vector<Double> (nCorrelations, 1);
1243 10957526 : if (doing_p.correctedData_p)
1244 : {
1245 48123933 : for (Int correlation = 0; correlation < nCorrelations; correlation ++)
1246 : {
1247 38417354 : adjustedWeight(correlation) = rowInput->weight(correlation);
1248 : }
1249 : }
1250 1250947 : else if (doing_p.observedData_p || doing_p.floatData_p)
1251 : {
1252 3920812 : for (Int correlation = 0; correlation < nCorrelations; correlation ++)
1253 : {
1254 2979222 : adjustedWeight(correlation) = AveragingTvi2::sigmaToWeight(rowInput->sigma (correlation));
1255 : }
1256 : }
1257 :
1258 21915052 : return std::make_pair (rowFlagged, adjustedWeight);
1259 10957526 : }
1260 :
1261 :
1262 : inline void
1263 776865919 : VbAvg::accumulateElementForCubes (AccumulationParameters & accumulationParameters,
1264 : Bool zeroAccumulation)
1265 : {
1266 :
1267 : // NOTE: The channelized flag check comes from the calling context (continue statement)
1268 776865919 : float weightCorrected = 1.0f;
1269 776865919 : float weightObserved = 1.0f;
1270 776865919 : const float One = 1.0f;
1271 :
1272 776865919 : if (doing_p.correctedData_p)
1273 : {
1274 : // The weight corresponding to CORRECTED_DATA is that stored in WEIGHT
1275 239440858 : weightCorrected = * accumulationParameters.weightSpectrumIn ();
1276 :
1277 :
1278 : // Accumulate weighted average contribution (normalization will come at the end)
1279 239440858 : accumulateElementForCube ( accumulationParameters.correctedIn (),
1280 : weightCorrected, zeroAccumulation,
1281 : accumulationParameters.correctedOut ());
1282 :
1283 : // The weight resulting from weighted average is the sum of the weights
1284 239440858 : accumulateElementForCube ( & weightCorrected,
1285 : One, zeroAccumulation,
1286 : accumulationParameters.weightSpectrumOut ());
1287 : }
1288 :
1289 776865919 : 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 335197513 : weightObserved = AveragingTvi2::sigmaToWeight(* accumulationParameters.sigmaSpectrumIn ());
1294 :
1295 : // Accumulate weighted average contribution (normalization will come at the end)
1296 :
1297 335197513 : accumulateElementForCube ( accumulationParameters.observedIn (),
1298 : weightObserved, zeroAccumulation,
1299 : accumulationParameters.observedOut ());
1300 :
1301 335197513 : if (not doing_p.correctedData_p)
1302 : {
1303 : // The weight resulting from weighted average is the sum of the weights
1304 266666226 : 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 335197513 : 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 776865919 : if (doing_p.modelData_p)
1325 : {
1326 234538916 : if (doing_p.correctedData_p)
1327 : {
1328 155852273 : accumulateElementForCube ( accumulationParameters.modelIn (),
1329 : weightCorrected, zeroAccumulation,
1330 : accumulationParameters.modelOut ());
1331 : }
1332 78686643 : else if (doing_p.observedData_p)
1333 : {
1334 0 : accumulateElementForCube ( accumulationParameters.modelIn (),
1335 : weightObserved, zeroAccumulation,
1336 : accumulationParameters.modelOut ());
1337 : }
1338 : else
1339 : {
1340 78686643 : 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 78686643 : accumulateElementForCube ( & One,
1346 : 1.0f, zeroAccumulation,
1347 : accumulationParameters.weightSpectrumOut ());
1348 : }
1349 : }
1350 :
1351 776865919 : if (doing_p.floatData_p)
1352 : {
1353 :
1354 : // The weight corresponding to FLOAT_DATA is that derived from the rms stored in SIGMA
1355 192072192 : weightObserved = AveragingTvi2::sigmaToWeight(* accumulationParameters.sigmaSpectrumIn ());
1356 :
1357 : // Accumulate weighted average contribution (normalization will come at the end)
1358 192072192 : accumulateElementForCube ( accumulationParameters.floatDataIn (),
1359 : weightObserved, zeroAccumulation,
1360 : accumulationParameters.floatDataOut ());
1361 :
1362 : // The weight resulting from weighted average is the sum of the weights
1363 192072192 : accumulateElementForCube ( & weightObserved,
1364 : 1.0f, zeroAccumulation,
1365 : accumulationParameters.weightSpectrumOut ());
1366 : }
1367 :
1368 1553731838 : return;
1369 : }
1370 :
1371 :
1372 :
1373 : template <typename T>
1374 : T
1375 43830104 : VbAvg::accumulateRowDatum (const T & averagedValue, const T & inputValue, Bool resetAverage)
1376 : {
1377 43830104 : if (resetAverage){
1378 1961276 : return inputValue;
1379 : }
1380 : else{
1381 41868828 : return inputValue + averagedValue;
1382 : }
1383 : }
1384 :
1385 : void
1386 10957526 : 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 10957526 : Bool accumulatorRowFlagged = rowAveraged->isRowFlagged();
1393 21454034 : Bool flagChange = accumulatorRowFlagged != rowFlagged || // actual change
1394 10496508 : rowAveraged->countsBaseline() == 0; // first time
1395 :
1396 10957526 : 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 10957526 : rowAveraged->setCountsBaseline (accumulateRowDatum (rowAveraged->countsBaseline(),
1405 10957526 : 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 10957526 : accumulatorRowFlagged = accumulatorRowFlagged && rowFlagged;
1413 10957526 : rowAveraged->setRowFlag (accumulatorRowFlagged);
1414 :
1415 10957526 : rowAveraged->setExposure (accumulateRowDatum (rowAveraged->exposure(),
1416 10957526 : 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 10957526 : if (accumulatorRowFlagged){
1425 661742 : weightToUse = 1;
1426 : }
1427 : else{
1428 10295784 : weightToUse = adjustedWeight;
1429 : }
1430 :
1431 10957526 : if (flagChange){
1432 490319 : rowAveraged->setNormalizationFactor(0.0);
1433 : }
1434 :
1435 10957526 : rowAveraged->accumulateNormalizationFactor(weightToUse);
1436 :
1437 10957526 : Double weightedTC = (rowInput->timeCentroid() - rowAveraged->timeFirst()) * weightToUse;
1438 10957526 : rowAveraged->setTimeCentroid (accumulateRowDatum (rowAveraged->timeCentroid(),
1439 : weightedTC,
1440 : flagChange));
1441 :
1442 21915052 : Vector<Double> weightedUvw = rowInput->uvw() * weightToUse;
1443 :
1444 10957526 : rowAveraged->setUvw (accumulateRowDatum (rowAveraged->uvw (),
1445 : weightedUvw,
1446 : flagChange));
1447 :
1448 : // Capture a couple pieces of data
1449 :
1450 10957526 : rowAveraged->setTimeLast (rowInput->time());
1451 :
1452 10957526 : rowAveraged->setIntervalLast (rowInput->interval());
1453 10957526 : }
1454 10957526 : }
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 10960983 : VbAvg::counts () const
1487 : {
1488 10960983 : return counts_p;
1489 : }
1490 :
1491 :
1492 : void
1493 490219 : VbAvg::copyIdValues (MsRow * rowInput, MsRowAvg * rowAveraged)
1494 : {
1495 490219 : rowAveraged->setAntenna1 (rowInput->antenna1());
1496 490219 : rowAveraged->setAntenna2 (rowInput->antenna2());
1497 490219 : rowAveraged->setArrayId (rowInput->arrayId());
1498 490219 : rowAveraged->setDataDescriptionId (rowInput->dataDescriptionId());
1499 490219 : rowAveraged->setFeed1 (rowInput->feed1());
1500 490219 : rowAveraged->setFeed2 (rowInput->feed2());
1501 490219 : rowAveraged->setFieldId (rowInput->fieldId());
1502 490219 : rowAveraged->setProcessorId (rowInput->processorId());
1503 490219 : rowAveraged->setScanNumber (rowInput->scanNumber());
1504 490219 : rowAveraged->setSpectralWindow (rowInput->spectralWindow());
1505 490219 : rowAveraged->setObservationId (rowInput->observationId());
1506 490219 : rowAveraged->setStateId (rowInput->stateId());
1507 490219 : }
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 10957526 : 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 10957526 : const Int antenna1 = msRow->antenna1 ();
1537 10957526 : const Int antenna2 = msRow->antenna2 ();
1538 10957526 : const Int spw = msRow->spectralWindow ();
1539 :
1540 10957526 : const Int index = baselineIndex_p (antenna1, antenna2, spw);
1541 :
1542 10957526 : 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 27424 : VbAvg::finalizeAverages ()
1555 : {
1556 27424 : if (empty_p){
1557 12733 : return; // nothing to finalize
1558 : }
1559 :
1560 14691 : MsRowAvg * msRowAvg = getRowMutable (0);
1561 :
1562 371193 : for (Int baseline = 0; baseline < nBaselines(); baseline ++){
1563 :
1564 356502 : msRowAvg->changeRow (baseline);
1565 :
1566 356502 : if (msRowAvg->baselinePresent()){
1567 278650 : finalizeBaseline (msRowAvg);
1568 : }
1569 :
1570 : }
1571 :
1572 14691 : delete msRowAvg;
1573 :
1574 14691 : empty_p = true;
1575 : }
1576 :
1577 : void
1578 490219 : 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 490219 : msRowAvg->setRowFlag(false);
1585 :
1586 490219 : finalizeCubeData (msRowAvg);
1587 :
1588 490219 : finalizeRowData (msRowAvg);
1589 :
1590 490219 : transferBaseline (msRowAvg);
1591 490219 : }
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 : public std::binary_function<L,R,RES>
1597 : {
1598 170654064 : RES operator() (const L& x, const R& y) const
1599 : {
1600 341308128 : { return y > 0? RES(x)/y : RES(x); }
1601 : }
1602 : };
1603 :
1604 :
1605 : void
1606 490219 : VbAvg::finalizeCubeData (MsRowAvg * msRowAvg)
1607 : {
1608 : // Divide each of the data cubes in use by the sum of the appropriate weights.
1609 :
1610 : typedef DividesNonZero <Complex, Float, Complex> DivideOp;
1611 : DivideOp op;
1612 :
1613 490219 : if (doing_p.correctedData_p)
1614 : {
1615 251047 : Matrix<Complex> corrected = msRowAvg->correctedMutable();
1616 251047 : arrayTransformInPlace<Complex, Float, DivideOp > (corrected,msRowAvg->weightSpectrum (), op);
1617 251047 : }
1618 :
1619 490219 : if (doing_p.observedData_p)
1620 : {
1621 237766 : Matrix<Complex> observed = msRowAvg->observedMutable();
1622 237766 : if (not doing_p.correctedData_p)
1623 209686 : arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->weightSpectrum (), op);
1624 : else
1625 28080 : arrayTransformInPlace<Complex, Float, DivideOp > (observed,msRowAvg->sigmaSpectrum (), op);
1626 237766 : }
1627 :
1628 490219 : if (doing_p.modelData_p)
1629 : {
1630 249053 : Matrix<Complex> model = msRowAvg->modelMutable();
1631 :
1632 249053 : if (doing_p.correctedData_p)
1633 220248 : arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->weightSpectrum (), op);
1634 28805 : else if (doing_p.observedData_p)
1635 0 : arrayTransformInPlace<Complex, Float, DivideOp > (model,msRowAvg->sigmaSpectrum (), op);
1636 : else
1637 28805 : arrayTransformInPlace<Complex, Int, DivideOp > (model,msRowAvg->counts (), op);
1638 249053 : }
1639 :
1640 490219 : if (doing_p.floatData_p)
1641 : {
1642 : typedef Divides <Float, Float, Float> DivideOpFloat;
1643 : DivideOpFloat opFloat;
1644 :
1645 681 : Matrix<Float> visCubeFloat = msRowAvg->singleDishDataMutable();
1646 681 : arrayTransformInPlace<Float, Float, DivideOpFloat > (visCubeFloat,msRowAvg->weightSpectrum (), opFloat);
1647 681 : }
1648 :
1649 :
1650 980438 : return;
1651 : }
1652 :
1653 : void
1654 490219 : VbAvg::finalizeRowData (MsRowAvg * msRow)
1655 : {
1656 490219 : Int n = msRow->countsBaseline ();
1657 :
1658 : // Obtain row-level WEIGHT by calculating the mean of WEIGHT_SPECTRUM
1659 : // msRow->setWeight(partialMedians(msRow->weightSpectrum(),IPosition(1,1),true));
1660 490219 : msRow->setWeight(AveragingTvi2::average(msRow->weightSpectrum(),msRow->flags()));
1661 :
1662 : // If doing both DATA and CORRECTED_DATA then SIGMA_SPECTRUM contains the weight
1663 : // (not sigma) accumulation for DATA, and we have to derive SIGMA from it
1664 490219 : if (doing_p.correctedData_p and doing_p.observedData_p)
1665 : {
1666 : // jagonzal: SIGMA is not derived from the mean of SIGMA_SPECTRUM but from the mean
1667 : // of the WEIGHT format of SIGMA_SPECTRUM turned into SIGMA by using 1/pow(weight,2)
1668 : // Vector<Float> weight = partialMedians(msRow->sigmaSpectrum(),IPosition(1,1),true);
1669 28080 : Vector<Float> weight = AveragingTvi2::average(msRow->sigmaSpectrum(),msRow->flags());
1670 28080 : arrayTransformInPlace (weight, AveragingTvi2::weightToSigma);
1671 28080 : msRow->setSigma (weight);
1672 :
1673 : // Now convert the DATA weight accumulation stored in sigmaSpectrum into the real SIGMA_SPECTRUM
1674 : // TODO: This should happen only if we are writing out SIGMA_SPECTRUM but
1675 : // multiple column operation is rare and might be forbidden in the future
1676 28080 : Matrix<Float> sigmaSpectrun = msRow->sigmaSpectrum(); // Reference copy
1677 28080 : arrayTransformInPlace (sigmaSpectrun, AveragingTvi2::weightToSigma);
1678 28080 : }
1679 : // Otherwise (doing only DATA/FLOAT_DATA or CORRECTED_DATA) we can derive SIGMA from WEIGHT directly
1680 462139 : else if ( not doing_p.onlymetadata_p)
1681 : {
1682 : // jagonzal: SIGMA is not derived from the mean of SIGMA_SPECTRUM
1683 : // but from WEIGHT turned into SIGMA by using 1/pow(weight,2)
1684 462139 : Vector<Float> sigma = msRow->sigma(); // Reference copy
1685 462139 : sigma = msRow->weight(); // Normal copy (transfer Weight values to Sigma)
1686 462139 : arrayTransformInPlace (sigma, AveragingTvi2::weightToSigma);
1687 :
1688 : // Derive SIGMA_SPECTRUM from computed WEIGHT_SPECTRUM
1689 462139 : Matrix<Float> sigmaSpectrun = msRow->sigmaSpectrum(); // Reference copy
1690 462139 : sigmaSpectrun = msRow->weightSpectrum(); // Normal copy (transfer WeightSpectrum values to SigmaSpectrum)
1691 462139 : arrayTransformInPlace (sigmaSpectrun, AveragingTvi2::weightToSigma);
1692 462139 : }
1693 :
1694 : // Get the normalization factor for this baseline, containing
1695 : // the accumulation of row-level) weight contributions
1696 490219 : Double weight = msRow->getNormalizationFactor();
1697 :
1698 490219 : if (n != 0){
1699 :
1700 490219 : if (weight == 0){
1701 :
1702 : // The weights are all zero so compute an arithmetic average
1703 : // so that a somewhat value can go into these columns (i.e. rather than NaN).
1704 :
1705 0 : weight = msRow->countsBaseline();
1706 : }
1707 :
1708 490219 : msRow->setTimeCentroid (msRow->timeCentroid() / weight + msRow->timeFirst());
1709 :
1710 490219 : msRow->setUvw (msRow->uvw() / weight);
1711 :
1712 : // Exposure is a simple sum, not an average so it is already
1713 : // done at this point.
1714 : }
1715 :
1716 : // Fill in the time and the interval
1717 : //
1718 : // The time of a sample is centered over the integration time period.
1719 : // Thus find the center of the averaged interval it is necessary to
1720 : // slide it back by 1/2 an interval.
1721 :
1722 490219 : Double dT = msRow->timeLast () - msRow->timeFirst();
1723 :
1724 490219 : Double centerOfInterval = msRow->timeFirst () + dT / 2;
1725 :
1726 490219 : msRow->setTime (centerOfInterval);
1727 :
1728 490219 : if (dT != 0){
1729 :
1730 : // The interval is the center-to-center time + half of the intervals of
1731 : // the first and the last sample (if dT == 0 then the interval is
1732 : // already the interval of the first sample and is correct)
1733 :
1734 440710 : Double interval = dT + msRow->interval() / 2 + msRow->intervalLast() / 2;
1735 440710 : msRow->setInterval (interval);
1736 : }
1737 490219 : }
1738 :
1739 : void
1740 10957526 : VbAvg::finalizeBaselineIfNeeded (MsRow * rowInput, MsRowAvg * rowAveraged, const Subchunk & /*subchunk*/)
1741 : {
1742 10957526 : if (! rowAveraged->baselinePresent()){
1743 278650 : return;
1744 : }
1745 :
1746 : // Finalization is needed if either the uvw distance or the time distance between the input
1747 : // baseline and the averaged baseline is above the maximum
1748 :
1749 10678876 : Bool needed = usingUvwDistance_p;
1750 :
1751 10678876 : if (needed) {
1752 16992 : Double deltaUvw = distanceSquared (rowInput->uvw(), rowAveraged->uvwFirst ());
1753 16992 : needed = deltaUvw > maxUvwDistanceSquared_p;
1754 : }
1755 :
1756 10678876 : needed = needed || (rowInput->time() - rowAveraged->timeFirst()) > maxTimeDistance_p;
1757 :
1758 10678876 : if (needed){
1759 :
1760 : // Finalize the data so that the final averaging products and then move them to
1761 : // output buffer.
1762 :
1763 211569 : finalizeBaseline (rowAveraged);
1764 : }
1765 : }
1766 :
1767 : MsRowAvg *
1768 0 : VbAvg::getRow (Int row) const
1769 : {
1770 0 : return new MsRowAvg (row, this);
1771 : }
1772 :
1773 : MsRowAvg *
1774 607269 : VbAvg::getRowMutable (Int row)
1775 : {
1776 607269 : return new MsRowAvg (row, this);
1777 : }
1778 :
1779 : void
1780 490219 : VbAvg::initializeBaseline (MsRow * rowInput, MsRowAvg * rowAveraged,
1781 : const Subchunk & /*subchunk*/)
1782 : {
1783 490219 : copyIdValues (rowInput, rowAveraged);
1784 :
1785 : // Size and zero out the counters
1786 :
1787 490219 : rowAveraged->setInterval (rowInput->interval()); // capture first one
1788 490219 : rowAveraged->setTimeFirst (rowInput->time());
1789 490219 : rowAveraged->setTimeLast (rowInput->time());
1790 490219 : rowAveraged->uvwFirst () = rowInput->uvw ();
1791 :
1792 490219 : rowAveraged->setCountsBaseline (0);
1793 :
1794 490219 : IPosition shape = rowInput->flags().shape();
1795 490219 : Int nCorrelations = shape (0);
1796 490219 : Int nChannels = shape (1);
1797 :
1798 490219 : rowAveraged->setCounts (zeroInt_p.getMatrix (nCorrelations, nChannels, 0));
1799 490219 : rowAveraged->setWeight (Vector<Float> (nCorrelations, 0));
1800 490219 : rowAveraged->setTimeCentroid (0.0);
1801 :
1802 490219 : if (doing_p.weightSpectrumOut_p){
1803 490219 : rowAveraged->setWeightSpectrum (zeroFloat_p.getMatrix (nCorrelations, nChannels, 0));
1804 : }
1805 :
1806 490219 : if (doing_p.sigmaSpectrumOut_p){
1807 490219 : rowAveraged->setSigmaSpectrum (zeroFloat_p.getMatrix (nCorrelations, nChannels, 0));
1808 : }
1809 :
1810 : // VisBufferComponents2 exclusions =
1811 : // VisBufferComponents2::these(VisibilityObserved, VisibilityCorrected,
1812 : // VisibilityModel, CorrType, JonesC, Unknown);
1813 : //
1814 : // cacheResizeAndZero(exclusions);
1815 :
1816 : // Flag everything to start with
1817 :
1818 490219 : rowAveraged->setRowFlag (true); // only for use during row-value accumulation
1819 490219 : rowAveraged->setFlags(trueBool_p.getMatrix (nCorrelations, nChannels, true));
1820 490219 : rowAveraged->correlationFlagsMutable() = Vector<Bool> (nCorrelations, true);
1821 :
1822 490219 : rowAveraged->setBaselinePresent(true);
1823 :
1824 490219 : rowAveraged->setNormalizationFactor(0.0);
1825 490219 : }
1826 :
1827 :
1828 : Bool
1829 0 : VbAvg::isComplete () const
1830 : {
1831 0 : return complete_p;
1832 : }
1833 :
1834 : Bool
1835 592501 : VbAvg::isUsingUvwDistance () const
1836 : {
1837 592501 : return usingUvwDistance_p;
1838 : }
1839 :
1840 :
1841 : //void
1842 : //VbAvg::markEmpty ()
1843 : //{
1844 : // empty_p = true;
1845 : // complete_p = false;
1846 : //}
1847 :
1848 : Int
1849 900934 : VbAvg::nBaselines () const
1850 : {
1851 900934 : return countsBaseline_p.nelements();
1852 : }
1853 :
1854 : Int
1855 592501 : VbAvg::nSpectralWindowsInBuffer () const
1856 : {
1857 592501 : const Vector<Int> & windows = spectralWindows();
1858 592501 : set <Int> s;
1859 :
1860 15842321 : for (uInt i = 0; i < windows.nelements(); i ++){
1861 15249820 : s.insert (windows(i));
1862 : }
1863 :
1864 1185002 : return (Int) s.size();
1865 :
1866 592501 : }
1867 :
1868 :
1869 : void
1870 39522 : VbAvg::captureIterationInfo (VisBufferImpl2 * dstVb, const VisBuffer2 * srcVb,
1871 : const Subchunk & subchunk)
1872 : {
1873 39522 : dstVb->setIterationInfo (srcVb->msId(),
1874 39522 : srcVb->msName(),
1875 39522 : srcVb->isNewMs(),
1876 39522 : srcVb->isNewArrayId(),
1877 39522 : srcVb->isNewFieldId(),
1878 39522 : srcVb->isNewSpectralWindow(),
1879 : subchunk,
1880 79044 : srcVb->getCorrelationTypes (),
1881 79044 : srcVb->getCorrelationTypesDefined(),
1882 79044 : srcVb->getCorrelationTypesSelected(),
1883 79044 : CountedPtr <WeightScaling> ( ));
1884 :
1885 : // Request info from the VB which will cause it to be filled
1886 : // into cache from the input VII at its current position.
1887 :
1888 39522 : dstVb->setRekeyable(true);
1889 39522 : dstVb->setShape(srcVb->nCorrelations(), srcVb->nChannels(), nBaselines(), false);
1890 : // Do not clear the cache since we're resuing the storage
1891 :
1892 39522 : dstVb->phaseCenter();
1893 39522 : dstVb->nAntennas();
1894 39522 : dstVb->correlationTypes();
1895 39522 : dstVb->polarizationFrame();
1896 39522 : dstVb->polarizationId();
1897 39522 : }
1898 :
1899 : //void
1900 : //VbAvg::prepareForFirstData (const VisBuffer2 * vb, const Subchunk & subchunk)
1901 : //{
1902 : // startTime_p = vb->time() (0);
1903 : // sampleInterval_p = vb->timeInterval() (0);
1904 : //
1905 : // Int nAntennas = vb->nAntennas();
1906 : // Int nSpw = vb->getVi()->nSpectralWindows();
1907 : // Int nBaselines = ((nAntennas * (nAntennas + 1)) / 2) * nSpw;
1908 : //
1909 : // // Size and zero out the counters
1910 : //
1911 : // timeFirst_p = Vector<Double> (nBaselines, vb->time() (0));
1912 : // timeLast_p = Vector<Double> (nBaselines, vb->time() (0));
1913 : // uvwFirst_p = Vector<Double> (nBaselines, vb->uvw().column(0));
1914 : //
1915 : // countsBaseline_p = Vector<Int> (nBaselines, 0);
1916 : // counts_p = Cube<Int> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1917 : // weightSum_p = Cube<Float> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1918 : // if (doing_p.sigmaSpectrum_p){
1919 : // weightCorrectedSum_p = Cube<Float> (vb->nCorrelations(), vb->nChannels(), nBaselines, 0);
1920 : // }
1921 : //
1922 : // baselineIndex_p.configure (nAntennas, nSpw, vb);
1923 : //
1924 : // // Reshape the inherited members from VisBuffer2
1925 : //
1926 : // captureIterationInfo (vb, subchunk);
1927 : //
1928 : // setShape (vb->nCorrelations(), vb->nChannels(), nBaselines, false);
1929 : //
1930 : // VisBufferComponents2 exclusions =
1931 : // VisBufferComponents2::these(VisibilityObserved, VisibilityCorrected,
1932 : // VisibilityModel, CorrType, JonesC, Unknown);
1933 : // cacheResizeAndZero(exclusions);
1934 : //
1935 : // prepareIds (vb);
1936 : //
1937 : // // Flag everything to start with
1938 : //
1939 : // setFlagCube (Cube<Bool> (vb->nCorrelations(), vb->nChannels(), nBaselines, true));
1940 : // setFlagRow (Vector<Bool> (nBaselines, true));
1941 : //
1942 : // complete_p = false;
1943 : //}
1944 :
1945 : void
1946 0 : VbAvg::prepareIds (const VisBuffer2 * vb)
1947 : {
1948 : // Set these row ID values to indicate they are unknown
1949 :
1950 0 : Vector<Int> minusOne (nBaselines(), -1);
1951 :
1952 0 : setAntenna1 (minusOne);
1953 0 : setAntenna2 (minusOne);
1954 0 : setDataDescriptionIds (minusOne);
1955 0 : setFeed1 (minusOne);
1956 0 : setFeed2 (minusOne);
1957 0 : setProcessorId (minusOne);
1958 0 : setScan (minusOne);
1959 0 : setObservationId (minusOne);
1960 0 : setSpectralWindows (minusOne);
1961 0 : setStateId (minusOne);
1962 :
1963 : // Copy the value from the input VB
1964 :
1965 0 : Vector<Int> tmp (nBaselines(), vb->arrayId()(0));
1966 :
1967 0 : setArrayId (tmp);
1968 :
1969 0 : tmp = vb->dataDescriptionIds()(0);
1970 0 : setDataDescriptionIds (tmp);
1971 :
1972 0 : tmp = vb->fieldId()(0);
1973 0 : setFieldId (tmp);
1974 0 : }
1975 :
1976 : //void
1977 : //VbAvg::removeMissingBaselines ()
1978 : //{
1979 : // // Some baselines may not be present in the portion of the input data
1980 : // // that made up this average. However, this is not known until after
1981 : // // all of the data is seen. Thus at finalization time these bogus
1982 : // // baselines should be removed from the average so as not to pass
1983 : // // flagged but zero-exposure baselines into the output.
1984 : //
1985 : //
1986 : // Vector<Int> rowsToDelete (nBaselines());
1987 : //
1988 : // Int nBaselinesDeleted = 0;
1989 : //
1990 : // for (Int baseline = 0; baseline < nBaselines(); baseline ++){
1991 : //
1992 : // if (countsBaseline_p (baseline) == 0){
1993 : // rowsToDelete (nBaselinesDeleted) = baseline;
1994 : // nBaselinesDeleted ++;
1995 : // }
1996 : // }
1997 : //
1998 : // rowsToDelete.resize (nBaselinesDeleted, true);
1999 : //
2000 : // deleteRows (rowsToDelete);
2001 : //}
2002 :
2003 : void
2004 52351 : VbAvg::setBufferToFill(VisBufferImpl2 * vb)
2005 : {
2006 52351 : bufferToFill_p = vb;
2007 :
2008 : // Set flag so that iteration information will be captured into
2009 : // this VB from the first input VB.
2010 :
2011 52351 : needIterationInfo_p = true;
2012 52351 : }
2013 :
2014 : void
2015 14691 : VbAvg::setupVbAvg (const VisBuffer2 * vb)
2016 : {
2017 : // Configure the index
2018 :
2019 14691 : Int nAntennas = vb->nAntennas();
2020 :
2021 : // This is a kluge to allow multiple spectral windows (of the same shape)
2022 : // to be combined into a single VB. This really shouldn't be allowed!!!
2023 :
2024 14691 : set<uInt> spwInVb;
2025 :
2026 313615 : for (rownr_t i = 0; i < vb->nRows(); i++){
2027 298924 : spwInVb.insert (vb->dataDescriptionIds()(i));
2028 : }
2029 :
2030 14691 : uInt nSpwInVb = spwInVb.size();
2031 :
2032 14691 : Int nSpw = averagingVii_p->nSpectralWindows();
2033 :
2034 14691 : baselineIndex_p.configure (nAntennas, nSpw, vb);
2035 :
2036 14691 : Int nBaselines = ((nAntennas * (nAntennas + 1)) / 2)* nSpwInVb;
2037 :
2038 14691 : setShape (vb->nCorrelations(), vb->nChannels(), nBaselines);
2039 :
2040 14691 : setupArrays (vb->nCorrelations(), vb->nChannels(), nBaselines);
2041 :
2042 14691 : baselinePresent_p.resize(nBaselines);
2043 14691 : baselinePresent_p = False;
2044 :
2045 14691 : normalizationFactor_p.resize(nBaselines);
2046 14691 : normalizationFactor_p = 0.0;
2047 :
2048 14691 : empty_p = false;
2049 14691 : }
2050 :
2051 : void
2052 14691 : VbAvg::setupArrays (Int nCorrelations, Int nChannels, Int nBaselines)
2053 : {
2054 :
2055 14691 : setShape (nCorrelations, nChannels, nBaselines);
2056 :
2057 : // Start out with all of the array-valued components except the
2058 : // optional ones.
2059 :
2060 : VisBufferComponents2 including =
2061 : VisBufferComponents2::these ({VisBufferComponent2::Antenna1,
2062 : VisBufferComponent2::Antenna2,
2063 : VisBufferComponent2::ArrayId,
2064 : VisBufferComponent2::CorrType,
2065 : VisBufferComponent2::DataDescriptionIds,
2066 : VisBufferComponent2::Exposure,
2067 : VisBufferComponent2::Feed1,
2068 : VisBufferComponent2::Feed2,
2069 : VisBufferComponent2::FieldId,
2070 : VisBufferComponent2::FlagCategory,
2071 : VisBufferComponent2::FlagCube,
2072 : VisBufferComponent2::FlagRow,
2073 : VisBufferComponent2::ObservationId,
2074 : VisBufferComponent2::ProcessorId,
2075 : VisBufferComponent2::RowIds,
2076 : VisBufferComponent2::Scan,
2077 : VisBufferComponent2::Sigma,
2078 : VisBufferComponent2::SpectralWindows,
2079 : VisBufferComponent2::StateId,
2080 : VisBufferComponent2::Time,
2081 : VisBufferComponent2::TimeCentroid,
2082 : VisBufferComponent2::TimeInterval,
2083 : VisBufferComponent2::Weight,
2084 14691 : VisBufferComponent2::Uvw});
2085 :
2086 14691 : if (doing_p.modelData_p){
2087 4307 : including += VisBufferComponent2::VisibilityCubeModel;
2088 : }
2089 :
2090 14691 : if (doing_p.correctedData_p){
2091 4346 : including += VisBufferComponent2::VisibilityCubeCorrected;
2092 : }
2093 :
2094 14691 : if (doing_p.observedData_p){
2095 9788 : including += VisBufferComponent2::VisibilityCubeObserved;
2096 : }
2097 :
2098 14691 : if (doing_p.floatData_p){
2099 498 : including += VisBufferComponent2::VisibilityCubeFloat;
2100 : }
2101 :
2102 14691 : if (doing_p.weightSpectrumOut_p){
2103 14691 : including += VisBufferComponent2::WeightSpectrum;
2104 : }
2105 :
2106 14691 : if (doing_p.sigmaSpectrumOut_p){
2107 14691 : including += VisBufferComponent2::SigmaSpectrum;
2108 : }
2109 :
2110 14691 : cacheResizeAndZero (including);
2111 :
2112 14691 : correlationFlags_p.reference (Matrix<Bool> (IPosition (2, nCorrelations, nBaselines), true));
2113 14691 : counts_p.reference (Cube<Int> (IPosition (3, nCorrelations, nChannels, nBaselines), 0));
2114 14691 : countsBaseline_p.reference (Vector<Int> (nBaselines, 0)); // number of items summed together for each baseline.
2115 14691 : intervalLast_p.reference (Vector<Double> (nBaselines, 0));
2116 14691 : timeFirst_p.reference (Vector<Double> (nBaselines, 0));
2117 14691 : timeLast_p.reference (Vector<Double> (nBaselines, 0));
2118 14691 : uvwFirst_p.reference (Matrix<Double> (IPosition (2, 3, nBaselines), 0.0));
2119 14691 : }
2120 :
2121 : void
2122 12530 : VbAvg::startChunk (ViImplementation2 * vi)
2123 : {
2124 12530 : empty_p = true;
2125 :
2126 12530 : rowIdGenerator_p = 0;
2127 12530 : row2AvgRow_p.resize(0);
2128 :
2129 : // See if the new MS has corrected and/or model data columns
2130 :
2131 12530 : doing_p.observedData_p = averagingOptions_p.contains (AveragingOptions::AverageObserved);
2132 22013 : doing_p.correctedData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeCorrected) &&
2133 9483 : averagingOptions_p.contains (AveragingOptions::AverageCorrected);
2134 25060 : doing_p.modelData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeModel) &&
2135 12530 : averagingOptions_p.contains (AveragingOptions::AverageModel);
2136 13066 : doing_p.floatData_p = vi->existsColumn (VisBufferComponent2::VisibilityCubeFloat) &&
2137 536 : averagingOptions_p.contains (AveragingOptions::AverageFloat);
2138 :
2139 12530 : doing_p.weightSpectrumIn_p = doing_p.correctedData_p;
2140 12530 : doing_p.sigmaSpectrumIn_p = doing_p.observedData_p || doing_p.floatData_p;
2141 12530 : doing_p.weightSpectrumOut_p = true; // We always use the output WeightSpectrum
2142 12530 : doing_p.sigmaSpectrumOut_p = true; // We always use the output SigmaSpectrum
2143 :
2144 12530 : if (doing_p.observedData_p or doing_p.correctedData_p or doing_p.modelData_p or doing_p.floatData_p)
2145 : {
2146 12530 : doing_p.onlymetadata_p = false;
2147 : }
2148 :
2149 : // Set up the flags for row copying
2150 :
2151 12530 : optionalComponentsToCopy_p = VisBufferComponents2::none();
2152 :
2153 12530 : if (doing_p.observedData_p){
2154 4082 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeObserved;
2155 : }
2156 :
2157 12530 : if (doing_p.correctedData_p){
2158 7851 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeCorrected;
2159 : }
2160 :
2161 12530 : if (doing_p.modelData_p){
2162 7806 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeModel;
2163 : }
2164 :
2165 12530 : if (doing_p.floatData_p){
2166 536 : optionalComponentsToCopy_p += VisBufferComponent2::VisibilityCubeFloat;
2167 : }
2168 :
2169 12530 : if (doing_p.weightSpectrumOut_p){
2170 12530 : optionalComponentsToCopy_p += VisBufferComponent2::WeightSpectrum;
2171 : }
2172 :
2173 12530 : if (doing_p.sigmaSpectrumOut_p){
2174 12530 : optionalComponentsToCopy_p += VisBufferComponent2::SigmaSpectrum;
2175 : }
2176 12530 : }
2177 :
2178 : void
2179 490219 : VbAvg::transferBaseline (MsRowAvg * rowAveraged)
2180 : {
2181 490219 : rowAveraged->setRowId (rowIdGenerator_p ++);
2182 490219 : bufferToFill_p->appendRow (rowAveraged, nBaselines(), optionalComponentsToCopy_p);
2183 :
2184 490219 : rowAveraged->setBaselinePresent(false);
2185 490219 : }
2186 :
2187 :
2188 : //VbSet::VbSet (const AveragingParameters & averagingParameters)
2189 : //: averagingInterval_p (averagingParameters.getAveragingInterval ()),
2190 : // averagingOptions_p (averagingParameters.getOptions()),
2191 : // averagingParameters_p (averagingParameters),
2192 : // doingCorrectedData_p (false),
2193 : // doingModelData_p (false),
2194 : // doingObservedData_p (false),
2195 : // doingWeightSpectrum_p (false),
2196 : // doingsigmaSpectrum_p (false),
2197 : // vbAveragers_p ()
2198 : //{}
2199 :
2200 : //VbSet::~VbSet ()
2201 : //{
2202 : // flush (true); // allow killing nonempty buffers
2203 : //}
2204 : //
2205 : //void
2206 : //VbSet::accumulate (const VisBuffer2 * input, const Subchunk & subchunk)
2207 : //{
2208 : // Int ddId = input->dataDescriptionIds()(0);
2209 : //
2210 : // if (! utilj::containsKey (ddId, vbAveragers_p)){ // Need a new averager
2211 : // add (ddId);
2212 : // }
2213 : //
2214 : // VbAvg * vba = vbAveragers_p [ddId];
2215 : //
2216 : // vba->accumulate (input, subchunk);
2217 : //}
2218 : //
2219 : //VbAvg *
2220 : //VbSet::add (Int ddId)
2221 : //{
2222 : // VbAvg::Doing doing;
2223 : //
2224 : // doing.correctedData_p = doingCorrectedData_p;
2225 : // doing.modelData_p = doingModelData_p;
2226 : // doing.observedData_p = doingObservedData_p;
2227 : // doing.weightSpectrum_p = doingWeightSpectrum_p;
2228 : // doing.sigmaSpectrum_p = doingsigmaSpectrum_p;
2229 : //
2230 : // VbAvg * newAverager = new VbAvg (doing, averagingParameters_p);
2231 : //
2232 : // vbAveragers_p [ddId] = newAverager;
2233 : //
2234 : // return newAverager;
2235 : //}
2236 : //
2237 : //Bool
2238 : //VbSet::anyAveragesReady(Int ddid) const
2239 : //{
2240 : // Bool any = false;
2241 : //
2242 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2243 : // a != vbAveragers_p.end();
2244 : // a++){
2245 : //
2246 : // if (a->second->isComplete() &&
2247 : // (ddid < 0 || ddid == a->second->dataDescriptionIds()(0))){
2248 : //
2249 : // any = true;
2250 : // break;
2251 : // }
2252 : // }
2253 : //
2254 : // return any;
2255 : //}
2256 : //
2257 : //Bool
2258 : //VbSet::contains (Int ddId) const
2259 : //{
2260 : // return vbAveragers_p.find (ddId) != vbAveragers_p.end();
2261 : //}
2262 : //
2263 : ////void
2264 : ////VbSet::finalizeAverage (Int ddId)
2265 : ////{
2266 : //// vbAveragers_p [ddId]->finalizeAverage();
2267 : ////}
2268 : //
2269 : //void
2270 : //VbSet::finalizeAverages ()
2271 : //{
2272 : //// for (Averagers::iterator a = vbAveragers_p.begin();
2273 : //// a != vbAveragers_p.end();
2274 : //// a ++){
2275 : ////
2276 : //// finalizeAverage (a->first);
2277 : //// }
2278 : //}
2279 : //
2280 : //void
2281 : //VbSet::flush (Bool okIfNonempty, ViImplementation2 * vi)
2282 : //{
2283 : // // Get rid of all of the averagers. This is done at
2284 : // // destruction time and when a sweeping into a new MS.
2285 : //
2286 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2287 : // a != vbAveragers_p.end();
2288 : // a ++){
2289 : //
2290 : // Assert (okIfNonempty || (a->second)->empty());
2291 : // // should have been emptied before calling this
2292 : //
2293 : // delete a->second;
2294 : // }
2295 : //
2296 : // vbAveragers_p.clear();
2297 : //
2298 : // seeIfCubeColumnsExist (vi);
2299 : //}
2300 : //
2301 : //void
2302 : //VbSet::seeIfCubeColumnsExist (ViImplementation2 * vi)
2303 : //{
2304 : // if (vi != 0){
2305 : //
2306 : // // See if the new MS has corrected and/or model data columns
2307 : //
2308 : // doingObservedData_p = averagingOptions_p.contains (AveragingOptions::AverageObserved);
2309 : // doingCorrectedData_p = vi->existsColumn (VisibilityCubeCorrected) &&
2310 : // averagingOptions_p.contains (AveragingOptions::AverageCorrected);
2311 : // doingModelData_p = vi->existsColumn (VisibilityCubeModel) &&
2312 : // averagingOptions_p.contains (AveragingOptions::AverageModel);
2313 : // doingWeightSpectrum_p = vi->existsColumn (WeightSpectrum);
2314 : //
2315 : // // If the use of corrected weights were specified for one of the averages, it's an
2316 : // // error if the column does not exist. Also set the doing flag for corrected weights
2317 : // // if it's being used in some way.
2318 : //
2319 : // Bool needCorrectedWeights =
2320 : // averagingOptions_p.contains (AveragingOptions::ModelUseCorrectedWeights) ||
2321 : // averagingOptions_p.contains (AveragingOptions::CorrectedUseCorrectedWeights);
2322 : //
2323 : // Bool correctedWeightsMissing = needCorrectedWeights &&
2324 : // ! vi->existsColumn (sigmaSpectrum);
2325 : //
2326 : // ThrowIf (correctedWeightsMissing,
2327 : // "Corrected_weight_spectrum not present but required by provided averaging options");
2328 : //
2329 : // doingsigmaSpectrum_p = needCorrectedWeights;
2330 : // }
2331 : //}
2332 : //
2333 : //Int
2334 : //VbSet::getFirstReadyDdid () const
2335 : //{
2336 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2337 : // a != vbAveragers_p.end();
2338 : // a ++){
2339 : //
2340 : // if (a->second->isComplete()){
2341 : // return a->first;
2342 : // }
2343 : // }
2344 : //
2345 : // Assert (false); // ought to have been one that's ready
2346 : //
2347 : // return -1; // shouldn't be called
2348 : //}
2349 : //
2350 : ////void
2351 : ////VbSet::transferAverage (Int ddId, VisBuffer2 * vb)
2352 : ////{
2353 : //// Assert (utilj::containsKey (ddId, vbAveragers_p));
2354 : ////
2355 : //// VbAvg * vba = vbAveragers_p [ddId];
2356 : ////
2357 : //// Assert (vba != 0 && ! vba->empty ());
2358 : ////
2359 : //// // Copy< the completed average into the provided VisBuffer, but
2360 : //// // first reshape the VB if it's shape is different.
2361 : ////
2362 : //// vba->transferAverage (vb);
2363 : ////
2364 : ////}
2365 : //
2366 : //void
2367 : //VbSet::zero ()
2368 : //{
2369 : // for (Averagers::const_iterator a = vbAveragers_p.begin();
2370 : // a != vbAveragers_p.end();
2371 : // a ++){
2372 : //
2373 : // a->second->markEmpty();
2374 : // }
2375 : //}
2376 :
2377 : ///////////////////////
2378 : // //
2379 : // End Namespace avg //
2380 : // //
2381 : ///////////////////////
2382 :
2383 : } // end avg
2384 :
2385 : using namespace avg;
2386 :
2387 363 : AveragingTvi2::AveragingTvi2 (ViImplementation2 * inputVi,
2388 363 : const AveragingParameters & averagingParameters)
2389 : : TransformingVi2 (inputVi),
2390 726 : averagingInterval_p (averagingParameters.getAveragingInterval()),
2391 363 : averagingOptions_p (averagingParameters.getOptions()),
2392 363 : averagingParameters_p (averagingParameters),
2393 363 : ddidLastUsed_p (-1),
2394 363 : inputViiAdvanced_p (false),
2395 726 : vbAvg_p (new VbAvg (averagingParameters, this))
2396 : {
2397 363 : validateInputVi (inputVi);
2398 :
2399 : // Position input Vi to the first subchunk
2400 :
2401 363 : getVii()->originChunks();
2402 363 : getVii()->origin ();
2403 :
2404 363 : setVisBuffer (createAttachedVisBuffer (VbNoOptions));
2405 363 : }
2406 :
2407 726 : AveragingTvi2::~AveragingTvi2 ()
2408 : {
2409 363 : delete vbAvg_p;
2410 726 : }
2411 :
2412 : //void
2413 : //AveragingTvi2::advanceInputVii ()
2414 : //{
2415 : // Assert (false);
2416 : //
2417 : // // Attempt to advance to the next subchunk
2418 : //
2419 : // getVii()->next ();
2420 : //
2421 : // if (! getVii()->more()){
2422 : //
2423 : // // No more subchunks so advance to the next chunk
2424 : //
2425 : // getVii()->nextChunk();
2426 : //
2427 : // if (! getVii()->moreChunks()){
2428 : // return; // no more chunks
2429 : // }
2430 : //
2431 : // // Position to the first subchunk
2432 : //
2433 : // getVii()->origin();
2434 : // }
2435 : //}
2436 :
2437 : //Int
2438 : //AveragingTvi2::determineDdidToUse() const
2439 : //{
2440 : // if (ddidLastUsed_p >= 0 && vbSet_p->anyAveragesReady (ddidLastUsed_p)){
2441 : // return ddidLastUsed_p;
2442 : // }
2443 : //
2444 : // return vbSet_p->getFirstReadyDdid();
2445 : //}
2446 :
2447 : Bool
2448 56208 : AveragingTvi2::more () const
2449 : {
2450 56208 : return more_p;
2451 : }
2452 :
2453 : Bool
2454 25285 : AveragingTvi2::moreChunks () const
2455 : {
2456 25285 : return getVii()->moreChunks();
2457 : }
2458 :
2459 : void
2460 37660 : AveragingTvi2::next ()
2461 : {
2462 37660 : subchunkExists_p = false;
2463 :
2464 37660 : startBuffer_p = endBuffer_p + 1;
2465 37660 : endBuffer_p = startBuffer_p - 1;
2466 :
2467 37660 : if (getVii()->more()){
2468 24927 : getVii()->next();
2469 24927 : endBuffer_p++;
2470 : }
2471 :
2472 37660 : produceSubchunk ();
2473 :
2474 37660 : subchunk_p.incrementSubChunk();
2475 37660 : }
2476 :
2477 : void
2478 12023 : AveragingTvi2::nextChunk ()
2479 : {
2480 : // New chunk, so toss all of the accumulators
2481 :
2482 12023 : vbAvg_p->startChunk (getVii());
2483 :
2484 : // Advance the input to the next chunk as well.
2485 :
2486 12023 : getVii()->nextChunk ();
2487 :
2488 12023 : subchunk_p.incrementChunk();
2489 :
2490 12023 : more_p = false;
2491 12023 : }
2492 :
2493 : void
2494 14691 : AveragingTvi2::origin ()
2495 : {
2496 : // Position input VI to the start of the chunk
2497 :
2498 14691 : subchunk_p.resetSubChunk();
2499 :
2500 14691 : getVii()->origin();
2501 :
2502 14691 : startBuffer_p = 0;
2503 14691 : endBuffer_p = -1;
2504 :
2505 : // Get the first subchunk ready.
2506 :
2507 14691 : produceSubchunk ();
2508 14691 : }
2509 :
2510 : void
2511 507 : AveragingTvi2::originChunks (Bool forceRewind)
2512 : {
2513 : // Ensure that the underlying VI is in a state where some metadata
2514 : // can be retrieved
2515 507 : getVii()->originChunks(forceRewind);
2516 :
2517 : // Initialize the chunk
2518 507 : vbAvg_p->startChunk (getVii());
2519 :
2520 507 : more_p = false;
2521 :
2522 507 : subchunk_p.resetToOrigin();
2523 507 : }
2524 :
2525 : /**
2526 : * Configure a VisBuffer with given averagingOptions related to phase shifting
2527 : *
2528 : * @param vb a VisBuffer to set up in terms of phase shifting
2529 : * @param averagingOptions averaging options enabled
2530 : * @param avgPars AveragingParmeters object to set into the buffer
2531 : */
2532 : void
2533 592578 : setPhaseShiftingOptions(VisBuffer2 * vb, const AveragingOptions &averagingOptions,
2534 : const AveragingParameters &avgPars)
2535 : {
2536 592578 : if (averagingOptions.contains (AveragingOptions::phaseShifting))
2537 : {
2538 0 : if (averagingOptions.contains (AveragingOptions::AverageObserved))
2539 : {
2540 0 : vb->visCube();
2541 : }
2542 :
2543 0 : if (averagingOptions.contains (AveragingOptions::AverageCorrected))
2544 : {
2545 0 : vb->visCubeCorrected();
2546 : }
2547 :
2548 0 : if (averagingOptions.contains (AveragingOptions::AverageModel))
2549 : {
2550 0 : vb->visCubeModel();
2551 : }
2552 :
2553 0 : vb->phaseCenterShift(avgPars.getXpcOffset(),
2554 : avgPars.getYpcOffset());
2555 : }
2556 592578 : }
2557 :
2558 : /**
2559 : * The iteration of this method is different depending on whether "row blocking" is used or
2560 : * not. The reason is that the while loop already had enough complexity embedded when fixes
2561 : * were done to get flagdata+time-average+row-blocking working (CAS-11910). Hopefully in the
2562 : * near future we can get rid of the hacky "row blocking" feature. For the time being, it is
2563 : * not clear how it could possibly work together with the "uvwdistance" feature. So better
2564 : * to keep those features separate.
2565 : *
2566 : * So the "if (block > 0)" separates iteration when using row blocking. That implies that
2567 : * row blocking takes precedence over (and disables) other features like
2568 : * "isUsingUvwDistance()".
2569 : * An alternative would be to add comparisons between block and vbToFill->appendSize() in
2570 : * the ifs below. Something like:
2571 : * if (! vbAvg_p->isUsingUvwDistance()
2572 : * && (block == 0 && vbToFill->appendSize() > 0
2573 : * || (block > 0 && vbToFill->appendSize() >= block)
2574 : * )
2575 : * ){
2576 : * ...
2577 : * else if ((block > 0 && vbToFill->appendSize() < block) ||
2578 : * vbToFill->appendSize() < nBaselines * nWindows){
2579 : * ...
2580 : * } else {
2581 : *
2582 : * But I prefer not adding more complexity to those ifs. The potential combinations would
2583 : * be too many to handle in a handful of if branches, and they are not well understood let
2584 : * alone well tested.
2585 : */
2586 : void
2587 52351 : AveragingTvi2::produceSubchunk ()
2588 : {
2589 52351 : VisBufferImpl2 * vbToFill = dynamic_cast<VisBufferImpl2 *> (getVisBuffer());
2590 52351 : assert (vbToFill != 0);
2591 :
2592 52351 : vbToFill->setFillable (true); // New subchunk, so it's fillable
2593 :
2594 52351 : vbAvg_p->setBufferToFill (vbToFill);
2595 :
2596 52351 : Int nBaselines = nAntennas() * (nAntennas() -1) / 2;
2597 : // This is just a heuristic to keep output VBs from being too small
2598 :
2599 : // jagonzal: Handle nBaselines for SD case
2600 52351 : if (nBaselines == 0) nBaselines = 1;
2601 :
2602 52351 : auto block = getVii()->getRowBlocking();
2603 620002 : while (getVii()->more()){
2604 :
2605 592578 : VisBuffer2 * vb = getVii()->getVisBuffer();
2606 :
2607 592578 : setPhaseShiftingOptions(vb, averagingOptions_p, averagingParameters_p);
2608 :
2609 592578 : bool rowBlocking = block > 0;
2610 592578 : vbAvg_p->accumulate (vb, subchunk_p);
2611 :
2612 592578 : if (rowBlocking) {
2613 77 : auto app = vbToFill->appendSize();
2614 77 : if (app <= block) {
2615 77 : getVii()->next();
2616 77 : endBuffer_p++;
2617 : } else {
2618 0 : break;
2619 : }
2620 : } else {
2621 592501 : Int nWindows = vbAvg_p->nSpectralWindowsInBuffer ();
2622 592501 : if (! vbAvg_p->isUsingUvwDistance() && vbToFill->appendSize() > 0){
2623 : // Doing straight average and some data has been produced so
2624 : // output it to the user
2625 24791 : break;
2626 : }
2627 567710 : else if (vbToFill->appendSize() < nBaselines * nWindows) {
2628 567574 : getVii()->next();
2629 567574 : endBuffer_p++;
2630 : }
2631 : else {
2632 136 : break;
2633 : }
2634 : }
2635 : };
2636 :
2637 52351 : if (! getVii()->more()){
2638 27424 : vbAvg_p->finalizeAverages ();
2639 : }
2640 :
2641 52351 : vbAvg_p->finalizeBufferFilling ();
2642 :
2643 79775 : more_p = getVii()->more() || // more to read
2644 27424 : vbToFill->nRows() > 0; // some to process
2645 52351 : }
2646 :
2647 : // Bool
2648 : // AveragingTvi2::reachedAveragingBoundary()
2649 : // {
2650 : // // An average can be terminated for a variety of reasons:
2651 : // // o the time interval has elapsed
2652 : // // o the current MS is completed
2653 : // // o no more input data
2654 : // // o other (e.g, change of scan, etc.)
2655 :
2656 : // Bool reachedIt = false;
2657 : // VisBuffer2 * vb = getVii()->getVisBuffer();
2658 :
2659 : // if (! getVii()->more() && ! getVii ()->moreChunks()){
2660 :
2661 : // reachedIt = true; // no more input data
2662 : // }
2663 : // else if (vb->isNewMs()){
2664 : // reachedIt = true; // new MS
2665 : // }
2666 :
2667 : // return reachedIt;
2668 : // }
2669 :
2670 : //Bool
2671 : //AveragingTvi2::subchunksReady() const
2672 : //{
2673 : // Bool ready = vbSet_p->anyAveragesReady();
2674 : //
2675 : // return ready;
2676 : //}
2677 :
2678 : void
2679 363 : AveragingTvi2::validateInputVi (ViImplementation2 *)
2680 : {
2681 : // Validate that the input VI is compatible with this VI.
2682 :
2683 : // No implmented right now
2684 363 : }
2685 :
2686 402248632 : Float AveragingTvi2::weightToSigma (Float weight)
2687 : {
2688 402248632 : return weight > FLT_MIN ? 1.0 / std::sqrt (weight) : -1.0; // bogosity indicator
2689 : }
2690 :
2691 : Vector<Float>
2692 5333923 : AveragingTvi2::average (const Matrix<Float> &data, const Matrix<Bool> &flags)
2693 : {
2694 5333923 : IPosition shape = data.shape();
2695 5333923 : Vector<Float> result(shape(0),0);
2696 5333923 : Vector<uInt> samples(shape(0),0);
2697 5333923 : uInt nCorrelations = shape (0);
2698 5333923 : uInt nChannels = shape (1);
2699 :
2700 20680899 : for (uInt correlation = 0; correlation < nCorrelations; correlation++)
2701 : {
2702 15346976 : int nSamples = 0;
2703 15346976 : float sum = 0;
2704 15346976 : bool accumulatorFlag = true;
2705 :
2706 654318022 : for (uInt channel=0; channel< nChannels; channel++)
2707 : {
2708 638971046 : Bool inputFlag = flags(correlation,channel);
2709 : // true/true or false/false
2710 638971046 : if (accumulatorFlag == inputFlag)
2711 : {
2712 611577896 : nSamples ++;
2713 611577896 : sum += data (correlation, channel);
2714 : }
2715 : // true/false: Reset accumulation when accumulator switches from flagged to unflagged
2716 27393150 : else if ( accumulatorFlag and ! inputFlag )
2717 : {
2718 14125890 : accumulatorFlag = false;
2719 14125890 : nSamples = 1;
2720 14125890 : sum = data (correlation, channel);
2721 : }
2722 : // else ignore case where accumulator is valid and data is not
2723 :
2724 : }
2725 :
2726 15346976 : result (correlation) = sum / (nSamples != 0 ? nSamples : 1);
2727 : }
2728 :
2729 10667846 : return result;
2730 5333923 : }
2731 :
2732 : Matrix<Float>
2733 300114 : AveragingTvi2::average (const Cube<Float> &data, const Cube<Bool> &flags)
2734 : {
2735 300114 : IPosition shape = data.shape();
2736 300114 : uInt nRows = shape(2);
2737 300114 : uInt nCorrelations = shape (0);
2738 :
2739 300114 : Matrix<Float> result(nCorrelations, nRows, 0);
2740 :
2741 5115738 : for (uInt row=0; row < nRows; row++)
2742 : {
2743 4815624 : result.column(row) = AveragingTvi2::average (data.xyPlane(row), flags.xyPlane(row));
2744 : }
2745 :
2746 600228 : return result;
2747 300114 : }
2748 :
2749 : /**
2750 : * Strategy to support different ways of propagating flags from the 'transformed' cube to
2751 : * the original ('propagated') cube. Iterates through rows, channels, correlations.
2752 : *
2753 : * This is meant to be used from propagateChanAvgFlags with at least two alternative
2754 : * functors. One to propagate flags as required by flagdata (preserve pre-existing flags
2755 : * in the original data cube), and a second one to propagate flags as required by plotms.
2756 : * CAS-12737, CAS-9985, CAS-12205.
2757 : *
2758 : * @param flagrow per row FLAG_ROW value
2759 : * @param flagMapped propagated FLAG_ROW
2760 : * @param flagCubeMapped Cube of flags in which flags are to be propagated
2761 : * @param row2AvgRow map of input_buffer_row_index->output_buffer_row_index (this is pre-
2762 : * calculated by the "VbAvg" when averaging rows and is needed here).
2763 : */
2764 : template <typename Functor>
2765 529 : void cubePropagateFlags(const Vector<Bool> &flagRow,
2766 : Vector<Bool> &flagMapped,
2767 : Cube<Bool> &flagCubeMapped,
2768 : std::vector<size_t> row2AvgRow,
2769 : Functor propagate)
2770 : {
2771 529 : Int nRowsMapped = flagCubeMapped.shape()[2];
2772 17169 : for(Int rowMapped=0; rowMapped < nRowsMapped; ++rowMapped) {
2773 16640 : size_t index = row2AvgRow[rowMapped];
2774 16640 : flagMapped(rowMapped) = flagRow(index);
2775 :
2776 1081600 : for (Int chan_i=0; chan_i < flagCubeMapped.shape()[1]; ++chan_i) {
2777 5324800 : for (Int corr_i=0; corr_i<flagCubeMapped.shape()[0]; ++corr_i) {
2778 4259840 : propagate(rowMapped, chan_i, corr_i, index);
2779 : }
2780 : }
2781 : }
2782 529 : }
2783 :
2784 : // -----------------------------------------------------------------------
2785 : //
2786 : // -----------------------------------------------------------------------
2787 208 : void AveragingTvi2::writeFlag (const Cube<Bool> & flag)
2788 : {
2789 : // Calculate FLAG_ROW from flag
2790 208 : Vector<Bool> flagRow;
2791 208 : TransformingVi2::calculateFlagRowFromFlagCube(flag, flagRow);
2792 :
2793 : const auto flagdataFlagPropagation =
2794 208 : averagingOptions_p.contains(AveragingOptions::flagdataFlagPropagation);
2795 :
2796 : // Propagate transformed flags
2797 208 : getVii()->origin();
2798 208 : Int currentBuffer = 0;
2799 7285 : while (getVii()->more())
2800 : {
2801 7285 : if ((currentBuffer >= startBuffer_p) and (currentBuffer <= endBuffer_p))
2802 : {
2803 : // Allocated propagated flag vector/cube
2804 529 : uInt nOriginalRows = getVii()->getVisBuffer()->nRows();
2805 529 : Vector<Bool> flagMapped(nOriginalRows,false);
2806 529 : Cube<Bool> flagCubeMapped;
2807 529 : flagCubeMapped = getVii()->getVisBuffer()->flagCube();
2808 :
2809 : // Keeping two separate blocks for 'flagdataFlagPropagation' (CAS-12737,
2810 : // CAS-12205, CAS-9985) until this issue is better settled.
2811 : // Fill propagated flag vector/cube
2812 529 : if (flagdataFlagPropagation)
2813 : {
2814 529 : cubePropagateFlags(flagRow, flagMapped, flagCubeMapped, vbAvg_p->row2AvgRow(),
2815 4259840 : [&flag, &flagCubeMapped]
2816 4902708 : (uInt rowMapped, uInt chan_i, uInt corr_i, uInt index) {
2817 4259840 : if (flag(corr_i,chan_i,index))
2818 642868 : flagCubeMapped(corr_i,chan_i,rowMapped) = true;
2819 4259840 : });
2820 : }
2821 : else
2822 : {
2823 0 : cubePropagateFlags(flagRow, flagMapped, flagCubeMapped, vbAvg_p->row2AvgRow(),
2824 0 : [&flag, &flagCubeMapped]
2825 0 : (uInt rowMapped, uInt chan_i, uInt corr_i, uInt index) {
2826 0 : flagCubeMapped(corr_i, chan_i, rowMapped) =
2827 0 : flag(corr_i, chan_i, index);
2828 0 : });
2829 : }
2830 :
2831 : // Write propagated flag vector/cube
2832 529 : getVii()->writeFlag(flagCubeMapped);
2833 529 : getVii()->writeFlagRow(flagMapped);
2834 529 : }
2835 :
2836 7285 : currentBuffer++;
2837 7285 : getVii()->next();
2838 7285 : if (currentBuffer > endBuffer_p) break;
2839 : }
2840 208 : }
2841 :
2842 : // -----------------------------------------------------------------------
2843 : //
2844 : // -----------------------------------------------------------------------
2845 14 : void AveragingTvi2::writeFlagRow (const Vector<Bool> & flagRow)
2846 : {
2847 : // Create index map for averaged data
2848 14 : VisBuffer2 *avgVB = getVisBuffer();
2849 14 : Vector<Int> avgAnt1 = avgVB->antenna1();
2850 14 : Vector<Int> avgAnt2 = avgVB->antenna2();
2851 14 : Vector<Int> avgSPW = avgVB->spectralWindows();
2852 :
2853 14 : std::map< Int, std::map <Int, std::map< Int, uInt> > > spwAnt1Ant2IndexMap;
2854 578 : for (uInt avgRow=0;avgRow<avgAnt1.size();avgRow++)
2855 : {
2856 564 : spwAnt1Ant2IndexMap[avgSPW(avgRow)][avgAnt1(avgRow)][avgAnt2(avgRow)] = avgRow;
2857 : }
2858 :
2859 : // Propagate transformed flags
2860 14 : getVii()->origin();
2861 14 : Int currentBuffer = 0;
2862 110 : while (getVii()->more())
2863 : {
2864 110 : if ((currentBuffer >= startBuffer_p) and (currentBuffer <= endBuffer_p))
2865 : {
2866 : // Allocated propagated flag vector/cube
2867 110 : uInt nOriginalRows = getVii()->getVisBuffer()->nRows();
2868 110 : Vector<Bool> flagMapped(nOriginalRows,false);
2869 :
2870 : // Get original ant1/ant2/spw cols. to determine the mapped index
2871 110 : Vector<Int> orgAnt1 = getVii()->getVisBuffer()->antenna1();
2872 110 : Vector<Int> orgAnt2 = getVii()->getVisBuffer()->antenna2();
2873 110 : Vector<Int> orgSPW = getVii()->getVisBuffer()->spectralWindows();
2874 :
2875 : // Fill propagated flag vector/cube
2876 1474 : for (uInt row=0;row<nOriginalRows;row++)
2877 : {
2878 1364 : uInt index = spwAnt1Ant2IndexMap[orgSPW(row)][orgAnt1(row)][orgAnt2(row)];
2879 1364 : flagMapped(row) = flagRow(index);
2880 : }
2881 :
2882 : // Write propagated flag vector/cube
2883 110 : getVii()->writeFlagRow(flagMapped);
2884 :
2885 110 : }
2886 :
2887 110 : currentBuffer += 1;
2888 110 : getVii()->next();
2889 110 : if (currentBuffer > endBuffer_p) break;
2890 : }
2891 :
2892 28 : return;
2893 14 : }
2894 :
2895 0 : void AveragingTvi2::visibilityObserved(casacore::Cube<casacore::Complex>& vis) const
2896 : {
2897 0 : if(!averagingOptions_p.contains(AveragingOptions::AverageObserved))
2898 0 : throw AipsError("Requesting visibilityCorrected but column has not been averaged");
2899 0 : VisBuffer2* vb = getVisBuffer();
2900 0 : vis = vb->visCube();
2901 0 : return;
2902 : }
2903 :
2904 0 : void AveragingTvi2::visibilityCorrected(casacore::Cube<casacore::Complex>& vis) const
2905 : {
2906 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeCorrected) ||
2907 0 : !averagingOptions_p.contains(AveragingOptions::AverageCorrected))
2908 0 : throw AipsError("Requesting visibilityCorrected but column not "
2909 0 : "provided by underlying VI or column has not been averaged");
2910 0 : VisBuffer2* vb = getVisBuffer();
2911 0 : vis = vb->visCubeCorrected();
2912 0 : return;
2913 : }
2914 :
2915 0 : void AveragingTvi2::visibilityModel(casacore::Cube<casacore::Complex>& vis) const
2916 : {
2917 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeModel) ||
2918 0 : !averagingOptions_p.contains(AveragingOptions::AverageModel))
2919 0 : throw AipsError("Requesting visibilityModel but column not "
2920 0 : "provided by underlying VI or column has not been averaged");
2921 0 : VisBuffer2* vb = getVisBuffer();
2922 0 : vis = vb->visCubeModel();
2923 0 : return;
2924 : }
2925 :
2926 0 : void AveragingTvi2::floatData(casacore::Cube<casacore::Float>& fcube) const
2927 : {
2928 0 : if(!getVii()->getVisBuffer()->existsColumn (VisBufferComponent2::VisibilityCubeFloat) ||
2929 0 : !averagingOptions_p.contains(AveragingOptions::AverageFloat))
2930 0 : throw AipsError("Requesting floatData but column not "
2931 0 : "provided by underlying VI or column has not been averaged");
2932 0 : VisBuffer2* vb = getVisBuffer();
2933 0 : fcube = vb->visCubeFloat();
2934 0 : return;
2935 : }
2936 :
2937 0 : void AveragingTvi2::flag(casacore::Cube<casacore::Bool>& flags) const
2938 : {
2939 0 : VisBuffer2* vb = getVisBuffer();
2940 0 : flags = vb->flagCube();
2941 0 : return;
2942 : }
2943 :
2944 0 : void AveragingTvi2::flagRow(casacore::Vector<casacore::Bool>& rowflags) const
2945 : {
2946 0 : VisBuffer2* vb = getVisBuffer();
2947 0 : rowflags = vb->flagRow();
2948 0 : return;
2949 : }
2950 :
2951 0 : void AveragingTvi2::sigma(casacore::Matrix<casacore::Float> & sigmat) const
2952 : {
2953 0 : VisBuffer2* vb = getVisBuffer();
2954 0 : sigmat = vb->sigma();
2955 0 : return;
2956 : }
2957 :
2958 0 : void AveragingTvi2::weight (casacore::Matrix<casacore::Float> & wtmat) const
2959 : {
2960 0 : VisBuffer2* vb = getVisBuffer();
2961 0 : wtmat = vb->weight();
2962 0 : return;
2963 : }
2964 :
2965 0 : void AveragingTvi2::weightSpectrum (casacore::Cube<casacore::Float> & wtsp) const
2966 : {
2967 0 : VisBuffer2* vb = getVisBuffer();
2968 0 : wtsp = vb->weightSpectrum();
2969 0 : return;
2970 : }
2971 :
2972 0 : void AveragingTvi2::sigmaSpectrum (casacore::Cube<casacore::Float> & sigsp) const
2973 : {
2974 0 : VisBuffer2* vb = getVisBuffer();
2975 0 : sigsp = vb->sigmaSpectrum();
2976 0 : return;
2977 : }
2978 :
2979 0 : void AveragingTvi2::exposure (casacore::Vector<double> & expo) const
2980 : {
2981 0 : VisBuffer2* vb = getVisBuffer();
2982 0 : expo = vb->exposure();
2983 0 : return;
2984 : }
2985 :
2986 0 : void AveragingTvi2::getRowIds (Vector<rownr_t> & rowids) const
2987 : {
2988 0 : VisBuffer2* vb = getVisBuffer();
2989 0 : rowids = vb->rowIds();
2990 0 : return;
2991 : }
2992 :
2993 33 : void AveragingTvi2::time (casacore::Vector<double> & t) const
2994 : {
2995 33 : VisBuffer2* vb = getVisBuffer();
2996 33 : t = vb->time();
2997 33 : return;
2998 : }
2999 :
3000 0 : void AveragingTvi2::timeInterval (casacore::Vector<double> & ti) const
3001 : {
3002 0 : VisBuffer2* vb = getVisBuffer();
3003 0 : ti = vb->timeInterval();
3004 0 : return;
3005 : }
3006 :
3007 0 : void AveragingTvi2::timeCentroid (casacore::Vector<double> & t) const
3008 : {
3009 0 : VisBuffer2* vb = getVisBuffer();
3010 0 : t = vb->timeCentroid();
3011 0 : return;
3012 : }
3013 :
3014 0 : void AveragingTvi2::uvw (casacore::Matrix<double> & uvwmat) const
3015 : {
3016 0 : VisBuffer2* vb = getVisBuffer();
3017 0 : uvwmat = vb->uvw();
3018 0 : return;
3019 : }
3020 :
3021 34 : void AveragingTvi2::antenna1 (casacore::Vector<casacore::Int> & ant1) const
3022 : {
3023 34 : VisBuffer2* vb = getVisBuffer();
3024 34 : ant1 = vb->antenna1();
3025 34 : return;
3026 : }
3027 :
3028 34 : void AveragingTvi2::antenna2 (casacore::Vector<casacore::Int> & ant2) const
3029 : {
3030 34 : VisBuffer2* vb = getVisBuffer();
3031 34 : ant2 = vb->antenna2();
3032 34 : return;
3033 : }
3034 :
3035 0 : casacore::Bool AveragingTvi2::weightSpectrumExists () const
3036 : {
3037 : //According to VbAvg::startChunk code comments,
3038 : //there is always an output weightSpectrum. See also CAS-11559.
3039 0 : return true;
3040 : }
3041 :
3042 0 : casacore::Bool AveragingTvi2::sigmaSpectrumExists () const
3043 : {
3044 : //According to VbAvg::startChunk code comments,
3045 : //there is always an output sigmaSpectrum. See also CAS-11559.
3046 0 : return true;
3047 : }
3048 :
3049 : } // end namespace vi
3050 :
3051 : using namespace casacore;
3052 : } // end namespace casa
|
