Line data Source code
1 : //# VisibilityIterator.cc: Step through MeasurementEquation by visibility
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
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: VisIterator.cc,v 19.15 2006/02/01 01:25:14 kgolap Exp $
27 :
28 : #include <msvis/MSVis/VisIterator.h>
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <stdcasa/UtilJ.h>
31 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
32 : #include <casacore/casa/Arrays/ArrayLogical.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/MaskedArray.h>
35 : #include <casacore/casa/Exceptions/Error.h>
36 : #include <casacore/casa/Utilities/Assert.h>
37 : #include <casacore/casa/Utilities/Sort.h>
38 : #include <casacore/ms/MeasurementSets/MSColumns.h>
39 : #include <casacore/casa/Quanta/MVTime.h>
40 : #include <casacore/tables/Tables/TableDesc.h>
41 : #include <casacore/tables/Tables/ColDescSet.h>
42 : #include <casacore/tables/Tables/TableRecord.h>
43 : #include <casacore/tables/DataMan/TiledColumnStMan.h>
44 : #include <casacore/tables/DataMan/TiledStManAccessor.h>
45 : #include <msvis/MSVis/VisibilityIteratorImpl.h>
46 : #include <cmath>
47 :
48 : using std::ceil;
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 : class ROVisIteratorImpl : public VisibilityIteratorReadImpl {
54 :
55 : friend class ROVisIterator;
56 : friend class VisIterator;
57 :
58 : public:
59 :
60 : // Default constructor - useful only to assign another iterator later
61 : ROVisIteratorImpl();
62 :
63 : // Construct from MS and a Block of MS column enums specifying the iteration
64 : // order, if none are specified, time iteration is implicit. An optional
65 : // timeInterval can be given to iterate through chunks of time. The default
66 : // interval of 0 groups all times together. Every 'chunk' of data contains
67 : // all data within a certain time interval (in seconds) and with identical
68 : // values of the other iteration columns (e.g. DATA_DESC_ID and FIELD_ID).
69 : // Using selectChannel(), a number of groups of channels can be requested.
70 : // At present the channel group iteration will always occur before the
71 : // interval iteration.
72 : ROVisIteratorImpl(ROVisIterator * rovi,
73 : const MeasurementSet& ms,
74 : const Block<Int>& sortColumns,
75 : Double timeInterval=0);
76 :
77 : // Copy construct. This calls the assignment operator.
78 : ROVisIteratorImpl(const ROVisIteratorImpl & other);
79 :
80 : // Destructor
81 : virtual ~ROVisIteratorImpl();
82 :
83 : // Assignment. Any attached VisBuffers are lost in the assign.
84 : ROVisIteratorImpl & operator=(const ROVisIteratorImpl &other);
85 :
86 : // Members
87 :
88 : // Advance iterator through data
89 : ROVisIteratorImpl & operator++(int);
90 : ROVisIteratorImpl & operator++();
91 :
92 : // Selected spws and channel counts
93 : virtual void allSelectedSpectralWindows (Vector<Int> & spws, Vector<Int> & nvischan);
94 :
95 : //
96 : virtual void lsrFrequency(const Int& spw, Vector<Double>& freq, Bool& convert, const Bool ignoreconv=false);
97 :
98 :
99 : // The following throws an exception, because this isn't the
100 : // language of channel selection in VisIterator
101 0 : virtual void getChannelSelection(Block< Vector<Int> >&,
102 : Block< Vector<Int> >&,
103 : Block< Vector<Int> >&,
104 : Block< Vector<Int> >&,
105 : Block< Vector<Int> >&)
106 0 : { throw(AipsError("ROVisIteratorImpl::getChannelSelection: you can't do that!")); };
107 :
108 : // Return channel numbers in selected VisSet spectrum
109 : // (i.e. disregarding possible selection on the iterator, but
110 : // including the selection set when creating the VisSet)
111 : Vector<Int>& channel(Vector<Int>& chan) const;
112 : Vector<Int>& chanIds(Vector<Int>& chanids) const;
113 : Vector<Int>& chanIds(Vector<Int>& chanids,Int spw) const;
114 :
115 : // Return selected correlation indices
116 : Vector<Int>& corrIds(Vector<Int>& corrids) const;
117 :
118 : // Return flag for each polarization, channel and row
119 : Cube<Bool>& flag(Cube<Bool>& flags) const;
120 :
121 : // Return current frequencies
122 : Vector<Double>& frequency(Vector<Double>& freq) const;
123 :
124 : // Return the correlation type (returns Stokes enums)
125 : Vector<Int>& corrType(Vector<Int>& corrTypes) const;
126 :
127 : // Return sigma matrix (pol-dep)
128 : Matrix<Float>& sigmaMat(Matrix<Float>& sigmat) const;
129 :
130 : // Return the visibilities as found in the MS, Cube(npol,nchan,nrow).
131 : Cube<Complex>& visibility(Cube<Complex>& vis,
132 : DataColumn whichOne) const;
133 : // Return weight matrix
134 : Matrix<Float>& weightMat(Matrix<Float>& wtmat) const;
135 :
136 : // Return weightspectrum (a weight for each corr & channel)
137 : Cube<Float>& weightSpectrum(Cube<Float>& wtsp) const;
138 :
139 : // Set up new chan/corr selection via Vector<Slice>
140 : void selectChannel(const Vector<Vector<Slice> >& chansel);
141 : void selectCorrelation(const Vector<Vector<Slice> >& corrsel);
142 :
143 : // Set up/return channel averaging bounds
144 : Vector<Matrix<Int> >& setChanAveBounds(Float factor, Vector<Matrix<Int> >& bounds);
145 :
146 : // Return number of chans/corrs per spw/pol
147 : Int numberChan(Int spw) const;
148 : Int numberCorr(Int pol) const;
149 :
150 : // Return the row ids as from the original root table. This is useful
151 : // to find correspondance between a given row in this iteration to the
152 : // original ms row
153 : Vector<rownr_t>& rowIds(Vector<rownr_t>& rowids) const;
154 :
155 : // Need to override this and not use getColArray
156 : Vector<RigidVector<Double,3> >& uvw(Vector<RigidVector<Double,3> >& uvwvec) const;
157 :
158 : protected:
159 :
160 : void setSelTable();
161 :
162 : virtual const Table attachTable() const;
163 :
164 : // update the DATA slicer
165 : void updateSlicer();
166 : // attach the column objects to the currently selected table
167 :
168 : // The ROVisibilityIterator version of this function sets the tile cache to 1
169 : // because of a feature in sliced data access that grows memory dramatically in
170 : // some cases. However, ROVisibilityIterator, because it uses
171 : // ArrayColumn::getColumn(Vector<Vector<Slice> >&), is (1/28/2011) incredibly
172 : // slow if the tile cache does not span all the selected channels, and it will
173 : // read orders of magnitude more data than it needs to. This sets the tile
174 : // cache to the minimum number of tiles required to span the selected channels.
175 : // Unlike ROVisibilityIterator, it does it for each hypercube, not just the
176 : // first one, and it does its work when the DDID has changed.
177 : void setTileCache();
178 :
179 : void getDataColumn(DataColumn whichOne, const Vector<Vector<Slice> >& slices,
180 : Cube<Complex>& data) const;
181 :
182 : // Column access functions
183 : void getCol(const ScalarColumn<Bool> &column, Vector<Bool> &array, Bool resize = false) const;
184 : void getCol(const ScalarColumn<Int> &column, Vector<Int> &array, Bool resize = false) const;
185 : void getCol(const ScalarColumn<Double> &column, Vector<Double> &array, Bool resize = false) const;
186 :
187 : void getCol(const ArrayColumn<Bool> &column, Array<Bool> &array, Bool resize = false) const;
188 : void getCol(const ArrayColumn<Float> &column, Array<Float> &array, Bool resize = false) const;
189 : void getCol(const ArrayColumn<Double> &column, Array<Double> &array, Bool resize = false) const;
190 : void getCol(const ArrayColumn<Complex> &column, Array<Complex> &array, Bool resize = false) const;
191 :
192 : void getCol(const ArrayColumn<Bool> &column, const Slicer &slicer, Array<Bool> &array, Bool resize = false) const;
193 : void getCol(const ArrayColumn<Float> &column, const Slicer &slicer, Array<Float> &array, Bool resize = false) const;
194 : void getCol(const ArrayColumn<Complex> &column, const Slicer &slicer, Array<Complex> &array, Bool resize = false) const;
195 :
196 :
197 : // New slicer supports multiple Slices in channel and correlation
198 :
199 : Vector<Matrix<Int> > chanAveBounds_p;
200 : Vector<Vector<Slice> > chanSlices_p;
201 : Vector<Vector<Slice> > corrSlices_p;
202 : Vector<Vector<Slice> > newSlicer_p;
203 : Vector<Vector<Slice> > newWtSlicer_p;
204 : Table selTable_p;
205 : Bool useNewSlicer_p;
206 :
207 : };
208 :
209 : class VisIteratorImpl : public ROVisIteratorImpl {
210 :
211 : friend class ROVisIterator;
212 : friend class VisIterator;
213 :
214 : public:
215 :
216 : // Default constructor - useful only to assign another iterator later
217 : VisIteratorImpl();
218 :
219 : // Construct from MS and a Block of MS column enums specifying the iteration
220 : // order, if none are specified, time iteration is implicit. An optional
221 : // timeInterval can be given to iterate through chunks of time. The default
222 : // interval of 0 groups all times together. Every 'chunk' of data contains
223 : // all data within a certain time interval (in seconds) and with identical
224 : // values of the other iteration columns (e.g. DATA_DESC_ID and FIELD_ID).
225 : // Using selectChannel(), a number of groups of channels can be requested.
226 : // At present the channel group iteration will always occur before the
227 : // interval iteration.
228 : VisIteratorImpl(ROVisIterator * rovi,
229 : const MeasurementSet& ms,
230 : const Block<Int>& sortColumns,
231 : Double timeInterval=0);
232 :
233 : // Copy construct. This calls the assignment operator.
234 : VisIteratorImpl(const ROVisIteratorImpl & other);
235 :
236 : // Destructor
237 : ~VisIteratorImpl();
238 :
239 : // Assignment. Any attached VisBuffers are lost in the assign.
240 : VisIteratorImpl & operator=(const VisIteratorImpl &other);
241 :
242 : // Members
243 :
244 : // Advance iterator through data
245 : VisIteratorImpl & operator++(int);
246 : VisIteratorImpl & operator++();
247 :
248 : // Return channel numbers in selected VisSet spectrum
249 : // (i.e. disregarding possible selection on the iterator, but
250 : // including the selection set when creating the VisSet)
251 : //Vector<Int>& channel(Vector<Int>& chan) const;
252 : //Vector<Int>& chanIds(Vector<Int>& chanids) const;
253 : //Vector<Int>& chanIds(Vector<Int>& chanids,Int spw) const;
254 :
255 : // Return selected correlation indices
256 : //Vector<Int>& corrIds(Vector<Int>& corrids) const;
257 :
258 : // Return flag for each polarization, channel and row
259 : //Cube<Bool>& flag(Cube<Bool>& flags) const;
260 :
261 : // Return current frequencies
262 : //Vector<Double>& frequency(Vector<Double>& freq) const;
263 :
264 : // Return the correlation type (returns Stokes enums)
265 : //Vector<Int>& corrType(Vector<Int>& corrTypes) const;
266 :
267 : // Return sigma matrix (pol-dep)
268 : //Matrix<Float>& sigmaMat(Matrix<Float>& sigmat) const;
269 :
270 : // Return the visibilities as found in the MS, Cube(npol,nchan,nrow).
271 : //Cube<Complex>& visibility(Cube<Complex>& vis,
272 : //DataColumn whichOne) const;
273 : // Return weight matrix
274 : //Matrix<Float>& weightMat(Matrix<Float>& wtmat) const;
275 :
276 : // Return weightspectrum (a weight for each corr & channel)
277 : //Cube<Float>& weightSpectrum(Cube<Float>& wtsp) const;
278 :
279 : // Set up new chan/corr selection via Vector<Slice>
280 : //void selectChannel(const Vector<Vector<Slice> >& chansel);
281 : //void selectCorrelation(const Vector<Vector<Slice> >& corrsel);
282 :
283 : // Set up/return channel averaging bounds
284 : //Vector<Matrix<Int> >& setChanAveBounds(Float factor, Vector<Matrix<Int> >& bounds);
285 :
286 : // Return number of chans/corrs per spw/pol
287 : //Int numberChan(Int spw) const;
288 : //Int numberCorr(Int pol) const;
289 :
290 : // Return the row ids as from the original root table. This is useful
291 : // to find correspondance between a given row in this iteration to the
292 : // original ms row
293 : //Vector<uInt>& rowIds(Vector<uInt>& rowids) const;
294 :
295 : // Need to override this and not use getColArray
296 : //Vector<RigidVector<Double,3> >& uvw(Vector<RigidVector<Double,3> >& uvwvec) const;
297 :
298 : // Set/modify the flag row column; dimension Vector(nrow)
299 : void setFlagRow(const Vector<Bool>& rowflags);
300 :
301 : // Set/modify the flags in the data.
302 : // This sets the flags as found in the MS, Cube(npol,nchan,nrow),
303 : // where nrow is the number of rows in the current iteration (given by
304 : // nRow()).
305 : void setFlag(const Cube<Bool>& flag);
306 :
307 : // Set/modify the visibilities
308 : // This sets the data as found in the MS, Cube(npol,nchan,nrow).
309 : void setVis(const Cube<Complex>& vis, DataColumn whichOne);
310 :
311 : // Set the visibility and flags, and interpolate from velocities if needed
312 : void setVisAndFlag(const Cube<Complex>& vis, const Cube<Bool>& flag,
313 : DataColumn whichOne);
314 :
315 : // Set/modify the weightMat
316 : void setWeightMat(const Matrix<Float>& wtmat);
317 :
318 : // Set/modify the weightSpectrum
319 : void setWeightSpectrum(const Cube<Float>& wtsp);
320 :
321 : protected:
322 :
323 : //void setSelTable();
324 :
325 : void attachColumns (const Table & t);
326 :
327 : // update the DATA slicer
328 : //void updateSlicer();
329 : // attach the column objects to the currently selected table
330 :
331 : // The ROVisibilityIterator version of this function sets the tile cache to 1
332 : // because of a feature in sliced data access that grows memory dramatically in
333 : // some cases. However, ROVisibilityIterator, because it uses
334 : // ArrayColumn::getColumn(Vector<Vector<Slice> >&), is (1/28/2011) incredibly
335 : // slow if the tile cache does not span all the selected channels, and it will
336 : // read orders of magnitude more data than it needs to. This sets the tile
337 : // cache to the minimum number of tiles required to span the selected channels.
338 : // Unlike ROVisibilityIterator, it does it for each hypercube, not just the
339 : // first one, and it does its work when the DDID has changed.
340 :
341 : void getDataColumn(DataColumn whichOne, const Vector<Vector<Slice> >& slices,
342 : Cube<Complex>& data) const;
343 :
344 : // Column access functions
345 : // void getCol(const ScalarColumn<Bool> &column, Vector<Bool> &array, Bool resize = false) const;
346 : // void getCol(const ScalarColumn<Int> &column, Vector<Int> &array, Bool resize = false) const;
347 : // void getCol(const ScalarColumn<Double> &column, Vector<Double> &array, Bool resize = false) const;
348 : //
349 : // void getCol(const ArrayColumn<Bool> &column, Array<Bool> &array, Bool resize = false) const;
350 : // void getCol(const ArrayColumn<Float> &column, Array<Float> &array, Bool resize = false) const;
351 : // void getCol(const ArrayColumn<Double> &column, Array<Double> &array, Bool resize = false) const;
352 : // void getCol(const ArrayColumn<Complex> &column, Array<Complex> &array, Bool resize = false) const;
353 : //
354 : // void getCol(const ArrayColumn<Bool> &column, const Slicer &slicer, Array<Bool> &array, Bool resize = false) const;
355 : // void getCol(const ArrayColumn<Float> &column, const Slicer &slicer, Array<Float> &array, Bool resize = false) const;
356 : // void getCol(const ArrayColumn<Complex> &column, const Slicer &slicer, Array<Complex> &array, Bool resize = false) const;
357 :
358 : void putDataColumn(DataColumn whichOne,
359 : const Cube<Complex>& data);
360 :
361 : void putDataColumn(DataColumn whichOne,
362 : const Vector<Vector<Slice> >& slices,
363 : const Cube<Complex>& data);
364 :
365 : // New slicer supports multiple Slices in channel and correlation
366 :
367 : // Vector<Matrix<Int> > chanAveBounds_p;
368 : // Vector<Vector<Slice> > chanSlices_p;
369 : // Vector<Vector<Slice> > corrSlices_p;
370 : // Vector<Vector<Slice> > newSlicer_p;
371 : // Vector<Vector<Slice> > newWtSlicer_p;
372 : // Table selTable_p;
373 : // Bool useNewSlicer_p;
374 :
375 : ArrayColumn<Complex> rwColVis_p;
376 : ArrayColumn<Float> rwColFloatVis_p;
377 : ArrayColumn<Complex> rwColModelVis_p;
378 : ArrayColumn<Complex> rwColCorrVis_p;
379 : ArrayColumn<Float> rwColWeight_p;
380 : ArrayColumn<Float> rwColWeightSpectrum_p;
381 : ArrayColumn<Float> rwColSigma_p;
382 : ArrayColumn<Bool> rwColFlag_p;
383 : ScalarColumn<Bool> rwColFlagRow_p;
384 : };
385 :
386 : // ********************
387 : // ********************
388 : // ** **
389 : // ** ROVisIterator **
390 : // ** **
391 : // ********************
392 : // ********************
393 :
394 0 : ROVisIterator::ROVisIterator()
395 : {
396 0 : }
397 :
398 3 : ROVisIterator::ROVisIterator(const MeasurementSet & ms, const Block<Int>& sortColumns, Double timeInterval)
399 3 : : ROVisibilityIterator (ms, sortColumns, timeInterval, Factory (this))
400 3 : {}
401 :
402 0 : ROVisIterator::ROVisIterator(const MeasurementSet & ms, const Block<Int>& sortColumns,
403 0 : Double timeInterval, const ROVisibilityIterator::Factory & factory)
404 0 : : ROVisibilityIterator (ms, sortColumns, timeInterval, factory)
405 0 : {}
406 :
407 :
408 0 : ROVisIterator::ROVisIterator(const ROVisIterator & other)
409 0 : : ROVisibilityIterator (other)
410 0 : {}
411 :
412 3 : ROVisIterator::~ROVisIterator()
413 3 : {}
414 :
415 : ROVisIterator &
416 0 : ROVisIterator::operator=(const ROVisIterator & other)
417 : {
418 0 : ROVisibilityIterator::operator= (other);
419 :
420 0 : return * this;
421 : }
422 :
423 : ROVisIterator &
424 0 : ROVisIterator::operator++(int)
425 : {
426 0 : advance();
427 0 : return * this;
428 : }
429 :
430 : ROVisIterator &
431 12 : ROVisIterator::operator++()
432 : {
433 12 : advance();
434 12 : return * this;
435 : }
436 :
437 : Vector<Int>&
438 0 : ROVisIterator::chanIds(Vector<Int>& chanids) const
439 : {
440 0 : return getReadImpl()->chanIds(chanids);
441 : }
442 :
443 : Vector<Int>&
444 0 : ROVisIterator::chanIds(Vector<Int>& chanids, Int spw) const
445 : {
446 0 : return getReadImpl()->chanIds(chanids, spw);
447 : }
448 :
449 : Vector<Int>&
450 0 : ROVisIterator::corrIds(Vector<Int>& corrids) const
451 : {
452 0 : return getReadImpl()->corrIds(corrids);
453 : }
454 :
455 :
456 : Vector<Int>&
457 12 : ROVisIterator::corrType(Vector<Int>& corrTypes) const
458 : {
459 12 : return getReadImpl()->corrType(corrTypes);
460 : }
461 :
462 : ROVisIteratorImpl *
463 18 : ROVisIterator::getReadImpl () const
464 : {
465 18 : ROVisIteratorImpl * impl = dynamic_cast<ROVisIteratorImpl *> (ROVisibilityIterator::getReadImpl ());
466 :
467 18 : Assert (impl != NULL);
468 :
469 18 : return impl;
470 : }
471 :
472 : void
473 3 : ROVisIterator::selectChannel(const Vector<Vector<Slice> >& chansel)
474 : {
475 3 : return getReadImpl()->selectChannel(chansel);
476 : }
477 :
478 : void
479 3 : ROVisIterator::selectCorrelation(const Vector<Vector<Slice> >& corrsel)
480 : {
481 3 : return getReadImpl()->selectCorrelation(corrsel);
482 : }
483 :
484 : Vector<Matrix<Int> >&
485 0 : ROVisIterator::setChanAveBounds(Float factor, Vector<Matrix<Int> >& bounds)
486 : {
487 0 : return getReadImpl()->setChanAveBounds(factor, bounds);
488 : }
489 :
490 : void
491 30 : ROVisIteratorImpl::setTileCache()
492 : {
493 : // Set the cache size when the DDID changes (as opposed to MS) to avoid
494 : // overreading in a case like:
495 : // hcubes: [2, 256], [4, 64]
496 : // tileshape: [4, 64]
497 : // spws (ddids): [4,64], [2, 256], [2, 256], [4,64], [2, 256], [4,64]
498 : // and spw = '0:1~7,1:1~7,2:100~200,3:20~40,4:200~230,5:40~50'
499 : //
500 : // For hypercube 0, spw 2 needs 3 tiles, but spw 4 only needs 1, AND the last
501 : // tile at that. So if hypercube 0 used a cache of 3 tiles throughout, every
502 : // read of 4:200~230 would likely also read the unwanted channels 0~127 of
503 : // the next row.
504 : //
505 30 : if (curStartRow_p != 0 && ! msIter_p.newDataDescriptionId()){
506 0 : return;
507 : }
508 :
509 30 : const MeasurementSet& thems = msIter_p.ms();
510 :
511 30 : if(thems.tableType() == Table::Memory){
512 0 : return;
513 : }
514 :
515 30 : const ColumnDescSet& cds=thems.tableDesc().columnDescSet();
516 :
517 : // Get the first row number for this DDID.
518 60 : Vector<rownr_t> rownums;
519 30 : rowIds(rownums);
520 30 : uInt startrow = rownums[0];
521 :
522 60 : Vector<String> columns(8);
523 : // complex
524 30 : columns(0)=MS::columnName(MS::DATA);
525 30 : columns(1)=MS::columnName(MS::CORRECTED_DATA);
526 30 : columns(2)=MS::columnName(MS::MODEL_DATA);
527 : // boolean
528 30 : columns(3)=MS::columnName(MS::FLAG);
529 : // float
530 30 : columns(4)=MS::columnName(MS::WEIGHT_SPECTRUM);
531 30 : columns(5)=MS::columnName(MS::WEIGHT);
532 30 : columns(6)=MS::columnName(MS::SIGMA);
533 : // double
534 30 : columns(7)=MS::columnName(MS::UVW);
535 : //
536 270 : for(uInt k = 0; k < columns.nelements(); ++k){
537 240 : if(cds.isDefined(columns(k))){
538 240 : const ColumnDesc& cdesc = cds[columns(k)];
539 240 : String dataManType="";
540 :
541 240 : dataManType = cdesc.dataManagerType();
542 : // We have to check WEIGHT_SPECTRUM as it tends to exist but not have
543 : // valid data.
544 270 : if(columns[k] == MS::columnName(MS::WEIGHT_SPECTRUM) &&
545 30 : !existsWeightSpectrum())
546 30 : dataManType="";
547 :
548 : // Sometimes columns may not contain anything yet
549 270 : if((columns[k]==MS::columnName(MS::DATA) && (columns_p.vis_p.isNull() ||
550 30 : !columns_p.vis_p.isDefined(0))) ||
551 240 : (columns[k]==MS::columnName(MS::MODEL_DATA) && (columns_p.modelVis_p.isNull() ||
552 30 : !columns_p.modelVis_p.isDefined(0))) ||
553 240 : (columns[k]==MS::columnName(MS::CORRECTED_DATA) && (columns_p.corrVis_p.isNull() ||
554 30 : !columns_p.corrVis_p.isDefined(0))) ||
555 240 : (columns[k]==MS::columnName(MS::FLAG) && (columns_p.flag_p.isNull() ||
556 30 : !columns_p.flag_p.isDefined(0))) ||
557 240 : (columns[k]==MS::columnName(MS::WEIGHT) && (columns_p.weight_p.isNull() ||
558 30 : !columns_p.weight_p.isDefined(0))) ||
559 240 : (columns[k]==MS::columnName(MS::SIGMA) && (columns_p.sigma_p.isNull() ||
560 510 : !columns_p.sigma_p.isDefined(0))) ||
561 240 : (columns[k]==MS::columnName(MS::UVW) && (columns_p.uvw_p.isNull() ||
562 30 : !columns_p.uvw_p.isDefined(0))) ){
563 0 : dataManType="";
564 : }
565 :
566 390 : if(dataManType.contains("Tiled") &&
567 390 : !String(cdesc.dataManagerGroup()).empty()){
568 : try {
569 : ROTiledStManAccessor tacc=ROTiledStManAccessor(thems,
570 300 : cdesc.dataManagerGroup());
571 :
572 : // This is for the data columns, WEIGHT_SPECTRUM and FLAG only.
573 270 : if((columns[k] != MS::columnName(MS::WEIGHT)) &&
574 120 : (columns[k] != MS::columnName(MS::UVW))){
575 : // Figure out how many tiles are needed to span the selected channels.
576 180 : const IPosition tileShape(tacc.tileShape(startrow));
577 180 : Vector<Int> ids;
578 90 : chanIds(ids);
579 90 : uInt startTile = ids[0] / tileShape[1];
580 90 : uInt endTile = ids[ids.nelements() - 1] / tileShape[1];
581 90 : uInt cachesize = endTile - startTile + 1;
582 :
583 : // and the selected correlations.
584 90 : corrIds(ids);
585 90 : startTile = ids[0] / tileShape[0];
586 90 : endTile = ids[ids.nelements() - 1] / tileShape[0];
587 90 : cachesize *= endTile - startTile + 1;
588 :
589 : // The cache needs to hold enough tiles to provide the data
590 : // for a single row for all channels and all correlations.
591 : // Determine the number of tiles in the 0 and 1 directions,
592 : // round them up to the nearest integer and take the product
593 : // of the two to yield the number of tiles needed.
594 : //
595 : // This could be reduced a little bit if necessary by using the
596 : // actual frequency selection since the current calculation uses
597 : // all of the channels in the row, thus reserves space in the cache
598 : // for tiles which contain no selected data (however the logic doesn't
599 : // appear to read those into the cache so there is no performance
600 : // penalty).
601 :
602 180 : const IPosition hShape(tacc.hypercubeShape(startrow));
603 :
604 90 : float nTiles0 = hShape [0] / (float) (tileShape [0]);
605 90 : float nTiles1 = hShape [1] / (float) (tileShape [1]);
606 :
607 90 : cachesize = (uInt) (ceil (nTiles0) * ceil (nTiles1));
608 :
609 90 : tacc.setCacheSize(startrow, cachesize);
610 : }
611 : else
612 60 : tacc.setCacheSize(startrow, 1);
613 : }
614 0 : catch (AipsError x) {
615 : //It failed so leave the caching as is.
616 0 : continue;
617 : }
618 : }
619 : }
620 : }
621 : }
622 :
623 :
624 : // Selected spws and channel counts
625 0 : void ROVisIterator::allSelectedSpectralWindows (Vector<Int> & spws, Vector<Int> & nvischan)
626 : {
627 0 : getReadImpl()->allSelectedSpectralWindows(spws,nvischan);
628 0 : }
629 :
630 0 : void ROVisIterator::lsrFrequency(const Int& spw, Vector<Double>& freq, Bool& convert, const Bool ignoreconv)
631 : {
632 0 : getReadImpl()->lsrFrequency(spw,freq,convert,ignoreconv);
633 0 : }
634 :
635 :
636 : Int
637 0 : ROVisIterator::numberChan(Int spw) const
638 : {
639 0 : return getReadImpl()->numberChan(spw);
640 : }
641 :
642 : Int
643 0 : ROVisIterator::numberCorr(Int pol) const
644 : {
645 0 : return getReadImpl()->numberCorr(pol);
646 : }
647 :
648 : void
649 0 : ROVisIterator::getDataColumn(DataColumn whichOne, const Vector<Vector<Slice> >& slices, Cube<Complex>& data) const
650 : {
651 0 : getReadImpl()->getDataColumn(whichOne, slices, data);
652 0 : }
653 :
654 : VisibilityIteratorReadImpl *
655 3 : ROVisIterator::Factory::operator() (const asyncio::PrefetchColumns * /*prefetchColumns*/,
656 : const Block<MeasurementSet>& mss,
657 : const Block<Int>& sortColumns,
658 : const Bool /*addDefaultSortCols*/,
659 : Double timeInterval) const
660 : {
661 3 : return new ROVisIteratorImpl (vi_p, mss[0], sortColumns, timeInterval);
662 : }
663 :
664 : // *************************
665 : // *************************
666 : // ** **
667 : // ** ROVisIteratorImpl **
668 : // ** **
669 : // *************************
670 : // *************************
671 :
672 0 : ROVisIteratorImpl::ROVisIteratorImpl() {}
673 :
674 3 : ROVisIteratorImpl::ROVisIteratorImpl(ROVisIterator * rovi,
675 : const MeasurementSet &ms,
676 : const Block<Int>& sortColumns,
677 3 : Double timeInterval)
678 3 : : VisibilityIteratorReadImpl(rovi, Block<MeasurementSet> (1, ms), sortColumns, true, timeInterval),
679 6 : useNewSlicer_p(false)
680 : {
681 :
682 : // Initialize multi-slicers with empty slices
683 3 : chanSlices_p.resize(numberSpw());
684 3 : chanSlices_p.set(Vector<Slice>());
685 3 : corrSlices_p.resize(numberPol());
686 3 : corrSlices_p.set(Vector<Slice>());
687 :
688 3 : }
689 :
690 0 : ROVisIteratorImpl::ROVisIteratorImpl(const ROVisIteratorImpl& other)
691 0 : : VisibilityIteratorReadImpl (other)
692 : {
693 0 : operator=(other);
694 0 : }
695 :
696 6 : ROVisIteratorImpl::~ROVisIteratorImpl() {}
697 :
698 : ROVisIteratorImpl&
699 0 : ROVisIteratorImpl::operator=(const ROVisIteratorImpl& other)
700 : {
701 0 : if (this!=&other) {
702 :
703 0 : VisibilityIteratorReadImpl::operator=(other);
704 :
705 0 : chanAveBounds_p = other.chanAveBounds_p;
706 0 : chanSlices_p = other.chanSlices_p;
707 0 : corrSlices_p = other.corrSlices_p;
708 0 : newSlicer_p = other.newSlicer_p;
709 0 : newWtSlicer_p = other.newWtSlicer_p;
710 0 : selTable_p = other.selTable_p;
711 0 : useNewSlicer_p = other.useNewSlicer_p;
712 :
713 : }
714 :
715 0 : return *this;
716 : }
717 :
718 0 : ROVisIteratorImpl & ROVisIteratorImpl::operator++(int)
719 : {
720 0 : if (! more ()){
721 0 : advance();
722 : }
723 0 : return *this;
724 : }
725 :
726 0 : ROVisIteratorImpl & ROVisIteratorImpl::operator++()
727 : {
728 0 : if (! more ()){
729 0 : advance();
730 : }
731 0 : return *this;
732 : }
733 :
734 :
735 :
736 21 : void ROVisIteratorImpl::updateSlicer()
737 : {
738 21 : useNewSlicer_p=true;
739 :
740 : // cout << "Using new slicer!..." << flush;
741 :
742 21 : newSlicer_p.resize(2);
743 21 : newWtSlicer_p.resize(1);
744 :
745 : // cout << "newSlicer_p.shape() = " << newSlicer_p.shape() << endl;
746 :
747 21 : useSlicer_p=false;
748 :
749 : // Refer to correct slices for the current spwid/polid
750 21 : newSlicer_p(0).reference(corrSlices_p(this->polarizationId()));
751 21 : newSlicer_p(1).reference(chanSlices_p(this->spectralWindow()));
752 :
753 21 : newWtSlicer_p(0).reference(corrSlices_p(this->polarizationId()));
754 :
755 21 : setTileCache();
756 21 : }
757 :
758 : //
759 :
760 : // Selected spws and channel counts
761 0 : void ROVisIteratorImpl::allSelectedSpectralWindows (Vector<Int> & spws, Vector<Int> & nvischan) {
762 :
763 0 : Vector<Int> ddids;
764 0 : msColumns().dataDescId().getColumn(ddids);
765 0 : Int ndd=genSort(ddids, Sort::Ascending, (Sort::QuickSort | Sort::NoDuplicates));
766 0 : ddids.resize(ndd,true);
767 : // cout << "ddids = " << ddids << endl;
768 0 : Vector<Int> spwperdd;
769 0 : msColumns().dataDescription().spectralWindowId().getColumn(spwperdd);
770 : // cout << "spwperdd = " << spwperdd << endl;
771 0 : spws.resize(ndd);
772 0 : nvischan.resize(numberSpw());
773 0 : nvischan.set(-1);
774 0 : for (Int idd=0;idd<ndd;++idd) {
775 0 : spws(idd)=spwperdd(ddids(idd));
776 0 : nvischan(spws(idd))=this->numberChan(spws(idd));
777 : }
778 :
779 0 : }
780 :
781 0 : void ROVisIteratorImpl::lsrFrequency(const Int& spw, Vector<Double>& freq,
782 : Bool& convert, const Bool ignoreconv) {
783 :
784 :
785 0 : MFrequency::Types obsMFreqType = (MFrequency::Types) msIter_p.msColumns ().spectralWindow ().measFreqRef ()(spw);
786 0 : convert = obsMFreqType != MFrequency::LSRK; // make this parameter write-only
787 0 : if (ignoreconv) convert=false; // override if user says so
788 :
789 : // Set up the frequency converter
790 0 : MEpoch ep;
791 0 : ROScalarMeasColumn<MEpoch>(msIter_p.table (), MS::columnName (MS::TIME)).get (curStartRow_p, ep); // Setting epoch to iteration's first one
792 0 : MPosition obsPos = msIter_p.telescopePosition ();
793 0 : MDirection dir = msIter_p.phaseCenter ();
794 0 : MeasFrame frame (ep, obsPos, dir);
795 : MFrequency::Convert tolsr (obsMFreqType,
796 0 : MFrequency::Ref (MFrequency::LSRK, frame));
797 :
798 :
799 : // Get the requested spw's nominal frequencies (all of them)
800 : // (we will index these with channel ids later)
801 0 : Vector<Double> chanFreq(0);
802 0 : chanFreq = msIter_p.msColumns ().spectralWindow ().chanFreq ()(spw);
803 :
804 : // The selected channel ids for the requested spw
805 0 : Vector<Int> chans(0);
806 0 : this->chanIds(chans,spw);
807 0 : Int nchan=chans.nelements();
808 :
809 : // Create output frequencies
810 0 : freq.resize(nchan);
811 0 : for (Int ich=0;ich<nchan;++ich) {
812 0 : if (convert)
813 0 : freq[ich]=tolsr(chanFreq(chans(ich))).getValue().getValue();
814 : else
815 0 : freq[ich]=chanFreq(chans(ich));
816 :
817 : }
818 :
819 0 : }
820 :
821 :
822 : // (Alternative syntax for ROVisIter::chanIds)
823 24 : Vector<Int>& ROVisIteratorImpl::channel(Vector<Int>& chan) const
824 24 : { return this->chanIds(chan); }
825 :
826 :
827 114 : Vector<Int>& ROVisIteratorImpl::chanIds(Vector<Int>& chanids) const
828 : {
829 114 : return chanIds(chanids,msIter_p.spectralWindowId());
830 : }
831 :
832 114 : Vector<Int>& ROVisIteratorImpl::chanIds(Vector<Int>& chanids, Int spw) const
833 : {
834 :
835 114 : chanids.resize(this->numberChan(spw));
836 :
837 : // if non-trivial slices available
838 114 : if (chanSlices_p(spw).nelements() > 0 ) {
839 :
840 192 : Vector<Slice> slices(chanSlices_p(spw));
841 96 : Int nslices=slices.nelements();
842 :
843 96 : Int ich0=0;
844 192 : for (Int i=0;i<nslices;++i) {
845 96 : Int nchan=slices(i).length();
846 96 : Int start=slices(i).start();
847 480 : for (Int ich=0;ich<nchan;++ich,++ich0)
848 384 : chanids(ich0)=start+ich;
849 : }
850 : }
851 : else {
852 : // all channels selected...
853 18 : indgen(chanids);
854 : }
855 114 : return chanids;
856 : }
857 :
858 30 : void ROVisIteratorImpl::setSelTable()
859 : {
860 30 : VisibilityIteratorReadImpl::setSelTable();
861 :
862 : // The following code (which uses Table::operator() to create
863 : // a new RefTable) is computationally expensive. This could
864 : // be optimized by using the same method as in the
865 : // VisibilityIterator base class (which is to not create
866 : // RefTables but instead access the table column directly in
867 : // getReadImpl()->msIter_p.table() using VisibilityIterator::getReadImpl()->selRows_p).
868 :
869 : // Doing so would mean replacing calls like
870 : // colSigma.getColumn(newWtSlicer_p,sigmat,true);
871 : // with
872 : // colSigma.getColumnCells(getReadImpl()->selRows_p,newWtSlicer_p,sigmat,true);
873 : //
874 : // However the ArrayColumn class does allow passing both a
875 : // Vector<Vector<Slice> > and a RefRows parameter to get-/putColumn.
876 : // Such put/get functions must be first implemented.
877 :
878 30 : Vector<rownr_t> rows(curNumRow_p);
879 30 : indgen(rows,curStartRow_p);
880 30 : selTable_p=msIter_p.table()(rows);
881 30 : attachColumns(attachTable());
882 30 : }
883 :
884 : // Return native correlation _indices_
885 102 : Vector<Int>& ROVisIteratorImpl::corrIds(Vector<Int>& corrids) const
886 : {
887 102 : Int pol = msIter_p.polarizationId();
888 :
889 102 : corrids.resize(this->numberCorr(pol));
890 :
891 102 : Vector<Slice> slices(corrSlices_p(pol));
892 102 : Int nslices=slices.nelements();
893 :
894 : // if non-trivial slices available
895 102 : if (nslices > 0 ) {
896 :
897 0 : Int icor0=0;
898 0 : for (Int i=0;i<nslices;++i) {
899 0 : Int ncorr=slices(i).length();
900 0 : Int start=slices(i).start();
901 0 : for (Int icor=0;icor<ncorr;++icor,++icor0)
902 0 : corrids(icor0)=start+icor;
903 : }
904 : }
905 : else {
906 : // all corrs selected...
907 102 : indgen(corrids);
908 : }
909 :
910 204 : return corrids;
911 : }
912 :
913 : Vector<Int>&
914 12 : ROVisIteratorImpl::corrType(Vector<Int>& corrTypes) const
915 : {
916 :
917 : // Get the nominal corrType list
918 12 : Int polId = msIter_p.polarizationId();
919 24 : Vector<Int> nomCorrTypes;
920 12 : msIter_p.msColumns().polarization().corrType().get(polId,nomCorrTypes,true);
921 :
922 : // Get the 0-based corr indices
923 12 : Vector<Int> corrids;
924 12 : corrIds(corrids);
925 :
926 : // Set the corrType values by the corrids
927 12 : Int nCor=corrids.nelements();
928 12 : corrTypes.resize(nCor);
929 60 : for (Int icor=0;icor<nCor;++icor)
930 48 : corrTypes[icor]=nomCorrTypes[corrids[icor]];
931 :
932 24 : return corrTypes;
933 : }
934 :
935 12 : Cube<Bool>& ROVisIteratorImpl::flag(Cube<Bool>& flags) const
936 : {
937 12 : if (useNewSlicer_p){
938 12 : columns_p.flag_p.getColumn(newSlicer_p,flags,true);
939 : }
940 : else {
941 0 : columns_p.flag_p.getColumn(flags,true);
942 : }
943 12 : return flags;
944 : }
945 :
946 12 : Vector<Double>& ROVisIteratorImpl::frequency(Vector<Double>& freq) const
947 : {
948 :
949 : // We need to change some internals transparently
950 : //ROVisIteratorImpl* self = const_cast<ROVisIteratorImpl*>(this);
951 :
952 12 : if (! cache_p.freqCacheOK_p) {
953 12 : cache_p.freqCacheOK_p=true;
954 12 : const Vector<Double>& chanFreq=msIter_p.frequency();
955 24 : Vector<Int> chan;
956 12 : channel(chan);
957 12 : Int nch=chan.nelements();
958 12 : cache_p.frequency_p.resize(nch);
959 :
960 60 : for (Int i=0;i<nch;++i){
961 48 : cache_p.frequency_p(i)=chanFreq(chan(i));
962 : }
963 : }
964 : // Always copy to output
965 12 : freq.resize(cache_p.frequency_p.nelements());
966 12 : freq=cache_p.frequency_p;
967 :
968 12 : return freq;
969 : }
970 :
971 : Cube<Complex>&
972 36 : ROVisIteratorImpl::visibility(Cube<Complex>& vis, DataColumn whichOne) const
973 : {
974 36 : if (useNewSlicer_p){
975 36 : getDataColumn(whichOne,newSlicer_p,vis);
976 : }
977 : else {
978 0 : VisibilityIteratorReadImpl::getDataColumn(whichOne,vis);
979 : }
980 36 : return vis;
981 : }
982 :
983 36 : void ROVisIteratorImpl::getDataColumn(DataColumn whichOne,
984 : const Vector<Vector<Slice> >& slices,
985 : Cube<Complex>& data) const
986 : {
987 : // Return the visibility (observed, model or corrected);
988 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
989 36 : switch (whichOne) {
990 12 : case ROVisIterator::Observed:
991 12 : if (floatDataFound_p) {
992 0 : Cube<Float> dataFloat;
993 0 : columns_p.floatVis_p.getColumn(slices,dataFloat,true);
994 0 : data.resize(dataFloat.shape());
995 0 : convertArray(data,dataFloat);
996 : } else {
997 12 : columns_p.vis_p.getColumn(slices,data,true);
998 : };
999 12 : break;
1000 12 : case ROVisIterator::Corrected:
1001 12 : columns_p.corrVis_p.getColumn(slices,data,true);
1002 12 : break;
1003 12 : case ROVisIterator::Model:
1004 12 : columns_p.modelVis_p.getColumn(slices,data,true);
1005 12 : break;
1006 : };
1007 :
1008 36 : };
1009 :
1010 12 : Matrix<Float>& ROVisIteratorImpl::sigmaMat(Matrix<Float>& sigmat) const
1011 : {
1012 12 : if (useNewSlicer_p) columns_p.sigma_p.getColumn(newWtSlicer_p,sigmat,true);
1013 : else {
1014 0 : sigmat.resize(nPol_p,curNumRow_p);
1015 0 : columns_p.sigma_p.getColumn(sigmat);
1016 : }
1017 12 : return sigmat;
1018 : }
1019 :
1020 12 : Matrix<Float>& ROVisIteratorImpl::weightMat(Matrix<Float>& wtmat) const
1021 : {
1022 12 : if (useNewSlicer_p) columns_p.weight_p.getColumn(newWtSlicer_p,wtmat,true);
1023 : else {
1024 0 : wtmat.resize(nPol_p,curNumRow_p);
1025 0 : columns_p.weight_p.getColumn(wtmat);
1026 : }
1027 12 : return wtmat;
1028 : }
1029 :
1030 0 : Cube<Float>& ROVisIteratorImpl::weightSpectrum(Cube<Float>& wtsp) const
1031 : {
1032 0 : if (this->existsWeightSpectrum()) {
1033 0 : if (useNewSlicer_p) columns_p.weightSpectrum_p.getColumn(newSlicer_p,wtsp,true);
1034 : else {
1035 0 : wtsp.resize(nPol_p,nChan_p,curNumRow_p);
1036 0 : columns_p.weightSpectrum_p.getColumn(wtsp);
1037 : }
1038 : } else {
1039 0 : wtsp.resize(0,0,0);
1040 : }
1041 0 : return wtsp;
1042 : }
1043 :
1044 3 : void ROVisIteratorImpl::selectChannel(const Vector<Vector<Slice> >& chansel) {
1045 : // cout << "selectChannel!...vss..." << flush;
1046 :
1047 3 : if (chansel.nelements() != uInt(numberSpw()))
1048 0 : throw(AipsError("Specified channel slices has incorrect number of spws."));
1049 :
1050 3 : chanSlices_p.resize(numberSpw(),false);
1051 3 : chanSlices_p=chansel;
1052 :
1053 : // Enforce use of the new slicer downstream
1054 3 : useNewSlicer_p=true;
1055 :
1056 : // cout << "done." << endl;
1057 :
1058 3 : }
1059 :
1060 3 : void ROVisIteratorImpl::selectCorrelation(const Vector<Vector<Slice> >& corrsel) {
1061 :
1062 : // cout << "selectCorrelation!...vvs..." << flush;
1063 :
1064 3 : if (corrsel.nelements() != uInt(numberPol()))
1065 0 : throw(AipsError("Specified correlation slices has incorrect number of polIds."));
1066 :
1067 3 : corrSlices_p.resize(numberPol(),false);
1068 3 : corrSlices_p=corrsel;
1069 :
1070 : // Enforce use of the new slicer downstream
1071 3 : useNewSlicer_p=true;
1072 :
1073 : // cout << "done." << endl;
1074 :
1075 3 : }
1076 :
1077 0 : Vector<Matrix<Int> >& ROVisIteratorImpl::setChanAveBounds(Float factor,
1078 : Vector<Matrix<Int> >& bounds)
1079 : {
1080 0 : if(!useNewSlicer_p)
1081 0 : throw(AipsError("Help!"));
1082 :
1083 : // For every spw...
1084 0 : bounds.resize(numberSpw());
1085 :
1086 : // If factor greater than zero, fill the bounds non-trivially
1087 0 : if (factor > 0) {
1088 :
1089 : // Decipher averaging factor
1090 0 : Int width = 1;
1091 0 : if(factor > 1.0)
1092 0 : width = Int(factor); // factor supplied in channel units
1093 :
1094 : // Calculate bounds for each spw
1095 0 : for(Int ispw = 0; ispw < numberSpw(); ++ispw){
1096 :
1097 : // Number of SELECTED channels PRIOR to averaging
1098 0 : Int nChan0 = numberChan(ispw);
1099 :
1100 : // factor might have been supplied as a fraction;
1101 : // width is then potentially spw-dependent
1102 0 : if(factor <= 1.0)
1103 0 : width = Int(factor * Float(nChan0));
1104 :
1105 : // Get the selected channel list
1106 0 : Vector<Int> chans;
1107 0 : chanIds(chans, ispw);
1108 :
1109 : // The nominal number of output channels depends on the full
1110 : // range of channels selected (not the number of them)
1111 : // (will be revised later, if nec.)
1112 0 : Int nChanOut0 = 1 + (chans[nChan0 - 1] - chans[0]) / width;
1113 :
1114 : // Initialize the bounds container for this spw
1115 0 : Matrix<Int>& currBounds(bounds[ispw]);
1116 0 : currBounds.resize(nChanOut0, 4);
1117 : //currBounds.set(0);
1118 :
1119 0 : Int outChan = 0;
1120 0 : Int firstchan = chans[0];
1121 0 : Int lastchan = chans[nChan0 - 1];
1122 :
1123 : // Index of input channel in SELECTED list, i.e.
1124 : // ich = vi.chanIds(chanids, spw)[selchanind].
1125 0 : uInt selchanind = 0;
1126 :
1127 0 : for(Int ich = firstchan; ich <= lastchan; ich += width){
1128 0 : Int w = width;
1129 :
1130 : // Use a clamped width in case (lastchan - firstchan + 1) % width != 0.
1131 0 : if(ich + w - 1 > lastchan)
1132 0 : w = lastchan + 1 - ich;
1133 :
1134 : // The start and end in MS channel #s.
1135 0 : currBounds(outChan, 0) = ich;
1136 0 : currBounds(outChan, 1) = ich + w - 1;
1137 :
1138 : // The start and end in selected reckoning.
1139 0 : currBounds(outChan, 2) = selchanind;
1140 0 : selchanind += w;
1141 0 : currBounds(outChan, 3) = selchanind - 1;
1142 :
1143 : // for(uInt ii = 0; ii < 4; ++ii)
1144 : // cerr << "currBounds(" << outChan << ", " << ii << ") = "
1145 : // << currBounds(outChan, ii) << endl;
1146 0 : ++outChan;
1147 : }
1148 : } // ispw
1149 :
1150 : } // factor > 0
1151 :
1152 : // Internal reference (needed?)
1153 0 : chanAveBounds_p.reference(bounds);
1154 :
1155 : // Return the bounds Vector reference
1156 0 : return bounds;
1157 : }
1158 :
1159 : // Vector<Matrix<Int> >& ROVisIteratorImpl::setChanAveBounds(Float factor,
1160 : // Vector<Matrix<Int> >& bounds) {
1161 :
1162 : // if (!useNewSlicer_p) throw(AipsError("Help!"));
1163 :
1164 : // // For every spw...
1165 : // bounds.resize(numberSpw());
1166 :
1167 : // // If factor greater than zero, fill the bounds non-trivially
1168 : // if (factor>0) {
1169 :
1170 : // // Decipher averaging factor
1171 : // Int width(1);
1172 : // if (factor>1.0) width=Int(factor); // factor supplied in channel units
1173 :
1174 : // // Calculate bounds for each spw
1175 : // for (Int ispw=0;ispw<numberSpw();++ispw) {
1176 :
1177 : // // Number of SELECTED channels PRIOR to averaging
1178 : // Int nChan0=numberChan(ispw);
1179 :
1180 : // // factor might have been supplied in factional units;
1181 : // // width is then potentially spw-dependent
1182 : // if (factor<=1.0)
1183 : // width=Int(factor*Float(nChan0));
1184 :
1185 : // // Get the selected channel list
1186 : // Vector<Int> chans;
1187 : // chanIds(chans,ispw);
1188 :
1189 : // // The nominal number of output channels depends on the full
1190 : // // range of channels selected (not the number of them)
1191 : // // (will be revised later, if nec.)
1192 : // Int nChanOut0((chans(nChan0-1)-chans(0)+1+width)/width);
1193 :
1194 : // // Initialize the bounds container for this spw
1195 : // Matrix<Int>& currBounds(bounds(ispw));
1196 : // currBounds.resize(nChanOut0,4);
1197 : // currBounds.set(0);
1198 :
1199 : // // Actual output channel count; initialization
1200 : // Int nChanOut=1;
1201 : // Int lo(chans(0));
1202 : // currBounds(0,0)=lo;
1203 : // currBounds(0,2)=0;
1204 : // for (Int ich=0;ich<nChan0;++ich)
1205 : // if ( (chans(ich)-lo+1)>width ) {
1206 : // currBounds(nChanOut-1,1)=chans(ich-1); // end of previous
1207 : // currBounds(nChanOut-1,3)=ich-1;
1208 : // lo=currBounds(nChanOut,0)=chans(ich); // start of next
1209 : // currBounds(nChanOut,2)=ich; // start of next
1210 : // ++nChanOut;
1211 : // }
1212 : // currBounds(nChanOut-1,1)=chans(nChan0-1); // end of final set
1213 : // currBounds(nChanOut-1,3)=(nChan0-1); // end of final set
1214 :
1215 : // // for(uInt ii = 0; ii < 4; ++ii)
1216 : // // cerr << "currBounds(" << nChanOut - 1 << ", " << ii << ") = "
1217 : // // << currBounds(nChanOut - 1, ii) << endl;
1218 :
1219 : // // contract bounds container, if necessary
1220 : // if (nChanOut<nChanOut0)
1221 : // currBounds.resize(nChanOut,4,true);
1222 :
1223 : // } // ispw
1224 :
1225 :
1226 : // } // factor > 0
1227 :
1228 : // // Internal reference (needed?)
1229 : // chanAveBounds_p.reference(bounds);
1230 :
1231 : // // Return the bounds Vector reference
1232 : // return bounds;
1233 :
1234 : // }
1235 :
1236 114 : Int ROVisIteratorImpl::numberChan(Int spw) const {
1237 :
1238 : // Nominally all channels this spw
1239 114 : Int nchan=msColumns().spectralWindow().numChan()(spw);
1240 :
1241 114 : if (useNewSlicer_p) {
1242 114 : Int nslices=chanSlices_p(spw).nelements();
1243 114 : if (nslices > 0 ) {
1244 96 : nchan=0;
1245 192 : for (Int isl=0;isl<nslices;++isl)
1246 96 : nchan+=chanSlices_p(spw)(isl).length();
1247 : }
1248 : }
1249 :
1250 114 : return nchan;
1251 :
1252 : }
1253 :
1254 :
1255 102 : Int ROVisIteratorImpl::numberCorr(Int pol) const {
1256 :
1257 : // Nominally all correlations this pol
1258 102 : Int ncorr=msColumns().polarization().numCorr()(pol);
1259 :
1260 102 : if (useNewSlicer_p) {
1261 102 : Int nslices=corrSlices_p(pol).nelements();
1262 102 : if (nslices > 0 ) {
1263 : // Accumulate from slice lengths
1264 0 : ncorr=0;
1265 0 : for (Int isl=0;isl<nslices;++isl)
1266 0 : ncorr+=corrSlices_p(pol)(isl).length();
1267 : }
1268 : }
1269 :
1270 102 : return ncorr;
1271 :
1272 : }
1273 :
1274 12 : void ROVisIteratorImpl::getCol(const ScalarColumn<Bool> &column, Vector<Bool> &array, Bool resize) const
1275 : {
1276 12 : column.getColumn(array, resize);
1277 12 : }
1278 :
1279 0 : void ROVisIteratorImpl::getCol(const ScalarColumn<Int> &column, Vector<Int> &array, Bool resize) const
1280 : {
1281 0 : column.getColumn(array, resize);
1282 0 : }
1283 :
1284 0 : void ROVisIteratorImpl::getCol(const ScalarColumn<Double> &column, Vector<Double> &array, Bool resize) const
1285 : {
1286 0 : column.getColumn(array, resize);
1287 0 : }
1288 :
1289 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Bool> &column, Array<Bool> &array, Bool resize) const
1290 : {
1291 0 : column.getColumn(array, resize);
1292 0 : }
1293 :
1294 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Float> &column, Array<Float> &array, Bool resize) const
1295 : {
1296 0 : column.getColumn(array, resize);
1297 0 : }
1298 :
1299 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Double> &column, Array<Double> &array, Bool resize) const
1300 : {
1301 0 : column.getColumn(array, resize);
1302 0 : }
1303 :
1304 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Complex> &column, Array<Complex> &array, Bool resize) const
1305 : {
1306 0 : column.getColumn(array, resize);
1307 0 : }
1308 :
1309 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Bool> &column, const Slicer &slicer, Array<Bool> &array, Bool resize) const
1310 : {
1311 0 : column.getColumn(slicer, array, resize);
1312 0 : }
1313 :
1314 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Float> &column, const Slicer &slicer, Array<Float> &array, Bool resize) const
1315 : {
1316 0 : column.getColumn(slicer, array, resize);
1317 0 : }
1318 :
1319 0 : void ROVisIteratorImpl::getCol(const ArrayColumn<Complex> &column, const Slicer &slicer, Array<Complex> &array, Bool resize) const
1320 : {
1321 0 : column.getColumn(slicer, array, resize);
1322 0 : }
1323 :
1324 :
1325 : Vector<RigidVector<Double,3> >&
1326 0 : ROVisIteratorImpl::uvw(Vector<RigidVector<Double,3> >& uvwvec) const
1327 : {
1328 0 : uvwvec.resize(curNumRow_p);
1329 0 : getCol(columns_p.uvw_p, cache_p.uvwMat_p,true);
1330 : // get a pointer to the raw storage for quick access
1331 : Bool deleteIt;
1332 0 : Double* pmat = cache_p.uvwMat_p.getStorage(deleteIt);
1333 0 : for (uInt row=0; row<curNumRow_p; row++, pmat+=3) uvwvec(row)=pmat;
1334 0 : return uvwvec;
1335 : }
1336 :
1337 :
1338 : const Table
1339 60 : ROVisIteratorImpl::attachTable() const
1340 : {
1341 60 : return selTable_p;
1342 : }
1343 :
1344 : // ******************
1345 : // ******************
1346 : // ** **
1347 : // ** VisIterator **
1348 : // ** **
1349 : // ******************
1350 : // ******************
1351 :
1352 :
1353 0 : VisIterator::VisIterator() {}
1354 :
1355 0 : VisIterator::VisIterator(MeasurementSet &MS,
1356 : const Block<Int>& sortColumns,
1357 0 : Double timeInterval)
1358 0 : : ROVisIterator (MS, sortColumns, timeInterval, Factory (this))
1359 0 : {}
1360 :
1361 0 : VisIterator::VisIterator(const VisIterator & other)
1362 0 : : ROVisIterator (other)
1363 : {
1364 0 : operator=(other);
1365 0 : }
1366 :
1367 0 : VisIterator::~VisIterator() {}
1368 :
1369 : VisIterator&
1370 0 : VisIterator::operator=(const VisIterator& other)
1371 : {
1372 0 : if (this!=&other) {
1373 :
1374 0 : ROVisIterator::operator= (other);
1375 : }
1376 :
1377 0 : return * this;
1378 : }
1379 :
1380 : VisIterator &
1381 0 : VisIterator::operator++(int)
1382 : {
1383 0 : advance();
1384 :
1385 0 : return *this;
1386 : }
1387 :
1388 : VisIterator &
1389 0 : VisIterator::operator++()
1390 : {
1391 0 : advance();
1392 :
1393 0 : return *this;
1394 : }
1395 :
1396 : void
1397 0 : VisIterator::attachColumns(const Table &table)
1398 : {
1399 0 : getImpl() -> attachColumns (table);
1400 0 : }
1401 :
1402 : VisIteratorImpl *
1403 0 : VisIterator::getImpl () const
1404 : {
1405 0 : VisIteratorImpl * impl = dynamic_cast<VisIteratorImpl *> (ROVisibilityIterator::getReadImpl ());
1406 :
1407 0 : Assert (impl != NULL);
1408 :
1409 0 : return impl;
1410 : }
1411 :
1412 0 : void VisIterator::setFlagRow(const Vector<Bool>& rowflags)
1413 : {
1414 0 : getImpl() -> setFlagRow (rowflags);
1415 0 : };
1416 :
1417 0 : void VisIterator::setVis(const Cube<Complex>& vis, DataColumn whichOne)
1418 : {
1419 0 : getImpl () -> setVis (vis, whichOne);
1420 0 : }
1421 :
1422 0 : void VisIterator::setFlag(const Cube<Bool>& flags)
1423 : {
1424 0 : getImpl () -> setFlag (flags);
1425 0 : }
1426 :
1427 0 : void VisIterator::setVisAndFlag(const Cube<Complex>& vis,
1428 : const Cube<Bool>& flag,
1429 : DataColumn whichOne)
1430 : {
1431 0 : getImpl () -> setVisAndFlag (vis, flag, whichOne);
1432 0 : }
1433 :
1434 0 : void VisIterator::setWeightMat(const Matrix<Float>& weightMat)
1435 : {
1436 0 : getImpl () -> setWeightMat (weightMat);
1437 0 : }
1438 :
1439 0 : void VisIterator::setWeightSpectrum(const Cube<Float>& weightSpectrum)
1440 : {
1441 0 : getImpl () -> setWeightSpectrum (weightSpectrum);
1442 0 : }
1443 :
1444 0 : void VisIterator::putDataColumn(DataColumn whichOne,
1445 : const Vector<Vector<Slice> >& slices,
1446 : const Cube<Complex>& data)
1447 : {
1448 0 : getImpl () -> putDataColumn (whichOne, slices, data);
1449 0 : };
1450 :
1451 0 : void VisIterator::putDataColumn(DataColumn whichOne,
1452 : const Cube<Complex>& data)
1453 : {
1454 0 : getImpl () -> putDataColumn (whichOne, data);
1455 0 : };
1456 :
1457 30 : Vector<rownr_t>& ROVisIteratorImpl::rowIds(Vector<rownr_t>& rowids) const
1458 : {
1459 30 : rowids.resize(curNumRow_p);
1460 30 : rowids=selTable_p.rowNumbers();
1461 30 : return rowids;
1462 : }
1463 :
1464 :
1465 0 : void VisIterator::putCol(ScalarColumn<Bool> &column, const Vector<Bool> &array)
1466 : {
1467 0 : column.putColumn(array);
1468 0 : }
1469 :
1470 0 : void VisIterator::putCol(ArrayColumn<Bool> &column, const Array<Bool> &array)
1471 : {
1472 0 : column.putColumn(array);
1473 0 : }
1474 :
1475 0 : void VisIterator::putCol(ArrayColumn<Float> &column, const Array<Float> &array)
1476 : {
1477 0 : column.putColumn(array);
1478 0 : }
1479 :
1480 0 : void VisIterator::putCol(ArrayColumn<Complex> &column, const Array<Complex> &array)
1481 : {
1482 0 : column.putColumn(array);
1483 0 : }
1484 :
1485 0 : void VisIterator::putCol(ArrayColumn<Bool> &column, const Slicer &slicer, const Array<Bool> &array)
1486 : {
1487 0 : column.putColumn(slicer, array);
1488 0 : }
1489 :
1490 0 : void VisIterator::putCol(ArrayColumn<Float> &column, const Slicer &slicer, const Array<Float> &array)
1491 : {
1492 0 : column.putColumn(slicer, array);
1493 0 : }
1494 :
1495 0 : void VisIterator::putCol(ArrayColumn<Complex> &column, const Slicer &slicer, const Array<Complex> &array)
1496 : {
1497 0 : column.putColumn(slicer, array);
1498 0 : }
1499 :
1500 : VisibilityIteratorReadImpl *
1501 0 : VisIterator::Factory::operator() (const asyncio::PrefetchColumns * /*prefetchColumns*/,
1502 : const Block<MeasurementSet>& mss,
1503 : const Block<Int>& sortColumns,
1504 : const Bool /*addDefaultSortCols*/,
1505 : Double timeInterval) const
1506 : {
1507 0 : return new VisIteratorImpl (vi_p, mss[0], sortColumns, timeInterval);
1508 : }
1509 :
1510 : // ***********************
1511 : // ***********************
1512 : // ** **
1513 : // ** VisIteratorImpl **
1514 : // ** **
1515 : // ***********************
1516 : // ***********************
1517 :
1518 0 : VisIteratorImpl::VisIteratorImpl(ROVisIterator * rovi,
1519 : const MeasurementSet& ms,
1520 : const Block<Int>& sortColumns,
1521 0 : Double timeInterval)
1522 0 : : ROVisIteratorImpl (rovi, ms, sortColumns, timeInterval)
1523 : {
1524 0 : }
1525 :
1526 0 : VisIteratorImpl::~VisIteratorImpl()
1527 0 : {}
1528 :
1529 : VisIteratorImpl &
1530 0 : VisIteratorImpl::operator=(const VisIteratorImpl &other){
1531 :
1532 0 : if (this != & other){
1533 :
1534 0 : ROVisIteratorImpl::operator=(other);
1535 :
1536 0 : rwColFlag_p.reference(other.rwColFlag_p);
1537 0 : rwColFlagRow_p.reference(other.rwColFlagRow_p);
1538 0 : rwColVis_p.reference(other.rwColVis_p);
1539 0 : rwColFloatVis_p.reference(other.rwColFloatVis_p);
1540 0 : rwColModelVis_p.reference(other.rwColModelVis_p);
1541 0 : rwColCorrVis_p.reference(other.rwColCorrVis_p);
1542 0 : rwColWeight_p.reference(other.rwColWeight_p);
1543 0 : rwColWeightSpectrum_p.reference(other.rwColWeightSpectrum_p);
1544 0 : rwColSigma_p.reference(other.rwColSigma_p);
1545 :
1546 : }
1547 :
1548 0 : return * this;
1549 : }
1550 :
1551 :
1552 : void
1553 0 : VisIteratorImpl::attachColumns(const Table &table)
1554 : {
1555 :
1556 0 : ROVisIteratorImpl::attachColumns(table);
1557 :
1558 0 : const ColumnDescSet& cds = table.tableDesc().columnDescSet();
1559 :
1560 0 : if (cds.isDefined("CORRECTED_DATA")){
1561 0 : rwColCorrVis_p.attach(table,"CORRECTED_DATA");
1562 : }
1563 :
1564 0 : if (cds.isDefined(MS::columnName(MS::DATA))) {
1565 0 : rwColVis_p.attach(table,MS::columnName(MS::DATA));
1566 : }
1567 :
1568 0 : rwColFlag_p.attach(table,MS::columnName(MS::FLAG));
1569 :
1570 0 : rwColFlagRow_p.attach(table,MS::columnName(MS::FLAG_ROW));
1571 :
1572 0 : if (cds.isDefined(MS::columnName(MS::FLOAT_DATA))) {
1573 0 : floatDataFound_p=true;
1574 0 : rwColFloatVis_p.attach(table,MS::columnName(MS::FLOAT_DATA));
1575 : }
1576 : else {
1577 0 : floatDataFound_p=false;
1578 : }
1579 :
1580 0 : if (cds.isDefined("MODEL_DATA")){
1581 0 : rwColModelVis_p.attach(table,"MODEL_DATA");
1582 : }
1583 :
1584 0 : rwColSigma_p.attach(table,MS::columnName(MS::SIGMA));
1585 :
1586 0 : rwColWeight_p.attach(table,MS::columnName(MS::WEIGHT));
1587 :
1588 0 : if (cds.isDefined("WEIGHT_SPECTRUM")){
1589 0 : rwColWeightSpectrum_p.attach(table,"WEIGHT_SPECTRUM");
1590 : }
1591 0 : }
1592 : void
1593 0 : VisIteratorImpl::putDataColumn(DataColumn whichOne,
1594 : const Vector<Vector<Slice> >& slices,
1595 : const Cube<Complex>& data)
1596 : {
1597 : // Set the visibility (observed, model or corrected);
1598 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
1599 :
1600 0 : switch (whichOne) {
1601 :
1602 0 : case VisibilityIterator::Observed:
1603 :
1604 0 : if (floatDataFound_p) {
1605 0 : Cube<Float> dataFloat=real(data);
1606 0 : rwColFloatVis_p.putColumn(slices,dataFloat);
1607 : } else {
1608 0 : rwColVis_p.putColumn(slices,data);
1609 : };
1610 0 : break;
1611 :
1612 0 : case VisibilityIterator::Corrected:
1613 :
1614 0 : rwColCorrVis_p.putColumn(slices,data);
1615 0 : break;
1616 :
1617 0 : case VisibilityIterator::Model:
1618 :
1619 0 : rwColModelVis_p.putColumn(slices,data);
1620 0 : break;
1621 : };
1622 0 : };
1623 :
1624 : void
1625 0 : VisIteratorImpl::putDataColumn(DataColumn whichOne,
1626 : const Cube<Complex>& data)
1627 : {
1628 : // Set the visibility (observed, model or corrected);
1629 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
1630 :
1631 0 : switch (whichOne) {
1632 :
1633 0 : case VisibilityIterator::Observed:
1634 :
1635 0 : if (floatDataFound_p) {
1636 0 : Cube<Float> dataFloat=real(data);
1637 0 : rwColFloatVis_p.putColumn(dataFloat);
1638 : } else {
1639 0 : rwColVis_p.putColumn(data);
1640 : };
1641 0 : break;
1642 :
1643 0 : case VisibilityIterator::Corrected:
1644 :
1645 0 : rwColCorrVis_p.putColumn(data);
1646 0 : break;
1647 :
1648 0 : case VisibilityIterator::Model:
1649 :
1650 0 : rwColModelVis_p.putColumn(data);
1651 0 : break;
1652 : };
1653 0 : };
1654 :
1655 :
1656 : void
1657 0 : VisIteratorImpl::setFlag(const Cube<Bool>& flags)
1658 : {
1659 0 : if (useNewSlicer_p){
1660 0 : rwColFlag_p.putColumn (newSlicer_p, flags);
1661 : }
1662 : else{
1663 0 : rwColFlag_p.putColumn(flags);
1664 : }
1665 0 : }
1666 :
1667 : void
1668 0 : VisIteratorImpl::setFlagRow(const Vector<Bool>& rowflags)
1669 : {
1670 0 : rwColFlagRow_p.putColumn(rowflags);
1671 0 : };
1672 :
1673 : void
1674 0 : VisIteratorImpl::setVis(const Cube<Complex>& vis, DataColumn whichOne)
1675 : {
1676 :
1677 0 : if (useNewSlicer_p){
1678 0 : putDataColumn (whichOne, newSlicer_p, vis);
1679 : }
1680 : else {
1681 0 : putDataColumn (whichOne, vis);
1682 : }
1683 :
1684 0 : }
1685 :
1686 : void
1687 0 : VisIteratorImpl::setVisAndFlag(const Cube<Complex>& vis,
1688 : const Cube<Bool>& flag,
1689 : DataColumn whichOne)
1690 : {
1691 0 : this->setFlag(flag);
1692 0 : this->setVis(vis,whichOne);
1693 0 : }
1694 :
1695 : void
1696 0 : VisIteratorImpl::setWeightMat(const Matrix<Float>& weightMat)
1697 : {
1698 0 : if (useNewSlicer_p){
1699 0 : rwColWeight_p.putColumn (newWtSlicer_p, weightMat);
1700 : }
1701 : else{
1702 0 : rwColWeight_p.putColumn (weightMat);
1703 : }
1704 0 : }
1705 :
1706 : void
1707 0 : VisIteratorImpl::setWeightSpectrum(const Cube<Float>& weightSpectrum)
1708 : {
1709 0 : if (existsWeightSpectrum()) {
1710 :
1711 0 : if (useNewSlicer_p){
1712 0 : rwColWeightSpectrum_p.putColumn(newSlicer_p,weightSpectrum);
1713 : }
1714 : else{
1715 0 : rwColWeightSpectrum_p.putColumn(weightSpectrum);
1716 : }
1717 : }
1718 : else {
1719 0 : throw(AipsError("Can't set WEIGHT_SPECTRUM -- it doesn't exist!"));
1720 : }
1721 0 : }
1722 :
1723 : } //# NAMESPACE CASA - END
1724 :
|
