Line data Source code
1 : //# VisibilityIterator2.cc: Step through MeasurementEquation by visibility
2 : //# Copyright (C) 1996-2012
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the !GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: VisibilityIterator2.cc,v 19.15 2006/02/01 01:25:14 kgolap Exp $
27 :
28 : #include <tuple>
29 : #include <casacore/casa/Arrays.h>
30 : #include <casacore/casa/BasicSL/Constants.h>
31 : #include <casacore/casa/Containers/Record.h>
32 : #include <casacore/casa/Exceptions.h>
33 : #include <casacore/casa/Quanta/MVTime.h>
34 : #include <casacore/casa/System/AipsrcValue.h>
35 : #include <casacore/casa/Utilities.h>
36 : #include <casacore/ms/MeasurementSets.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/ms/MSSel/MSSelection.h>
39 : #include <casacore/ms/MSSel/MSSpwIndex.h>
40 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
41 : //#include <msvis/MSVis/StokesVector.h>
42 : #include <stdcasa/UtilJ.h>
43 : #include <msvis/MSVis/MeasurementSet2.h>
44 : #include <msvis/MSVis/MSUtil.h>
45 : #include <msvis/MSVis/MSIter2.h>
46 : #include <msvis/MSVis/SpectralWindow.h>
47 : #include <msvis/MSVis/ViFrequencySelection.h>
48 : #include <msvis/MSVis/VisBuffer2.h>
49 : #include <msvis/MSVis/VisBufferComponents2.h>
50 : #include <msvis/MSVis/VisibilityIterator2.h>
51 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
52 : #include <msvis/MSVis/PointingDirectionCache.h>
53 : #include <msvis/MSVis/VisModelDataI.h>
54 : #include <casacore/tables/Tables/ColDescSet.h>
55 : #include <casacore/tables/Tables/ArrayColumn.h>
56 : #include <casacore/tables/DataMan/IncrStManAccessor.h>
57 : #include <casacore/tables/DataMan/StandardStManAccessor.h>
58 : #include <casacore/tables/Tables/TableDesc.h>
59 : #include <casacore/tables/Tables/TableRecord.h>
60 : #include <casacore/tables/DataMan/TiledStManAccessor.h>
61 :
62 : #include <cassert>
63 : #include <algorithm>
64 : #include <limits>
65 : #include <memory>
66 : #include <map>
67 : #include <vector>
68 :
69 : using std::make_pair;
70 : using namespace casacore;
71 : using namespace casa::vi;
72 : using namespace std;
73 :
74 : using namespace casacore;
75 : namespace casa { //# NAMESPACE CASA - BEGIN
76 :
77 : namespace vi {
78 :
79 : Bool
80 0 : operator!=(const Slice & a, const Slice & b)
81 : {
82 0 : Bool result = a.start() != b.start()||
83 0 : a.length() != b.length() ||
84 0 : a.inc() != b.inc();
85 :
86 0 : return result;
87 : }
88 :
89 : // ChannelSubslicer - represents a single selection from a DATA matrix
90 : //
91 : // Typically the ChannelSubslicer uses selects elements from the DATA matrix in
92 : // two dimensions: correlation and channel. Two slice objects are used, the
93 : // first one for correlation and the other for channel. The slice object is a
94 : // triple of starting channel, increment and the number of channels.
95 :
96 : class ChannelSubslicer
97 : {
98 :
99 : public:
100 :
101 : enum {Correlation, Channel};
102 :
103 2143704 : ChannelSubslicer()
104 2143704 : : subslicer_p()
105 2143704 : {}
106 :
107 714568 : ChannelSubslicer(Int n)
108 714568 : : subslicer_p(n)
109 714568 : {}
110 :
111 714568 : ChannelSubslicer(const Vector<Slice> & axis)
112 714568 : : subslicer_p(axis.nelements())
113 : {
114 1429231 : for (uInt i = 0; i < axis.nelements(); i++) {
115 714663 : subslicer_p[i] = axis(i);
116 : }
117 714568 : }
118 :
119 : Bool
120 : operator==(const ChannelSubslicer & other) const
121 : {
122 : if (other.nelements() != nelements()) {
123 : return false;
124 : }
125 :
126 : for (uInt i = 0; i < nelements(); i++) {
127 :
128 : if (!slicesEqual(subslicer_p[i], other.subslicer_p[i])) {
129 : return false;
130 : }
131 : }
132 :
133 : return true;
134 : }
135 :
136 : Bool
137 : operator!=(const ChannelSubslicer & other) const
138 : {
139 : return !(* this == other);
140 : }
141 :
142 : const Slice &
143 13330342 : getSlice(Int i) const
144 : {
145 13330342 : return subslicer_p[i];
146 : }
147 :
148 : Vector<Slice>
149 13200836 : getSlices() const
150 : {
151 13200836 : return Vector<Slice>(subslicer_p);
152 : }
153 :
154 26634096 : size_t nelements() const
155 : {
156 26634096 : return subslicer_p.size();
157 : }
158 :
159 : void
160 725347 : setSlice(Int i, const Slice & slice)
161 : {
162 725347 : subslicer_p[i] = slice;
163 725347 : }
164 :
165 : protected:
166 :
167 : static Bool
168 : slicesEqual(const Slice & a, const Slice & b) {
169 :
170 : return a.start() == b.start() &&
171 : a.length() == b.length() &&
172 : a.inc() == b.inc();
173 :
174 : }
175 :
176 : private:
177 :
178 : vector<Slice> subslicer_p;
179 : };
180 :
181 : // ChannelSlicer - represents the channels selected in a window.
182 : //
183 : // The ChannelSlicer is a collection of channel and correlation selections.
184 : // Each selection, a ChannelSubslicer, specifies a slice for each dimension of
185 : // the data cube: correlation, channel (row selection is not managed by this
186 : // data structure).
187 :
188 : class ChannelSlicer
189 : {
190 :
191 : public:
192 :
193 : typedef vector<ChannelSubslicer> Rep;
194 : typedef Vector<Vector <Slice> > CoreRep;
195 :
196 714568 : ChannelSlicer()
197 714568 : : slicer_p()
198 714568 : {}
199 :
200 714568 : ChannelSlicer(Int nAxes)
201 714568 : : slicer_p(nAxes)
202 714568 : {}
203 :
204 : bool
205 : operator==(const ChannelSlicer & other) const
206 : {
207 :
208 : if (nelements() != other.nelements()) {
209 : return false;
210 : }
211 :
212 :
213 : for (uInt i = 0; i < nelements(); i++) {
214 :
215 : if (slicer_p[i] != other.slicer_p[i]) {
216 : return false;
217 : }
218 :
219 : }
220 :
221 : return true;
222 : }
223 :
224 : void
225 714568 : appendDimension()
226 : {
227 714568 : slicer_p.push_back(ChannelSubslicer());
228 714568 : }
229 :
230 : CoreRep
231 5385197 : getSlicerInCoreRep() const
232 : {
233 : // Convert to Vector<Vector <Slice> > for use by
234 : // casacore methods
235 :
236 5385197 : CoreRep rep(nelements());
237 :
238 16155591 : for (uInt i = 0; i < nelements(); i++) {
239 :
240 10770394 : const ChannelSubslicer & subslicer = slicer_p[i];
241 :
242 10770394 : rep[i] = subslicer.getSlices();
243 :
244 : }
245 :
246 5385197 : return rep;
247 : }
248 :
249 : ColumnSlicer
250 5385197 : getColumnSlicer() const
251 : {
252 :
253 : // The goal is to create a 2D array of Slicers. Two types of
254 : // Slicers are required: one to slice into the data as read and
255 : // another to put it into the target array. Both are computed and
256 : // returned here. We are effectively computing a cross-product
257 : // between the correlation and channel slices.
258 :
259 : // might be able to replace this?
260 10770394 : CoreRep slicing = getSlicerInCoreRep();
261 :
262 10770394 : Vector<Slice> correlationSlices = slicing(0);
263 10770394 : Vector<Slice> channelSlices = slicing(1);
264 :
265 10770394 : IPosition start(2), length(2), increment(2);
266 :
267 5385197 : uInt nCorrelationSlices = slicing(0).size();
268 5385197 : uInt nChannelSlices = channelSlices.size();
269 5385197 : uInt nSlices = nCorrelationSlices * nChannelSlices;
270 :
271 10770394 : Vector<Slicer *> dataSlices(nSlices, 0);
272 10770394 : Vector<Slicer *> destinationSlices(nSlices, 0);
273 10770394 : IPosition shape(2, 0);
274 :
275 5385197 : uInt outputSlice = 0;
276 5385197 : const uInt Channel = 1;
277 5385197 : const uInt Correlation = 0;
278 :
279 5385197 : uInt correlationDestination = 0;
280 :
281 5385197 : for (uInt correlationSlice = 0;
282 10774402 : correlationSlice < nCorrelationSlices;
283 : correlationSlice++) {
284 :
285 5389205 : const Slice & aSlice = correlationSlices(correlationSlice);
286 :
287 :
288 5389205 : uInt channelDestination = 0;
289 :
290 5389205 : for (uInt channelSlice = 0;
291 10809885 : channelSlice < nChannelSlices; channelSlice++) {
292 :
293 5420680 : const Slice & bSlice = channelSlices(channelSlice);
294 :
295 5420680 : start(Channel) = bSlice.start();
296 5420680 : length(Channel) = bSlice.length();
297 5420680 : increment(Channel) = bSlice.inc();
298 :
299 : // Can't move outside loop because of of mutation during
300 : // destination logic below
301 5420680 : start(Correlation) = aSlice.start();
302 5420680 : length(Correlation) = aSlice.length();
303 5420680 : increment(Correlation) = aSlice.inc();
304 :
305 5420680 : dataSlices[outputSlice] =
306 5420680 : new Slicer(start, length, increment);
307 :
308 : // The destination slicer will always have increment one and
309 : // the same length as the data slicer. However, it's
310 : // starting location depends on how many slices (both axes)
311 : // it is away from the origin.
312 :
313 5420680 : start(Channel) = channelDestination;
314 5420680 : increment(Channel) = 1;
315 5420680 : channelDestination += length(Channel);
316 :
317 5420680 : start(Correlation) = correlationDestination;
318 5420680 : increment(Correlation) = 1;
319 :
320 5420680 : destinationSlices[outputSlice] =
321 5420680 : new Slicer(start, length, increment);
322 :
323 5420680 : outputSlice++;
324 5420680 : shape(Channel) = max(shape(Channel), channelDestination);
325 :
326 : }
327 :
328 5389205 : correlationDestination += length(Correlation);
329 5389205 : shape(Correlation) = correlationDestination;
330 : }
331 :
332 10770394 : return ColumnSlicer(shape, dataSlices, destinationSlices);
333 : }
334 :
335 : const ChannelSubslicer &
336 15734196 : getSubslicer(Int i) const
337 : {
338 15734196 : return slicer_p[i];
339 : }
340 :
341 21540788 : size_t nelements() const
342 : {
343 21540788 : return slicer_p.size();
344 : }
345 :
346 : void
347 1429136 : setSubslicer(Int i, const ChannelSubslicer & subslice)
348 : {
349 1429136 : slicer_p[i] = subslice;
350 1429136 : }
351 :
352 :
353 : String
354 : toString() const
355 : {
356 : String result = "{";
357 :
358 : for (Rep::const_iterator i = slicer_p.begin();
359 : i != slicer_p.end();
360 : i++) {
361 :
362 : result += String((i == slicer_p.begin()) ? "" : ", ") + "{ ";
363 :
364 : const ChannelSubslicer & subslicer = * i;
365 :
366 : for (uInt j = 0; j < subslicer.nelements(); j++) {
367 :
368 : const Slice & slice = subslicer.getSlice(j);
369 : result += String::format(
370 : "(st=%d, len=%d, inc=%d)",
371 : slice.start(), slice.length(), slice.inc());
372 : }
373 :
374 : result += " }";
375 : }
376 :
377 : result += " }";
378 :
379 : return result;
380 : }
381 :
382 : private:
383 :
384 : Rep slicer_p;
385 : };
386 :
387 : // ChannelSelector - class that provides the channels that need to be extracted
388 : // from a spectral window at a given time for a given MS.
389 : //
390 :
391 : class ChannelSelector
392 : {
393 :
394 : public:
395 :
396 714568 : ChannelSelector(
397 : Double time,
398 : Int _msId,
399 : Int _spectralWindowId,
400 : Int _polarizationId,
401 : const ChannelSlicer & slicer)
402 714568 : : msId(_msId)
403 : , polarizationId(_polarizationId)
404 : , spectralWindowId(_spectralWindowId)
405 : , timeStamp(time)
406 714568 : , slicer_p(slicer)
407 : {
408 : // Count up the number of frequencies selected
409 :
410 714568 : const ChannelSubslicer & frequencySlicer = slicer_p.getSubslicer(1);
411 :
412 714568 : nFrequencies_p = 0;
413 :
414 1439915 : for (Int i = 0; i < (int) frequencySlicer.nelements(); i++) {
415 :
416 725347 : nFrequencies_p += frequencySlicer.getSlice(i).length();
417 : }
418 :
419 : // Create the slicer for FlagCategory data which can't use the
420 : // normal slicer.
421 :
422 714568 : createFlagCategorySlicer();
423 714568 : }
424 :
425 : Bool
426 : equals(const ChannelSelector & other) const
427 : {
428 : // To be equal the other ChannelSelector must be for the same
429 : // MS and spectral window
430 :
431 : Bool equal = msId == other.msId &&
432 : spectralWindowId == other.spectralWindowId;
433 :
434 : // If the timestamps match, then they're equal
435 :
436 : if (timeStamp == other.timeStamp) {
437 : return true;
438 : }
439 :
440 : // They differed on timestamps, but if they select the same channels
441 : // then they're equivalent.
442 :
443 : equal = equal && slicer_p == other.slicer_p;
444 :
445 : return equal;
446 : }
447 :
448 : void
449 714568 : createFlagCategorySlicer()
450 : {
451 : // The shape of the flag categories column cell is [nC, nF,
452 : //nCategories] whereas every other sliced column cell is[nC, nF].
453 : //Use the normal slicer to create the flag category slicer.
454 :
455 714568 : slicerFlagCategories_p = slicer_p;
456 714568 : slicerFlagCategories_p.appendDimension();
457 :
458 714568 : }
459 :
460 : // Returns a list of channels to be extracted.
461 :
462 : Vector<Int>
463 1413082 : getChannels() const
464 : {
465 :
466 : // create result of appropriate size
467 1413082 : Vector<Int> frequencies(nFrequencies_p);
468 :
469 : // get channel axis of slicer
470 1413082 : const ChannelSubslicer & channelSlices = slicer_p.getSubslicer(1);
471 :
472 : // Iterator over all of the slices contained in the channel portion
473 : // of the slicer. For each slice, use each channel number to fill
474 : // in the appropriate index of the result. The result will be that
475 : // frequencies[x] will contain the actual channel number that
476 : // resulted from the frequency selection process.
477 :
478 1413082 : int k = 0;
479 :
480 2835943 : for (int i = 0; i < (int) channelSlices.nelements(); i++) {
481 :
482 1422861 : const Slice & slice = channelSlices.getSlice(i);
483 1422861 : Int channel = slice.start();
484 1422861 : Int increment = slice.inc();
485 1422861 : Int nChannels = slice.length();
486 :
487 1422861 : assert(k + nChannels - 1 <= nFrequencies_p);
488 :
489 135370341 : for (int j = 0; j < (int) nChannels; j++) {
490 :
491 133947480 : frequencies[k++] = channel;
492 133947480 : channel += increment;
493 : }
494 :
495 : }
496 :
497 1413082 : return frequencies;
498 : }
499 :
500 : Vector<Int>
501 11176104 : getCorrelations() const {
502 :
503 11176104 : const ChannelSubslicer & correlationAxis = slicer_p.getSubslicer(0);
504 :
505 22352208 : vector<Int> correlations;
506 :
507 22358238 : for (uInt i = 0; i < correlationAxis.nelements(); i++) {
508 :
509 11182134 : const Slice & slice = correlationAxis.getSlice(i);
510 :
511 11182134 : for (uInt j = slice.start(), n = 0;
512 45424221 : n < slice.length();
513 34242087 : j += slice.inc()) {
514 :
515 34242087 : correlations.push_back(j);
516 34242087 : n++;
517 :
518 : }
519 : }
520 :
521 11176104 : Vector<Int> result(correlations.size());
522 :
523 45418191 : for (uInt i = 0; i < correlations.size(); i++) {
524 34242087 : result[i] = correlations[i];
525 : }
526 :
527 22352208 : return result;
528 : }
529 :
530 : Int
531 5051174 : getNFrequencies() const
532 : {
533 5051174 : return nFrequencies_p;
534 : }
535 :
536 : // Returns the ChannelSlicer object which contains the actual
537 : // channelselection for the current time, window and MS.
538 :
539 : const ChannelSlicer &
540 7815639 : getSlicer() const
541 : {
542 7815639 : return slicer_p;
543 : }
544 :
545 : const ChannelSlicer &
546 0 : getSlicerForFlagCategories() const
547 : {
548 0 : return slicerFlagCategories_p;
549 : }
550 :
551 : const Int msId;
552 : const Int polarizationId;
553 : const Int spectralWindowId;
554 : const Double timeStamp;
555 :
556 : private:
557 :
558 : Int nFrequencies_p;
559 : ChannelSlicer slicer_p;
560 : ChannelSlicer slicerFlagCategories_p;
561 : };
562 :
563 : class ChannelSelectorCache
564 : {
565 : public:
566 :
567 5434 : ChannelSelectorCache(Int maxEntries = 200)
568 5434 : : frameOfReference_p(FrequencySelection::Unknown)
569 : , maxEntries_p(maxEntries)
570 5434 : , msId_p(-1)
571 5434 : {}
572 :
573 5434 : ~ChannelSelectorCache() {};
574 :
575 714568 : void add(std::shared_ptr<ChannelSelector> entry, Int frameOfReference)
576 : {
577 714568 : if (entry->msId != msId_p
578 714568 : || frameOfReference != frameOfReference_p) {
579 :
580 : // Cache only holds values for a single MS and a single frame of
581 : // reference at a time.
582 :
583 5583 : flush();
584 :
585 5583 : msId_p = entry->msId;
586 5583 : frameOfReference_p = frameOfReference;
587 : }
588 :
589 714568 : if (cache_p.size() >= maxEntries_p) {
590 :
591 : // Boot the first entry out of the cache.
592 : // In most situations would have the
593 : // lowest timestamp and lowest spectral window id.
594 :
595 520925 : cache_p.erase(cache_p.begin());
596 : }
597 :
598 : Double time =
599 714568 : ((frameOfReference_p == FrequencySelection::ByChannel)
600 1416679 : ? -1 // channel selection not function of time
601 702111 : : entry->timeStamp);
602 :
603 : // take ownership of entry
604 714568 : cache_p[Key(time, entry->spectralWindowId)] = entry;
605 714568 : }
606 :
607 : std::shared_ptr<ChannelSelector>
608 256127026 : find(Double time, Int msId, Int frameOfReference,
609 : Int spectralWindowId) const {
610 :
611 256127026 : std::shared_ptr<ChannelSelector> result = nullptr;
612 :
613 256127026 : if (msId == msId_p && frameOfReference == frameOfReference_p) {
614 :
615 256121443 : if (frameOfReference_p == FrequencySelection::ByChannel) {
616 33293018 : time = -1; // channel selection is not function of time
617 : }
618 :
619 256121443 : Cache::const_iterator i = cache_p.find(Key(time, spectralWindowId));
620 :
621 256121443 : if (i != cache_p.end()) {
622 255412458 : result = i->second;
623 : }
624 : }
625 :
626 256127026 : return result;
627 : }
628 :
629 9107 : void flush() {
630 9107 : cache_p.clear();
631 9107 : }
632 :
633 : private:
634 :
635 : typedef pair<Double, Int> Key; // (time, spectralWindowId)
636 : typedef map <Key, std::shared_ptr<ChannelSelector>> Cache;
637 :
638 : Cache cache_p; // the cache itself
639 : Int frameOfReference_p; // cache's frame of reference
640 : const uInt maxEntries_p; // max # of entries to keep in the cache
641 : Int msId_p; // cache's MS ID
642 :
643 : };
644 :
645 : class SpectralWindowChannel
646 : {
647 :
648 : public:
649 :
650 : SpectralWindowChannel() // for use of vector only
651 : : channel_p(-1)
652 : , frequency_p(-1)
653 : , width_p(-1)
654 : {}
655 :
656 996160 : SpectralWindowChannel(Int channel, Double frequency, Double width)
657 996160 : : channel_p(channel)
658 : , frequency_p(frequency)
659 996160 : , width_p(width)
660 996160 : {}
661 :
662 : Bool
663 : operator<(const SpectralWindowChannel other) const
664 : {
665 : return frequency_p < other.frequency_p;
666 : }
667 :
668 : Bool
669 : operator<(Double other) const
670 : {
671 : return frequency_p < other;
672 : }
673 :
674 : Int
675 2414920 : getChannel() const
676 : {
677 2414920 : return channel_p;
678 : }
679 :
680 : Double
681 131290687 : getFrequency() const
682 : {
683 131290687 : return frequency_p;
684 : }
685 :
686 : Double
687 3580117 : getWidth() const
688 : {
689 3580117 : return width_p;
690 : }
691 :
692 : private:
693 :
694 : Int channel_p;
695 : Double frequency_p;
696 : Double width_p;
697 : };
698 :
699 : Bool
700 0 : operator<(Double frequency, const SpectralWindowChannel & swc)
701 : {
702 0 : return frequency < swc.getFrequency();
703 : }
704 :
705 : Bool
706 0 : operator>(Double frequency, const SpectralWindowChannel & swc)
707 : {
708 0 : return frequency > swc.getFrequency();
709 : }
710 :
711 : class SpectralWindowChannels : public std::vector <SpectralWindowChannel>
712 : {
713 : public:
714 :
715 3441 : SpectralWindowChannels(const Vector<Double> & frequencies,
716 : const Vector<Double> widths)
717 3441 : : std::vector <SpectralWindowChannel>()
718 : {
719 3441 : reserve(frequencies.nelements());
720 :
721 : // Use the passed frequencies and widths to fill out the values.
722 : // The first indices of frequency are the channel numbers.
723 :
724 3441 : int channel = 0;
725 6882 : Array<Double>::const_iterator frequency = frequencies.begin();
726 6882 : Array<Double>::const_iterator width = widths.begin();
727 3441 : int nChannels = frequencies.size();
728 :
729 501521 : for (; channel < nChannels; ++channel, ++frequency, ++width) {
730 498080 : push_back(SpectralWindowChannel(channel, * frequency, * width));
731 : }
732 :
733 : // Remember the positions of the highest and lowest entries and
734 : // whether the channels are in order of increasing frequency.
735 :
736 3441 : increasingOrder_p =
737 3441 : begin()->getFrequency() < rbegin()->getFrequency();
738 :
739 3441 : if (increasingOrder_p) {
740 2049 : lowest_p = 0;
741 2049 : highest_p = ((int) size()) - 1;
742 : } else {
743 1392 : lowest_p = ((int) size()) - 1;
744 1392 : highest_p = 0;
745 : }
746 3441 : }
747 :
748 : double
749 127740836 : getFrequency(int channel) const
750 : {
751 :
752 127740836 : ThrowIf(channel < 0 || channel >= (int) size(),
753 : String::format("Channel %d not in range [0,%d])", channel,
754 : size()-1));
755 :
756 127740836 : double result = at(channel).getFrequency();
757 :
758 127740836 : return result;
759 : }
760 :
761 : double
762 37148 : getWidth(int channel) const
763 : {
764 :
765 37148 : ThrowIf(channel < 0 || channel >= (int) size(),
766 : String::format("Channel %d not in range [0,%d])", channel,
767 : size()-1));
768 :
769 37148 : double result = at(channel).getWidth();
770 :
771 37148 : return result;
772 : }
773 :
774 : Slice
775 712821 : getIntersection(double lowerFrequency, double upperFrequency) const
776 : {
777 : bool noIntersection =
778 712821 : (at(lowest_p).getFrequency() - 0.5 * at(lowest_p).getWidth()
779 1425642 : > upperFrequency) ||
780 712821 : (at(highest_p).getFrequency() + 0.5 * at(highest_p).getWidth()
781 712821 : < lowerFrequency);
782 :
783 712821 : if (noIntersection) {
784 0 : return Slice(0,0); // sentinel value
785 : }
786 :
787 : // Search for the left and right ends of the intersection which need
788 : // to run. The above check for intersection should guarantee that
789 : // one will be found so there's no need to check the iterator
790 : // limits.
791 :
792 712821 : int increment = increasingOrder_p ? 1 : -1;
793 :
794 : int left, right;
795 :
796 712821 : for (left = lowest_p;
797 1019133 : at(left).getFrequency() + 0.5 * at(left).getWidth()
798 1019133 : < lowerFrequency;
799 306312 : left += increment) {
800 : }
801 :
802 712821 : for (right = highest_p;
803 1098194 : at(right).getFrequency() - 0.5 * at(right).getWidth()
804 1098194 : > upperFrequency;
805 385373 : right -= increment) {
806 : }
807 :
808 : // Slices can only go from upward so if the channel frequencies are
809 : // in reverse order, swap them (i.e., (128, 1) --> (1, 128))
810 :
811 : int a, b;
812 712821 : if (increasingOrder_p) {
813 580555 : a = at(left).getChannel();
814 580555 : b = at(right).getChannel();
815 : } else {
816 132266 : a = at(right).getChannel();
817 132266 : b = at(left).getChannel();
818 : }
819 :
820 712821 : Slice s(a, b - a + 1); // (start, length)
821 :
822 712821 : return s;
823 : }
824 :
825 : private:
826 :
827 : int highest_p;
828 : bool increasingOrder_p;
829 : int lowest_p;
830 :
831 : };
832 :
833 :
834 : class SpectralWindowChannelsCache
835 : {
836 :
837 : public:
838 :
839 5434 : SpectralWindowChannelsCache()
840 5434 : : msId_p(-1)
841 5434 : , nBytes_p(0)
842 5434 : {}
843 :
844 5434 : ~SpectralWindowChannelsCache()
845 5434 : {
846 5434 : flush();
847 5434 : }
848 :
849 : void
850 3441 : add(const SpectralWindowChannels * entry, Int msId, Int spectralWindowId)
851 : {
852 3441 : if (msId != msId_p) {
853 2357 : flush();
854 2357 : nBytes_p = 0;
855 : }
856 :
857 : // If necessary, insert code here to put an upper limit on the size
858 : // of the cache.
859 :
860 : // Add the entry to the cache
861 :
862 3441 : map_p[spectralWindowId] = entry;
863 3441 : msId_p = msId;
864 3441 : nBytes_p += entry->size() * sizeof(SpectralWindowChannel);
865 3441 : }
866 :
867 : const SpectralWindowChannels *
868 1841987 : find(Int msId, Int spectralWindowId)
869 : {
870 1841987 : const SpectralWindowChannels * result = 0;
871 :
872 1841987 : if (msId == msId_p) {
873 :
874 1839630 : Map::const_iterator i = map_p.find(spectralWindowId);
875 :
876 1839630 : if (i != map_p.end()) {
877 1838546 : result = i->second;
878 : }
879 : }
880 :
881 1841987 : return result;
882 : }
883 :
884 11315 : void flush()
885 : {
886 : // Delete all of the contained SpectralWindowChannels objects
887 :
888 14756 : for (Map::iterator i = map_p.begin(); i != map_p.end(); i++) {
889 3441 : delete (i->second);
890 : }
891 :
892 11315 : map_p.clear();
893 11315 : }
894 :
895 :
896 : private:
897 :
898 : typedef map<Int, const SpectralWindowChannels *> Map;
899 : // spectralWindowId->SpectralWindowChannels *
900 :
901 : Map map_p;
902 : Int msId_p;
903 : Int nBytes_p;
904 :
905 : };
906 :
907 : Subchunk
908 0 : Subchunk::noMoreData()
909 : {
910 0 : Int maxInt = std::numeric_limits<Int>::max();
911 0 : return Subchunk(maxInt, maxInt);
912 : }
913 :
914 : String
915 0 : Subchunk::toString() const
916 : {
917 0 : return String::format("(%d,%d)", first, second);
918 : }
919 :
920 : template <typename T>
921 : void
922 4599557 : VisibilityIteratorImpl2::getColumnRows(const ArrayColumn<T> & column,
923 : Array<T> & array) const
924 : {
925 9199114 : ColumnSlicer columnSlicer =
926 4599557 : channelSelectors_p[0]->getSlicer().getColumnSlicer();
927 :
928 4599557 : column.getColumnCells(rowBounds_p.subchunkRows_p,
929 : columnSlicer,
930 : array,
931 : true);
932 4599557 : }
933 :
934 : template <typename T>
935 : void
936 0 : VisibilityIteratorImpl2::getColumnRows(const ArrayColumn<T> & column,
937 : Vector<Cube<T>> & cubeVector) const
938 : {
939 0 : cubeVector.resize(channelSelectors_p.size());
940 0 : for (size_t iVector=0; iVector < channelSelectors_p.size(); iVector++)
941 : {
942 0 : ColumnSlicer columnSlicer =
943 0 : channelSelectors_p[iVector]->getSlicer().getColumnSlicer();
944 :
945 0 : column.getColumnCells(rowBounds_p.subchunkEqChanSelRows_p[iVector],
946 : columnSlicer,
947 0 : cubeVector[iVector],
948 : true);
949 : }
950 0 : }
951 :
952 : template <typename T>
953 : void
954 3834304 : VisibilityIteratorImpl2::getColumnRowsMatrix(const ArrayColumn<T> & column,
955 : Matrix<T> & array,
956 : Bool correlationSlicing) const
957 : {
958 3834304 : if (correlationSlicing) {
959 :
960 : // Extract the correlations slices and repackage them for use by
961 : // getColumnCells
962 :
963 2430442 : const ChannelSlicer & slicer = channelSelectors_p[0]->getSlicer();
964 : // has to be at least one
965 4860884 : const ChannelSubslicer subslicer = slicer.getSubslicer(0);
966 :
967 4860884 : Vector<Slice> correlationSlices = subslicer.getSlices();
968 :
969 4860884 : Vector<Slicer *> dataSlicers(correlationSlices.size(), 0);
970 4860884 : Vector<Slicer *> destinationSlicers(correlationSlices.size(), 0);
971 :
972 4860884 : IPosition start(1, 0), length(1, 0), increment(1, 0);
973 2430442 : uInt sliceStart = 0;
974 :
975 4861688 : for (uInt i = 0; i < correlationSlices.size(); i++) {
976 :
977 2431246 : start(0) = correlationSlices(i).start();
978 2431246 : length(0) = correlationSlices(i).length();
979 2431246 : increment(0) = correlationSlices(i).inc();
980 2431246 : dataSlicers(i) = new Slicer(start, length, increment);
981 :
982 2431246 : start(0) = sliceStart;
983 2431246 : increment(0) = 1;
984 2431246 : destinationSlicers(i) = new Slicer(start, length, increment);
985 :
986 2431246 : sliceStart += length(0);
987 : }
988 :
989 4860884 : IPosition shape(1, sliceStart);
990 :
991 4860884 : ColumnSlicer columnSlicer(shape, dataSlicers, destinationSlicers);
992 :
993 2430442 : column.getColumnCells(rowBounds_p.subchunkRows_p, columnSlicer, array,
994 : true);
995 : }
996 : else{
997 :
998 1403862 : column.getColumnCells(rowBounds_p.subchunkRows_p, array, true);
999 : }
1000 3834304 : }
1001 :
1002 : template <typename T>
1003 : void
1004 0 : VisibilityIteratorImpl2::getColumnRowsMatrix(const ArrayColumn<T> & column,
1005 : Vector<Matrix<T>> & matrixVector) const
1006 : {
1007 0 : matrixVector.resize(channelSelectors_p.size());
1008 0 : for (size_t iVector=0; iVector < channelSelectors_p.size(); iVector++)
1009 : {
1010 :
1011 : // Extract the correlations slices and repackage them for use by
1012 : // getColumnCells
1013 :
1014 0 : const ChannelSlicer & slicer = channelSelectors_p[iVector]->getSlicer();
1015 : // has to be at least one
1016 0 : const ChannelSubslicer subslicer = slicer.getSubslicer(0);
1017 :
1018 0 : Vector<Slice> correlationSlices = subslicer.getSlices();
1019 :
1020 0 : Vector<Slicer *> dataSlicers(correlationSlices.size(), 0);
1021 0 : Vector<Slicer *> destinationSlicers(correlationSlices.size(), 0);
1022 :
1023 0 : IPosition start(1, 0), length(1, 0), increment(1, 0);
1024 0 : uInt sliceStart = 0;
1025 :
1026 0 : for (uInt i = 0; i < correlationSlices.size(); i++) {
1027 :
1028 0 : start(0) = correlationSlices(i).start();
1029 0 : length(0) = correlationSlices(i).length();
1030 0 : increment(0) = correlationSlices(i).inc();
1031 0 : dataSlicers(i) = new Slicer(start, length, increment);
1032 :
1033 0 : start(0) = sliceStart;
1034 0 : increment(0) = 1;
1035 0 : destinationSlicers(i) = new Slicer(start, length, increment);
1036 :
1037 0 : sliceStart += length(0);
1038 : }
1039 :
1040 0 : IPosition shape(1, sliceStart);
1041 :
1042 0 : ColumnSlicer columnSlicer(shape, dataSlicers, destinationSlicers);
1043 :
1044 0 : column.getColumnCells(rowBounds_p.subchunkEqChanSelRows_p[iVector],
1045 : columnSlicer,
1046 0 : matrixVector[iVector],
1047 : true);
1048 : }
1049 0 : }
1050 :
1051 : template <typename T>
1052 : void
1053 26826866 : VisibilityIteratorImpl2::getColumnRows(const ScalarColumn<T> & column,
1054 : Vector<T> & array) const
1055 : {
1056 26826866 : column.getColumnCells(rowBounds_p.subchunkRows_p, array, true);
1057 26826866 : }
1058 :
1059 : template <typename T>
1060 : void
1061 785640 : VisibilityIteratorImpl2::putColumnRows(
1062 : ArrayColumn<T> & column,
1063 : const Array<T> & array)
1064 : {
1065 1571280 : ColumnSlicer columnSlicer =
1066 785640 : channelSelectors_p[0]->getSlicer().getColumnSlicer();
1067 :
1068 785640 : column.putColumnCells(rowBounds_p.subchunkRows_p,
1069 : columnSlicer,
1070 : array);
1071 785640 : }
1072 :
1073 : template <typename T>
1074 : void
1075 184695 : VisibilityIteratorImpl2::putColumnRows(
1076 : ArrayColumn<T> & column,
1077 : const Matrix<T> & array)
1078 : {
1079 184695 : RefRows & rows = rowBounds_p.subchunkRows_p;
1080 :
1081 184695 : column.putColumnCells(rows, array);
1082 184695 : }
1083 :
1084 : template <typename T>
1085 : void
1086 383243 : VisibilityIteratorImpl2::putColumnRows(
1087 : ScalarColumn<T> & column,
1088 : const Vector <T> & array)
1089 : {
1090 383243 : RefRows & rows = rowBounds_p.subchunkRows_p;
1091 :
1092 383243 : column.putColumnCells(rows, array);
1093 383243 : }
1094 :
1095 5434 : VisibilityIteratorImpl2::VisibilityIteratorImpl2(
1096 : const Block<const MeasurementSet *> &mss,
1097 : const SortColumns & sortColumns,
1098 : Double timeInterval,
1099 : Bool writable,
1100 5434 : Bool useMSIter2)
1101 : : ViImplementation2(),
1102 : channelSelectors_p(),
1103 5434 : channelSelectorCache_p(new ChannelSelectorCache()),
1104 : columns_p(),
1105 : floatDataFound_p(false),
1106 : frequencySelections_p(nullptr),
1107 : measurementFrame_p(VisBuffer2::FrameNotSpecified),
1108 5434 : modelDataGenerator_p(VisModelDataI::create2()),
1109 : more_p(false),
1110 : msIndex_p(0),
1111 : msIterAtOrigin_p(false),
1112 : msIter_p(),
1113 : nCorrelations_p(-1),
1114 : nRowBlocking_p(0),
1115 5434 : pendingChanges_p(new PendingChanges()),
1116 : reportingFrame_p(VisBuffer2::FrameNotSpecified),
1117 5434 : spectralWindowChannelsCache_p(new SpectralWindowChannelsCache()),
1118 : subtableColumns_p(nullptr),
1119 5434 : tileCacheModMtx_p(new std::mutex()),
1120 5434 : tileCacheIsSet_p(new std::vector<bool>()),
1121 : timeInterval_p(timeInterval),
1122 : vb_p(nullptr),
1123 : weightScaling_p(),
1124 : writable_p(writable),
1125 : ddIdScope_p(UnknownScope),
1126 : timeScope_p(UnknownScope),
1127 : sortColumns_p(sortColumns),
1128 38038 : subchunkSortColumns_p(false)
1129 : {
1130 : // Set the default subchunk iteration sorting scheme, i.e.
1131 : // unique timestamps in each subchunk.
1132 5434 : auto singleTimeCompare = std::make_shared<ObjCompare<Double> >();
1133 5434 : subchunkSortColumns_p.addSortingColumn(MS::TIME, singleTimeCompare);
1134 5434 : timeScope_p = SubchunkScope;
1135 :
1136 5434 : initialize(mss, useMSIter2);
1137 :
1138 5434 : VisBufferOptions options = writable_p ? VbWritable : VbNoOptions;
1139 :
1140 5434 : vb_p = createAttachedVisBuffer(options);
1141 5434 : }
1142 :
1143 0 : VisibilityIteratorImpl2::VisibilityIteratorImpl2(
1144 : const casacore::Block<const casacore::MeasurementSet *> & mss,
1145 : const SortColumns & chunkSortColumns,
1146 : const SortColumns & subchunkSortColumns,
1147 0 : bool isWritable)
1148 : : ViImplementation2(),
1149 : channelSelectors_p(),
1150 0 : channelSelectorCache_p(new ChannelSelectorCache()),
1151 : columns_p(),
1152 : floatDataFound_p(false),
1153 : frequencySelections_p(nullptr),
1154 : measurementFrame_p(VisBuffer2::FrameNotSpecified),
1155 0 : modelDataGenerator_p(VisModelDataI::create2()),
1156 : more_p(false),
1157 : msIndex_p(0),
1158 : msIterAtOrigin_p(false),
1159 : msIter_p(),
1160 : nCorrelations_p(-1),
1161 : nRowBlocking_p(0),
1162 0 : pendingChanges_p(new PendingChanges()),
1163 : reportingFrame_p(VisBuffer2::FrameNotSpecified),
1164 0 : spectralWindowChannelsCache_p(new SpectralWindowChannelsCache()),
1165 : subtableColumns_p(nullptr),
1166 0 : tileCacheModMtx_p(new std::mutex()),
1167 0 : tileCacheIsSet_p(new std::vector<bool>()),
1168 : timeInterval_p(0),
1169 : vb_p(nullptr),
1170 : weightScaling_p(),
1171 : writable_p(isWritable),
1172 : ddIdScope_p(UnknownScope),
1173 : timeScope_p(UnknownScope),
1174 : sortColumns_p(chunkSortColumns),
1175 0 : subchunkSortColumns_p(subchunkSortColumns)
1176 : {
1177 0 : initialize(mss);
1178 :
1179 0 : VisBufferOptions options = writable_p ? VbWritable : VbNoOptions;
1180 :
1181 0 : vb_p = createAttachedVisBuffer(options);
1182 0 : }
1183 :
1184 0 : VisibilityIteratorImpl2::VisibilityIteratorImpl2(
1185 0 : const VisibilityIteratorImpl2& vii)
1186 : : ViImplementation2()
1187 : , channelSelectors_p()
1188 : , channelSelectorCache_p(nullptr)
1189 : , frequencySelections_p(nullptr)
1190 : , modelDataGenerator_p(nullptr)
1191 : , pendingChanges_p(nullptr)
1192 : , spectralWindowChannelsCache_p(nullptr)
1193 0 : , vb_p(nullptr)
1194 : {
1195 0 : *this = vii;
1196 0 : }
1197 :
1198 : VisibilityIteratorImpl2 &
1199 0 : VisibilityIteratorImpl2::operator=(const VisibilityIteratorImpl2& vii)
1200 : {
1201 : // clone msIter_p
1202 0 : msIter_p = vii.msIter_p->clone();
1203 :
1204 : // copy cache
1205 0 : cache_p = vii.cache_p;
1206 :
1207 0 : tileCacheModMtx_p = vii.tileCacheModMtx_p;
1208 0 : tileCacheIsSet_p = vii.tileCacheIsSet_p;
1209 :
1210 : // clone frequencySelections_p
1211 0 : if (frequencySelections_p)
1212 0 : delete frequencySelections_p;
1213 0 : frequencySelections_p = vii.frequencySelections_p->clone();
1214 :
1215 : // copy channelSelectors_p.
1216 0 : channelSelectors_p.clear();
1217 0 : for (auto chSel : vii.channelSelectors_p)
1218 : {
1219 0 : channelSelectors_p.push_back(std::shared_ptr<ChannelSelector>
1220 0 : (new ChannelSelector(chSel->timeStamp, chSel->msId,
1221 0 : chSel->spectralWindowId, chSel->polarizationId,
1222 0 : chSel->getSlicer())));
1223 : }
1224 :
1225 : // get frame of reference for current MS
1226 : const FrequencySelection &selection =
1227 0 : vii.frequencySelections_p->get(vii.msId());
1228 0 : Int frameOfReference = selection.getFrameOfReference();
1229 :
1230 : // initialize channelSelectorCache_p with current channelSelectors_p
1231 0 : channelSelectorCache_p->flush();
1232 0 : for (auto chSel : channelSelectors_p)
1233 0 : channelSelectorCache_p->add(chSel, frameOfReference);
1234 :
1235 : // copy assign some values
1236 0 : columns_p = vii.columns_p;
1237 0 : floatDataFound_p = vii.floatDataFound_p;
1238 0 : imwgt_p = vii.imwgt_p;
1239 0 : measurementFrame_p = vii.measurementFrame_p;
1240 0 : more_p = vii.more_p;
1241 0 : msIndex_p = vii.msIndex_p;
1242 0 : msIterAtOrigin_p = vii.msIterAtOrigin_p;
1243 0 : nCorrelations_p = vii.nCorrelations_p;
1244 0 : nRowBlocking_p = vii.nRowBlocking_p;
1245 0 : reportingFrame_p = vii.reportingFrame_p;
1246 0 : rowBounds_p = vii.rowBounds_p;
1247 0 : subchunk_p = vii.subchunk_p;
1248 0 : timeFrameOfReference_p = vii.timeFrameOfReference_p;
1249 0 : timeInterval_p = vii.timeInterval_p;
1250 0 : weightScaling_p = vii.weightScaling_p;
1251 0 : writable_p = vii.writable_p;
1252 0 : ddIdScope_p = vii.ddIdScope_p;
1253 0 : timeScope_p = vii.timeScope_p;
1254 :
1255 : // clone modelDataGenerator_p
1256 0 : if (modelDataGenerator_p) delete modelDataGenerator_p;
1257 0 : modelDataGenerator_p = vii.modelDataGenerator_p->clone();
1258 :
1259 : // initialize MSDerivedValues as done in configureNewChunk()
1260 0 : msd_p.setAntennas(msIter_p->msColumns().antenna());
1261 0 : msd_p.setFieldCenter(msIter_p->phaseCenter());
1262 :
1263 : // clone pendingChanges_p
1264 0 : pendingChanges_p.reset(vii.pendingChanges_p->clone());
1265 :
1266 : // initialize subtableColumns_p
1267 0 : if (subtableColumns_p) delete subtableColumns_p;
1268 0 : subtableColumns_p = new SubtableColumns(msIter_p);
1269 :
1270 : // initialize attached VisBuffer...vii.vb_p does *not* get copied
1271 : // TODO: it would be better to use a shared_ptr to the attached VisBuffer,
1272 : // since we don't know if there are any outstanding references, as they can
1273 : // escape from VisibilityIteratorImpl2...for now, just delete it
1274 0 : if (vb_p) delete vb_p;
1275 0 : vb_p = createAttachedVisBuffer(writable_p ? VbWritable : VbNoOptions);
1276 :
1277 0 : return *this;
1278 : }
1279 :
1280 0 : VisibilityIteratorImpl2::VisibilityIteratorImpl2(VisibilityIteratorImpl2&& vii)
1281 : : ViImplementation2()
1282 : , channelSelectors_p()
1283 : , channelSelectorCache_p(nullptr)
1284 : , frequencySelections_p(nullptr)
1285 : , modelDataGenerator_p(nullptr)
1286 : , pendingChanges_p(nullptr)
1287 : , spectralWindowChannelsCache_p(nullptr)
1288 : , vb_p(nullptr)
1289 0 : , writable_p(false)
1290 : {
1291 0 : *this = std::move(vii);
1292 0 : }
1293 :
1294 : VisibilityIteratorImpl2 &
1295 0 : VisibilityIteratorImpl2::operator=(VisibilityIteratorImpl2&& vii)
1296 : {
1297 : // copy msIter_p
1298 0 : msIter_p = vii.msIter_p;
1299 :
1300 : // copy cache
1301 0 : cache_p = vii.cache_p;
1302 :
1303 0 : tileCacheModMtx_p = std::move(vii.tileCacheModMtx_p);
1304 0 : tileCacheIsSet_p = std::move(vii.tileCacheIsSet_p);
1305 :
1306 : // move frequencySelections_p
1307 0 : if (frequencySelections_p)
1308 0 : delete frequencySelections_p;
1309 0 : frequencySelections_p = vii.frequencySelections_p;
1310 0 : vii.frequencySelections_p = nullptr;
1311 :
1312 : // move channelSelectors_p. Owned by channelSelectorCache_p
1313 : // so don't delete initial destination value or source value
1314 0 : channelSelectors_p = std::move(vii.channelSelectors_p);
1315 :
1316 : // move channelSelectorCache_p
1317 0 : if (channelSelectorCache_p)
1318 0 : delete channelSelectorCache_p;
1319 0 : channelSelectorCache_p = vii.channelSelectorCache_p;
1320 0 : vii.channelSelectorCache_p = nullptr;
1321 :
1322 : // move backWriters_p
1323 0 : backWriters_p = std::move(vii.backWriters_p);
1324 :
1325 : // copy assign some values
1326 0 : autoTileCacheSizing_p = vii.autoTileCacheSizing_p;
1327 0 : columns_p = vii.columns_p;
1328 0 : floatDataFound_p = vii.floatDataFound_p;
1329 0 : imwgt_p = vii.imwgt_p;
1330 0 : measurementFrame_p = vii.measurementFrame_p;
1331 0 : measurementSets_p = vii.measurementSets_p;
1332 0 : more_p = vii.more_p;
1333 0 : msIndex_p = vii.msIndex_p;
1334 0 : msIterAtOrigin_p = vii.msIterAtOrigin_p;
1335 0 : nCorrelations_p = vii.nCorrelations_p;
1336 0 : nRowBlocking_p = vii.nRowBlocking_p;
1337 0 : reportingFrame_p = vii.reportingFrame_p;
1338 0 : rowBounds_p = vii.rowBounds_p;
1339 0 : subchunk_p = vii.subchunk_p;
1340 0 : timeFrameOfReference_p = vii.timeFrameOfReference_p;
1341 0 : timeInterval_p = vii.timeInterval_p;
1342 0 : weightScaling_p = vii.weightScaling_p;
1343 0 : writable_p = vii.writable_p;
1344 0 : ddIdScope_p = vii.ddIdScope_p;
1345 0 : timeScope_p = vii.timeScope_p;
1346 0 : sortColumns_p = vii.sortColumns_p;
1347 0 : subchunkSortColumns_p = vii.subchunkSortColumns_p;
1348 :
1349 : // move modelDataGenerator_p
1350 0 : if (modelDataGenerator_p)
1351 0 : delete modelDataGenerator_p;
1352 0 : modelDataGenerator_p = vii.modelDataGenerator_p;
1353 0 : vii.modelDataGenerator_p = nullptr;
1354 :
1355 : // initialize MSDerivedValues as done in configureNewChunk()...moving or
1356 : // copying the value is not well supported
1357 0 : msd_p.setAntennas(msIter_p->msColumns().antenna());
1358 0 : msd_p.setFieldCenter(msIter_p->phaseCenter());
1359 :
1360 : // move pendingChanges_p
1361 0 : pendingChanges_p = std::move(vii.pendingChanges_p);
1362 :
1363 : // initialize subtableColumns_p
1364 0 : if (subtableColumns_p)
1365 0 : delete subtableColumns_p;
1366 0 : subtableColumns_p = vii.subtableColumns_p;
1367 0 : vii.subtableColumns_p = nullptr;
1368 :
1369 : // initialize vb_p...can't steal VisBuffer2 instance since it is attached to
1370 : // vii
1371 0 : if (vb_p)
1372 0 : delete vb_p;
1373 0 : vb_p = createAttachedVisBuffer(writable_p ? VbWritable : VbNoOptions);
1374 : // TODO: again, it would be better were vb_p a shared_ptr
1375 0 : delete vii.vb_p;
1376 0 : vii.vb_p = nullptr;
1377 :
1378 0 : return *this;
1379 : }
1380 :
1381 : void
1382 5503 : VisibilityIteratorImpl2::addDataSelection(const MeasurementSet & ms)
1383 : {
1384 5503 : const MeasurementSet2 * ms2 = dynamic_cast <const MeasurementSet2 *>(&ms);
1385 :
1386 5503 : if (ms2 == 0) {
1387 :
1388 : // A normal MS was used so there is no frequency selection attached to
1389 : // it. Simply add in an empty selection which will select everything.
1390 :
1391 5503 : frequencySelections_p->add(FrequencySelectionUsingChannels());
1392 :
1393 5503 : return;
1394 :
1395 : }
1396 :
1397 : // Get the channel and correlation selectin.
1398 : //
1399 : // Channel slices are indexed by spectral window ID and correlation slices
1400 : // by polarization ID
1401 :
1402 : MSSelection * msSelection =
1403 0 : const_cast<MSSelection *>(ms2->getMSSelection());
1404 : // *KLUGE* -- MSSelection is somewhat sloppy about making methods const
1405 : // so simply getting the slices requires a non-const object ;-(
1406 :
1407 0 : Vector <Vector <Slice> > channelSlices;
1408 0 : msSelection->getChanSlices(channelSlices, ms2);
1409 0 : Vector <Vector <Slice> > correlationSlices;
1410 0 : msSelection->getCorrSlices(correlationSlices, ms2);
1411 :
1412 0 : FrequencySelectionUsingChannels selection;
1413 :
1414 0 : for (uInt spectralWindow = 0;
1415 0 : spectralWindow < channelSlices.nelements();
1416 : spectralWindow++) {
1417 :
1418 : // Get the frequency slices for this spectral window
1419 :
1420 0 : Vector<Slice> & slices = channelSlices[spectralWindow];
1421 :
1422 0 : for (uInt s = 0; s < slices.nelements(); s++) {
1423 :
1424 : // Add each frequency slice to the selection for this spectral
1425 : // window
1426 :
1427 0 : Slice & slice = slices[s];
1428 :
1429 0 : selection.add(spectralWindow, slice.start(), slice.length(),
1430 0 : slice.inc());
1431 : }
1432 : }
1433 :
1434 0 : selection.addCorrelationSlices(correlationSlices);
1435 :
1436 0 : frequencySelections_p->add(selection);
1437 : }
1438 :
1439 :
1440 : void
1441 5434 : VisibilityIteratorImpl2::initialize(const Block<const MeasurementSet *> &mss,
1442 : Bool useMSIter2)
1443 : {
1444 :
1445 5434 : ThrowIf(!sortColumns_p.usingDefaultSortingFunctions(),
1446 : "Sorting definition for chunks doesn't support generic functions yet");
1447 :
1448 5434 : cache_p.flush();
1449 :
1450 5434 : msIndex_p = 0;
1451 :
1452 5434 : frequencySelections_p = new FrequencySelections();
1453 :
1454 5434 : Int nMs = mss.nelements();
1455 5434 : measurementSets_p.resize(nMs);
1456 5434 : tileCacheIsSet_p->resize(nMs);
1457 :
1458 10937 : for (Int k = 0; k < nMs; ++k) {
1459 5503 : measurementSets_p[k] = * mss[k];
1460 5503 : addDataSelection(measurementSets_p[k]);
1461 5503 : (*tileCacheIsSet_p)[k] = false;
1462 : }
1463 :
1464 5434 : if (useMSIter2)
1465 :
1466 : // This version uses the MSSmartInterval for time comparisons in the
1467 : // Table sort/iteration
1468 167 : msIter_p = new MSIter2(measurementSets_p,
1469 167 : sortColumns_p.getColumnIds(),
1470 : timeInterval_p,
1471 167 : sortColumns_p.shouldAddDefaultColumns(),
1472 167 : false);
1473 : else
1474 : // The old-fashioned version
1475 5267 : msIter_p = new MSIter(measurementSets_p,
1476 5267 : sortColumns_p.getColumnIds(),
1477 : timeInterval_p,
1478 5267 : sortColumns_p.shouldAddDefaultColumns(),
1479 5267 : false);
1480 :
1481 5434 : subtableColumns_p = new SubtableColumns(msIter_p);
1482 :
1483 : // Check whether DDID is unique within each chunk. Otherwise assume
1484 : // that it can change for every row.
1485 5434 : ddIdScope_p = RowScope;
1486 5434 : freqSelScope_p = RowScope;
1487 5434 : if (sortColumns_p.shouldAddDefaultColumns())
1488 3124 : ddIdScope_p = ChunkScope;
1489 : else
1490 : {
1491 13488 : for (auto sortCol : sortColumns_p.getColumnIds())
1492 : {
1493 11178 : if(sortCol == MS::DATA_DESC_ID)
1494 2218 : ddIdScope_p = ChunkScope;
1495 : }
1496 : }
1497 :
1498 5434 : freqSelScope_p = ddIdScope_p;
1499 : // If frequency/channel selection also depends on time (selection based on
1500 : // frequencies), then the scope can be further limited by timestamp scope
1501 5434 : if (!(frequencySelections_p->getFrameOfReference() == FrequencySelection::ByChannel))
1502 : {
1503 0 : if(freqSelScope_p == ChunkScope)
1504 0 : freqSelScope_p = timeScope_p;
1505 : }
1506 :
1507 :
1508 5434 : casacore::AipsrcValue<Bool>::find(
1509 5434 : autoTileCacheSizing_p,
1510 10868 : VisibilityIterator2::getAipsRcBase() + ".AutoTileCacheSizing", false);
1511 5434 : }
1512 :
1513 : void
1514 0 : VisibilityIteratorImpl2::initialize(const Block<const MeasurementSet *> &mss)
1515 : {
1516 0 : cache_p.flush();
1517 :
1518 0 : msIndex_p = 0;
1519 :
1520 0 : frequencySelections_p = new FrequencySelections();
1521 :
1522 0 : Int nMs = mss.nelements();
1523 0 : measurementSets_p.resize(nMs);
1524 0 : tileCacheIsSet_p->resize(nMs);
1525 :
1526 0 : for (Int k = 0; k < nMs; ++k) {
1527 0 : measurementSets_p[k] = * mss[k];
1528 0 : addDataSelection(measurementSets_p[k]);
1529 0 : (*tileCacheIsSet_p)[k] = false;
1530 : }
1531 :
1532 0 : msIter_p = new MSIter(measurementSets_p,
1533 0 : sortColumns_p.sortingDefinition());
1534 :
1535 0 : subtableColumns_p = new SubtableColumns(msIter_p);
1536 :
1537 : // Set the scope of each of the metadata to track
1538 0 : setMetadataScope();
1539 :
1540 0 : casacore::AipsrcValue<Bool>::find(
1541 0 : autoTileCacheSizing_p,
1542 0 : VisibilityIterator2::getAipsRcBase() + ".AutoTileCacheSizing", false);
1543 0 : }
1544 :
1545 10868 : VisibilityIteratorImpl2::~VisibilityIteratorImpl2()
1546 : {
1547 5434 : if (channelSelectorCache_p) delete channelSelectorCache_p;
1548 5434 : if (frequencySelections_p) delete frequencySelections_p;
1549 5434 : if (modelDataGenerator_p) delete modelDataGenerator_p;
1550 5434 : if (spectralWindowChannelsCache_p) delete spectralWindowChannelsCache_p;
1551 5434 : if (subtableColumns_p) delete subtableColumns_p;
1552 5434 : if (vb_p) delete vb_p;
1553 10868 : }
1554 :
1555 : std::unique_ptr<VisibilityIteratorImpl2>
1556 0 : VisibilityIteratorImpl2::clone() const
1557 : {
1558 0 : unsigned nms = measurementSets_p.nelements();
1559 0 : Block<const MeasurementSet *> msps(nms);
1560 0 : for (unsigned i = 0; i < nms; ++i)
1561 0 : msps[i] = &measurementSets_p[i];
1562 :
1563 : std::unique_ptr<VisibilityIteratorImpl2> result(
1564 : new VisibilityIteratorImpl2(
1565 : msps,
1566 0 : sortColumns_p,
1567 0 : timeInterval_p,
1568 0 : writable_p,
1569 0 : false));
1570 0 : *result = *this;
1571 0 : return result;
1572 : }
1573 :
1574 7048 : void VisibilityIteratorImpl2::setMetadataScope()
1575 : {
1576 : // Check if each chunk will receive a unique DDId.
1577 : // For that it must be a sorting column and comparison function should
1578 : // be of the type ObjCompare<Int>
1579 7048 : bool uniqueDDIdInChunk = false;
1580 7048 : bool uniqueDDIdInSubchunk = false;
1581 21712 : for(auto& sortDef : sortColumns_p.sortingDefinition())
1582 18170 : if(sortDef.first == MS::columnName(MS::DATA_DESC_ID) &&
1583 3506 : dynamic_cast<ObjCompare<Int>*>(sortDef.second.get()))
1584 0 : uniqueDDIdInChunk = true;
1585 :
1586 14096 : for(auto& sortDef : subchunkSortColumns_p.sortingDefinition())
1587 7048 : if(sortDef.first == MS::columnName(MS::DATA_DESC_ID) &&
1588 0 : dynamic_cast<ObjCompare<Int>*>(sortDef.second.get()))
1589 0 : uniqueDDIdInSubchunk = true;
1590 :
1591 7048 : if(uniqueDDIdInChunk)
1592 0 : ddIdScope_p = ChunkScope;
1593 7048 : else if(uniqueDDIdInSubchunk)
1594 0 : ddIdScope_p = SubchunkScope;
1595 : else
1596 7048 : ddIdScope_p = RowScope;
1597 :
1598 : // Similar for time
1599 7048 : bool uniqueTimeInChunk = false, uniqueTimeInSubchunk = false;
1600 21712 : for(auto& sortDef : sortColumns_p.sortingDefinition())
1601 17088 : if(sortDef.first == MS::columnName(MS::TIME) &&
1602 2424 : dynamic_cast<ObjCompare<Int>*>(sortDef.second.get()))
1603 0 : uniqueTimeInChunk = true;
1604 :
1605 14096 : for(auto& sortDef : subchunkSortColumns_p.sortingDefinition())
1606 21144 : if(sortDef.first == MS::columnName(MS::TIME) &&
1607 14096 : dynamic_cast<ObjCompare<Int>*>(sortDef.second.get()))
1608 0 : uniqueTimeInSubchunk = true;
1609 :
1610 7048 : if(uniqueTimeInChunk)
1611 0 : timeScope_p = ChunkScope;
1612 7048 : else if(uniqueTimeInSubchunk)
1613 0 : timeScope_p = SubchunkScope;
1614 : else
1615 7048 : timeScope_p = RowScope;
1616 :
1617 : // Determine the scope of the frequency/channel selections
1618 : // The scope of the frequency selections is at most the same as the DDID
1619 7048 : freqSelScope_p = ddIdScope_p;
1620 : // If frequency/channel selection also depends on time (selection based on
1621 : // frequencies), then the scope can be further limited by row (timestamp) scope
1622 7048 : if (!(frequencySelections_p->getFrameOfReference() == FrequencySelection::ByChannel))
1623 : {
1624 1724 : if(!(freqSelScope_p == RowScope)) // Only if scope is broader than Row can be further restricted
1625 : {
1626 0 : if(freqSelScope_p == SubchunkScope && timeScope_p == RowScope)
1627 0 : freqSelScope_p = RowScope;
1628 0 : if(freqSelScope_p == ChunkScope)
1629 0 : freqSelScope_p = timeScope_p;
1630 : }
1631 : }
1632 :
1633 7048 : }
1634 :
1635 5434 : VisibilityIteratorImpl2::Cache::Cache()
1636 : :
1637 : azel0Time_p(-1),
1638 : azelTime_p(-1),
1639 : feedpaTime_p(-1),
1640 : hourang_p(0),
1641 : hourangTime_p(-1),
1642 : msHasFlagCategory_p(false),
1643 : msHasWeightSpectrum_p(false),
1644 : msHasSigmaSpectrum_p(false),
1645 : parang0_p(0),
1646 : parang0Time_p(-1),
1647 5434 : parangTime_p(-1)
1648 5434 : {}
1649 :
1650 : void
1651 23385 : VisibilityIteratorImpl2::Cache::flush()
1652 : {
1653 23385 : azel0Time_p = -1;
1654 23385 : azelTime_p = -1;
1655 23385 : feedpaTime_p = -1;
1656 23385 : hourangTime_p = -1;
1657 23385 : parangTime_p = -1;
1658 23385 : parang0Time_p = -1;
1659 23385 : }
1660 :
1661 : const Cube<RigidVector<Double, 2> > &
1662 0 : VisibilityIteratorImpl2::getBeamOffsets() const
1663 : {
1664 0 : return msIter_p->getBeamOffsets();
1665 : }
1666 :
1667 : std::pair <bool, casacore::MDirection>
1668 7123259 : VisibilityIteratorImpl2::getPointingAngle(int antenna, double time) const
1669 : {
1670 7123259 : if (!pointingDirectionCache_p) {
1671 64 : pointingSource_p.reset(
1672 64 : new PointingColumns(subtableColumns_p->pointing()));
1673 64 : pointingDirectionCache_p.reset(
1674 64 : new PointingDirectionCache(nAntennas(), *pointingSource_p.get()));
1675 : }
1676 :
1677 : return pointingDirectionCache_p->getPointingDirection(
1678 7123259 : antenna, time, phaseCenter());
1679 : }
1680 :
1681 :
1682 : Vector<Double>
1683 1139398 : VisibilityIteratorImpl2::getFrequencies(Double time, Int frameOfReference,
1684 : Int spectralWindowId, Int msId) const
1685 : {
1686 : const SpectralWindowChannels & spectralWindowChannels =
1687 1139398 : getSpectralWindowChannels(msId, spectralWindowId);
1688 :
1689 : // Get the channel numbers selected for this time (the spectral window and
1690 : // MS index are assumed to be the same as those currently pointed to by the
1691 : // MSIter).
1692 :
1693 : Vector<Int> channels =
1694 2278796 : getChannels(time, frameOfReference, spectralWindowId, msId);
1695 :
1696 1139398 : Vector<Double> frequencies(channels.nelements());
1697 : // MFrequency::Types observatoryFrequencyType =
1698 : // getObservatoryFrequencyType();
1699 :
1700 :
1701 : MFrequency::Types measurementFrequencyType =
1702 1139398 : static_cast<MFrequency::Types>(getMeasurementFrame(spectralWindowId));
1703 :
1704 1139398 : if (frameOfReference == measurementFrequencyType) {
1705 :
1706 : // Since the observed frequency is being asked for, no conversion is
1707 : // necessary. Just copy each selected channel's observed frequency over
1708 : // to the result.
1709 :
1710 123473936 : for (Int i = 0; i < (int) channels.nelements(); i++) {
1711 :
1712 122914122 : Int channelNumber = channels[i];
1713 :
1714 122914122 : frequencies[i] = spectralWindowChannels.getFrequency(channelNumber);
1715 : }
1716 :
1717 559814 : return frequencies;
1718 : }
1719 :
1720 : // Get the converter from the observed to the requested frame.
1721 :
1722 : MFrequency::Convert fromObserved =
1723 : makeFrequencyConverter(time, spectralWindowId, frameOfReference,
1724 1738752 : false, Unit(String("Hz")));
1725 :
1726 : // For each selected channel, get its observed frequency and convert it into
1727 : // the frequency in the requested frame. Put the result into the
1728 : // corresponding slot of "frequencies".
1729 :
1730 5406298 : for (Int i = 0; i < (int) channels.nelements(); i++) {
1731 :
1732 4826714 : Int channelNumber = channels[i];
1733 :
1734 4826714 : Double fO = spectralWindowChannels.getFrequency(channelNumber);
1735 : // Observed frequency
1736 :
1737 4826714 : Double fF = fromObserved(fO).getValue();
1738 : // Frame frequency
1739 :
1740 4826714 : frequencies[i] = fF;
1741 : }
1742 :
1743 579584 : return frequencies;
1744 : }
1745 :
1746 : Vector<Double>
1747 478 : VisibilityIteratorImpl2::getChanWidths(Double time, Int frameOfReference,
1748 : Int spectralWindowId, Int msId) const
1749 : {
1750 : // This gets the native frame channel widths (no frame conversions performed, for now)
1751 :
1752 : const SpectralWindowChannels & spectralWindowChannels =
1753 478 : getSpectralWindowChannels(msId, spectralWindowId);
1754 :
1755 : // Get the channel numbers selected for this time (the spectral window and
1756 : // MS index are assumed to be the same as those currently pointed to by the
1757 : // MSIter).
1758 :
1759 : Vector<Int> channels =
1760 956 : getChannels(time, frameOfReference, spectralWindowId, msId);
1761 :
1762 478 : Vector<Double> widths(channels.nelements());
1763 :
1764 37626 : for (Int i = 0; i < (int) channels.nelements(); i++) {
1765 :
1766 37148 : Int channelNumber = channels[i];
1767 :
1768 37148 : widths[i] = spectralWindowChannels.getWidth(channelNumber);
1769 :
1770 : }
1771 :
1772 956 : return widths;
1773 : }
1774 :
1775 : Vector<Int>
1776 1413082 : VisibilityIteratorImpl2::getChannels(Double time, Int /*frameOfReference*/,
1777 : Int spectralWindowId, Int msId) const
1778 : {
1779 : std::shared_ptr<ChannelSelector> channelSelector =
1780 2826164 : determineChannelSelection(time, spectralWindowId, -1, msId);
1781 :
1782 2826164 : return channelSelector->getChannels();
1783 : }
1784 :
1785 : Vector<Int>
1786 4129040 : VisibilityIteratorImpl2::getCorrelations() const
1787 : {
1788 4129040 : assert(!channelSelectors_p.empty());
1789 :
1790 4129040 : return channelSelectors_p[0]->getCorrelations();
1791 : }
1792 :
1793 : Vector<Stokes::StokesTypes>
1794 7222964 : VisibilityIteratorImpl2::getCorrelationTypesDefined() const
1795 : {
1796 7222964 : assert(!channelSelectors_p.empty());
1797 :
1798 14445928 : Vector<Int> typesAsInt;
1799 7222964 : Int polarizationId = channelSelectors_p[0]->polarizationId;
1800 7222964 : subtableColumns_p->polarization().corrType().get(
1801 : polarizationId, typesAsInt, true);
1802 7222964 : Vector<Stokes::StokesTypes> correlationTypesDefined(typesAsInt.size());
1803 :
1804 29168504 : for (uInt i = 0; i < typesAsInt.size(); i++) {
1805 65836620 : correlationTypesDefined(i) =
1806 21945540 : static_cast<Stokes::StokesTypes>(typesAsInt(i));
1807 : }
1808 :
1809 14445928 : return correlationTypesDefined;
1810 : }
1811 :
1812 : Vector<Stokes::StokesTypes>
1813 3611482 : VisibilityIteratorImpl2::getCorrelationTypesSelected() const
1814 : {
1815 3611482 : assert(!channelSelectors_p.empty());
1816 :
1817 7222964 : Vector<Int> correlationIndices = getCorrelations();
1818 : Vector<Stokes::StokesTypes> correlationTypesDefined =
1819 7222964 : getCorrelationTypesDefined();
1820 : Vector<Stokes::StokesTypes> correlationTypesSelected(
1821 3611482 : correlationIndices.size());
1822 :
1823 14549140 : for (uInt i = 0; i < correlationIndices.size(); i++) {
1824 32812974 : correlationTypesSelected(i) =
1825 10937658 : correlationTypesDefined(correlationIndices(i));
1826 : }
1827 :
1828 7222964 : return correlationTypesSelected;
1829 : }
1830 :
1831 : Double
1832 0 : VisibilityIteratorImpl2::getInterval() const
1833 : {
1834 0 : return timeInterval_p;
1835 : }
1836 :
1837 : int
1838 4073846 : VisibilityIteratorImpl2::getMeasurementFrame(int spectralWindowId) const
1839 : {
1840 4073846 : int frame = subtableColumns_p->spectralWindow().measFreqRef()(
1841 : spectralWindowId);
1842 :
1843 4073846 : return frame;
1844 : }
1845 :
1846 :
1847 : Bool
1848 3615194 : VisibilityIteratorImpl2::isNewArrayId() const
1849 : {
1850 3615194 : return msIter_p->newArray();
1851 : }
1852 :
1853 : Bool
1854 3615194 : VisibilityIteratorImpl2::isNewFieldId() const
1855 : {
1856 3615194 : return msIter_p->newField();
1857 : }
1858 :
1859 : Bool
1860 3885471 : VisibilityIteratorImpl2::isNewMs() const
1861 : {
1862 3885471 : return msIter_p->newMS();
1863 : }
1864 :
1865 : Bool
1866 3615194 : VisibilityIteratorImpl2::isNewSpectralWindow() const
1867 : {
1868 3615194 : return msIter_p->newSpectralWindow();
1869 : }
1870 :
1871 : Bool
1872 7143628 : VisibilityIteratorImpl2::allBeamOffsetsZero() const
1873 : {
1874 7143628 : return msIter_p->allBeamOffsetsZero();
1875 : }
1876 :
1877 : rownr_t
1878 6998 : VisibilityIteratorImpl2::nRows() const
1879 : {
1880 6998 : return rowBounds_p.subchunkNRows_p;
1881 : }
1882 :
1883 : rownr_t
1884 0 : VisibilityIteratorImpl2::nShapes() const
1885 : {
1886 0 : return 1;
1887 : }
1888 :
1889 : const casacore::Vector<casacore::rownr_t>&
1890 521270 : VisibilityIteratorImpl2::nRowsPerShape () const
1891 : {
1892 521270 : return nRowsPerShape_p;
1893 : }
1894 :
1895 : const casacore::Vector<casacore::Int>&
1896 293298 : VisibilityIteratorImpl2::nChannelsPerShape () const
1897 : {
1898 293298 : return nChannPerShape_p;
1899 : }
1900 :
1901 : const casacore::Vector<casacore::Int>&
1902 517558 : VisibilityIteratorImpl2::nCorrelationsPerShape () const
1903 : {
1904 517558 : return nCorrsPerShape_p;
1905 : }
1906 :
1907 25735 : rownr_t VisibilityIteratorImpl2::nRowsInChunk() const
1908 : {
1909 25735 : return msIter_p->table().nrow();
1910 : }
1911 :
1912 251 : Int VisibilityIteratorImpl2::nTimes() const {
1913 251 : static const auto timeName = MeasurementSet::columnName(MSMainEnums::TIME);
1914 753 : auto times = ScalarColumn<Double>(msIter_p->table(), timeName).getColumn();
1915 251 : std::set<Double> uniqueTimes(times.cbegin(), times.cend());
1916 502 : return uniqueTimes.size();
1917 : }
1918 :
1919 : Bool
1920 3815760 : VisibilityIteratorImpl2::more() const
1921 : {
1922 3815760 : return more_p;
1923 : }
1924 :
1925 : Bool
1926 307302 : VisibilityIteratorImpl2::moreChunks() const
1927 : {
1928 307302 : return msIter_p->more();
1929 : }
1930 :
1931 : void
1932 1062 : VisibilityIteratorImpl2::result(casacore::Record& res) const
1933 : {
1934 1062 : if (moreChunks()) {
1935 : throw AipsError("TransformingVi2::result(Record&) can only be called at the end of "
1936 : "the iteration. It has been called while there are still "
1937 0 : "moreChunks(). Please check and/or revisit this condition.");
1938 : }
1939 : // For now nothing to add to result record from here
1940 1062 : }
1941 :
1942 : const MSColumns *
1943 0 : VisibilityIteratorImpl2::msColumnsKluge() const
1944 : {
1945 0 : return & msIter_p->msColumns();
1946 : }
1947 :
1948 : Int
1949 259096630 : VisibilityIteratorImpl2::msId() const
1950 : {
1951 259096630 : return msIter_p->msId();
1952 : }
1953 :
1954 : const MeasurementSet &
1955 4852124 : VisibilityIteratorImpl2::ms() const
1956 : {
1957 4852124 : return msIter_p->ms();
1958 : }
1959 :
1960 : String
1961 521306 : VisibilityIteratorImpl2::msName() const
1962 : {
1963 : // Name of current MS
1964 521306 : return ms().tableName();
1965 : }
1966 :
1967 : void
1968 2492372 : VisibilityIteratorImpl2::fieldIds(Vector<Int> & fieldIds) const
1969 : {
1970 2492372 : getColumnRows(columns_p.field_p, fieldIds);
1971 2492372 : }
1972 :
1973 : // Return the current ArrayId
1974 :
1975 : void
1976 976229 : VisibilityIteratorImpl2::arrayIds(Vector<Int> & arrayIds) const
1977 : {
1978 976229 : getColumnRows(columns_p.array_p, arrayIds);
1979 976229 : }
1980 :
1981 : // Return the current Field Name
1982 : String
1983 0 : VisibilityIteratorImpl2::fieldName() const
1984 : {
1985 0 : return msIter_p->fieldName();
1986 : }
1987 :
1988 : // Return the current Source Name
1989 : String
1990 0 : VisibilityIteratorImpl2::sourceName() const
1991 : {
1992 0 : return msIter_p->sourceName();
1993 : }
1994 :
1995 : const Vector<String> &
1996 0 : VisibilityIteratorImpl2::antennaMounts() const
1997 : {
1998 0 : return msIter_p->antennaMounts();
1999 : }
2000 :
2001 : void
2002 0 : VisibilityIteratorImpl2::setInterval(Double timeInterval)
2003 : {
2004 0 : pendingChanges_p->setInterval(timeInterval);
2005 0 : }
2006 :
2007 : void
2008 824 : VisibilityIteratorImpl2::setRowBlocking(rownr_t nRow)
2009 : {
2010 824 : pendingChanges_p->setNRowBlocking(nRow);
2011 824 : }
2012 :
2013 : rownr_t
2014 52303 : VisibilityIteratorImpl2::getRowBlocking() const
2015 : {
2016 52303 : return nRowBlocking_p;
2017 : }
2018 :
2019 : const MDirection &
2020 7131622 : VisibilityIteratorImpl2::phaseCenter() const
2021 : {
2022 7131622 : return msIter_p->phaseCenter();
2023 : }
2024 :
2025 : Int
2026 93626 : VisibilityIteratorImpl2::polFrame() const
2027 : {
2028 93626 : return msIter_p->polFrame();
2029 : }
2030 :
2031 : void
2032 5118806 : VisibilityIteratorImpl2::spectralWindows(Vector<Int> & spws) const
2033 : {
2034 : // Get's the list of spectral windows for each row in the VB window
2035 :
2036 10237612 : Vector<Int> ddis;
2037 5118806 : dataDescriptionIds(ddis);
2038 5118806 : spws.resize(ddis.size());
2039 :
2040 507472815 : for (uInt idx = 0; idx < ddis.size(); idx++) {
2041 1004708018 : spws(idx) = subtableColumns_p->dataDescription().spectralWindowId()(
2042 502354009 : ddis(idx));
2043 : }
2044 :
2045 10237612 : return;
2046 : }
2047 :
2048 : void
2049 1776704 : VisibilityIteratorImpl2::polarizationIds(Vector<Int> & polIds) const
2050 : {
2051 : // Get's the list of polarization Ids for each row in the VB window
2052 :
2053 3553408 : Vector<Int> ddis;
2054 1776704 : dataDescriptionIds(ddis);
2055 1776704 : polIds.resize(ddis.size());
2056 :
2057 254605369 : for (uInt idx = 0; idx < ddis.size(); idx++) {
2058 505657330 : polIds(idx) = subtableColumns_p->dataDescription().polarizationId()(
2059 252828665 : ddis(idx));
2060 : }
2061 :
2062 3553408 : return;
2063 : }
2064 :
2065 : // Return current Polarization Id
2066 : Int
2067 17626 : VisibilityIteratorImpl2::polarizationId() const
2068 : {
2069 17626 : return msIter_p->polarizationId();
2070 : }
2071 :
2072 : // Return current DataDescription Id
2073 : Int
2074 47949 : VisibilityIteratorImpl2::dataDescriptionId() const
2075 : {
2076 47949 : return msIter_p->dataDescriptionId();
2077 : }
2078 :
2079 : void
2080 7284775 : VisibilityIteratorImpl2::dataDescriptionIds(Vector<Int> & ddis) const
2081 : {
2082 7284775 : getColumnRows(columns_p.dataDescription_p, ddis);
2083 7284775 : }
2084 :
2085 : Bool
2086 0 : VisibilityIteratorImpl2::newFieldId() const
2087 : {
2088 0 : return (rowBounds_p.subchunkBegin_p == 0 && msIter_p->newField());
2089 : }
2090 :
2091 : Bool
2092 0 : VisibilityIteratorImpl2::newArrayId() const
2093 : {
2094 0 : return (rowBounds_p.subchunkBegin_p == 0 && msIter_p->newArray());
2095 : }
2096 :
2097 : Bool
2098 0 : VisibilityIteratorImpl2::newSpectralWindow() const
2099 : {
2100 0 : return (rowBounds_p.subchunkBegin_p == 0 && msIter_p->newSpectralWindow());
2101 : }
2102 :
2103 : Bool
2104 1217436 : VisibilityIteratorImpl2::existsColumn(VisBufferComponent2 id) const
2105 : {
2106 : Bool result;
2107 1217436 : switch (id) {
2108 :
2109 747066 : case VisBufferComponent2::VisibilityCorrected:
2110 : case VisBufferComponent2::VisibilityCubeCorrected:
2111 :
2112 747066 : result =
2113 747066 : !columns_p.corrVis_p.isNull() && columns_p.corrVis_p.isDefined(0);
2114 747066 : break;
2115 :
2116 9608 : case VisBufferComponent2::VisibilityModel:
2117 : case VisBufferComponent2::VisibilityCubeModel:
2118 :
2119 9608 : result =
2120 16348 : (!columns_p.modelVis_p.isNull() && columns_p.modelVis_p.isDefined(0)) ||
2121 6740 : modelDataGenerator_p != nullptr;
2122 9608 : break;
2123 :
2124 273 : case VisBufferComponent2::VisibilityObserved:
2125 : case VisBufferComponent2::VisibilityCubeObserved:
2126 :
2127 273 : result = (!columns_p.vis_p.isNull() && columns_p.vis_p.isDefined(0)) ||
2128 0 : (columns_p.floatVis_p.isNull() && columns_p.floatVis_p.isNull());
2129 :
2130 273 : break;
2131 :
2132 276542 : case VisBufferComponent2::VisibilityCubeFloat:
2133 :
2134 276542 : result =
2135 276542 : !columns_p.floatVis_p.isNull() && columns_p.floatVis_p.isDefined(0);
2136 :
2137 276542 : break;
2138 :
2139 39 : case VisBufferComponent2::WeightSpectrum:
2140 :
2141 39 : result =
2142 39 : !columns_p.weightSpectrum_p.isNull()
2143 39 : && columns_p.weightSpectrum_p.isDefined(0);
2144 39 : break;
2145 :
2146 0 : case VisBufferComponent2::SigmaSpectrum:
2147 :
2148 0 : result =
2149 0 : !columns_p.sigmaSpectrum_p.isNull()
2150 0 : && columns_p.sigmaSpectrum_p.isDefined(0);
2151 0 : break;
2152 :
2153 183908 : default:
2154 183908 : result = true; // required columns
2155 183908 : break;
2156 : }
2157 :
2158 1217436 : return result;
2159 : }
2160 :
2161 : const SubtableColumns &
2162 2799644 : VisibilityIteratorImpl2::subtableColumns() const
2163 : {
2164 2799644 : return *subtableColumns_p;
2165 : }
2166 :
2167 : void
2168 0 : VisibilityIteratorImpl2::allSpectralWindowsSelected(
2169 : Vector<Int> & selectedWindows,
2170 : Vector<Int> & nChannels) const
2171 : {
2172 :
2173 0 : Vector<Int> firstChannels; // ignored
2174 0 : Vector<Int> channelIncrement; // ignored
2175 :
2176 : // info generation should not use time as input
2177 0 : std::tie(selectedWindows, nChannels, firstChannels, channelIncrement) =
2178 0 : getChannelInformation();
2179 0 : }
2180 :
2181 : void
2182 3252 : VisibilityIteratorImpl2::useImagingWeight(const VisImagingWeight & imWgt)
2183 : {
2184 3252 : imwgt_p = imWgt;
2185 3252 : }
2186 :
2187 : void
2188 261717 : VisibilityIteratorImpl2::origin()
2189 : {
2190 261717 : ThrowIf(rowBounds_p.chunkNRows_p < 0,
2191 : "Call to origin without first initializing chunk");
2192 :
2193 261717 : throwIfPendingChanges();
2194 :
2195 261717 : rowBounds_p.subchunkBegin_p = 0; // begin at the beginning
2196 261717 : more_p = true;
2197 261717 : subchunk_p.resetSubChunk();
2198 :
2199 :
2200 261717 : if( ! (nRowBlocking_p > 0) )
2201 : {
2202 : // Create a MeasurementSet which points
2203 : // to the current iteration with msIter
2204 519534 : msSubchunk_p.reset(new casacore::MeasurementSet(msIter_p->table(),
2205 259767 : &(msIter_p->ms())));
2206 :
2207 : // Create a MSIter for the subchunk loop which iterates the
2208 : // the MS created before.
2209 519534 : msIterSubchunk_p.reset(new casacore::MSIter(*msSubchunk_p,
2210 259767 : subchunkSortColumns_p.sortingDefinition()));
2211 259767 : msIterSubchunk_p->origin();
2212 : }
2213 :
2214 261717 : configureNewSubchunk();
2215 261717 : }
2216 :
2217 : void
2218 0 : VisibilityIteratorImpl2::originChunks()
2219 : {
2220 0 : originChunks(false);
2221 0 : }
2222 :
2223 : void
2224 17847 : VisibilityIteratorImpl2::applyPendingChanges()
2225 : {
2226 17847 : if (!pendingChanges_p->empty()) {
2227 :
2228 : Bool exists;
2229 :
2230 : // Handle a pending frequency selection if it exists.
2231 :
2232 : FrequencySelections * newSelection;
2233 3605 : std::tie(exists, newSelection) =
2234 3605 : pendingChanges_p->popFrequencySelections();
2235 :
2236 3605 : if (exists) {
2237 :
2238 3524 : delete frequencySelections_p; // out with the old
2239 :
2240 3524 : frequencySelections_p = newSelection; // in with the new
2241 3524 : setMetadataScope();
2242 : }
2243 :
2244 : // Handle any pending interval change
2245 :
2246 : Double newInterval;
2247 3605 : std::tie(exists, newInterval) = pendingChanges_p->popInterval();
2248 :
2249 3605 : if (exists) {
2250 :
2251 0 : msIter_p->setInterval(newInterval);
2252 0 : timeInterval_p = newInterval;
2253 : }
2254 :
2255 : // Handle any row-blocking change
2256 :
2257 : Int newBlocking;
2258 3605 : std::tie(exists, newBlocking) = pendingChanges_p->popNRowBlocking();
2259 :
2260 3605 : if (exists) {
2261 :
2262 824 : nRowBlocking_p = newBlocking;
2263 :
2264 : }
2265 :
2266 : // force rewind since window selections may have changed
2267 3605 : msIterAtOrigin_p = false;
2268 : }
2269 17847 : }
2270 :
2271 : void
2272 17847 : VisibilityIteratorImpl2::originChunks(Bool forceRewind)
2273 : {
2274 17847 : subchunk_p.resetToOrigin();
2275 :
2276 17847 : applyPendingChanges();
2277 :
2278 17847 : if (!msIterAtOrigin_p || forceRewind) {
2279 :
2280 7606 : msIter_p->origin();
2281 7606 : msIterAtOrigin_p = true;
2282 :
2283 7606 : positionMsIterToASelectedSpectralWindow();
2284 :
2285 7606 : msIndex_p = msId();
2286 : }
2287 :
2288 17847 : configureNewChunk();
2289 17847 : }
2290 :
2291 : void
2292 259931 : VisibilityIteratorImpl2::positionMsIterToASelectedSpectralWindow()
2293 : {
2294 512912 : while (msIter_p->more() &&
2295 505962 : !frequencySelections_p->isSpectralWindowSelected(
2296 252981 : msIter_p->msId(), msIter_p->spectralWindowId())) {
2297 :
2298 0 : (* msIter_p)++;
2299 : }
2300 259931 : }
2301 :
2302 : void
2303 3077590 : VisibilityIteratorImpl2::next()
2304 : {
2305 3077590 : ThrowIf(!more_p, "Attempt to advance subchunk past end of chunk");
2306 :
2307 3077590 : throwIfPendingChanges(); // throw if unapplied changes exist
2308 :
2309 3077590 : measurementFrame_p = VisBuffer2::FrameNotSpecified; // flush cached value
2310 :
2311 3077590 : if(nRowBlocking_p > 0)
2312 : {
2313 : // Attempt to advance to the next subchunk
2314 3316 : rowBounds_p.subchunkBegin_p = rowBounds_p.subchunkEnd_p + 1;
2315 3316 : more_p = rowBounds_p.subchunkBegin_p < rowBounds_p.chunkNRows_p;
2316 : }
2317 : else
2318 : {
2319 : // Increment the subchunk MSIter
2320 3074274 : (*msIterSubchunk_p)++;
2321 3074274 : more_p = msIterSubchunk_p->more();
2322 : }
2323 :
2324 3077590 : if (more_p) {
2325 :
2326 2832207 : subchunk_p.incrementSubChunk();
2327 :
2328 2832207 : configureNewSubchunk();
2329 : }
2330 : else
2331 : {
2332 : // Leave the columns referencing a valid table. This ensures that some
2333 : // TVIs can still get some valid metadata when they are not at the end
2334 : // of iteration (even if the underlying VI2 is already at the end).
2335 245383 : attachColumns(msIter_p->table());
2336 : }
2337 3077590 : }
2338 :
2339 : Subchunk
2340 685438 : VisibilityIteratorImpl2::getSubchunkId() const
2341 : {
2342 685438 : return subchunk_p;
2343 : }
2344 :
2345 : const SortColumns &
2346 0 : VisibilityIteratorImpl2::getSortColumns() const
2347 : {
2348 0 : return sortColumns_p;
2349 : }
2350 :
2351 : void
2352 3591632 : VisibilityIteratorImpl2::throwIfPendingChanges()
2353 : {
2354 : // Throw an exception if there are any pending changes to the operation of
2355 : // the visibility iterator. The user must call originChunks to cause the
2356 : // changes to take effect; it is an error to try to advance the VI if there
2357 : // are unapplied changes pending.
2358 :
2359 3591632 : ThrowIf(!pendingChanges_p->empty(),
2360 : "Call to originChunks required after applying frequencySelection");
2361 :
2362 3591632 : }
2363 :
2364 : Bool
2365 0 : VisibilityIteratorImpl2::isInASelectedSpectralWindow() const
2366 : {
2367 0 : return frequencySelections_p->isSpectralWindowSelected(
2368 0 : msIter_p->msId(), msIter_p->spectralWindowId());
2369 : }
2370 :
2371 : Bool
2372 1559980 : VisibilityIteratorImpl2::isWritable() const
2373 : {
2374 1559980 : return writable_p;
2375 : }
2376 :
2377 : void
2378 252325 : VisibilityIteratorImpl2::nextChunk()
2379 : {
2380 252325 : ThrowIf(!msIter_p->more(),
2381 : "Attempt to advance past end of data using nextChunk");
2382 :
2383 252325 : throwIfPendingChanges(); // error if unapplied changes exist
2384 :
2385 : // Advance the MS Iterator until either there's no more data or it points to
2386 : // a selected spectral window.
2387 :
2388 252325 : (* msIter_p)++;
2389 :
2390 252325 : positionMsIterToASelectedSpectralWindow();
2391 :
2392 252325 : msIterAtOrigin_p = false;
2393 :
2394 : // If the MS Iterator was successfully advanced then
2395 : // set up for a new chunk
2396 :
2397 252325 : if (msIter_p->more()) {
2398 :
2399 245375 : subchunk_p.incrementChunk();
2400 :
2401 245375 : configureNewChunk();
2402 :
2403 245375 : vb_p->invalidate(); // flush the cache ??????????
2404 : }
2405 :
2406 252325 : more_p = msIter_p->more();
2407 252325 : }
2408 :
2409 : void
2410 3093924 : VisibilityIteratorImpl2::configureNewSubchunk()
2411 : {
2412 :
2413 : // Only for rowBlocking: work out how many rows to return for the moment
2414 : // we return all rows with
2415 : // the same value for time unless row blocking is set, in which case we
2416 : // return more rows at once.
2417 :
2418 3093924 : rowBounds_p.subchunkEqChanSelRows_p.clear();
2419 3093924 : if (nRowBlocking_p > 0) {
2420 3322 : rowBounds_p.subchunkEnd_p =
2421 3322 : rowBounds_p.subchunkBegin_p + nRowBlocking_p;
2422 :
2423 3322 : if (rowBounds_p.subchunkEnd_p >= rowBounds_p.chunkNRows_p) {
2424 1950 : rowBounds_p.subchunkEnd_p = rowBounds_p.chunkNRows_p - 1;
2425 : }
2426 : // This is needed because the call to spectralWindows() needs to
2427 : // have rowBounds_p.subchunkRows_p properly initialized
2428 : rowBounds_p.subchunkRows_p =
2429 3322 : RefRows(rowBounds_p.subchunkBegin_p, rowBounds_p.subchunkEnd_p);
2430 :
2431 :
2432 : // Scan the subchunk to see if the same channels are selected in each
2433 : // row. End the subchunk when a row using different channels is
2434 : // encountered.
2435 : Double previousRowTime =
2436 3322 : rowBounds_p.times_p(rowBounds_p.subchunkBegin_p);
2437 3322 : channelSelectors_p.clear();
2438 3322 : channelSelectorsNrows_p.clear();
2439 :
2440 3322 : channelSelectors_p.push_back(determineChannelSelection(previousRowTime,
2441 3322 : -1, polarizationId(), msId()));
2442 :
2443 1857354 : for (Int i = rowBounds_p.subchunkBegin_p + 1;
2444 1857354 : i <= rowBounds_p.subchunkEnd_p;
2445 : i++) {
2446 :
2447 1854032 : Double rowTime = rowBounds_p.times_p(i);
2448 :
2449 1854032 : if (rowTime == previousRowTime) {
2450 5914 : continue; // Same time means same rows.
2451 : }
2452 :
2453 : // Compute the channel selector for this row so it can be compared
2454 : // with the previous row's channel selector.
2455 :
2456 : std::shared_ptr<ChannelSelector> newSelector =
2457 : determineChannelSelection(rowTime, msIter_p->spectralWindowId(),
2458 3696236 : msIter_p->polarizationId(), msId());
2459 :
2460 1848118 : if (newSelector.get() != channelSelectors_p[0].get()) {
2461 :
2462 : // This row uses different channels than the previous row and so
2463 : // it cannot be included in this subchunk. Make the previous
2464 : // row the end of the subchunk.
2465 :
2466 0 : rowBounds_p.subchunkEnd_p = i - 1;
2467 : }
2468 : }
2469 : // Set the number of rows that use this channelSelector
2470 3322 : channelSelectorsNrows_p.push_back(rowBounds_p.subchunkEnd_p - rowBounds_p.subchunkBegin_p + 1);
2471 :
2472 3322 : rowBounds_p.subchunkNRows_p =
2473 3322 : rowBounds_p.subchunkEnd_p - rowBounds_p.subchunkBegin_p + 1;
2474 : // Reset this in case rowBounds_p.subchunkEnd_p has changed
2475 : rowBounds_p.subchunkRows_p =
2476 3322 : RefRows(rowBounds_p.subchunkBegin_p, rowBounds_p.subchunkEnd_p);
2477 3322 : rowBounds_p.subchunkEqChanSelRows_p.push_back(rowBounds_p.subchunkRows_p);
2478 : }
2479 : else {
2480 : // All the information is in the subchunk MSIter
2481 3090602 : rowBounds_p.subchunkNRows_p = msIterSubchunk_p->table().nrow();
2482 :
2483 3090602 : attachColumns(attachTable());
2484 :
2485 : // Fetch all of the times in this chunk and get the min/max
2486 : // of those times
2487 :
2488 3090602 : rowBounds_p.times_p.resize(rowBounds_p.subchunkNRows_p);
2489 3090602 : columns_p.time_p.getColumn(rowBounds_p.times_p);
2490 :
2491 : // The subchunk rows refer to the subchunk iterator
2492 : // and therefore are consecutive.
2493 3090602 : rowBounds_p.subchunkBegin_p = 0;
2494 3090602 : rowBounds_p.subchunkEnd_p = msIterSubchunk_p->table().nrow() - 1;
2495 3090602 : rowBounds_p.subchunkNRows_p =
2496 3090602 : rowBounds_p.subchunkEnd_p - rowBounds_p.subchunkBegin_p + 1;
2497 : rowBounds_p.subchunkRows_p =
2498 3090602 : RefRows(rowBounds_p.subchunkBegin_p, rowBounds_p.subchunkEnd_p);
2499 :
2500 :
2501 : // Under some circumstances, there is only one channel selector per chunk:
2502 : // 1. The selection doesn't depend on time (it is based only on channel number)
2503 : // and DDId (and consequently SPW, polID) is the same for the whole subchunk.
2504 : // 2. The selection might depend on time but DDid *and* time are the same for
2505 : // the whole subchunk.
2506 3090602 : if(freqSelScope_p == SubchunkScope)
2507 : {
2508 0 : channelSelectors_p.clear();
2509 0 : channelSelectorsNrows_p.clear();
2510 0 : double timeStamp = -1;
2511 0 : if(frequencySelections_p->getFrameOfReference() != FrequencySelection::ByChannel)
2512 0 : timeStamp = columns_p.time_p.asdouble(0);
2513 0 : channelSelectors_p.push_back(
2514 0 : determineChannelSelection(timeStamp,
2515 0 : msIterSubchunk_p->spectralWindowId(),
2516 0 : msIterSubchunk_p->polarizationId(), msId()));
2517 0 : channelSelectorsNrows_p.push_back(rowBounds_p.subchunkEnd_p - rowBounds_p.subchunkBegin_p + 1);
2518 0 : rowBounds_p.subchunkEqChanSelRows_p.push_back(rowBounds_p.subchunkRows_p);
2519 : }
2520 : // In all other cases the channel selector needs to be computed
2521 : // for each row. Each channel selector will then apply to a set
2522 : // of n consecutive rows (as defined in channelSelectorsNrows_p).
2523 3090602 : else if(freqSelScope_p == RowScope)
2524 : {
2525 1776704 : channelSelectors_p.clear();
2526 1776704 : channelSelectorsNrows_p.clear();
2527 3553408 : Vector<Int> spws, polIds;
2528 1776704 : spectralWindows(spws);
2529 1776704 : polarizationIds(polIds);
2530 1776704 : double timeStamp = -1;
2531 254605369 : for(Int irow = 0 ; irow < rowBounds_p.subchunkNRows_p; ++irow)
2532 : {
2533 252828665 : if(frequencySelections_p->getFrameOfReference() != FrequencySelection::ByChannel)
2534 222017779 : timeStamp = columns_p.time_p.asdouble(0);
2535 : auto newChannelSelector = determineChannelSelection(
2536 : timeStamp,
2537 505657330 : spws[irow], polIds[irow], msId());
2538 252828665 : if(irow == 0 || newChannelSelector != channelSelectors_p.back())
2539 : {
2540 2635920 : channelSelectors_p.push_back(newChannelSelector);
2541 2635920 : channelSelectorsNrows_p.push_back(1);
2542 : }
2543 : else
2544 250192745 : channelSelectorsNrows_p.back()++;
2545 : }
2546 1776704 : size_t beginRefRowIdx = rowBounds_p.subchunkBegin_p;
2547 4412624 : for (auto nrows : channelSelectorsNrows_p)
2548 : {
2549 2635920 : rowBounds_p.subchunkEqChanSelRows_p.push_back(RefRows(beginRefRowIdx, beginRefRowIdx + nrows - 1));
2550 2635920 : beginRefRowIdx += nrows;
2551 : }
2552 : }
2553 : // The remaining case is that scope of frequency selections is chunk.
2554 : // In this case the channelSelector is constant for a chunk
2555 : // and has already been computed in configureNewChunk.
2556 : // The number of rows still needs to be updated
2557 : // to account for the the number of rows in this subchunk
2558 : else
2559 : {
2560 1313898 : channelSelectorsNrows_p.clear();
2561 1313898 : channelSelectorsNrows_p.push_back(rowBounds_p.subchunkEnd_p - rowBounds_p.subchunkBegin_p + 1);
2562 1313898 : rowBounds_p.subchunkEqChanSelRows_p.push_back(rowBounds_p.subchunkRows_p);
2563 : }
2564 : }
2565 :
2566 : // Set flags for current subchunk
2567 :
2568 6187848 : Vector<Int> correlations = channelSelectors_p[0]->getCorrelations();
2569 3093924 : nCorrelations_p = correlations.nelements();
2570 :
2571 : Vector<Stokes::StokesTypes> correlationsDefined =
2572 6187848 : getCorrelationTypesDefined();
2573 : Vector<Stokes::StokesTypes> correlationsSelected =
2574 6187848 : getCorrelationTypesSelected();
2575 :
2576 3093924 : String msName = ms().tableName();
2577 :
2578 3093924 : auto nShapes = channelSelectors_p.size();
2579 3093924 : nRowsPerShape_p = Vector<uInt64>(channelSelectorsNrows_p);
2580 3093924 : nChannPerShape_p.resize(nShapes);
2581 3093924 : nCorrsPerShape_p.resize(nShapes);
2582 :
2583 3093924 : rownr_t ishape = 0;
2584 7047064 : for (auto channelSelector : channelSelectors_p)
2585 : {
2586 3953140 : nChannPerShape_p[ishape] = channelSelector->getNFrequencies();
2587 3953140 : nCorrsPerShape_p[ishape] = channelSelector->getCorrelations().nelements();
2588 3953140 : ++ishape;
2589 : }
2590 :
2591 3093924 : vb_p->configureNewSubchunk(
2592 3093924 : msId(), msName, isNewMs(), isNewArrayId(), isNewFieldId(),
2593 3093924 : isNewSpectralWindow(), subchunk_p,
2594 3093924 : nRowsPerShape_p,
2595 3093924 : nChannPerShape_p,
2596 3093924 : nCorrsPerShape_p,
2597 : correlations, correlationsDefined, correlationsSelected,
2598 3093924 : weightScaling_p);
2599 3093924 : }
2600 :
2601 : std::shared_ptr<ChannelSelector>
2602 256127026 : VisibilityIteratorImpl2::determineChannelSelection(
2603 : Double time,
2604 : Int spectralWindowId,
2605 : Int polarizationId,
2606 : Int msId) const
2607 : {
2608 : // --> The channels selected will need to be identical for all members of the
2609 : // subchunk; otherwise the underlying fetch method won't work since it takes
2610 : // a single Vector<Vector <Slice> > to specify the channels.
2611 :
2612 256127026 : assert(frequencySelections_p != 0);
2613 :
2614 256127026 : if (spectralWindowId == -1) {
2615 3322 : Vector<Int> spws;
2616 3322 : this->spectralWindows(spws);
2617 3322 : spectralWindowId = spws[0];
2618 : }
2619 :
2620 256127026 : if (msId < 0) {
2621 0 : msId = this->msId();
2622 : }
2623 :
2624 256127026 : const FrequencySelection & selection = frequencySelections_p->get(msId);
2625 256127026 : Int frameOfReference = selection.getFrameOfReference();
2626 :
2627 : // See if the appropriate channel selector is in the cache.
2628 :
2629 : std::shared_ptr<ChannelSelector> cachedSelector =
2630 256127026 : channelSelectorCache_p->find(time, msId, frameOfReference,
2631 512254052 : spectralWindowId);
2632 :
2633 256127026 : if (cachedSelector != nullptr) {
2634 255412458 : return cachedSelector;
2635 : }
2636 :
2637 : // Create the appropriate frequency selection for the current
2638 : // MS and Spectral Window
2639 :
2640 714568 : selection.filterByWindow(spectralWindowId);
2641 :
2642 : // Find(or create) the appropriate channel selection.
2643 :
2644 714568 : std::shared_ptr<ChannelSelector> newSelector;
2645 :
2646 714568 : if (polarizationId < 0) {
2647 831 : polarizationId = getPolarizationId(spectralWindowId, msId);
2648 : }
2649 :
2650 714568 : if (selection.getFrameOfReference() == FrequencySelection::ByChannel) {
2651 24914 : newSelector = makeChannelSelectorC(selection, time, msId,
2652 12457 : spectralWindowId, polarizationId);
2653 : }
2654 : else{
2655 1404222 : newSelector = makeChannelSelectorF(selection, time, msId,
2656 702111 : spectralWindowId, polarizationId);
2657 : }
2658 :
2659 : // Cache it for possible future use. The cache may hold multiple equivalent
2660 : // selectors, each having a different timestamp. Since selectors are small
2661 : // and there are not expected to be many equivalent selectors in the cache
2662 : // at a time, this shouldn't be a problem (the special case of selection by
2663 : // channel number is already handled).
2664 :
2665 714568 : channelSelectorCache_p->add(newSelector, frameOfReference);
2666 :
2667 714568 : return newSelector;
2668 : }
2669 :
2670 : Int
2671 831 : VisibilityIteratorImpl2::getPolarizationId(Int spectralWindowId, Int msId) const
2672 : {
2673 831 : ThrowIf(msId != this->msId(),
2674 : String::format("Requested MS number is %d but current is %d", msId,
2675 : this->msId()));
2676 :
2677 831 : auto & ddCols = subtableColumns_p->dataDescription();
2678 831 : auto nSpw = subtableColumns_p->spectralWindow().nrow();
2679 :
2680 : // This will break if the same spectral window is referenced by two
2681 : // different data_descrption IDs. Ideally, this whole thing should be
2682 : // reworked to used DDIDs with spectral window ID only used internally.
2683 831 : Int polID = -1;
2684 27184 : for (uInt idd = 0; idd < ddCols.spectralWindowId().nrow(); idd++) {
2685 26353 : if(ddCols.spectralWindowId()(idd) == spectralWindowId)
2686 831 : polID = ddCols.polarizationId()(idd);
2687 : }
2688 :
2689 : // If the SPW is not found in the DD it will return -1, rather than failing.
2690 : // This can happen for the so-called phantom SPWs. See CAS-11734
2691 831 : if(spectralWindowId < (Int)nSpw)
2692 831 : return polID;
2693 :
2694 : // spectralWindowId is not present in subtables
2695 0 : ThrowIf(true, String::format(
2696 : "Could not find entry for spectral window id"
2697 : "%d in spectral_window in MS #%d", spectralWindowId, msId));
2698 :
2699 0 : return -1; // Can't get here so make the compiler happy
2700 : }
2701 :
2702 :
2703 : std::shared_ptr<vi::ChannelSelector>
2704 12457 : VisibilityIteratorImpl2::makeChannelSelectorC(
2705 : const FrequencySelection & selectionIn,
2706 : Double time,
2707 : Int msId,
2708 : Int spectralWindowId,
2709 : Int polarizationId) const
2710 : {
2711 : const FrequencySelectionUsingChannels & selection =
2712 12457 : dynamic_cast<const FrequencySelectionUsingChannels &>(selectionIn);
2713 :
2714 12457 : if (selection.refinementNeeded())
2715 : {
2716 : auto spwcFetcher =
2717 : [this, msId]
2718 0 : (int spwId, double lowerFrequency, double upperFrequency)
2719 : {
2720 : const SpectralWindowChannels & spwChannels =
2721 0 : getSpectralWindowChannels(msId, spwId);
2722 : return spwChannels.getIntersection(
2723 0 : lowerFrequency, upperFrequency);
2724 0 : };
2725 0 : selection.applyRefinement(spwcFetcher);
2726 : }
2727 :
2728 24914 : vector<Slice> frequencySlices;
2729 :
2730 : // Iterate over all of the frequency selections for the specified spectral
2731 : // window and collect them into a vector of Slice objects.
2732 14775 : for (FrequencySelectionUsingChannels::const_iterator i = selection.begin();
2733 14775 : i != selection.end(); i++)
2734 : {
2735 :
2736 2318 : frequencySlices.push_back(i->getSlice());
2737 : }
2738 :
2739 12457 : if (frequencySlices.empty())
2740 : {
2741 : // And empty selection implies all channels
2742 : Int nChannels =
2743 10208 : subtableColumns_p->spectralWindow().numChan()(spectralWindowId);
2744 10208 : frequencySlices.push_back(Slice(0, nChannels, 1));
2745 : }
2746 :
2747 24914 : ChannelSlicer slices(2);
2748 :
2749 : // Install the polarization selections
2750 12457 : if(polarizationId != -1)
2751 : {
2752 : Vector<Slice> correlationSlices =
2753 12457 : selection.getCorrelationSlices(polarizationId);
2754 :
2755 12457 : if (correlationSlices.empty())
2756 : {
2757 : Int nCorrelations =
2758 12053 : subtableColumns_p->polarization().numCorr().get(polarizationId);
2759 12053 : correlationSlices = Vector<Slice>(1, Slice(0, nCorrelations));
2760 : }
2761 :
2762 12457 : slices.setSubslicer(0, ChannelSubslicer(correlationSlices));
2763 :
2764 : }
2765 :
2766 : // Create and install the frequency selections
2767 24914 : ChannelSubslicer frequencyAxis(frequencySlices.size());
2768 24983 : for (Int i = 0; i <(int) frequencySlices.size(); i++)
2769 : {
2770 12526 : frequencyAxis.setSlice(i, frequencySlices[i]);
2771 : }
2772 :
2773 12457 : slices.setSubslicer(1, frequencyAxis);
2774 :
2775 : // Package up the result and return it.
2776 : std::shared_ptr<ChannelSelector> result(new ChannelSelector
2777 12457 : (time, msId, spectralWindowId, polarizationId, slices));
2778 :
2779 24914 : return result;
2780 : }
2781 :
2782 : std::shared_ptr<ChannelSelector>
2783 702111 : VisibilityIteratorImpl2::makeChannelSelectorF(
2784 : const FrequencySelection & selectionIn,
2785 : Double time, Int msId, Int spectralWindowId,
2786 : Int polarizationId) const
2787 : {
2788 : // Make a ChannelSelector from a frame-based frequency selection.
2789 :
2790 : const FrequencySelectionUsingFrame & selection =
2791 702111 : dynamic_cast<const FrequencySelectionUsingFrame &>(selectionIn);
2792 :
2793 1404222 : vector<Slice> frequencySlices;
2794 :
2795 702111 : selection.filterByWindow(spectralWindowId);
2796 :
2797 : // Set up frequency converter so that we can convert to the measured
2798 : // frequency
2799 :
2800 : MFrequency::Convert convertToObservedFrame =
2801 : makeFrequencyConverter(
2802 : time, spectralWindowId, selection.getFrameOfReference(),
2803 1404222 : true, Unit("Hz"));
2804 :
2805 : // Convert each frequency selection into a Slice(interval) of channels.
2806 :
2807 : const SpectralWindowChannels & spectralWindowChannels =
2808 702111 : getSpectralWindowChannels(msId, spectralWindowId);
2809 :
2810 1414932 : for (FrequencySelectionUsingFrame::const_iterator i = selection.begin();
2811 1414932 : i != selection.end();
2812 712821 : i++) {
2813 :
2814 712821 : Double f = i->getBeginFrequency();
2815 712821 : Double lowerFrequency = convertToObservedFrame(f).getValue();
2816 :
2817 712821 : f = i->getEndFrequency();
2818 712821 : Double upperFrequency = convertToObservedFrame(f).getValue();
2819 :
2820 : Slice s =
2821 : findChannelsInRange(lowerFrequency, upperFrequency,
2822 712821 : spectralWindowChannels);
2823 :
2824 712821 : if (s.length() > 0) {
2825 712821 : frequencySlices.push_back(s);
2826 : }
2827 : }
2828 :
2829 : // Convert the STL vector of slices into the expected Casa Vector<Vector
2830 : // <Slice>> form. Element one of the Vector is empty indicating that the
2831 : // entire correlations axis is desired. The second element of the outer
2832 : // array specifies different channel intervals along the channel axis.
2833 :
2834 1404222 : ChannelSlicer slices(2);
2835 :
2836 : // Install the polarization selections
2837 :
2838 : Vector<Slice> correlationSlices =
2839 1404222 : selection.getCorrelationSlices(polarizationId);
2840 702111 : if (correlationSlices.empty()) {
2841 :
2842 : Int nCorrelations;
2843 702111 : subtableColumns_p->polarization().numCorr().get(
2844 : polarizationId, nCorrelations);
2845 :
2846 702111 : correlationSlices = Vector<Slice>(1, Slice(0, nCorrelations));
2847 : }
2848 :
2849 702111 : slices.setSubslicer(0, ChannelSubslicer(correlationSlices));
2850 :
2851 1404222 : ChannelSubslicer frequencyAxis(frequencySlices.size());
2852 :
2853 1414932 : for (Int i = 0; i <(int) frequencySlices.size(); i++) {
2854 712821 : frequencyAxis.setSlice(i, frequencySlices[i]);
2855 : }
2856 :
2857 702111 : slices.setSubslicer(1, frequencyAxis);
2858 :
2859 : // Package up result and return it.
2860 :
2861 : std::shared_ptr<ChannelSelector> result(new ChannelSelector
2862 702111 : (time, msId, spectralWindowId, polarizationId, slices));
2863 :
2864 1404222 : return result;
2865 : }
2866 :
2867 : MFrequency::Convert
2868 1281695 : VisibilityIteratorImpl2::makeFrequencyConverter(
2869 : Double time,
2870 : int spectralWindowId,
2871 : Int otherFrameOfReference,
2872 : Bool toObservedFrame,
2873 : Unit unit) const
2874 : {
2875 :
2876 : // Time is from the "time" field of the data row and is potentially in
2877 : // MFrequency::Types observatoryFrequencyType =
2878 1281695 : getObservatoryFrequencyType();
2879 :
2880 : MFrequency::Types measurementFrequencyType =
2881 1281695 : static_cast<MFrequency::Types>(getMeasurementFrame(spectralWindowId));
2882 :
2883 : // Set up frequency converter so that we can convert to the
2884 : // measured frequency
2885 : // =========================================================
2886 :
2887 : // Take the provided time in units native to the MS and attach the frame of
2888 : // reference so that the time users won't be confused(most MSs have time in
2889 : // UTC, but there are some that use different time standards).
2890 :
2891 3845085 : MEpoch epoch(MVEpoch(Quantity(time, "s")), timeFrameOfReference_p);
2892 :
2893 1281695 : const auto &telescopePosition = msIter_p->telescopePosition();
2894 1281695 : const auto ¤tDirection = msIter_p->phaseCenter();
2895 :
2896 2563390 : MeasFrame measFrame(epoch, telescopePosition, currentDirection);
2897 :
2898 2563390 : MFrequency::Ref observedFrame(measurementFrequencyType, measFrame);
2899 :
2900 1281695 : MFrequency::Types selectionFrame =
2901 : static_cast<MFrequency::Types>(otherFrameOfReference);
2902 :
2903 1281695 : MFrequency::Convert result;
2904 :
2905 1281695 : if (toObservedFrame) {
2906 :
2907 702111 : result = MFrequency::Convert(unit, selectionFrame, observedFrame);
2908 : }
2909 : else{
2910 :
2911 579584 : result = MFrequency::Convert(unit, observedFrame, selectionFrame);
2912 : }
2913 :
2914 2563390 : return result;
2915 : }
2916 :
2917 : Slice
2918 712821 : VisibilityIteratorImpl2::findChannelsInRange(
2919 : Double lowerFrequency, Double upperFrequency,
2920 : const vi::SpectralWindowChannels & spectralWindowChannels) const
2921 : {
2922 712821 : ThrowIf(spectralWindowChannels.empty(),
2923 : String::format(
2924 : "No spectral window channel info for ms=%d", msId()));
2925 :
2926 712821 : return spectralWindowChannels.getIntersection(lowerFrequency, upperFrequency);
2927 : }
2928 :
2929 : Int
2930 3361 : VisibilityIteratorImpl2::getNMs() const
2931 : {
2932 3361 : return measurementSets_p.nelements();
2933 : }
2934 :
2935 :
2936 : MFrequency::Types
2937 1281695 : VisibilityIteratorImpl2::getObservatoryFrequencyType() const
2938 : {
2939 1281695 : Int const spwId = msIter_p->spectralWindowId();
2940 1281695 : Int const measFreqRef = getMeasurementFrame(spwId);
2941 1281695 : return (measFreqRef >= 0) ? MFrequency::castType((uInt)measFreqRef)
2942 1281695 : : MFrequency::DEFAULT;
2943 : }
2944 :
2945 : MPosition
2946 0 : VisibilityIteratorImpl2::getObservatoryPosition() const
2947 : {
2948 0 : return msIter_p->telescopePosition();
2949 : }
2950 :
2951 : const SpectralWindowChannels &
2952 1841987 : VisibilityIteratorImpl2::getSpectralWindowChannels(
2953 : Int msId,
2954 : Int spectralWindowId) const
2955 : {
2956 : const SpectralWindowChannels * cached =
2957 1841987 : spectralWindowChannelsCache_p->find(msId, spectralWindowId);
2958 :
2959 1841987 : if (cached != 0) {
2960 1838546 : return * cached;
2961 : }
2962 :
2963 : // Get the columns for the spectral window subtable and then select out the
2964 : // frequency and width columns. Use those columns to extract the frequency
2965 : // and width lists for the specified spectral window.
2966 :
2967 : const MSSpWindowColumns& spectralWindow =
2968 3441 : subtableColumns_p->spectralWindow();
2969 :
2970 3441 : const ArrayColumn<Double>& frequenciesColumn = spectralWindow.chanFreq();
2971 6882 : Vector<Double> frequencies;
2972 3441 : frequenciesColumn.get(spectralWindowId, frequencies, true);
2973 :
2974 3441 : const ArrayColumn<Double>& widthsColumn = spectralWindow.chanWidth();
2975 3441 : Vector<Double> widths;
2976 3441 : widthsColumn.get(spectralWindowId, widths, true);
2977 :
2978 3441 : Assert(!frequencies.empty());
2979 3441 : Assert(frequencies.size() == widths.size());
2980 :
2981 : // Use the frequencies and widths to fill out a vector of
2982 : // SpectralWindowChannel objects.(No: This array needs to be in order of
2983 : // increasing frequency.) N.B.: If frequencies are in random order(i.e.,
2984 : // rather than reverse order) then all sorts of things will break elsewhere.
2985 : // Width also needs to be positive.
2986 :
2987 : SpectralWindowChannels *result =
2988 3441 : new SpectralWindowChannels(frequencies, widths);
2989 3441 : bool inOrder = true;
2990 :
2991 501521 : for (Int i = 0; i <(int) frequencies.nelements(); i++) {
2992 498080 : (*result)[i] = SpectralWindowChannel(i, frequencies[i], abs(widths[i]));
2993 498080 : inOrder = inOrder &&(i == 0 || frequencies[i] > frequencies[i - 1]);
2994 : }
2995 :
2996 : // Sanity check for frequencies that aren't monotonically
2997 : // increasing/decreasing.
2998 :
2999 498080 : for (Int i = 1; i <(int) frequencies.nelements(); i++) {
3000 494639 : ThrowIf(
3001 : abs((*result)[i].getChannel() -(*result)[i-1].getChannel()) != 1,
3002 : String::format(
3003 : "Spectral window %d in MS #%d has random ordered frequencies",
3004 : spectralWindowId, msId));
3005 : }
3006 :
3007 3441 : spectralWindowChannelsCache_p->add(result, msId, spectralWindowId);
3008 :
3009 3441 : return * result;
3010 : }
3011 :
3012 : void
3013 263222 : VisibilityIteratorImpl2::configureNewChunk()
3014 : {
3015 263222 : rowBounds_p.chunkNRows_p = msIter_p->table().nrow();
3016 263222 : rowBounds_p.subchunkBegin_p = -1; // undefined value
3017 263222 : rowBounds_p.subchunkEnd_p = -1; // will increment to 1st row
3018 :
3019 263222 : cache_p.chunkRowIds_p.resize(0); // flush cached row number map.
3020 :
3021 : // If this is a new MeasurementSet then set up the antenna locations, etc.
3022 :
3023 263222 : if (msIter_p->newMS()) {
3024 :
3025 : // Flush some cache flag values
3026 :
3027 17951 : cache_p.flush();
3028 :
3029 17951 : msd_p.setAntennas(msIter_p->msColumns().antenna());
3030 :
3031 : // Grab the time frame of reference so that it can be converted to UTC
3032 : // for use in other frame of reference conversions.
3033 :
3034 17951 : timeFrameOfReference_p = msIter_p->msColumns().timeMeas()(0).getRef();
3035 :
3036 : }
3037 :
3038 263222 : if (isNewMs()) { // New ms so flush pointing caches(if they exist).
3039 17951 : pointingDirectionCache_p.reset();
3040 17951 : pointingSource_p.reset();
3041 : }
3042 :
3043 263222 : if (msIter_p->newField() || msIterAtOrigin_p) {
3044 33156 : msd_p.setFieldCenter(msIter_p->phaseCenter());
3045 : }
3046 :
3047 263222 : if(nRowBlocking_p >0)
3048 : {
3049 1905 : attachColumns(msIter_p->table());
3050 :
3051 : // Fetch all of the times in this chunk
3052 :
3053 1905 : rowBounds_p.times_p.resize(rowBounds_p.chunkNRows_p);
3054 1905 : columns_p.time_p.getColumn(rowBounds_p.times_p);
3055 : }
3056 : else
3057 : {
3058 : // Columns are attached to the msIter chunk iteration.
3059 : // This is needed for the call of setTileCache() below, which
3060 : // performs some tests on the attached columns
3061 : // Later, in configureNewSubchunk the columns are reset to
3062 : // the subchunk msIterSubchunk_p columns.
3063 261317 : attachColumns(msIter_p->table());
3064 : }
3065 :
3066 : // Reset channel selectors vector in each chunk
3067 263222 : channelSelectors_p.clear();
3068 :
3069 : // If frequency selections are constant for the whole chunk
3070 263222 : if (freqSelScope_p == ChunkScope)
3071 : {
3072 : // If selection does not depend on time, assign a time value of -1.
3073 : // Since there is only one DDId we can ask for a single SPWId, PolID
3074 : // in msIter.
3075 33839 : double timeStamp = -1;
3076 33839 : if(frequencySelections_p->getFrameOfReference() == FrequencySelection::ByChannel)
3077 33839 : timeStamp = msIter_p->msColumns().time().asdouble(0);
3078 : // Note that channelSelectorNRow is not initialized since it should
3079 : // refer to nunber of subchunk rows. It is properly set in configureNewSubchunk
3080 33839 : channelSelectors_p.push_back(
3081 67678 : determineChannelSelection(timeStamp,
3082 : msIter_p->spectralWindowId(),
3083 33839 : msIter_p->polarizationId(), msId()));
3084 : }
3085 :
3086 263222 : setTileCache();
3087 263222 : }
3088 :
3089 : const MSDerivedValues &
3090 0 : VisibilityIteratorImpl2::getMsd() const
3091 : {
3092 0 : return msd_p;
3093 : }
3094 :
3095 : void
3096 263222 : VisibilityIteratorImpl2::setTileCache()
3097 : {
3098 263222 : std::lock_guard<std::mutex> lck(*tileCacheModMtx_p);
3099 :
3100 263222 : if ((*tileCacheIsSet_p)[msId()]
3101 268714 : || !msIter_p->newMS()
3102 268714 : || autoTileCacheSizing_p) {
3103 257730 : return; // Only useful when at start of an MS
3104 : }
3105 :
3106 : // This function sets the tile cache because of a feature in sliced data
3107 : // access that grows memory dramatically in some cases if (useSlicer_p) {
3108 :
3109 : // if (!(msIter_p->newDataDescriptionId() || msIter_p->newMS()) ) {
3110 : // return;
3111 : // }
3112 :
3113 5492 : const MeasurementSet & theMs = msIter_p->ms();
3114 5492 : if (theMs.tableType() == Table::Memory) {
3115 0 : return;
3116 : }
3117 :
3118 : // if (autoTileCacheSizing_p) {
3119 : // return; // take the default behavior
3120 : // }
3121 :
3122 : vector<MSMainEnums::PredefinedColumns> columnIds =
3123 : { MS::CORRECTED_DATA, MS::DATA, MS::FLAG, MS::MODEL_DATA, MS::SIGMA,
3124 10984 : MS::SIGMA_SPECTRUM, MS::UVW, MS::WEIGHT, MS::WEIGHT_SPECTRUM };
3125 :
3126 5492 : vector<String> msNames;
3127 :
3128 5492 : if (theMs.tableInfo().subType() == "CONCATENATED") {
3129 :
3130 126 : Block<String> names = theMs.getPartNames(false);
3131 63 : msNames.assign(names.begin(), names.end());
3132 :
3133 338 : for (String msName : msNames) {
3134 275 : MeasurementSet anMs(msName);
3135 :
3136 275 : setMsCacheSizes(anMs, columnIds);
3137 : }
3138 :
3139 : } else {
3140 5429 : setMsCacheSizes(theMs, columnIds);
3141 : }
3142 5492 : (*tileCacheIsSet_p)[msId()] = true;
3143 : }
3144 :
3145 5704 : void VisibilityIteratorImpl2::setMsCacheSizes(
3146 : const MeasurementSet & ms,
3147 : vector<MSMainEnums::PredefinedColumns> columnIds)
3148 : {
3149 5704 : const ColumnDescSet & cds = ms.tableDesc().columnDescSet();
3150 :
3151 57040 : for (MSMainEnums::PredefinedColumns columnId : columnIds) {
3152 :
3153 51336 : String column = MS::columnName(columnId);
3154 :
3155 : try {
3156 :
3157 51336 : if (!cds.isDefined(column) || !usesTiledDataManager(column, ms)) {
3158 20521 : continue; // skip if column not in MS or not using tiles
3159 : }
3160 :
3161 30815 : if (columnId == MS::WEIGHT_SPECTRUM
3162 29611 : || columnId == MS::SIGMA_SPECTRUM) {
3163 :
3164 : // These two columns are frequently present in an MS but
3165 : // uninitialized.
3166 :
3167 1550 : TableColumn tableColumn(ms, column);
3168 1550 : if (!tableColumn.hasContent()) {
3169 0 : continue; // Skip
3170 : }
3171 : }
3172 :
3173 30815 : setMsColumnCacheSizes(ms, column);
3174 :
3175 0 : } catch(AipsError & e) {
3176 0 : continue; // It failed so leave the caching as is
3177 : }
3178 : }
3179 5704 : }
3180 :
3181 : void
3182 30815 : VisibilityIteratorImpl2::setMsColumnCacheSizes(const MeasurementSet & partMs,
3183 : const string & column)
3184 : {
3185 : // For the column in the provided MS, loop over the hypercubes and
3186 : // set the cache size appropriately.
3187 :
3188 61630 : ROTiledStManAccessor accessor(partMs, column, true);
3189 30815 : uInt nHypercubes = accessor.nhypercubes();
3190 :
3191 90698 : for (uInt cube = 0; cube != nHypercubes; cube++) {
3192 :
3193 : // Get hypercube shape(includes row axis) and tile shape(does not
3194 : // include the row axis).
3195 :
3196 59883 : const IPosition tileShape(accessor.getTileShape(cube));
3197 59883 : const IPosition hypercubeShape(accessor.getHypercubeShape(cube));
3198 :
3199 59883 : uInt nAxes = hypercubeShape.size(); // how many axes
3200 59883 : if (nAxes < 1) {
3201 :
3202 : // Empty hypercube so skip it. Can't rely on loop below to handle
3203 : // this case because nAxes is unsigned so nAxes-1 is going to wrap
3204 : // around and become yuuge!
3205 :
3206 27006 : continue;
3207 : }
3208 :
3209 : // Compute the appropriate cache size which will hold at least a single
3210 : // row's worth of tiles in the cache. Use the factor of 2 as the
3211 : // initial value since some large Alma data sets cause the row to span
3212 : // tiles along the baseline axis due to the autocorrelations not
3213 : // occurring near the corresponding cross correlation baselines.
3214 :
3215 : // Doubling to handle case where baselines span tiles in the row
3216 : // direction
3217 32877 : uInt cacheSize = 2;
3218 :
3219 : // Compute the number of tiles required to span the non-row axes of the
3220 : // hypercube.
3221 :
3222 83994 : for (uInt axis = 0; axis < nAxes - 1; ++axis) {
3223 51117 : cacheSize *=
3224 51117 : (uInt)ceil(hypercubeShape[axis] /(Float)(tileShape[axis]));
3225 : }
3226 :
3227 : // Apply the cache size(in tiles).
3228 :
3229 32877 : accessor.setHypercubeCacheSize(cube, cacheSize);
3230 : }
3231 30815 : }
3232 :
3233 :
3234 : Bool
3235 35987 : VisibilityIteratorImpl2::usesTiledDataManager(
3236 : const String & columnName,
3237 : const MeasurementSet & theMs) const
3238 : {
3239 35987 : Bool noData = false;
3240 :
3241 : // Have to do something special about weight_spectrum as it tend to exist
3242 : // but has no valid data.
3243 :
3244 71974 : noData = noData ||
3245 35987 : (columnName == MS::columnName(MS::WEIGHT_SPECTRUM)
3246 2636 : && !weightSpectrumExists());
3247 :
3248 70542 : noData = noData ||
3249 34555 : (columnName == MS::columnName(MS::SIGMA_SPECTRUM)
3250 447 : && !sigmaSpectrumExists());
3251 :
3252 : // Check to see if the column exist and have valid data
3253 :
3254 70466 : noData = noData ||
3255 39461 : (columnName == MS::columnName(MS::DATA) &&
3256 9964 : (columns_p.vis_p.isNull() || !columns_p.vis_p.isDefined(0)));
3257 :
3258 70466 : noData = noData ||
3259 36972 : (columnName == MS::columnName(MS::MODEL_DATA) &&
3260 4986 : (columns_p.modelVis_p.isNull() || !columns_p.modelVis_p.isDefined(0)));
3261 :
3262 70466 : noData = noData ||
3263 37092 : (columnName == MS::columnName(MS::CORRECTED_DATA) &&
3264 5226 : (columns_p.corrVis_p.isNull() || !columns_p.corrVis_p.isDefined(0)));
3265 :
3266 70466 : noData = noData ||
3267 40183 : (columnName == MS::columnName(MS::FLAG) &&
3268 11408 : (columns_p.flag_p.isNull() || !columns_p.flag_p.isDefined(0)));
3269 :
3270 70466 : noData = noData ||
3271 40183 : (columnName == MS::columnName(MS::WEIGHT) &&
3272 11408 : (columns_p.weight_p.isNull() || !columns_p.weight_p.isDefined(0)));
3273 :
3274 70466 : noData = noData ||
3275 40183 : (columnName == MS::columnName(MS::SIGMA) &&
3276 11408 : (columns_p.sigma_p.isNull() || !columns_p.sigma_p.isDefined(0)));
3277 :
3278 70466 : noData = noData ||
3279 40183 : (columnName == MS::columnName(MS::UVW) &&
3280 11408 : (columns_p.uvw_p.isNull() || !columns_p.uvw_p.isDefined(0)));
3281 :
3282 35987 : Bool usesTiles = false;
3283 :
3284 35987 : if (!noData) {
3285 : String dataManType =
3286 68958 : RODataManAccessor(theMs, columnName, true).dataManagerType();
3287 :
3288 34479 : usesTiles = dataManType.contains("Tiled");
3289 : }
3290 :
3291 35987 : return usesTiles;
3292 : }
3293 :
3294 : void
3295 3599207 : VisibilityIteratorImpl2::attachColumns(const Table & t)
3296 : {
3297 3599207 : columns_p.attachColumns(t);
3298 :
3299 3599207 : floatDataFound_p = columns_p.isFloatDataPresent();
3300 3599207 : }
3301 :
3302 : MEpoch
3303 131749 : VisibilityIteratorImpl2::getEpoch() const
3304 : {
3305 131749 : MEpoch mEpoch = msIter_p->msColumns().timeMeas()(0);
3306 :
3307 131749 : return mEpoch;
3308 : }
3309 :
3310 : Vector<Float>
3311 0 : VisibilityIteratorImpl2::getReceptor0Angle()
3312 : {
3313 0 : Int nAntennas = this->nAntennas();
3314 :
3315 0 : Vector<Float> receptor0Angle(nAntennas);
3316 :
3317 0 : for (int i = 0; i < nAntennas; i++) {
3318 0 : receptor0Angle[i] = msIter_p->receptorAngle()(0, i);
3319 : }
3320 :
3321 0 : return receptor0Angle;
3322 : }
3323 :
3324 : void
3325 26116 : VisibilityIteratorImpl2::getRowIds(Vector<rownr_t> & rowIds) const
3326 : {
3327 26116 : if(nRowBlocking_p >0)
3328 : {
3329 : // Resize the rowIds vector and fill it with the row numbers contained in
3330 : // the current subchunk. These row numbers are relative to the reference
3331 : // table used by MSIter to define a chunk; thus rowId 0 is the first row in
3332 : // the chunk.
3333 219 : rowIds.resize(rowBounds_p.subchunkNRows_p);
3334 219 : rowIds = rowBounds_p.subchunkRows_p.convert();
3335 :
3336 219 : if (cache_p.chunkRowIds_p.nelements() == 0) {
3337 : // Create chunkRowIds_p as a "map" from chunk rows to MS rows. This
3338 : // needs to be created once per chunk since a new reference table is
3339 : // created each time the MSIter moves to the next chunk.
3340 219 : cache_p.chunkRowIds_p = msIter_p->table().rowNumbers(msIter_p->ms());
3341 : }
3342 :
3343 : // Using chunkRowIds_p as a map from chunk rows to MS rows replace the
3344 : // chunk-relative row numbers with the actual row number from the MS.
3345 2814 : for (uInt i = 0; i < rowIds.nelements(); i++) {
3346 2595 : rowIds(i) = cache_p.chunkRowIds_p(rowIds(i));
3347 : }
3348 : }
3349 : else
3350 : {
3351 : // Resize the rowIds vector and fill it with the row numbers contained in
3352 : // the current subchunk.
3353 25897 : rowIds.resize(rowBounds_p.subchunkNRows_p);
3354 :
3355 : // Initialize the cache it if not yet done
3356 : // (it is reset each time nextChunk() is called).
3357 : // The cache contains the mapping between chunk rows and MS rows.
3358 : // This needs to be created once per chunk since a new reference table is
3359 : // created each time the MSIter moves to the next chunk.
3360 25897 : if (cache_p.chunkRowIds_p.size() == 0)
3361 1166 : cache_p.chunkRowIds_p = msIter_p->table().rowNumbers(msIter_p->ms());
3362 :
3363 : // Now create the map from subchunk rows to chunk rows. This
3364 : // needs to be created for each subchunk since a new reference table
3365 : // in the msIterInner_p iterator is created each time the MSIter moves.
3366 : // Note that what we get are row Ids for msIter_p->table(), which is itself
3367 : // a reference table.
3368 77691 : auto subchunkRowIds = msIterSubchunk_p->table().rowNumbers(msIter_p->table(), true);
3369 :
3370 : // Now, for each row in the subchunk (i), get the row in the outer loop
3371 : // table (subchunkRowId(i)) and use cache_p.chunkRowIds_p to get the row
3372 : // in the original MS.
3373 5381974 : for (uInt i = 0; i < rowIds.size(); i++)
3374 5356077 : rowIds(i) = cache_p.chunkRowIds_p(subchunkRowIds(i));
3375 : }
3376 26116 : }
3377 :
3378 : void
3379 2198944 : VisibilityIteratorImpl2::antenna1(Vector<Int> & ant1) const
3380 : {
3381 2198944 : getColumnRows(columns_p.antenna1_p, ant1);
3382 2198944 : }
3383 :
3384 : void
3385 1865002 : VisibilityIteratorImpl2::antenna2(Vector<Int> & ant2) const
3386 : {
3387 1865002 : getColumnRows(columns_p.antenna2_p, ant2);
3388 1865002 : }
3389 :
3390 : void
3391 518936 : VisibilityIteratorImpl2::feed1(Vector<Int> & fd1) const
3392 : {
3393 518936 : getColumnRows(columns_p.feed1_p, fd1);
3394 518936 : }
3395 :
3396 : void
3397 353815 : VisibilityIteratorImpl2::feed2(Vector<Int> & fd2) const
3398 : {
3399 353815 : getColumnRows(columns_p.feed2_p, fd2);
3400 353815 : }
3401 :
3402 : void
3403 740803 : VisibilityIteratorImpl2::corrType(Vector<Int> & corrTypes) const
3404 : {
3405 740803 : Int polId = msIter_p->polarizationId();
3406 :
3407 740803 : subtableColumns_p->polarization().corrType().get(polId, corrTypes, true);
3408 740803 : }
3409 :
3410 : void
3411 2586206 : VisibilityIteratorImpl2::flag(Cube<Bool> & flags) const
3412 : {
3413 2586206 : getColumnRows(columns_p.flag_p, flags);
3414 2586206 : }
3415 :
3416 : void
3417 0 : VisibilityIteratorImpl2::flag(Vector<Cube<Bool>> & flags) const
3418 : {
3419 0 : getColumnRows(columns_p.flag_p, flags);
3420 0 : }
3421 :
3422 : void
3423 0 : VisibilityIteratorImpl2::flag(Matrix<Bool> & flags) const
3424 : {
3425 0 : Cube<Bool> flagCube;
3426 :
3427 0 : flag(flagCube);
3428 :
3429 0 : getColumnRows(columns_p.flag_p, flagCube);
3430 :
3431 0 : uInt nChannels = flagCube.shape()(1);
3432 :
3433 0 : flags.resize(nChannels, rowBounds_p.subchunkNRows_p);
3434 :
3435 0 : for (Int row = 0; row < rowBounds_p.subchunkNRows_p; row++) {
3436 :
3437 0 : for (uInt channel = 0; channel < nChannels; channel++) {
3438 :
3439 0 : Bool flagIt = flagCube(0, channel, row);
3440 :
3441 0 : for (Int correlation = 1;
3442 0 : correlation < nCorrelations_p && not flagIt;
3443 : correlation++) {
3444 :
3445 0 : flagIt = flagCube(correlation, channel, row);
3446 : }
3447 :
3448 0 : flags(channel, row) = flagIt;
3449 : }
3450 : }
3451 0 : }
3452 :
3453 : Bool
3454 0 : VisibilityIteratorImpl2::flagCategoryExists() const
3455 : {
3456 0 : if (msIter_p->newMS()) { // Cache to avoid testing unnecessarily.
3457 0 : cache_p.msHasFlagCategory_p = columns_p.flagCategory_p.hasContent();
3458 : }
3459 0 : return cache_p.msHasFlagCategory_p;
3460 : }
3461 :
3462 : void
3463 0 : VisibilityIteratorImpl2::flagCategory(Array<Bool> & /*flagCategories*/) const
3464 : {
3465 0 : ThrowIf(true, "The flag_category column is not supported.");
3466 : // if (columns_p.flagCategory_p.isNull() ||
3467 : // !columns_p.flagCategory_p.isDefined(0)) { // It often is.
3468 : //
3469 : // flagCategories.resize(); // Zap it.
3470 : // }
3471 : // else {
3472 : //
3473 : // // Since flag category is shaped [nC, nF, nCategories] it requires a
3474 : // // slightly different slicer and cannot use the usual getColumns
3475 : // method.
3476 : //
3477 : // const ChannelSlicer & channelSlicer =
3478 : // channelSelectors_p[0]->getSlicerForFlagCategories();
3479 : //
3480 : // columns_p.flagCategory_p.getSliceForRows(
3481 : // rowBounds_p.subchunkRows_p,
3482 : // channelSlicer.getSlicerInCoreRep(),
3483 : // flagCategories);
3484 : // }
3485 0 : }
3486 :
3487 : void
3488 2040771 : VisibilityIteratorImpl2::flagRow(Vector<Bool> & rowflags) const
3489 : {
3490 2040771 : getColumnRows(columns_p.flagRow_p, rowflags);
3491 2040771 : }
3492 :
3493 : void
3494 1346936 : VisibilityIteratorImpl2::observationId(Vector<Int> & obsIDs) const
3495 : {
3496 1346936 : getColumnRows(columns_p.observation_p, obsIDs);
3497 1346936 : }
3498 :
3499 : void
3500 339090 : VisibilityIteratorImpl2::processorId(Vector<Int> & procIDs) const
3501 : {
3502 339090 : getColumnRows(columns_p.processor_p, procIDs);
3503 339090 : }
3504 :
3505 : void
3506 1562516 : VisibilityIteratorImpl2::scan(Vector<Int> & scans) const
3507 : {
3508 1562516 : getColumnRows(columns_p.scan_p, scans);
3509 1562516 : }
3510 :
3511 : void
3512 805065 : VisibilityIteratorImpl2::stateId(Vector<Int> & stateIds) const
3513 : {
3514 805065 : getColumnRows(columns_p.state_p, stateIds);
3515 805065 : }
3516 :
3517 : void
3518 2131647 : VisibilityIteratorImpl2::time(Vector<Double> & t) const
3519 : {
3520 2131647 : getColumnRows(columns_p.time_p, t);
3521 2131647 : }
3522 :
3523 : void
3524 935433 : VisibilityIteratorImpl2::timeCentroid(Vector<Double> & t) const
3525 : {
3526 935433 : getColumnRows(columns_p.timeCentroid_p, t);
3527 935433 : }
3528 :
3529 : void
3530 1067799 : VisibilityIteratorImpl2::timeInterval(Vector<Double> & t) const
3531 : {
3532 1067799 : getColumnRows(columns_p.timeInterval_p, t);
3533 1067799 : }
3534 :
3535 : void
3536 907536 : VisibilityIteratorImpl2::exposure(Vector<Double> & expo) const
3537 : {
3538 907536 : getColumnRows(columns_p.exposure_p, expo);
3539 907536 : }
3540 :
3541 : void
3542 432466 : VisibilityIteratorImpl2::visibilityCorrected(Cube<Complex> & vis) const
3543 : {
3544 432466 : if(columns_p.corrVis_p.isNull())
3545 0 : throw AipsError("Requesting visibilityCorrected but column is null");
3546 432466 : getColumnRows(columns_p.corrVis_p, vis);
3547 432466 : }
3548 :
3549 : void
3550 0 : VisibilityIteratorImpl2::visibilityCorrected(Vector<Cube<Complex>> & vis) const
3551 : {
3552 0 : if(columns_p.corrVis_p.isNull())
3553 0 : throw AipsError("Requesting visibilityCorrected but column is null");
3554 0 : getColumnRows(columns_p.corrVis_p, vis);
3555 0 : }
3556 :
3557 : void
3558 449122 : VisibilityIteratorImpl2::visibilityModel(Cube<Complex> & vis) const
3559 : {
3560 : // See if the data can be filled from a virtual model column; if not then
3561 : // get it from the model column.
3562 :
3563 449122 : if (!fillFromVirtualModel(vis)) {
3564 186560 : getColumnRows(columns_p.modelVis_p, vis);
3565 : }
3566 449122 : }
3567 :
3568 : void
3569 0 : VisibilityIteratorImpl2::visibilityModel(Vector<Cube<Complex>> & vis) const
3570 : {
3571 0 : if (!fillFromVirtualModel(vis[0])) {
3572 0 : getColumnRows(columns_p.modelVis_p, vis);
3573 : }
3574 : else
3575 0 : throw AipsError("VisibilityIteratorImpl2::visibilityModel(Vector<Cube<Complex>> & vis) from model not yet implemented");
3576 0 : }
3577 :
3578 : void
3579 1205021 : VisibilityIteratorImpl2::visibilityObserved(Cube<Complex> & vis) const
3580 : {
3581 1205021 : if (floatDataFound_p) {
3582 :
3583 : // Since there is a floating data column, read that and convert it into
3584 : // the expected Complex form.
3585 :
3586 219578 : Cube<Float> dataFloat;
3587 :
3588 109789 : getColumnRows(columns_p.floatVis_p, dataFloat);
3589 :
3590 109789 : vis.resize(dataFloat.shape());
3591 :
3592 109789 : convertArray(vis, dataFloat);
3593 : }
3594 : else {
3595 1095232 : if(columns_p.vis_p.isNull())
3596 0 : throw AipsError("Requesting visibilityObserved but column is null");
3597 1095232 : getColumnRows(columns_p.vis_p, vis);
3598 : }
3599 1205021 : }
3600 :
3601 : void
3602 0 : VisibilityIteratorImpl2::visibilityObserved(Vector<Cube<Complex>> & vis) const
3603 : {
3604 0 : if (floatDataFound_p) {
3605 :
3606 : // Since there is a floating data column, read that and convert it into
3607 : // the expected Complex form.
3608 :
3609 0 : Vector<Cube<Float>> dataFloat;
3610 :
3611 0 : getColumnRows(columns_p.floatVis_p, dataFloat);
3612 :
3613 0 : vis.resize(dataFloat.size());
3614 :
3615 0 : size_t iVec = 0;
3616 0 : for (iVec= 0; iVec < vis.size(); iVec++)
3617 : {
3618 0 : vis[iVec].resize(dataFloat[iVec].shape());
3619 :
3620 0 : convertArray(vis[iVec], dataFloat[iVec]);
3621 : }
3622 : }
3623 : else {
3624 0 : if(columns_p.vis_p.isNull())
3625 0 : throw AipsError("Requesting visibilityObserved but column is null");
3626 0 : getColumnRows(columns_p.vis_p, vis);
3627 : }
3628 0 : }
3629 :
3630 : void
3631 114083 : VisibilityIteratorImpl2::floatData(Cube<Float> & fcube) const
3632 : {
3633 114083 : if (floatDataFound_p) {
3634 114083 : getColumnRows(columns_p.floatVis_p, fcube);
3635 : }
3636 : else{
3637 0 : fcube.resize();
3638 : }
3639 114083 : }
3640 :
3641 : void
3642 0 : VisibilityIteratorImpl2::floatData(Vector<Cube<Float>> & fcubes) const
3643 : {
3644 0 : if (floatDataFound_p) {
3645 0 : getColumnRows(columns_p.floatVis_p, fcubes);
3646 : }
3647 : else{
3648 0 : fcubes.resize(1);
3649 0 : fcubes[0].resize();
3650 : }
3651 0 : }
3652 :
3653 : void
3654 1403862 : VisibilityIteratorImpl2::uvw(Matrix<Double> & uvwmat) const
3655 : {
3656 1403862 : getColumnRowsMatrix(columns_p.uvw_p, uvwmat, false);
3657 1403862 : }
3658 :
3659 : // Fill in parallactic angle.
3660 : const Vector<Float> &
3661 7415896 : VisibilityIteratorImpl2::feed_pa(Double time) const
3662 : {
3663 : // Absolute UT
3664 :
3665 7415896 : Double ut = time;
3666 :
3667 7415896 : if (ut != cache_p.feedpaTime_p) {
3668 :
3669 : // Set up the Epoch using the absolute MJD in seconds get the Epoch
3670 : // reference from the column
3671 :
3672 194666 : MEpoch mEpoch = getEpoch();
3673 :
3674 194666 : const Matrix<Double> & angles = receptorAngles().xyPlane(0);
3675 97333 : Int nAnt = angles.shape()(1);
3676 :
3677 97333 : Vector<Float> receptor0Angle(nAnt, 0);
3678 :
3679 1262454 : for (int i = 0; i < nAnt; i++) {
3680 1165121 : receptor0Angle[i] = angles(0, i);
3681 : }
3682 :
3683 97333 : cache_p.feedpa_p.assign(
3684 194666 : feed_paCalculate(time, msd_p, nAnt, mEpoch, receptor0Angle));
3685 :
3686 97333 : cache_p.feedpaTime_p = ut;
3687 : }
3688 7415896 : return cache_p.feedpa_p;
3689 : }
3690 :
3691 : // Fill in parallactic angle.
3692 : const Float &
3693 0 : VisibilityIteratorImpl2::parang0(Double time) const
3694 : {
3695 0 : if (time != cache_p.parang0Time_p) {
3696 :
3697 0 : cache_p.parang0Time_p = time;
3698 :
3699 : // Set up the Epoch using the absolute MJD in seconds get the Epoch
3700 : // reference from the column
3701 0 : MEpoch mEpoch = getEpoch();
3702 0 : cache_p.parang0_p = parang0Calculate(time, msd_p, mEpoch);
3703 : }
3704 0 : return cache_p.parang0_p;
3705 : }
3706 :
3707 : // Fill in parallactic angle(NO FEED PA offset!).
3708 : const Vector<Float> &
3709 0 : VisibilityIteratorImpl2::parang(Double time) const
3710 : {
3711 0 : if (time != cache_p.parangTime_p) {
3712 :
3713 0 : cache_p.parangTime_p = time;
3714 :
3715 : // Set up the Epoch using the absolute MJD in seconds get the Epoch
3716 : // reference from the column
3717 :
3718 0 : MEpoch mEpoch = getEpoch();
3719 0 : Int nAnt = msIter_p->receptorAngle().shape()(1);
3720 :
3721 0 : cache_p.parang_p = parangCalculate(time, msd_p, nAnt, mEpoch);
3722 :
3723 : }
3724 0 : return cache_p.parang_p;
3725 : }
3726 :
3727 : // Fill in azimuth/elevation of the antennas. Cloned from feed_pa, we need to
3728 : // check that this is all correct!
3729 : const Vector<MDirection> &
3730 159282 : VisibilityIteratorImpl2::azel(Double ut) const
3731 : {
3732 159282 : if (ut != cache_p.azelTime_p) {
3733 :
3734 1131 : cache_p.azelTime_p = ut;
3735 :
3736 1131 : Int nAnt = msIter_p->receptorAngle().shape()(1);
3737 :
3738 2262 : MEpoch mEpoch = getEpoch();
3739 :
3740 1131 : azelCalculate(ut, msd_p, cache_p.azel_p, nAnt, mEpoch);
3741 :
3742 : }
3743 159282 : return cache_p.azel_p;
3744 : }
3745 :
3746 :
3747 : // Fill in azimuth/elevation of the antennas. Cloned from feed_pa, we need to
3748 : // check that this is all correct!
3749 : MDirection
3750 0 : VisibilityIteratorImpl2::azel0(Double time) const
3751 : {
3752 : // Absolute UT
3753 :
3754 0 : Double ut = time;
3755 :
3756 0 : if (ut != cache_p.azel0Time_p) {
3757 :
3758 0 : cache_p.azel0Time_p = ut;
3759 :
3760 0 : MEpoch mEpoch = getEpoch();
3761 :
3762 0 : azel0Calculate(time, msd_p, cache_p.azel0_p, mEpoch);
3763 :
3764 : }
3765 0 : return cache_p.azel0_p;
3766 : }
3767 :
3768 :
3769 : // Hour angle at specified time.
3770 : Double
3771 0 : VisibilityIteratorImpl2::hourang(Double time) const
3772 : {
3773 : // Absolute UT
3774 :
3775 0 : Double ut = time;
3776 :
3777 0 : if (ut != cache_p.hourangTime_p) {
3778 :
3779 0 : cache_p.hourangTime_p = ut;
3780 :
3781 : // Set up the Epoch using the absolute MJD in seconds get the Epoch
3782 : // reference from the column keyword
3783 :
3784 0 : MEpoch mEpoch = getEpoch();
3785 :
3786 0 : cache_p.hourang_p = hourangCalculate(time, msd_p, mEpoch);
3787 :
3788 : }
3789 0 : return cache_p.hourang_p;
3790 : }
3791 :
3792 :
3793 : void
3794 1242959 : VisibilityIteratorImpl2::sigma(Matrix<Float> & sigma) const
3795 : {
3796 1242959 : getColumnRowsMatrix(columns_p.sigma_p, sigma, true);
3797 1242959 : }
3798 :
3799 : void
3800 0 : VisibilityIteratorImpl2::sigma(Vector<Matrix<Float>> & sigma) const
3801 : {
3802 0 : getColumnRowsMatrix(columns_p.sigma_p, sigma);
3803 0 : }
3804 :
3805 : void
3806 1187483 : VisibilityIteratorImpl2::weight(Matrix<Float> & wt) const
3807 : {
3808 1187483 : getColumnRowsMatrix(columns_p.weight_p, wt, true);
3809 1187483 : }
3810 :
3811 : void
3812 0 : VisibilityIteratorImpl2::weight(Vector<Matrix<Float>> & wt) const
3813 : {
3814 0 : getColumnRowsMatrix(columns_p.weight_p, wt);
3815 0 : }
3816 :
3817 : Bool
3818 1675289 : VisibilityIteratorImpl2::weightSpectrumExists() const
3819 : {
3820 1675289 : if (msIter_p->newSpectralWindow()) {
3821 : // Cache to avoid testing unnecessarily.
3822 1364761 : cache_p.msHasWeightSpectrum_p = columns_p.weightSpectrum_p.hasContent();
3823 : }
3824 :
3825 1675289 : return cache_p.msHasWeightSpectrum_p;
3826 : }
3827 :
3828 : Bool
3829 627349 : VisibilityIteratorImpl2::sigmaSpectrumExists() const
3830 : {
3831 627349 : if (msIter_p->newMS()) { // Cache to avoid testing unnecessarily.
3832 :
3833 140859 : cache_p.msHasSigmaSpectrum_p = columns_p.sigmaSpectrum_p.hasContent();
3834 :
3835 : }
3836 :
3837 627349 : return cache_p.msHasSigmaSpectrum_p;
3838 : }
3839 :
3840 : void
3841 36444 : VisibilityIteratorImpl2::weightSpectrum(Cube<Float> & spectrum) const
3842 : {
3843 36444 : if (weightSpectrumExists())
3844 36444 : getColumnRows(columns_p.weightSpectrum_p, spectrum);
3845 : else
3846 0 : spectrum.resize(0, 0, 0);
3847 36444 : }
3848 :
3849 : void
3850 0 : VisibilityIteratorImpl2::weightSpectrum(Vector<Cube<Float>> & spectrum) const
3851 : {
3852 0 : if (weightSpectrumExists())
3853 0 : getColumnRows(columns_p.weightSpectrum_p, spectrum);
3854 : else
3855 : {
3856 0 : spectrum.resize(1);
3857 0 : spectrum[0].resize(0, 0, 0);
3858 : }
3859 0 : }
3860 :
3861 : void
3862 39652 : VisibilityIteratorImpl2::sigmaSpectrum(Cube<Float> & spectrum) const
3863 : {
3864 39652 : if (sigmaSpectrumExists())
3865 38777 : getColumnRows(columns_p.sigmaSpectrum_p, spectrum);
3866 : else
3867 875 : spectrum.resize(0, 0, 0);
3868 39652 : }
3869 :
3870 : void
3871 0 : VisibilityIteratorImpl2::sigmaSpectrum(Vector<Cube<Float>> & spectrum) const
3872 : {
3873 0 : if (sigmaSpectrumExists())
3874 0 : getColumnRows(columns_p.sigmaSpectrum_p, spectrum);
3875 : else
3876 : {
3877 0 : spectrum.resize(1);
3878 0 : spectrum[0].resize(0, 0, 0);
3879 : }
3880 0 : }
3881 :
3882 : void
3883 427 : VisibilityIteratorImpl2::setWeightScaling(
3884 : CountedPtr<WeightScaling> weightScaling)
3885 : {
3886 427 : weightScaling_p = weightScaling;
3887 427 : }
3888 :
3889 : Bool
3890 0 : VisibilityIteratorImpl2::hasWeightScaling() const
3891 : {
3892 0 : return !weightScaling_p.null();
3893 : }
3894 :
3895 : CountedPtr<WeightScaling>
3896 521270 : VisibilityIteratorImpl2::getWeightScaling() const
3897 : {
3898 521270 : return weightScaling_p;
3899 : }
3900 :
3901 :
3902 : const VisImagingWeight &
3903 471125 : VisibilityIteratorImpl2::getImagingWeightGenerator() const
3904 : {
3905 471125 : return imwgt_p;
3906 : }
3907 :
3908 : Block<MeasurementSet>
3909 0 : VisibilityIteratorImpl2::getMeasurementSets() const
3910 : {
3911 0 : return measurementSets_p;
3912 : }
3913 :
3914 : Int
3915 31297816 : VisibilityIteratorImpl2::getReportingFrameOfReference() const
3916 : {
3917 : Int frame;
3918 31297816 : if (reportingFrame_p == VisBuffer2::FrameNotSpecified) {
3919 :
3920 31297816 : if (frequencySelections_p != 0) {
3921 :
3922 31297816 : frame = frequencySelections_p->getFrameOfReference();
3923 :
3924 31297816 : if (frame == FrequencySelection::ByChannel) {
3925 :
3926 : // Since selection was done by channels, the frequencies are
3927 : // native.
3928 :
3929 371058 : Vector<Int> spws;
3930 371058 : spectralWindows(spws);
3931 371058 : measurementFrame_p = getMeasurementFrame(spws[0]);
3932 371058 : frame = measurementFrame_p;
3933 : }
3934 : }
3935 : else{
3936 0 : frame = VisBuffer2::FrameNotSpecified;
3937 : }
3938 : }
3939 : else{
3940 0 : frame = reportingFrame_p;
3941 : }
3942 :
3943 31297816 : return frame;
3944 : }
3945 :
3946 : void
3947 0 : VisibilityIteratorImpl2::setReportingFrameOfReference(Int frame)
3948 : {
3949 0 : ThrowIf(frame < 0 || frame >= MFrequency::N_Types,
3950 : String::format("Unknown frame: id=%d", frame));
3951 :
3952 0 : reportingFrame_p = frame;
3953 0 : }
3954 :
3955 : VisBuffer2 *
3956 3477543 : VisibilityIteratorImpl2::getVisBuffer() const
3957 : {
3958 3477543 : return vb_p;
3959 : }
3960 :
3961 : VisBuffer2 *
3962 0 : VisibilityIteratorImpl2::getVisBuffer(const VisibilityIterator2 * vi) const
3963 : {
3964 0 : ThrowIf(vb_p == nullptr, "VI Implementation has no VisBuffer.");
3965 0 : vb_p->associateWithVi2(vi);
3966 0 : return vb_p;
3967 : }
3968 :
3969 :
3970 : Int
3971 163965 : VisibilityIteratorImpl2::nAntennas() const
3972 : {
3973 163965 : return subtableColumns_p->antenna().nrow(); // for single(sub)array only..
3974 : }
3975 :
3976 : Int
3977 15912 : VisibilityIteratorImpl2::nSpectralWindows() const
3978 : {
3979 15912 : return subtableColumns_p->spectralWindow().nrow();
3980 : }
3981 :
3982 : Int
3983 0 : VisibilityIteratorImpl2::nDataDescriptionIds() const
3984 : {
3985 0 : return subtableColumns_p->dataDescription().nrow();
3986 : }
3987 :
3988 : Int
3989 164 : VisibilityIteratorImpl2::nPolarizationIds() const
3990 : {
3991 164 : return subtableColumns_p->polarization().nrow();
3992 : }
3993 :
3994 : rownr_t
3995 0 : VisibilityIteratorImpl2::nRowsViWillSweep() const
3996 : {
3997 0 : Int numcoh = 0;
3998 0 : for (uInt k = 0; k < uInt(msIter_p->numMS()) ; ++k) {
3999 0 : numcoh += msIter_p->ms(k).nrow();
4000 : }
4001 0 : return numcoh;
4002 :
4003 : }
4004 :
4005 : const Table
4006 3090602 : VisibilityIteratorImpl2::attachTable() const
4007 : {
4008 3090602 : return msIterSubchunk_p->table();
4009 : }
4010 :
4011 : void
4012 1039 : VisibilityIteratorImpl2::slurp() const
4013 : {
4014 : // Set the table data manager(ISM and SSM) cache size to the full column
4015 : // size, for the columns ANTENNA1, ANTENNA2, FEED1, FEED2, TIME, INTERVAL,
4016 : // FLAG_ROW, SCAN_NUMBER and UVW
4017 :
4018 2078 : Record dmInfo(msIter_p->ms().dataManagerInfo());
4019 :
4020 2078 : RecordDesc desc = dmInfo.description();
4021 :
4022 19629 : for (unsigned i = 0; i < dmInfo.nfields(); i++) {
4023 :
4024 18590 : if (desc.isSubRecord(i)) {
4025 :
4026 37180 : Record sub = dmInfo.subRecord(i);
4027 :
4028 18590 : if (sub.fieldNumber("NAME") >= 0 &&
4029 37180 : sub.fieldNumber("TYPE") >= 0 &&
4030 37180 : sub.fieldNumber("COLUMNS") >= 0 &&
4031 37180 : sub.type(sub.fieldNumber("NAME")) == TpString &&
4032 74360 : sub.type(sub.fieldNumber("TYPE")) == TpString &&
4033 37180 : sub.type(sub.fieldNumber("COLUMNS")) == TpArrayString) {
4034 :
4035 37180 : Array<String> columns;
4036 18590 : dmInfo.subRecord(i).get("COLUMNS", columns);
4037 :
4038 18590 : bool match = false;
4039 :
4040 42869 : for (unsigned j = 0; j < columns.nelements(); j++) {
4041 :
4042 24279 : String column = columns(IPosition(1, j));
4043 :
4044 24279 : match |=(column == MS::columnName(MS::ANTENNA1) ||
4045 23240 : column == MS::columnName(MS::ANTENNA2) ||
4046 22201 : column == MS::columnName(MS::FEED1) ||
4047 21162 : column == MS::columnName(MS::FEED2) ||
4048 20123 : column == MS::columnName(MS::TIME) ||
4049 19084 : column == MS::columnName(MS::INTERVAL) ||
4050 18045 : column == MS::columnName(MS::FLAG_ROW) ||
4051 63486 : column == MS::columnName(MS::SCAN_NUMBER) ||
4052 15967 : column == MS::columnName(MS::UVW));
4053 : }
4054 :
4055 18590 : if (match) {
4056 :
4057 13832 : String dm_name;
4058 6916 : dmInfo.subRecord(i).get("NAME", dm_name);
4059 :
4060 13832 : String dm_type;
4061 6916 : dmInfo.subRecord(i).get("TYPE", dm_type);
4062 :
4063 6916 : Bool can_exceed_nr_buckets = false;
4064 6916 : uInt num_buckets = msIter_p->ms().nrow();
4065 : // One bucket is at least one row, so this is enough
4066 :
4067 6916 : if (dm_type == "IncrementalStMan") {
4068 :
4069 7880 : ROIncrementalStManAccessor acc(msIter_p->ms(), dm_name);
4070 3940 : acc.setCacheSize(num_buckets, can_exceed_nr_buckets);
4071 :
4072 2976 : } else if (dm_type == "StandardStMan") {
4073 :
4074 3952 : ROStandardStManAccessor acc(msIter_p->ms(), dm_name);
4075 1976 : acc.setCacheSize(num_buckets, can_exceed_nr_buckets);
4076 : }
4077 : /* These are the only storage managers which use the
4078 : BucketCache (and therefore are slow for random access and
4079 : small cache sizes)
4080 : */
4081 : }
4082 : else {
4083 :
4084 11674 : String dm_name;
4085 11674 : dmInfo.subRecord(i).get("NAME", dm_name);
4086 :
4087 : }
4088 : } else {
4089 0 : cerr << "Data manager info has unexpected shape! "
4090 0 : << sub << endl;
4091 : }
4092 : }
4093 : }
4094 2078 : return;
4095 : }
4096 :
4097 : const Cube<Double> &
4098 97333 : VisibilityIteratorImpl2::receptorAngles() const
4099 : {
4100 97333 : return msIter_p->receptorAngles();
4101 : }
4102 :
4103 : IPosition
4104 1098034 : VisibilityIteratorImpl2::visibilityShape() const
4105 : {
4106 :
4107 : IPosition result(3,
4108 1098034 : nCorrelations_p,
4109 1098034 : channelSelectors_p[0]->getNFrequencies(),
4110 2196068 : rowBounds_p.subchunkNRows_p);
4111 :
4112 1098034 : return result;
4113 : }
4114 :
4115 : void
4116 3524 : VisibilityIteratorImpl2::setFrequencySelections(
4117 : FrequencySelections const& frequencySelections)
4118 : {
4119 3524 : pendingChanges_p->setFrequencySelections(frequencySelections.clone());
4120 :
4121 3524 : channelSelectorCache_p->flush();
4122 3524 : spectralWindowChannelsCache_p->flush();
4123 3524 : channelSelectors_p.clear();
4124 3524 : channelSelectorsNrows_p.clear();
4125 3524 : setMetadataScope();
4126 3524 : }
4127 :
4128 : void
4129 0 : VisibilityIteratorImpl2::jonesC(
4130 : Vector<SquareMatrix<complex<float>, 2> >& cjones) const
4131 : {
4132 0 : cjones.resize(msIter_p->CJones().nelements());
4133 0 : cjones = msIter_p->CJones();
4134 0 : }
4135 :
4136 : void
4137 589286 : VisibilityIteratorImpl2::writeFlag(const Cube<Bool> & flags)
4138 : {
4139 589286 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4140 :
4141 589286 : putColumnRows(columns_p.flag_p, flags);
4142 589286 : }
4143 :
4144 : void
4145 0 : VisibilityIteratorImpl2::writeFlagCategory(const Array<Bool>& flagCategory)
4146 : {
4147 0 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4148 :
4149 : // Flag categories are [nC, nF, nCategories] and therefore must use a
4150 : // different slicer which also prevents use of more usual putColumn method.
4151 :
4152 0 : RefRows & rows = rowBounds_p.subchunkRows_p;
4153 : const ChannelSlicer & channelSlicer =
4154 0 : channelSelectors_p[0]->getSlicerForFlagCategories();
4155 :
4156 0 : columns_p.flagCategory_p.putSliceFromRows(
4157 0 : rows, channelSlicer.getSlicerInCoreRep(), flagCategory);
4158 0 : }
4159 :
4160 : void
4161 383243 : VisibilityIteratorImpl2::writeFlagRow(const Vector<Bool> & rowflags)
4162 : {
4163 383243 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4164 :
4165 383243 : putColumnRows(columns_p.flagRow_p, rowflags);
4166 383243 : }
4167 :
4168 : void
4169 178303 : VisibilityIteratorImpl2::writeVisCorrected(const Cube<Complex> & vis)
4170 : {
4171 178303 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4172 :
4173 178303 : putColumnRows(columns_p.corrVis_p, vis);
4174 178303 : }
4175 :
4176 : void
4177 14320 : VisibilityIteratorImpl2::writeVisModel(const Cube<Complex> & vis)
4178 : {
4179 14320 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4180 :
4181 14320 : putColumnRows(columns_p.modelVis_p, vis);
4182 14320 : }
4183 :
4184 : void
4185 0 : VisibilityIteratorImpl2::writeVisObserved(const Cube<Complex> & vis)
4186 : {
4187 0 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4188 :
4189 0 : if (floatDataFound_p) {
4190 :
4191 : // This MS has float data; convert the cube to float
4192 : // and right it out
4193 :
4194 0 : Cube<Float> dataFloat = real(vis);
4195 0 : putColumnRows(columns_p.floatVis_p, dataFloat);
4196 :
4197 : }
4198 : else {
4199 0 : putColumnRows(columns_p.vis_p, vis);
4200 : }
4201 :
4202 0 : }
4203 :
4204 : void
4205 178753 : VisibilityIteratorImpl2::writeWeight(const Matrix<Float> & weight)
4206 : {
4207 178753 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4208 :
4209 178753 : putColumnRows(columns_p.weight_p, weight);
4210 178753 : }
4211 :
4212 : void
4213 3059 : VisibilityIteratorImpl2::writeWeightSpectrum(const Cube<Float> & weightSpectrum)
4214 : {
4215 3059 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4216 :
4217 3059 : if (!columns_p.weightSpectrum_p.isNull()) {
4218 3059 : putColumnRows(columns_p.weightSpectrum_p, weightSpectrum);
4219 : }
4220 3059 : }
4221 :
4222 : void
4223 2064 : VisibilityIteratorImpl2::initWeightSpectrum(const Cube<Float> & weightSpectrum)
4224 : {
4225 2064 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4226 :
4227 2064 : if (!columns_p.weightSpectrum_p.isNull()) {
4228 1900 : RefRows & rows = rowBounds_p.subchunkRows_p;
4229 1900 : columns_p.weightSpectrum_p.putColumnCells(rows, weightSpectrum);
4230 : }
4231 2064 : }
4232 :
4233 : void
4234 672 : VisibilityIteratorImpl2::writeSigmaSpectrum(const Cube<Float> & sigmaSpectrum)
4235 : {
4236 672 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4237 :
4238 672 : if (!columns_p.sigmaSpectrum_p.isNull()) {
4239 672 : putColumnRows(columns_p.sigmaSpectrum_p, sigmaSpectrum);
4240 : }
4241 672 : }
4242 :
4243 : void
4244 1556 : VisibilityIteratorImpl2::initSigmaSpectrum(const Cube<Float> & sigmaSpectrum)
4245 : {
4246 1556 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4247 :
4248 1556 : if (!columns_p.sigmaSpectrum_p.isNull() ) {
4249 1392 : RefRows & rows = rowBounds_p.subchunkRows_p;
4250 1392 : columns_p.sigmaSpectrum_p.putColumnCells(rows, sigmaSpectrum);
4251 : }
4252 1556 : }
4253 :
4254 : void
4255 5942 : VisibilityIteratorImpl2::writeSigma(const Matrix<Float> & sigma)
4256 : {
4257 5942 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4258 :
4259 5942 : putColumnRows(columns_p.sigma_p, sigma);
4260 5942 : }
4261 :
4262 : void
4263 17 : VisibilityIteratorImpl2::writeModel(
4264 : const RecordInterface& rec, Bool iscomponentlist, Bool incremental)
4265 : {
4266 :
4267 17 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4268 : /* Version 1 stuff
4269 : Vector<Int> fields = columns_p.field_p.getColumn();
4270 :
4271 : const Int option = Sort::HeapSort | Sort::NoDuplicates;
4272 : const Sort::Order order = Sort::Ascending;
4273 :
4274 : Int nFields = GenSort<Int>::sort(fields, order, option);
4275 :
4276 : // Make sure we have the right size
4277 :
4278 : fields.resize(nFields, true);
4279 : */
4280 :
4281 34 : Matrix<Int> combiIndex;
4282 17 : MSUtil::getIndexCombination(MSColumns(ms()), combiIndex);
4283 34 : Vector<Int> selectedWindows;
4284 34 : Vector<Int> nChannels;
4285 34 : Vector<Int> firstChannels;
4286 34 : Vector<Int> channelIncrement;
4287 :
4288 0 : std::tie(selectedWindows, nChannels, firstChannels, channelIncrement) =
4289 17 : getChannelInformation();
4290 34 : Matrix<Int> chansel(selectedWindows.nelements(),4);
4291 17 : chansel.column(0)=selectedWindows;
4292 17 : chansel.column(1)=firstChannels;
4293 17 : chansel.column(2)=nChannels;
4294 17 : chansel.column(3)=channelIncrement;
4295 34 : CountedPtr<VisModelDataI> visModelData = VisModelDataI::create();
4296 :
4297 : /*visModelData->putModelI(
4298 : ms(), rec, fields, selectedWindows, firstChannels, nChannels,
4299 : channelIncrement, iscomponentlist, incremental);*/
4300 : //Version 2 interface to keep state and scan number in track
4301 17 : visModelData->putModelI (ms(), rec, combiIndex, chansel, iscomponentlist, incremental);
4302 :
4303 17 : }
4304 :
4305 : VisibilityIteratorImpl2::ChannelInfo
4306 17 : VisibilityIteratorImpl2::getChannelInformationUsingFrequency() const
4307 : {
4308 : const FrequencySelectionUsingFrame *frequencySelection =
4309 0 : dynamic_cast<const FrequencySelectionUsingFrame*>(
4310 17 : &frequencySelections_p->get(msId()));
4311 17 : if (!frequencySelection)
4312 : throw(AipsError(
4313 0 : "Programmer Error channel info with wrong object called"));
4314 34 : set<Int> windows = frequencySelection->getSelectedWindows();
4315 :
4316 34 : Vector<Int> spectralWindow(windows.size());
4317 34 : Vector<Int> nChannels(windows.size(), -1);
4318 34 : Vector<Int> firstChannel(windows.size(), -1);
4319 34 : Vector<Int> channelIncrement(windows.size(), -1);
4320 :
4321 :
4322 17 : Int i = 0;
4323 34 : map<int, pair<int, int> > spwRanges=frequencySelection->getChannelRange ( measurementSets_p [msId()]) ;
4324 :
4325 34 : for (set<Int>::iterator j = windows.begin(); j != windows.end(); j++){
4326 :
4327 : //spectralWindow [i] = * j;
4328 17 : auto sel = spwRanges.find(*j);
4329 :
4330 17 : if(sel != spwRanges.end()){
4331 17 : spectralWindow.resize(i+1, True);
4332 17 : nChannels.resize(i+1,True);
4333 17 : firstChannel.resize(i+1, True);
4334 17 : channelIncrement.resize(i+1,True);
4335 17 : spectralWindow [i] = * j;
4336 17 : nChannels [i] = (sel->second).first;
4337 17 : firstChannel [i] =(sel->second).second;
4338 17 : channelIncrement[i] = 1;
4339 :
4340 17 : ++i;
4341 : }
4342 :
4343 : }
4344 :
4345 : return std::make_tuple(spectralWindow, nChannels, firstChannel,
4346 34 : channelIncrement);
4347 : }
4348 :
4349 :
4350 : VisibilityIteratorImpl2::ChannelInfo
4351 17 : VisibilityIteratorImpl2::getChannelInformation() const
4352 : {
4353 : const FrequencySelectionUsingChannels * frequencySelection =
4354 0 : dynamic_cast<const FrequencySelectionUsingChannels *>(
4355 17 : &frequencySelections_p->get(msId()));
4356 :
4357 17 : if (frequencySelection == 0) {
4358 17 : return getChannelInformationUsingFrequency();
4359 : }
4360 :
4361 0 : Vector<Int> spectralWindow;
4362 0 : Vector<Int> nChannels;
4363 0 : Vector<Int> firstChannel;
4364 0 : Vector<Int> channelIncrement;
4365 :
4366 0 : if (frequencySelection->empty()) {
4367 :
4368 : // No explicit selection, so everything is selected.
4369 :
4370 0 : casa::ms::SpectralWindows spectralWindows(& measurementSets_p[msId()]);
4371 :
4372 0 : spectralWindow.resize(spectralWindows.size());
4373 0 : nChannels.resize(spectralWindows.size());
4374 0 : firstChannel.resize(spectralWindows.size());
4375 0 : channelIncrement.resize(spectralWindows.size());
4376 :
4377 0 : Int i = 0;
4378 :
4379 0 : for(casa::ms::SpectralWindows::const_iterator s =
4380 0 : spectralWindows.begin();
4381 0 : s != spectralWindows.end();
4382 0 : s++) {
4383 :
4384 0 : spectralWindow(i) = s->id();
4385 0 : nChannels(i) = s->nChannels();
4386 0 : firstChannel(i) = 0;
4387 0 : channelIncrement(i) = 1;
4388 :
4389 0 : i++;
4390 : }
4391 : }
4392 : else {
4393 :
4394 : // Use the explicit channel-based selection to compute the result.
4395 :
4396 0 : spectralWindow.resize(frequencySelection->size());
4397 0 : nChannels.resize(frequencySelection->size());
4398 0 : firstChannel.resize(frequencySelection->size());
4399 0 : channelIncrement.resize(frequencySelection->size());
4400 :
4401 0 : Int i = 0;
4402 0 : for (FrequencySelectionUsingChannels::const_iterator j =
4403 0 : frequencySelection->begin();
4404 0 : j != frequencySelection->end();
4405 0 : ++j) {
4406 :
4407 0 : spectralWindow(i) = j->spectralWindow_p;
4408 0 : nChannels(i) = j->nChannels_p;
4409 0 : firstChannel(i) = j->firstChannel_p;
4410 0 : channelIncrement(i) = j->increment_p;
4411 :
4412 0 : i++;
4413 : }
4414 : }
4415 :
4416 : return std::make_tuple(spectralWindow, nChannels, firstChannel,
4417 0 : channelIncrement);
4418 : }
4419 :
4420 0 : Vector<casacore::Vector<Int> > VisibilityIteratorImpl2::getAllSelectedSpws() const{
4421 :
4422 0 : Vector<Vector<Int> > retval( frequencySelections_p->size());
4423 0 : for (uInt k=0; k < retval.nelements(); ++k){
4424 0 : std::set<Int> spw=(frequencySelections_p->get(k)).getSelectedWindows();
4425 0 : retval[k]=Vector<Int>(std::vector<Int>(spw.begin(), spw.end()));
4426 : }
4427 0 : return retval;
4428 : }
4429 :
4430 : void
4431 202765 : VisibilityIteratorImpl2::writeBackChanges(VisBuffer2 * vb)
4432 : {
4433 202765 : ThrowIf(!isWritable(), "This visibility iterator is not writable");
4434 :
4435 202765 : if (backWriters_p.empty()) {
4436 182 : initializeBackWriters();
4437 : }
4438 :
4439 405530 : VisBufferComponents2 dirtyComponents = vb->dirtyComponentsGet();
4440 :
4441 556198 : for (VisBufferComponents2::const_iterator dirtyComponent =
4442 202765 : dirtyComponents.begin();
4443 758963 : dirtyComponent != dirtyComponents.end();
4444 556198 : dirtyComponent++) {
4445 :
4446 556198 : ThrowIf(backWriters_p.find(* dirtyComponent) == backWriters_p.end(),
4447 : String::format(
4448 : "No writer defined for VisBuffer component %d",
4449 : *dirtyComponent));
4450 556198 : BackWriter * backWriter = backWriters_p[ * dirtyComponent];
4451 :
4452 : try {
4453 556198 : (* backWriter)(this, vb);
4454 0 : } catch(AipsError & e) {
4455 0 : Rethrow(
4456 : e,
4457 : String::format(
4458 : "Error while writing back VisBuffer component %d",
4459 : *dirtyComponent));
4460 : }
4461 : }
4462 202765 : }
4463 :
4464 : void
4465 182 : VisibilityIteratorImpl2::initializeBackWriters()
4466 : {
4467 0 : backWriters_p[VisBufferComponent2::FlagCube] =
4468 182 : makeBackWriter(& VisibilityIteratorImpl2::writeFlag, & VisBuffer2::flagCube);
4469 0 : backWriters_p[VisBufferComponent2::FlagRow] =
4470 182 : makeBackWriter(& VisibilityIteratorImpl2::writeFlagRow, & VisBuffer2::flagRow);
4471 0 : backWriters_p[VisBufferComponent2::FlagCategory] =
4472 182 : makeBackWriter(& VisibilityIteratorImpl2::writeFlagCategory, & VisBuffer2::flagCategory);
4473 0 : backWriters_p[VisBufferComponent2::Sigma] =
4474 182 : makeBackWriter(& VisibilityIteratorImpl2::writeSigma, & VisBuffer2::sigma);
4475 0 : backWriters_p[VisBufferComponent2::Weight] =
4476 182 : makeBackWriter(& VisibilityIteratorImpl2::writeWeight, & VisBuffer2::weight);
4477 0 : backWriters_p[VisBufferComponent2::WeightSpectrum] =
4478 182 : makeBackWriter(& VisibilityIteratorImpl2::writeWeightSpectrum, & VisBuffer2::weightSpectrum);
4479 0 : backWriters_p[VisBufferComponent2::SigmaSpectrum] =
4480 182 : makeBackWriter(& VisibilityIteratorImpl2::writeSigmaSpectrum, & VisBuffer2::sigmaSpectrum);
4481 :
4482 : // Now do the visibilities.
4483 :
4484 0 : backWriters_p[VisBufferComponent2::VisibilityCubeObserved] =
4485 182 : makeBackWriter(& VisibilityIteratorImpl2::writeVisObserved, & VisBuffer2::visCube);
4486 0 : backWriters_p[VisBufferComponent2::VisibilityCubeCorrected] =
4487 182 : makeBackWriter(& VisibilityIteratorImpl2::writeVisCorrected, & VisBuffer2::visCubeCorrected);
4488 546 : backWriters_p[VisBufferComponent2::VisibilityCubeModel] =
4489 182 : makeBackWriter(& VisibilityIteratorImpl2::writeVisModel, & VisBuffer2::visCubeModel);
4490 :
4491 182 : }
4492 :
4493 5434 : VisibilityIteratorImpl2::PendingChanges::PendingChanges()
4494 : : frequencySelections_p(0)
4495 : , frequencySelectionsPending_p(false)
4496 : , interval_p(Empty)
4497 5434 : , nRowBlocking_p(Empty)
4498 5434 : {}
4499 :
4500 10868 : VisibilityIteratorImpl2::PendingChanges::~PendingChanges()
4501 : {
4502 5434 : delete frequencySelections_p;
4503 5434 : }
4504 :
4505 : VisibilityIteratorImpl2::PendingChanges *
4506 0 : VisibilityIteratorImpl2::PendingChanges::clone() const
4507 : {
4508 0 : PendingChanges * theClone = new PendingChanges();
4509 :
4510 0 : theClone->frequencySelections_p =
4511 0 : new FrequencySelections(* frequencySelections_p);
4512 0 : theClone->frequencySelectionsPending_p = frequencySelectionsPending_p;
4513 0 : theClone->interval_p = interval_p;
4514 0 : theClone->nRowBlocking_p = nRowBlocking_p;
4515 :
4516 0 : return theClone;
4517 : }
4518 :
4519 :
4520 :
4521 : Bool
4522 3609479 : VisibilityIteratorImpl2::PendingChanges::empty() const
4523 : {
4524 10824913 : Bool result = frequencySelections_p == 0 &&
4525 7215434 : interval_p == Empty &&
4526 3605955 : nRowBlocking_p == Empty;
4527 :
4528 3609479 : return result;
4529 : }
4530 :
4531 : pair<Bool, FrequencySelections *>
4532 3605 : VisibilityIteratorImpl2::PendingChanges::popFrequencySelections()
4533 : {
4534 : // yields ownership
4535 3605 : FrequencySelections * result = frequencySelections_p;
4536 3605 : Bool wasPending = frequencySelectionsPending_p;
4537 :
4538 3605 : frequencySelections_p = 0;
4539 3605 : frequencySelectionsPending_p = false;
4540 :
4541 7210 : return make_pair(wasPending, result);
4542 : }
4543 :
4544 : pair<Bool, Double>
4545 3605 : VisibilityIteratorImpl2::PendingChanges::popInterval()
4546 : {
4547 3605 : pair<Bool,Double> result = make_pair(interval_p != Empty, interval_p);
4548 :
4549 3605 : interval_p = Empty;
4550 :
4551 3605 : return result;
4552 : }
4553 :
4554 : pair<Bool, Int>
4555 3605 : VisibilityIteratorImpl2::PendingChanges::popNRowBlocking()
4556 : {
4557 3605 : pair<Bool,Int> result = make_pair(nRowBlocking_p != Empty, nRowBlocking_p);
4558 :
4559 3605 : nRowBlocking_p = Empty;
4560 :
4561 3605 : return result;
4562 : }
4563 :
4564 : void
4565 3524 : VisibilityIteratorImpl2::PendingChanges::setFrequencySelections(
4566 : FrequencySelections * fs)
4567 : {
4568 : // takes ownership
4569 3524 : Assert(!frequencySelectionsPending_p);
4570 :
4571 3524 : frequencySelections_p = fs;
4572 3524 : frequencySelectionsPending_p = true;
4573 3524 : }
4574 :
4575 : void
4576 0 : VisibilityIteratorImpl2::PendingChanges::setInterval(Double interval)
4577 : {
4578 0 : Assert(interval_p == Empty);
4579 :
4580 0 : interval_p = interval;
4581 0 : }
4582 :
4583 : void
4584 824 : VisibilityIteratorImpl2::PendingChanges::setNRowBlocking(Int nRowBlocking)
4585 : {
4586 824 : Assert(nRowBlocking_p == Empty);
4587 :
4588 824 : nRowBlocking_p = nRowBlocking;
4589 824 : }
4590 :
4591 : bool
4592 449122 : VisibilityIteratorImpl2::fillFromVirtualModel(Cube <Complex> & value) const
4593 : {
4594 : // The model is virtual if there is no model column or there is a virtual
4595 : // model defined in the MS.
4596 :
4597 898244 : String modelKey = "definedmodel_field_" + String::toString(vb_p->fieldId());
4598 : Int sourceRow;
4599 898244 : Bool hasModelKey= !(modelDataGenerator_p == 0) &&
4600 449122 : modelDataGenerator_p->isModelDefinedI(vb_p->fieldId()(0), ms(),
4601 449122 : modelKey, sourceRow);
4602 :
4603 449122 : Bool isVirtual = hasModelKey || !(ms().tableDesc().isColumn("MODEL_DATA"));
4604 :
4605 449122 : if (isVirtual) {
4606 262562 : modelDataGenerator_p->init(*vb_p);
4607 :
4608 :
4609 : //////////This bit can be removed once version 1 is no longer read
4610 262562 : if(!(modelDataGenerator_p->isVersion2())){
4611 :
4612 262318 : auto field = vb_p->fieldId()(0);
4613 262318 : auto spw = vb_p->spectralWindows()(0);
4614 262318 : if (modelDataGenerator_p->hasModel(msId(),field , spw) == -1) {
4615 :
4616 : // If the model generator does not have a model for this(ms,field,
4617 : // spectralWindow) then try to add it.
4618 :
4619 262318 : if (hasModelKey) {
4620 :
4621 : // Read the record of model information and if found add it to
4622 : // the model generator.
4623 :
4624 0 : TableRecord modelRecord;
4625 0 : if (modelDataGenerator_p->getModelRecordI(
4626 0 : modelKey, modelRecord, ms())) {
4627 :
4628 0 : modelDataGenerator_p->addModel(
4629 0 : modelRecord, Vector<Int>(1, msId()), * vb_p);
4630 : }
4631 : }
4632 : }
4633 : }
4634 : /////////////////////////
4635 : // Now use the model data generator to fill in the model data.
4636 : // Temporarily make the VisBuffer writable, if it wasn't already.
4637 :
4638 262562 : Bool wasWritable = vb_p->setWritability(true);
4639 262562 : modelDataGenerator_p->getModelVis(const_cast <VisBuffer2 &>(* vb_p));
4640 262562 : vb_p->setWritability(wasWritable);
4641 :
4642 : // Put the model cube into the provided container. If that turns out to
4643 : // be the model component of the VIIs VB2 then this will be a no-op.
4644 :
4645 262562 : value = vb_p->visCubeModel();
4646 262562 : return true; // filled it
4647 : }
4648 :
4649 186560 : return false; // Did not fill
4650 : }
4651 :
4652 : //**********************************************************************
4653 : // Methods to access the subtables.
4654 : //**********************************************************************
4655 :
4656 : // Access to antenna subtable
4657 0 : const casacore::MSAntennaColumns& VisibilityIteratorImpl2::antennaSubtablecols() const
4658 : {
4659 0 : return msIter_p->msColumns().antenna();
4660 : }
4661 :
4662 : // Access to dataDescription subtable
4663 123 : const casacore::MSDataDescColumns& VisibilityIteratorImpl2::dataDescriptionSubtablecols() const
4664 : {
4665 123 : return msIter_p->msColumns().dataDescription();
4666 : }
4667 :
4668 : // Access to feed subtable
4669 0 : const casacore::MSFeedColumns& VisibilityIteratorImpl2::feedSubtablecols() const
4670 : {
4671 0 : return msIter_p->msColumns().feed();
4672 : }
4673 :
4674 : // Access to field subtable
4675 21 : const casacore::MSFieldColumns& VisibilityIteratorImpl2::fieldSubtablecols() const
4676 : {
4677 21 : return msIter_p->msColumns().field();
4678 : }
4679 :
4680 : // Access to flagCmd subtable
4681 0 : const casacore::MSFlagCmdColumns& VisibilityIteratorImpl2::flagCmdSubtablecols() const
4682 : {
4683 0 : return msIter_p->msColumns().flagCmd();
4684 : }
4685 :
4686 : // Access to history subtable
4687 0 : const casacore::MSHistoryColumns& VisibilityIteratorImpl2::historySubtablecols() const
4688 : {
4689 0 : return msIter_p->msColumns().history();
4690 : }
4691 :
4692 : // Access to observation subtable
4693 0 : const casacore::MSObservationColumns& VisibilityIteratorImpl2::observationSubtablecols() const
4694 : {
4695 0 : return msIter_p->msColumns().observation();
4696 : }
4697 :
4698 : // Access to pointing subtable
4699 0 : const casacore::MSPointingColumns& VisibilityIteratorImpl2::pointingSubtablecols() const
4700 : {
4701 0 : return msIter_p->msColumns().pointing();
4702 : }
4703 :
4704 : // Access to polarization subtable
4705 47 : const casacore::MSPolarizationColumns& VisibilityIteratorImpl2::polarizationSubtablecols() const
4706 : {
4707 47 : return msIter_p->msColumns().polarization();
4708 : }
4709 :
4710 : // Access to processor subtable
4711 0 : const casacore::MSProcessorColumns& VisibilityIteratorImpl2::processorSubtablecols() const
4712 : {
4713 0 : return msIter_p->msColumns().processor();
4714 : }
4715 :
4716 : // Access to spectralWindow subtable
4717 311 : const casacore::MSSpWindowColumns& VisibilityIteratorImpl2::spectralWindowSubtablecols() const
4718 : {
4719 311 : return msIter_p->msColumns().spectralWindow();
4720 : }
4721 :
4722 : // Access to state subtable
4723 40014 : const casacore::MSStateColumns& VisibilityIteratorImpl2::stateSubtablecols() const
4724 : {
4725 40014 : return msIter_p->msColumns().state();
4726 : }
4727 :
4728 : // Access to doppler subtable
4729 0 : const casacore::MSDopplerColumns& VisibilityIteratorImpl2::dopplerSubtablecols() const
4730 : {
4731 0 : return msIter_p->msColumns().doppler();
4732 : }
4733 :
4734 : // Access to freqOffset subtable
4735 0 : const casacore::MSFreqOffsetColumns& VisibilityIteratorImpl2::freqOffsetSubtablecols() const
4736 : {
4737 0 : return msIter_p->msColumns().freqOffset();
4738 : }
4739 :
4740 : // Access to source subtable
4741 0 : const casacore::MSSourceColumns& VisibilityIteratorImpl2::sourceSubtablecols() const
4742 : {
4743 0 : return msIter_p->msColumns().source();
4744 : }
4745 :
4746 : // Access to sysCal subtable
4747 0 : const casacore::MSSysCalColumns& VisibilityIteratorImpl2::sysCalSubtablecols() const
4748 : {
4749 0 : return msIter_p->msColumns().sysCal();
4750 : }
4751 :
4752 : // Access to weather subtable
4753 0 : const casacore::MSWeatherColumns& VisibilityIteratorImpl2::weatherSubtablecols() const
4754 : {
4755 0 : return msIter_p->msColumns().weather();
4756 : }
4757 :
4758 : } // end namespace vi
4759 :
4760 : } //# NAMESPACE CASA - END
|
