Changeset 2466

Timestamp:

04/13/12 21:33:47 (13 years ago)

Author:

Kana Sugimoto

Message:

New Development: No

JIRA Issue: No (a bug fix)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: unit test of sdcal, using averageall=True

Put in Release Notes: No

Module(s):

Description:

Fixed a bug recently found by Malte in asapmath.new_average, with compel=True, or
sdcal with averageall=True.
The bug caused the method/task crash, when the spectra to be averaged over have
IFs with negative and positive frequency increments mixed.

File:

• trunk/src/STMath.cpp (modified) (16 diffs)

trunk/src/STMath.cpp

@@ -3045 +3045 @@
       uInt rows = tmpin[itable]->nrow() ;
       uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
+      freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
+      //ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
         if ( freqid[itable].size() == freqnrows ) {
@@ -3050 +3052 @@
         }
         else {
-          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
-          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
           uInt id = freqIDCol( irow ) ;
           if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
             //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
             vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
+            double lbedge, rbedge;
             freqid[itable].push_back( id ) ;
-            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
-            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
+            lbedge = abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] );
+            rbedge = abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] );
+            iffreq[itable].push_back( min(lbedge, rbedge) ) ;
+            iffreq[itable].push_back( max(lbedge, rbedge) ) ;
           }
         }
@@ -3143 +3146 @@
     // Grouping continuous IF groups (without frequency gap)
     // freqgrp: list of IF group indexes in each frequency group
-    // freqrange: list of minimum and maximum frequency in each frequency group
     // freqgrp[numgrp][nummember]
     // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
@@ -3149 +3151 @@
     //           ...
     //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
-    // freqrange[numgrp*2]
-    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
     vector< vector<uInt> > freqgrp ;
     double freqrange = 0.0 ;
@@ -3388 +3388 @@
         double resol = abcissa[1] - abcissa[0] ;
         //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
-        if ( minfreq <= abcissa[0] ) 
-          nminchan = 0 ;
+        uInt imin, imax;
+        int sigres;
+        if ( resol >= 0. ) {
+          imin = 0;
+          imax = nchan - 1 ;
+          sigres = 1 ;
+        } else {
+          imin = nchan - 1 ;
+          imax = 0 ;
+          sigres = -1 ;
+          resol = abs(resol) ;
+        }
+        if ( minfreq <= abcissa[imin] ) 
+          nminchan = imin ;
         else {
           //double cfreq = ( minfreq - abcissa[0] ) / resol ;
-          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
-          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
+          double cfreq = ( minfreq - abcissa[imin] + 0.5 * resol ) / resol ;
+          nminchan = imin + sigres * (int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 )) ;
         }
-        if ( maxfreq >= abcissa[abcissa.size()-1] )
-          nmaxchan = abcissa.size() - 1 ;
+        if ( maxfreq >= abcissa[imax] )
+          nmaxchan = imax;
         else {
           //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
-          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
-          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
+          double cfreq = ( abcissa[imax] - maxfreq + 0.5 * resol ) / resol ;
+          nmaxchan = imax - sigres * (int(cfreq) +( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 )) ;
         }
-        //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
+        //os << "freqrange = [" << abcissa[nminchan] << ", " << abcissa[nmaxchan] << "]"<< LogIO::POST; 
+        if ( nmaxchan < nminchan ) {
+          int tmp = nmaxchan ;
+          nmaxchan = nminchan ;
+          nminchan = tmp ;
+        }
         if ( nmaxchan > nminchan ) {
           newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
           int newchan = nmaxchan - nminchan + 1 ;
-          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
+          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan ) {
             freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
+
+            // Update before nminchan is lost
+            uInt freqId = freqIdVec[itable][irow] ;
+            Double refpix ;
+            Double refval ;
+            Double increment ;
+            newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
+            //refval = refval - ( refpix - nminchan ) * increment ;
+            refval = abcissa[nminchan] ;
+            refpix = 0 ;
+            newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
+          }
         }
         else {
@@ -3413 +3442 @@
         }
       }
-      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
-        uInt freqId = freqIdUpdate[i] ;
-        Double refpix ;
-        Double refval ;
-        Double increment ;
+//       for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
+//      uInt freqId = freqIdUpdate[i] ;
+//      Double refpix ;
+//      Double refval ;
+//      Double increment ;
         
-        // update row
-        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
-        refval = refval - ( refpix - nminchan ) * increment ;
-        refpix = 0 ;
-        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
-      }    
+//      // update row
+//      newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ; 
+//      refval = refval - ( refpix - nminchan ) * increment ;
+//      refpix = 0 ;
+//      newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
+//       }    
     }
 
@@ -3461 +3490 @@
             if ( nchan < minchan ) {
               minchan = nchan ;
-              maxdnu = dnu ;
+              maxdnu = abs(dnu) ;
               newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
               refreqid = rowid ;
               reftable = tableid ;
+              // Spectra are reversed when dnu < 0
+              if (dnu < 0)
+                refpixref = minchan - 1 -refpixref ;
             }
           }
@@ -3479 +3511 @@
           //os << "   regridChannel applied to " ;
           //if ( tableid != reftable ) 
-          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
+          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, maxdnu ) ;
           for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
             uInt tfreqid = freqIdVec[tableid][irow] ;
@@ -3500 +3532 @@
           vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
           minchan = abcissa.size() ;
-          maxdnu = abcissa[1] - abcissa[0] ;
+          maxdnu = abs( abcissa[1] - abcissa[0]) ;
         }
       }
@@ -3538 +3570 @@
     // combine frequency group
     os << "Combine spectra based on frequency grouping" << LogIO::POST ;
-    os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
+    //os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
+    os << "IFNO is renumbered as IF group ID (see above)" << LogIO::POST ;
     vector<string> coordinfo = tmpout->getCoordInfo() ;
     oldinfo[0] = coordinfo[0] ;
@@ -3570 +3603 @@
           vector<double> abcissa = tmpout->getAbcissa( irow ) ;
           double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
-          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
+          int nchan = (int) abs( bw / gmaxdnu[igrp] ) ;
+          // All spectra will have positive frequency increments
           tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
           break ;
@@ -3593 +3627 @@
     cols[0] = String("POLNO") ;
     TableIterator iter( tab, cols ) ;
-    bool done = false ;
     vector< vector<uInt> > sizes( freqgrp.size() ) ;
     while( !iter.pastEnd() ) {
@@ -3644 +3677 @@
 //       }
       // get a list of number of channels for each frequency group member
-      if ( !done ) {
-        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
-          sizes[igrp].resize( freqgrp[igrp].size() ) ;
-          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
-            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
-              uInt ifno = ifnoCol( irow ) ;
-              if ( ifno == freqgrp[igrp][imem] ) {
-                Vector<Float> spec = specCols( irow ) ;
-                sizes[igrp][imem] = spec.nelements() ;
-                break ;
-              }                
-            }
-          }
-        }
-        done = true ;
+      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+        sizes[igrp].resize( freqgrp[igrp].size() ) ;
+        for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
+          for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
+            uInt ifno = ifnoCol( irow ) ;
+            if ( ifno == freqgrp[igrp][imem] ) {
+              Vector<Float> spec = specCols( irow ) ;
+              sizes[igrp][imem] = spec.nelements() ;
+              break ;
+            }
+          }
+        }
       }
       // combine spectra
@@ -3705 +3735 @@
           specColOut.put( irow, newspec ) ;
           flagColOut.put( irow, newflag ) ;
-          // IFNO renumbering
+          // IFNO renumbering (renumbered as frequency group ID)
           ifnoColOut.put( irow, igrp ) ; 
         }
@@ -3716 +3746 @@
       uInt index = 0 ;
       uInt pixShift = 0 ;
+
       while ( freqgrp[igrp][index] != gmemid[igrp] ) {
         pixShift += sizes[igrp][index++] ;
       }
       for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
-        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
+        //if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
+        if ( ifnoColOut( irow ) == igrp && !updated[igrp] ) {
           uInt freqidOut = freqidColOut( irow ) ;
           //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
@@ -3727 +3759 @@
           double increm ;
           out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
+          if (increm < 0)
+            refpix = sizes[igrp][index] - 1 - refpix ; // reversed
           refpix += pixShift ;
-          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
+          out->frequencies().setEntry( refpix, refval, gmaxdnu[igrp], freqidOut ) ;
           updated[igrp] = true ;
         }