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$
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1