casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
MSIter.h
Go to the documentation of this file.
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: MSIter.h 20998 2010-11-17 07:10:12Z gervandiepen $
00027 
00028 #ifndef MS_MSITER_H
00029 #define MS_MSITER_H
00030 
00031 #include <casa/aips.h>
00032 #include <casa/Arrays/Matrix.h>
00033 #include <casa/Arrays/Cube.h>
00034 #include <ms/MeasurementSets/MeasurementSet.h>
00035 #include <measures/Measures/MFrequency.h>
00036 #include <measures/Measures/MDirection.h>
00037 #include <measures/Measures/MPosition.h>
00038 #include <tables/Tables/ScalarColumn.h>
00039 #include <casa/Utilities/Compare.h>
00040 #include <casa/BasicSL/String.h>
00041 #include <scimath/Mathematics/SquareMatrix.h>
00042 #include <scimath/Mathematics/RigidVector.h>
00043 
00044 namespace casa { //# NAMESPACE CASA - 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   MSIter(const MeasurementSet& ms, const Block<Int>& sortColumns, 
00186          Double timeInterval=0, Bool addDefaultSortColumns=True);
00187 
00188   // Same as above with multiple MSs as input.
00189   MSIter(const Block<MeasurementSet>& mss, const Block<Int>& sortColumns, 
00190          Double timeInterval=0, Bool addDefaultSortColumns=True);
00191 
00192   // Copy construct. This calls the assigment operator.
00193   MSIter(const MSIter & other);
00194 
00195   // Destructor
00196   virtual ~MSIter();
00197   
00198   // Assigment. This will reset the iterator to the origin.
00199   MSIter & operator=(const MSIter &other);
00200 
00201   //# Members
00202  
00203   // Set or reset the time interval to use for iteration.
00204   // You should call origin() to reset the iteration after 
00205   // calling this.
00206   void setInterval(Double timeInterval);
00207  
00208   // Reset iterator to start of data
00209   void origin();
00210  
00211   // Return False if there is no more data
00212   Bool more() const;
00213 
00214   // Advance iterator through data
00215   MSIter & operator++(int);
00216   MSIter & operator++();
00217 
00218   // Return the current Table iteration
00219   Table table() const;
00220 
00221   // Return reference to the current MS
00222   const MS& ms() const;
00223 
00224   // Return reference to the current ROMSColumns
00225   const ROMSColumns& msColumns() const;
00226 
00227   // Return the current MS Id (according to the order in which 
00228   // they appeared in the constructor)
00229   Int msId() const;
00230 
00231   // Return true if msId has changed since last iteration
00232   Bool newMS() const;
00233 
00234   // Return the current ArrayId
00235   Int arrayId() const;
00236 
00237   // Return True if ArrayId has changed since last iteration
00238   Bool newArray() const;
00239 
00240   // Return the current FieldId
00241   Int fieldId() const;
00242 
00243   // Return the current Field Name
00244   const String& fieldName() const;
00245 
00246   // Return the current Source Name
00247   const String& sourceName() const;
00248 
00249   // Return True if FieldId/Source has changed since last iteration
00250   Bool newField() const;
00251 
00252   // Return current SpectralWindow
00253   Int spectralWindowId() const;
00254 
00255   // Return True if SpectralWindow has changed since last iteration
00256   Bool newSpectralWindow() const;
00257 
00258   // Return current DataDescriptionId
00259   Int dataDescriptionId() const;
00260 
00261   // Return True if DataDescriptionId has changed since last iteration
00262   Bool newDataDescriptionId() const;
00263 
00264   // Return current PolarizationId
00265   Int polarizationId() const;
00266 
00267   // Return True if polarization has changed since last iteration
00268   Bool newPolarizationId() const;
00269 
00270   // Return the current phase center as MDirection
00271   const MDirection& phaseCenter() const;
00272 
00273   // Return frame for polarization (returns PolFrame enum)
00274   Int polFrame() const;
00275 
00276   // Return the frequencies corresponding to the DATA matrix.
00277   const Vector<Double>& frequency() const;
00278 
00279   // Return frequency of first channel with reference frame as a Measure.
00280   // The reference frame Epoch is that of the first row, reset it as needed
00281   // for each row.
00282   // The reference frame Position is the average of the antenna positions.
00283   const MFrequency& frequency0() const;
00284 
00285   // Return the rest frequency of the specified line as a Measure
00286   const MFrequency& restFrequency(Int line=0) const;
00287 
00288   // Return the telescope position (if a known telescope) or the 
00289   // position of the first antenna (if unknown)
00290   const MPosition& telescopePosition() const;
00291 
00292   // Return the feed configuration/leakage matrix for feed 0 on each antenna
00293   // TODO: CJonesAll can be used instead of this method in all instances
00294   const Vector<SquareMatrix<Complex,2> >& CJones() const;
00295 
00296   // Return the feed configuration/leakage matrix for all feeds and antennae
00297   // First axis is antennaId, 2nd axis is feedId. Result of CJones() is
00298   // a reference to the first column of the matrix returned by this method
00299   const Matrix<SquareMatrix<Complex,2> >& CJonesAll() const;
00300 
00301   // Return the receptor angle for feed 0 on each antenna.
00302   // First axis is receptor number, 2nd axis is antennaId.
00303   // TODO: receptorAngles() can be used instead of this method
00304   const Matrix<Double>& receptorAngle() const;
00305 
00306   // Return the receptor angles for all feeds and antennae
00307   // First axis is a receptor number, 2nd axis is antennaId,
00308   // 3rd axis is feedId. Result of receptorAngle() is just a reference
00309   // to the first plane of the cube returned by this method
00310   const Cube<Double>& receptorAngles() const;
00311 
00312   // Return the channel number of the first channel in the DATA.
00313   // (non-zero for reference MS created by VisSet with channel selection)
00314   Int startChan() const;
00315 
00316   // Return a string mount identifier for each antenna
00317   const Vector<String>& antennaMounts() const;
00318 
00319   // Return a cube containing pairs of coordinate offset for each receptor
00320   // of each feed (values are in radians, coordinate system is fixed with 
00321   // antenna and is the same as used to define the BEAM_OFFSET parameter 
00322   // in the feed table). The cube axes are receptor, antenna, feed. 
00323   const Cube<RigidVector<Double, 2> >& getBeamOffsets() const;
00324 
00325   // True if all elements of the cube returned by getBeamOffsets are zero
00326   Bool allBeamOffsetsZero() const;
00327 
00328   // Get the spw, start  and nchan for all the ms's is this msiter that 
00329   // match the frequecy "freqstart-freqStep" and "freqEnd+freqStep" range
00330   
00331   void getSpwInFreqRange(Block<Vector<Int> >& spw, 
00332                          Block<Vector<Int> >& start, 
00333                          Block<Vector<Int> >& nchan, 
00334                          Double freqStart, Double freqEnd, 
00335                          Double freqStep);
00336 
00337   //Get the number of actual ms's associated wth this iterator
00338   Int numMS() const;
00339 
00340   //Get a reference to the nth ms in the list of ms associated with this 
00341   // iterator. If larger than the list of ms's current ms is returned
00342   // So better check wth numMS() before making the call
00343   const MS& ms(const uInt n) const;
00344 
00345 protected:
00346   // handle the construction details
00347   void construct(const Block<Int>& sortColumns, Bool addDefaultSortColumns);
00348   // advance the iteration
00349   void advance();
00350   // set the iteration state
00351   void setState();
00352   void setMSInfo();
00353   void setArrayInfo();
00354   void setFeedInfo();
00355   void setDataDescInfo();
00356   void setFieldInfo();
00357 
00358 // Determine if the numbers in r1 are a sorted subset of those in r2
00359   Bool isSubSet(const Vector<uInt>& r1, const Vector<uInt>& r2);
00360 
00361   MSIter* This;
00362   Block<MeasurementSet> bms_p;
00363   PtrBlock<TableIterator* > tabIter_p;
00364   Block<Bool> tabIterAtStart_p;
00365 
00366   Int nMS_p;
00367   ROMSColumns* msc_p;
00368   Table curTable_p;
00369   Int curMS_p, lastMS_p, curArray_p, lastArray_p, curSource_p;
00370   String curFieldName_p, curSourceName_p;
00371   Int curField_p, lastField_p, curSpectralWindow_p, lastSpectralWindow_p;
00372   Int curPolarizationId_p, lastPolarizationId_p;
00373   Int curDataDescId_p, lastDataDescId_p;
00374   Bool more_p, newMS_p, newArray_p, newField_p, newSpectralWindow_p, 
00375     newPolarizationId_p, newDataDescId_p, preselected_p,
00376     timeDepFeed_p, spwDepFeed_p, checkFeed_p;
00377   Int startChan_p;
00378 
00379   // time selection
00380   Double interval_p;
00381   // channel selection
00382   Block<Int> preselectedChanStart_p,preselectednChan_p;
00383   
00384   // columns
00385   ROScalarColumn<Int> colArray_p, colDataDesc_p, colField_p;
00386 
00387   //cache for access functions
00388   MDirection phaseCenter_p;
00389   Matrix<Double> receptorAnglesFeed0_p; // former receptorAngle_p,
00390                                    // temporary retained for compatibility
00391                                    // contain actually a reference to the
00392                                    // first plane of receptorAngles_p
00393   Cube<Double> receptorAngles_p;
00394   Vector<SquareMatrix<Complex,2> > CJonesFeed0_p; // a temporary reference
00395                                    // similar to receptorAngle_p
00396   Matrix<SquareMatrix<Complex,2> > CJones_p;
00397   Vector<String>  antennaMounts_p; // a string mount identifier for each
00398                                    // antenna (e.g. EQUATORIAL, ALT-AZ,...)
00399   Cube<RigidVector<Double, 2> > beamOffsets_p;// angular offsets (two values for
00400                                    // each element of the cube in radians) 
00401                                    // in the antenna coordinate system. 
00402                                    // Cube axes are: receptor, antenna, feed.
00403   Bool allBeamOffsetsZero_p;       // True if all elements of beamOffsets_p
00404                                    // are zero (to speed things up in a 
00405                                    // single beam case)
00406   PolFrame polFrame_p;
00407   Bool freqCacheOK_p;
00408   Vector<Double> frequency_p;
00409   MFrequency frequency0_p;
00410   MFrequency restFrequency_p;
00411   MPosition telescopePosition_p;
00412 
00413   MSInterval *timeComp_p;          // Points to the time comparator.
00414                                    // NULL if not using a time interval.
00415 };
00416 
00417 inline Bool MSIter::more() const { return more_p;}
00418 inline Table MSIter::table() const {return curTable_p;}
00419 inline const MS& MSIter::ms() const {return bms_p[curMS_p];}
00420 inline const ROMSColumns& MSIter::msColumns() const { return *msc_p;}
00421 inline Bool MSIter::newMS() const { return newMS_p;}
00422 inline Bool MSIter::newArray() const {return newArray_p;}
00423 inline Bool MSIter::newField() const { return newField_p;}
00424 inline Bool MSIter::newSpectralWindow() const 
00425 { return newSpectralWindow_p;}
00426 inline Int MSIter::msId() const { return curMS_p;}
00427 inline Int MSIter::numMS() const { return nMS_p;}
00428 inline Int MSIter::arrayId() const {return curArray_p;}
00429 inline Int MSIter::fieldId() const { return curField_p;}
00430 inline const String& MSIter::fieldName() const { return curFieldName_p;}
00431 inline const String& MSIter::sourceName() const { return curSourceName_p;}
00432 inline Int MSIter::spectralWindowId() const 
00433 { return curSpectralWindow_p;}
00434 inline Int MSIter::polarizationId() const {return curPolarizationId_p;}
00435 inline Int MSIter::dataDescriptionId() const {return curDataDescId_p;}
00436 inline Bool MSIter::newPolarizationId() const { return newPolarizationId_p;}
00437 inline Bool MSIter::newDataDescriptionId() const { return newDataDescId_p;}
00438 inline Int MSIter::polFrame() const { return polFrame_p;}
00439 inline const MDirection& MSIter::phaseCenter() const 
00440 { return phaseCenter_p; }
00441 inline const MPosition& MSIter::telescopePosition() const
00442 { return telescopePosition_p;}
00443 inline const Vector<SquareMatrix<Complex,2> >& MSIter::CJones() const  
00444 { return CJonesFeed0_p;}
00445 inline const Matrix<SquareMatrix<Complex,2> >& MSIter::CJonesAll() const
00446 { return CJones_p;}
00447 inline const Matrix<Double>& MSIter::receptorAngle() const 
00448 {return receptorAnglesFeed0_p;}
00449 inline const Cube<Double>& MSIter::receptorAngles() const
00450 {return receptorAngles_p;}
00451 inline const Vector<String>& MSIter::antennaMounts() const
00452 {return antennaMounts_p;}
00453 inline const Cube<RigidVector<Double, 2> >& MSIter::getBeamOffsets() const
00454 {return beamOffsets_p;}
00455 inline Int MSIter::startChan() const {return startChan_p;}
00456 inline Bool MSIter::allBeamOffsetsZero() const {return allBeamOffsetsZero_p;}
00457 
00458 } //# NAMESPACE CASA - END
00459 
00460 #endif