- Timestamp:
- 12/28/11 16:51:43 (13 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/asapgrid.py
r2387 r2390 9 9 class asapgrid: 10 10 def __init__( self, infile ): 11 self.infile = infile12 11 self.outfile = None 13 self.gridder = stgrid( self.infile )14 12 self.ifno = None 13 self.gridder = stgrid() 14 self.setData( infile ) 15 15 16 16 def setData( self, infile ): 17 self.gridder._setin( infile ) 17 if isinstance( infile, str ): 18 self.gridder._setin( infile ) 19 else: 20 self.gridder._setfiles( infile ) 21 self.infile = infile 18 22 19 23 def setIF( self, ifno ): … … 58 62 class _SDGridPlotter: 59 63 def __init__( self, infile, outfile=None, ifno=-1 ): 60 self.infile = infile 64 if isinstance( infile, str ): 65 self.infile = [infile] 66 else: 67 self.infile = infile 61 68 self.outfile = outfile 62 69 if self.outfile is None: 63 self.outfile = self.infile .rstrip('/')+'.grid'70 self.outfile = self.infile[0].rstrip('/')+'.grid' 64 71 self.nx = -1 65 72 self.ny = -1 … … 138 145 # plot observed position 139 146 if plotobs: 140 self.createTableIn() 141 irow = 0 142 while ( irow < self.nrow ): 143 chunk = self.getPointingChunk( irow ) 144 #print chunk 145 pl.plot(chunk[0],chunk[1],',',color='green') 146 irow += chunk.shape[1] 147 #print irow 147 for i in xrange(len(self.infile)): 148 self.createTableIn( self.infile[i] ) 149 irow = 0 150 while ( irow < self.nrow ): 151 chunk = self.getPointingChunk( irow ) 152 #print chunk 153 pl.plot(chunk[0],chunk[1],',',color='green') 154 irow += chunk.shape[1] 155 #print irow 148 156 # show image 149 157 extent=[self.blc[0]-0.5*self.cellx, … … 157 165 pl.title( title ) 158 166 159 def createTableIn( self ): 160 self.tablein = scantable( self.infile, average=False ) 167 def createTableIn( self, tab ): 168 del self.tablein 169 self.tablein = scantable( tab, average=False ) 161 170 if self.ifno < 0: 162 171 ifno = self.tablein.getif(0) … … 166 175 sel = selector() 167 176 sel.set_ifs( ifno ) 168 self.tablein.set_selection( sel ) 177 self.tablein.set_selection( sel ) 169 178 self.nchan = len(self.tablein._getspectrum(0)) 170 self.nrow = self.tablein.nrow() 179 self.nrow = self.tablein.nrow() 171 180 del sel 172 181 … … 195 204 else: 196 205 retval = data[pol] 197 #retval[0][self.nx-1] = -1.0198 206 return retval 199 207 -
trunk/src/STGrid.cpp
r2389 r2390 47 47 48 48 setFileIn( infile ) ; 49 } 50 51 STGrid::STGrid( const vector<string> infile ) 52 { 53 init() ; 54 55 setFileList( infile ) ; 49 56 } 50 57 … … 82 89 infileList_.resize( 1 ) ; 83 90 infileList_[0] = String(infile) ; 91 nfile_ = 1 ; 92 } 93 } 94 95 void STGrid::setFileList( const vector<string> infile ) 96 { 97 nfile_ = infile.size() ; 98 infileList_.resize( nfile_ ) ; 99 for ( uInt i = 0 ; i < nfile_ ; i++ ) { 100 infileList_[i] = infile[i] ; 84 101 } 85 102 } … … 343 360 double eInitPol = 0.0 ; 344 361 345 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) { 346 t0 = mathutil::gettimeofday_sec() ; 347 initPol( ipol ) ; 348 t1 = mathutil::gettimeofday_sec() ; 349 eInitPol += t1-t0 ; 350 351 polMap[0] = ipol ; 352 353 os << "start pol " << ipol << LogIO::POST ; 354 355 while( !pastEnd() ) { 356 // data storage 357 IPosition cshape( 3, npol_, nchan_, nchunk_ ) ; 358 IPosition mshape( 2, npol_, nchunk_ ) ; 359 IPosition vshape( 1, nchunk_ ) ; 360 IPosition dshape( 2, 2, nchunk_ ) ; 361 IPosition wshape( 2, nchan_, nchunk_ ) ; 362 Array<Complex> spectra( wshape ) ; 363 Array<Int> flagtra( wshape ) ; 364 Array<Int> rflag( vshape ) ; 365 Array<Float> weight( wshape ) ; 366 Array<Double> direction( dshape ) ; 367 Array<Double> xypos( dshape ) ; 362 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) { 363 initTable( ifile ) ; 364 365 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 366 367 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) { 368 t0 = mathutil::gettimeofday_sec() ; 369 initPol( ipol ) ; 370 t1 = mathutil::gettimeofday_sec() ; 371 eInitPol += t1-t0 ; 368 372 369 spectraF_.resize( wshape ) ; 370 flagtraUC_.resize( wshape ) ; 371 rflagUI_.resize( vshape ) ; 373 polMap[0] = ipol ; 374 375 os << "start pol " << ipol << LogIO::POST ; 376 377 while( !pastEnd() ) { 378 // data storage 379 IPosition cshape( 3, npol_, nchan_, nchunk_ ) ; 380 IPosition mshape( 2, npol_, nchunk_ ) ; 381 IPosition vshape( 1, nchunk_ ) ; 382 IPosition dshape( 2, 2, nchunk_ ) ; 383 IPosition wshape( 2, nchan_, nchunk_ ) ; 384 Array<Complex> spectra( wshape ) ; 385 Array<Int> flagtra( wshape ) ; 386 Array<Int> rflag( vshape ) ; 387 Array<Float> weight( wshape ) ; 388 Array<Double> direction( dshape ) ; 389 Array<Double> xypos( dshape ) ; 390 391 spectraF_.resize( wshape ) ; 392 flagtraUC_.resize( wshape ) ; 393 rflagUI_.resize( vshape ) ; 394 395 // retrieve data 396 t0 = mathutil::gettimeofday_sec() ; 397 Int nrow = getDataChunk( spectra, direction, flagtra, rflag, weight ) ; 398 t1 = mathutil::gettimeofday_sec() ; 399 eGetData += t1-t0 ; 400 401 // world -> pixel 402 t0 = mathutil::gettimeofday_sec() ; 403 toPixel( direction, xypos ) ; 404 t1 = mathutil::gettimeofday_sec() ; 405 eToPixel += t1-t0 ; 406 407 // call ggridsd 408 t0 = mathutil::gettimeofday_sec() ; 409 call_ggridsd( xypos, 410 spectra, 411 nvispol, 412 nchan_, 413 flagtra, 414 rflag, 415 weight, 416 nrow, 417 irow, 418 gdataArrC, 419 gwgtArr, 420 gnx, 421 gny, 422 npol_, 423 nchan_, 424 convSupport_, 425 convSampling_, 426 convFunc, 427 chanMap, 428 polMap ) ; 429 t1 = mathutil::gettimeofday_sec() ; 430 eGGridSD += t1-t0 ; 431 432 } 433 434 os << "end pol " << ipol << LogIO::POST ; 435 436 nprocessed_ = 0 ; 437 } 438 } 439 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; 440 os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ; 441 os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ; 442 os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ; 443 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ; 444 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ; 445 446 delete polMap ; 447 delete chanMap ; 448 449 // set data 450 setData( gdataArrC, gwgtArr ) ; 451 452 } 453 454 void STGrid::gridPerPol() 455 { 456 LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ; 457 double t0, t1 ; 458 459 // convolution kernel 460 Vector<Float> convFunc ; 461 t0 = mathutil::gettimeofday_sec() ; 462 setConvFunc( convFunc ) ; 463 t1 = mathutil::gettimeofday_sec() ; 464 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 465 466 // prepare grid data storage 467 Int gnx = nx_ ; 468 Int gny = ny_ ; 469 // // Extend grid plane with convSupport_ 470 // Int gnx = nx_+convSupport_*2 ; 471 // Int gny = ny_+convSupport_*2 ; 472 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ; 473 Array<Complex> gdataArrC( gshape, 0.0 ) ; 474 //Array<Float> gwgtArr( gshape, 0.0 ) ; 475 // 2011/12/20 TN 476 // data_ and weight array shares storage 477 data_.resize( gshape ) ; 478 data_ = 0.0 ; 479 Array<Float> gwgtArr( data_ ) ; 480 481 // maps 482 Int *chanMap = new Int[nchan_] ; 483 { 484 Int *work_p = chanMap ; 485 for ( Int i = 0 ; i < nchan_ ; i++ ) { 486 *work_p = i ; 487 work_p++ ; 488 } 489 } 490 Int *polMap = new Int[1] ; 491 492 // some parameters for ggridsd 493 Int nvispol = 1 ; 494 Int irow = -1 ; 495 496 // for performance check 497 double eInitPol = 0.0 ; 498 double eGetData = 0.0 ; 499 double eToPixel = 0.0 ; 500 double eGGridSD = 0.0 ; 501 double eToInt = 0.0 ; 502 503 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) { 504 initTable( ifile ) ; 505 506 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 507 // to read data from the table 508 IPosition mshape( 2, nchan_, nrow_ ) ; 509 IPosition dshape( 2, 2, nrow_ ) ; 510 Array<Complex> spectra( mshape ) ; 511 Array<Double> direction( dshape ) ; 512 Array<Int> flagtra( mshape ) ; 513 Array<Int> rflag( IPosition(1,nrow_) ) ; 514 Array<Float> weight( mshape ) ; 515 Array<Double> xypos( dshape ) ; 516 517 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) { 518 t0 = mathutil::gettimeofday_sec() ; 519 initPol( ipol ) ; 520 t1 = mathutil::gettimeofday_sec() ; 521 eInitPol += t1-t0 ; 372 522 373 523 // retrieve data 374 524 t0 = mathutil::gettimeofday_sec() ; 375 Int nrow = getDataChunk( spectra, direction, flagtra, rflag, weight ) ;525 getDataPerPol( spectra, direction, flagtra, rflag, weight ) ; 376 526 t1 = mathutil::gettimeofday_sec() ; 377 527 eGetData += t1-t0 ; … … 379 529 // world -> pixel 380 530 t0 = mathutil::gettimeofday_sec() ; 381 toPixel( direction, xypos ) ; 531 toPixel( direction, xypos ) ; 382 532 t1 = mathutil::gettimeofday_sec() ; 383 533 eToPixel += t1-t0 ; 384 534 385 535 // call ggridsd 536 polMap[0] = ipol ; 386 537 t0 = mathutil::gettimeofday_sec() ; 387 538 call_ggridsd( xypos, … … 392 543 rflag, 393 544 weight, 394 nrow ,545 nrow_, 395 546 irow, 396 547 gdataArrC, … … 407 558 t1 = mathutil::gettimeofday_sec() ; 408 559 eGGridSD += t1-t0 ; 409 410 } 411 412 os << "end pol " << ipol << LogIO::POST ; 413 414 nprocessed_ = 0 ; 560 } 415 561 } 416 562 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; … … 420 566 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ; 421 567 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ; 422 423 delete polMap ;424 delete chanMap ;425 426 // set data427 setData( gdataArrC, gwgtArr ) ;428 429 }430 431 void STGrid::gridPerPol()432 {433 LogIO os( LogOrigin("STGrid", "gridPerPol", WHERE) ) ;434 double t0, t1 ;435 436 // convolution kernel437 Vector<Float> convFunc ;438 t0 = mathutil::gettimeofday_sec() ;439 setConvFunc( convFunc ) ;440 t1 = mathutil::gettimeofday_sec() ;441 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;442 443 // prepare grid data storage444 Int gnx = nx_ ;445 Int gny = ny_ ;446 // // Extend grid plane with convSupport_447 // Int gnx = nx_+convSupport_*2 ;448 // Int gny = ny_+convSupport_*2 ;449 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;450 Array<Complex> gdataArrC( gshape, 0.0 ) ;451 //Array<Float> gwgtArr( gshape, 0.0 ) ;452 // 2011/12/20 TN453 // data_ and weight array shares storage454 data_.resize( gshape ) ;455 data_ = 0.0 ;456 Array<Float> gwgtArr( data_ ) ;457 458 // to read data from the table459 IPosition mshape( 2, nchan_, nrow_ ) ;460 IPosition dshape( 2, 2, nrow_ ) ;461 Array<Complex> spectra( mshape ) ;462 Array<Double> direction( dshape ) ;463 Array<Int> flagtra( mshape ) ;464 Array<Int> rflag( IPosition(1,nrow_) ) ;465 Array<Float> weight( mshape ) ;466 Array<Double> xypos( dshape ) ;467 468 // maps469 Int *chanMap = new Int[nchan_] ;470 {471 Int *work_p = chanMap ;472 for ( Int i = 0 ; i < nchan_ ; i++ ) {473 *work_p = i ;474 work_p++ ;475 }476 }477 Int *polMap = new Int[1] ;478 479 // some parameters for ggridsd480 Int nvispol = 1 ;481 Int irow = -1 ;482 483 // for performance check484 double eInitPol = 0.0 ;485 double eGetData = 0.0 ;486 double eToPixel = 0.0 ;487 double eGGridSD = 0.0 ;488 double eToInt = 0.0 ;489 490 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {491 t0 = mathutil::gettimeofday_sec() ;492 initPol( ipol ) ;493 t1 = mathutil::gettimeofday_sec() ;494 eInitPol += t1-t0 ;495 496 // retrieve data497 t0 = mathutil::gettimeofday_sec() ;498 getDataPerPol( spectra, direction, flagtra, rflag, weight ) ;499 t1 = mathutil::gettimeofday_sec() ;500 eGetData += t1-t0 ;501 502 // world -> pixel503 t0 = mathutil::gettimeofday_sec() ;504 toPixel( direction, xypos ) ;505 t1 = mathutil::gettimeofday_sec() ;506 eToPixel += t1-t0 ;507 508 // call ggridsd509 polMap[0] = ipol ;510 t0 = mathutil::gettimeofday_sec() ;511 call_ggridsd( xypos,512 spectra,513 nvispol,514 nchan_,515 flagtra,516 rflag,517 weight,518 nrow_,519 irow,520 gdataArrC,521 gwgtArr,522 gnx,523 gny,524 npol_,525 nchan_,526 convSupport_,527 convSampling_,528 convFunc,529 chanMap,530 polMap ) ;531 t1 = mathutil::gettimeofday_sec() ;532 eGGridSD += t1-t0 ;533 }534 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;535 os << "getData: elapsed time is " << eGetData-eToInt-eGetWeight << " sec." << LogIO::POST ;536 os << "toPixel: elapsed time is " << eToPixel << " sec." << LogIO::POST ;537 os << "ggridsd: elapsed time is " << eGGridSD << " sec." << LogIO::POST ;538 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;539 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;540 568 541 569 // delete maps … … 560 588 561 589 attach( ptab_ ) ; 590 } 591 592 void STGrid::initTable( uInt idx ) 593 { 594 tab_ = tableList_[idx] ; 595 nrow_ = rows_[idx] ; 562 596 } 563 597 … … 720 754 { 721 755 //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ; 722 ROArrayColumn<Double> dirCol( tab_, "DIRECTION" ) ;723 Matrix<Double> direction = dir Col.getColumn() ;756 directionCol_.attach( tableList_[0], "DIRECTION" ) ; 757 Matrix<Double> direction = directionCol_.getColumn() ; 724 758 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ; 725 759 minMax( xmin, xmax, direction.row( 0 ) ) ; 726 760 minMax( ymin, ymax, direction.row( 1 ) ) ; 761 Double amin, amax, bmin, bmax ; 762 for ( uInt i = 1 ; i < nfile_ ; i++ ) { 763 directionCol_.attach( tableList_[i], "DIRECTION" ) ; 764 direction.assign( directionCol_.getColumn() ) ; 765 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ; 766 minMax( amin, amax, direction.row( 0 ) ) ; 767 minMax( bmin, bmax, direction.row( 1 ) ) ; 768 xmin = min( xmin, amin ) ; 769 xmax = max( xmax, amax ) ; 770 ymin = min( ymin, bmin ) ; 771 ymax = max( ymax, bmax ) ; 772 } 727 773 //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ; 728 774 //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ; … … 733 779 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ; 734 780 Int ifno = ifno_ ; 735 uInt nfile = infileList_.size() ; 736 tableList_.resize( nfile ) ; 737 for ( uInt i = 0 ; i < nfile ; i++ ) { 781 tableList_.resize( nfile_ ) ; 782 for ( uInt i = 0 ; i < nfile_ ; i++ ) { 738 783 Table taborg( infileList_[i] ) ; 739 784 if ( ifno == -1 ) { … … 766 811 } 767 812 } 768 tab_ = tableList_[0] ;769 813 } 770 814 … … 785 829 tsysCol_.attach( tab, "TSYS" ) ; 786 830 intervalCol_.attach( tab, "INTERVAL" ) ; 787 // Vector<String> colnames( 6 ) ;788 // colnames[0] = "SPECTRA" ;789 // colnames[1] = "FLAGTRA" ;790 // colnames[2] = "DIRECTION" ;791 // colnames[3] = "FLAGROW" ;792 // colnames[4] = "TSYS" ;793 // colnames[5] = "INTERVAL" ;794 // row_ = ROTableRow( tab_, colnames ) ;795 // const TableRecord &rec = row_.record() ;796 // spectraRF_.attachToRecord( rec, colnames[0] ) ;797 // flagtraRF_.attachToRecord( rec, colnames[1] ) ;798 // directionRF_.attachToRecord( rec, colnames[2] ) ;799 // flagRowRF_.attachToRecord( rec, colnames[3] ) ;800 // tsysRF_.attachToRecord( rec, colnames[4] ) ;801 // intervalRF_.attachToRecord( rec, colnames[5] ) ;802 831 } 803 832 … … 939 968 { 940 969 LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ; 941 ROScalarColumn<uInt> polnoCol( tab _, "POLNO" ) ;970 ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ; 942 971 Vector<uInt> pols = polnoCol.getColumn() ; 943 972 //os << pols << LogIO::POST ; … … 946 975 uInt polno ; 947 976 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) { 948 //polno = polnoCol( i ) ;949 977 polno = pols( i ) ; 950 978 if ( allNE( pollistOrg, polno ) ) { … … 972 1000 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ; 973 1001 } 974 nrow_ = tab_.nrow() / npolOrg_ ; 975 ROArrayColumn<uChar> tmpCol( tab_, "FLAGTRA" ) ; 976 nchan_ = tmpCol( 0 ).nelements() ; 1002 rows_.resize( nfile_ ) ; 1003 for ( uInt i = 0 ; i < nfile_ ; i++ ) { 1004 rows_[i] = tableList_[i].nrow() / npolOrg_ ; 1005 } 1006 flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ; 1007 nchan_ = flagtraCol_( 0 ).nelements() ; 977 1008 // os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl 978 1009 // << "nchan_ = " << nchan_ << endl … … 1214 1245 t0 = mathutil::gettimeofday_sec() ; 1215 1246 1216 //Int polno = 0 ;1217 1247 String outfile_ ; 1218 1248 if ( outfile.size() == 0 ) { -
trunk/src/STGrid.h
r2389 r2390 41 41 STGrid() ; 42 42 STGrid( const string infile ) ; 43 STGrid( const vector<string> infile ) ; 43 44 virtual ~STGrid() {} ; 44 45 45 46 void setFileIn( const string infile ) ; 47 void setFileList( const vector<string> infile ) ; 46 48 47 49 void setIF( unsigned int ifno ) { ifno_ = ifno ; } ; … … 156 158 157 159 void initPol( Int ipol ) ; 160 void initTable( uInt idx ) ; 158 161 Bool isMultiIF( Table &tab ) ; 159 162 … … 167 170 168 171 Block<String> infileList_ ; 172 uInt nfile_ ; 169 173 Int ifno_ ; 170 174 Int nx_ ; … … 174 178 Int nchan_ ; 175 179 Int nrow_ ; 180 Block<Int> rows_ ; 176 181 Double cellx_ ; 177 182 Double celly_ ; … … 196 201 ROScalarColumn<Double> intervalCol_ ; 197 202 Int nprocessed_ ; 198 Vector<uInt> rows_ ;199 203 Int nchunk_ ; 200 204 -
trunk/src/python_STGrid.cpp
r2364 r2390 21 21 .def( init <> () ) 22 22 .def( init < const std::string > () ) 23 .def( init < const std::vector<std::string> > () ) 23 24 .def("_setif", &STGrid::setIF) 24 25 .def("_setpollist", &STGrid::setPolList) … … 28 29 .def("_grid", &STGrid::grid) 29 30 .def("_setin", &STGrid::setFileIn) 31 .def("_setfiles", &STGrid::setFileList) 30 32 .def("_setweight", &STGrid::setWeight) 31 33 .def("_save", &STGrid::saveData)
Note:
See TracChangeset
for help on using the changeset viewer.