Changeset 2519

Timestamp:

05/15/12 15:23:48 (13 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: Defined STIdxIter and its python interface

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Defined STIdxIter module that contains STIdxIterAcc class.
Python interface for STIdxIterAcc class also defined.
STMath::almacal is rewritten using STIdxIterAcc.

Location:

branches/hpc33/src

Files:

• . (modified) (1 prop)
• CMakeLists.txt (modified) (2 diffs)
• FillerWrapper.h (modified) (1 prop)
• SConscript (modified) (1 prop)
• STIdxIter.cpp (added)
• STIdxIter.h (added)
• STMath.cpp (modified) (8 diffs)
• STMath.h (modified) (1 diff)
• python_asap.cpp (modified) (1 diff)
• python_asap.h (modified) (1 diff)

branches/hpc33/src/CMakeLists.txt

@@ -60 +60 @@
      ${SRCDIR}/STUpgrade.cpp 
      ${SRCDIR}/STGrid.cpp
+     ${SRCDIR}/STIdxIter.cpp
      ${SRCDIR}/Templates.cpp )
 
@@ -81 +82 @@
      ${SRCDIR}/python_LogSink.cpp
      ${SRCDIR}/python_STGrid.cpp
+     ${SRCDIR}/python_Iterator.cpp
      ${SRCDIR}/python_asap.cpp )
 

branches/hpc33/src/STMath.cpp

@@ -54 +54 @@
 #include "STSelector.h"
 #include "Accelerator.h"
+#include "STIdxIter.h"
 
 using namespace casa;
@@ -3860 +3861 @@
   }
   else {
+//     double t0 = mathutil::gettimeofday_sec() ;
     vector<bool> masks = s->getMask( 0 ) ;
     
     // off scan
     STSelector sel = STSelector() ;
-    vector<int> types ;
-    types.push_back( SrcType::PSOFF ) ;
+    vector<int> types( 1, SrcType::PSOFF ) ;
+    //types.push_back( SrcType::PSOFF ) ;
     sel.setTypes( types ) ;
     s->setSelection( sel ) ;
@@ -3874 +3876 @@
     CountedPtr<Scantable> soff = getScantable( s, false ) ;
     Table ttab = soff->table() ;
-    ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
-    uInt nrow = timeCol.nrow() ;
-    Vector<Double> timeSep( nrow - 1 ) ;
+//     String tmpname = File::newUniqueName( "./", "temp" ).baseName() ;
+//     Table ttab = s->table().copyToMemoryTable( tmpname, False ) ;
+    ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
+    uInt nrow = timeCol->nrow() ;
+    Vector<Double> timeSep = timeCol->getColumn() ;
+    delete timeCol ;
+    for ( uInt i = nrow-2 ; i > 0 ; i-- ) {
+      timeSep[i] -= timeSep[i-1] ;
+    }
+    Vector<Double> interval = soff->integrCol_.getColumn() ;
+    interval /= 86400.0 ;
+    Block<uInt> gaplist( 100 ) ;
+    uInt glidx = 0 ;
+    uInt glsize = gaplist.size() ;
     for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
-      timeSep[i] = timeCol(i+1) - timeCol(i) ;
-    }
-    ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
-    Vector<Double> interval = intervalCol.getColumn() ;
-    interval /= 86400.0 ;
-    ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
-    vector<uInt> glist ;
-    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
-      double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
+      double gap = 2.0 * timeSep[i+1] / ( interval[i] + interval[i+1] ) ;
       //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
       if ( gap > 1.1 ) {
-        glist.push_back( i ) ;
-      }
-    }
-    Vector<uInt> gaplist( glist ) ;
+//         cout << "detected gap between " << i << " and " << i+1 << endl ;
+        gaplist[glidx] = i ;
+        glidx++ ;
+      }
+      if ( glidx >= glsize ) {
+//         cout << "resize gaplist" << endl ;
+        glsize += 100 ;
+        gaplist.resize( glsize ) ;
+      }
+    }
     //cout << "gaplist = " << gaplist << endl ;
     uInt newid = 0 ;
     for ( uInt i = 0 ; i < nrow ; i++ ) {
-      scanCol.put( i, newid ) ;
+      soff->scanCol_.put( i, newid ) ;
       if ( i == gaplist[newid] ) {
         newid++ ;
@@ -3907 +3918 @@
     s->unsetSelection() ;
     sel.reset() ;
-    types.clear() ;
+    //types.clear() ;
+//     double t1 = mathutil::gettimeofday_sec() ;
+//     cout << "elapsed time for off averaging: " << t1-t0 << " sec" << endl ;
     
     // on scan
+//     t0 = mathutil::gettimeofday_sec() ;
     bool insitu = insitu_ ;
     insitu_ = false ;
     CountedPtr<Scantable> out = getScantable( s, true ) ;
     insitu_ = insitu ;
-    types.push_back( SrcType::PSON ) ;
+    //types.push_back( SrcType::PSON ) ;
+    types[0] = SrcType::PSON ;
     sel.setTypes( types ) ;
     s->setSelection( sel ) ;
@@ -3920 +3935 @@
     s->unsetSelection() ;
     sel.reset() ;
-    types.clear() ;
-    
+    //types.clear() ;
+//     t1 = mathutil::gettimeofday_sec() ;
+//     cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
+
     // process each on scan
+//     t0 = mathutil::gettimeofday_sec() ;
 //     for ( int i = 0 ; i < out->nrow() ; i++ ) {
 //       vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
@@ -3928 +3946 @@
 //     }
 
-    // selection first, then process
-    vector<unsigned int> ifs = out->getIFNos() ;
-    vector<unsigned int> beams = out->getBeamNos() ;
-    vector<unsigned int> pols = out->getPolNos() ;
-    vector<int> v( 1 ) ;
-    for ( vector<unsigned int>::iterator i = ifs.begin() ; \
-          i != ifs.end() ; i++ ) {
-      v[0] = *i ;
-      sel.setIFs( v ) ;
-      for ( vector<unsigned int>::iterator j = pols.begin() ; \
-            j != pols.end() ; j++ ) {
-        v[0] = *j ;
-        sel.setPolarizations( v ) ;
-        for ( vector<unsigned int>::iterator k = beams.begin() ; \
-              k != beams.end() ; k++ ) {
-          v[0] = *k ;
-          sel.setBeams( v ) ;
-          out->setSelection( sel ) ;
-          aoff->setSelection( sel ) ;
-
-          // calibrate data
-          calibrateALMA( out, aoff ) ;
-
-          out->unsetSelection() ;
-          aoff->unsetSelection() ;
-          sel.reset() ;
-        }
-      }
-    }
+    // using STIdxIterAcc 
+    vector<string> cols( 3 ) ;
+    cols[0] = "BEAMNO" ;
+    cols[1] = "POLNO" ;
+    cols[2] = "IFNO" ;
+    STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
+    while ( !iter->pastEnd() ) {
+      Vector<uInt> ids = iter->current() ;
+      stringstream ss ;
+      ss << "SELECT FROM $1 WHERE "
+         << "BEAMNO==" << ids[0] << "&&"
+         << "POLNO==" << ids[1] << "&&"
+         << "IFNO==" << ids[2] ;
+      sel.setTaQL( ss.str() ) ;
+      aoff->setSelection( sel ) ;
+      Vector<uInt> rows = iter->getRows() ;
+      calibrateALMA( out, aoff, rows ) ;
+      aoff->unsetSelection() ;
+      sel.reset() ;
+      iter->next() ;
+    }
+    delete iter ;
+
+//     t1 = mathutil::gettimeofday_sec() ;
+//     cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
 
     // flux unit
@@ -5087 +5101 @@
 
 void STMath::calibrateALMA( CountedPtr<Scantable>& on,
-                                        CountedPtr<Scantable>& off )
-{
-//   string reftime = on->getTime( index ) ;
+                            CountedPtr<Scantable>& off )
+{
+  //   string reftime = on->getTime( index ) ;
   ROTableColumn timeCol( on->table(), "TIME" ) ;
   ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
@@ -5115 +5129 @@
     }
     on->setSpectrum( sp, index ) ;
+  }
+}
+
+void STMath::calibrateALMA( CountedPtr<Scantable>& on,
+                            CountedPtr<Scantable>& off,
+                            Vector<uInt>& rows )
+{
+//   string reftime = on->getTime( index ) ;
+  ROTableColumn timeCol( on->table(), "TIME" ) ;
+  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
+  vector<float> sp( on->nchan( on->getIF(rows[0]) ) ) ;
+  unsigned int spsize = sp.size() ;
+  // I know that the data is contiguous
+  const uInt *p = rows.data() ;
+  for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
+    //int index = rows[irow] ;
+    //double reftime = timeCol.asdouble(index) ;
+    double reftime = timeCol.asdouble(*p) ;
+    vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
+    //vector<float> spec = on->getSpectrum( index ) ;
+    //Vector<Float> tsys = tsysCol( index ) ;
+    vector<float> spec = on->getSpectrum( *p ) ;
+    Vector<Float> tsys = tsysCol( *p ) ;
+    // ALMA Calibration
+    // 
+    // Ta* = Tsys * ( ON - OFF ) / OFF
+    //
+    // 2010/01/07 Takeshi Nakazato
+    unsigned int tsyssize = tsys.nelements() ;
+    for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
+      float tscale = 0.0 ;
+      if ( tsyssize == spsize ) 
+        tscale = tsys[j] ;
+      else 
+        tscale = tsys[0] ;
+      float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
+      sp[j] = v ;
+    }
+    //on->setSpectrum( sp, index ) ;
+    on->setSpectrum( sp, *p ) ;
+    p++ ;
   }
 }

branches/hpc33/src/STMath.h

@@ -423 +423 @@
   void calibrateALMA( casa::CountedPtr<Scantable>& on,
                       casa::CountedPtr<Scantable>& off ) ;
+  void calibrateALMA( casa::CountedPtr<Scantable>& on,
+                      casa::CountedPtr<Scantable>& off,
+                      casa::Vector<casa::uInt>& rows ) ;
   vector<float> getFSCalibratedSpectra( casa::CountedPtr<Scantable>& sig,
                                         casa::CountedPtr<Scantable>& ref,

branches/hpc33/src/python_asap.cpp

@@ -85 +85 @@
   asap::python::python_SrcType();
   asap::python::python_STGrid();
+  asap::python::python_Iterator();
 
 #ifndef HAVE_LIBPYRAP

branches/hpc33/src/python_asap.h

@@ -53 +53 @@
     void python_SrcType();
     void python_STGrid();
+    void python_Iterator();
 
   } // python