Changeset 2396


Ignore:
Timestamp:
01/11/12 13:51:25 (13 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...

Added min/max clipping functionality.
Introduced control parameter and methods for clipping.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asapgrid.py

    r2391 r2396  
    158158        weightType -- weight type ('uniform', 'tint', 'tsys', 'tintsys')
    159159        """
    160         self.gridder._setweight( weightType )
     160        self.gridder._setweight( weightType )
     161
     162    def enableClip( self ):
     163        """
     164        Enable min/max clipping.
     165
     166        By default, min/max clipping is disabled so that you should
     167        call this method before actual gridding if you want to do
     168        clipping.
     169        """
     170        self.gridder._enableclip()
     171
     172    def disableClip( self ):
     173        """
     174        Disable min/max clipping.
     175        """
     176        self.gridder._disableclip()
    161177
    162178    def grid( self ):

  • trunk/src/STGrid.cpp

    r2393 r2396  
    11#include <iostream>
    22#include <fstream>
     3#include <cfloat>
    34
    45#include <casa/BasicSL/String.h>
    … …  
    8485  cellyUI_ = "" ;
    8586  centerUI_ = "" ;
     87  doclip_ = False ;
    8688}
    8789
    … …  
    251253}
    252254
     255#define NEED_UNDERSCORES
     256#if defined(NEED_UNDERSCORES)
     257#define ggridsd2 ggridsd2_
     258#endif
     259extern "C" {
     260   void ggridsd2(Double*,
     261                 const Complex*,
     262                 Int*,
     263                 Int*,
     264                 Int*,
     265                 const Int*,
     266                 const Int*,
     267                 const Float*,
     268                 Int*,
     269                 Int*,
     270                 Complex*,
     271                 Float*,
     272                 Int*,
     273                 Complex*,
     274                 Float*,
     275                 Float*,
     276                 Complex*,
     277                 Float*,
     278                 Float*,
     279                 Int*,
     280                 Int*,
     281                 Int *,
     282                 Int *,
     283                 Int*,
     284                 Int*,
     285                 Float*,
     286                 Int*,
     287                 Int*,
     288                 Double*);
     289}
     290void STGrid::call_ggridsd2( Array<Double> &xypos,
     291                            Array<Complex> &spectra,
     292                            Int &nvispol,
     293                            Int &nvischan,
     294                            Array<Int> &flagtra,
     295                            Array<Int> &flagrow,
     296                            Array<Float> &weight,
     297                            Int &nrow,
     298                            Int &irow,
     299                            Array<Complex> &gdata,
     300                            Array<Float> &gwgt,
     301                            Array<Int> &npoints,
     302                            Array<Complex> &clipmin,
     303                            Array<Float> &clipwmin,
     304                            Array<Float> &clipcmin,
     305                            Array<Complex> &clipmax,
     306                            Array<Float> &clipwmax,
     307                            Array<Float> &clipcmax,
     308                            Int &nx,
     309                            Int &ny,
     310                            Int &npol,
     311                            Int &nchan,
     312                            Int &support,
     313                            Int &sampling,
     314                            Vector<Float> &convFunc,
     315                            Int *chanMap,
     316                            Int *polMap )
     317{
     318  // parameters for gridding
     319  Int idopsf = 0 ;
     320  Int len = npol*nchan ;
     321  Double *sumw_p = new Double[len] ;
     322  {
     323    Double *work_p = sumw_p ;
     324    for ( Int i = 0 ; i < len ; i++ ) {
     325      *work_p = 0.0 ;
     326      work_p++ ;
     327    }
     328  }
     329
     330  // prepare pointer
     331  Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG, deleteNpts, deleteCMin, deleteCWMin, deleteCCMin, deleteCMax, deleteCWMax, deleteCCMax ;
     332  Double *xy_p = xypos.getStorage( deletePos ) ;
     333  const Complex *values_p = spectra.getStorage( deleteData ) ;
     334  const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
     335  const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
     336  const Float *wgt_p = weight.getStorage( deleteWgt ) ;
     337  Complex *grid_p = gdata.getStorage( deleteDataG ) ;
     338  Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
     339  Float *conv_p = convFunc.getStorage( deleteConv ) ;
     340  Int *npts_p = npoints.getStorage( deleteNpts ) ;
     341  Complex *cmin_p = clipmin.getStorage( deleteCMin ) ;
     342  Float *cwmin_p = clipwmin.getStorage( deleteCWMin ) ;
     343  Float *ccmin_p = clipcmin.getStorage( deleteCCMin ) ;
     344  Complex *cmax_p = clipmax.getStorage( deleteCMax ) ;
     345  Float *cwmax_p = clipwmax.getStorage( deleteCWMax ) ;
     346  Float *ccmax_p = clipcmax.getStorage( deleteCCMax ) ;
     347
     348  // pass copy of irow to ggridsd since it will be modified in theroutine
     349  Int irowCopy = irow ;
     350     
     351  // call ggridsd
     352  ggridsd2( xy_p,
     353            values_p,
     354            &nvispol,
     355            &nvischan,
     356            &idopsf,
     357            flag_p,
     358            rflag_p,
     359            wgt_p,
     360            &nrow,
     361            &irowCopy,
     362            grid_p,
     363            wgrid_p,
     364            npts_p,
     365            cmin_p,
     366            cwmin_p,
     367            ccmin_p,
     368            cmax_p,
     369            cwmax_p,
     370            ccmax_p,
     371            &nx,
     372            &ny,
     373            &npol,
     374            &nchan,
     375            &support,
     376            &sampling,
     377            conv_p,
     378            chanMap,
     379            polMap,
     380            sumw_p ) ;
     381
     382  // finalization
     383  xypos.putStorage( xy_p, deletePos ) ;
     384  spectra.freeStorage( values_p, deleteData ) ;
     385  flagtra.freeStorage( flag_p, deleteFlag ) ;
     386  flagrow.freeStorage( rflag_p, deleteFlagR ) ;
     387  weight.freeStorage( wgt_p, deleteWgt ) ;
     388  gdata.putStorage( grid_p, deleteDataG ) ;
     389  gwgt.putStorage( wgrid_p, deleteWgtG ) ;
     390  convFunc.putStorage( conv_p, deleteConv ) ;
     391  clipmin.putStorage( cmin_p, deleteCMin ) ;
     392  clipwmin.putStorage( cwmin_p, deleteCWMin ) ;
     393  clipcmin.putStorage( ccmin_p, deleteCCMin ) ;
     394  clipmax.putStorage( cmax_p, deleteCMax ) ;
     395  clipwmax.putStorage( cwmax_p, deleteCWMax ) ;
     396  clipcmax.putStorage( ccmax_p, deleteCCMax ) ;
     397  delete sumw_p ;
     398}
    253399
    254400void STGrid::grid()
    … …  
    289435  os << "   weighting = " << wtype_ << endl ;
    290436  os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ;
     437  os << "   doclip = " << (doclip_?"True":"False") << endl ;
    291438  os << "----------" << LogIO::POST ;
    292439  os << LogIO::NORMAL ;
    … …  
    296443  if ( doAll )
    297444    gridPerPol() ;
     445  else if ( doclip_ )
     446    gridPerRowWithClipping() ;
    298447  else
    299448    gridPerRow() ;
     449
    300450}
    301451
    … …  
    338488};
    339489
     490struct STCommonDataWithClipping {
     491  Int gnx;
     492  Int gny;
     493  Int *chanMap;
     494  Vector<Float> convFunc ;
     495  Array<Complex> gdataArrC;
     496  Array<Float> gwgtArr;
     497  Array<Int> npoints ;
     498  Array<Complex> clipMin ;
     499  Array<Float> clipWMin ;
     500  Array<Float> clipCMin ;
     501  Array<Complex> clipMax ;
     502  Array<Float> clipWMax ;
     503  Array<Float> clipCMax ; 
     504  STCommonDataWithClipping(IPosition const &gshape,
     505                           IPosition const &pshape,
     506                           Array<Float> const &data)
     507    : gdataArrC(gshape, 0.0),
     508      gwgtArr(data),
     509      npoints(pshape, 0),
     510      clipMin(gshape, Complex(FLT_MAX,0.0)),
     511      clipWMin(gshape, 0.0),
     512      clipCMin(gshape, 0.0),
     513      clipMax(gshape, Complex(-FLT_MAX,0.0)),
     514      clipWMax(gshape, 0.0),
     515      clipCMax(gshape, 0.0)
     516  {}
     517};
     518
    340519#define DO_AHEAD 3
    341520
    … …  
    346525  const Int pol;
    347526  STContext(STGrid *obj, STCommonData &common, Int pol)
     527    : self(obj), common(common), pol(pol) {}
     528};
     529
     530struct STContextWithClipping {
     531  STCommonDataWithClipping &common;
     532  FIFO<STGChunk *, DO_AHEAD> queue;
     533  STGrid *const self;
     534  const Int pol;
     535  STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
    348536    : self(obj), common(common), pol(pol) {}
    349537};
    … …  
    512700  setData( common.gdataArrC, common.gwgtArr ) ;
    513701
     702}
     703
     704void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
     705{
     706  STContextWithClipping &context = *(STContextWithClipping *)ctx;
     707  STGChunk *chunk = NULL;
     708  try {
     709    context.queue.lock();
     710    chunk = context.queue.get();
     711    context.queue.unlock();
     712  } catch (FullException &e) {
     713    context.queue.unlock();
     714    // TODO: log error
     715    throw PCException();
     716  }
     717
     718  double t0, t1 ;
     719  // world -> pixel
     720  Array<Double> xypos( context.self->dshape_ ) ;
     721  t0 = mathutil::gettimeofday_sec() ;
     722  context.self->toPixel( chunk->direction, xypos ) ;
     723  t1 = mathutil::gettimeofday_sec() ;
     724  context.self->eToPixel_ += t1-t0 ;
     725   
     726  // call ggridsd
     727  Int nvispol = 1 ;
     728  Int irow = -1 ;
     729  t0 = mathutil::gettimeofday_sec() ;
     730  context.self->call_ggridsd2( xypos,
     731                chunk->spectra,
     732                nvispol,
     733                context.self->nchan_,
     734                chunk->flagtra,
     735                chunk->rflag,
     736                chunk->weight,
     737                chunk->nrow,
     738                irow,
     739                context.common.gdataArrC,
     740                context.common.gwgtArr,
     741                context.common.npoints,
     742                context.common.clipMin,
     743                context.common.clipWMin,
     744                context.common.clipCMin,
     745                context.common.clipMax,
     746                context.common.clipWMax,
     747                context.common.clipCMax,
     748                context.common.gnx,
     749                context.common.gny,
     750                context.self->npol_,
     751                context.self->nchan_,
     752                context.self->convSupport_,
     753                context.self->convSampling_,
     754                context.common.convFunc,
     755                context.common.chanMap,
     756                (Int*)&context.pol ) ;
     757  t1 = mathutil::gettimeofday_sec() ;
     758  context.self->eGGridSD_ += t1-t0 ;
     759 
     760  delete chunk;
     761}
     762
     763void STGrid::gridPerRowWithClipping()
     764{
     765  LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
     766  double t0, t1 ;
     767
     768
     769  // grid data
     770  // Extend grid plane with convSupport_
     771  //   Int gnx = nx_+convSupport_*2 ;
     772  //   Int gny = ny_+convSupport_*2 ;
     773  Int gnx = nx_;
     774  Int gny = ny_;
     775
     776  IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
     777  IPosition pshape( 3, gnx, gny, npol_ ) ;
     778  // 2011/12/20 TN
     779  // data_ and gwgtArr share storage
     780  data_.resize( gshape ) ;
     781  data_ = 0.0 ;
     782  //STCommonData common = STCommonData(gshape, data_);
     783  STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
     784                                                              pshape,
     785                                                              data_ ) ;
     786  common.gnx = gnx ;
     787  common.gny = gny ;
     788
     789  // parameters for gridding
     790  Int *chanMap = new Int[nchan_] ;
     791  for ( Int i = 0 ; i < nchan_ ; i++ ) {
     792    chanMap[i] = i ;
     793  }
     794  common.chanMap = chanMap;
     795
     796  // convolution kernel
     797  t0 = mathutil::gettimeofday_sec() ;
     798  setConvFunc( common.convFunc ) ;
     799  t1 = mathutil::gettimeofday_sec() ;
     800  os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
     801
     802  // for performance check
     803  eGetData_ = 0.0 ;
     804  eToPixel_ = 0.0 ;
     805  eGGridSD_ = 0.0 ;
     806  double eInitPol = 0.0 ;
     807
     808  //Broker broker = Broker(produceChunk, consumeChunk);
     809  Broker broker = Broker(produceChunk, consumeChunkWithClipping);
     810  for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
     811    initTable( ifile ) ;
     812
     813    os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;   
     814    for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
     815      t0 = mathutil::gettimeofday_sec() ;
     816      initPol( ipol ) ; // set ptab_ and attach()
     817      t1 = mathutil::gettimeofday_sec() ;
     818      eInitPol += t1-t0 ;
     819     
     820      //STContext context(this, common, ipol);
     821      STContextWithClipping context(this, common, ipol);
     822     
     823      os << "start pol " << ipol << LogIO::POST ;
     824     
     825      nprocessed_ = 0 ;
     826#if 1
     827      broker.runProducerAsMasterThread(&context, DO_AHEAD);
     828#else
     829      for (;;) {
     830        bool produced = produceChunk(&context);
     831        if (! produced) {
     832          break;
     833        }
     834        //consumeChunk(&context);
     835        consumeChunkWithClipping(&context);
     836      }
     837#endif
     838
     839      os << "end pol " << ipol << LogIO::POST ;
     840
     841    }
     842    os << "end table " << ifile << LogIO::POST ;   
     843  }
     844  os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
     845  os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
     846  os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
     847  os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
     848  os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
     849  os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
     850 
     851  delete chanMap ;
     852
     853  // clip min and max in each grid
     854//   os << "BEFORE CLIPPING" << LogIO::POST ;
     855//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
     856//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
     857  t0 = mathutil::gettimeofday_sec() ;
     858  clipMinMax( common.gdataArrC,
     859              common.gwgtArr,
     860              common.npoints,
     861              common.clipMin,
     862              common.clipWMin,
     863              common.clipCMin,
     864              common.clipMax,
     865              common.clipWMax,
     866              common.clipCMax ) ;
     867  t1 = mathutil::gettimeofday_sec() ;
     868  os << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
     869//   os << "AFTER CLIPPING" << LogIO::POST ;
     870//   os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
     871//   os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
     872
     873  // set data
     874  setData( common.gdataArrC, common.gwgtArr ) ;
     875
     876}
     877
     878void STGrid::clipMinMax( Array<Complex> &grid,
     879                         Array<Float> &weight,
     880                         Array<Int> &npoints,
     881                         Array<Complex> &clipmin,
     882                         Array<Float> &clipwmin,
     883                         Array<Float> &clipcmin,
     884                         Array<Complex> &clipmax,
     885                         Array<Float> &clipwmax,
     886                         Array<Float> &clipcmax )
     887{
     888  //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
     889
     890  // prepare pointers
     891  Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
     892  Complex *grid_p = grid.getStorage( delG ) ;
     893  Float *wgt_p = weight.getStorage( delW ) ;
     894  const Int *npts_p = npoints.getStorage( delNP ) ;
     895  const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
     896  const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
     897  const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
     898  const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
     899  const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
     900  const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
     901
     902  const IPosition &gshape = grid.shape() ;
     903  long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
     904  Int nchan = gshape[3] ;
     905  long origin = nchan * offset ;
     906  for ( long i = 0 ; i < offset ; i++ ) {
     907    if ( *npts_p > 2 ) {
     908      for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
     909        // clip minimum and maximum
     910        *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
     911          + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
     912        *wgt_p -= (*cwmin_p)*(*ccmin_p)
     913          + (*cwmax_p)*(*ccmax_p) ;
     914       
     915        grid_p += offset ;
     916        wgt_p += offset ;
     917        cmin_p += offset ;
     918        cwmin_p += offset ;
     919        ccmin_p += offset ;
     920        cmax_p += offset ;
     921        cwmax_p += offset ;
     922        ccmax_p += offset ;
     923      }
     924      grid_p -= origin ;
     925      wgt_p -= origin ;
     926      cmin_p -= origin ;
     927      cwmin_p -= origin ;
     928      ccmin_p -= origin ;
     929      cmax_p -= origin ;
     930      cwmax_p -= origin ;
     931      ccmax_p -= origin ;
     932    }
     933    grid_p++ ;
     934    wgt_p++ ;
     935    npts_p++ ;
     936    cmin_p++ ;
     937    cwmin_p++ ;
     938    ccmin_p++ ;
     939    cmax_p++ ;
     940    cwmax_p++ ;
     941    ccmax_p++ ;
     942  }
     943  grid_p -= offset ;
     944  wgt_p -= offset ;
     945  npts_p -= offset ;
     946  cmin_p -= offset ;
     947  cwmin_p -= offset ;
     948  ccmin_p -= offset ;
     949  cmax_p -= offset ;
     950  cwmax_p -= offset ;
     951  ccmax_p -= offset ; 
     952
     953  // finalization
     954  grid.putStorage( grid_p, delG ) ;
     955  weight.putStorage( wgt_p, delW ) ;
     956  npoints.freeStorage( npts_p, delNP ) ;
     957  clipmin.freeStorage( cmin_p, delCMin ) ;
     958  clipwmin.freeStorage( cwmin_p, delCWMin ) ;
     959  clipcmin.freeStorage( ccmin_p, delCCMin ) ;
     960  clipmax.freeStorage( cmax_p, delCMax ) ;
     961  clipwmax.freeStorage( cwmax_p, delCWMax ) ;
     962  clipcmax.freeStorage( ccmax_p, delCCMax ) ;
    514963}
    515964

  • trunk/src/STGrid.h

    r2393 r2396  
    6565  void setWeight( const string wType="uniform" ) ;
    6666
     67  void enableClip() { doclip_ = True ; } ;
     68  void disableClip() { doclip_ = False ; } ;
     69
    6770  void grid() ;
    6871 
    … …  
    7477  // actual gridding
    7578  void gridPerRow() ;
     79  void gridPerRowWithClipping() ;
    7680  void gridPerPol() ;
     81
     82  // clipping
     83  void clipMinMax( Array<Complex> &data,
     84                   Array<Float> &weight,
     85                   Array<Int> &npoints,
     86                   Array<Complex> &clipmin,
     87                   Array<Float> &clipwmin,
     88                   Array<Float> &clipcmin,
     89                   Array<Complex> &clipmax,
     90                   Array<Float> &clipwmax,
     91                   Array<Float> &clipcmax ) ;
     92                   
    7793
    7894  void setupGrid() ;
    … …  
    166182                     Int *chanMap,
    167183                     Int *polMap ) ;
     184  void call_ggridsd2( Array<Double> &xy,
     185                      Array<Complex> &values,
     186                      Int &nvispol,
     187                      Int &nvischan,
     188                      Array<Int> &flag,
     189                      Array<Int> &rflag,
     190                      Array<Float> &weight,
     191                      Int &nrow,
     192                      Int &irow,
     193                      Array<Complex> &grid,
     194                      Array<Float> &wgrid,
     195                      Array<Int> &npoints,
     196                      Array<Complex> &clipmin,
     197                      Array<Float> &clipwmin,
     198                      Array<Float> &clipcmin,
     199                      Array<Complex> &clipmax,
     200                      Array<Float> &clipwmax,
     201                      Array<Float> &clipcmax,
     202                      Int &nx,
     203                      Int &ny,
     204                      Int &npol,
     205                      Int &nchan,
     206                      Int &support,
     207                      Int &sampling,
     208                      Vector<Float> &convFunc,
     209                      Int *chanMap,
     210                      Int *polMap ) ;
    168211
    169212  void initPol( Int ipol ) ;
    … …  
    172215  static bool produceChunk(void *ctx) throw(concurrent::PCException);
    173216  static void consumeChunk(void *ctx) throw(concurrent::PCException);
     217  static void consumeChunkWithClipping(void *ctx) throw(concurrent::PCException);
    174218
    175219
    … …  
    184228  uInt nfile_ ;
    185229  Int ifno_ ;
     230  Bool doclip_ ;
    186231
    187232  Int nx_ ;

  • trunk/src/python_STGrid.cpp

    r2390 r2396  
    3131    .def("_setfiles", &STGrid::setFileList)
    3232    .def("_setweight", &STGrid::setWeight)
     33    .def("_enableclip", &STGrid::enableClip)
     34    .def("_disableclip", &STGrid::disableClip)
    3335    .def("_save", &STGrid::saveData)
    3436    ;
Note: See TracChangeset for help on using the changeset viewer.