Changeset 2542

Timestamp:

05/22/12 10:17:06 (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...

More speedup of STMath::average().

• Reduced number of call of RowAccumulator::replaceNaN().
• Performance of RowAccumulator::add() is improved.

Location:

branches/hpc33/src

Files:

• RowAccumulator.cpp (modified) (2 diffs)
• RowAccumulator.h (modified) (1 diff)
• STMath.cpp (modified) (7 diffs)

branches/hpc33/src/RowAccumulator.cpp

@@ -81 +81 @@
                                  const Double time)
 {
-  doAddSpectrum(v, m, tsys, interval, time, False);
-  doAddSpectrum(v, m, tsys, interval, time, True);  // CAS-2776
+//   doAddSpectrum(v, m, tsys, interval, time, False);
+//   doAddSpectrum(v, m, tsys, interval, time, True);  // CAS-2776
+  doAddSpectrum2( v, m, tsys, interval, time ) ;
 }
 
@@ -112 +113 @@
 }
 
+void RowAccumulator::doAddSpectrum2(const Vector<Float>& v,
+                                    const Vector<Bool>& m,
+                                    const Vector<Float>& tsys,
+                                    const Double interval,
+                                    const Double time)
+{
+  const MaskedArray<Float> vadd(v, m);
+  const MaskedArray<Float> vaddN(v, !m);
+  Float totalWeight = getTotalWeight( vadd, tsys, interval, time, False ) ;
+  Float totalWeightNoMask = getTotalWeight( vaddN, tsys, interval, time, True ) ;
+  // process spectrum with input mask and spectrum with inverted mask 
+  // together
+  // prepare pointers
+  Bool vD, mD ;
+  const Float *v_p = v.getStorage( vD ) ;
+  const Float *v_wp = v_p ;
+  const Bool *m_p = m.getStorage( mD ) ;
+  const Bool *m_wp = m_p ;
+
+  Bool spD, spND, wgtD, wgtND, nD, nND ;
+  Float *sp_p = spectrum_.getRWArrayStorage( spD ) ;
+  Float *sp_wp = sp_p ;
+  Float *wgt_p = weightSum_.getRWArrayStorage( wgtD ) ;
+  Float *wgt_wp = wgt_p ;
+  Float *spN_p = spectrumNoMask_.getRWArrayStorage( spND ) ;
+  Float *spN_wp = spN_p ;
+  Float *wgtN_p = weightSumNoMask_.getRWArrayStorage( wgtND ) ;
+  Float *wgtN_wp = wgtN_p ;
+  uInt *n_p = n_.getRWArrayStorage( nD ) ;
+  uInt *n_wp = n_p ;
+  uInt *nN_p = nNoMask_.getRWArrayStorage( nND ) ;
+  uInt *nN_wp = nN_p ;
+
+  // actual processing
+  uInt len = m.nelements() ;
+  // loop over channels
+  for ( uInt i = 0 ; i < len ; i++ ) {
+    // masks for spectrum_ and spectrumNoMask_ are not needed since 
+    // it is initialized as True for all channels and those are constant
+    if ( *m_wp ) {
+      // mask is True
+      // add v * totalWeight to spectrum_
+      // add totalWeight to weightSum_
+      // increment n_
+      *sp_wp += *v_wp * totalWeight ;
+      *wgt_wp += totalWeight ;
+      *n_wp += 1 ;
+    }
+    else {
+      // mask is False
+      // add v * totalWeightNoMask to spectrumNoMask_
+      // add totalWeightNoMask to weightSumNoMask_
+      // increment nNoMask_
+      *spN_wp += *v_wp * totalWeightNoMask ;
+      *wgtN_wp += totalWeightNoMask ;
+      *nN_wp += 1 ;
+    }
+    sp_wp++ ;
+    wgt_wp++ ;
+    n_wp++ ;
+    spN_wp++ ;
+    wgtN_wp++ ;
+    nN_wp++ ;
+    v_wp++ ;
+    m_wp++ ;
+  }
+
+  // free storage
+  spectrum_.putArrayStorage( sp_p, spD ) ;
+  weightSum_.putArrayStorage( wgt_p, wgtD ) ;
+  spectrumNoMask_.putArrayStorage( spN_p, spND ) ;
+  weightSumNoMask_.putArrayStorage( wgtN_p, wgtND ) ;
+  n_.putArrayStorage( n_p, nD ) ;
+  nNoMask_.putArrayStorage( nN_p, nND ) ;
+
+  v.freeStorage( v_p, vD ) ; 
+  m.freeStorage( m_p, mD ) ;
+}
+
 Float RowAccumulator::getTotalWeight(const MaskedArray<Float>& data,
                                      const Vector<Float>& tsys,

branches/hpc33/src/RowAccumulator.h

@@ -109 +109 @@
                      const casa::Double time,
                      const casa::Bool inverseMask);
+  void doAddSpectrum2(const casa::Vector<casa::Float>& v,
+                      const casa::Vector<casa::Bool>& m,
+                      const casa::Vector<casa::Float>& tsys,
+                      const casa::Double interval,
+                      const casa::Double time);
   casa::Float getTotalWeight(const casa::MaskedArray<casa::Float>& data,
                              const casa::Vector<casa::Float>& tsys,

branches/hpc33/src/STMath.cpp

@@ -164 +164 @@
   // use STIdxIterExAcc instead of TableIterator
   STIdxIterExAcc iter( in[0], cols ) ;
+//   double t2 = 0 ;
+//   double t3 = 0 ;
 //   TableIterator iter(baset, cols);
 //   int count = 0 ;
@@ -229 +231 @@
 
     // in[0] is already selected by TableItertor
-//     specCol.attach(subt,"SPECTRA");
-//     flagCol.attach(subt,"FLAGTRA");
-//     tsysCol.attach(subt,"TSYS");
-//     intCol.attach(subt,"INTERVAL");
-//     mjdCol.attach(subt,"TIME");
     specCol.attach(baset,"SPECTRA");
     flagCol.attach(baset,"FLAGTRA");
@@ -258 +255 @@
       mjdCol.get(k, time);
       // spectrum has to be added last to enable weighting by the other values                             
+//       t2 = mathutil::gettimeofday_sec() ;
       acc.add(spec, !bflag, tsys, inter, time);
-    }
-
-    // If there exists a channel at which all the input spectra are masked,                                
-    // spec has 'nan' values for that channel and it may affect the following                              
-    // processes. To avoid this, replacing 'nan' values in spec with                                       
-    // weighted-mean of all spectra in the following line.                                                 
-    // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)                                                  
-    acc.replaceNaN();
+//       t3 += mathutil::gettimeofday_sec() - t2 ;
+      
+    }
+
 
     // in[0] is already selected by TableIterator so that index is 
@@ -337 +331 @@
       intCol.attach(subt,"INTERVAL");
       mjdCol.attach(subt,"TIME");
-//       Vector<Float> spec,tsys;
-//       Vector<uChar> flag;
-//       Double inter,time;
       for (uInt k = 0; k < subt.nrow(); ++k ) {
         flagCol.get(k, flag);
@@ -354 +345 @@
         mjdCol.get(k, time);
         // spectrum has to be added last to enable weighting by the other values
+//         t2 = mathutil::gettimeofday_sec() ;
         acc.add(spec, !bflag, tsys, inter, time);
-      }
-
+//         t3 += mathutil::gettimeofday_sec() - t2 ;
+      }
+
+    }
+    const Vector<Bool>& msk = acc.getMask();
+    if ( allEQ(msk, False) ) {
+      uint n = rowstodelete.nelements();
+      rowstodelete.resize(n+1, True);
+      rowstodelete[n] = i;
+      continue;
+    }
+    //write out
+    if (acc.state()) {
       // If there exists a channel at which all the input spectra are masked, 
       // spec has 'nan' values for that channel and it may affect the following 
@@ -363 +366 @@
       // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
       acc.replaceNaN();
-    }
-    const Vector<Bool>& msk = acc.getMask();
-    if ( allEQ(msk, False) ) {
-      uint n = rowstodelete.nelements();
-      rowstodelete.resize(n+1, True);
-      rowstodelete[n] = i;
-      continue;
-    }
-    //write out
-    if (acc.state()) {
+
       Vector<uChar> flg(msk.shape());
       convertArray(flg, !msk);
@@ -413 +407 @@
 //    t1 = mathutil::gettimeofday_sec() ;
 //    cout << "elapsed time for average(): " << t1-t0 << " sec" << endl ;
+//    cout << "   elapsed time for acc.add(): " << t3 << " sec" << endl ;
 
   return out;