- Timestamp:
- 01/05/12 12:10:37 (13 years ago)
- Location:
- trunk/src
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/CMakeLists.txt
r2370 r2393 21 21 # source files for libpyrap 22 22 set( ASAP_SRCS 23 ${SRCDIR}/concurrent.cpp 23 24 ${SRCDIR}/MathUtils.cpp 24 25 ${SRCDIR}/TableTraverse.cpp -
trunk/src/STGrid.cpp
r2390 r2393 27 27 28 28 using namespace std ; 29 using namespace concurrent ; 29 30 using namespace casa ; 30 31 using namespace asap ; … … 38 39 // constructor 39 40 STGrid::STGrid() 41 : vshape_( 1 ), wshape_( 2 ), dshape_( 2 ) 40 42 { 41 43 init() ; … … 43 45 44 46 STGrid::STGrid( const string infile ) 47 : vshape_( 1 ), wshape_( 2 ), dshape_( 2 ) 45 48 { 46 49 init() ; … … 85 88 void STGrid::setFileIn( const string infile ) 86 89 { 90 nfile_ = 1 ; 87 91 String name( infile ) ; 88 if ( infileList_.size() == 0 || infileList_[0].compare( name ) != 0 ) { 89 infileList_.resize( 1 ) ; 90 infileList_[0] = String(infile) ; 91 nfile_ = 1 ; 92 } 92 infileList_.resize( nfile_ ) ; 93 infileList_[0] = String(infile) ; 93 94 } 94 95 … … 250 251 } 251 252 252 Bool STGrid::pastEnd()253 {254 LogIO os( LogOrigin("STGrid","pastEnd",WHERE) ) ;255 Bool b = nprocessed_ >= nrow_ ;256 return b ;257 }258 253 259 254 void STGrid::grid() … … 312 307 Bool b = nchunk_ >= nrow_ ; 313 308 nchunk_ = min( nchunk_, nrow_ ) ; 309 vshape_ = IPosition( 1, nchunk_ ) ; 310 wshape_ = IPosition( 2, nchan_, nchunk_ ) ; 311 dshape_ = IPosition( 2, 2, nchunk_ ) ; 314 312 return b ; 313 } 314 315 struct STGChunk { 316 Int nrow ; 317 Array<Complex> spectra; 318 Array<Int> flagtra; 319 Array<Int> rflag; 320 Array<Float> weight; 321 Array<Double> direction; 322 STGChunk(IPosition const &wshape, IPosition const &vshape, 323 IPosition const &dshape) 324 : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape), 325 direction(dshape) 326 { } 327 }; 328 329 struct STCommonData { 330 Int gnx; 331 Int gny; 332 Int *chanMap; 333 Vector<Float> convFunc ; 334 Array<Complex> gdataArrC; 335 Array<Float> gwgtArr; 336 STCommonData(IPosition const &gshape, Array<Float> const &data) 337 : gdataArrC(gshape, 0.0), gwgtArr(data) {} 338 }; 339 340 #define DO_AHEAD 3 341 342 struct STContext { 343 STCommonData &common; 344 FIFO<STGChunk *, DO_AHEAD> queue; 345 STGrid *const self; 346 const Int pol; 347 STContext(STGrid *obj, STCommonData &common, Int pol) 348 : self(obj), common(common), pol(pol) {} 349 }; 350 351 352 bool STGrid::produceChunk(void *ctx) throw(PCException) 353 { 354 STContext &context = *(STContext *)ctx; 355 if ( context.self->nprocessed_ >= context.self->nrow_ ) { 356 return false; 357 } 358 STGChunk *chunk = new STGChunk(context.self->wshape_, 359 context.self->vshape_, 360 context.self->dshape_); 361 362 double t0 = mathutil::gettimeofday_sec() ; 363 chunk->nrow = context.self->getDataChunk( 364 context.self->wshape_, context.self->vshape_, context.self->dshape_, 365 chunk->spectra, chunk->direction, 366 chunk->flagtra, chunk->rflag, chunk->weight); 367 double t1 = mathutil::gettimeofday_sec() ; 368 context.self->eGetData_ += t1-t0 ; 369 370 context.queue.lock(); 371 context.queue.put(chunk); 372 context.queue.unlock(); 373 return true; 374 } 375 376 void STGrid::consumeChunk(void *ctx) throw(PCException) 377 { 378 STContext &context = *(STContext *)ctx; 379 STGChunk *chunk = NULL; 380 try { 381 context.queue.lock(); 382 chunk = context.queue.get(); 383 context.queue.unlock(); 384 } catch (FullException &e) { 385 context.queue.unlock(); 386 // TODO: log error 387 throw PCException(); 388 } 389 390 double t0, t1 ; 391 // world -> pixel 392 Array<Double> xypos( context.self->dshape_ ) ; 393 t0 = mathutil::gettimeofday_sec() ; 394 context.self->toPixel( chunk->direction, xypos ) ; 395 t1 = mathutil::gettimeofday_sec() ; 396 context.self->eToPixel_ += t1-t0 ; 397 398 // call ggridsd 399 Int nvispol = 1 ; 400 Int irow = -1 ; 401 t0 = mathutil::gettimeofday_sec() ; 402 context.self->call_ggridsd( xypos, 403 chunk->spectra, 404 nvispol, 405 context.self->nchan_, 406 chunk->flagtra, 407 chunk->rflag, 408 chunk->weight, 409 chunk->nrow, 410 irow, 411 context.common.gdataArrC, 412 context.common.gwgtArr, 413 context.common.gnx, 414 context.common.gny, 415 context.self->npol_, 416 context.self->nchan_, 417 context.self->convSupport_, 418 context.self->convSampling_, 419 context.common.convFunc, 420 context.common.chanMap, 421 (Int*)&context.pol ) ; 422 t1 = mathutil::gettimeofday_sec() ; 423 context.self->eGGridSD_ += t1-t0 ; 424 425 delete chunk; 315 426 } 316 427 … … 320 431 double t0, t1 ; 321 432 322 // convolution kernel323 Vector<Float> convFunc ;324 t0 = mathutil::gettimeofday_sec() ;325 setConvFunc( convFunc ) ;326 t1 = mathutil::gettimeofday_sec() ;327 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;328 433 329 434 // grid data 330 Int gnx = nx_ ;331 Int gny = ny_ ;332 435 // Extend grid plane with convSupport_ 333 // Int gnx = nx_+convSupport_*2 ; 334 // Int gny = ny_+convSupport_*2 ; 436 // Int gnx = nx_+convSupport_*2 ; 437 // Int gny = ny_+convSupport_*2 ; 438 Int gnx = nx_; 439 Int gny = ny_; 440 335 441 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ; 336 Array<Complex> gdataArrC( gshape, 0.0 ) ;337 442 // 2011/12/20 TN 338 443 // data_ and gwgtArr share storage 339 444 data_.resize( gshape ) ; 340 445 data_ = 0.0 ; 341 Array<Float> gwgtArr( data_ ) ; 446 STCommonData common = STCommonData(gshape, data_); 447 common.gnx = gnx ; 448 common.gny = gny ; 342 449 343 450 // parameters for gridding 344 451 Int *chanMap = new Int[nchan_] ; 345 {346 Int *work_p = chanMap;347 for ( Int i = 0 ; i < nchan_ ; i++ ) {348 *work_p = i;349 work_p++ ; 350 }351 }352 Int *polMap = new Int[1];353 Int nvispol = 1;354 Int irow = -1 ;452 for ( Int i = 0 ; i < nchan_ ; i++ ) { 453 chanMap[i] = i ; 454 } 455 common.chanMap = chanMap; 456 457 // convolution kernel 458 t0 = mathutil::gettimeofday_sec() ; 459 setConvFunc( common.convFunc ) ; 460 t1 = mathutil::gettimeofday_sec() ; 461 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ; 355 462 356 463 // for performance check 357 double eGetData= 0.0 ;358 double eToPixel= 0.0 ;359 double eGGridSD= 0.0 ;464 eGetData_ = 0.0 ; 465 eToPixel_ = 0.0 ; 466 eGGridSD_ = 0.0 ; 360 467 double eInitPol = 0.0 ; 361 468 469 Broker broker = Broker(produceChunk, consumeChunk); 362 470 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) { 363 471 initTable( ifile ) ; 364 472 365 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 366 473 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 367 474 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) { 368 475 t0 = mathutil::gettimeofday_sec() ; 369 initPol( ipol ) ; 476 initPol( ipol ) ; // set ptab_ and attach() 370 477 t1 = mathutil::gettimeofday_sec() ; 371 478 eInitPol += t1-t0 ; 372 479 373 polMap[0] = ipol;480 STContext context(this, common, ipol); 374 481 375 482 os << "start pol " << ipol << LogIO::POST ; 376 483 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 484 nprocessed_ = 0 ; 485 #if 1 486 broker.runProducerAsMasterThread(&context, DO_AHEAD); 487 #else 488 for (;;) { 489 bool produced = produceChunk(&context); 490 if (! produced) { 491 break; 492 } 493 consumeChunk(&context); 432 494 } 433 495 #endif 496 434 497 os << "end pol " << ipol << LogIO::POST ; 435 436 nprocessed_ = 0 ;437 }498 499 } 500 os << "end table " << ifile << LogIO::POST ; 438 501 } 439 502 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 ;503 os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ; 504 os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ; 505 os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ; 443 506 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ; 444 507 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ; 445 508 446 delete polMap ;447 509 delete chanMap ; 448 510 449 511 // set data 450 setData( gdataArrC,gwgtArr ) ;512 setData( common.gdataArrC, common.gwgtArr ) ; 451 513 452 514 } … … 488 550 } 489 551 } 490 Int *polMap = new Int[1] ;552 Int polMap[1] ; 491 553 492 554 // some parameters for ggridsd … … 504 566 initTable( ifile ) ; 505 567 506 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 568 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ; 507 569 // to read data from the table 508 570 IPosition mshape( 2, nchan_, nrow_ ) ; … … 520 582 t1 = mathutil::gettimeofday_sec() ; 521 583 eInitPol += t1-t0 ; 522 584 585 os << "start pol " << ipol << LogIO::POST ; 586 523 587 // retrieve data 524 588 t0 = mathutil::gettimeofday_sec() ; … … 558 622 t1 = mathutil::gettimeofday_sec() ; 559 623 eGGridSD += t1-t0 ; 560 } 624 625 os << "end pol " << ipol << LogIO::POST ; 626 627 } 628 os << "end table " << ifile << LogIO::POST ; 561 629 } 562 630 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ; … … 568 636 569 637 // delete maps 570 delete polMap ;571 638 delete chanMap ; 572 639 … … 780 847 Int ifno = ifno_ ; 781 848 tableList_.resize( nfile_ ) ; 849 if ( ifno == -1 ) { 850 Table taborg( infileList_[0] ) ; 851 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ; 852 ifno = ifnoCol( 0 ) ; 853 os << LogIO::WARN 854 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ; 855 } 782 856 for ( uInt i = 0 ; i < nfile_ ; i++ ) { 783 857 Table taborg( infileList_[i] ) ; 784 if ( ifno == -1 ) {785 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;786 ifno = ifnoCol( 0 ) ;787 os << LogIO::WARN788 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;789 }790 858 TableExprNode node ; 791 859 if ( isMultiIF( taborg ) ) { … … 803 871 tableList_[i] = taborg( node ) ; 804 872 } 873 os << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ; 805 874 if ( tableList_[i].nrow() == 0 ) { 806 875 os << LogIO::SEVERE … … 890 959 } 891 960 892 Int STGrid::getDataChunk( Array<Complex> &spectra, 893 Array<Double> &direction, 894 Array<Int> &flagtra, 895 Array<Int> &rflag, 896 Array<Float> &weight ) 961 Int STGrid::getDataChunk( 962 IPosition const &wshape, 963 IPosition const &vshape, 964 IPosition const &dshape, 965 Array<Complex> &spectra, 966 Array<Double> &direction, 967 Array<Int> &flagtra, 968 Array<Int> &rflag, 969 Array<Float> &weight ) 897 970 { 898 971 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ; 972 973 Array<Float> spectraF_(wshape); 974 Array<uChar> flagtraUC_(wshape); 975 Array<uInt> rflagUI_(vshape); 899 976 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ; 900 977 if ( nrow < nchunk_ ) { … … 913 990 return nrow ; 914 991 } 992 993 #if 0 994 Int STGrid::getDataChunk( Array<Complex> &spectra, 995 Array<Double> &direction, 996 Array<Int> &flagtra, 997 Array<Int> &rflag, 998 Array<Float> &weight ) 999 { 1000 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ; 1001 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ; 1002 if ( nrow < nchunk_ ) { 1003 spectra.resize( spectraF_.shape() ) ; 1004 flagtra.resize( flagtraUC_.shape() ) ; 1005 rflag.resize( rflagUI_.shape() ) ; 1006 } 1007 double t0, t1 ; 1008 t0 = mathutil::gettimeofday_sec() ; 1009 convertArray( spectra, spectraF_ ) ; 1010 toInt( flagtraUC_, flagtra ) ; 1011 toInt( rflagUI_, rflag ) ; 1012 t1 = mathutil::gettimeofday_sec() ; 1013 eToInt = t1 - t0 ; 1014 1015 return nrow ; 1016 } 1017 #endif 915 1018 916 1019 Int STGrid::getDataChunk( Array<Float> &spectra, … … 975 1078 uInt polno ; 976 1079 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) { 1080 //polno = polnoCol( i ) ; 977 1081 polno = pols( i ) ; 978 1082 if ( allNE( pollistOrg, polno ) ) { … … 1003 1107 for ( uInt i = 0 ; i < nfile_ ; i++ ) { 1004 1108 rows_[i] = tableList_[i].nrow() / npolOrg_ ; 1109 if ( nrow_ < rows_[i] ) 1110 nrow_ = rows_[i] ; 1005 1111 } 1006 1112 flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ; … … 1245 1351 t0 = mathutil::gettimeofday_sec() ; 1246 1352 1353 //Int polno = 0 ; 1247 1354 String outfile_ ; 1248 1355 if ( outfile.size() == 0 ) { -
trunk/src/STGrid.h
r2390 r2393 26 26 #include <tables/Tables/ScalarColumn.h> 27 27 #include <tables/Tables/ArrayColumn.h> 28 29 #include "concurrent.h" 28 30 //#include <tables/Tables/TableRow.h> 29 31 … … 100 102 Array<uInt> &rflag, 101 103 Array<Float> &weight ) ; 104 Int getDataChunk( IPosition const &wshape, 105 IPosition const &vshape, 106 IPosition const &dshape, 107 Array<Complex> &spectra, 108 Array<Double> &direction, 109 Array<Int> &flagtra, 110 Array<Int> &rflag, 111 Array<Float> &weight ) ; 102 112 Int getDataChunk( Array<Complex> &spectra, 103 113 Array<Double> &direction, … … 128 138 void prepareTable( Table &tab, String &name ) ; 129 139 130 Bool pastEnd() ;140 // Bool pastEnd() ; 131 141 132 142 void selectData() ; … … 160 170 void initTable( uInt idx ) ; 161 171 Bool isMultiIF( Table &tab ) ; 172 static bool produceChunk(void *ctx) throw(concurrent::PCException); 173 static void consumeChunk(void *ctx) throw(concurrent::PCException); 162 174 163 175 … … 172 184 uInt nfile_ ; 173 185 Int ifno_ ; 186 174 187 Int nx_ ; 175 188 Int ny_ ; … … 177 190 Int npolOrg_ ; 178 191 Int nchan_ ; 179 Int nrow_ ;180 Block<Int> rows_ ;181 192 Double cellx_ ; 182 193 Double celly_ ; … … 186 197 Int userSupport_ ; 187 198 Int convSampling_ ; 188 Array<Float> data_ ;189 199 Vector<uInt> pollist_ ; 190 200 Vector<uInt> scanlist_ ; 191 201 String wtype_ ; 202 Block<Table> tableList_ ; 203 Vector<uInt> rows_ ; 204 Int nchunk_ ; 205 206 /////////////// gridPerRow variable 207 IPosition vshape_; 208 IPosition wshape_; 209 IPosition dshape_; 210 // loop variable 211 Int nrow_ ; 212 Array<Float> data_ ; 192 213 193 214 Table tab_ ; 194 Block<Table> tableList_ ;215 // per pol 195 216 Table ptab_ ; 196 217 ROArrayColumn<Float> spectraCol_ ; … … 200 221 ROArrayColumn<Float> tsysCol_ ; 201 222 ROScalarColumn<Double> intervalCol_ ; 223 202 224 Int nprocessed_ ; 203 Int nchunk_ ; 204 205 Array<Float> spectraF_ ; 206 Array<uChar> flagtraUC_ ; 207 Array<uInt> rflagUI_ ; 208 225 226 227 double eGetData_; 228 double eToPixel_; 229 double eGGridSD_; 209 230 }; 210 231 }
Note:
See TracChangeset
for help on using the changeset viewer.