Changeset 2361 for trunk/src


Ignore:
Timestamp:
12/06/11 19:18:40 (13 years ago)
Author:
Takeshi Nakazato
Message:

New Development: No

JIRA Issue: Yes CAS-2816

Ready for Test: No

Interface Changes: Yes

What Interface Changed: Added method to set weight type

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Weighting type can be specified.
Supported types are:

UNIFORM
TINT
TSYS
TINTSYS


Location:
trunk/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STGrid.cpp

    r2360 r2361  
    77#include <casa/Arrays/Cube.h>
    88#include <casa/Arrays/ArrayMath.h>
    9 // #include <casa/Inputs/Input.h>
     9#include <casa/Arrays/ArrayPartMath.h>
    1010#include <casa/Quanta/Quantum.h>
    1111#include <casa/Quanta/QuantumHolder.h>
    1212#include <casa/Utilities/CountedPtr.h>
     13#include <casa/Logging/LogIO.h>
    1314
    1415#include <tables/Tables/Table.h>
    … …  
    4647  nx_ = -1 ;
    4748  ny_ = -1 ;
     49  npol_ = 0 ;
     50  nchan_ = 0 ;
     51  nrow_ = 0 ;
     52  ngrid_ = 0 ;
    4853  cellx_ = 0.0 ;
    4954  celly_ = 0.0 ;
    5055  center_ = Vector<Double> ( 2, 0.0 ) ;
    5156  convType_ = "BOX" ;
     57  wtype_ = "UNIFORM" ;
    5258  convSupport_ = -1 ;
    5359  userSupport_ = -1 ;
    … …  
    6975  pollist_.assign( Vector<uInt>( pols ) ) ;
    7076  cout << "pollist_ = " << pollist_ << endl ;
     77}
     78
     79void STGrid::setWeight( const string wType )
     80{
     81  wtype_ = String( wType ) ;
     82  wtype_.upcase() ;
     83  cout << "wtype_ = " << wtype_ << endl ;
    7184}
    7285
    … …  
    133146void STGrid::grid()
    134147{
     148  LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
    135149
    136150  // retrieve data
    … …  
    139153  Cube<uChar> flagtra ;
    140154  Matrix<uInt> rflag ;
    141   getData( infile_, spectra, direction, flagtra, rflag ) ;
     155  Matrix<Float> weight ;
     156  getData( infile_, spectra, direction, flagtra, rflag, weight ) ;
    142157  IPosition sshape = spectra.shape() ;
    143   Int npol = sshape[0] ;
    144   Int nchan = sshape[1] ;
    145   Int nrow = sshape[2] ;
    146   cout << "spectra.shape()=" << spectra.shape() << endl ;
    147   cout << "max(spectra) = " << max(spectra) << endl ;
     158  //os << "spectra.shape()=" << spectra.shape() << LogIO::POST ;
     159  //os << "max(spectra) = " << max(spectra) << LogIO::POST ;
     160  //os << "weight = " << weight << LogIO::POST ;
    148161
    149162  // flagtra: uChar -> Int
    150163  // rflag: uInt -> Int
    151 //   Matrix<Int> flagI ;
    152 //   Vector<Int> rflagI ;
    153164  Cube<Int> flagI ;
    154165  Matrix<Int> rflagI ;
    … …  
    156167  toInt( &rflag, &rflagI ) ;
    157168 
    158   cout << "flagI.shape() = " << flagI.shape() << endl ;
    159   cout << "rflagI.shape() = " << rflagI.shape() << endl ;
    160  
    161169  // grid parameter
    162 //   cout << "----------" << endl ;
    163 //   cout << "Grid parameter summary" << endl ;
    164 //   cout << "   (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
    165 //   cout << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
    166 //   cout << "   center = " << center_ << endl ;
    167 //   cout << "----------" << endl ;
     170//   os << "----------" << endl ;
     171//   os << "Grid parameter summary" << endl ;
     172//   os << "   (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
     173//   os << "   (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
     174//   os << "   center = " << center_ << endl ;
     175//   os << "----------" << LogIO::POST ;
    168176
    169177  // convolution kernel
    … …  
    176184  Matrix<Double> xypos( direction.shape(), 0.0 ) ;
    177185  toPixel( direction, xypos ) ; 
    178   //cout << "max(xypos.row(0))=" << max(xypos.row(0)) << endl ;
    179   //cout << "min(xypos.row(0))=" << min(xypos.row(0)) << endl ;
    180   //cout << "max(xypos.row(1))=" << max(xypos.row(1)) << endl ;
    181   //cout << "min(xypos.row(1))=" << min(xypos.row(1)) << endl ;
    182 //   for ( Int irow = 0 ; irow < nrow ; irow++ ) {
    183 //     if ( spectra.column( irow )(0) > 0.0 ) {
    184 //       cout << irow << ": xypos=" << xypos.column( irow )
    185 //            << " data = " << spectra.column( irow ) << endl ;
    186 //     }
    187 //   }
    188186 
    189   // weighting factor
    190   Matrix<Float> weight( IPosition( 2, nchan, nrow ), 1.0 ) ;
    191 
    192 //   // call ggridsd
     187  // call ggridsd
    193188  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
    194189  Double *xypos_p = xypos.getStorage( deletePos ) ;
    195   //Matrix<Complex> dataC( spectra.shape(), 0.0 ) ;
    196190  Cube<Complex> dataC( spectra.shape(), 0.0 ) ;
    197191  setReal( dataC, spectra ) ;
    … …  
    201195  const Int *rflag_p = rflagI.getStorage( deleteFlagR ) ;
    202196  Float *conv_p = convFunc.getStorage( deleteConv ) ;
    203   //Int npol = 1 ;
    204197  // Extend grid plane with convSupport_
    205198  //IPosition gshape( 4, nx_, ny_, npol, nchan ) ;
    206199  Int gnx = nx_+convSupport_*2 ;
    207200  Int gny = ny_+convSupport_*2 ;
    208   IPosition gshape( 4, gnx, gny, npol, nchan ) ;
     201  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
    209202  Array<Complex> gdataArrC( gshape, 0.0 ) ;
    210203  Array<Float> gwgtArr( gshape, 0.0 ) ;
    … …  
    213206  Int idopsf = 0 ;
    214207  Int irow = -1 ;
    215   Int *chanMap = new Int[nchan] ;
     208  Int *chanMap = new Int[nchan_] ;
    216209  {
    217210    Int *work_p = chanMap ;
    218     for ( Int i = 0 ; i < nchan ; i++ ) {
     211    for ( Int i = 0 ; i < nchan_ ; i++ ) {
    219212      *work_p = i ;
    220213      work_p++ ;
    221214    }
    222215  }
    223   Int *polMap = new Int[npol] ;
     216  Int *polMap = new Int[npol_] ;
    224217  {
    225218    Int *work_p = polMap ;
    226     for ( Int i = 0 ; i < npol ; i++ ) {
     219    for ( Int i = 0 ; i < npol_ ; i++ ) {
    227220      *work_p = i ;
    228221      work_p++ ;
    229222    }
    230223  }
    231   Double *sumw_p = new Double[npol*nchan] ;
     224  Double *sumw_p = new Double[npol_*nchan_] ;
    232225  {
    233226    Double *work_p = sumw_p ;
    234     for ( Int i = 0 ; i < npol*nchan ; i++ ) {
     227    for ( Int i = 0 ; i < npol_*nchan_ ; i++ ) {
    235228      *work_p = 0.0 ;
    236229      work_p++ ;
    … …  
    239232  ggridsd( xypos_p,
    240233           data_p,
    241            &npol,
    242            &nchan,
     234           &npol_,
     235           &nchan_,
    243236           &idopsf,
    244237           flag_p,
    245238           rflag_p,
    246239           wgt_p,
    247            &nrow,
     240           &nrow_,
    248241           &irow,
    249242           gdata_p,
    … …  
    251244           &gnx,
    252245           &gny,
    253            &npol,
    254            &nchan,
     246           &npol_,
     247           &nchan_,
    255248           &convSupport_,
    256249           &convSampling_,
    … …  
    274267  for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
    275268    for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
    276       for ( Int ip = 0 ; ip < npol ; ip++ ) {
    277         for ( Int ic = 0 ; ic < nchan ; ic++ ) {
     269      for ( Int ip = 0 ; ip < npol_ ; ip++ ) {
     270        for ( Int ic = 0 ; ic < nchan_ ; ic++ ) {
    278271          IPosition pos( 4, ix, iy, ip, ic ) ;
    279           //if ( gwgtArr( pos ) > 0.0 )
    280           //  data_( pos ) = gdataArr( pos ) / gwgtArr( pos ) ;
    281272          IPosition gpos( 4, ix+convSupport_, iy+convSupport_, ip, ic ) ;
    282273          if ( gwgtArr( gpos ) > 0.0 )
    … …  
    286277    }
    287278  }
    288   Matrix<Double> sumWeight( IPosition( 2, npol, nchan ), sumw_p, TAKE_OVER ) ;
     279  //Matrix<Double> sumWeight( IPosition( 2, npol_, nchan_ ), sumw_p, TAKE_OVER ) ;
     280  delete sumw_p ;
    289281  //cout << "sumWeight = " << sumWeight << endl ;
    290282  //cout << "gdataArr = " << gdataArr << endl ;
    … …  
    396388                      Matrix<Double> &direction,
    397389                      Cube<uChar> &flagtra,
    398                       Matrix<uInt> &rflag )
     390                      Matrix<uInt> &rflag,
     391                      Matrix<Float> &weight )
    399392{
    400393  Table tab( infile ) ;
    … …  
    426419    pollist_ = newlist ;
    427420  }
    428   uInt npol = pollist_.size() ;
     421  npol_ = pollist_.size() ;
    429422  ROArrayColumn<uChar> tmpCol( tab, "FLAGTRA" ) ;
    430   uInt nchan = tmpCol( 0 ).nelements() ;
    431   uInt nrow = tab.nrow() / npolOrg ;
    432   cout << "npol = " << npol << endl ;
    433   cout << "nchan = " << nchan << endl ;
    434   cout << "nrow = " << nrow << endl ;
    435   spectra.resize( npol, nchan, nrow ) ;
    436   flagtra.resize( npol, nchan, nrow ) ;
    437   rflag.resize( npol, nrow ) ;
    438   cout << "pollist_=" << pollist_ << endl ;
    439   for ( uInt ipol = 0 ; ipol < npol ; ipol++ ) {
     423  nchan_ = tmpCol( 0 ).nelements() ;
     424  nrow_ = tab.nrow() / npolOrg ;
     425//   cout << "npol_ = " << npol_ << endl ;
     426//   cout << "nchan_ = " << nchan_ << endl ;
     427//   cout << "nrow_ = " << nrow_ << endl ;
     428  spectra.resize( npol_, nchan_, nrow_ ) ;
     429  flagtra.resize( npol_, nchan_, nrow_ ) ;
     430  rflag.resize( npol_, nrow_ ) ;
     431  Cube<Float> tsys( npol_, nchan_, nrow_ ) ;
     432  Matrix<Double> tint( npol_, nrow_ ) ;
     433  for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
    440434    Table subt = tab( tab.col("POLNO") == pollist_[ipol] ) ;
    441435    ROArrayColumn<Float> spectraCol( subt, "SPECTRA" ) ;
    … …  
    443437    ROArrayColumn<uChar> flagtraCol( subt, "FLAGTRA" ) ;
    444438    ROScalarColumn<uInt> rflagCol( subt, "FLAGROW" ) ;
     439    ROArrayColumn<Float> tsysCol( subt, "TSYS" ) ;
     440    ROScalarColumn<Double> tintCol( subt, "INTERVAL" ) ;
    445441    Matrix<Float> tmpF = spectra.yzPlane( ipol ) ;
    446442    Matrix<uChar> tmpUC = flagtra.yzPlane( ipol ) ;
    447443    Vector<uInt> tmpUI = rflag.row( ipol ) ;
    448     cout << "tmpF.shape()=" << tmpF.shape() << endl ;
    449     cout << "tmpUC.shape()=" << tmpUC.shape() << endl ;
    450     cout << "tmpUI.shape()=" << tmpUI.shape() << endl ;
    451444    spectraCol.getColumn( tmpF ) ;
    452445    flagtraCol.getColumn( tmpUC ) ;
    … …  
    454447    if ( ipol == 0 )
    455448      directionCol.getColumn( direction ) ;
    456   }
    457   cout << "getData end" << endl ;
     449    Matrix<Float> tmpF2 = tsysCol.getColumn() ;
     450    Vector<Double> tmpD = tint.row( ipol ) ;
     451    if ( tmpF2.shape()(0) == nchan_ ) {
     452      tsys.yzPlane( ipol ) = tmpF2 ;
     453    }
     454    else {
     455      tsys.yzPlane( ipol ) = tmpF2(0,0) ;
     456    }
     457    tintCol.getColumn( tmpD ) ;
     458  }
     459
     460  getWeight( weight, tsys, tint ) ;
     461}
     462
     463void STGrid::getWeight( Matrix<Float> &w,
     464                        Cube<Float> &tsys,
     465                        Matrix<Double> &tint )
     466{
     467  // resize
     468  w.resize( nchan_, nrow_ ) ;
     469
     470  // set weight
     471  w = 1.0 ;
     472  Bool warn = False ;
     473  if ( wtype_.compare( "UNIFORM" ) == 0 ) {
     474    // do nothing
     475  }
     476  else if ( wtype_.compare( "TINT" ) == 0 ) {
     477    if ( npol_ > 1 ) warn = True ;
     478    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
     479      Float val = mean( tint.column( irow ) ) ;
     480      w.column( irow ) = w.column( irow ) *  val ;
     481    }
     482  }
     483  else if ( wtype_.compare( "TSYS" ) == 0 ) {
     484    if ( npol_ > 1 ) warn = True ;
     485    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
     486      Matrix<Float> arr = tsys.xyPlane( irow ) ;
     487      for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
     488        Float val = mean( arr.column( ichan ) ) ;
     489        w(ichan,irow) = w(ichan,irow) / ( val * val ) ;
     490      }
     491    }
     492  }
     493  else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
     494    if ( npol_ > 1 ) warn = True ;
     495    for ( Int irow = 0 ; irow < nrow_ ; irow++ ) {
     496      Float interval = mean( tint.column( irow ) ) ;
     497      Matrix<Float> arr = tsys.xyPlane( irow ) ;
     498      for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
     499        Float temp = mean( arr.column( ichan ) ) ;
     500        w(ichan,irow) = w(ichan,irow) * interval / ( temp * temp ) ;
     501      }
     502    }
     503  }
     504  else {
     505    LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
     506    os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
     507  }
     508
     509  if ( npol_ > 1 ) {
     510    LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
     511    os << LogIO::WARN << "STGrid doesn't support assigning independent polarization-dependent weight. Use averaged weight over polarization." << LogIO::POST ;
     512  }
    458513}
    459514
    … …  
    612667  Table tab = out->table() ;
    613668  IPosition dshape = data_.shape() ;
    614   Int npol = dshape[2] ;
    615   Int nchan = dshape[3] ;
    616   Int nrow = nx_ * ny_ * npol ;
    617   tab.rwKeywordSet().define( "nPol", npol ) ;
     669  //Int npol = dshape[2] ;
     670  //Int nchan = dshape[3] ;
     671  Int nrow = nx_ * ny_ * npol_ ;
     672  tab.rwKeywordSet().define( "nPol", npol_ ) ;
    618673  tab.addRow( nrow ) ;
    619674  Vector<Double> cpix( 2 ) ;
    … …  
    627682  for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
    628683    for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
    629       for ( Int ipol = 0 ; ipol < npol ; ipol++ ) {
     684      for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
    630685        IPosition start( 4, ix, iy, ipol, 0 ) ;
    631         IPosition end( 4, ix, iy, ipol, nchan-1 ) ;
     686        IPosition end( 4, ix, iy, ipol, nchan_-1 ) ;
    632687        IPosition inc( 4, 1, 1, 1, 1 ) ;
    633688        Vector<Float> sp = data_( start, end, inc ) ;

  • trunk/src/STGrid.h

    r2360 r2361  
    5858  void setOption( string convtype="box",
    5959                  int convsupport=-1 ) ;
     60
     61  void setWeight( const string wType="uniform" ) ;
     62
    6063  void grid() ;
    6164 
    … …  
    7982                Matrix<Double> &direction,
    8083                Cube<uChar> &flagtra,
    81                 Matrix<uInt> &rflag ) ;
     84                Matrix<uInt> &rflag,
     85                Matrix<Float> &weight ) ;
     86
     87  void getWeight( Matrix<Float> &w,
     88                  Cube<Float> &tsys,
     89                  Matrix<Double> &tint ) ;
    8290
    8391  void toInt( Array<uChar> *u, Array<Int> *v ) ;
    … …  
    95103  Int nx_ ;
    96104  Int ny_ ;
     105  Int npol_ ;
     106  Int nchan_ ;
     107  Int nrow_ ;
     108  Int ngrid_ ;
    97109  Double cellx_ ;
    98110  Double celly_ ;
    … …  
    104116  Array<Float> data_ ;
    105117  Vector<uInt> pollist_ ;
     118  String wtype_ ;
    106119
    107120  Table tab_ ;

  • trunk/src/python_STGrid.cpp

    r2360 r2361  
    2121    .def( init <> () )
    2222    .def( init < const std::string > () )
    23     //.def( init < const STGrid& > () )
    24 //     .def("_defineimage", &STGridWrapper::defineImage)
    25 //     .def("_setoption", &STGridWrapper::setOption)
    26 //     .def("_grid", &STGridWrapper::grid)
    2723    .def("_setpollist", &STGrid::setPolList)
    2824    .def("_defineimage", &STGrid::defineImage)
    … …  
    3026    .def("_grid", &STGrid::grid)
    3127    .def("_setin", &STGrid::setFileIn)
    32 //     .def("_setout", &STGrid::setFileOut)
     28    .def("_setweight", &STGrid::setWeight)
    3329    .def("_save", &STGrid::saveData)
    3430    ;
Note: See TracChangeset for help on using the changeset viewer.