Changeset 1446 for branches/alma/src/STMath.cpp

Timestamp:

11/12/08 17:04:01 (16 years ago)

Author:

TakTsutsumi

Message:

Merged recent updates (since 2007) from nrao-asap

File:

• branches/alma/src/STMath.cpp (modified) (5 diffs)

branches/alma/src/STMath.cpp

@@ -191 +191 @@
     }
     //write out
-    Vector<uChar> flg(msk.shape());
-    convertArray(flg, !msk);
-    flagColOut.put(i, flg);
-    specColOut.put(i, acc.getSpectrum());
-    tsysColOut.put(i, acc.getTsys());
-    intColOut.put(i, acc.getInterval());
-    mjdColOut.put(i, acc.getTime());
-    // we should only have one cycle now -> reset it to be 0
-    // frequency switched data has different CYCLENO for different IFNO
-    // which requires resetting this value
-    cycColOut.put(i, uInt(0));
+    if (acc.state()) {
+      Vector<uChar> flg(msk.shape());
+      convertArray(flg, !msk);
+      flagColOut.put(i, flg);
+      specColOut.put(i, acc.getSpectrum());
+      tsysColOut.put(i, acc.getTsys());
+      intColOut.put(i, acc.getInterval());
+      mjdColOut.put(i, acc.getTime());
+      // we should only have one cycle now -> reset it to be 0
+      // frequency switched data has different CYCLENO for different IFNO
+      // which requires resetting this value
+      cycColOut.put(i, uInt(0));
+    } else {
+      ostringstream oss;
+      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
+      pushLog(String(oss));
+    }
     acc.reset();
   }
@@ -947 +953 @@
 {
 
-
+  
   STSelector sel;
   CountedPtr< Scantable > ws = getScantable(s, false);
   CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
-  CountedPtr< Scantable > calsig, calref, out;
+  CountedPtr< Scantable > calsig, calref, out, out1, out2;
+  Bool nofold=False;
 
   //split the data
@@ -973 +980 @@
   calref = dototalpower(refwcal, ref, tcal=tcal);
 
-  out=dosigref(calsig,calref,smoothref,tsysv,tau);
-
+  out1=dosigref(calsig,calref,smoothref,tsysv,tau); 
+  out2=dosigref(calref,calsig,smoothref,tsysv,tau); 
+
+  Table& tabout1=out1->table();
+  Table& tabout2=out2->table();
+  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
+  ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
+  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
+  Vector<Float> spec; specCol.get(0, spec);
+  uInt nchan = spec.nelements();
+  uInt freqid1; freqidCol1.get(0,freqid1);
+  uInt freqid2; freqidCol2.get(0,freqid2);
+  Double rp1, rp2, rv1, rv2, inc1, inc2;
+  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
+  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
+  if (rp1==rp2) {
+    Double foffset = rv1 - rv2;
+    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
+    if (choffset >= nchan) {
+      cerr<<"out-band frequency switching, no folding"<<endl;
+      nofold = True;
+    }
+  }
+ 
+  if (nofold) {
+    std::vector< CountedPtr< Scantable > > tabs; 
+    tabs.push_back(out1);
+    tabs.push_back(out2);
+    out = merge(tabs);
+  }
+  else { //folding is not implemented yet
+    out = out1;
+  }
+    
   return out;
 }
-
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
@@ -1565 +1603 @@
           id = out->frequencies().addEntry(rp, rv, inc);
           freqidcol.put(k,id);
-          String name,fname;Double rf;
+          //String name,fname;Double rf;
+          Vector<String> name,fname;Vector<Double> rf;
           (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
           id = out->molecules().addEntry(rf, name, fname);
@@ -1981 +2020 @@
   return out;
 }
+
+// Averaging spectra with different channel/resolution
+CountedPtr<Scantable>
+STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
+                     const bool& compel,
+                     const std::vector<bool>& mask,
+                     const std::string& weight,
+                     const std::string& avmode )
+  throw ( casa::AipsError )
+{
+  if ( avmode == "SCAN" && in.size() != 1 )
+    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
+                    "Use merge first."));
+  
+  CountedPtr<Scantable> out ;     // processed result 
+  if ( compel ) {
+    std::vector< CountedPtr<Scantable> > newin ; // input for average process
+    uInt insize = in.size() ;    // number of scantables 
+
+    // TEST: do normal average in each table before IF grouping
+    cout << "Do preliminary averaging" << endl ;
+    vector< CountedPtr<Scantable> > tmpin( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
+      tmpin[itable] = average( v, mask, weight, avmode ) ;
+    }
+
+    // warning
+    cout << "Average spectra with different spectral resolution" << endl ;
+    cout << endl ;
+
+    // temporarily set coordinfo
+    vector<string> oldinfo( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = in[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // columns
+    ScalarColumn<uInt> freqIDCol ;
+    ScalarColumn<uInt> ifnoCol ;
+    ScalarColumn<uInt> scannoCol ;
+
+
+    // check IF frequency coverage
+    //cout << "Check IF settings in each table" << endl ;
+    vector< vector<uInt> > freqid( insize );
+    vector< vector<double> > iffreq( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = tmpin[itable]->nrow() ;
+      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        if ( freqid[itable].size() == freqnrows ) {
+          break ;
+        }
+        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 ) {
+            //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
+            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
+            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] ) ) ;
+          }
+        }
+      }
+    }
+
+    // debug
+    //cout << "IF settings summary:" << endl ;
+    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
+    //cout << "   Table" << i << endl ;
+    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
+    //cout << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
+    //}
+    //}
+    //cout << endl ;
+
+    // IF grouping based on their frequency coverage
+    //cout << "IF grouping based on their frequency coverage" << endl ;
+
+    // IF group 
+    // ifgrp[numgrp][nummember*2]
+    // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
+    //         [table11, freqrow11, table11, freqrow11, ...],
+    //         ...
+    // ifgfreq[numgrp*2] 
+    // ifgfreq = [min0, max0, min1, max1, ...]
+    vector< vector<uInt> > ifgrp ;
+    vector<double> ifgfreq ;
+
+    // parameter for IF grouping
+    // groupmode = OR    retrieve all region 
+    //             AND   only retrieve overlaped region
+    //string groupmode = "AND" ;
+    string groupmode = "OR" ;
+    uInt sizecr = 0 ;
+    if ( groupmode == "AND" ) 
+      sizecr = 2 ;
+    else if ( groupmode == "OR" ) 
+      sizecr = 0 ;
+
+    vector<double> sortedfreq ;
+    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
+      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
+        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
+          sortedfreq.push_back( iffreq[i][j] ) ;
+      }
+    }
+    sort( sortedfreq.begin(), sortedfreq.end() ) ;
+    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
+      ifgfreq.push_back( *i ) ;
+      ifgfreq.push_back( *(i+1) ) ;
+    }
+    ifgrp.resize( ifgfreq.size()/2 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
+        double range0 = iffreq[itable][2*iif] ;
+        double range1 = iffreq[itable][2*iif+1] ;
+        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
+          double fmin = max( range0, ifgfreq[2*j] ) ;
+          double fmax = min( range1, ifgfreq[2*j+1] ) ;
+          if ( fmin < fmax ) {
+            ifgrp[j].push_back( itable ) ;
+            ifgrp[j].push_back( freqid[itable][iif] ) ;
+          }
+        }
+      }
+    }
+    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ; 
+    vector<double>::iterator giter = ifgfreq.begin() ;
+    while( fiter != ifgrp.end() ) {
+      if ( fiter->size() <= sizecr ) {
+        fiter = ifgrp.erase( fiter ) ;
+        giter = ifgfreq.erase( giter ) ;
+        giter = ifgfreq.erase( giter ) ;
+      }
+      else {
+        fiter++ ;
+        advance( giter, 2 ) ;
+      }
+    }
+
+    // freqgrp[numgrp][nummember]
+    // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
+    //           [ifgrp10, ifgrp11, ifgrp12, ...],
+    //           ...
+    // freqrange[numgrp*2]
+    // freqrange = [min0, max0, min1, max1, ...]
+    vector< vector<uInt> > freqgrp ;
+    double freqrange = 0.0 ;
+    uInt grpnum = 0 ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      // Assumed that ifgfreq was sorted
+      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
+        freqgrp[grpnum-1].push_back( i ) ;
+      }
+      else {
+        vector<uInt> grp0( 1, i ) ;
+        freqgrp.push_back( grp0 ) ;
+        grpnum++ ;
+      }
+      freqrange = ifgfreq[2*i+1] ;
+    }
+        
+
+    // print IF groups
+    cout << "IF Group summary: " << endl ;
+    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+        cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+      }
+      cout << endl ;
+    }
+    cout << endl ;
+    
+    // print frequency group
+    cout << "Frequency Group summary: " << endl ;
+    cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
+    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+      cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
+      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
+        cout << freqgrp[i][j] << " " ;
+      }
+      cout << endl ;
+    }
+    cout << endl ;
+
+    // membership check
+    // groups[numtable][numIF][nummembership]
+    vector< vector< vector<uInt> > > groups( insize ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      groups[i].resize( freqid[i].size() ) ;
+    }
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+        uInt tableid = ifgrp[igrp][2*imem] ;
+        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
+        if ( iter != freqid[tableid].end() ) {
+          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
+          groups[tableid][rowid].push_back( igrp ) ;
+        }
+      }
+    }
+
+    // print membership
+    //for ( uInt i = 0 ; i < insize ; i++ ) {
+    //cout << "Table " << i << endl ;
+    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
+    //cout << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
+    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
+    //cout << setw( 2 ) << groups[i][j][k] << " " ;
+    //}
+    //cout << endl ;
+    //}
+    //}
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      tmpin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // Create additional table if needed
+    bool oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    vector< vector<uInt> > addrow( insize ) ;
+    vector<uInt> addtable( insize, 0 ) ;
+    vector<uInt> newtableids( insize ) ;
+    vector<uInt> newifids( insize, 0 ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      //cout << "Table " << setw(2) << itable << ": " ;
+      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
+        //cout << addrow[itable][ifrow] << " " ;
+      }
+      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
+      //cout << "(" << addtable[itable] << ")" << endl ;
+    }
+    newin.resize( insize ) ;
+    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
+    for ( uInt i = 0 ; i < insize ; i++ ) {
+      newtableids[i] = i ;
+    }
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
+        vector<int> freqidlist ;
+        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
+          if ( groups[itable][i].size() > iadd + 1 ) {
+            freqidlist.push_back( freqid[itable][i] ) ;
+          }
+        }
+        stringstream taqlstream ;
+        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
+        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
+          taqlstream << i ;
+          if ( i < freqidlist.size() - 1 )
+            taqlstream << "," ;
+          else
+            taqlstream << "]" ;
+        }
+        string taql = taqlstream.str() ;
+        //cout << "taql = " << taql << endl ;
+        STSelector selector = STSelector() ;
+        selector.setTaQL( taql ) ;
+        add->setSelection( selector ) ;
+        newin.push_back( add ) ;
+        newtableids.push_back( itable ) ;
+        newifids.push_back( iadd + 1 ) ;
+      }
+    }
+
+    // udpate ifgrp
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
+        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
+          if ( groups[itable][ifrow].size() > iadd + 1 ) {
+            uInt igrp = groups[itable][ifrow][iadd+1] ;
+            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
+                ifgrp[igrp][2*imem] = insize + iadd ;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    // print IF groups again for debug
+    //cout << "IF Group summary: " << endl ;
+    //cout << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
+    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
+    //cout << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
+    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
+    //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ; 
+    //}
+    //cout << endl ;
+    //}
+    //cout << endl ;
+
+    // reset SCANNO and IF: IF number is reset by the result of sortation 
+    cout << "All scan number is set to 0" << endl ;
+    //cout << "All IF number is set to IF group index" << endl ;
+    cout << endl ;
+    insize = newin.size() ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      Table &tmpt = newin[itable]->table() ;
+      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
+      scannoCol.attach( tmpt, "SCANNO" ) ;
+      ifnoCol.attach( tmpt, "IFNO" ) ;
+      for ( uInt irow=0 ; irow < rows ; irow++ ) {
+        scannoCol.put( irow, 0 ) ;
+        uInt freqID = freqIDCol( irow ) ;
+        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
+        if ( iter != freqid[newtableids[itable]].end() ) {
+          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
+          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
+        }
+        else {
+          throw(AipsError("IF grouping was wrong in additional tables.")) ;
+        } 
+      }
+    }
+    oldinfo.resize( insize ) ;
+    setInsitu( oldInsitu ) ;
+
+    // temporarily set coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      oldinfo[itable] = coordinfo[0] ;
+      coordinfo[0] = "Hz" ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }
+
+    // save column values in the vector
+    vector< vector<uInt> > freqTableIdVec( insize ) ;
+    vector< vector<uInt> > freqIdVec( insize ) ;
+    vector< vector<uInt> > ifNoVec( insize ) ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      ScalarColumn<uInt> freqIDs ;
+      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
+      ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
+      freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
+      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
+        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
+      }
+      for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
+        freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
+        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
+      }
+    }
+
+    // reset spectra and flagtra: pick up common part of frequency coverage
+    //cout << "Pick common frequency range and align resolution" << endl ;
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      uInt rows = newin[itable]->nrow() ;
+      int nminchan = -1 ;
+      int nmaxchan = -1 ;
+      vector<uInt> freqIdUpdate ; 
+      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
+        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
+        double minfreq = ifgfreq[2*ifno] ;
+        double maxfreq = ifgfreq[2*ifno+1] ;
+        //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
+        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
+        int nchan = abcissa.size() ;
+        double resol = abcissa[1] - abcissa[0] ;
+        //cout << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
+        if ( minfreq <= abcissa[0] ) 
+          nminchan = 0 ;
+        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 ) ;
+        }
+        if ( maxfreq >= abcissa[abcissa.size()-1] )
+          nmaxchan = abcissa.size() - 1 ;
+        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 ) ;
+        }
+        //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
+        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 )
+            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
+        }
+        else {
+          throw(AipsError("Failed to pick up common part of frequency range.")) ;
+        }
+      }
+      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 ) ;
+      }    
+    }
+
+    
+    // reset spectra and flagtra: align spectral resolution
+    //cout << "Align spectral resolution" << endl ;
+    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
+    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
+    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
+      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution 
+      int minchan = INT_MAX ;     // minimum channel number
+      Double refpixref = -1 ;     // reference of 'reference pixel'
+      Double refvalref = -1 ;     // reference of 'reference frequency'
+      Double refinc = -1 ;        // reference frequency resolution 
+      uInt refreqid ;
+      uInt reftable = INT_MAX;
+      // process only if group member > 1
+      if ( ifgrp[igrp].size() > 2 ) {
+        // find minchan and maxdnu in each group
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+          if ( iter != freqIdVec[tableid].end() ) {
+            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+            int nchan = abcissa.size() ;
+            double dnu = abcissa[1] - abcissa[0] ;
+            //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ; 
+            if ( nchan < minchan ) {
+              minchan = nchan ;
+              maxdnu = dnu ;
+              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
+              refreqid = rowid ;
+              reftable = tableid ;
+            }
+          }
+        }
+        // regrid spectra in each group
+        cout << "GROUP " << igrp << endl ;
+        cout << "   Channel number is adjusted to " << minchan << endl ;
+        cout << "   Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
+        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
+          uInt tableid = ifgrp[igrp][2*imem] ;
+          uInt rowid = ifgrp[igrp][2*imem+1] ;
+          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
+          //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
+          //cout << "   regridChannel applied to " ;
+          if ( tableid != reftable ) 
+            refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
+          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
+            uInt tfreqid = freqIdVec[tableid][irow] ;
+            if ( tfreqid == rowid ) {     
+              //cout << irow << " " ;
+              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
+              freqIDCol.put( irow, refreqid ) ;
+              freqIdVec[tableid][irow] = refreqid ;
+            }
+          }
+          //cout << endl ;
+        } 
+      }
+      else {
+        uInt tableid = ifgrp[igrp][0] ;
+        uInt rowid = ifgrp[igrp][1] ;
+        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
+        if ( iter != freqIdVec[tableid].end() ) {
+          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
+          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
+          minchan = abcissa.size() ;
+          maxdnu = abcissa[1] - abcissa[0] ;
+        }
+      }
+      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
+        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
+          if ( maxdnu > gmaxdnu[i] ) {
+            gmaxdnu[i] = maxdnu ;
+            gmemid[i] = igrp ;
+          }
+          break ;
+        }
+      }
+    }
+    cout << endl ;
+
+    // set back coordinfo
+    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
+      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
+      coordinfo[0] = oldinfo[itable] ;
+      newin[itable]->setCoordInfo( coordinfo ) ;
+    }      
+
+    // accumulate all rows into the first table 
+    // NOTE: assumed in.size() = 1
+    vector< CountedPtr<Scantable> > tmp( 1 ) ;
+    if ( newin.size() == 1 )
+      tmp[0] = newin[0] ;
+    else 
+      tmp[0] = merge( newin ) ;
+
+    //return tmp[0] ;
+
+    // average
+    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
+
+    //return tmpout ;
+
+    // combine frequency group
+    cout << "Combine spectra based on frequency grouping" << endl ;
+    cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
+    vector<string> coordinfo = tmpout->getCoordInfo() ;
+    oldinfo[0] = coordinfo[0] ;
+    coordinfo[0] = "Hz" ;
+    tmpout->setCoordInfo( coordinfo ) ;
+    // create proformas of output table
+    stringstream taqlstream ;
+    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
+    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
+      taqlstream << gmemid[i] ;
+      if ( i < gmemid.size() - 1 )
+        taqlstream << "," ;
+      else
+        taqlstream << "]" ;
+    }
+    string taql = taqlstream.str() ;
+    //cout << "taql = " << taql << endl ;
+    STSelector selector = STSelector() ;
+    selector.setTaQL( taql ) ;
+    oldInsitu = insitu_ ;
+    setInsitu( false ) ;
+    out = getScantable( tmpout, false ) ;
+    setInsitu( oldInsitu ) ;
+    out->setSelection( selector ) ;
+    // regrid rows
+    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
+    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
+      uInt ifno = ifnoCol( irow ) ;
+      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
+          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
+          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
+          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
+          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
+          break ;
+        }
+      }
+    }
+    // combine spectra
+    ArrayColumn<Float> specColOut ;
+    specColOut.attach( out->table(), "SPECTRA" ) ;
+    ArrayColumn<uChar> flagColOut ;
+    flagColOut.attach( out->table(), "FLAGTRA" ) ;
+    ScalarColumn<uInt> ifnoColOut ;
+    ifnoColOut.attach( out->table(), "IFNO" ) ;
+    ScalarColumn<uInt> polnoColOut ;
+    polnoColOut.attach( out->table(), "POLNO" ) ;
+    ScalarColumn<uInt> freqidColOut ;
+    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
+    Table &tab = tmpout->table() ;
+    TableIterator iter( tab, "POLNO" ) ;
+    vector< vector<uInt> > sizes( freqgrp.size() ) ;
+    while( !iter.pastEnd() ) {
+      vector< vector<Float> > specout( freqgrp.size() ) ;
+      vector< vector<uChar> > flagout( freqgrp.size() ) ;
+      ArrayColumn<Float> specCols ;
+      specCols.attach( iter.table(), "SPECTRA" ) ;
+      ArrayColumn<uChar> flagCols ;
+      flagCols.attach( iter.table(), "FLAGTRA" ) ;
+      ifnoCol.attach( iter.table(), "IFNO" ) ;
+      ScalarColumn<uInt> polnos ;
+      polnos.attach( iter.table(), "POLNO" ) ;
+      uInt polno = polnos( 0 ) ;
+      //cout << "POLNO iteration: " << polno << endl ;
+      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 ) ;
+              Vector<uChar> flag = flagCols( irow ) ;
+              vector<Float> svec ;
+              spec.tovector( svec ) ;
+              vector<uChar> fvec ;
+              flag.tovector( fvec ) ;
+              //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
+              specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
+              flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
+              //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
+              sizes[igrp][imem] = spec.nelements() ;
+            }
+          }
+        }
+        for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
+          uInt ifout = ifnoColOut( irow ) ;
+          uInt polout = polnoColOut( irow ) ;
+          if ( ifout == gmemid[igrp] && polout == polno ) {
+            // set SPECTRA and FRAGTRA
+            Vector<Float> newspec( specout[igrp] ) ;
+            Vector<uChar> newflag( flagout[igrp] ) ;
+            specColOut.put( irow, newspec ) ;
+            flagColOut.put( irow, newflag ) ;
+            // IFNO renumbering
+            ifnoColOut.put( irow, igrp ) ;
+          }
+        } 
+      }
+      iter++ ;
+    }
+    // update FREQUENCIES subtable
+    vector<bool> updated( freqgrp.size(), false ) ;
+    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
+      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] ) {
+          uInt freqidOut = freqidColOut( irow ) ;
+          //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
+          double refpix ;
+          double refval ;
+          double increm ;
+          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
+          refpix += pixShift ;
+          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
+          updated[igrp] = true ;
+        }
+      }
+    }
+
+    //out = tmpout ;
+
+    coordinfo = tmpout->getCoordInfo() ;
+    coordinfo[0] = oldinfo[0] ;
+    tmpout->setCoordInfo( coordinfo ) ;
+  }
+  else {
+    // simple average
+    out =  average( in, mask, weight, avmode ) ;
+  }
+  
+  return out ;
+}