Changeset 1607 for branches/alma/src


Ignore:
Timestamp:
07/21/09 13:05:36 (15 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-1423

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: run sdaverage with OTF data

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Previous averaging methods could not deal with OTF data since these methods
does not refer DIRECTION information so that dumped spectra with different
DIRECTION are unwillingly accumulated.
Currently, these methods refer DIRECTION column to support OTF data.
Therefore, spectra are accumulated only when DIRECTION is exactly same.
I have defined a tolerance for the determination of the coincidence of
DIRECTION, but it is very small value (order of significant digits for
double) right now.


Location:
branches/alma/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/alma/src/STMath.cpp

    r1603 r1607  
    5454using namespace asap;
    5555
     56// tolerance for direction comparison (rad)
     57Double tol = 1.0e-15 ;
     58//Double tol = 1.0 ;
     59
    5660STMath::STMath(bool insitu) :
    5761  insitu_(insitu)
    … …  
    7478                    "Use merge first."));
    7579  WeightType wtype = stringToWeight(weight);
     80
     81  // check if OTF observation
     82  String obstype = in[0]->getHeader().obstype ;
     83  bool otfscan = false ;
     84  if ( obstype.find( "OTF" ) != String::npos ) {
     85    //cout << "OTF scan" << endl ;
     86    otfscan = true ;
     87  }
    7688
    7789  // output
    … …  
    121133  uInt outrowCount = 0;
    122134  TableIterator iter(baset, cols);
     135  int count = 0 ;
    123136  while (!iter.pastEnd()) {
    124137    Table subt = iter.table();
    125     // copy the first row of this selection into the new table
    126     tout.addRow();
    127     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    128     // re-index to 0
    129     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
    130       scanColOut.put(outrowCount, uInt(0));
    131     }
    132     ++outrowCount;
     138//     // copy the first row of this selection into the new table
     139//     tout.addRow();
     140//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     141//     // re-index to 0
     142//     if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     143//       scanColOut.put(outrowCount, uInt(0));
     144//     }
     145//     ++outrowCount;
     146    if ( otfscan ) {
     147      MDirection::ScalarColumn dircol ;
     148      dircol.attach( subt, "DIRECTION" ) ;
     149      Int length = subt.nrow() ;
     150      vector< Vector<Double> > dirs ;
     151      vector<int> indexes ;
     152      for ( Int i = 0 ; i < length ; i++ ) {
     153        Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     154        //cout << setw( 3 ) << count++ << ": " ;
     155        //cout << "[" << t[0] << "," << t[1] << "]" <<  endl ;
     156        bool adddir = true ;
     157        for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     158          //if ( allTrue( t == dirs[j] ) ) {
     159          if ( allNearAbs( t, dirs[j], tol ) ) {
     160            adddir = false ;
     161            break ;
     162          }
     163        }
     164        if ( adddir ) {
     165          dirs.push_back( t ) ;
     166          indexes.push_back( i ) ;
     167        }
     168      }
     169      uInt rowNum = dirs.size() ;
     170      //cout << "dirs.size()=" << dirs.size() << endl ;
     171      tout.addRow( rowNum ) ;
     172      for ( uInt i = 0 ; i < rowNum ; i++ ) {
     173        TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
     174        // re-index to 0
     175        if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     176          scanColOut.put(outrowCount+i, uInt(0));
     177        }       
     178      }
     179      outrowCount += rowNum ;
     180    }
     181    else {
     182      // copy the first row of this selection into the new table
     183      tout.addRow();
     184      TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     185      // re-index to 0
     186      if ( avmode != "SCAN" && avmode != "SOURCE" ) {
     187        scanColOut.put(outrowCount, uInt(0));
     188      }
     189      ++outrowCount;
     190    }
    133191    ++iter;
    134192  }
    … …  
    162220      } else {
    163221        subt = basesubt;
     222      }
     223      // for OTF
     224      if ( otfscan ) {
     225        vector<uInt> removeRows ;
     226        uInt nrsubt = subt.nrow() ;
     227        for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
     228          //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
     229          if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
     230            removeRows.push_back( irow ) ;
     231          }
     232        }
     233        //cout << "removeRows.size()=" << removeRows.size() << endl ;
     234        if ( removeRows.size() != 0 ) {
     235          //cout << "[" ;
     236          //for ( uInt irow=0 ; irow<removeRows.size()-1 ; irow++ )
     237          //cout << removeRows[irow] << "," ;
     238          //cout << removeRows[removeRows.size()-1] << "]" << endl ;
     239          subt.removeRow( removeRows ) ;
     240        }
     241       
     242        if ( nrsubt == removeRows.size() )
     243          throw(AipsError("Averaging data is empty.")) ;
    164244      }
    165245      specCol.attach(subt,"SPECTRA");
    … …  
    230310                          const std::string& avmode )
    231311{
     312  // check if OTF observation
     313  String obstype = in->getHeader().obstype ;
     314  bool otfscan = false ;
     315  if ( obstype.find( "OTF" ) != String::npos ) {
     316    //cout << "OTF scan" << endl ;
     317    otfscan = true ;
     318  }
     319
    232320  // clone as this is non insitu
    233321  bool insitu = insitu_;
    … …  
    261349    flagCol.attach(subt,"FLAGTRA");
    262350    tsysCol.attach(subt,"TSYS");
    263     tout.addRow();
    264     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
    265     if ( avmode != "SCAN") {
    266       scanColOut.put(outrowCount, uInt(0));
    267     }
    268     Vector<Float> tmp;
    269     specCol.get(0, tmp);
    270     uInt nchan = tmp.nelements();
    271     // have to do channel by channel here as MaskedArrMath
    272     // doesn't have partialMedians
    273     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
    274     Vector<Float> outspec(nchan);
    275     Vector<uChar> outflag(nchan,0);
    276     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
    277     for (uInt i=0; i<nchan; ++i) {
    278       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
    279       MaskedArray<Float> ma = maskedArray(specs,flags);
    280       outspec[i] = median(ma);
    281       if ( allEQ(ma.getMask(), False) )
    282         outflag[i] = userflag;// flag data
    283     }
    284     outtsys[0] = median(tsysCol.getColumn());
    285     specColOut.put(outrowCount, outspec);
    286     flagColOut.put(outrowCount, outflag);
    287     tsysColOut.put(outrowCount, outtsys);
    288     Double intsum = sum(intCol.getColumn());
    289     intColOut.put(outrowCount, intsum);
    290     ++outrowCount;
     351//     tout.addRow();
     352//     TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     353//     if ( avmode != "SCAN") {
     354//       scanColOut.put(outrowCount, uInt(0));
     355//     }
     356//     Vector<Float> tmp;
     357//     specCol.get(0, tmp);
     358//     uInt nchan = tmp.nelements();
     359//     // have to do channel by channel here as MaskedArrMath
     360//     // doesn't have partialMedians
     361//     Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     362//     Vector<Float> outspec(nchan);
     363//     Vector<uChar> outflag(nchan,0);
     364//     Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     365//     for (uInt i=0; i<nchan; ++i) {
     366//       Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     367//       MaskedArray<Float> ma = maskedArray(specs,flags);
     368//       outspec[i] = median(ma);
     369//       if ( allEQ(ma.getMask(), False) )
     370//         outflag[i] = userflag;// flag data
     371//     }
     372//     outtsys[0] = median(tsysCol.getColumn());
     373//     specColOut.put(outrowCount, outspec);
     374//     flagColOut.put(outrowCount, outflag);
     375//     tsysColOut.put(outrowCount, outtsys);
     376//     Double intsum = sum(intCol.getColumn());
     377//     intColOut.put(outrowCount, intsum);
     378//     ++outrowCount;
     379//     ++iter;
     380    if ( otfscan ) {
     381      MDirection::ScalarColumn dircol ;
     382      dircol.attach( subt, "DIRECTION" ) ;
     383      Int length = subt.nrow() ;
     384      vector< Vector<Double> > dirs ;
     385      vector<int> indexes ;
     386      for ( Int i = 0 ; i < length ; i++ ) {
     387        Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
     388        bool adddir = true ;
     389        for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
     390          //if ( allTrue( t == dirs[j] ) ) {
     391          if ( allNearAbs( t, dirs[j], tol ) ) {
     392            adddir = false ;
     393            break ;
     394          }
     395        }
     396        if ( adddir ) {
     397          dirs.push_back( t ) ;
     398          indexes.push_back( i ) ;
     399        }
     400      }
     401      uInt rowNum = dirs.size() ;
     402      tout.addRow( rowNum );
     403      for ( uInt i = 0 ; i < rowNum ; i++ ) {
     404        TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
     405        if ( avmode != "SCAN") {
     406          //scanColOut.put(outrowCount+i, uInt(0));
     407        }
     408      }
     409      MDirection::ScalarColumn dircolOut ;
     410      dircolOut.attach( tout, "DIRECTION" ) ;
     411      for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
     412        Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
     413        Vector<Float> tmp;
     414        specCol.get(0, tmp);
     415        uInt nchan = tmp.nelements();
     416        // have to do channel by channel here as MaskedArrMath
     417        // doesn't have partialMedians
     418        Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     419        // mask spectra for different DIRECTION
     420        //cout << "irow=" << outrowCount+irow << ": flagged [" ;
     421        for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
     422          Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     423          //if ( t[0] != direction[0] || t[1] != direction[1] ) {
     424          if ( !allNearAbs( t, direction, tol ) ) {
     425            //cout << jrow << " " ;
     426            flags[jrow] = userflag ;
     427          }
     428        }
     429        //cout << "]" << endl ;
     430        //cout << "flags=" << flags << endl ;
     431        Vector<Float> outspec(nchan);
     432        Vector<uChar> outflag(nchan,0);
     433        Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     434        for (uInt i=0; i<nchan; ++i) {
     435          Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     436          MaskedArray<Float> ma = maskedArray(specs,flags);
     437          outspec[i] = median(ma);
     438          if ( allEQ(ma.getMask(), False) )
     439            outflag[i] = userflag;// flag data
     440        }
     441        outtsys[0] = median(tsysCol.getColumn());
     442        specColOut.put(outrowCount+irow, outspec);
     443        flagColOut.put(outrowCount+irow, outflag);
     444        tsysColOut.put(outrowCount+irow, outtsys);
     445        Vector<Double> integ = intCol.getColumn() ;
     446        MaskedArray<Double> mi = maskedArray( integ, flags ) ;
     447        Double intsum = sum(mi);
     448        intColOut.put(outrowCount+irow, intsum);
     449      }
     450      outrowCount += rowNum ;
     451    }
     452    else {
     453      tout.addRow();
     454      TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
     455      if ( avmode != "SCAN") {
     456        scanColOut.put(outrowCount, uInt(0));
     457      }
     458      Vector<Float> tmp;
     459      specCol.get(0, tmp);
     460      uInt nchan = tmp.nelements();
     461      // have to do channel by channel here as MaskedArrMath
     462      // doesn't have partialMedians
     463      Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
     464      Vector<Float> outspec(nchan);
     465      Vector<uChar> outflag(nchan,0);
     466      Vector<Float> outtsys(1);/// @fixme when tsys is channel based
     467      for (uInt i=0; i<nchan; ++i) {
     468        Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
     469        MaskedArray<Float> ma = maskedArray(specs,flags);
     470        outspec[i] = median(ma);
     471        if ( allEQ(ma.getMask(), False) )
     472          outflag[i] = userflag;// flag data
     473      }
     474      outtsys[0] = median(tsysCol.getColumn());
     475      specColOut.put(outrowCount, outspec);
     476      flagColOut.put(outrowCount, outflag);
     477      tsysColOut.put(outrowCount, outtsys);
     478      Double intsum = sum(intCol.getColumn());
     479      intColOut.put(outrowCount, intsum);
     480      ++outrowCount;
     481    }
    291482    ++iter;
    292483  }
    … …  
    395586  convertArray(mask, f);
    396587  return MaskedArray<Float>(s,!mask);
     588}
     589
     590MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
     591                                         const Vector<uChar>& f)
     592{
     593  Vector<Bool> mask;
     594  mask.resize(f.shape());
     595  convertArray(mask, f);
     596  return MaskedArray<Double>(s,!mask);
    397597}
    398598
    … …  
    20912291                    "Use merge first."));
    20922292 
     2293  // check if OTF observation
     2294  String obstype = in[0]->getHeader().obstype ;
     2295  bool otfscan = false ;
     2296  if ( obstype.find( "OTF" ) != String::npos ) {
     2297    //cout << "OTF scan" << endl ;
     2298    otfscan = true ;
     2299  }
     2300
    20932301  CountedPtr<Scantable> out ;     // processed result
    20942302  if ( compel ) {
    20952303    std::vector< CountedPtr<Scantable> > newin ; // input for average process
    2096     uInt insize = in.size() ;    // number of scantables
     2304    uInt insize = in.size() ;    // number of input scantables
    20972305
    20982306    // TEST: do normal average in each table before IF grouping
    … …  
    21242332
    21252333    // check IF frequency coverage
     2334    // freqid: list of FREQ_ID, which is used, in each table 
     2335    // iffreq: list of minimum and maximum frequency for each FREQ_ID in
     2336    //         each table
     2337    // freqid[insize][numIF]
     2338    // freqid: [[id00, id01, ...],
     2339    //          [id10, id11, ...],
     2340    //          ...
     2341    //          [idn0, idn1, ...]]
     2342    // iffreq[insize][numIF*2]
     2343    // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
     2344    //          [min_id10, max_id10, min_id11, max_id11, ...],
     2345    //          ...
     2346    //          [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
    21262347    //cout << "Check IF settings in each table" << endl ;
    21272348    vector< vector<uInt> > freqid( insize );
    … …  
    21602381
    21612382    // IF grouping based on their frequency coverage
     2383    // ifgrp: list of table index and FREQ_ID for all members in each IF group
     2384    // ifgfreq: list of minimum and maximum frequency in each IF group
     2385    // ifgrp[numgrp][nummember*2]
     2386    // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
     2387    //         [table10, freqrow10, table11, freqrow11, ...],
     2388    //         ...
     2389    //         [tablen0, freqrown0, tablen1, freqrown1, ...]]
     2390    // ifgfreq[numgrp*2]
     2391    // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
    21622392    //cout << "IF grouping based on their frequency coverage" << endl ;
    2163 
    2164     // IF group
    2165     // ifgrp[numgrp][nummember*2]
    2166     // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
    2167     //         [table11, freqrow11, table11, freqrow11, ...],
    2168     //         ...
    2169     // ifgfreq[numgrp*2]
    2170     // ifgfreq = [min0, max0, min1, max1, ...]
    21712393    vector< vector<uInt> > ifgrp ;
    21722394    vector<double> ifgfreq ;
    … …  
    22242446    }
    22252447
     2448    // Grouping continuous IF groups (without frequency gap)
     2449    // freqgrp: list of IF group indexes in each frequency group
     2450    // freqrange: list of minimum and maximum frequency in each frequency group
    22262451    // freqgrp[numgrp][nummember]
    2227     // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
     2452    // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
    22282453    //           [ifgrp10, ifgrp11, ifgrp12, ...],
    22292454    //           ...
     2455    //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
    22302456    // freqrange[numgrp*2]
    2231     // freqrange = [min0, max0, min1, max1, ...]
     2457    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
    22322458    vector< vector<uInt> > freqgrp ;
    22332459    double freqrange = 0.0 ;
    … …  
    22722498
    22732499    // membership check
     2500    // groups: list of IF group indexes whose frequency range overlaps with
     2501    //         that of each table and IF
    22742502    // groups[numtable][numIF][nummembership]
     2503    // groups: [[[grp, grp,...], [grp, grp,...],...],
     2504    //          [[grp, grp,...], [grp, grp,...],...],
     2505    //          ...
     2506    //          [[grp, grp,...], [grp, grp,...],...]]
    22752507    vector< vector< vector<uInt> > > groups( insize ) ;
    22762508    for ( uInt i = 0 ; i < insize ; i++ ) {
    … …  
    23852617    //cout << endl ;
    23862618
    2387     // reset SCANNO and IF: IF number is reset by the result of sortation
     2619    // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
    23882620    cout << "All scan number is set to 0" << endl ;
    23892621    //cout << "All IF number is set to IF group index" << endl ;
    … …  
    24962728    // reset spectra and flagtra: align spectral resolution
    24972729    //cout << "Align spectral resolution" << endl ;
     2730    // gmaxdnu: the coarsest frequency resolution in the frequency group
     2731    // gmemid: member index that have a resolution equal to gmaxdnu
     2732    // gmaxdnu[numfreqgrp]
     2733    // gmaxdnu: [dnu0, dnu1, ...]
     2734    // gmemid[numfreqgrp]
     2735    // gmemid: [id0, id1, ...]
    24982736    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
    24992737    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
    … …  
    26482886    ScalarColumn<uInt> freqidColOut ;
    26492887    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
     2888    MDirection::ScalarColumn dirColOut ;
     2889    dirColOut.attach( out->table(), "DIRECTION" ) ;
    26502890    Table &tab = tmpout->table() ;
    2651     TableIterator iter( tab, "POLNO" ) ;
     2891    Block<String> cols(1);
     2892    cols[0] = String("POLNO") ;
     2893    TableIterator iter( tab, cols ) ;
     2894    bool done = false ;
    26522895    vector< vector<uInt> > sizes( freqgrp.size() ) ;
    26532896    while( !iter.pastEnd() ) {
    … …  
    26612904      ScalarColumn<uInt> polnos ;
    26622905      polnos.attach( iter.table(), "POLNO" ) ;
     2906      MDirection::ScalarColumn dircol ;
     2907      dircol.attach( iter.table(), "DIRECTION" ) ;
    26632908      uInt polno = polnos( 0 ) ;
    26642909      //cout << "POLNO iteration: " << polno << endl ;
    2665       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
    2666         sizes[igrp].resize( freqgrp[igrp].size() ) ;
    2667         for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
    2668           for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
    2669             uInt ifno = ifnoCol( irow ) ;
    2670             if ( ifno == freqgrp[igrp][imem] ) {
    2671               Vector<Float> spec = specCols( irow ) ;
    2672               Vector<uChar> flag = flagCols( irow ) ;
    2673               vector<Float> svec ;
    2674               spec.tovector( svec ) ;
    2675               vector<uChar> fvec ;
    2676               flag.tovector( fvec ) ;
    2677               //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
    2678               specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
    2679               flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
    2680               //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
    2681               sizes[igrp][imem] = spec.nelements() ;
    2682             }
    2683           }
    2684         }
    2685         for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
    2686           uInt ifout = ifnoColOut( irow ) ;
    2687           uInt polout = polnoColOut( irow ) ;
    2688           if ( ifout == gmemid[igrp] && polout == polno ) {
    2689             // set SPECTRA and FRAGTRA
    2690             Vector<Float> newspec( specout[igrp] ) ;
    2691             Vector<uChar> newflag( flagout[igrp] ) ;
    2692             specColOut.put( irow, newspec ) ;
    2693             flagColOut.put( irow, newflag ) ;
    2694             // IFNO renumbering
    2695             ifnoColOut.put( irow, igrp ) ;
    2696           }
    2697         }
     2910//       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     2911//      sizes[igrp].resize( freqgrp[igrp].size() ) ;
     2912//      for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     2913//        for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
     2914//          uInt ifno = ifnoCol( irow ) ;
     2915//          if ( ifno == freqgrp[igrp][imem] ) {
     2916//            Vector<Float> spec = specCols( irow ) ;
     2917//            Vector<uChar> flag = flagCols( irow ) ;
     2918//            vector<Float> svec ;
     2919//            spec.tovector( svec ) ;
     2920//            vector<uChar> fvec ;
     2921//            flag.tovector( fvec ) ;
     2922//            //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
     2923//            specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
     2924//            flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
     2925//            //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
     2926//            sizes[igrp][imem] = spec.nelements() ;
     2927//          }
     2928//        }
     2929//      }
     2930//      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     2931//        uInt ifout = ifnoColOut( irow ) ;
     2932//        uInt polout = polnoColOut( irow ) ;
     2933//        if ( ifout == gmemid[igrp] && polout == polno ) {
     2934//          // set SPECTRA and FRAGTRA
     2935//          Vector<Float> newspec( specout[igrp] ) ;
     2936//          Vector<uChar> newflag( flagout[igrp] ) ;
     2937//          specColOut.put( irow, newspec ) ;
     2938//          flagColOut.put( irow, newflag ) ;
     2939//          // IFNO renumbering
     2940//          ifnoColOut.put( irow, igrp ) ;
     2941//        }
     2942//      }
     2943//       }
     2944      // get a list of number of channels for each frequency group member
     2945      if ( !done ) {
     2946        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
     2947          sizes[igrp].resize( freqgrp[igrp].size() ) ;
     2948          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     2949            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
     2950              uInt ifno = ifnoCol( irow ) ;
     2951              if ( ifno == freqgrp[igrp][imem] ) {
     2952                Vector<Float> spec = specCols( irow ) ;
     2953                sizes[igrp][imem] = spec.nelements() ;
     2954                break ;
     2955              }               
     2956            }
     2957          }
     2958        }
     2959        done = true ;
     2960      }
     2961      // combine spectra
     2962      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
     2963        uInt polout = polnoColOut( irow ) ;
     2964        if ( polout == polno ) {
     2965          uInt ifout = ifnoColOut( irow ) ;
     2966          Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
     2967          uInt igrp ;
     2968          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
     2969            if ( ifout == gmemid[jgrp] ) {
     2970              igrp = jgrp ;
     2971              break ;
     2972            }
     2973          }
     2974          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
     2975            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
     2976              uInt ifno = ifnoCol( jrow ) ;
     2977              Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
     2978              //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
     2979              if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
     2980                Vector<Float> spec = specCols( jrow ) ;
     2981                Vector<uChar> flag = flagCols( jrow ) ;
     2982                vector<Float> svec ;
     2983                spec.tovector( svec ) ;
     2984                vector<uChar> fvec ;
     2985                flag.tovector( fvec ) ;
     2986                //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
     2987                specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
     2988                flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
     2989                //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
     2990              }
     2991            }
     2992          }
     2993          // set SPECTRA and FRAGTRA
     2994          Vector<Float> newspec( specout[igrp] ) ;
     2995          Vector<uChar> newflag( flagout[igrp] ) ;
     2996          specColOut.put( irow, newspec ) ;
     2997          flagColOut.put( irow, newflag ) ;
     2998          // IFNO renumbering
     2999          ifnoColOut.put( irow, igrp ) ;
     3000        }
    26983001      }
    26993002      iter++ ;

  • branches/alma/src/STMath.h

    r1603 r1607  
    315315    maskedArray( const casa::Vector<casa::Float>& s,
    316316                 const casa::Vector<casa::uChar>& f );
     317  casa::MaskedArray<casa::Double>
     318    maskedArray( const casa::Vector<casa::Double>& s,
     319                 const casa::Vector<casa::uChar>& f );
    317320  casa::Vector<casa::uChar>
    318321    flagsFromMA(const casa::MaskedArray<casa::Float>& ma);
Note: See TracChangeset for help on using the changeset viewer.