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() const {return offset_p;}
64  virtual void setOffset(Double offset) {offset_p=offset;}
65  Double getInterval() const {return interval_p;}
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  virtual void origin();
220 
221  // Return False if there is no more data
222  virtual Bool more() const;
223 
224  // Advance iterator through data
225  virtual MSIter & operator++(int);
226  virtual MSIter & operator++();
227 
228  // Report Name of slowest column that changes at end of current iteration
229  const String& keyChange() const;
230 
231  // Return the current Table iteration
232  Table table() const;
233 
234  // Return reference to the current MS
235  const MS& ms() const;
236 
237  // Return reference to the current ROMSColumns
238  const ROMSColumns& msColumns() const;
239 
240  // Return the current MS Id (according to the order in which
241  // they appeared in the constructor)
242  Int msId() const;
243 
244  // Return true if msId has changed since last iteration
245  Bool newMS() const;
246 
247  // Return the current ArrayId
248  Int arrayId() const;
249 
250  // Return True if ArrayId has changed since last iteration
251  Bool newArray() const;
252 
253  // Return the current FieldId
254  Int fieldId() const;
255 
256  // Return the current Field Name
257  const String& fieldName() const;
258 
259  // Return the current Source Name
260  const String& sourceName() const;
261 
262  // Return True if FieldId/Source has changed since last iteration
263  Bool newField() const;
264 
265  // Return current SpectralWindow
266  Int spectralWindowId() const;
267 
268  // Return True if SpectralWindow has changed since last iteration
269  Bool newSpectralWindow() const;
270 
271  // Return current DataDescriptionId
272  Int dataDescriptionId() const;
273 
274  // Return True if DataDescriptionId has changed since last iteration
275  Bool newDataDescriptionId() const;
276 
277  // Return current PolarizationId
278  Int polarizationId() const;
279 
280  // Return True if polarization has changed since last iteration
281  Bool newPolarizationId() const;
282 
283  // Return the current phase center as MDirection
284  const MDirection& phaseCenter() const;
285 
286  // Return frame for polarization (returns PolFrame enum)
287  Int polFrame() const;
288 
289  // Return the frequencies corresponding to the DATA matrix.
290  const Vector<Double>& frequency() const;
291 
292  // Return frequency of first channel with reference frame as a Measure.
293  // The reference frame Epoch is that of the first row, reset it as needed
294  // for each row.
295  // The reference frame Position is the average of the antenna positions.
296  const MFrequency& frequency0() const;
297 
298  // Return the rest frequency of the specified line as a Measure
299  const MFrequency& restFrequency(Int line=0) const;
300 
301  // Return the telescope position (if a known telescope) or the
302  // position of the first antenna (if unknown)
303  const MPosition& telescopePosition() const;
304 
305  // Return the feed configuration/leakage matrix for feed 0 on each antenna
306  // TODO: CJonesAll can be used instead of this method in all instances
307  const Vector<SquareMatrix<Complex,2> >& CJones() const;
308 
309  // Return the feed configuration/leakage matrix for all feeds and antennae
310  // First axis is antennaId, 2nd axis is feedId. Result of CJones() is
311  // a reference to the first column of the matrix returned by this method
312  const Matrix<SquareMatrix<Complex,2> >& CJonesAll() const;
313 
314  // Return the receptor angle for feed 0 on each antenna.
315  // First axis is receptor number, 2nd axis is antennaId.
316  // TODO: receptorAngles() can be used instead of this method
317  const Matrix<Double>& receptorAngle() const;
318 
319  // Return the receptor angles for all feeds and antennae
320  // First axis is a receptor number, 2nd axis is antennaId,
321  // 3rd axis is feedId. Result of receptorAngle() is just a reference
322  // to the first plane of the cube returned by this method
323  const Cube<Double>& receptorAngles() const;
324 
325  // Return the channel number of the first channel in the DATA.
326  // (non-zero for reference MS created by VisSet with channel selection)
327  Int startChan() const;
328 
329  // Return a string mount identifier for each antenna
330  const Vector<String>& antennaMounts() const;
331 
332  // Return a cube containing pairs of coordinate offset for each receptor
333  // of each feed (values are in radians, coordinate system is fixed with
334  // antenna and is the same as used to define the BEAM_OFFSET parameter
335  // in the feed table). The cube axes are receptor, antenna, feed.
336  const Cube<RigidVector<Double, 2> >& getBeamOffsets() const;
337 
338  // True if all elements of the cube returned by getBeamOffsets are zero
339  Bool allBeamOffsetsZero() const;
340 
341  // Get the spw, start and nchan for all the ms's is this msiter that
342  // match the frequecy "freqstart-freqStep" and "freqEnd+freqStep" range
343 
344  void getSpwInFreqRange(Block<Vector<Int> >& spw,
345  Block<Vector<Int> >& start,
346  Block<Vector<Int> >& nchan,
347  Double freqStart, Double freqEnd,
348  Double freqStep);
349 
350  //Get the number of actual ms's associated wth this iterator
351  Int numMS() const;
352 
353  //Get a reference to the nth ms in the list of ms associated with this
354  // iterator. If larger than the list of ms's current ms is returned
355  // So better check wth numMS() before making the call
356  const MS& ms(const uInt n) const;
357 
358 protected:
359  // handle the construction details
360  void construct(const Block<Int>& sortColumns, Bool addDefaultSortColumns);
361  // advance the iteration
362  void advance();
363  // set the iteration state
364  virtual void setState();
365  void setMSInfo();
366  void setArrayInfo();
367  void setFeedInfo();
368  void setDataDescInfo();
369  void setFieldInfo();
370 
371 // Determine if the numbers in r1 are a sorted subset of those in r2
372  Bool isSubSet(const Vector<uInt>& r1, const Vector<uInt>& r2);
373 
374  MSIter* This;
375  Block<MeasurementSet> bms_p;
376  PtrBlock<TableIterator* > tabIter_p;
377  Block<Bool> tabIterAtStart_p;
378 
379  Int nMS_p;
380  ROMSColumns* msc_p;
381  Table curTable_p;
382  Int curMS_p, lastMS_p, curArray_p, lastArray_p, curSource_p;
383  String curFieldName_p, curSourceName_p;
384  Int curField_p, lastField_p, curSpectralWindow_p, lastSpectralWindow_p;
385  Int curPolarizationId_p, lastPolarizationId_p;
386  Int curDataDescId_p, lastDataDescId_p;
387  Bool more_p, newMS_p, newArray_p, newField_p, newSpectralWindow_p,
388  newPolarizationId_p, newDataDescId_p, preselected_p,
389  timeDepFeed_p, spwDepFeed_p, checkFeed_p;
390  Int startChan_p;
391 
392  // Globally control disk storage of SORTED_TABLE
393  Bool storeSorted_p;
394 
395  // time selection
397  // channel selection
398  Block<Int> preselectedChanStart_p,preselectednChan_p;
399 
400  // columns
401  ScalarColumn<Int> colArray_p, colDataDesc_p, colField_p;
402 
403  //cache for access functions
404  MDirection phaseCenter_p;
405  Matrix<Double> receptorAnglesFeed0_p; // former receptorAngle_p,
406  // temporary retained for compatibility
407  // contain actually a reference to the
408  // first plane of receptorAngles_p
409  Cube<Double> receptorAngles_p;
410  Vector<SquareMatrix<Complex,2> > CJonesFeed0_p; // a temporary reference
411  // similar to receptorAngle_p
413  Vector<String> antennaMounts_p; // a string mount identifier for each
414  // antenna (e.g. EQUATORIAL, ALT-AZ,...)
415  Cube<RigidVector<Double, 2> > beamOffsets_p;// angular offsets (two values for
416  // each element of the cube in radians)
417  // in the antenna coordinate system.
418  // Cube axes are: receptor, antenna, feed.
419  Bool allBeamOffsetsZero_p; // True if all elements of beamOffsets_p
420  // are zero (to speed things up in a
421  // single beam case)
422  PolFrame polFrame_p;
423  Bool freqCacheOK_p;
424  Vector<Double> frequency_p;
425  MFrequency frequency0_p;
426  MFrequency restFrequency_p;
427  MPosition telescopePosition_p;
428 
429  MSInterval *timeComp_p; // Points to the time comparator.
430  // 0 if not using a time interval.
431 };
433 inline Bool MSIter::more() const { return more_p;}
434 inline Table MSIter::table() const {return curTable_p;}
435 inline const MS& MSIter::ms() const {return bms_p[curMS_p];}
436 inline const ROMSColumns& MSIter::msColumns() const { return *msc_p;}
437 inline Bool MSIter::newMS() const { return newMS_p;}
438 inline Bool MSIter::newArray() const {return newArray_p;}
439 inline Bool MSIter::newField() const { return newField_p;}
441 { return newSpectralWindow_p;}
442 inline Int MSIter::msId() const { return curMS_p;}
443 inline Int MSIter::numMS() const { return nMS_p;}
444 inline Int MSIter::arrayId() const {return curArray_p;}
445 inline Int MSIter::fieldId() const { return curField_p;}
446 inline const String& MSIter::fieldName() const { return curFieldName_p;}
447 inline const String& MSIter::sourceName() const { return curSourceName_p;}
449 { return curSpectralWindow_p;}
450 inline Int MSIter::polarizationId() const {return curPolarizationId_p;}
451 inline Int MSIter::dataDescriptionId() const {return curDataDescId_p;}
452 inline Bool MSIter::newPolarizationId() const { return newPolarizationId_p;}
453 inline Bool MSIter::newDataDescriptionId() const { return newDataDescId_p;}
454 inline Int MSIter::polFrame() const { return polFrame_p;}
455 inline const MDirection& MSIter::phaseCenter() const
456 { return phaseCenter_p; }
457 inline const MPosition& MSIter::telescopePosition() const
458 { return telescopePosition_p;}
460 { return CJonesFeed0_p;}
462 { return CJones_p;}
463 inline const Matrix<Double>& MSIter::receptorAngle() const
464 {return receptorAnglesFeed0_p;}
466 {return receptorAngles_p;}
467 inline const Vector<String>& MSIter::antennaMounts() const
468 {return antennaMounts_p;}
470 {return beamOffsets_p;}
471 inline Int MSIter::startChan() const {return startChan_p;}
472 inline Bool MSIter::allBeamOffsetsZero() const {return allBeamOffsetsZero_p;}
474 } //# NAMESPACE CASACORE - END
475 
476 #endif
const MS & ms() const
Return reference to the current MS.
Definition: MSIter.h:496
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:506
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
virtual Bool more() const
Return False if there is no more data.
Definition: MSIter.h:494
Bool allBeamOffsetsZero() const
True if all elements of the cube returned by getBeamOffsets are zero.
Definition: MSIter.h:533
Int polFrame() const
Return frame for polarization (returns PolFrame enum)
Definition: MSIter.h:515
const MDirection & phaseCenter() const
Return the current phase center as MDirection.
Definition: MSIter.h:516
Int numMS() const
Get the number of actual ms&#39;s associated wth this iterator.
Definition: MSIter.h:504
PtrHolder< T > & operator=(const PtrHolder< T > &other)
Int polarizationId() const
Return current PolarizationId.
Definition: MSIter.h:511
virtual void setOffset(Double offset)
Definition: MSIter.h:65
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:503
Bool newPolarizationId() const
Return True if polarization has changed since last iteration.
Definition: MSIter.h:513
Table table() const
Return the current Table iteration.
Definition: MSIter.h:495
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:524
Int startChan() const
Return the channel number of the first channel in the DATA.
Definition: MSIter.h:532
Bool newMS() const
Return true if msId has changed since last iteration.
Definition: MSIter.h:498
Bool newDataDescriptionId() const
Return True if DataDescriptionId has changed since last iteration.
Definition: MSIter.h:514
Double getInterval() const
Definition: MSIter.h:66
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:520
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:522
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:501
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:530
const Vector< String > & antennaMounts() const
Return a string mount identifier for each antenna.
Definition: MSIter.h:528
Bool newArray() const
Return True if ArrayId has changed since last iteration.
Definition: MSIter.h:499
A drop-in replacement for Block<T*>.
Definition: Block.h:861
Int arrayId() const
Return the current ArrayId.
Definition: MSIter.h:505
const String & fieldName() const
Return the current Field Name.
Definition: MSIter.h:507
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:500
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:526
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:508
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:512
const ROMSColumns & msColumns() const
Return reference to the current ROMSColumns.
Definition: MSIter.h:497
String: the storage and methods of handling collections of characters.
Definition: String.h:223
const MPosition & telescopePosition() const
Return the telescope position (if a known telescope) or the position of the first antenna (if unknown...
Definition: MSIter.h:518
Int spectralWindowId() const
Return current SpectralWindow.
Definition: MSIter.h:509
Double getOffset() const
Definition: MSIter.h:64
const Bool True
Definition: aipstype.h:43
this file contains all the compiler specific defines
Definition: mainpage.dox:28
unsigned int uInt
Definition: aipstype.h:51