Changeset 2571

Timestamp:

06/15/12 15:35:25 (13 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

STMath::averageWithinSession is rewritten such that input table is
purely read-only.

File:

• branches/hpc33/src/STMath.cpp (modified) (2 diffs)

branches/hpc33/src/STMath.cpp

@@ -3955 +3955 @@
                                                        masks,
                                                        "TINT" ) ;
+
     // on scan
 //     t0 = mathutil::gettimeofday_sec() ;
@@ -5115 +5116 @@
                                                     string weight )
 {
-    Table ttab = s->table() ;
-    ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
-    uInt nrow = timeCol->nrow() ;
-    Vector<Double> timeSep = timeCol->getColumn() ;
-    Vector<uInt> scannoOrg = s->scanCol_.getColumn() ;
-    delete timeCol ;
-    for ( uInt i = nrow-2 ; i > 0 ; i-- ) {
-      timeSep[i] -= timeSep[i-1] ;
-    }
-    Vector<Double> interval = s->integrCol_.getColumn() ;
-    interval /= 86400.0 ;
-    uInt *newscan = new uInt[nrow] ;
-    Vector<uInt> newscanno( IPosition(1,nrow), newscan, TAKE_OVER ) ; 
-    uInt *p = newscan ;
-    uInt newid = 0 ;
-    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
-      *p = newid ;
-      double gap = 2.0 * timeSep[i+1] / ( interval[i] + interval[i+1] ) ;
+  // prepare output table
+  bool insitu = insitu_ ;
+  insitu_ = false ;
+  CountedPtr<Scantable> a = getScantable( s, true ) ;
+  insitu_ = insitu ;
+  Table &atab = a->table() ;
+  ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
+
+  if ( s->nrow() == 0 ) 
+    return a ;
+
+  // setup RowAccumulator
+  WeightType wtype = stringToWeight( weight ) ;
+  RowAccumulator acc( wtype ) ;
+  Vector<Bool> cmask( mask ) ;
+  acc.setUserMask( cmask ) ;
+
+  vector<string> cols( 1, "POLNO" ) ;
+  STIdxIterAcc iter( s, cols ) ;
+
+  uInt nchan = s->nchan(s->getIF(0)) ;
+  Table ttab = s->table() ;
+  ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
+  Vector<Double> timeVec = timeCol->getColumn() ;
+  delete timeCol ;
+  Vector<Double> interval = s->integrCol_.getColumn() ;
+  uInt nrow = timeVec.nelements() ;
+  uInt outrow = 0 ;
+
+  Vector<uChar> flag( nchan ) ;
+  Vector<Bool> bflag( nchan ) ;
+  Vector<Float> spec( nchan ) ;
+  Vector<Float> tsys( nchan ) ;
+
+  while( !iter.pastEnd() ) {
+
+    Vector<uInt> rows = iter.getRows( SHARE ) ;
+    uInt len = rows.nelements() ;
+
+    Vector<Double> timeSep( len-1 ) ; 
+    for ( uInt i = 0 ; i < len-1 ; i++ ) {
+      timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
+    }
+
+    uInt irow ;
+    uInt jrow ;
+    for ( uInt i = 0 ; i < len-1 ; i++ ) {
+      irow = rows[i] ;
+      jrow = rows[i+1] ;
+      // accumulate data
+      s->flagsCol_.get( irow, flag ) ;
+      convertArray( bflag, flag ) ;
+      s->specCol_.get( irow, spec ) ;
+      tsys.assign( s->tsysCol_( irow ) ) ;
+      if ( !allEQ(bflag,True) ) 
+        acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
+      double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
       //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
       if ( gap > 1.1 ) {
-//         cout << "detected gap between " << i << " and " << i+1 << endl ;
-        newid++ ;
-      }
-      p++ ;
-    }
-    *p = newid ;
-    s->scanCol_.putColumn( newscanno ) ;
-//     double t2 = mathutil::gettimeofday_sec() ;
-    vector< CountedPtr<Scantable> > tmp( 1, s ) ;
-    CountedPtr<Scantable> a = average( tmp, mask, weight, "SCAN" ) ;
-    s->scanCol_.putColumn( scannoOrg ) ;
-//     double t3 = mathutil::gettimeofday_sec() ;
-//     cout << "a.nrow = " << a->nrow() << endl ;
-//     t1 = mathutil::gettimeofday_sec() ;
-//     cout << "   elapsed time for average(): " << t3-t2 << " sec" << endl ;
-    return a ;
-}
+        //cout << "detected gap between " << i << " and " << i+1 << endl ;
+        // put data to output table
+        // reset RowAccumulator
+        if ( acc.state() ) {
+          atab.addRow() ;
+          copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ; 
+          acc.replaceNaN() ;
+          const Vector<Bool> &msk = acc.getMask() ;
+          convertArray( flag, !msk ) ;
+          for (uInt k = 0; k < nchan; ++k) {
+            uChar userFlag = 1 << 7;
+            if (msk[k]==True) userFlag = 0 << 7;
+            flag(k) = userFlag;
+          }
+          a->flagsCol_.put( outrow, flag ) ;
+          a->specCol_.put( outrow, acc.getSpectrum() ) ;
+          a->tsysCol_.put( outrow, acc.getTsys() ) ;
+          a->integrCol_.put( outrow, acc.getInterval() ) ;
+          timeColOut.put( outrow, acc.getTime() ) ;
+          a->cycleCol_.put( outrow, 0 ) ;
+        }
+        acc.reset() ;
+        outrow++ ;
+      }
+    }
+
+    // accumulate and add last data
+    irow = rows[len-1] ;
+    s->flagsCol_.get( irow, flag ) ;
+    convertArray( bflag, flag ) ;
+    s->specCol_.get( irow, spec ) ;
+    tsys.assign( s->tsysCol_( irow ) ) ;
+    if (!allEQ(bflag,True) ) 
+      acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
+    if ( acc.state() ) {
+      atab.addRow() ;
+      copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
+      acc.replaceNaN() ;
+      const Vector<Bool> &msk = acc.getMask() ;
+      convertArray( flag, !msk ) ;
+      for (uInt k = 0; k < nchan; ++k) {
+        uChar userFlag = 1 << 7;
+        if (msk[k]==True) userFlag = 0 << 7;
+        flag(k) = userFlag;
+      }
+      a->flagsCol_.put( outrow, flag ) ;
+      a->specCol_.put( outrow, acc.getSpectrum() ) ;
+      a->tsysCol_.put( outrow, acc.getTsys() ) ;
+      a->integrCol_.put( outrow, acc.getInterval() ) ;
+      timeColOut.put( outrow, acc.getTime() ) ;
+      a->cycleCol_.put( outrow, 0 ) ;
+    }
+    acc.reset() ;
+    outrow++ ;
+
+    iter.next() ;
+  }
+
+  return a ;
+}