casacore
MSIter.h
Go to the documentation of this file.
1 //# MSIter.h: Step through the MeasurementEquation by table
2 //# Copyright (C) 1996,1999,2000,2001,2002
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$
27 
28 #ifndef MS_MSITER_H
29 #define MS_MSITER_H
30 
31 #include <casacore/casa/aips.h>
32 #include <casacore/casa/Arrays/Matrix.h>
33 #include <casacore/casa/Arrays/Cube.h>
34 #include <casacore/ms/MeasurementSets/MeasurementSet.h>
35 #include <casacore/measures/Measures/MFrequency.h>
36 #include <casacore/measures/Measures/MDirection.h>
37 #include <casacore/measures/Measures/MPosition.h>
38 #include <casacore/tables/Tables/ScalarColumn.h>
39 #include <casacore/casa/Utilities/Compare.h>
40 #include <casacore/casa/BasicSL/String.h>
41 #include <casacore/scimath/Mathematics/SquareMatrix.h>
42 #include <casacore/scimath/Mathematics/RigidVector.h>
43 
44 namespace casacore { //# NAMESPACE CASACORE - BEGIN
45 
46 //# forward decl
47 class ROMSColumns;
48 class TableIterator;
49 
50 // <summary>
51 // Small helper class to specify an 'interval' comparison
52 // </summary>
53 // <synopsis>
54 // Small helper class to specify an 'interval' comparison for table iteration
55 // by time interval.
56 // </synopsis>
57 class MSInterval : public BaseCompare
58 {
59 public:
60  explicit MSInterval(Double interval) : interval_p(interval), offset_p(0) {}
61  virtual ~MSInterval() {}
62  virtual int comp(const void * obj1, const void * obj2) const;
63  Double getOffset() {return offset_p;}
64  void setOffset(Double offset) {offset_p=offset;}
66  void setInterval(Double interval) {interval_p=interval;}
67 private:
69  mutable Double offset_p;
70 };
71 
72 // <summary>
73 // An iterator class for MeasurementSets
74 // </summary>
75 
76 // <use visibility=export>
77 
78 // <prerequisite>
79 // <li> <linkto class="MeasurementSet:description">MeasurementSet</linkto>
80 // </prerequisite>
81 //
82 // <etymology>
83 // MSIter stands for the MeasurementSet Iterator class.
84 // </etymology>
85 //
86 // <synopsis>
87 // An MSIter is a class to traverse a MeasurementSet in various orders. It
88 // automatically adds four predefined sort columns to your selection of sort
89 // columns (see constructor) so that it can keep track of changes in frequency
90 // or polarization setup, field position and sub-array. Note that this can
91 // cause iterations to occur in a different way from what you would expect, see
92 // examples below. MSIter implements iteration by time interval for the use of
93 // e.g., calibration tasks that want to calculate solutions over some interval
94 // of time. You can iterate over multiple MeasurementSets with this class.
95 // </synopsis>
96 //
97 // <example>
98 // <srcblock>
99 // // The following code iterates by by ARRAY_ID, FIELD_ID, DATA_DESC_ID and
100 // // TIME (all implicitly added columns) and then by baseline (antenna pair),
101 // // in 3000s intervals.
102 // MeasurementSet ms("3C273XC1.ms");
103 // Block<int> sort(2);
104 // sort[0] = MS::ANTENNA1;
105 // sort[1] = MS::ANTENNA2;
106 // Double timeInterval = 3000;
107 // MSIter msIter(ms,sort,timeInteval);
108 // for (msIter.origin(); msIter.more(); msIter++) {
109 // // print out some of the iteration state
110 // cout << msIter.fieldId() << endl;
111 // cout << msIter.fieldName() << endl;
112 // cout << msIter.dataDescriptionId() << endl;
113 // cout << msIter.frequency0() << endl;
114 // cout << msIter.table().nrow() << endl;
115 // process(msIter.table()); // process the data in the current iteration
116 // }
117 // // Output shows only 1 row at a time because the table is sorted on TIME
118 // // first and ANTENNA1, ANTENNA2 next and each baseline occurs only once per
119 // // TIME stamp. The interval has no effect in this case.
120 // </srcblock>
121 // </example>
122 
123 // <example>
124 // <srcblock>
125 // // The following code iterates by baseline (antenna pair), TIME, and,
126 // // implicitly, by ARRAY_ID, FIELD_ID and DATA_DESC_ID in 3000s
127 // // intervals.
128 // MeasurementSet ms("3C273XC1.ms");
129 // Block<int> sort(3);
130 // sort[0] = MS::ANTENNA1;
131 // sort[1] = MS::ANTENNA2;
132 // sort[2] = MS::TIME;
133 // Double timeInterval = 3000;
134 // MSIter msIter(ms,sort,timeInteval);
135 // for (msIter.origin(); msIter.more(); msIter++) {
136 // // print out some of the iteration state
137 // cout << msIter.fieldId() << endl;
138 // cout << msIter.fieldName() << endl;
139 // cout << msIter.dataDescriptionId() << endl;
140 // cout << msIter.frequency0() << endl;
141 // cout << msIter.table().nrow() << endl;
142 // process(msIter.table()); // process the data in the current iteration
143 // // Now the output shows 7 rows at a time, all with identical ANTENNA1
144 // // and ANTENNA2 values and TIME values within a 3000s interval.
145 // }
146 // </srcblock>
147 // </example>
148 //
149 // <motivation>
150 // This class was originally part of the VisibilityIterator class, but that
151 // class was getting too large and complicated. By splitting out the toplevel
152 // iteration into this class the code is much easier to understand. It is now
153 // also available through the ms tool.
154 // </motivation>
155 //
156 // <todo>
157 // multiple observatories in a single MS are not handled correctly (need to
158 // sort on observation id and check observatory name to set position frame)
159 // </todo>
160 
161 class MSIter
162 {
163 public:
164  enum PolFrame {
165  // Circular polarization
166  Circular=0,
167  // Linear polarization
168  Linear=1
169  };
170 
171  // Default constructor - useful only to assign another iterator later.
172  // Use of other member functions on this object is likely to dump core.
173  MSIter();
174 
175  // Construct from MS and a Block of MS column enums specifying the
176  // iteration order, if none are specified, ARRAY_ID, FIELD_ID, DATA_DESC_ID,
177  // and TIME iteration is implicit (unless addDefaultSortColumns=False)
178  // These columns will be added first if they are not specified.
179  // An optional timeInterval can be given to iterate through chunks of time.
180  // The default interval of 0 groups all times together.
181  // Every 'chunk' of data contains all data within a certain time interval
182  // and with identical values of the other iteration columns (e.g.
183  // DATA_DESCRIPTION_ID and FIELD_ID).
184  // See the examples above for the effect of different sort orders.
185  //
186  // The storeSorted parameter determines how the resulting SORT_TABLE is
187  // managed. If storeSorted is true then the table will be stored on disk;
188  // this potentially allows its reuse in the future but has also been shown
189  // to be a problem when the MS is being read in parallel. If storeSorted is
190  // false then the SORTED_TABLE is constructed and used in memory which keeps
191  // concurrent readers from interfering with each other.
192 
193  MSIter(const MeasurementSet& ms, const Block<Int>& sortColumns,
194  Double timeInterval=0, Bool addDefaultSortColumns=True,
195  Bool storeSorted=True);
196 
197  // Same as above with multiple MSs as input.
198  MSIter(const Block<MeasurementSet>& mss, const Block<Int>& sortColumns,
199  Double timeInterval=0, Bool addDefaultSortColumns=True,
200  Bool storeSorted=True);
201 
202  // Copy construct. This calls the assigment operator.
203  MSIter(const MSIter & other);
204 
205  // Destructor
206  virtual ~MSIter();
207 
208  // Assigment. This will reset the iterator to the origin.
209  MSIter & operator=(const MSIter &other);
210 
211  //# Members
212 
213  // Set or reset the time interval to use for iteration.
214  // You should call origin() to reset the iteration after
215  // calling this.
216  void setInterval(Double timeInterval);
217 
218  // Reset iterator to start of data
219  void origin();
220 
221  // Return False if there is no more data
222  Bool more() const;
223 
224  // Advance iterator through data
225  MSIter & operator++(int);
226  MSIter & operator++();
227 
228  // Return the current Table iteration
229  Table table() const;
230 
231  // Return reference to the current MS
232  const MS& ms() const;
233 
234  // Return reference to the current ROMSColumns
235  const ROMSColumns& msColumns() const;
236 
237  // Return the current MS Id (according to the order in which
238  // they appeared in the constructor)
239  Int msId() const;
240 
241  // Return true if msId has changed since last iteration
242  Bool newMS() const;
243 
244  // Return the current ArrayId
245  Int arrayId() const;
246 
247  // Return True if ArrayId has changed since last iteration
248  Bool newArray() const;
249 
250  // Return the current FieldId
251  Int fieldId() const;
252 
253  // Return the current Field Name
254  const String& fieldName() const;
255 
256  // Return the current Source Name
257  const String& sourceName() const;
258 
259  // Return True if FieldId/Source has changed since last iteration
260  Bool newField() const;
261 
262  // Return current SpectralWindow
263  Int spectralWindowId() const;
264 
265  // Return True if SpectralWindow has changed since last iteration
266  Bool newSpectralWindow() const;
267 
268  // Return current DataDescriptionId
269  Int dataDescriptionId() const;
270 
271  // Return True if DataDescriptionId has changed since last iteration
272  Bool newDataDescriptionId() const;
273 
274  // Return current PolarizationId
275  Int polarizationId() const;
276 
277  // Return True if polarization has changed since last iteration
278  Bool newPolarizationId() const;
279 
280  // Return the current phase center as MDirection
281  const MDirection& phaseCenter() const;
282 
283  // Return frame for polarization (returns PolFrame enum)
284  Int polFrame() const;
285 
286  // Return the frequencies corresponding to the DATA matrix.
287  const Vector<Double>& frequency() const;
288 
289  // Return frequency of first channel with reference frame as a Measure.
290  // The reference frame Epoch is that of the first row, reset it as needed
291  // for each row.
292  // The reference frame Position is the average of the antenna positions.
293  const MFrequency& frequency0() const;
294 
295  // Return the rest frequency of the specified line as a Measure
296  const MFrequency& restFrequency(Int line=0) const;
297 
298  // Return the telescope position (if a known telescope) or the
299  // position of the first antenna (if unknown)
300  const MPosition& telescopePosition() const;
301 
302  // Return the feed configuration/leakage matrix for feed 0 on each antenna
303  // TODO: CJonesAll can be used instead of this method in all instances
304  const Vector<SquareMatrix<Complex,2> >& CJones() const;
305 
306  // Return the feed configuration/leakage matrix for all feeds and antennae
307  // First axis is antennaId, 2nd axis is feedId. Result of CJones() is
308  // a reference to the first column of the matrix returned by this method
309  const Matrix<SquareMatrix<Complex,2> >& CJonesAll() const;
310 
311  // Return the receptor angle for feed 0 on each antenna.
312  // First axis is receptor number, 2nd axis is antennaId.
313  // TODO: receptorAngles() can be used instead of this method
314  const Matrix<Double>& receptorAngle() const;
315 
316  // Return the receptor angles for all feeds and antennae
317  // First axis is a receptor number, 2nd axis is antennaId,
318  // 3rd axis is feedId. Result of receptorAngle() is just a reference
319  // to the first plane of the cube returned by this method
320  const Cube<Double>& receptorAngles() const;
321 
322  // Return the channel number of the first channel in the DATA.
323  // (non-zero for reference MS created by VisSet with channel selection)
324  Int startChan() const;
325 
326  // Return a string mount identifier for each antenna
327  const Vector<String>& antennaMounts() const;
328 
329  // Return a cube containing pairs of coordinate offset for each receptor
330  // of each feed (values are in radians, coordinate system is fixed with
331  // antenna and is the same as used to define the BEAM_OFFSET parameter
332  // in the feed table). The cube axes are receptor, antenna, feed.
333  const Cube<RigidVector<Double, 2> >& getBeamOffsets() const;
334 
335  // True if all elements of the cube returned by getBeamOffsets are zero
336  Bool allBeamOffsetsZero() const;
337 
338  // Get the spw, start and nchan for all the ms's is this msiter that
339  // match the frequecy "freqstart-freqStep" and "freqEnd+freqStep" range
340 
341  void getSpwInFreqRange(Block<Vector<Int> >& spw,
342  Block<Vector<Int> >& start,
343  Block<Vector<Int> >& nchan,
344  Double freqStart, Double freqEnd,
345  Double freqStep);
346 
347  //Get the number of actual ms's associated wth this iterator
348  Int numMS() const;
349 
350  //Get a reference to the nth ms in the list of ms associated with this
351  // iterator. If larger than the list of ms's current ms is returned
352  // So better check wth numMS() before making the call
353  const MS& ms(const uInt n) const;
354 
355 protected:
356  // handle the construction details
357  void construct(const Block<Int>& sortColumns, Bool addDefaultSortColumns);
358  // advance the iteration
359  void advance();
360  // set the iteration state
361  void setState();
362  void setMSInfo();
363  void setArrayInfo();
364  void setFeedInfo();
365  void setDataDescInfo();
366  void setFieldInfo();
367 
368 // Determine if the numbers in r1 are a sorted subset of those in r2
369  Bool isSubSet(const Vector<uInt>& r1, const Vector<uInt>& r2);
370 
371  MSIter* This;
372  Block<MeasurementSet> bms_p;
373  PtrBlock<TableIterator* > tabIter_p;
374  Block<Bool> tabIterAtStart_p;
375 
376  Int nMS_p;
377  ROMSColumns* msc_p;
378  Table curTable_p;
379  Int curMS_p, lastMS_p, curArray_p, lastArray_p, curSource_p;
380  String curFieldName_p, curSourceName_p;
381  Int curField_p, lastField_p, curSpectralWindow_p, lastSpectralWindow_p;
382  Int curPolarizationId_p, lastPolarizationId_p;
383  Int curDataDescId_p, lastDataDescId_p;
384  Bool more_p, newMS_p, newArray_p, newField_p, newSpectralWindow_p,
385  newPolarizationId_p, newDataDescId_p, preselected_p,
386  timeDepFeed_p, spwDepFeed_p, checkFeed_p;
387  Int startChan_p;
388 
389  // Globally control disk storage of SORTED_TABLE
390  Bool storeSorted_p;
391 
392  // time selection
394  // channel selection
395  Block<Int> preselectedChanStart_p,preselectednChan_p;
396 
397  // columns
398  ScalarColumn<Int> colArray_p, colDataDesc_p, colField_p;
399 
400  //cache for access functions
401  MDirection phaseCenter_p;
402  Matrix<Double> receptorAnglesFeed0_p; // former receptorAngle_p,
403  // temporary retained for compatibility
404  // contain actually a reference to the
405  // first plane of receptorAngles_p
406  Cube<Double> receptorAngles_p;
407  Vector<SquareMatrix<Complex,2> > CJonesFeed0_p; // a temporary reference
408  // similar to receptorAngle_p
410  Vector<String> antennaMounts_p; // a string mount identifier for each
411  // antenna (e.g. EQUATORIAL, ALT-AZ,...)
412  Cube<RigidVector<Double, 2> > beamOffsets_p;// angular offsets (two values for
413  // each element of the cube in radians)
414  // in the antenna coordinate system.
415  // Cube axes are: receptor, antenna, feed.
416  Bool allBeamOffsetsZero_p; // True if all elements of beamOffsets_p
417  // are zero (to speed things up in a
418  // single beam case)
419  PolFrame polFrame_p;
420  Bool freqCacheOK_p;
421  Vector<Double> frequency_p;
422  MFrequency frequency0_p;
423  MFrequency restFrequency_p;
424  MPosition telescopePosition_p;
426  MSInterval *timeComp_p; // Points to the time comparator.
427  // 0 if not using a time interval.
428 };
430 inline Bool MSIter::more() const { return more_p;}
431 inline Table MSIter::table() const {return curTable_p;}
432 inline const MS& MSIter::ms() const {return bms_p[curMS_p];}
433 inline const ROMSColumns& MSIter::msColumns() const { return *msc_p;}
434 inline Bool MSIter::newMS() const { return newMS_p;}
435 inline Bool MSIter::newArray() const {return newArray_p;}
436 inline Bool MSIter::newField() const { return newField_p;}
437 inline Bool MSIter::newSpectralWindow() const
438 { return newSpectralWindow_p;}
439 inline Int MSIter::msId() const { return curMS_p;}
440 inline Int MSIter::numMS() const { return nMS_p;}
441 inline Int MSIter::arrayId() const {return curArray_p;}
442 inline Int MSIter::fieldId() const { return curField_p;}
443 inline const String& MSIter::fieldName() const { return curFieldName_p;}
444 inline const String& MSIter::sourceName() const { return curSourceName_p;}
445 inline Int MSIter::spectralWindowId() const
446 { return curSpectralWindow_p;}
447 inline Int MSIter::polarizationId() const {return curPolarizationId_p;}
448 inline Int MSIter::dataDescriptionId() const {return curDataDescId_p;}
449 inline Bool MSIter::newPolarizationId() const { return newPolarizationId_p;}
450 inline Bool MSIter::newDataDescriptionId() const { return newDataDescId_p;}
451 inline Int MSIter::polFrame() const { return polFrame_p;}
452 inline const MDirection& MSIter::phaseCenter() const
453 { return phaseCenter_p; }
454 inline const MPosition& MSIter::telescopePosition() const
455 { return telescopePosition_p;}
457 { return CJonesFeed0_p;}
459 { return CJones_p;}
460 inline const Matrix<Double>& MSIter::receptorAngle() const
461 {return receptorAnglesFeed0_p;}
463 {return receptorAngles_p;}
464 inline const Vector<String>& MSIter::antennaMounts() const
465 {return antennaMounts_p;}
467 {return beamOffsets_p;}
468 inline Int MSIter::startChan() const {return startChan_p;}
469 inline Bool MSIter::allBeamOffsetsZero() const {return allBeamOffsetsZero_p;}
470 
471 } //# NAMESPACE CASACORE - END
472 
473 #endif
const MS & ms() const
Return reference to the current MS.
Definition: MSIter.h:492
A Measure: astronomical direction.
Definition: MDirection.h:174
A Measure: position on Earth.
Definition: MPosition.h:79
int Int
Definition: aipstype.h:50
Int fieldId() const
Return the current FieldId.
Definition: MSIter.h:502
Main interface class to a read/write table.
Definition: Table.h:149
void setInterval(Double interval)
Definition: MSIter.h:67
MSInterval(Double interval)
Definition: MSIter.h:61
Bool more() const
Return False if there is no more data.
Definition: MSIter.h:490
Bool allBeamOffsetsZero() const
True if all elements of the cube returned by getBeamOffsets are zero.
Definition: MSIter.h:529
Int polFrame() const
Return frame for polarization (returns PolFrame enum)
Definition: MSIter.h:511
const MDirection & phaseCenter() const
Return the current phase center as MDirection.
Definition: MSIter.h:512
Int numMS() const
Get the number of actual ms&#39;s associated wth this iterator.
Definition: MSIter.h:500
PtrHolder< T > & operator=(const PtrHolder< T > &other)
void setOffset(Double offset)
Definition: MSIter.h:65
Int polarizationId() const
Return current PolarizationId.
Definition: MSIter.h:507
A 2-D Specialization of the Array class.
Definition: Array.h:53
Int msId() const
Return the current MS Id (according to the order in which they appeared in the constructor) ...
Definition: MSIter.h:499
Bool newPolarizationId() const
Return True if polarization has changed since last iteration.
Definition: MSIter.h:509
Table table() const
Return the current Table iteration.
Definition: MSIter.h:491
A Measure: wave characteristics.
Definition: MFrequency.h:161
const Matrix< Double > & receptorAngle() const
Return the receptor angle for feed 0 on each antenna.
Definition: MSIter.h:520
Int startChan() const
Return the channel number of the first channel in the DATA.
Definition: MSIter.h:528
Bool newMS() const
Return true if msId has changed since last iteration.
Definition: MSIter.h:494
Bool newDataDescriptionId() const
Return True if DataDescriptionId has changed since last iteration.
Definition: MSIter.h:510
Double interval_p
Definition: MSIter.h:69
double Double
Definition: aipstype.h:55
A class to provide easy read-only access to MeasurementSet columns.
Definition: MSColumns.h:111
const Vector< SquareMatrix< Complex, 2 > > & CJones() const
Return the feed configuration/leakage matrix for feed 0 on each antenna TODO: CJonesAll can be used i...
Definition: MSIter.h:516
const Matrix< SquareMatrix< Complex, 2 > > & CJonesAll() const
Return the feed configuration/leakage matrix for all feeds and antennae First axis is antennaId...
Definition: MSIter.h:518
bool Bool
Define the standard types used by Casacore.
Definition: aipstype.h:42
Bool newSpectralWindow() const
Return True if SpectralWindow has changed since last iteration.
Definition: MSIter.h:497
const Cube< RigidVector< Double, 2 > > & getBeamOffsets() const
Return a cube containing pairs of coordinate offset for each receptor of each feed (values are in rad...
Definition: MSIter.h:526
const Vector< String > & antennaMounts() const
Return a string mount identifier for each antenna.
Definition: MSIter.h:524
Bool newArray() const
Return True if ArrayId has changed since last iteration.
Definition: MSIter.h:495
A drop-in replacement for Block<T*>.
Definition: Block.h:861
Int arrayId() const
Return the current ArrayId.
Definition: MSIter.h:501
const String & fieldName() const
Return the current Field Name.
Definition: MSIter.h:503
A Table intended to hold astronomical data (a set of Measurements).
Bool newField() const
Return True if FieldId/Source has changed since last iteration.
Definition: MSIter.h:496
const Cube< Double > & receptorAngles() const
Return the receptor angles for all feeds and antennae First axis is a receptor number, 2nd axis is antennaId, 3rd axis is feedId.
Definition: MSIter.h:522
virtual ~MSInterval()
Definition: MSIter.h:62
Small helper class to specify an &#39;interval&#39; comparison.
Definition: MSIter.h:58
const String & sourceName() const
Return the current Source Name.
Definition: MSIter.h:504
virtual int comp(const void *obj1, const void *obj2) const
Compare two objects, and return.
An iterator class for MeasurementSets.
Definition: MSIter.h:162
Int dataDescriptionId() const
Return current DataDescriptionId.
Definition: MSIter.h:508
const ROMSColumns & msColumns() const
Return reference to the current ROMSColumns.
Definition: MSIter.h:493
String: the storage and methods of handling collections of characters.
Definition: String.h:223
Double getOffset()
Definition: MSIter.h:64
const MPosition & telescopePosition() const
Return the telescope position (if a known telescope) or the position of the first antenna (if unknown...
Definition: MSIter.h:514
Int spectralWindowId() const
Return current SpectralWindow.
Definition: MSIter.h:505
const Bool True
Definition: aipstype.h:43
this file contains all the compiler specific defines
Definition: mainpage.dox:28
Double getInterval()
Definition: MSIter.h:66
unsigned int uInt
Definition: aipstype.h:51