Changeset 2964


Ignore:
Timestamp:
06/30/14 17:13:19 (10 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-6585, CAS-6571

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Code refactoring.


Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STApplyCal.cpp

    r2963 r2964  
    4646
    4747namespace {
    48 // utility functions
    49 Vector<uInt> indexSort(Vector<Double> &t)
    50 {
    51   Sort sort;
    52   sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
    53   Vector<uInt> idx;
    54   sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
    55   return idx;
    56 }
    57 
    58 void getSortedData(const Vector<Double> key, const Vector<Float> data,
    59                    const Vector<uChar> flag,
    60                    Vector<Double> sortedKey, Vector<Float> sortedData,
    61                    Vector<uChar> sortedFlag)
    62 {
    63 }
     48template<class Accessor, class Type>
     49class AccessorInterface
     50{
     51public:
     52  typedef Type TableType;
     53  static void GetSortedData(const vector<TableType *> &tablelist,
     54                           const Vector<uInt> &tableIndex,
     55                           const uInt nrow,
     56                           const uInt nchan,
     57                           Vector<Double> &time,
     58                           Matrix<Float> &data,
     59                           Matrix<uChar> &flag)
     60  {
     61    Vector<Double> timeUnsorted;
     62    Matrix<Float> dataUnsorted;
     63    Matrix<uChar> flagUnsorted;
     64    GetFromTable(tablelist, tableIndex, nrow, nchan,
     65                 timeUnsorted, dataUnsorted, flagUnsorted);
     66    SortData(timeUnsorted, dataUnsorted, flagUnsorted,
     67             time, data, flag);
     68  }
     69private:
     70  static void GetFromTable(const vector<TableType *> &tableList,
     71                           const Vector<uInt> &tableIndex,
     72                           const uInt nrow,
     73                           const uInt nchan,
     74                           Vector<Double> &time,
     75                           Matrix<Float> &data,
     76                           Matrix<uChar> &flag)
     77  {
     78    time.resize(nrow, False);
     79    const IPosition shape(2, nrow, nchan);
     80    data.resize(shape, False);
     81    flag.resize(shape, False);
     82    uInt rowIndex = 0;
     83    for (uInt i = 0 ; i < tableIndex.nelements(); i++) {
     84      TableType *p = tableList[tableIndex[i]];
     85      Vector<Double> t = Accessor::GetTime(p);
     86      Matrix<Float> dt = Accessor::GetData(p);
     87      Matrix<uChar> fl = Accessor::GetFlag(p);
     88      for (uInt j = 0; j < t.nelements(); j++) {
     89        time[rowIndex] = t[j];
     90        data.row(rowIndex) = dt.column(j);
     91        flag.row(rowIndex) = fl.column(j);
     92        rowIndex++;
     93      }
     94    }   
     95  }
     96
     97  static Vector<uInt> IndexSort(const Vector<Double> &t)
     98  {
     99    Sort sort;
     100    sort.sortKey(&t[0], TpDouble, 0, Sort::Ascending);
     101    Vector<uInt> idx;
     102    sort.sort(idx, t.nelements(), Sort::QuickSort|Sort::NoDuplicates);
     103    return idx;
     104  }
     105
     106  static void SortData(const Vector<Double> &key, const Matrix<Float> &data,
     107                         const Matrix<uChar> &flag,
     108                         Vector<Double> &sortedKey, Matrix<Float> &sortedData,
     109                         Matrix<uChar> &sortedFlag)
     110  {
     111    Vector<uInt> sortIndex = IndexSort(key);
     112    uInt len = sortIndex.nelements();
     113    IPosition shape = data.shape();
     114    shape[0] = len;
     115    Int64 nelements = shape.product();
     116    sortedKey.takeStorage(IPosition(1, len), new Double[len], TAKE_OVER);
     117    sortedData.takeStorage(shape, new Float[nelements], TAKE_OVER);
     118    sortedFlag.takeStorage(shape, new uChar[nelements], TAKE_OVER);
     119    for (uInt i = 0 ; i < len; i++) {
     120      sortedKey[i] = key[sortIndex[i]];
     121    }
     122    for (uInt i = 0; i < len; ++i) {
     123      sortedData.row(i) = data.row(sortIndex[i]);
     124      sortedFlag.row(i) = flag.row(sortIndex[i]);
     125    }
     126  }
     127
     128};
     129
     130class SkyTableAccessor : public AccessorInterface<SkyTableAccessor, asap::STCalSkyTable>
     131{
     132public:
     133  static Vector<Double> GetTime(const TableType *t) {return t->getTime();}
     134  static Matrix<Float> GetData(const TableType *t) {return t->getSpectra();}
     135  static Matrix<uChar> GetFlag(const TableType *t) {return t->getFlagtra();}
     136};
     137
     138class TsysTableAccessor : public AccessorInterface<TsysTableAccessor, asap::STCalTsysTable>
     139{
     140public:
     141  static Vector<Double> GetTime(const TableType *t) {return t->getTime();}
     142  static Matrix<Float> GetData(const TableType *t) {return t->getTsys();}
     143  static Matrix<uChar> GetFlag(const TableType *t) {return t->getFlagtra();}
     144};
    64145}
    65146
    … …  
    277358
    278359  // apply selection to apply tables
    279   uInt nrowSky = 0;
    280   uInt nrowTsys = 0;
     360  uInt nrowSkyTotal = 0;
     361  uInt nrowTsysTotal = 0;
    281362  for (uInt i = 0; i < skylist.nelements(); i++) {
    282363    skytable_[skylist[i]]->setSelection(sel);
    283     nrowSky += skytable_[skylist[i]]->nrow();
    284     os_ << "nrowSky=" << nrowSky << LogIO::POST;
     364    nrowSkyTotal += skytable_[skylist[i]]->nrow();
     365    os_ << "nrowSkyTotal=" << nrowSkyTotal << LogIO::POST;
    285366  }
    286367
    287368  // Skip IFNO without sky data
    288   if (nrowSky == 0)
     369  if (nrowSkyTotal == 0)
    289370    return;
    290371
    … …  
    309390      tsystable_[i]->setSelection(sel);
    310391      uInt nrowThisTsys = tsystable_[i]->nrow();
    311       nrowTsys += nrowThisTsys;
     392      nrowTsysTotal += nrowThisTsys;
    312393      if (nrowThisTsys > 0 && nchanTsys == 0) {
    313394        nchanTsys = tsystable_[i]->nchan(tsysifno);
    … …  
    315396      }
    316397    }
    317     interpolatorF_->setX(ftsys.data(), nchanTsys);
    318398  }
    319399
    320400  uInt nchanSp = skytable_[skylist[0]]->nchan(ifno);
    321   uInt nrowSkySorted = nrowSky;
    322   Vector<Double> timeSkySorted;
    323   Matrix<Float> spoffSorted;
    324   Matrix<uChar> flagoffSorted;
    325   {
    326     Vector<Double> timeSky(nrowSky);
    327     Matrix<Float> spoff(nrowSky, nchanSp);
    328     Matrix<uChar> flagoff(nrowSky, nchanSp);
    329     nrowSky = 0;
    330     for (uInt i = 0 ; i < skylist.nelements(); i++) {
    331       STCalSkyTable *p = skytable_[skylist[i]];
    332       Vector<Double> t = p->getTime();
    333       Matrix<Float> sp = p->getSpectra();
    334       Matrix<uChar> fl = p->getFlagtra();
    335       for (uInt j = 0; j < t.nelements(); j++) {
    336         timeSky[nrowSky] = t[j];
    337         spoff.row(nrowSky) = sp.column(j);
    338         flagoff.row(nrowSky) = fl.column(j);
    339         nrowSky++;
    340       }
    341     }
    342    
    343     Vector<uInt> skyIdx = indexSort(timeSky);
    344     nrowSkySorted = skyIdx.nelements();
    345    
    346     timeSkySorted.takeStorage(IPosition(1, nrowSkySorted),
    347                               new Double[nrowSkySorted],
    348                               TAKE_OVER);
    349     for (uInt i = 0 ; i < nrowSkySorted; i++) {
    350       timeSkySorted[i] = timeSky[skyIdx[i]];
    351     }
    352     interpolatorS_->setX(timeSkySorted.data(), nrowSkySorted);
    353    
    354     spoffSorted.takeStorage(IPosition(2, nrowSky, nchanSp),
    355                             new Float[nrowSky * nchanSp],
    356                             TAKE_OVER);
    357     flagoffSorted.takeStorage(IPosition(2, nrowSkySorted, nchanSp),
    358                               new uChar[nrowSkySorted * nchanSp],
    359                               TAKE_OVER);
    360     for (uInt i = 0 ; i < nrowSky; i++) {
    361       spoffSorted.row(i) = spoff.row(skyIdx[i]);
    362       flagoffSorted.row(i) = flagoff.row(skyIdx[i]);
    363     }
    364   }
    365 
    366   uInt nrowTsysSorted = nrowTsys;
    367   Matrix<Float> tsysSorted;
    368   Matrix<uChar> flagtsysSorted;
    369   Vector<Double> timeTsysSorted;
     401  uInt nrowSky = nrowSkyTotal;
     402  Vector<Double> timeSky;
     403  Matrix<Float> spoff;
     404  Matrix<uChar> flagoff;
     405  SkyTableAccessor::GetSortedData(skytable_, skylist,
     406                                  nrowSkyTotal, nchanSp,
     407                                  timeSky, spoff, flagoff);
     408  nrowSky = timeSky.nelements();
     409
     410  uInt nrowTsys = nrowTsysTotal;
     411  Vector<Double> timeTsys;
     412  Matrix<Float> tsys;
     413  Matrix<uChar> flagtsys;
    370414  if (doTsys) {
    371415    //os_ << "doTsys" << LogIO::POST;
    372     Vector<Double> timeTsys(nrowTsys);
    373     Matrix<Float> tsys(nrowTsys, nchanTsys);
    374     Matrix<uChar> flagtsys(nrowTsys, nchanTsys);
    375     tsysSorted.takeStorage(IPosition(2, nrowTsys, nchanTsys),
    376                            new Float[nrowTsys * nchanTsys],
    377                            TAKE_OVER);
    378     nrowTsys = 0;
    379     for (uInt i = 0 ; i < tsystable_.size(); i++) {
    380       STCalTsysTable *p = tsystable_[i];
    381       Vector<Double> t = p->getTime();
    382       Matrix<Float> ts = p->getTsys();
    383       Matrix<uChar> fl = p->getFlagtra();
    384       for (uInt j = 0; j < t.nelements(); j++) {
    385         timeTsys[nrowTsys] = t[j];
    386         tsys.row(nrowTsys) = ts.column(j);
    387         flagtsys.row(nrowTsys) = fl.column(j);
    388         nrowTsys++;
    389       }
    390     }
    391     Vector<uInt> tsysIdx = indexSort(timeTsys);
    392     nrowTsysSorted = tsysIdx.nelements();
    393 
    394     timeTsysSorted.takeStorage(IPosition(1, nrowTsysSorted),
    395                                new Double[nrowTsysSorted],
    396                                TAKE_OVER);
    397     flagtsysSorted.takeStorage(IPosition(2, nrowTsysSorted, nchanTsys),
    398                                new uChar[nrowTsysSorted * nchanTsys],
    399                                TAKE_OVER);
    400     for (uInt i = 0 ; i < nrowTsysSorted; i++) {
    401       timeTsysSorted[i] = timeTsys[tsysIdx[i]];
    402     }
    403     interpolatorT_->setX(timeTsysSorted.data(), nrowTsysSorted);
    404 
    405     for (uInt i = 0; i < nrowTsys; ++i) {
    406       tsysSorted.row(i) = tsys.row(tsysIdx[i]);
    407       flagtsysSorted.row(i) = flagtsys.row(tsysIdx[i]);
    408     }
     416    Vector<uInt> tsyslist(tsystable_.size());
     417    indgen(tsyslist);
     418    TsysTableAccessor::GetSortedData(tsystable_, tsyslist,
     419                                     nrowTsysTotal, nchanTsys,
     420                                     timeTsys, tsys, flagtsys);
     421    nrowTsys = timeTsys.nelements();
    409422  }
    410423
    … …  
    414427  ArrayColumn<Float> tsysCol(tab, "TSYS");
    415428  ScalarColumn<Double> timeCol(tab, "TIME");
    416   //Vector<Float> on;
    417429
    418430  // Array for scaling factor (aka Tsys)
    … …  
    427439
    428440  // working array for interpolation with flags
    429   uInt arraySize = max(max(nrowTsysSorted, nchanTsys), nrowSkySorted);
     441  uInt arraySize = max(max(nrowTsys, nchanTsys), nrowSky);
    430442  Vector<Double> xwork(IPosition(1, arraySize), new Double[arraySize], TAKE_OVER);
    431443  Vector<Float> ywork(IPosition(1, arraySize), new Float[arraySize], TAKE_OVER);
    … …  
    447459    Float *ywork_p = ywork.data();
    448460    for (uInt ichan = 0; ichan < nchanSp; ichan++) {
    449       Float *tmpY = &(spoffSorted.data()[ichan * nrowSkySorted]);
    450       uChar *tmpF = &(flagoffSorted.data()[ichan * nrowSkySorted]);
     461      Float *tmpY = &(spoff.data()[ichan * nrowSky]);
     462      uChar *tmpF = &(flagoff.data()[ichan * nrowSky]);
    451463      uInt wnrow = 0;
    452       for (uInt ir = 0; ir < nrowSkySorted; ++ir) {
     464      for (uInt ir = 0; ir < nrowSky; ++ir) {
    453465        if (tmpF[ir] == 0) {
    454           xwork_p[wnrow] = timeSkySorted.data()[ir];
     466          xwork_p[wnrow] = timeSky.data()[ir];
    455467          ywork_p[wnrow] = tmpY[ir];
    456468          wnrow++;
    457469        }
    458470      }
    459       //if (allNE(flagoffSorted.column(ichan), (uChar)0)) {
    460471      if (wnrow > 0) {
    461472        // any valid reference data
    … …  
    465476        // no valid reference data
    466477        // interpolate data regardless of flag
    467         interpolatorS_->setData(timeSkySorted.data(), tmpY, nrowSkySorted);
     478        interpolatorS_->setData(timeSky.data(), tmpY, nrowSky);
    468479        // flag this channel for calibrated data
    469480        flag[ichan] = 1 << 7; // user flag
    470481      }
    471       //interpolatorS_->setY(tmpY, nrowSkySorted);
    472482      iOff[ichan] = interpolatorS_->interpolate(t0);
    473483    }
    … …  
    481491      uChar *fwork_p = fwork.data();
    482492      for (uInt ichan = 0; ichan < nchanTsys; ichan++) {
    483         Float *tmpY = &(tsysSorted.data()[ichan * nrowTsysSorted]);
    484         uChar *tmpF = &(flagtsysSorted.data()[ichan * nrowTsysSorted]);
     493        Float *tmpY = &(tsys.data()[ichan * nrowTsys]);
     494        uChar *tmpF = &(flagtsys.data()[ichan * nrowTsys]);
    485495        uInt wnrow = 0;
    486         for (uInt ir = 0; ir < nrowTsysSorted; ++ir) {
     496        for (uInt ir = 0; ir < nrowTsys; ++ir) {
    487497          if (tmpF[ir] == 0) {
    488             xwork_p[wnrow] = timeTsysSorted.data()[ir];
     498            xwork_p[wnrow] = timeTsys.data()[ir];
    489499            ywork_p[wnrow] = tmpY[ir];
    490500            wnrow++;
    … …  
    493503        if (wnrow > 0) {
    494504          // any valid value exists
    495           //interpolatorT_->setY(tmpY, nrowTsysSorted);
    496505          interpolatorT_->setData(xwork_p, ywork_p, wnrow);
    497506          iTsysT[ichan] = interpolatorT_->interpolate(t0);

  • trunk/src/STApplyTable.h

    r2720 r2964  
    5959  casa::uInt nrow() {return table_.nrow();}
    6060
    61   casa::Vector<casa::uInt> getScan() {return scanCol_.getColumn();}
    62   casa::Vector<casa::uInt> getCycle() {return cycleCol_.getColumn();}
    63   casa::Vector<casa::uInt> getBeam() {return beamCol_.getColumn();}
    64   casa::Vector<casa::uInt> getIF() {return ifCol_.getColumn();}
    65   casa::Vector<casa::uInt> getPol() {return polCol_.getColumn();}
    66   casa::Vector<casa::Double> getTime() {return timeCol_.getColumn();}
     61  casa::Vector<casa::uInt> getScan() const {return scanCol_.getColumn();}
     62  casa::Vector<casa::uInt> getCycle() const {return cycleCol_.getColumn();}
     63  casa::Vector<casa::uInt> getBeam() const {return beamCol_.getColumn();}
     64  casa::Vector<casa::uInt> getIF() const {return ifCol_.getColumn();}
     65  casa::Vector<casa::uInt> getPol() const {return polCol_.getColumn();}
     66  casa::Vector<casa::Double> getTime() const {return timeCol_.getColumn();}
    6767
    6868  void setSelection(STSelector &sel, bool sortByTime=false);

  • trunk/src/STCalSkyTable.h

    r2958 r2964  
    5555                  const casa::Vector<casa::uChar> &flagtra);
    5656 
    57   casa::Vector<casa::Float> getElevation() {return elCol_.getColumn();}
    58   casa::Matrix<casa::Float> getSpectra() {return spectraCol_.getColumn();}
    59   casa::Matrix<casa::uChar> getFlagtra() {return flagtraCol_.getColumn();}
     57  casa::Vector<casa::Float> getElevation() const {return elCol_.getColumn();}
     58  casa::Matrix<casa::Float> getSpectra() const {return spectraCol_.getColumn();}
     59  casa::Matrix<casa::uChar> getFlagtra() const {return flagtraCol_.getColumn();}
    6060  casa::uInt nchan(casa::uInt ifno);
    6161

  • trunk/src/STCalTsysTable.h

    r2958 r2964  
    5353                  const casa::Vector<casa::uChar> &flagtra);
    5454 
    55   casa::Vector<casa::Float> getElevation() {return elCol_.getColumn();}
    56   casa::Matrix<casa::Float> getTsys() {return tsysCol_.getColumn();}
    57   casa::Matrix<casa::uChar> getFlagtra() {return flagtraCol_.getColumn();}
     55  casa::Vector<casa::Float> getElevation() const {return elCol_.getColumn();}
     56  casa::Matrix<casa::Float> getTsys() const {return tsysCol_.getColumn();}
     57  casa::Matrix<casa::uChar> getFlagtra() const {return flagtraCol_.getColumn();}
    5858  casa::uInt nchan(casa::uInt ifno);
    5959
Note: See TracChangeset for help on using the changeset viewer.