Changeset 2379


Ignore:
Timestamp:
12/21/11 15:42:01 (12 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: Yes

Interface Changes: 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...

Introduced row-by-row gridding.


Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STGrid.cpp

    r2378 r2379  
    1818#include <tables/Tables/ScalarColumn.h>
    1919#include <tables/Tables/ArrayColumn.h>
     20#include <tables/Tables/TableIter.h>
    2021
    2122#include <measures/Measures/MDirection.h>
    … …  
    6061  userSupport_ = -1 ;
    6162  convSampling_ = 100 ;
    62   irow_ = 0 ;
     63  nprocessed_ = 0 ;
     64  nchunk_ = 0 ;
    6365}
    6466
    … …  
    177179  //cout << "convSupport=" << convSupport_ << endl ;
    178180  //cout << "convFunc=" << convFunc << endl ;
    179   Float *conv_p = convFunc.getStorage( deleteConv ) ;
    180 
    181   // Extend grid plane with convSupport_
     181
     182  // grid data
    182183  Int gnx = nx_ ;
    183184  Int gny = ny_ ;
     185  // Extend grid plane with convSupport_
    184186//   Int gnx = nx_+convSupport_*2 ;
    185187//   Int gny = ny_+convSupport_*2 ;
    … …  
    191193  data_ = 0.0 ;
    192194  Array<Float> gwgtArr( data_ ) ;
    193   //Array<Float> gwgtArr( gshape, 0.0 ) ;
    194   Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
    195   Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
    196 
    197   // data selection
    198   selectData( tab_ ) ;
    199   setupArray( tab_ ) ;
    200195
    201196  // data storage
    202   IPosition mshape( 3, npol_, nchan_, 1 ) ;
    203   IPosition vshape( 2, npol_, 1 ) ;
    204   IPosition dshape( 2, 2, 1 ) ;
     197  Int irow = -1 ;
     198  IPosition mshape( 3, npol_, nchan_, nchunk_ ) ;
     199  IPosition vshape( 2, npol_, nchunk_ ) ;
     200  IPosition dshape( 2, 2, nchunk_ ) ;
    205201  Array<Complex> spectra( mshape ) ;
    206202  Array<Double> direction( dshape ) ;
    … …  
    210206  Array<Double> xypos( dshape ) ;
    211207
     208  // parameters for gridding
     209  Int idopsf = 0 ;
     210  Int *chanMap = new Int[nchan_] ;
     211  {
     212    Int *work_p = chanMap ;
     213    for ( Int i = 0 ; i < nchan_ ; i++ ) {
     214      *work_p = i ;
     215      work_p++ ;
     216    }
     217  }
     218  Int *polMap = new Int[npol_] ;
     219  {
     220    Int *work_p = polMap ;
     221    for ( Int i = 0 ; i < npol_ ; i++ ) {
     222      *work_p = i ;
     223      work_p++ ;
     224    }
     225  }
     226  Double *sumw_p = new Double[npol_*nchan_] ;
     227  {
     228    Double *work_p = sumw_p ;
     229    for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
     230      *work_p = 0.0 ;
     231      work_p++ ;
     232    }
     233  }
     234
     235  // for performance check
     236  double eGetDataChunk = 0.0 ;
     237  double eToPixel = 0.0 ;
     238  double eGGridSD = 0.0 ;
     239
    212240  while( !pastEnd() ) {
    213241    // retrieve data
    214     getDataChunk( spectra, direction, flagtra, rflag, weight ) ;
     242    t0 = mathutil::gettimeofday_sec() ;
     243    Int nrow = getDataChunk( spectra, direction, flagtra, rflag, weight ) ;
     244    //os << "nrow = " << nrow << LogIO::POST ;
     245    t1 = mathutil::gettimeofday_sec() ;
     246    eGetDataChunk += t1-t0 ;
    215247   
    216248    // world -> pixel
     249    t0 = mathutil::gettimeofday_sec() ;
    217250    toPixel( direction, xypos ) ;
    218    
     251    t1 = mathutil::gettimeofday_sec() ;
     252    eToPixel += t1-t0 ;
     253   
     254    // prepare pointer
     255    Float *conv_p = convFunc.getStorage( deleteConv ) ;
     256    Complex *gdata_p = gdataArrC.getStorage( deleteDataG ) ;
     257    Float *wdata_p = gwgtArr.getStorage( deleteWgtG ) ;
     258    const Complex *data_p = spectra.getStorage( deleteData ) ;
     259    const Float *wgt_p = weight.getStorage( deleteWgt ) ;
     260    const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
     261    const Int *rflag_p = rflag.getStorage( deleteFlagR ) ;
     262    Double *xypos_p = xypos.getStorage( deletePos ) ;
     263
    219264    // call ggridsd
    220   }
    221  
     265    irow = -1 ;
     266    t0 = mathutil::gettimeofday_sec() ;
     267    ggridsd( xypos_p,
     268             data_p,
     269             &npol_,
     270             &nchan_,
     271             &idopsf,
     272             flag_p,
     273             rflag_p,
     274             wgt_p,
     275             &nrow,
     276             &irow,
     277             gdata_p,
     278             wdata_p,
     279             &gnx,
     280             &gny,
     281             &npol_,
     282             &nchan_,
     283             &convSupport_,
     284             &convSampling_,
     285             conv_p,
     286             chanMap,
     287             polMap,
     288             sumw_p ) ;
     289    t1 = mathutil::gettimeofday_sec() ;
     290    eGGridSD += t1-t0 ;
     291
     292    // finalization
     293    convFunc.putStorage( conv_p, deleteConv ) ;
     294    gdataArrC.putStorage( gdata_p, deleteDataG ) ;
     295    gwgtArr.putStorage( wdata_p, deleteWgtG ) ;
     296    xypos.putStorage( xypos_p, deletePos ) ;
     297    spectra.freeStorage( data_p, deleteData ) ;
     298    weight.freeStorage( wgt_p, deleteWgt ) ;
     299    flagtra.freeStorage( flag_p, deleteFlag ) ;
     300    rflag.freeStorage( rflag_p, deleteFlagR ) ;
     301  }
     302  os << "getDataChunk: elapsed time is " << eGetDataChunk << " sec." << LogIO::POST ;
     303  os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;
     304  os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;
     305
     306
     307  // finalization
     308  delete polMap ;
     309  delete chanMap ;
     310  delete sumw_p ;
     311
    222312  // set data
    223313  setData( gdataArrC, gwgtArr ) ;
     314
    224315}
    225316
    226317Bool STGrid::pastEnd()
    227318{
    228   Bool b = irow_ >= nrow_ ;
     319  LogIO os( LogOrigin("STGrid","pastEnd",WHERE) ) ;
     320  Bool b = nprocessed_ >= nrow_ ;
    229321  return b ;
    230322}
    231323
    232 void STGrid::grid()
     324void STGrid::grid()
     325{
     326  // data selection
     327  selectData() ;
     328  setupArray() ;
     329  sortData() ;
     330
     331  if ( npol_ > 1 ) {
     332    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
     333    os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
     334  }
     335 
     336  if ( wtype_.compare("UNIFORM") != 0 &&
     337       wtype_.compare("TINT") != 0 &&
     338       wtype_.compare("TINTSYS") != 0 ) {
     339    LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
     340    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
     341    wtype_ = "UNIFORM" ;
     342  }
     343
     344  Bool doAll = examine() ;
     345
     346  if ( doAll )
     347    gridAll() ;
     348  else
     349    gridPerRow() ;
     350}
     351
     352Bool STGrid::examine()
     353{
     354  // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
     355  //       by considering data size to be allocated for ggridsd input/output
     356  nchunk_ = 100 ;
     357  Bool b = nchunk_ >= nrow_ ;
     358  return b ;
     359}
     360
     361void STGrid::gridAll()
    233362{
    234363  LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
    … …  
    380509  double t0, t1 ;
    381510  t0 = mathutil::gettimeofday_sec() ;
    382   //data.resize( gdata.shape() ) ;
    383511  uInt len = data_.nelements() ;
    384512  const Complex *w1_p ;
    385513  Float *w2_p ;
    386   Bool b0, b1, b2 ;
     514  Bool b1, b2 ;
    387515  const Complex *gdata_p = gdata.getStorage( b1 ) ;
    388516  Float *gwgt_p = data_.getStorage( b2 ) ;
    … …  
    390518  w2_p = gwgt_p ;
    391519  for ( uInt i = 0 ; i < len ; i++ ) {
    392     //*w2_p = (*w2_p > 0.0) ? ((*w1_p).real() / *w2_p) : 0.0 ;
    393520    if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
    394521    w1_p++ ;
    … …  
    519646}
    520647
    521 void STGrid::selectData( Table &tab )
     648void STGrid::selectData()
    522649{
    523650  Int ifno = ifno_ ;
    … …  
    539666    node = node && taborg.col("SCANNO").in( scanlist_ ) ;
    540667  }
    541 //   Block<String> cols( 3 ) ;
    542 //   cols[0] = "TIME" ;
    543 //   cols[1] = "BEAMNO" ;
    544 //   cols[2] = "POLNO" ;
    545 //   Block<Int> order( 3, Sort::Ascending ) ;
    546 //   tab = taborg( node ).sort( cols, order ) ;
    547   tab = taborg( node ) ;
    548   if ( tab.nrow() == 0 ) {
     668  tab_ = taborg( node ) ;
     669  if ( tab_.nrow() == 0 ) {
    549670    LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
    550671    os << LogIO::SEVERE
    … …  
    554675    os << LogIO::EXCEPTION ;
    555676  }
     677  attach() ;
     678}
     679
     680void STGrid::attach()
     681{
     682  // attach to table
     683  spectraCol_.attach( tab_, "SPECTRA" ) ;
     684  flagtraCol_.attach( tab_, "FLAGTRA" ) ;
     685  directionCol_.attach( tab_, "DIRECTION" ) ;
     686  flagRowCol_.attach( tab_, "FLAGROW" ) ;
     687  tsysCol_.attach( tab_, "TSYS" ) ;
     688  intervalCol_.attach( tab_, "INTERVAL" ) ;
    556689}
    557690
    … …  
    564697  LogIO os( LogOrigin("STGrid","getData",WHERE) ) ;
    565698//   os << "start" << LogIO::POST ;
    566   Table tab ;
    567   selectData( tab ) ;
    568   setupArray( tab ) ;
    569699//   os << "npol_ = " << npol_ << LogIO::POST ;
    570700//   os << "nchan_ = " << nchan_ << LogIO::POST ;
    … …  
    598728  long rincr = npol_ * nchan_ ;
    599729  for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
    600     Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
     730    Table subt = tab_( tab_.col("POLNO") == pollist_[ipol] ) ;
    601731    ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
    602732    ROArrayColumn<Double> directionCol( subt, "DIRECTION" ) ;
    … …  
    670800}
    671801
    672 void STGrid::getDataChunk( Array<Complex> &spectra,
    673                            Array<Double> &direction,
    674                            Array<Int> &flagtra,
    675                            Array<Int> &rflag,
    676                            Array<Float> &weight )
    677 {
    678 }
    679 
    680 void STGrid::setupArray( Table &tab )
     802Int STGrid::getDataChunk( Array<Complex> &spectra,
     803                          Array<Double> &direction,
     804                          Array<Int> &flagtra,
     805                          Array<Int> &rflag,
     806                          Array<Float> &weight )
     807{
     808  Array<Float> spFloat( spectra.shape() ) ;
     809  Array<uChar> fluChar( flagtra.shape() ) ;
     810  Array<uInt> fruInt( rflag.shape() ) ;
     811  Int nrow = getDataChunk( spFloat, direction, fluChar, fruInt, weight ) ;
     812  convertArray( spectra, spFloat ) ;
     813  convertArray( flagtra, fluChar ) ;
     814  convertArray( rflag, fruInt ) ;
     815
     816  return nrow ;
     817}
     818
     819Int STGrid::getDataChunk( Array<Float> &spectra,
     820                          Array<Double> &direction,
     821                          Array<uChar> &flagtra,
     822                          Array<uInt> &rflag,
     823                          Array<Float> &weight )
     824{
     825  LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
     826  Int nrow = min( spectra.shape()[2], nrow_-nprocessed_ ) ;
     827  Array<Float> tsys( spectra.shape() ) ;
     828  Array<Double> tint( rflag.shape() ) ;
     829  IPosition m( 2, 0, 2 ) ;
     830  IPosition v( 1, 0 ) ;
     831  ArrayIterator<Float> spi( spectra, m, False ) ;
     832  ArrayIterator<uChar> fli( flagtra, m, False ) ;
     833  ArrayIterator<Float> tsi( tsys, m, False ) ;
     834  ArrayIterator<Double> di( direction, v ) ;
     835  Bool bfr, bti ;
     836  uInt *fr_p = rflag.getStorage( bfr ) ;
     837  Double *ti_p = tint.getStorage( bti ) ;
     838  uInt *wfr_p = fr_p ;
     839  Double *wti_p = ti_p ;
     840  Int offset = nprocessed_ * npol_ ;
     841  for ( Int irow = 0 ; irow < npol_ * nrow ; irow++ ) {
     842    Array<Float> spSlice = spi.array() ;
     843    Array<uChar> flSlice = fli.array() ;
     844    Array<Float> tsSlice = tsi.array() ;
     845
     846    uInt idx = rows_[offset+irow] ;
     847    spectraCol_.get( idx, spSlice ) ;
     848
     849    flagtraCol_.get( idx, flSlice ) ;
     850    Vector<Float> tmpF = tsysCol_( idx ) ;
     851    if ( tmpF.nelements() == (uInt)nchan_ ) {
     852      tsSlice = tmpF ;
     853    }
     854    else {
     855      tsSlice = tmpF[0] ;
     856    }
     857    *wfr_p = flagRowCol_( idx ) ;
     858    *wti_p = intervalCol_( idx ) ;
     859
     860    spi.next() ;
     861    fli.next() ;
     862    tsi.next() ;
     863    wfr_p++ ;
     864    wti_p++ ;
     865  }
     866  rflag.putStorage( fr_p, bfr ) ;
     867  tint.putStorage( ti_p, bti ) ;
     868
     869  for ( Int irow = 0 ; irow < nrow ; irow += npol_ ) {
     870    uInt idx = rows_[offset+irow] ;
     871    directionCol_.get( idx, di.array() ) ;
     872    di.next() ;
     873  }
     874
     875  getWeight( weight, tsys, tint ) ;
     876
     877  nprocessed_ += nrow ;
     878 
     879  return nrow ;
     880}
     881
     882void STGrid::setupArray()
    681883{
    682884  LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
    683   ROScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
     885  ROScalarColumn<uInt> polnoCol( tab_, "POLNO" ) ;
    684886  Vector<uInt> pols = polnoCol.getColumn() ;
    685887  //os << pols << LogIO::POST ;
    … …  
    714916    os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
    715917  }
    716   nrow_ = tab.nrow() / npolOrg ;
    717   ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
     918  nrow_ = tab_.nrow() / npolOrg ;
     919  ROArrayColumn<uChar> tmpCol( tab_, "FLAGTRA" ) ;
    718920  nchan_ = tmpCol( 0 ).nelements() ;
    719921//   os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
    … …  
    727929{
    728930  LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
    729   double t0, t1 ;
    730   t0 = mathutil::gettimeofday_sec() ;
     931//   double t0, t1 ;
     932//   t0 = mathutil::gettimeofday_sec() ;
    731933  // resize
    732934  IPosition refShape = tsys.shape() ;
    … …  
    797999  else {
    7981000    //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
    799     os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
     1001    //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
    8001002    w = 1.0 ;
    8011003  }
    8021004
    803   if ( npol_ > 1 ) {
    804     //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
    805     os << LogIO::WARN << "STGrid doesn't support assigning polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
    806   }
    807   t1 = mathutil::gettimeofday_sec() ;
    808   os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
     1005//   t1 = mathutil::gettimeofday_sec() ;
     1006//   os << "getWeight: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
    8091007}
    8101008
    … …  
    10631261  tab = Table( name, Table::Update ) ;
    10641262}
    1065 }
     1263
     1264void STGrid::sortData()
     1265{
     1266  LogIO os( LogOrigin("STGRid","sortData",WHERE) ) ;
     1267  double t0, t1 ;
     1268  t0 = mathutil::gettimeofday_sec() ;
     1269  rows_.resize( npol_*nrow_ ) ;
     1270  Table tab = tab_( tab_.col("POLNO").in(pollist_) ) ;
     1271  Block<String> cols( 2 ) ;
     1272  cols[0] = "TIME" ;
     1273  cols[1] = "BEAMNO" ;
     1274  TableIterator iter( tab, cols ) ;
     1275  uInt idx = 0 ;
     1276  while( !iter.pastEnd() ) {
     1277    Table t = iter.table().sort( "POLNO" ) ;
     1278    Vector<uInt> rows = t.rowNumbers() ;
     1279    //os << "rows=" << rows << LogIO::POST ;
     1280    for ( uInt i = 0 ; i < rows.nelements() ; i++ ) {
     1281      rows_[idx] = rows[i] ;
     1282      idx++ ;
     1283    }
     1284    iter.next() ;
     1285  }
     1286  t1 = mathutil::gettimeofday_sec() ;
     1287  os << "sortData: elapsed time is " << t1-t0 << " sec" << LogIO::POST ;
     1288  //os << "rows_ = " << rows_ << LogIO::POST ;
     1289}
     1290}

  • trunk/src/STGrid.h

    r2378 r2379  
    2929
    3030#include <tables/Tables/Table.h>
    31 // #include <tables/Tables/ScalarColumn.h>
    32 // #include <tables/Tables/ArrayColumn.h>
     31#include <tables/Tables/ScalarColumn.h>
     32#include <tables/Tables/ArrayColumn.h>
    3333
    3434// #include <measures/Measures/MDirection.h>
    … …  
    6666
    6767  void grid() ;
     68  void gridAll() ;
    6869  void gridPerRow() ;
    6970 
    … …  
    9697                Array<uInt> &rflag,
    9798                Array<Float> &weight ) ;
    98   void getDataChunk( Array<Complex> &spectra,
    99                      Array<Double> &direction,
    100                      Array<Int> &flagtra,
    101                      Array<Int> &rflag,
    102                      Array<Float> &weight ) ;
     99  Int getDataChunk( Array<Complex> &spectra,
     100                    Array<Double> &direction,
     101                    Array<Int> &flagtra,
     102                    Array<Int> &rflag,
     103                    Array<Float> &weight ) ;
     104  Int getDataChunk( Array<Float> &spectra,
     105                    Array<Double> &direction,
     106                    Array<uChar> &flagtra,
     107                    Array<uInt> &rflag,
     108                    Array<Float> &weight ) ;
    103109
    104110  void getWeight( Array<Float> &w,
    … …  
    116122  void pbFunc( Vector<Float> &convFunc ) ;
    117123  void setConvFunc( Vector<Float> &convFunc ) ;
    118   void selectData( Table &tab ) ;
    119124
    120125  Float polMean( const Float *p ) ;
    121126  Double polMean( const Double *p ) ;
    122127
    123   void setupArray( Table &tab ) ;
    124 
    125128  void prepareTable( Table &tab, String &name ) ;
    126129
    127130  Bool pastEnd() ;
     131
     132  void selectData() ;
     133  void setupArray() ;
     134  void sortData() ;
     135
     136  Bool examine() ;
     137  void attach() ;
    128138
    129139
    … …  
    148158
    149159  Table tab_ ;
    150   Int irow_ ;
     160  ROArrayColumn<Float> spectraCol_ ;
     161  ROArrayColumn<uChar> flagtraCol_ ;
     162  ROArrayColumn<Double> directionCol_ ;
     163  ROScalarColumn<uInt> flagRowCol_ ;
     164  ROArrayColumn<Float> tsysCol_ ;
     165  ROScalarColumn<Double> intervalCol_ ;
     166  Int nprocessed_ ;
     167  Vector<uInt> rows_ ;
     168  Int nchunk_ ;
    151169};
    152170}
Note: See TracChangeset for help on using the changeset viewer.