Changeset 2550

Timestamp:

05/29/12 11:24:15 (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...

Clean up code.

Location:

branches/hpc33/src

Files:

• MSFiller.cpp (modified) (4 diffs)
• STIdxIter.cpp (modified) (17 diffs)
• STMath.cpp (modified) (9 diffs)
• STMath.h (modified) (1 diff)

branches/hpc33/src/MSFiller.cpp

@@ -304 +304 @@
     // set dummy epoch
     mf.set( currentTime ) ;
-
-    // initialize dirType
-    //dirType = MDirection::N_Types ;
 
     //
@@ -957 +954 @@
                     pointingDirection.xyPlane( idx-1 ), pointingDirection.xyPlane( idx ) ) ;
     }
-    //mf.resetEpoch( currentTime ) ;
     mf.set( currentTime ) ;
     Quantum< Vector<Double> > tmp( d.column( 0 ), Unit( "rad" ) ) ;
     if ( dirType != MDirection::J2000 ) {
-      //MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
       dir = toj2000( tmp ).getAngle( "rad" ).getValue() ;
     }
@@ -968 +963 @@
     }
     if ( dirType != MDirection::AZELGEO ) {
-      //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
-      //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZEL, mf ) ) ;
       azel = toazel( tmp ).getAngle( "rad" ).getValue() ;
     }
@@ -981 +974 @@
   {
     dir = sourceDir.getAngle( "rad" ).getValue() ;
-    mf.resetEpoch( currentTime ) ;
-    //MDirection::Convert toazel( dirType, MDirection::Ref( MDirection::AZELGEO, mf ) ) ;
+    mf.set( currentTime ) ;
     azel = toazel( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
     if ( dirType != MDirection::J2000 ) {
-      //MDirection::Convert toj2000( dirType, MDirection::Ref( MDirection::J2000, mf ) ) ;
       dir = toj2000( Quantum< Vector<Double> >( dir, Unit("rad") ) ).getAngle( "rad" ).getValue() ;
     }

branches/hpc33/src/STIdxIter.cpp

@@ -16 +16 @@
   prod_m.resize( nfield_m ) ;
   idx_m.resize( nfield_m ) ;
-//   prod_m[0] = 1 ;
-//   idx_m[0] = 0 ;
-//   for ( uInt i = 1 ; i < nfield_m ; i++ ) {
-//     prod_m[i] = shape[i-1] * prod_m[i-1] ;
-//     idx_m[i] = 0 ;
-//   }
-//   maxiter_m = prod_m[nfield_m-1] * shape[nfield_m-1] ;
   prod_m[nfield_m-1] = 1 ;
   idx_m[nfield_m-1] = 0 ;
@@ -42 +35 @@
   niter_m++ ;
   uInt x = niter_m ;
-//   for ( Int i = nfield_m-1 ; i >= 0 ; i-- ) {
   for ( Int i = 0 ; i < nfield_m ; i++ ) {
     idx_m[i] = x / prod_m[i] ;
@@ -136 +128 @@
   // ...
   storage_m.resize( arr_m.nelements()+nrow_m ) ;
-//   uInt *p = storage_m.storage() + arr_m.nelements() ;
+  // initialize constant temporary storage
   uInt *p = storage_m.storage() ;
   for ( uInt i = 0 ; i < nrow_m ; i++ ) {
@@ -142 +134 @@
     p++ ;
   }
+  // len_m layout
+  // len[0]: length of temporary storage that incidates whole rows (constant)
+  // len[1]: number of matched row for column 0 selection
+  // len[2]: number of matched row for column 1 selection
+  // ...
   len_m.resize( ncol_m+1 ) ;
   p = len_m.storage() ;
@@ -149 +146 @@
 //     cout << "len[" << i << "]=" << len_m[i] << endl ;
   }
+  // skip_m layout
+  // skip_m[0]: True if column 0 is filled by unique value
+  // skip_m[1]: True if column 1 is filled by unique value
+  // ...
+  skip_m.resize( ncol_m ) ;
+  for ( uInt i = 0 ; i < ncol_m ; i++ ) {
+    skip_m[i] = (Bool)(idxlist_m[i].size()==1) ;
+//     cout << "skip_m[" << i << "]=" << skip_m[i] << endl ;
+  }
   prev_m = iter_m->current() ;
   p = prev_m.storage() ;
@@ -156 +162 @@
   }
 //   cout << "prev_m=" << Vector<uInt>(IPosition(1,ncol_m),prev_m.storage()) << endl ;
-  skip_m.resize( ncol_m ) ;
-  for ( uInt i = 0 ; i < ncol_m ; i++ ) {
-    skip_m[i] = (Bool)(idxlist_m[i].size()==1) ;
-//     cout << "skip_m[" << i << "]=" << skip_m[i] << endl ;
-  }
 }
 
@@ -175 +176 @@
 //   cout << "v=" << Vector<uInt>(IPosition(1,v.nelements()),v.storage(),SHARE) << endl ;
 //   cout << "c=" << c << endl ;
-//   if ( c < 0 ) {
   if ( c > ncol_m-1 ) {
-//     pos_m[0] = len_m[0] ;
-//     uInt offset = 0 ;
-//     while( offset < ncol_m && skip_m[offset] )
-//       offset++ ;
     pos_m[0] = len_m[ncol_m] ;
     Int offset = ncol_m - 1 ;
@@ -189 +185 @@
     return Vector<uInt>( pos_m, storage_m.storage()+offset*nrow_m, policy ) ;
   }
-//   Int offset = c + 1 ;
-//   while( offset < ncol_m && skip_m[offset] ) 
-//     offset++ ;
   Int offset = c - 1 ;
   while( offset >= 0 && skip_m[offset] )
@@ -200 +193 @@
 //   cout << "len_m[c+1]=" << len_m[c+1] << endl ;
 //   cout << "base=" << Vector<uInt>(IPosition(1,len_m[c+1]),base,SHARE) << endl ;
-//   for ( Int i = c ; i >= 0 ; i-- ) {
-//     base = updateStorage( i, base, idxlist_m[i][v[i]] ) ;
-  for ( Int i = c+1 ; i <= ncol_m ; i++ ) {
-    base = updateStorage( i, base, idxlist_m[i-1][v[i-1]] ) ;
+  for ( Int i = c ; i < ncol_m ; i++ ) {
+    base = updateStorage( i, base, idxlist_m[i][v[i]] ) ;
 //     cout << "len_m[" << i << "]=" << len_m[i] << endl ;
 //     cout << "base=" << Vector<uInt>(IPosition(1,len_m[i]),base,SHARE) << endl ;
   }
-//   pos_m[0] = len_m[0] ;
   pos_m[0] = len_m[ncol_m] ;
 //   cout << "pos_m=" << pos_m << endl ;
@@ -216 +206 @@
 Int ArrayIndexIteratorAcc::isChanged( Block<uInt> &idx )
 {
-//   Int i = ncol_m - 1 ;
-//   while( i >= 0 && idx[i] == prev_m[i] )
-//     i-- ;
   Int i = 0 ;
   while( i < ncol_m && idx[i] == prev_m[i] )
@@ -229 +216 @@
                                             uInt &v )
 {
+  // CAUTION: 
+  // indexes for storage_m and len_m differ from index for skip_m by 1
+  // (skip_m[0] corresponds to storage_m[1] and len_m[1]) since there 
+  // is additional temporary storage at first segment in storage_m 
   uInt *p ;
-//   if ( skip_m[icol] ) {
-  if ( skip_m[icol-1] ) {
+  if ( skip_m[icol] ) {
     // skip update, just update len_m[icol] and pass appropriate pointer
 //     cout << "skip " << icol << endl ;
     p = base ;
-//     len_m[icol] = len_m[icol+1] ;
-    len_m[icol] = len_m[icol-1] ;
+    len_m[icol+1] = len_m[icol] ;
   }
   else {
 //     cout << "update " << icol << endl ;
     uInt len = 0 ;
-    p = storage_m.storage() + icol * nrow_m ;
+    p = storage_m.storage() + (icol+1) * nrow_m ;
     uInt *work = p ;
     Bool b ;
@@ -247 +236 @@
     long offset = 0 ;
     // warr_p points a first element of (icol)-th column
-//     const uInt *warr_p = arr_p + icol * nrow_m ;
-    const uInt *warr_p = arr_p + (icol-1) * nrow_m ;
-//     for ( uInt i = 0 ; i < len_m[icol+1] ; i++ ) {
-    for ( uInt i = 0 ; i < len_m[icol-1] ; i++ ) {
+    const uInt *warr_p = arr_p + icol * nrow_m ;
+    for ( uInt i = 0 ; i < len_m[icol] ; i++ ) {
       // increment warr_p by (*(base)-*(base-1))
       warr_p += *base - offset ;
@@ -267 +254 @@
     }
     arr_m.freeStorage( arr_p, b ) ;
-    len_m[icol] = len ;
+    len_m[icol+1] = len ;
   }
   return p ;
@@ -464 +451 @@
 {
   convertArray( out, in ) ;
-  //uInt len = in.nelements() ;
-  //Vector<Int> tmp = in.copy() ;
-  //uInt n = genSort( tmp, Sort::Ascending, Sort::QuickSort|Sort::NoDuplicates ) ;
-  //val.resize( n ) ;
-  //for ( uInt i = 0 ; i < n ; i++ ) {
-  //  val[i] = tmp[i] ;
-//     cout << "val[" << i << "]=" << val[i] << endl ;
-  //}
-  //out.resize( len ) ;
-  //uInt idx = -1 ;
-  //cout << "idx=" << idx << endl ;
-  //int idx2 = (Int)idx ;
-  //cout << "idx2=" << idx2 << endl ;
-  //idx = (uInt) idx2 ;
-  //cout << "idx=" << idx << endl ;
-  //if ( n == 1 ) {
-  //  cout << "n=1" << endl ;
-  //  out = 0 ;
-  //}
-  //else if ( n == 2 ) {
-  //  cout << "n=2" << endl ;
-  //  for ( uInt i = 0 ; i < len ; i++ ) {
-  //    out[i] = (in[i] == val[0]) ? 0 : 1 ;
-  //  }
-  //}
-  //else {
-  //  cout << "n=" << n << endl ;
-  //  map<Int,uInt> m ;
-  // for( uInt i = 0 ; i < n ; i++ )
-  //    m[val[i]] = i ;
-  // for ( uInt i = 0 ; i < len ; i++ ) {
-     //for ( uInt j = 0 ; j < n ; j++ ) {
-      //  if ( in[i] == val[j] ) {
-      //    out[i] = j ;
-      //    break ;
-      //  }
-      //}
-  //    out[i] = m[in[i]] ;
-  //  }
-  // }
 } 
 
@@ -518 +465 @@
 //     cout << "val[" << i << "]=" << val[i] << endl ;
   }
-  //out.resize( len ) ;
   if ( n == 1 ) {
     //cout << "n=1" << endl ;
@@ -535 +481 @@
       m[val[i]] = i ;
     for ( uInt i = 0 ; i < len ; i++ ) {
-      //for ( uInt j = 0 ; j < n ; j++ ) {
-      //  if ( in[i] == val[j] ) {
-      //    out[i] = j ;
-      //    break ;
-      //  }
-      //}
       out[i] = m[in[i]] ;
     }
@@ -548 +488 @@
 Int STIdxIterExAcc::getSrcType()
 {
-  //if ( srctype_m.nelements() > 0 ) 
   if ( srctypeid_m >= 0 ) 
-    //return srctype_m[iter_m->current()[srctypeid_m]] ;
     return (Int)(iter_m->current()[srctypeid_m]) ;
   else

branches/hpc33/src/STMath.cpp

@@ -184 +184 @@
     tout.addRow();
 //     t4 = mathutil::gettimeofday_sec() ;
-    //TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
-    //TableCopy::copyRows(tout, baset, outrowCount, rows[0], 1);
     // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are 
     // overwritten in the following process
@@ -235 +233 @@
 //     outrowCount += rowNum ;
 
-// merge loop 
+    // merge loop 
     uInt i = outrowCount ;
-//  for (uInt i=0; i < tout.nrow(); ++i) {
-
-    // in[0] is already selected by TableItertor
+    // in[0] is already selected by iterator
     specCol.attach(baset,"SPECTRA");
     flagCol.attach(baset,"FLAGTRA");
@@ -248 +244 @@
     Vector<uChar> flag;
     Double inter,time;
-//     for (uInt k = 0; k < subt.nrow(); ++k ) {
+
     for (uInt l = 0; l < rows.nelements(); ++l ) {
       uInt k = rows[l] ;
@@ -273 +269 @@
     // in[0] is already selected by TableIterator so that index is 
     // started from 1
-    //for ( int j=0; j < int(in.size()); ++j ) {
     for ( int j=1; j < int(in.size()); ++j ) {
       const Table& tin = in[j]->table();
@@ -362 +357 @@
     const Vector<Bool>& msk = acc.getMask();
     if ( allEQ(msk, False) ) {
-      //uint n = rowstodelete.nelements();
-      //rowstodelete.resize(n+1, True);
-      //rowstodelete[n] = i;
       rowstodelB[nrowdel] = i ;
       nrowdel++ ;
@@ -408 +400 @@
   }
 
-  //if (rowstodelete.nelements() > 0) {
   if ( nrowdel > 0 ) {
     Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ;
@@ -5315 +5306 @@
 }
 
-void STMath::calibrateALMA( CountedPtr<Scantable>& on,
-                            CountedPtr<Scantable>& off )
-{
-  if ( on->nrow() == 0 )
-    return ;
-  //   string reftime = on->getTime( index ) ;
-  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
-  Vector<Double> timeVec = timeCol.getColumn() ;
-  //ROTableColumn timeCol( on->table(), "TIME" ) ;
-  timeCol.attach( on->table(), "TIME" ) ;
-  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
-  vector<float> sp( on->nchan( on->getIF(0) ) ) ;
-  unsigned int spsize = sp.size() ;
-  for ( int index = 0 ; index < on->nrow() ; index++ ) {
-    double reftime = timeCol.asdouble(index) ;
-    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
-    vector<float> spec = on->getSpectrum( index ) ;
-    Vector<Float> tsys = tsysCol( index ) ;
-    // 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 ) ;
-  }
-}
-
-void STMath::calibrateALMA( CountedPtr<Scantable>& on,
-                            CountedPtr<Scantable>& off,
-                            Vector<uInt>& rows )
-{
-  // if rows is empty, just return
-  if ( rows.nelements() == 0 )
-    return ;
-//   string reftime = on->getTime( index ) ;
-  ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
-  Vector<Double> timeVec = timeCol.getColumn() ;
-  timeCol.attach( 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++ ) {
-    double reftime = timeCol.asdouble(*p) ;
-    Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
-    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, *p ) ;
-    p++ ;
-  }
-}
-
 void STMath::calibrateALMA( CountedPtr<Scantable>& out,
                             const CountedPtr<Scantable>& on,
@@ -5409 +5322 @@
   ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
   unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
-  //Vector<Float> spoff( spsize ) ;
-  //Vector<Float> spec( spsize ) ;
   // I know that the data is contiguous
   const uInt *p = rows.data() ;
@@ -5417 +5328 @@
     Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
     Vector<Float> spec = on->specCol_( *p ) ;
-    //getSpectrumFromTime2( spoff, reftime, timeVec, off, "linear" ) ;
-    //on->specCol_.get( *p, spec ) ;
     Vector<Float> tsys = tsysCol( *p ) ;
     // ALMA Calibration

branches/hpc33/src/STMath.h

@@ -422 +422 @@
                                       casa::CountedPtr<Scantable>& off,
                                       int index ) ;
-  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 ) ;
   void calibrateALMA( casa::CountedPtr<Scantable>& out,
                       const casa::CountedPtr<Scantable>& on,