00001 //# MSIter.h: Step through the MeasurementEquation by table 00002 //# Copyright (C) 1996,1999,2000,2001,2002 00003 //# Associated Universities, Inc. Washington DC, USA. 00004 //# 00005 //# This library is free software; you can redistribute it and/or modify it 00006 //# under the terms of the GNU Library General Public License as published by 00007 //# the Free Software Foundation; either version 2 of the License, or (at your 00008 //# option) any later version. 00009 //# 00010 //# This library is distributed in the hope that it will be useful, but WITHOUT 00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 00013 //# License for more details. 00014 //# 00015 //# You should have received a copy of the GNU Library General Public License 00016 //# along with this library; if not, write to the Free Software Foundation, 00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 00018 //# 00019 //# Correspondence concerning AIPS++ should be addressed as follows: 00020 //# Internet email: aips2-request@nrao.edu. 00021 //# Postal address: AIPS++ Project Office 00022 //# National Radio Astronomy Observatory 00023 //# 520 Edgemont Road 00024 //# Charlottesville, VA 22903-2475 USA 00025 //# 00026 //# $Id$ 00027 00028 #ifndef MS_MSITER_H 00029 #define MS_MSITER_H 00030 00031 #include <casacore/casa/aips.h> 00032 #include <casacore/casa/Arrays/Matrix.h> 00033 #include <casacore/casa/Arrays/Cube.h> 00034 #include <casacore/ms/MeasurementSets/MeasurementSet.h> 00035 #include <casacore/measures/Measures/MFrequency.h> 00036 #include <casacore/measures/Measures/MDirection.h> 00037 #include <casacore/measures/Measures/MPosition.h> 00038 #include <casacore/tables/Tables/ScalarColumn.h> 00039 #include <casacore/casa/Utilities/Compare.h> 00040 #include <casacore/casa/BasicSL/String.h> 00041 #include <casacore/scimath/Mathematics/SquareMatrix.h> 00042 #include <casacore/scimath/Mathematics/RigidVector.h> 00043 00044 namespace casacore { //# NAMESPACE CASACORE - BEGIN 00045 00046 //# forward decl 00047 class ROMSColumns; 00048 class TableIterator; 00049 00050 // <summary> 00051 // Small helper class to specify an 'interval' comparison 00052 // </summary> 00053 // <synopsis> 00054 // Small helper class to specify an 'interval' comparison for table iteration 00055 // by time interval. 00056 // </synopsis> 00057 class MSInterval : public BaseCompare 00058 { 00059 public: 00060 explicit MSInterval(Double interval) : interval_p(interval), offset_p(0) {} 00061 virtual ~MSInterval() {} 00062 virtual int comp(const void * obj1, const void * obj2) const; 00063 Double getOffset() {return offset_p;} 00064 void setOffset(Double offset) {offset_p=offset;} 00065 Double getInterval() {return interval_p;} 00066 void setInterval(Double interval) {interval_p=interval;} 00067 private: 00068 Double interval_p; 00069 mutable Double offset_p; 00070 }; 00071 00072 // <summary> 00073 // An iterator class for MeasurementSets 00074 // </summary> 00075 00076 // <use visibility=export> 00077 00078 // <prerequisite> 00079 // <li> <linkto class="MeasurementSet:description">MeasurementSet</linkto> 00080 // </prerequisite> 00081 // 00082 // <etymology> 00083 // MSIter stands for the MeasurementSet Iterator class. 00084 // </etymology> 00085 // 00086 // <synopsis> 00087 // An MSIter is a class to traverse a MeasurementSet in various orders. It 00088 // automatically adds four predefined sort columns to your selection of sort 00089 // columns (see constructor) so that it can keep track of changes in frequency 00090 // or polarization setup, field position and sub-array. Note that this can 00091 // cause iterations to occur in a different way from what you would expect, see 00092 // examples below. MSIter implements iteration by time interval for the use of 00093 // e.g., calibration tasks that want to calculate solutions over some interval 00094 // of time. You can iterate over multiple MeasurementSets with this class. 00095 // </synopsis> 00096 // 00097 // <example> 00098 // <srcblock> 00099 // // The following code iterates by by ARRAY_ID, FIELD_ID, DATA_DESC_ID and 00100 // // TIME (all implicitly added columns) and then by baseline (antenna pair), 00101 // // in 3000s intervals. 00102 // MeasurementSet ms("3C273XC1.ms"); 00103 // Block<int> sort(2); 00104 // sort[0] = MS::ANTENNA1; 00105 // sort[1] = MS::ANTENNA2; 00106 // Double timeInterval = 3000; 00107 // MSIter msIter(ms,sort,timeInteval); 00108 // for (msIter.origin(); msIter.more(); msIter++) { 00109 // // print out some of the iteration state 00110 // cout << msIter.fieldId() << endl; 00111 // cout << msIter.fieldName() << endl; 00112 // cout << msIter.dataDescriptionId() << endl; 00113 // cout << msIter.frequency0() << endl; 00114 // cout << msIter.table().nrow() << endl; 00115 // process(msIter.table()); // process the data in the current iteration 00116 // } 00117 // // Output shows only 1 row at a time because the table is sorted on TIME 00118 // // first and ANTENNA1, ANTENNA2 next and each baseline occurs only once per 00119 // // TIME stamp. The interval has no effect in this case. 00120 // </srcblock> 00121 // </example> 00122 00123 // <example> 00124 // <srcblock> 00125 // // The following code iterates by baseline (antenna pair), TIME, and, 00126 // // implicitly, by ARRAY_ID, FIELD_ID and DATA_DESC_ID in 3000s 00127 // // intervals. 00128 // MeasurementSet ms("3C273XC1.ms"); 00129 // Block<int> sort(3); 00130 // sort[0] = MS::ANTENNA1; 00131 // sort[1] = MS::ANTENNA2; 00132 // sort[2] = MS::TIME; 00133 // Double timeInterval = 3000; 00134 // MSIter msIter(ms,sort,timeInteval); 00135 // for (msIter.origin(); msIter.more(); msIter++) { 00136 // // print out some of the iteration state 00137 // cout << msIter.fieldId() << endl; 00138 // cout << msIter.fieldName() << endl; 00139 // cout << msIter.dataDescriptionId() << endl; 00140 // cout << msIter.frequency0() << endl; 00141 // cout << msIter.table().nrow() << endl; 00142 // process(msIter.table()); // process the data in the current iteration 00143 // // Now the output shows 7 rows at a time, all with identical ANTENNA1 00144 // // and ANTENNA2 values and TIME values within a 3000s interval. 00145 // } 00146 // </srcblock> 00147 // </example> 00148 // 00149 // <motivation> 00150 // This class was originally part of the VisibilityIterator class, but that 00151 // class was getting too large and complicated. By splitting out the toplevel 00152 // iteration into this class the code is much easier to understand. It is now 00153 // also available through the ms tool. 00154 // </motivation> 00155 // 00156 // <todo> 00157 // multiple observatories in a single MS are not handled correctly (need to 00158 // sort on observation id and check observatory name to set position frame) 00159 // </todo> 00160 00161 class MSIter 00162 { 00163 public: 00164 enum PolFrame { 00165 // Circular polarization 00166 Circular=0, 00167 // Linear polarization 00168 Linear=1 00169 }; 00170 00171 // Default constructor - useful only to assign another iterator later. 00172 // Use of other member functions on this object is likely to dump core. 00173 MSIter(); 00174 00175 // Construct from MS and a Block of MS column enums specifying the 00176 // iteration order, if none are specified, ARRAY_ID, FIELD_ID, DATA_DESC_ID, 00177 // and TIME iteration is implicit (unless addDefaultSortColumns=False) 00178 // These columns will be added first if they are not specified. 00179 // An optional timeInterval can be given to iterate through chunks of time. 00180 // The default interval of 0 groups all times together. 00181 // Every 'chunk' of data contains all data within a certain time interval 00182 // and with identical values of the other iteration columns (e.g. 00183 // DATA_DESCRIPTION_ID and FIELD_ID). 00184 // See the examples above for the effect of different sort orders. 00185 // 00186 // The storeSorted parameter determines how the resulting SORT_TABLE is 00187 // managed. If storeSorted is true then the table will be stored on disk; 00188 // this potentially allows its reuse in the future but has also been shown 00189 // to be a problem when the MS is being read in parallel. If storeSorted is 00190 // false then the SORTED_TABLE is constructed and used in memory which keeps 00191 // concurrent readers from interfering with each other. 00192 00193 MSIter(const MeasurementSet& ms, const Block<Int>& sortColumns, 00194 Double timeInterval=0, Bool addDefaultSortColumns=True, 00195 Bool storeSorted=True); 00196 00197 // Same as above with multiple MSs as input. 00198 MSIter(const Block<MeasurementSet>& mss, const Block<Int>& sortColumns, 00199 Double timeInterval=0, Bool addDefaultSortColumns=True, 00200 Bool storeSorted=True); 00201 00202 // Copy construct. This calls the assigment operator. 00203 MSIter(const MSIter & other); 00204 00205 // Destructor 00206 virtual ~MSIter(); 00207 00208 // Assigment. This will reset the iterator to the origin. 00209 MSIter & operator=(const MSIter &other); 00210 00211 //# Members 00212 00213 // Set or reset the time interval to use for iteration. 00214 // You should call origin() to reset the iteration after 00215 // calling this. 00216 void setInterval(Double timeInterval); 00217 00218 // Reset iterator to start of data 00219 void origin(); 00220 00221 // Return False if there is no more data 00222 Bool more() const; 00223 00224 // Advance iterator through data 00225 MSIter & operator++(int); 00226 MSIter & operator++(); 00227 00228 // Return the current Table iteration 00229 Table table() const; 00230 00231 // Return reference to the current MS 00232 const MS& ms() const; 00233 00234 // Return reference to the current ROMSColumns 00235 const ROMSColumns& msColumns() const; 00236 00237 // Return the current MS Id (according to the order in which 00238 // they appeared in the constructor) 00239 Int msId() const; 00240 00241 // Return true if msId has changed since last iteration 00242 Bool newMS() const; 00243 00244 // Return the current ArrayId 00245 Int arrayId() const; 00246 00247 // Return True if ArrayId has changed since last iteration 00248 Bool newArray() const; 00249 00250 // Return the current FieldId 00251 Int fieldId() const; 00252 00253 // Return the current Field Name 00254 const String& fieldName() const; 00255 00256 // Return the current Source Name 00257 const String& sourceName() const; 00258 00259 // Return True if FieldId/Source has changed since last iteration 00260 Bool newField() const; 00261 00262 // Return current SpectralWindow 00263 Int spectralWindowId() const; 00264 00265 // Return True if SpectralWindow has changed since last iteration 00266 Bool newSpectralWindow() const; 00267 00268 // Return current DataDescriptionId 00269 Int dataDescriptionId() const; 00270 00271 // Return True if DataDescriptionId has changed since last iteration 00272 Bool newDataDescriptionId() const; 00273 00274 // Return current PolarizationId 00275 Int polarizationId() const; 00276 00277 // Return True if polarization has changed since last iteration 00278 Bool newPolarizationId() const; 00279 00280 // Return the current phase center as MDirection 00281 const MDirection& phaseCenter() const; 00282 00283 // Return frame for polarization (returns PolFrame enum) 00284 Int polFrame() const; 00285 00286 // Return the frequencies corresponding to the DATA matrix. 00287 const Vector<Double>& frequency() const; 00288 00289 // Return frequency of first channel with reference frame as a Measure. 00290 // The reference frame Epoch is that of the first row, reset it as needed 00291 // for each row. 00292 // The reference frame Position is the average of the antenna positions. 00293 const MFrequency& frequency0() const; 00294 00295 // Return the rest frequency of the specified line as a Measure 00296 const MFrequency& restFrequency(Int line=0) const; 00297 00298 // Return the telescope position (if a known telescope) or the 00299 // position of the first antenna (if unknown) 00300 const MPosition& telescopePosition() const; 00301 00302 // Return the feed configuration/leakage matrix for feed 0 on each antenna 00303 // TODO: CJonesAll can be used instead of this method in all instances 00304 const Vector<SquareMatrix<Complex,2> >& CJones() const; 00305 00306 // Return the feed configuration/leakage matrix for all feeds and antennae 00307 // First axis is antennaId, 2nd axis is feedId. Result of CJones() is 00308 // a reference to the first column of the matrix returned by this method 00309 const Matrix<SquareMatrix<Complex,2> >& CJonesAll() const; 00310 00311 // Return the receptor angle for feed 0 on each antenna. 00312 // First axis is receptor number, 2nd axis is antennaId. 00313 // TODO: receptorAngles() can be used instead of this method 00314 const Matrix<Double>& receptorAngle() const; 00315 00316 // Return the receptor angles for all feeds and antennae 00317 // First axis is a receptor number, 2nd axis is antennaId, 00318 // 3rd axis is feedId. Result of receptorAngle() is just a reference 00319 // to the first plane of the cube returned by this method 00320 const Cube<Double>& receptorAngles() const; 00321 00322 // Return the channel number of the first channel in the DATA. 00323 // (non-zero for reference MS created by VisSet with channel selection) 00324 Int startChan() const; 00325 00326 // Return a string mount identifier for each antenna 00327 const Vector<String>& antennaMounts() const; 00328 00329 // Return a cube containing pairs of coordinate offset for each receptor 00330 // of each feed (values are in radians, coordinate system is fixed with 00331 // antenna and is the same as used to define the BEAM_OFFSET parameter 00332 // in the feed table). The cube axes are receptor, antenna, feed. 00333 const Cube<RigidVector<Double, 2> >& getBeamOffsets() const; 00334 00335 // True if all elements of the cube returned by getBeamOffsets are zero 00336 Bool allBeamOffsetsZero() const; 00337 00338 // Get the spw, start and nchan for all the ms's is this msiter that 00339 // match the frequecy "freqstart-freqStep" and "freqEnd+freqStep" range 00340 00341 void getSpwInFreqRange(Block<Vector<Int> >& spw, 00342 Block<Vector<Int> >& start, 00343 Block<Vector<Int> >& nchan, 00344 Double freqStart, Double freqEnd, 00345 Double freqStep); 00346 00347 //Get the number of actual ms's associated wth this iterator 00348 Int numMS() const; 00349 00350 //Get a reference to the nth ms in the list of ms associated with this 00351 // iterator. If larger than the list of ms's current ms is returned 00352 // So better check wth numMS() before making the call 00353 const MS& ms(const uInt n) const; 00354 00355 protected: 00356 // handle the construction details 00357 void construct(const Block<Int>& sortColumns, Bool addDefaultSortColumns); 00358 // advance the iteration 00359 void advance(); 00360 // set the iteration state 00361 void setState(); 00362 void setMSInfo(); 00363 void setArrayInfo(); 00364 void setFeedInfo(); 00365 void setDataDescInfo(); 00366 void setFieldInfo(); 00367 00368 // Determine if the numbers in r1 are a sorted subset of those in r2 00369 Bool isSubSet(const Vector<uInt>& r1, const Vector<uInt>& r2); 00370 00371 MSIter* This; 00372 Block<MeasurementSet> bms_p; 00373 PtrBlock<TableIterator* > tabIter_p; 00374 Block<Bool> tabIterAtStart_p; 00375 00376 Int nMS_p; 00377 ROMSColumns* msc_p; 00378 Table curTable_p; 00379 Int curMS_p, lastMS_p, curArray_p, lastArray_p, curSource_p; 00380 String curFieldName_p, curSourceName_p; 00381 Int curField_p, lastField_p, curSpectralWindow_p, lastSpectralWindow_p; 00382 Int curPolarizationId_p, lastPolarizationId_p; 00383 Int curDataDescId_p, lastDataDescId_p; 00384 Bool more_p, newMS_p, newArray_p, newField_p, newSpectralWindow_p, 00385 newPolarizationId_p, newDataDescId_p, preselected_p, 00386 timeDepFeed_p, spwDepFeed_p, checkFeed_p; 00387 Int startChan_p; 00388 00389 // Globally control disk storage of SORTED_TABLE 00390 Bool storeSorted_p; 00391 00392 // time selection 00393 Double interval_p; 00394 // channel selection 00395 Block<Int> preselectedChanStart_p,preselectednChan_p; 00396 00397 // columns 00398 ScalarColumn<Int> colArray_p, colDataDesc_p, colField_p; 00399 00400 //cache for access functions 00401 MDirection phaseCenter_p; 00402 Matrix<Double> receptorAnglesFeed0_p; // former receptorAngle_p, 00403 // temporary retained for compatibility 00404 // contain actually a reference to the 00405 // first plane of receptorAngles_p 00406 Cube<Double> receptorAngles_p; 00407 Vector<SquareMatrix<Complex,2> > CJonesFeed0_p; // a temporary reference 00408 // similar to receptorAngle_p 00409 Matrix<SquareMatrix<Complex,2> > CJones_p; 00410 Vector<String> antennaMounts_p; // a string mount identifier for each 00411 // antenna (e.g. EQUATORIAL, ALT-AZ,...) 00412 Cube<RigidVector<Double, 2> > beamOffsets_p;// angular offsets (two values for 00413 // each element of the cube in radians) 00414 // in the antenna coordinate system. 00415 // Cube axes are: receptor, antenna, feed. 00416 Bool allBeamOffsetsZero_p; // True if all elements of beamOffsets_p 00417 // are zero (to speed things up in a 00418 // single beam case) 00419 PolFrame polFrame_p; 00420 Bool freqCacheOK_p; 00421 Vector<Double> frequency_p; 00422 MFrequency frequency0_p; 00423 MFrequency restFrequency_p; 00424 MPosition telescopePosition_p; 00425 00426 MSInterval *timeComp_p; // Points to the time comparator. 00427 // 0 if not using a time interval. 00428 }; 00429 00430 inline Bool MSIter::more() const { return more_p;} 00431 inline Table MSIter::table() const {return curTable_p;} 00432 inline const MS& MSIter::ms() const {return bms_p[curMS_p];} 00433 inline const ROMSColumns& MSIter::msColumns() const { return *msc_p;} 00434 inline Bool MSIter::newMS() const { return newMS_p;} 00435 inline Bool MSIter::newArray() const {return newArray_p;} 00436 inline Bool MSIter::newField() const { return newField_p;} 00437 inline Bool MSIter::newSpectralWindow() const 00438 { return newSpectralWindow_p;} 00439 inline Int MSIter::msId() const { return curMS_p;} 00440 inline Int MSIter::numMS() const { return nMS_p;} 00441 inline Int MSIter::arrayId() const {return curArray_p;} 00442 inline Int MSIter::fieldId() const { return curField_p;} 00443 inline const String& MSIter::fieldName() const { return curFieldName_p;} 00444 inline const String& MSIter::sourceName() const { return curSourceName_p;} 00445 inline Int MSIter::spectralWindowId() const 00446 { return curSpectralWindow_p;} 00447 inline Int MSIter::polarizationId() const {return curPolarizationId_p;} 00448 inline Int MSIter::dataDescriptionId() const {return curDataDescId_p;} 00449 inline Bool MSIter::newPolarizationId() const { return newPolarizationId_p;} 00450 inline Bool MSIter::newDataDescriptionId() const { return newDataDescId_p;} 00451 inline Int MSIter::polFrame() const { return polFrame_p;} 00452 inline const MDirection& MSIter::phaseCenter() const 00453 { return phaseCenter_p; } 00454 inline const MPosition& MSIter::telescopePosition() const 00455 { return telescopePosition_p;} 00456 inline const Vector<SquareMatrix<Complex,2> >& MSIter::CJones() const 00457 { return CJonesFeed0_p;} 00458 inline const Matrix<SquareMatrix<Complex,2> >& MSIter::CJonesAll() const 00459 { return CJones_p;} 00460 inline const Matrix<Double>& MSIter::receptorAngle() const 00461 {return receptorAnglesFeed0_p;} 00462 inline const Cube<Double>& MSIter::receptorAngles() const 00463 {return receptorAngles_p;} 00464 inline const Vector<String>& MSIter::antennaMounts() const 00465 {return antennaMounts_p;} 00466 inline const Cube<RigidVector<Double, 2> >& MSIter::getBeamOffsets() const 00467 {return beamOffsets_p;} 00468 inline Int MSIter::startChan() const {return startChan_p;} 00469 inline Bool MSIter::allBeamOffsetsZero() const {return allBeamOffsetsZero_p;} 00470 00471 } //# NAMESPACE CASACORE - END 00472 00473 #endif