Changeset 2595 for trunk/src


Ignore:
Timestamp:
07/11/12 12:29:28 (12 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-4010

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: sdcal unit test

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Rewrote STMath::new_average based on new algorithm. To do that,
Scantable::regridChannel is re-defined.


Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STMath.cpp

    r2580 r2595  
    30723072                    "Use merge first."));
    30733073 
    3074   // 2012/02/17 TN
    3075   // Since STGrid is implemented, average doesn't consider direction
    3076   // when accumulating
    3077   // check if OTF observation
    3078 //   String obstype = in[0]->getHeader().obstype ;
    3079 //   Double tol = 0.0 ;
    3080 //   if ( obstype.find( "OTF" ) != String::npos ) {
    3081 //     tol = TOL_OTF ;
    3082 //   }
    3083 //   else {
    3084 //     tol = TOL_POINT ;
    3085 //   }
    3086 
    30873074  CountedPtr<Scantable> out ;     // processed result
    30883075  if ( compel ) {
    … …  
    30903077    uInt insize = in.size() ;    // number of input scantables
    30913078
    3092     // TEST: do normal average in each table before IF grouping
    3093     os << "Do preliminary averaging" << LogIO::POST ;
    3094     vector< CountedPtr<Scantable> > tmpin( insize ) ;
     3079    // setup newin
     3080    bool oldInsitu = insitu_ ;
     3081    setInsitu( false ) ;
     3082    newin.resize( insize ) ;
    30953083    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3096       vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
    3097       tmpin[itable] = average( v, mask, weight, avmode ) ;
    3098     }
     3084      newin[itable] = getScantable( in[itable], false ) ;
     3085    }
     3086    setInsitu( oldInsitu ) ;
    30993087
    31003088    // warning
    … …  
    31073095      oldinfo[itable] = coordinfo[0] ;
    31083096      coordinfo[0] = "Hz" ;
    3109       tmpin[itable]->setCoordInfo( coordinfo ) ;
    3110     }
    3111 
    3112     // columns
    3113     ScalarColumn<uInt> freqIDCol ;
    3114     ScalarColumn<uInt> ifnoCol ;
    3115     ScalarColumn<uInt> scannoCol ;
    3116 
     3097      newin[itable]->setCoordInfo( coordinfo ) ;
     3098    }
     3099
     3100    ostringstream oss ;
    31173101
    31183102    // check IF frequency coverage
    … …  
    31343118    vector< vector<double> > iffreq( insize ) ;
    31353119    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3136       uInt rows = tmpin[itable]->nrow() ;
    3137       uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
    3138       freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
    3139       //ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
    3140       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
    3141         if ( freqid[itable].size() == freqnrows ) {
    3142           break ;
    3143         }
    3144         else {
    3145           uInt id = freqIDCol( irow ) ;
    3146           if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
    3147             //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
    3148             vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
    3149             double lbedge, rbedge;
    3150             freqid[itable].push_back( id ) ;
    3151             lbedge = abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] );
    3152             rbedge = abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] );
    3153             iffreq[itable].push_back( min(lbedge, rbedge) ) ;
    3154             iffreq[itable].push_back( max(lbedge, rbedge) ) ;
    3155           }
    3156         }
     3120      Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
     3121      vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
     3122      for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
     3123            i != uniqueFreqId.end() ; i++ ) {
     3124        //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
     3125        uInt target = 0 ;
     3126        while ( freqIds[target] != *i )
     3127          target++ ;
     3128        vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
     3129        freqid[itable].push_back( *i ) ;
     3130        double incr = abs( abcissa[1] - abcissa[0] ) ;
     3131        iffreq[itable].push_back( (*min_element(abcissa.begin(),abcissa.end()))-0.5*incr ) ;
     3132        iffreq[itable].push_back( (*max_element(abcissa.begin(),abcissa.end()))+0.5*incr ) ;
    31573133      }
    31583134    }
    31593135
    31603136    // debug
    3161     //os << "IF settings summary:" << endl ;
    3162     //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
    3163     //os << "   Table" << i << endl ;
    3164     //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
    3165     //os << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
    3166     //}
    3167     //}
    3168     //os << endl ;
    3169     //os.post() ;
     3137//     os << "IF settings summary:" << endl ;
     3138//     for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
     3139//       os << "   Table" << i << endl ;
     3140//       for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
     3141//         os << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
     3142//       }
     3143//     }
     3144//     os << endl ;
     3145//     os.post() ;
    31703146
    31713147    // IF grouping based on their frequency coverage
    3172     // ifgrp: list of table index and FREQ_ID for all members in each IF group
    3173     // ifgfreq: list of minimum and maximum frequency in each IF group
    3174     // ifgrp[numgrp][nummember*2]
    3175     // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
    3176     //         [table10, freqrow10, table11, freqrow11, ...],
    3177     //         ...
    3178     //         [tablen0, freqrown0, tablen1, freqrown1, ...]]
    3179     // ifgfreq[numgrp*2]
    3180     // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
     3148    // ifgrp: number of member in each IF group
     3149    // ifgrp[numgrp]
     3150    // ifgrp: [n0, n1, ...]
    31813151    //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
    3182     vector< vector<uInt> > ifgrp ;
    3183     vector<double> ifgfreq ;
    31843152
    31853153    // parameter for IF grouping
    … …  
    32023170    }
    32033171    sort( sortedfreq.begin(), sortedfreq.end() ) ;
    3204     for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
    3205       ifgfreq.push_back( *i ) ;
    3206       ifgfreq.push_back( *(i+1) ) ;
    3207     }
    3208     ifgrp.resize( ifgfreq.size()/2 ) ;
     3172    vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
    32093173    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    32103174      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
    32113175        double range0 = iffreq[itable][2*iif] ;
    32123176        double range1 = iffreq[itable][2*iif+1] ;
    3213         for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
    3214           double fmin = max( range0, ifgfreq[2*j] ) ;
    3215           double fmax = min( range1, ifgfreq[2*j+1] ) ;
    3216           if ( fmin < fmax ) {
    3217             ifgrp[j].push_back( itable ) ;
    3218             ifgrp[j].push_back( freqid[itable][iif] ) ;
    3219           }
    3220         }
    3221       }
    3222     }
    3223     vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
    3224     vector<double>::iterator giter = ifgfreq.begin() ;
    3225     while( fiter != ifgrp.end() ) {
    3226       if ( fiter->size() <= sizecr ) {
    3227         fiter = ifgrp.erase( fiter ) ;
    3228         giter = ifgfreq.erase( giter ) ;
    3229         giter = ifgfreq.erase( giter ) ;
    3230       }
    3231       else {
    3232         fiter++ ;
    3233         advance( giter, 2 ) ;
     3177        for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
     3178          double fmin = max( range0, sortedfreq[j] ) ;
     3179          double fmax = min( range1, sortedfreq[j+1] ) ;
     3180          if ( fmin < fmax ) {
     3181            ifgrp[j]++ ;
     3182          }
     3183        }
    32343184      }
    32353185    }
    … …  
    32423192    //           ...
    32433193    //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
     3194    // grprange[2*numgrp]
     3195    // grprange: [fmin0,fmax0,fmin1,fmax1,...]
    32443196    vector< vector<uInt> > freqgrp ;
    3245     double freqrange = 0.0 ;
    3246     uInt grpnum = 0 ;
    3247     for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
    3248       // Assumed that ifgfreq was sorted
    3249       if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
    3250         freqgrp[grpnum-1].push_back( i ) ;
    3251       }
    3252       else {
    3253         vector<uInt> grp0( 1, i ) ;
    3254         freqgrp.push_back( grp0 ) ;
    3255         grpnum++ ;
    3256       }
    3257       freqrange = ifgfreq[2*i+1] ;
    3258     }
    3259 
    3260 
    3261     // print IF groups
    3262     ostringstream oss ;
    3263     oss << "IF Group summary: " << endl ;
    3264     oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
    3265     for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
    3266       oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
    3267       for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
    3268         oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
    3269       }
    3270       oss << endl ;
    3271     }
    3272     oss << endl ;
    3273     os << oss.str() << LogIO::POST ;
    3274    
     3197    vector<double> grprange ;
     3198    vector<uInt> grpedge ;
     3199    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
     3200      if ( ifgrp[igrp] <= sizecr ) {
     3201        grpedge.push_back( igrp ) ;
     3202      }
     3203    }
     3204    grpedge.push_back( ifgrp.size() ) ;
     3205    uInt itmp = 0 ;
     3206    for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
     3207      int n = grpedge[i] - itmp ;
     3208      if ( n > 0 ) {
     3209        vector<uInt> members( n ) ;
     3210        for ( int j = 0 ; j < n ; j++ ) {
     3211          members[j] = itmp+j ;
     3212        }
     3213        freqgrp.push_back( members ) ;
     3214        grprange.push_back( sortedfreq[itmp] ) ;
     3215        grprange.push_back( sortedfreq[grpedge[i]] ) ;
     3216      }
     3217      itmp += n + 1 ;
     3218    }
     3219
    32753220    // print frequency group
    32763221    oss.str("") ;
    32773222    oss << "Frequency Group summary: " << endl ;
    3278     oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
     3223    oss << "   GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
    32793224    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
    3280       oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
    3281       for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
    3282         oss << freqgrp[i][j] << " " ;
    3283       }
     3225      oss << "   GROUP " << setw( 2 ) << i << ": [" << grprange[2*i] << "," << grprange[2*i+1] << "]" ;
    32843226      oss << endl ;
    32853227    }
    … …  
    32873229    os << oss.str() << LogIO::POST ;
    32883230
    3289     // membership check
    3290     // groups: list of IF group indexes whose frequency range overlaps with
    3291     //         that of each table and IF
    3292     // groups[numtable][numIF][nummembership]
    3293     // groups: [[[grp, grp,...], [grp, grp,...],...],
    3294     //          [[grp, grp,...], [grp, grp,...],...],
     3231    // groups: list of frequency group index whose frequency range overlaps
     3232    //         with that of each table and IF
     3233    // groups[numtable][numIF]
     3234    // groups: [[grpx, grpy,...],
     3235    //          [grpa, grpb,...],
    32953236    //          ...
    3296     //          [[grp, grp,...], [grp, grp,...],...]]
    3297     vector< vector< vector<uInt> > > groups( insize ) ;
     3237    //          [grpk, grpm,...]]
     3238    vector< vector<uInt> > groups( insize ) ;
    32983239    for ( uInt i = 0 ; i < insize ; i++ ) {
    32993240      groups[i].resize( freqid[i].size() ) ;
    33003241    }
    3301     for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
    3302       for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
    3303         uInt tableid = ifgrp[igrp][2*imem] ;
    3304         vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
    3305         if ( iter != freqid[tableid].end() ) {
    3306           uInt rowid = distance( freqid[tableid].begin(), iter ) ;
    3307           groups[tableid][rowid].push_back( igrp ) ;
    3308         }
    3309       }
    3310     }
     3242    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3243      for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
     3244        double minf = iffreq[itable][2*ifreq] ;
     3245        uInt groupid ;
     3246        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3247          vector<uInt> memberlist = freqgrp[igrp] ;
     3248          if ( (minf >= grprange[2*igrp]) && (minf <= grprange[2*igrp+1]) ) {
     3249            groupid = igrp ;
     3250            break ;
     3251          }
     3252        }
     3253        groups[itable][ifreq] = groupid ;
     3254      }
     3255    }
     3256                                                         
    33113257
    33123258    // print membership
    3313 //     oss.str("") ;
    3314 //     for ( uInt i = 0 ; i < insize ; i++ ) {
    3315 //       oss << "Table " << i << endl ;
    3316 //       for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
    3317 //         oss << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
    3318 //         for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
    3319 //           oss << setw( 2 ) << groups[i][j][k] << " " ;
    3320 //         }
    3321 //         oss << endl ;
    3322 //       }
    3323 //     }
    3324 //     os << oss.str() << LogIO::POST ;
    3325 
    3326     // set back coordinfo
    3327     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3328       vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
    3329       coordinfo[0] = oldinfo[itable] ;
    3330       tmpin[itable]->setCoordInfo( coordinfo ) ;
    3331     }
    3332 
    3333     // Create additional table if needed
    3334     bool oldInsitu = insitu_ ;
    3335     setInsitu( false ) ;
    3336     vector< vector<uInt> > addrow( insize ) ;
    3337     vector<uInt> addtable( insize, 0 ) ;
    3338     vector<uInt> newtableids( insize ) ;
    3339     vector<uInt> newifids( insize, 0 ) ;
    3340     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3341       //os << "Table " << itable << ": " ;
    3342       for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
    3343         addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
    3344         //os << addrow[itable][ifrow] << " " ;
    3345       }
    3346       addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
    3347       //os << "(" << addtable[itable] << ")" << LogIO::POST ;
    3348     }
    3349     newin.resize( insize ) ;
    3350     copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
     3259    oss.str("") ;
    33513260    for ( uInt i = 0 ; i < insize ; i++ ) {
    3352       newtableids[i] = i ;
    3353     }
    3354     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3355       for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
    3356         CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
    3357         vector<int> freqidlist ;
    3358         for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
    3359           if ( groups[itable][i].size() > iadd + 1 ) {
    3360             freqidlist.push_back( freqid[itable][i] ) ;
    3361           }
    3362         }
    3363         stringstream taqlstream ;
    3364         taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
    3365         for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
    3366           taqlstream << freqidlist[i] ;
    3367           if ( i < freqidlist.size() - 1 )
    3368             taqlstream << "," ;
    3369           else
    3370             taqlstream << "]" ;
    3371         }
    3372         string taql = taqlstream.str() ;
    3373         //os << "taql = " << taql << LogIO::POST ;
    3374         STSelector selector = STSelector() ;
    3375         selector.setTaQL( taql ) ;
    3376         add->setSelection( selector ) ;
    3377         newin.push_back( add ) ;
    3378         newtableids.push_back( itable ) ;
    3379         newifids.push_back( iadd + 1 ) ;
    3380       }
    3381     }
    3382 
    3383     // udpate ifgrp
    3384     uInt offset = insize ;
    3385     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3386       for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
    3387         for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
    3388           if ( groups[itable][ifrow].size() > iadd + 1 ) {
    3389             uInt igrp = groups[itable][ifrow][iadd+1] ;
    3390             for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
    3391               if ( ifgrp[igrp][2*imem] == newtableids[offset+iadd] && ifgrp[igrp][2*imem+1] == freqid[newtableids[offset+iadd]][ifrow] ) {
    3392                 ifgrp[igrp][2*imem] = offset + iadd ;
    3393               }
    3394             }
    3395           }
    3396         }
    3397       }
    3398       offset += addtable[itable] ;
    3399     }
    3400 
    3401     // print IF groups again for debug
    3402 //     oss.str( "" ) ;
    3403 //     oss << "IF Group summary: " << endl ;
    3404 //     oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
    3405 //     for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
    3406 //       oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
    3407 //       for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
    3408 //         oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
    3409 //       }
    3410 //       oss << endl ;
    3411 //     }
    3412 //     oss << endl ;
    3413 //     os << oss.str() << LogIO::POST ;
     3261      oss << "Table " << i << endl ;
     3262      for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
     3263        oss << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
     3264        oss << setw( 2 ) << groups[i][j] ;
     3265        oss << endl ;
     3266      }
     3267    }
     3268    os << oss.str() << LogIO::POST ;
    34143269
    34153270    // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
    34163271    //os << "All IF number is set to IF group index" << LogIO::POST ;
    3417     insize = newin.size() ;
    34183272    // reset SCANNO only when avmode != "SCAN"
    34193273    if ( avmode != "SCAN" ) {
    … …  
    34213275      for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    34223276        uInt nrow = newin[itable]->nrow() ;
    3423         Table &tmpt = newin[itable]->table() ;
    3424         scannoCol.attach( tmpt, "SCANNO" ) ;
    3425         for ( uInt irow = 0 ; irow < nrow ; irow++ )
    3426           scannoCol.put( irow, 0 ) ;
    3427       }
    3428     }
    3429     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3430       uInt rows = newin[itable]->nrow() ;
    3431       Table &tmpt = newin[itable]->table() ;
    3432       freqIDCol.attach( tmpt, "FREQ_ID" ) ;
    3433       ifnoCol.attach( tmpt, "IFNO" ) ;
    3434       for ( uInt irow=0 ; irow < rows ; irow++ ) {
    3435         uInt freqID = freqIDCol( irow ) ;
    3436         vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
    3437         if ( iter != freqid[newtableids[itable]].end() ) {
    3438           uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
    3439           ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
    3440         }
    3441         else {
    3442           throw(AipsError("IF grouping was wrong in additional tables.")) ;
    3443         }
    3444       }
    3445     }
    3446     oldinfo.resize( insize ) ;
    3447     setInsitu( oldInsitu ) ;
    3448 
    3449     // temporarily set coordinfo
    3450     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3451       vector<string> coordinfo = newin[itable]->getCoordInfo() ;
    3452       oldinfo[itable] = coordinfo[0] ;
    3453       coordinfo[0] = "Hz" ;
    3454       newin[itable]->setCoordInfo( coordinfo ) ;
    3455     }
    3456 
    3457     // save column values in the vector
    3458     vector< vector<uInt> > freqIdVec( insize ) ;
    3459     vector< vector<uInt> > ifNoVec( insize ) ;
    3460     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3461       ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
    3462       freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
    3463       for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
    3464         freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
    3465         ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
    3466       }
    3467     }
    3468 
    3469     // reset spectra and flagtra: pick up common part of frequency coverage
    3470     //os << "Pick common frequency range and align resolution" << LogIO::POST ;
    3471     for ( uInt itable = 0 ; itable < insize ; itable++ ) {
    3472       uInt rows = newin[itable]->nrow() ;
    3473       int nminchan = -1 ;
    3474       int nmaxchan = -1 ;
    3475       vector<uInt> freqIdUpdate ;
    3476       for ( uInt irow = 0 ; irow < rows ; irow++ ) {
    3477         uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index (IF group index)
    3478         double minfreq = ifgfreq[2*ifno] ;
    3479         double maxfreq = ifgfreq[2*ifno+1] ;
    3480         //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ;
    3481         vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
    3482         int nchan = abcissa.size() ;
    3483         double resol = abcissa[1] - abcissa[0] ;
    3484         //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
    3485         uInt imin, imax;
    3486         int sigres;
    3487         if ( resol >= 0. ) {
    3488           imin = 0;
    3489           imax = nchan - 1 ;
    3490           sigres = 1 ;
    3491         } else {
    3492           imin = nchan - 1 ;
    3493           imax = 0 ;
    3494           sigres = -1 ;
    3495           resol = abs(resol) ;
    3496         }
    3497         if ( minfreq <= abcissa[imin] )
    3498           nminchan = imin ;
    3499         else {
    3500           //double cfreq = ( minfreq - abcissa[0] ) / resol ;
    3501           double cfreq = ( minfreq - abcissa[imin] + 0.5 * resol ) / resol ;
    3502           nminchan = imin + sigres * (int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 )) ;
    3503         }
    3504         if ( maxfreq >= abcissa[imax] )
    3505           nmaxchan = imax;
    3506         else {
    3507           //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
    3508           double cfreq = ( abcissa[imax] - maxfreq + 0.5 * resol ) / resol ;
    3509           nmaxchan = imax - sigres * (int(cfreq) +( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 )) ;
    3510         }
    3511         //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
    3512         if ( nmaxchan < nminchan ) {
    3513           int tmp = nmaxchan ;
    3514           nmaxchan = nminchan ;
    3515           nminchan = tmp ;
    3516         }
    3517         if ( nmaxchan > nminchan ) {
    3518           newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
    3519           int newchan = nmaxchan - nminchan + 1 ;
    3520           if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan ) {
    3521             freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
    3522 
    3523             // Update before nminchan is lost
    3524             uInt freqId = freqIdVec[itable][irow] ;
    3525             Double refpix ;
    3526             Double refval ;
    3527             Double increment ;
    3528             newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
    3529             //refval = refval - ( refpix - nminchan ) * increment ;
    3530             refval = abcissa[nminchan] ;
    3531             refpix = 0 ;
    3532             newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
    3533           }
    3534         }
    3535         else {
    3536           throw(AipsError("Failed to pick up common part of frequency range.")) ;
    3537         }
    3538       }
    3539 //       for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
    3540 //      uInt freqId = freqIdUpdate[i] ;
    3541 //      Double refpix ;
    3542 //      Double refval ;
    3543 //      Double increment ;
    3544        
    3545 //      // update row
    3546 //      newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
    3547 //      refval = refval - ( refpix - nminchan ) * increment ;
    3548 //      refpix = 0 ;
    3549 //      newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
    3550 //       }   
    3551     }
    3552 
    3553    
     3277        Vector<uInt> resetScan( nrow, 0 ) ;
     3278        newin[itable]->scanCol_.putColumn( resetScan ) ;
     3279      }
     3280    }
     3281
    35543282    // reset spectra and flagtra: align spectral resolution
    35553283    //os << "Align spectral resolution" << LogIO::POST ;
    35563284    // gmaxdnu: the coarsest frequency resolution in the frequency group
    3557     // gmemid: member index that have a resolution equal to gmaxdnu
    3558     // gmaxdnu[numfreqgrp]
    3559     // gmaxdnu: [dnu0, dnu1, ...]
    3560     // gmemid[numfreqgrp]
    3561     // gmemid: [id0, id1, ...]
     3285    // gminfreq: lower frequency edge of the frequency group
     3286    // gnchan: number of channels for the frequency group
    35623287    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
    3563     vector<uInt> gmemid( freqgrp.size(), 0 ) ;
    3564     for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
    3565       double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution
    3566       int minchan = INT_MAX ;     // minimum channel number
    3567       Double refpixref = -1 ;     // reference of 'reference pixel'
    3568       Double refvalref = -1 ;     // reference of 'reference frequency'
    3569       Double refinc = -1 ;        // reference frequency resolution
    3570       uInt refreqid ;
    3571       uInt reftable = INT_MAX;
    3572       // process only if group member > 1
    3573       if ( ifgrp[igrp].size() > 2 ) {
    3574         // find minchan and maxdnu in each group
    3575         for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
    3576           uInt tableid = ifgrp[igrp][2*imem] ;
    3577           uInt rowid = ifgrp[igrp][2*imem+1] ;
    3578           vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
    3579           if ( iter != freqIdVec[tableid].end() ) {
    3580             uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
    3581             vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
    3582             int nchan = abcissa.size() ;
    3583             double dnu = abcissa[1] - abcissa[0] ;
    3584             //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ;
    3585             if ( nchan < minchan ) {
    3586               minchan = nchan ;
    3587               maxdnu = abs(dnu) ;
    3588               newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
    3589               refreqid = rowid ;
    3590               reftable = tableid ;
    3591               // Spectra are reversed when dnu < 0
    3592               if (dnu < 0)
    3593                 refpixref = minchan - 1 -refpixref ;
    3594             }
    3595           }
    3596         }
    3597         // regrid spectra in each group
    3598         os << "GROUP " << igrp << endl ;
    3599         os << "   Channel number is adjusted to " << minchan << endl ;
    3600         os << "   Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ;
    3601         for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
    3602           uInt tableid = ifgrp[igrp][2*imem] ;
    3603           uInt rowid = ifgrp[igrp][2*imem+1] ;
    3604           freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
    3605           //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ;
    3606           //os << "   regridChannel applied to " ;
    3607           //if ( tableid != reftable )
    3608           refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, maxdnu ) ;
    3609           for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
    3610             uInt tfreqid = freqIdVec[tableid][irow] ;
    3611             if ( tfreqid == rowid ) {     
    3612               //os << irow << " " ;
    3613               newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
    3614               freqIDCol.put( irow, refreqid ) ;
    3615               freqIdVec[tableid][irow] = refreqid ;
    3616             }
    3617           }
    3618           //os << LogIO::POST ;
    3619         }
    3620       }
    3621       else {
    3622         uInt tableid = ifgrp[igrp][0] ;
    3623         uInt rowid = ifgrp[igrp][1] ;
    3624         vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
    3625         if ( iter != freqIdVec[tableid].end() ) {
    3626           uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
    3627           vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
    3628           minchan = abcissa.size() ;
    3629           maxdnu = abs( abcissa[1] - abcissa[0]) ;
    3630         }
    3631       }
    3632       for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
    3633         if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
    3634           if ( maxdnu > gmaxdnu[i] ) {
    3635             gmaxdnu[i] = maxdnu ;
    3636             gmemid[i] = igrp ;
    3637           }
    3638           break ;
    3639         }
     3288    vector<double> gminfreq( freqgrp.size() ) ;
     3289    vector<double> gnchan( freqgrp.size() ) ;
     3290    for ( uInt i = 0 ; i < insize ; i++ ) {
     3291      vector<uInt> members = groups[i] ;
     3292      for ( uInt j = 0 ; j < members.size() ; j++ ) {
     3293        uInt groupid = members[j] ;
     3294        Double rp,rv,ic ;
     3295        newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
     3296        if ( abs(ic) > abs(gmaxdnu[groupid]) )
     3297          gmaxdnu[groupid] = ic ;
     3298      }
     3299    }
     3300    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     3301      gminfreq[igrp] = grprange[2*igrp] ;
     3302      double maxfreq = grprange[2*igrp+1] ;
     3303      gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
     3304    }
     3305     
     3306    // regrid spectral data and update frequency info
     3307    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
     3308      Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
     3309      Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
     3310
     3311      // update MAIN
     3312      for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
     3313        uInt groupid = groups[itable][oldFreqId[irow]] ;
     3314        newFreqId[irow] = groupid ;
     3315        newin[itable]->regridChannel( gnchan[groupid],
     3316                                      gmaxdnu[groupid],
     3317                                      gminfreq[groupid],
     3318                                      irow ) ;
     3319      }
     3320      newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
     3321      newin[itable]->ifCol_.putColumn( newFreqId ) ;
     3322
     3323      // update FREQUENCIES
     3324      Table tab = newin[itable]->frequencies().table() ;
     3325      ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
     3326      ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
     3327      ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
     3328      ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
     3329      if ( freqgrp.size() > tab.nrow() ) {
     3330        tab.addRow( freqgrp.size()-tab.nrow() ) ;
     3331      }
     3332      for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
     3333        Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
     3334        Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
     3335        Double increment = gmaxdnu[irow] ;
     3336        fIdCol.put( irow, irow ) ;
     3337        fRefPixCol.put( irow, refpix ) ;
     3338        fRefValCol.put( irow, refval ) ;
     3339        fIncrCol.put( irow, increment ) ;
    36403340      }
    36413341    }
    … …  
    36463346      coordinfo[0] = oldinfo[itable] ;
    36473347      newin[itable]->setCoordInfo( coordinfo ) ;
    3648     }     
    3649 
    3650     // accumulate all rows into the first table
    3651     // NOTE: assumed in.size() = 1
    3652     vector< CountedPtr<Scantable> > tmp( 1 ) ;
    3653     if ( newin.size() == 1 )
    3654       tmp[0] = newin[0] ;
    3655     else
    3656       tmp[0] = merge( newin ) ;
    3657 
    3658     //return tmp[0] ;
     3348    }
    36593349
    36603350    // average
    3661     CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
    3662 
    3663     //return tmpout ;
    3664 
    3665     // combine frequency group
    3666     os << "Combine spectra based on frequency grouping" << LogIO::POST ;
    3667     os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
    3668     vector<string> coordinfo = tmpout->getCoordInfo() ;
    3669     oldinfo[0] = coordinfo[0] ;
    3670     coordinfo[0] = "Hz" ;
    3671     tmpout->setCoordInfo( coordinfo ) ;
    3672     // create proformas of output table
    3673     stringstream taqlstream ;
    3674     taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
    3675     for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
    3676       taqlstream << gmemid[i] ;
    3677       if ( i < gmemid.size() - 1 )
    3678         taqlstream << "," ;
    3679       else
    3680         taqlstream << "]" ;
    3681     }
    3682     string taql = taqlstream.str() ;
    3683     //os << "taql = " << taql << LogIO::POST ;
    3684     STSelector selector = STSelector() ;
    3685     selector.setTaQL( taql ) ;
    3686     oldInsitu = insitu_ ;
    3687     setInsitu( false ) ;
    3688     out = getScantable( tmpout, false ) ;
    3689     setInsitu( oldInsitu ) ;
    3690     out->setSelection( selector ) ;
    3691     // regrid rows
    3692     ifnoCol.attach( tmpout->table(), "IFNO" ) ;
    3693     for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
    3694       uInt ifno = ifnoCol( irow ) ;
    3695       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
    3696         if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
    3697           vector<double> abcissa = tmpout->getAbcissa( irow ) ;
    3698           double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
    3699           int nchan = (int) abs( bw / gmaxdnu[igrp] ) ;
    3700           // All spectra will have positive frequency increments
    3701           tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
    3702           break ;
    3703         }
    3704       }
    3705     }
    3706     // combine spectra
    3707     ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
    3708     ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
    3709     ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
    3710     ScalarColumn<uInt> ifnoColOut( out->table(), "IFNO" ) ;
    3711     ScalarColumn<uInt> polnoColOut( out->table(), "POLNO" ) ;
    3712     ScalarColumn<uInt> freqidColOut( out->table(), "FREQ_ID" ) ;
    3713 //     MDirection::ScalarColumn dirColOut ;
    3714 //     dirColOut.attach( out->table(), "DIRECTION" ) ;
    3715     Table &tab = tmpout->table() ;
    3716 //     Block<String> cols(1);
    3717 //     cols[0] = String("POLNO") ;
    3718     Block<String> cols(3);
    3719     cols[0] = String("POLNO") ;
    3720     cols[1] = String("SCANNO") ;
    3721     cols[2] = String("BEAMNO") ;
    3722     TableIterator iter( tab, cols ) ;
    3723     vector< vector<uInt> > sizes( freqgrp.size() ) ;
    3724     vector<uInt> totalsizes( freqgrp.size(), 0 ) ;
    3725     ArrayColumn<Float> specCols ;
    3726     ArrayColumn<uChar> flagCols ;
    3727     ArrayColumn<Float> tsysCols ;
    3728     ScalarColumn<uInt> polnos ;
    3729     vector< Vector<Float> > specout( freqgrp.size() ) ;
    3730     vector< Vector<uChar> > flagout( freqgrp.size() ) ;
    3731     vector< Vector<Float> > tsysout( freqgrp.size() ) ;
    3732     while( !iter.pastEnd() ) {
    3733       specCols.attach( iter.table(), "SPECTRA" ) ;
    3734       flagCols.attach( iter.table(), "FLAGTRA" ) ;
    3735       tsysCols.attach( iter.table(), "TSYS" ) ;
    3736       ifnoCol.attach( iter.table(), "IFNO" ) ;
    3737       polnos.attach( iter.table(), "POLNO" ) ;
    3738 //       MDirection::ScalarColumn dircol ;
    3739 //       dircol.attach( iter.table(), "DIRECTION" ) ;
    3740       uInt polno = polnos( 0 ) ;
    3741       //os << "POLNO iteration: " << polno << LogIO::POST ;
    3742 //       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
    3743 //      sizes[igrp].resize( freqgrp[igrp].size() ) ;
    3744 //      for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
    3745 //        for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
    3746 //          uInt ifno = ifnoCol( irow ) ;
    3747 //          if ( ifno == freqgrp[igrp][imem] ) {
    3748 //            Vector<Float> spec = specCols( irow ) ;
    3749 //            Vector<uChar> flag = flagCols( irow ) ;
    3750 //            vector<Float> svec ;
    3751 //            spec.tovector( svec ) ;
    3752 //            vector<uChar> fvec ;
    3753 //            flag.tovector( fvec ) ;
    3754 //            //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
    3755 //            specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
    3756 //            flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
    3757 //            //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
    3758 //            sizes[igrp][imem] = spec.nelements() ;
    3759 //          }
    3760 //        }
    3761 //      }
    3762 //      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
    3763 //        uInt ifout = ifnoColOut( irow ) ;
    3764 //        uInt polout = polnoColOut( irow ) ;
    3765 //        if ( ifout == gmemid[igrp] && polout == polno ) {
    3766 //          // set SPECTRA and FRAGTRA
    3767 //          Vector<Float> newspec( specout[igrp] ) ;
    3768 //          Vector<uChar> newflag( flagout[igrp] ) ;
    3769 //          specColOut.put( irow, newspec ) ;
    3770 //          flagColOut.put( irow, newflag ) ;
    3771 //          // IFNO renumbering
    3772 //          ifnoColOut.put( irow, igrp ) ;
    3773 //        }
    3774 //      }
    3775 //       }
    3776       // get a list of number of channels for each frequency group member
    3777       if ( totalsizes[0] == 0 ) {
    3778         for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
    3779           sizes[igrp].resize( freqgrp[igrp].size() ) ;
    3780           for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
    3781             for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
    3782               uInt ifno = ifnoCol( irow ) ;
    3783               if ( ifno == freqgrp[igrp][imem] ) {
    3784                 Vector<Float> spec = specCols( irow ) ;
    3785                 sizes[igrp][imem] = spec.nelements() ;
    3786                 totalsizes[igrp] += sizes[igrp][imem] ;
    3787                 break ;
    3788               }
    3789             }
    3790           }
    3791           specout[igrp].resize( totalsizes[igrp] ) ;
    3792           flagout[igrp].resize( totalsizes[igrp] ) ;
    3793           tsysout[igrp].resize( totalsizes[igrp] ) ;
    3794         }
    3795       }
    3796 
    3797       // combine spectra
    3798       for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
    3799         uInt polout = polnoColOut( irow ) ;
    3800         if ( polout == polno ) {
    3801           uInt ifout = ifnoColOut( irow ) ;
    3802 //           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
    3803           uInt igrp ;
    3804           for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
    3805             if ( ifout == gmemid[jgrp] ) {
    3806               igrp = jgrp ;
    3807               break ;
    3808             }
    3809           }
    3810           IPosition startpos( 1, 0 ) ;
    3811           IPosition endpos( 1, 0 ) ;
    3812           for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
    3813             for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
    3814               uInt ifno = ifnoCol( jrow ) ;
    3815               // 2012/02/17 TN
    3816               // Since STGrid is implemented, average doesn't consider direction
    3817               // when accumulating
    3818 //               Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
    3819 //               //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
    3820 //               Double dx = tdir[0] - direction[0] ;
    3821 //               Double dy = tdir[1] - direction[1] ;
    3822 //               Double dd = sqrt( dx * dx + dy * dy ) ;
    3823               //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
    3824 //               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
    3825               if ( ifno == freqgrp[igrp][imem] ) {
    3826                 endpos[0] = startpos[0] + sizes[igrp][imem] - 1 ;
    3827                 Vector<Float> spec = specCols( jrow ) ;
    3828                 Vector<uChar> flag = flagCols( jrow ) ;
    3829                 Vector<Float> tsys = tsysCols( jrow ) ;
    3830                 //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
    3831                 specout[igrp]( startpos, endpos ) = spec ;
    3832                 flagout[igrp]( startpos, endpos ) = flag ;
    3833                 if ( spec.size() == tsys.size() ) {
    3834                   tsysout[igrp]( startpos, endpos ) = tsys ;
    3835                 }
    3836                 else {
    3837                   tsysout[igrp]( startpos, endpos ) = tsys[0] ;
    3838                 }
    3839                 //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
    3840                 startpos[0] += sizes[igrp][imem] ;
    3841               }
    3842             }
    3843           }
    3844           // set SPECTRA and FRAGTRA
    3845           specColOut.put( irow, specout[igrp] ) ;
    3846           flagColOut.put( irow, flagout[igrp] ) ;
    3847           tsysColOut.put( irow, tsysout[igrp] ) ;
    3848           // IFNO renumbering (renumbered as frequency group ID)
    3849           ifnoColOut.put( irow, igrp ) ;
    3850         }
    3851       }
    3852       iter++ ;
    3853     }
    3854     // update FREQUENCIES subtable
    3855     vector<bool> updated( freqgrp.size(), false ) ;
    3856     for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
    3857       uInt index = 0 ;
    3858       uInt pixShift = 0 ;
    3859 
    3860       while ( freqgrp[igrp][index] != gmemid[igrp] ) {
    3861         pixShift += sizes[igrp][index++] ;
    3862       }
    3863       for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
    3864         //if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
    3865         if ( ifnoColOut( irow ) == igrp && !updated[igrp] ) {
    3866           uInt freqidOut = freqidColOut( irow ) ;
    3867           //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
    3868           double refpix ;
    3869           double refval ;
    3870           double increm ;
    3871           out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
    3872           if (increm < 0)
    3873             refpix = sizes[igrp][index] - 1 - refpix ; // reversed
    3874           refpix += pixShift ;
    3875           out->frequencies().setEntry( refpix, refval, gmaxdnu[igrp], freqidOut ) ;
    3876           updated[igrp] = true ;
    3877         }
    3878       }
    3879     }
    3880 
    3881     //out = tmpout ;
    3882 
    3883     coordinfo = tmpout->getCoordInfo() ;
    3884     coordinfo[0] = oldinfo[0] ;
    3885     tmpout->setCoordInfo( coordinfo ) ;
     3351    out = average( newin, mask, weight, avmode ) ;
    38863352  }
    38873353  else {

  • trunk/src/Scantable.cpp

    r2591 r2595  
    23042304}
    23052305
     2306void Scantable::regridChannel( int nChan, double dnu, double fmin, int irow )
     2307{
     2308  Vector<Float> oldspec = specCol_( irow ) ;
     2309  Vector<uChar> oldflag = flagsCol_( irow ) ;
     2310  Vector<Float> oldtsys = tsysCol_( irow ) ;
     2311  Vector<Float> newspec( nChan, 0 ) ;
     2312  Vector<uChar> newflag( nChan, true ) ;
     2313  Vector<Float> newtsys ;
     2314  bool regridTsys = false ;
     2315  if (oldtsys.size() == oldspec.size()) {
     2316    regridTsys = true ;
     2317    newtsys.resize(nChan,false) ;
     2318    newtsys = 0 ;
     2319  }
     2320 
     2321  // regrid
     2322  vector<double> abcissa = getAbcissa( irow ) ;
     2323  int oldsize = abcissa.size() ;
     2324  double olddnu = abcissa[1] - abcissa[0] ;
     2325  //int ichan = 0 ;
     2326  double wsum = 0.0 ;
     2327  Vector<double> zi( nChan+1 ) ;
     2328  Vector<double> yi( oldsize + 1 ) ;
     2329  Block<uInt> count( nChan, 0 ) ;
     2330  yi[0] = abcissa[0] - 0.5 * olddnu ;
     2331  for ( int ii = 1 ; ii < oldsize ; ii++ )
     2332    yi[ii] = 0.5* (abcissa[ii-1] + abcissa[ii]) ;
     2333  yi[oldsize] = abcissa[oldsize-1] \
     2334    + 0.5 * (abcissa[oldsize-1] - abcissa[oldsize-2]) ;
     2335//   cout << "olddnu=" << olddnu << ", dnu=" << dnu << " (diff=" << olddnu-dnu << ")" << endl ;
     2336//   cout << "yi[0]=" << yi[0] << ", fmin=" << fmin << " (diff=" << yi[0]-fmin << ")" << endl ;
     2337//   cout << "oldsize=" << oldsize << ", nChan=" << nChan << endl ;
     2338
     2339  // do not regrid if input parameters are almost same as current
     2340  // spectral setup
     2341  double dnuDiff = abs( ( dnu - olddnu ) / olddnu ) ;
     2342  double oldfmin = min( yi[0], yi[oldsize] ) ;
     2343  double fminDiff = abs( ( fmin - oldfmin ) / oldfmin ) ;
     2344  double nChanDiff = nChan - oldsize ;
     2345  double eps = 1.0e-8 ;
     2346  if ( nChanDiff == 0 && dnuDiff < eps && fminDiff < eps )
     2347    return ;
     2348
     2349  //zi[0] = abcissa[0] - 0.5 * olddnu ;
     2350  //zi[0] = ((olddnu*dnu > 0) ? yi[0] : yi[oldsize]) ;
     2351  if ( dnu > 0 )
     2352    zi[0] = fmin - 0.5 * dnu ;
     2353  else
     2354    zi[0] = fmin + nChan * abs(dnu) ;
     2355  for ( int ii = 1 ; ii < nChan ; ii++ )
     2356    zi[ii] = zi[0] + dnu * ii ;
     2357  zi[nChan] = zi[nChan-1] + dnu ;
     2358  // Access zi and yi in ascending order
     2359  int izs = ((dnu > 0) ? 0 : nChan ) ;
     2360  int ize = ((dnu > 0) ? nChan : 0 ) ;
     2361  int izincr = ((dnu > 0) ? 1 : -1 ) ;
     2362  int ichan =  ((olddnu > 0) ? 0 : oldsize ) ;
     2363  int iye = ((olddnu > 0) ? oldsize : 0 ) ;
     2364  int iyincr = ((olddnu > 0) ? 1 : -1 ) ;
     2365  //for ( int ii = izs ; ii != ize ; ii+=izincr ){
     2366  int ii = izs ;
     2367  while (ii != ize) {
     2368    // always zl < zr
     2369    double zl = zi[ii] ;
     2370    double zr = zi[ii+izincr] ;
     2371    // Need to access smaller index for the new spec, flag, and tsys.
     2372    // Values between zi[k] and zi[k+1] should be stored in newspec[k], etc.
     2373    int i = min(ii, ii+izincr) ;
     2374    //for ( int jj = ichan ; jj != iye ; jj+=iyincr ) {
     2375    int jj = ichan ;
     2376    while (jj != iye) {
     2377      // always yl < yr
     2378      double yl = yi[jj] ;
     2379      double yr = yi[jj+iyincr] ;
     2380      // Need to access smaller index for the original spec, flag, and tsys.
     2381      // Values between yi[k] and yi[k+1] are stored in oldspec[k], etc.
     2382      int j = min(jj, jj+iyincr) ;
     2383      if ( yr <= zl ) {
     2384        jj += iyincr ;
     2385        continue ;
     2386      }
     2387      else if ( yl <= zl ) {
     2388        if ( yr < zr ) {
     2389          if (!oldflag[j]) {
     2390            newspec[i] += oldspec[j] * ( yr - zl ) ;
     2391            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - zl ) ;
     2392            wsum += ( yr - zl ) ;
     2393            count[i]++ ;
     2394          }
     2395          newflag[i] = newflag[i] && oldflag[j] ;
     2396        }
     2397        else {
     2398          if (!oldflag[j]) {
     2399            newspec[i] += oldspec[j] * abs(dnu) ;
     2400            if (regridTsys) newtsys[i] += oldtsys[j] * abs(dnu) ;
     2401            wsum += abs(dnu) ;
     2402            count[i]++ ;
     2403          }
     2404          newflag[i] = newflag[i] && oldflag[j] ;
     2405          ichan = jj ;
     2406          break ;
     2407        }
     2408      }
     2409      else if ( yl < zr ) {
     2410        if ( yr <= zr ) {
     2411          if (!oldflag[j]) {
     2412            newspec[i] += oldspec[j] * ( yr - yl ) ;
     2413            if (regridTsys) newtsys[i] += oldtsys[j] * ( yr - yl ) ;
     2414            wsum += ( yr - yl ) ;
     2415            count[i]++ ;
     2416          }
     2417          newflag[i] = newflag[i] && oldflag[j] ;
     2418        }
     2419        else {
     2420          if (!oldflag[j]) {
     2421            newspec[i] += oldspec[j] * ( zr - yl ) ;
     2422            if (regridTsys) newtsys[i] += oldtsys[j] * ( zr - yl ) ;
     2423            wsum += ( zr - yl ) ;
     2424            count[i]++ ;
     2425          }
     2426          newflag[i] = newflag[i] && oldflag[j] ;
     2427          ichan = jj ;
     2428          break ;
     2429        }
     2430      }
     2431      else {
     2432        //ichan = jj - iyincr ;
     2433        break ;
     2434      }
     2435      jj += iyincr ;
     2436    }
     2437    if ( wsum != 0.0 ) {
     2438      newspec[i] /= wsum ;
     2439      if (regridTsys) newtsys[i] /= wsum ;
     2440    }
     2441    wsum = 0.0 ;
     2442    ii += izincr ;
     2443  }
     2444
     2445  // flag out channels without data
     2446  // this is tentative since there is no specific definition
     2447  // on bit flag...
     2448  uChar noData = 1 << 7 ;
     2449  for ( Int i = 0 ; i < nChan ; i++ ) {
     2450    if ( count[i] == 0 )
     2451      newflag[i] = noData ;
     2452  }
     2453
     2454  specCol_.put( irow, newspec ) ;
     2455  flagsCol_.put( irow, newflag ) ;
     2456  if (regridTsys) tsysCol_.put( irow, newtsys );
     2457
     2458  return ;
     2459}
     2460
    23062461std::vector<float> Scantable::getWeather(int whichrow) const
    23072462{

  • trunk/src/Scantable.h

    r2591 r2595  
    501501  void regridChannel( int nchan, double dnu ) ;
    502502  void regridChannel( int nchan, double dnu, int irow ) ;
     503  void regridChannel( int nchan, double dnu, double fmin, int irow ) ;
    503504
    504505  void regridSpecChannel( double dnu, int nchan=-1 ) ;
Note: See TracChangeset for help on using the changeset viewer.