source: branches/hpc33/src/STGrid.cpp@ 2457

Last change on this file since 2457 was 2398, checked in by Takeshi Nakazato, 13 years ago

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...

Reorganization of procedure:

  • Dropped STGrid::gridPerPol
  • Broker is created per input table
  • STGrid::examine is redefined as STGrid::updateChunkShape
File size: 46.5 KB
Line 
1#include <iostream>
2#include <fstream>
3#include <cfloat>
4
5#include <casa/BasicSL/String.h>
6#include <casa/Arrays/Vector.h>
7#include <casa/Arrays/Matrix.h>
8#include <casa/Arrays/Cube.h>
9#include <casa/Arrays/ArrayMath.h>
10#include <casa/Arrays/ArrayIter.h>
11#include <casa/Quanta/Quantum.h>
12#include <casa/Quanta/QuantumHolder.h>
13#include <casa/Utilities/CountedPtr.h>
14#include <casa/Logging/LogIO.h>
15
16#include <tables/Tables/Table.h>
17#include <tables/Tables/TableRecord.h>
18#include <tables/Tables/ExprNode.h>
19#include <tables/Tables/ScalarColumn.h>
20#include <tables/Tables/ArrayColumn.h>
21#include <tables/Tables/TableIter.h>
22
23#include <measures/Measures/MDirection.h>
24
25#include <MathUtils.h>
26
27#include "STGrid.h"
28
29using namespace std ;
30using namespace concurrent ;
31using namespace casa ;
32using namespace asap ;
33
34namespace asap {
35
36// for performance check
37double eToInt = 0.0 ;
38double eGetWeight = 0.0 ;
39
40// constructor
41STGrid::STGrid()
42 : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
43{
44 init() ;
45}
46
47STGrid::STGrid( const string infile )
48 : vshape_( 1 ), wshape_( 2 ), dshape_( 2 )
49{
50 init() ;
51
52 setFileIn( infile ) ;
53}
54
55STGrid::STGrid( const vector<string> infile )
56{
57 init() ;
58
59 setFileList( infile ) ;
60}
61
62void STGrid::init()
63{
64 ifno_ = -1 ;
65 nx_ = -1 ;
66 ny_ = -1 ;
67 npol_ = 0 ;
68 nchan_ = 0 ;
69 nrow_ = 0 ;
70 cellx_ = 0.0 ;
71 celly_ = 0.0 ;
72 center_ = Vector<Double> ( 2, 0.0 ) ;
73 convType_ = "BOX" ;
74 wtype_ = "UNIFORM" ;
75 convSupport_ = -1 ;
76 userSupport_ = -1 ;
77 convSampling_ = 100 ;
78 nprocessed_ = 0 ;
79 nchunk_ = 0 ;
80
81 // initialize user input
82 nxUI_ = -1 ;
83 nyUI_ = -1 ;
84 cellxUI_ = "" ;
85 cellyUI_ = "" ;
86 centerUI_ = "" ;
87 doclip_ = False ;
88}
89
90void STGrid::setFileIn( const string infile )
91{
92 nfile_ = 1 ;
93 String name( infile ) ;
94 infileList_.resize( nfile_ ) ;
95 infileList_[0] = String(infile) ;
96}
97
98void STGrid::setFileList( const vector<string> infile )
99{
100 nfile_ = infile.size() ;
101 infileList_.resize( nfile_ ) ;
102 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
103 infileList_[i] = infile[i] ;
104 }
105}
106
107void STGrid::setPolList( vector<unsigned int> pols )
108{
109 pollist_.assign( Vector<uInt>( pols ) ) ;
110}
111
112void STGrid::setScanList( vector<unsigned int> scans )
113{
114 scanlist_.assign( Vector<uInt>( scans ) ) ;
115}
116
117void STGrid::setWeight( const string wType )
118{
119 wtype_ = String( wType ) ;
120 wtype_.upcase() ;
121}
122
123void STGrid::defineImage( int nx,
124 int ny,
125 string scellx,
126 string scelly,
127 string scenter )
128{
129 nxUI_ = (Int)nx ;
130 nyUI_ = (Int)ny ;
131 cellxUI_ = String( scellx ) ;
132 cellyUI_ = String( scelly ) ;
133 centerUI_ = String( scenter ) ;
134}
135
136void STGrid::setFunc( string convType,
137 int convSupport )
138{
139 convType_ = String( convType ) ;
140 convType_.upcase() ;
141 userSupport_ = (Int)convSupport ;
142}
143
144#define NEED_UNDERSCORES
145#if defined(NEED_UNDERSCORES)
146#define ggridsd ggridsd_
147#endif
148extern "C" {
149 void ggridsd(Double*,
150 const Complex*,
151 Int*,
152 Int*,
153 Int*,
154 const Int*,
155 const Int*,
156 const Float*,
157 Int*,
158 Int*,
159 Complex*,
160 Float*,
161 Int*,
162 Int*,
163 Int *,
164 Int *,
165 Int*,
166 Int*,
167 Float*,
168 Int*,
169 Int*,
170 Double*);
171}
172void STGrid::call_ggridsd( Array<Double> &xypos,
173 Array<Complex> &spectra,
174 Int &nvispol,
175 Int &nvischan,
176 Array<Int> &flagtra,
177 Array<Int> &flagrow,
178 Array<Float> &weight,
179 Int &nrow,
180 Int &irow,
181 Array<Complex> &gdata,
182 Array<Float> &gwgt,
183 Int &nx,
184 Int &ny,
185 Int &npol,
186 Int &nchan,
187 Int &support,
188 Int &sampling,
189 Vector<Float> &convFunc,
190 Int *chanMap,
191 Int *polMap )
192{
193 // parameters for gridding
194 Int idopsf = 0 ;
195 Int len = npol*nchan ;
196 Double *sumw_p = new Double[len] ;
197 {
198 Double *work_p = sumw_p ;
199 for ( Int i = 0 ; i < len ; i++ ) {
200 *work_p = 0.0 ;
201 work_p++ ;
202 }
203 }
204
205 // prepare pointer
206 Bool deletePos, deleteData, deleteWgt, deleteFlag, deleteFlagR, deleteConv, deleteDataG, deleteWgtG ;
207 Double *xy_p = xypos.getStorage( deletePos ) ;
208 const Complex *values_p = spectra.getStorage( deleteData ) ;
209 const Int *flag_p = flagtra.getStorage( deleteFlag ) ;
210 const Int *rflag_p = flagrow.getStorage( deleteFlagR ) ;
211 const Float *wgt_p = weight.getStorage( deleteWgt ) ;
212 Complex *grid_p = gdata.getStorage( deleteDataG ) ;
213 Float *wgrid_p = gwgt.getStorage( deleteWgtG ) ;
214 Float *conv_p = convFunc.getStorage( deleteConv ) ;
215
216 // pass copy of irow to ggridsd since it will be modified in theroutine
217 Int irowCopy = irow ;
218
219 // call ggridsd
220 ggridsd( xy_p,
221 values_p,
222 &nvispol,
223 &nvischan,
224 &idopsf,
225 flag_p,
226 rflag_p,
227 wgt_p,
228 &nrow,
229 &irowCopy,
230 grid_p,
231 wgrid_p,
232 &nx,
233 &ny,
234 &npol,
235 &nchan,
236 &support,
237 &sampling,
238 conv_p,
239 chanMap,
240 polMap,
241 sumw_p ) ;
242
243 // finalization
244 xypos.putStorage( xy_p, deletePos ) ;
245 spectra.freeStorage( values_p, deleteData ) ;
246 flagtra.freeStorage( flag_p, deleteFlag ) ;
247 flagrow.freeStorage( rflag_p, deleteFlagR ) ;
248 weight.freeStorage( wgt_p, deleteWgt ) ;
249 gdata.putStorage( grid_p, deleteDataG ) ;
250 gwgt.putStorage( wgrid_p, deleteWgtG ) ;
251 convFunc.putStorage( conv_p, deleteConv ) ;
252 delete sumw_p ;
253}
254
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}
399
400void STGrid::grid()
401{
402 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
403 double t0,t1 ;
404
405 // data selection
406 t0 = mathutil::gettimeofday_sec() ;
407 selectData() ;
408 t1 = mathutil::gettimeofday_sec() ;
409 os << "selectData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
410
411 setupGrid() ;
412 setupArray() ;
413
414 if ( wtype_.compare("UNIFORM") != 0 &&
415 wtype_.compare("TINT") != 0 &&
416 wtype_.compare("TSYS") != 0 &&
417 wtype_.compare("TINTSYS") != 0 ) {
418 LogIO os( LogOrigin("STGrid", "grid", WHERE) ) ;
419 os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
420 wtype_ = "UNIFORM" ;
421 }
422
423 // grid parameter
424 os << LogIO::DEBUGGING ;
425 os << "----------" << endl ;
426 os << "Data selection summary" << endl ;
427 os << " ifno = " << ifno_ << endl ;
428 os << " pollist = " << pollist_ << endl ;
429 os << " scanlist = " << scanlist_ << endl ;
430 os << "----------" << endl ;
431 os << "Grid parameter summary" << endl ;
432 os << " (nx,ny) = (" << nx_ << "," << ny_ << ")" << endl ;
433 os << " (cellx,celly) = (" << cellx_ << "," << celly_ << ")" << endl ;
434 os << " center = " << center_ << endl ;
435 os << " weighting = " << wtype_ << endl ;
436 os << " convfunc = " << convType_ << " with support " << convSupport_ << endl ;
437 os << " doclip = " << (doclip_?"True":"False") << endl ;
438 os << "----------" << LogIO::POST ;
439 os << LogIO::NORMAL ;
440
441 if ( doclip_ )
442 gridPerRowWithClipping() ;
443 else
444 gridPerRow() ;
445}
446
447void STGrid::updateChunkShape()
448{
449 // TODO: nchunk_ must be determined from nchan_, npol_, and (nx_,ny_)
450 // by considering data size to be allocated for ggridsd input/output
451 nchunk_ = 400 ;
452 nchunk_ = min( nchunk_, nrow_ ) ;
453 vshape_ = IPosition( 1, nchunk_ ) ;
454 wshape_ = IPosition( 2, nchan_, nchunk_ ) ;
455 dshape_ = IPosition( 2, 2, nchunk_ ) ;
456}
457
458struct STGChunk {
459 Int nrow ;
460 Array<Complex> spectra;
461 Array<Int> flagtra;
462 Array<Int> rflag;
463 Array<Float> weight;
464 Array<Double> direction;
465 STGChunk(IPosition const &wshape, IPosition const &vshape,
466 IPosition const &dshape)
467 : spectra(wshape), flagtra(wshape), rflag(vshape), weight(wshape),
468 direction(dshape)
469 { }
470};
471
472struct STCommonData {
473 Int gnx;
474 Int gny;
475 Int *chanMap;
476 Vector<Float> convFunc ;
477 Array<Complex> gdataArrC;
478 Array<Float> gwgtArr;
479 STCommonData(IPosition const &gshape, Array<Float> const &data)
480 : gdataArrC(gshape, 0.0), gwgtArr(data) {}
481};
482
483struct STCommonDataWithClipping {
484 Int gnx;
485 Int gny;
486 Int *chanMap;
487 Vector<Float> convFunc ;
488 Array<Complex> gdataArrC;
489 Array<Float> gwgtArr;
490 Array<Int> npoints ;
491 Array<Complex> clipMin ;
492 Array<Float> clipWMin ;
493 Array<Float> clipCMin ;
494 Array<Complex> clipMax ;
495 Array<Float> clipWMax ;
496 Array<Float> clipCMax ;
497 STCommonDataWithClipping(IPosition const &gshape,
498 IPosition const &pshape,
499 Array<Float> const &data)
500 : gdataArrC(gshape, 0.0),
501 gwgtArr(data),
502 npoints(pshape, 0),
503 clipMin(gshape, Complex(FLT_MAX,0.0)),
504 clipWMin(gshape, 0.0),
505 clipCMin(gshape, 0.0),
506 clipMax(gshape, Complex(-FLT_MAX,0.0)),
507 clipWMax(gshape, 0.0),
508 clipCMax(gshape, 0.0)
509 {}
510};
511
512#define DO_AHEAD 3
513
514struct STContext {
515 STCommonData &common;
516 FIFO<STGChunk *, DO_AHEAD> queue;
517 STGrid *const self;
518 const Int pol;
519 STContext(STGrid *obj, STCommonData &common, Int pol)
520 : self(obj), common(common), pol(pol) {}
521};
522
523struct STContextWithClipping {
524 STCommonDataWithClipping &common;
525 FIFO<STGChunk *, DO_AHEAD> queue;
526 STGrid *const self;
527 const Int pol;
528 STContextWithClipping(STGrid *obj, STCommonDataWithClipping &common, Int pol)
529 : self(obj), common(common), pol(pol) {}
530};
531
532
533bool STGrid::produceChunk(void *ctx) throw(PCException)
534{
535 STContext &context = *(STContext *)ctx;
536 if ( context.self->nprocessed_ >= context.self->nrow_ ) {
537 return false;
538 }
539 STGChunk *chunk = new STGChunk(context.self->wshape_,
540 context.self->vshape_,
541 context.self->dshape_);
542
543 double t0 = mathutil::gettimeofday_sec() ;
544 chunk->nrow = context.self->getDataChunk(
545 context.self->wshape_, context.self->vshape_, context.self->dshape_,
546 chunk->spectra, chunk->direction,
547 chunk->flagtra, chunk->rflag, chunk->weight);
548 double t1 = mathutil::gettimeofday_sec() ;
549 context.self->eGetData_ += t1-t0 ;
550
551 context.queue.lock();
552 context.queue.put(chunk);
553 context.queue.unlock();
554 return true;
555}
556
557void STGrid::consumeChunk(void *ctx) throw(PCException)
558{
559 STContext &context = *(STContext *)ctx;
560 STGChunk *chunk = NULL;
561 try {
562 context.queue.lock();
563 chunk = context.queue.get();
564 context.queue.unlock();
565 } catch (FullException &e) {
566 context.queue.unlock();
567 // TODO: log error
568 throw PCException();
569 }
570
571 double t0, t1 ;
572 // world -> pixel
573 Array<Double> xypos( context.self->dshape_ ) ;
574 t0 = mathutil::gettimeofday_sec() ;
575 context.self->toPixel( chunk->direction, xypos ) ;
576 t1 = mathutil::gettimeofday_sec() ;
577 context.self->eToPixel_ += t1-t0 ;
578
579 // call ggridsd
580 Int nvispol = 1 ;
581 Int irow = -1 ;
582 t0 = mathutil::gettimeofday_sec() ;
583 context.self->call_ggridsd( xypos,
584 chunk->spectra,
585 nvispol,
586 context.self->nchan_,
587 chunk->flagtra,
588 chunk->rflag,
589 chunk->weight,
590 chunk->nrow,
591 irow,
592 context.common.gdataArrC,
593 context.common.gwgtArr,
594 context.common.gnx,
595 context.common.gny,
596 context.self->npol_,
597 context.self->nchan_,
598 context.self->convSupport_,
599 context.self->convSampling_,
600 context.common.convFunc,
601 context.common.chanMap,
602 (Int*)&context.pol ) ;
603 t1 = mathutil::gettimeofday_sec() ;
604 context.self->eGGridSD_ += t1-t0 ;
605
606 delete chunk;
607}
608
609void STGrid::gridPerRow()
610{
611 LogIO os( LogOrigin("STGrid", "gridPerRow", WHERE) ) ;
612 double t0, t1 ;
613
614
615 // grid data
616 // Extend grid plane with convSupport_
617 // Int gnx = nx_+convSupport_*2 ;
618 // Int gny = ny_+convSupport_*2 ;
619 Int gnx = nx_;
620 Int gny = ny_;
621
622 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
623 // 2011/12/20 TN
624 // data_ and gwgtArr share storage
625 data_.resize( gshape ) ;
626 data_ = 0.0 ;
627 STCommonData common = STCommonData(gshape, data_);
628 common.gnx = gnx ;
629 common.gny = gny ;
630
631 // parameters for gridding
632 Int *chanMap = new Int[nchan_] ;
633 for ( Int i = 0 ; i < nchan_ ; i++ ) {
634 chanMap[i] = i ;
635 }
636 common.chanMap = chanMap;
637
638 // convolution kernel
639 t0 = mathutil::gettimeofday_sec() ;
640 setConvFunc( common.convFunc ) ;
641 t1 = mathutil::gettimeofday_sec() ;
642 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
643
644 // for performance check
645 eGetData_ = 0.0 ;
646 eToPixel_ = 0.0 ;
647 eGGridSD_ = 0.0 ;
648 double eInitPol = 0.0 ;
649
650 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
651 initTable( ifile ) ;
652
653 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
654 Broker broker = Broker(produceChunk, consumeChunk);
655 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
656 t0 = mathutil::gettimeofday_sec() ;
657 initPol( ipol ) ; // set ptab_ and attach()
658 t1 = mathutil::gettimeofday_sec() ;
659 eInitPol += t1-t0 ;
660
661 STContext context(this, common, ipol);
662
663 os << "start pol " << ipol << LogIO::POST ;
664
665 nprocessed_ = 0 ;
666#if 1
667 broker.runProducerAsMasterThread(&context, DO_AHEAD);
668#else
669 for (;;) {
670 bool produced = produceChunk(&context);
671 if (! produced) {
672 break;
673 }
674 consumeChunk(&context);
675 }
676#endif
677
678 os << "end pol " << ipol << LogIO::POST ;
679
680 }
681 os << "end table " << ifile << LogIO::POST ;
682 }
683 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
684 os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
685 os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
686 os << "ggridsd: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
687 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
688 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
689
690 delete chanMap ;
691
692 // set data
693 setData( common.gdataArrC, common.gwgtArr ) ;
694
695}
696
697void STGrid::consumeChunkWithClipping(void *ctx) throw(PCException)
698{
699 STContextWithClipping &context = *(STContextWithClipping *)ctx;
700 STGChunk *chunk = NULL;
701 try {
702 context.queue.lock();
703 chunk = context.queue.get();
704 context.queue.unlock();
705 } catch (FullException &e) {
706 context.queue.unlock();
707 // TODO: log error
708 throw PCException();
709 }
710
711 double t0, t1 ;
712 // world -> pixel
713 Array<Double> xypos( context.self->dshape_ ) ;
714 t0 = mathutil::gettimeofday_sec() ;
715 context.self->toPixel( chunk->direction, xypos ) ;
716 t1 = mathutil::gettimeofday_sec() ;
717 context.self->eToPixel_ += t1-t0 ;
718
719 // call ggridsd
720 Int nvispol = 1 ;
721 Int irow = -1 ;
722 t0 = mathutil::gettimeofday_sec() ;
723 context.self->call_ggridsd2( xypos,
724 chunk->spectra,
725 nvispol,
726 context.self->nchan_,
727 chunk->flagtra,
728 chunk->rflag,
729 chunk->weight,
730 chunk->nrow,
731 irow,
732 context.common.gdataArrC,
733 context.common.gwgtArr,
734 context.common.npoints,
735 context.common.clipMin,
736 context.common.clipWMin,
737 context.common.clipCMin,
738 context.common.clipMax,
739 context.common.clipWMax,
740 context.common.clipCMax,
741 context.common.gnx,
742 context.common.gny,
743 context.self->npol_,
744 context.self->nchan_,
745 context.self->convSupport_,
746 context.self->convSampling_,
747 context.common.convFunc,
748 context.common.chanMap,
749 (Int*)&context.pol ) ;
750 t1 = mathutil::gettimeofday_sec() ;
751 context.self->eGGridSD_ += t1-t0 ;
752
753 delete chunk;
754}
755
756void STGrid::gridPerRowWithClipping()
757{
758 LogIO os( LogOrigin("STGrid", "gridPerRowWithClipping", WHERE) ) ;
759 double t0, t1 ;
760
761
762 // grid data
763 // Extend grid plane with convSupport_
764 // Int gnx = nx_+convSupport_*2 ;
765 // Int gny = ny_+convSupport_*2 ;
766 Int gnx = nx_;
767 Int gny = ny_;
768
769 IPosition gshape( 4, gnx, gny, npol_, nchan_ ) ;
770 IPosition pshape( 3, gnx, gny, npol_ ) ;
771 // 2011/12/20 TN
772 // data_ and gwgtArr share storage
773 data_.resize( gshape ) ;
774 data_ = 0.0 ;
775 STCommonDataWithClipping common = STCommonDataWithClipping( gshape,
776 pshape,
777 data_ ) ;
778 common.gnx = gnx ;
779 common.gny = gny ;
780
781 // parameters for gridding
782 Int *chanMap = new Int[nchan_] ;
783 for ( Int i = 0 ; i < nchan_ ; i++ ) {
784 chanMap[i] = i ;
785 }
786 common.chanMap = chanMap;
787
788 // convolution kernel
789 t0 = mathutil::gettimeofday_sec() ;
790 setConvFunc( common.convFunc ) ;
791 t1 = mathutil::gettimeofday_sec() ;
792 os << "setConvFunc: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
793
794 // for performance check
795 eGetData_ = 0.0 ;
796 eToPixel_ = 0.0 ;
797 eGGridSD_ = 0.0 ;
798 double eInitPol = 0.0 ;
799
800 for ( uInt ifile = 0 ; ifile < nfile_ ; ifile++ ) {
801 initTable( ifile ) ;
802
803 os << "start table " << ifile << ": " << infileList_[ifile] << LogIO::POST ;
804 Broker broker = Broker(produceChunk, consumeChunkWithClipping);
805 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
806 t0 = mathutil::gettimeofday_sec() ;
807 initPol( ipol ) ; // set ptab_ and attach()
808 t1 = mathutil::gettimeofday_sec() ;
809 eInitPol += t1-t0 ;
810
811 STContextWithClipping context(this, common, ipol);
812
813 os << "start pol " << ipol << LogIO::POST ;
814
815 nprocessed_ = 0 ;
816#if 1
817 broker.runProducerAsMasterThread(&context, DO_AHEAD);
818#else
819 for (;;) {
820 bool produced = produceChunk(&context);
821 if (! produced) {
822 break;
823 }
824 consumeChunkWithClipping(&context);
825 }
826#endif
827
828 os << "end pol " << ipol << LogIO::POST ;
829
830 }
831 os << "end table " << ifile << LogIO::POST ;
832 }
833 os << "initPol: elapsed time is " << eInitPol << " sec." << LogIO::POST ;
834 os << "getData: elapsed time is " << eGetData_-eToInt-eGetWeight << " sec." << LogIO::POST ;
835 os << "toPixel: elapsed time is " << eToPixel_ << " sec." << LogIO::POST ;
836 os << "ggridsd2: elapsed time is " << eGGridSD_ << " sec." << LogIO::POST ;
837 os << "toInt: elapsed time is " << eToInt << " sec." << LogIO::POST ;
838 os << "getWeight: elapsed time is " << eGetWeight << " sec." << LogIO::POST ;
839
840 delete chanMap ;
841
842 // clip min and max in each grid
843// os << "BEFORE CLIPPING" << LogIO::POST ;
844// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
845// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
846 t0 = mathutil::gettimeofday_sec() ;
847 clipMinMax( common.gdataArrC,
848 common.gwgtArr,
849 common.npoints,
850 common.clipMin,
851 common.clipWMin,
852 common.clipCMin,
853 common.clipMax,
854 common.clipWMax,
855 common.clipCMax ) ;
856 t1 = mathutil::gettimeofday_sec() ;
857 os << "clipMinMax: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
858// os << "AFTER CLIPPING" << LogIO::POST ;
859// os << "gdataArrC=" << common.gdataArrC << LogIO::POST ;
860// os << "gwgtArr=" << common.gwgtArr << LogIO::POST ;
861
862 // set data
863 setData( common.gdataArrC, common.gwgtArr ) ;
864
865}
866
867void STGrid::clipMinMax( Array<Complex> &grid,
868 Array<Float> &weight,
869 Array<Int> &npoints,
870 Array<Complex> &clipmin,
871 Array<Float> &clipwmin,
872 Array<Float> &clipcmin,
873 Array<Complex> &clipmax,
874 Array<Float> &clipwmax,
875 Array<Float> &clipcmax )
876{
877 //LogIO os( LogOrigin("STGrid","clipMinMax",WHERE) ) ;
878
879 // prepare pointers
880 Bool delG, delW, delNP, delCMin, delCWMin, delCCMin, delCMax, delCWMax, delCCMax ;
881 Complex *grid_p = grid.getStorage( delG ) ;
882 Float *wgt_p = weight.getStorage( delW ) ;
883 const Int *npts_p = npoints.getStorage( delNP ) ;
884 const Complex *cmin_p = clipmin.getStorage( delCMin ) ;
885 const Float *cwmin_p = clipwmin.getStorage( delCWMin ) ;
886 const Float *ccmin_p = clipcmin.getStorage( delCCMin ) ;
887 const Complex *cmax_p = clipmax.getStorage( delCMax ) ;
888 const Float *cwmax_p = clipwmax.getStorage( delCWMax ) ;
889 const Float *ccmax_p = clipcmax.getStorage( delCCMax ) ;
890
891 const IPosition &gshape = grid.shape() ;
892 long offset = gshape[0] * gshape[1] * gshape[2] ; // nx * ny * npol
893 Int nchan = gshape[3] ;
894 long origin = nchan * offset ;
895 for ( long i = 0 ; i < offset ; i++ ) {
896 if ( *npts_p > 2 ) {
897 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
898 // clip minimum and maximum
899 *grid_p -= (*cmin_p)*(*cwmin_p)*(*ccmin_p)
900 + (*cmax_p)*(*cwmax_p)*(*ccmax_p) ;
901 *wgt_p -= (*cwmin_p)*(*ccmin_p)
902 + (*cwmax_p)*(*ccmax_p) ;
903
904 grid_p += offset ;
905 wgt_p += offset ;
906 cmin_p += offset ;
907 cwmin_p += offset ;
908 ccmin_p += offset ;
909 cmax_p += offset ;
910 cwmax_p += offset ;
911 ccmax_p += offset ;
912 }
913 grid_p -= origin ;
914 wgt_p -= origin ;
915 cmin_p -= origin ;
916 cwmin_p -= origin ;
917 ccmin_p -= origin ;
918 cmax_p -= origin ;
919 cwmax_p -= origin ;
920 ccmax_p -= origin ;
921 }
922 grid_p++ ;
923 wgt_p++ ;
924 npts_p++ ;
925 cmin_p++ ;
926 cwmin_p++ ;
927 ccmin_p++ ;
928 cmax_p++ ;
929 cwmax_p++ ;
930 ccmax_p++ ;
931 }
932 grid_p -= offset ;
933 wgt_p -= offset ;
934 npts_p -= offset ;
935 cmin_p -= offset ;
936 cwmin_p -= offset ;
937 ccmin_p -= offset ;
938 cmax_p -= offset ;
939 cwmax_p -= offset ;
940 ccmax_p -= offset ;
941
942 // finalization
943 grid.putStorage( grid_p, delG ) ;
944 weight.putStorage( wgt_p, delW ) ;
945 npoints.freeStorage( npts_p, delNP ) ;
946 clipmin.freeStorage( cmin_p, delCMin ) ;
947 clipwmin.freeStorage( cwmin_p, delCWMin ) ;
948 clipcmin.freeStorage( ccmin_p, delCCMin ) ;
949 clipmax.freeStorage( cmax_p, delCMax ) ;
950 clipwmax.freeStorage( cwmax_p, delCWMax ) ;
951 clipcmax.freeStorage( ccmax_p, delCCMax ) ;
952}
953
954void STGrid::initPol( Int ipol )
955{
956 LogIO os( LogOrigin("STGrid","initPol",WHERE) ) ;
957 if ( npolOrg_ == 1 ) {
958 os << "single polarization data." << LogIO::POST ;
959 ptab_ = tab_ ;
960 }
961 else
962 ptab_ = tab_( tab_.col("POLNO") == (uInt)ipol ) ;
963
964 attach( ptab_ ) ;
965}
966
967void STGrid::initTable( uInt idx )
968{
969 tab_ = tableList_[idx] ;
970 nrow_ = rows_[idx] ;
971 updateChunkShape() ;
972}
973
974void STGrid::setData( Array<Complex> &gdata,
975 Array<Float> &gwgt )
976{
977 // 2011/12/20 TN
978 // gwgt and data_ share storage
979 LogIO os( LogOrigin("STGrid","setData",WHERE) ) ;
980 double t0, t1 ;
981 t0 = mathutil::gettimeofday_sec() ;
982 uInt len = data_.nelements() ;
983 const Complex *w1_p ;
984 Float *w2_p ;
985 Bool b1, b2 ;
986 const Complex *gdata_p = gdata.getStorage( b1 ) ;
987 Float *gwgt_p = data_.getStorage( b2 ) ;
988 w1_p = gdata_p ;
989 w2_p = gwgt_p ;
990 for ( uInt i = 0 ; i < len ; i++ ) {
991 if ( *w2_p > 0.0 ) *w2_p = (*w1_p).real() / *w2_p ;
992 w1_p++ ;
993 w2_p++ ;
994 }
995 gdata.freeStorage( gdata_p, b1 ) ;
996 data_.putStorage( gwgt_p, b2 ) ;
997 t1 = mathutil::gettimeofday_sec() ;
998 os << "setData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
999}
1000
1001void STGrid::setupGrid()
1002{
1003 Double xmin,xmax,ymin,ymax ;
1004 mapExtent( xmin, xmax, ymin, ymax ) ;
1005
1006 setupGrid( nxUI_, nyUI_, cellxUI_, cellyUI_,
1007 xmin, xmax, ymin, ymax, centerUI_ ) ;
1008}
1009
1010void STGrid::setupGrid( Int &nx,
1011 Int &ny,
1012 String &cellx,
1013 String &celly,
1014 Double &xmin,
1015 Double &xmax,
1016 Double &ymin,
1017 Double &ymax,
1018 String &center )
1019{
1020 LogIO os( LogOrigin("STGrid","setupGrid",WHERE) ) ;
1021 //cout << "nx=" << nx << ", ny=" << ny << endl ;
1022
1023 // center position
1024 if ( center.size() == 0 ) {
1025 center_(0) = 0.5 * ( xmin + xmax ) ;
1026 center_(1) = 0.5 * ( ymin + ymax ) ;
1027 }
1028 else {
1029 String::size_type pos0 = center.find( " " ) ;
1030 if ( pos0 == String::npos ) {
1031 throw AipsError( "bad string format in parameter center" ) ;
1032 }
1033 String::size_type pos1 = center.find( " ", pos0+1 ) ;
1034 String typestr, xstr, ystr ;
1035 if ( pos1 != String::npos ) {
1036 typestr = center.substr( 0, pos0 ) ;
1037 xstr = center.substr( pos0+1, pos1-pos0 ) ;
1038 ystr = center.substr( pos1+1 ) ;
1039 // todo: convert to J2000 (or direction ref for DIRECTION column)
1040 }
1041 else {
1042 typestr = "J2000" ;
1043 xstr = center.substr( 0, pos0 ) ;
1044 ystr = center.substr( pos0+1 ) ;
1045 }
1046 QuantumHolder qh ;
1047 String err ;
1048 qh.fromString( err, xstr ) ;
1049 Quantum<Double> xcen = qh.asQuantumDouble() ;
1050 qh.fromString( err, ystr ) ;
1051 Quantum<Double> ycen = qh.asQuantumDouble() ;
1052 center_(0) = xcen.getValue( "rad" ) ;
1053 center_(1) = ycen.getValue( "rad" ) ;
1054 }
1055
1056
1057 nx_ = nx ;
1058 ny_ = ny ;
1059 if ( nx < 0 && ny > 0 ) {
1060 nx_ = ny ;
1061 ny_ = ny ;
1062 }
1063 if ( ny < 0 && nx > 0 ) {
1064 nx_ = nx ;
1065 ny_ = nx ;
1066 }
1067
1068 //Double wx = xmax - xmin ;
1069 //Double wy = ymax - ymin ;
1070 Double wx = max( abs(xmax-center_(0)), abs(xmin-center_(0)) ) * 2 ;
1071 Double wy = max( abs(ymax-center_(1)), abs(ymin-center_(1)) ) * 2 ;
1072 // take 10% margin
1073 wx *= 1.10 ;
1074 wy *= 1.10 ;
1075
1076 Quantum<Double> qcellx ;
1077 Quantum<Double> qcelly ;
1078 //cout << "nx_ = " << nx_ << ", ny_ = " << ny_ << endl ;
1079 if ( cellx.size() != 0 && celly.size() != 0 ) {
1080 readQuantity( qcellx, cellx ) ;
1081 readQuantity( qcelly, celly ) ;
1082 }
1083 else if ( celly.size() != 0 ) {
1084 os << "Using celly to x-axis..." << LogIO::POST ;
1085 readQuantity( qcelly, celly ) ;
1086 qcellx = qcelly ;
1087 }
1088 else if ( cellx.size() != 0 ) {
1089 os << "Using cellx to y-axis..." << LogIO::POST ;
1090 readQuantity( qcellx, cellx ) ;
1091 qcelly = qcellx ;
1092 }
1093 else {
1094 if ( nx_ < 0 ) {
1095 os << "No user preference in grid setting. Using default..." << LogIO::POST ;
1096 readQuantity( qcellx, "1.0arcmin" ) ;
1097 qcelly = qcellx ;
1098 }
1099 else {
1100 if ( wx == 0.0 ) {
1101 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1102 wx = 0.00290888 ;
1103 }
1104 if ( wy == 0.0 ) {
1105 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1106 wy = 0.00290888 ;
1107 }
1108 qcellx = Quantum<Double>( wx/nx_, "rad" ) ;
1109 qcelly = Quantum<Double>( wy/ny_, "rad" ) ;
1110 }
1111 }
1112 cellx_ = qcellx.getValue( "rad" ) ;
1113 celly_ = qcelly.getValue( "rad" ) ;
1114 if ( nx_ < 0 ) {
1115 if ( wx == 0.0 ) {
1116 os << "Using default spatial extent (10arcmin) in x" << LogIO::POST ;
1117 wx = 0.00290888 ;
1118 }
1119 if ( wy == 0.0 ) {
1120 os << "Using default spatial extent (10arcmin) in y" << LogIO::POST ;
1121 wy = 0.00290888 ;
1122 }
1123 nx_ = Int( ceil( wx/cellx_ ) ) ;
1124 ny_ = Int( ceil( wy/celly_ ) ) ;
1125 }
1126}
1127
1128void STGrid::mapExtent( Double &xmin, Double &xmax,
1129 Double &ymin, Double &ymax )
1130{
1131 //LogIO os( LogOrigin("STGrid","mapExtent",WHERE) ) ;
1132 directionCol_.attach( tableList_[0], "DIRECTION" ) ;
1133 Matrix<Double> direction = directionCol_.getColumn() ;
1134 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1135 minMax( xmin, xmax, direction.row( 0 ) ) ;
1136 minMax( ymin, ymax, direction.row( 1 ) ) ;
1137 Double amin, amax, bmin, bmax ;
1138 for ( uInt i = 1 ; i < nfile_ ; i++ ) {
1139 directionCol_.attach( tableList_[i], "DIRECTION" ) ;
1140 direction.assign( directionCol_.getColumn() ) ;
1141 //os << "dirCol.nrow() = " << dirCol.nrow() << LogIO::POST ;
1142 minMax( amin, amax, direction.row( 0 ) ) ;
1143 minMax( bmin, bmax, direction.row( 1 ) ) ;
1144 xmin = min( xmin, amin ) ;
1145 xmax = max( xmax, amax ) ;
1146 ymin = min( ymin, bmin ) ;
1147 ymax = max( ymax, bmax ) ;
1148 }
1149 //os << "(xmin,xmax)=(" << xmin << "," << xmax << ")" << LogIO::POST ;
1150 //os << "(ymin,ymax)=(" << ymin << "," << ymax << ")" << LogIO::POST ;
1151}
1152
1153void STGrid::selectData()
1154{
1155 LogIO os( LogOrigin("STGrid","selectData",WHERE) ) ;
1156 Int ifno = ifno_ ;
1157 tableList_.resize( nfile_ ) ;
1158 if ( ifno == -1 ) {
1159 Table taborg( infileList_[0] ) ;
1160 ROScalarColumn<uInt> ifnoCol( taborg, "IFNO" ) ;
1161 ifno = ifnoCol( 0 ) ;
1162 os << LogIO::WARN
1163 << "IFNO is not given. Using default IFNO: " << ifno << LogIO::POST ;
1164 }
1165 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1166 Table taborg( infileList_[i] ) ;
1167 TableExprNode node ;
1168 if ( isMultiIF( taborg ) ) {
1169 os << "apply selection on IFNO" << LogIO::POST ;
1170 node = taborg.col("IFNO") == ifno ;
1171 }
1172 if ( scanlist_.size() > 0 ) {
1173 os << "apply selection on SCANNO" << LogIO::POST ;
1174 node = node && taborg.col("SCANNO").in( scanlist_ ) ;
1175 }
1176 if ( node.isNull() ) {
1177 tableList_[i] = taborg ;
1178 }
1179 else {
1180 tableList_[i] = taborg( node ) ;
1181 }
1182 os << "tableList_[" << i << "].nrow()=" << tableList_[i].nrow() << LogIO::POST ;
1183 if ( tableList_[i].nrow() == 0 ) {
1184 os << LogIO::SEVERE
1185 << "No corresponding rows for given selection: IFNO " << ifno ;
1186 if ( scanlist_.size() > 0 )
1187 os << " SCANNO " << scanlist_ ;
1188 os << LogIO::EXCEPTION ;
1189 }
1190 }
1191}
1192
1193Bool STGrid::isMultiIF( Table &tab )
1194{
1195 ROScalarColumn<uInt> ifnoCol( tab, "IFNO" ) ;
1196 Vector<uInt> ifnos = ifnoCol.getColumn() ;
1197 return anyNE( ifnos, ifnos[0] ) ;
1198}
1199
1200void STGrid::attach( Table &tab )
1201{
1202 // attach to table
1203 spectraCol_.attach( tab, "SPECTRA" ) ;
1204 flagtraCol_.attach( tab, "FLAGTRA" ) ;
1205 directionCol_.attach( tab, "DIRECTION" ) ;
1206 flagRowCol_.attach( tab, "FLAGROW" ) ;
1207 tsysCol_.attach( tab, "TSYS" ) ;
1208 intervalCol_.attach( tab, "INTERVAL" ) ;
1209}
1210
1211Int STGrid::getDataChunk(
1212 IPosition const &wshape,
1213 IPosition const &vshape,
1214 IPosition const &dshape,
1215 Array<Complex> &spectra,
1216 Array<Double> &direction,
1217 Array<Int> &flagtra,
1218 Array<Int> &rflag,
1219 Array<Float> &weight )
1220{
1221 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1222
1223 Array<Float> spectraF_(wshape);
1224 Array<uChar> flagtraUC_(wshape);
1225 Array<uInt> rflagUI_(vshape);
1226 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1227 if ( nrow < nchunk_ ) {
1228 spectra.resize( spectraF_.shape() ) ;
1229 flagtra.resize( flagtraUC_.shape() ) ;
1230 rflag.resize( rflagUI_.shape() ) ;
1231 }
1232 double t0, t1 ;
1233 t0 = mathutil::gettimeofday_sec() ;
1234 convertArray( spectra, spectraF_ ) ;
1235 toInt( flagtraUC_, flagtra ) ;
1236 toInt( rflagUI_, rflag ) ;
1237 t1 = mathutil::gettimeofday_sec() ;
1238 eToInt = t1 - t0 ;
1239
1240 return nrow ;
1241}
1242
1243#if 0
1244Int STGrid::getDataChunk( Array<Complex> &spectra,
1245 Array<Double> &direction,
1246 Array<Int> &flagtra,
1247 Array<Int> &rflag,
1248 Array<Float> &weight )
1249{
1250 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1251 Int nrow = getDataChunk( spectraF_, direction, flagtraUC_, rflagUI_, weight ) ;
1252 if ( nrow < nchunk_ ) {
1253 spectra.resize( spectraF_.shape() ) ;
1254 flagtra.resize( flagtraUC_.shape() ) ;
1255 rflag.resize( rflagUI_.shape() ) ;
1256 }
1257 double t0, t1 ;
1258 t0 = mathutil::gettimeofday_sec() ;
1259 convertArray( spectra, spectraF_ ) ;
1260 toInt( flagtraUC_, flagtra ) ;
1261 toInt( rflagUI_, rflag ) ;
1262 t1 = mathutil::gettimeofday_sec() ;
1263 eToInt = t1 - t0 ;
1264
1265 return nrow ;
1266}
1267#endif
1268
1269Int STGrid::getDataChunk( Array<Float> &spectra,
1270 Array<Double> &direction,
1271 Array<uChar> &flagtra,
1272 Array<uInt> &rflag,
1273 Array<Float> &weight )
1274{
1275 LogIO os( LogOrigin("STGrid","getDataChunk",WHERE) ) ;
1276 Int nrow = spectra.shape()[1] ;
1277 Int remainingRow = nrow_ - nprocessed_ ;
1278 if ( remainingRow < nrow ) {
1279 nrow = remainingRow ;
1280 IPosition mshape( 2, nchan_, nrow ) ;
1281 IPosition vshape( 1, nrow ) ;
1282 spectra.resize( mshape ) ;
1283 flagtra.resize( mshape ) ;
1284 direction.resize( IPosition(2,2,nrow) ) ;
1285 rflag.resize( vshape ) ;
1286 weight.resize( mshape ) ;
1287 }
1288 // 2011/12/22 TN
1289 // tsys shares its storage with weight
1290 Array<Float> tsys( weight ) ;
1291 Array<Double> tint( rflag.shape() ) ;
1292
1293 Vector<uInt> rflagVec( rflag ) ;
1294 Vector<Double> tintVec( tint ) ;
1295
1296 RefRows rows( nprocessed_, nprocessed_+nrow-1, 1 ) ;
1297 //os<<LogIO::DEBUGGING<<"nprocessed_="<<nprocessed_<<": rows.nrows()="<<rows.nrows()<<LogIO::POST ;
1298 spectraCol_.getColumnCells( rows, spectra ) ;
1299 flagtraCol_.getColumnCells( rows, flagtra ) ;
1300 directionCol_.getColumnCells( rows, direction ) ;
1301 flagRowCol_.getColumnCells( rows, rflagVec ) ;
1302 intervalCol_.getColumnCells( rows, tintVec ) ;
1303 Vector<Float> tsysTemp = tsysCol_( nprocessed_ ) ;
1304 if ( tsysTemp.nelements() == (uInt)nchan_ )
1305 tsysCol_.getColumnCells( rows, tsys ) ;
1306 else
1307 tsys = tsysTemp[0] ;
1308
1309 double t0,t1 ;
1310 t0 = mathutil::gettimeofday_sec() ;
1311 getWeight( weight, tsys, tint ) ;
1312 t1 = mathutil::gettimeofday_sec() ;
1313 eGetWeight += t1-t0 ;
1314
1315 nprocessed_ += nrow ;
1316
1317 return nrow ;
1318}
1319
1320void STGrid::setupArray()
1321{
1322 LogIO os( LogOrigin("STGrid","setupArray",WHERE) ) ;
1323 ROScalarColumn<uInt> polnoCol( tableList_[0], "POLNO" ) ;
1324 Vector<uInt> pols = polnoCol.getColumn() ;
1325 //os << pols << LogIO::POST ;
1326 Vector<uInt> pollistOrg ;
1327 npolOrg_ = 0 ;
1328 uInt polno ;
1329 for ( uInt i = 0 ; i < polnoCol.nrow() ; i++ ) {
1330 //polno = polnoCol( i ) ;
1331 polno = pols( i ) ;
1332 if ( allNE( pollistOrg, polno ) ) {
1333 pollistOrg.resize( npolOrg_+1, True ) ;
1334 pollistOrg[npolOrg_] = polno ;
1335 npolOrg_++ ;
1336 }
1337 }
1338 if ( pollist_.size() == 0 )
1339 pollist_ = pollistOrg ;
1340 else {
1341 Vector<uInt> newlist ;
1342 uInt newsize = 0 ;
1343 for ( uInt i = 0 ; i < pollist_.size() ; i++ ) {
1344 if ( anyEQ( pollistOrg, pollist_[i] ) ) {
1345 newlist.resize( newsize+1, True ) ;
1346 newlist[newsize] = pollist_[i] ;
1347 newsize++ ;
1348 }
1349 }
1350 pollist_.assign( newlist ) ;
1351 }
1352 npol_ = pollist_.size() ;
1353 if ( npol_ == 0 ) {
1354 os << LogIO::SEVERE << "Empty pollist" << LogIO::EXCEPTION ;
1355 }
1356 rows_.resize( nfile_ ) ;
1357 for ( uInt i = 0 ; i < nfile_ ; i++ ) {
1358 rows_[i] = tableList_[i].nrow() / npolOrg_ ;
1359 //if ( nrow_ < rows_[i] )
1360 // nrow_ = rows_[i] ;
1361 }
1362 flagtraCol_.attach( tableList_[0], "FLAGTRA" ) ;
1363 nchan_ = flagtraCol_( 0 ).nelements() ;
1364// os << "npol_ = " << npol_ << "(" << pollist_ << ")" << endl
1365// << "nchan_ = " << nchan_ << endl
1366// << "nrow_ = " << nrow_ << LogIO::POST ;
1367}
1368
1369void STGrid::getWeight( Array<Float> &w,
1370 Array<Float> &tsys,
1371 Array<Double> &tint )
1372{
1373 LogIO os( LogOrigin("STGrid","getWeight",WHERE) ) ;
1374
1375 // 2011/12/22 TN
1376 // w (weight) and tsys share storage
1377 IPosition refShape = tsys.shape() ;
1378 Int nchan = refShape[0] ;
1379 Int nrow = refShape[1] ;
1380// os << "nchan=" << nchan << ", nrow=" << nrow << LogIO::POST ;
1381// os << "w.shape()=" << w.shape() << endl
1382// << "tsys.shape()=" << tsys.shape() << endl
1383// << "tint.shape()=" << tint.shape() << LogIO::POST ;
1384
1385 // set weight
1386 if ( wtype_.compare( "UNIFORM" ) == 0 ) {
1387 w = 1.0 ;
1388 }
1389 else if ( wtype_.compare( "TINT" ) == 0 ) {
1390 Bool b0, b1 ;
1391 Float *w_p = w.getStorage( b0 ) ;
1392 Float *w0_p = w_p ;
1393 const Double *ti_p = tint.getStorage( b1 ) ;
1394 const Double *w1_p = ti_p ;
1395 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1396 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1397 *w0_p = *w1_p ;
1398 w0_p++ ;
1399 }
1400 w1_p++ ;
1401 }
1402 w.putStorage( w_p, b0 ) ;
1403 tint.freeStorage( ti_p, b1 ) ;
1404 }
1405 else if ( wtype_.compare( "TSYS" ) == 0 ) {
1406 Bool b0 ;
1407 Float *w_p = w.getStorage( b0 ) ;
1408 Float *w0_p = w_p ;
1409 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1410 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1411 Float temp = *w0_p ;
1412 *w0_p = 1.0 / ( temp * temp ) ;
1413 w0_p++ ;
1414 }
1415 }
1416 w.putStorage( w_p, b0 ) ;
1417 }
1418 else if ( wtype_.compare( "TINTSYS" ) == 0 ) {
1419 Bool b0, b1 ;
1420 Float *w_p = w.getStorage( b0 ) ;
1421 Float *w0_p = w_p ;
1422 const Double *ti_p = tint.getStorage( b1 ) ;
1423 const Double *w1_p = ti_p ;
1424 for ( Int irow = 0 ; irow < nrow ; irow++ ) {
1425 Float interval = *w1_p ;
1426 for ( Int ichan = 0 ; ichan < nchan ; ichan++ ) {
1427 Float temp = *w0_p ;
1428 *w0_p = interval / ( temp * temp ) ;
1429 w0_p++ ;
1430 }
1431 w1_p++ ;
1432 }
1433 w.putStorage( w_p, b0 ) ;
1434 tint.freeStorage( ti_p, b1 ) ;
1435 }
1436 else {
1437 //LogIO os( LogOrigin("STGrid", "getWeight", WHERE) ) ;
1438 //os << LogIO::WARN << "Unsupported weight type '" << wtype_ << "', apply UNIFORM weight" << LogIO::POST ;
1439 w = 1.0 ;
1440 }
1441}
1442
1443void STGrid::toInt( Array<uChar> &u, Array<Int> &v )
1444{
1445 uInt len = u.nelements() ;
1446 Int *int_p = new Int[len] ;
1447 Bool deleteIt ;
1448 const uChar *data_p = u.getStorage( deleteIt ) ;
1449 Int *i_p = int_p ;
1450 const uChar *u_p = data_p ;
1451 for ( uInt i = 0 ; i < len ; i++ ) {
1452 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1453 i_p++ ;
1454 u_p++ ;
1455 }
1456 u.freeStorage( data_p, deleteIt ) ;
1457 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1458}
1459
1460void STGrid::toInt( Array<uInt> &u, Array<Int> &v )
1461{
1462 uInt len = u.nelements() ;
1463 Int *int_p = new Int[len] ;
1464 Bool deleteIt ;
1465 const uInt *data_p = u.getStorage( deleteIt ) ;
1466 Int *i_p = int_p ;
1467 const uInt *u_p = data_p ;
1468 for ( uInt i = 0 ; i < len ; i++ ) {
1469 *i_p = ( *u_p == 0 ) ? 0 : 1 ;
1470 i_p++ ;
1471 u_p++ ;
1472 }
1473 u.freeStorage( data_p, deleteIt ) ;
1474 v.takeStorage( u.shape(), int_p, TAKE_OVER ) ;
1475}
1476
1477void STGrid::toPixel( Array<Double> &world, Array<Double> &pixel )
1478{
1479 // gridding will be done on (nx_+2*convSupport_) x (ny_+2*convSupport_)
1480 // grid plane to avoid unexpected behavior on grid edge
1481 Block<Double> pixc( 2 ) ;
1482 pixc[0] = Double( nx_-1 ) * 0.5 ;
1483 pixc[1] = Double( ny_-1 ) * 0.5 ;
1484// pixc[0] = Double( nx_+2*convSupport_-1 ) * 0.5 ;
1485// pixc[1] = Double( ny_+2*convSupport_-1 ) * 0.5 ;
1486 uInt nrow = world.shape()[1] ;
1487 Bool bw, bp ;
1488 const Double *w_p = world.getStorage( bw ) ;
1489 Double *p_p = pixel.getStorage( bp ) ;
1490 const Double *ww_p = w_p ;
1491 Double *wp_p = p_p ;
1492 for ( uInt i = 0 ; i < nrow ; i++ ) {
1493 *wp_p = pixc[0] + ( *ww_p - center_[0] ) / cellx_ ;
1494 wp_p++ ;
1495 ww_p++ ;
1496 *wp_p = pixc[1] + ( *ww_p - center_[1] ) / celly_ ;
1497 wp_p++ ;
1498 ww_p++ ;
1499 }
1500 world.freeStorage( w_p, bw ) ;
1501 pixel.putStorage( p_p, bp ) ;
1502}
1503
1504void STGrid::boxFunc( Vector<Float> &convFunc, Int &convSize )
1505{
1506 convFunc = 0.0 ;
1507 for ( Int i = 0 ; i < convSize/2 ; i++ )
1508 convFunc(i) = 1.0 ;
1509}
1510
1511#define NEED_UNDERSCORES
1512#if defined(NEED_UNDERSCORES)
1513#define grdsf grdsf_
1514#endif
1515extern "C" {
1516 void grdsf(Double*, Double*);
1517}
1518void STGrid::spheroidalFunc( Vector<Float> &convFunc )
1519{
1520 convFunc = 0.0 ;
1521 for ( Int i = 0 ; i < convSampling_*convSupport_ ; i++ ) {
1522 Double nu = Double(i) / Double(convSupport_*convSampling_) ;
1523 Double val ;
1524 grdsf( &nu, &val ) ;
1525 convFunc(i) = ( 1.0 - nu * nu ) * val ;
1526 }
1527}
1528
1529void STGrid::gaussFunc( Vector<Float> &convFunc )
1530{
1531 convFunc = 0.0 ;
1532 // HWHM of the Gaussian is convSupport_ / 4
1533 // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
1534 Int len = convSampling_ * convSupport_ ;
1535 Double hwhm = len * 0.25 ;
1536 for ( Int i = 0 ; i < len ; i++ ) {
1537 Double val = Double(i) / hwhm ;
1538 convFunc(i) = exp( -log(2)*val*val ) ;
1539 }
1540}
1541
1542void STGrid::pbFunc( Vector<Float> &convFunc )
1543{
1544 convFunc = 0.0 ;
1545}
1546
1547void STGrid::setConvFunc( Vector<Float> &convFunc )
1548{
1549 convSupport_ = userSupport_ ;
1550 if ( convType_ == "BOX" ) {
1551 if ( convSupport_ < 0 )
1552 convSupport_ = 0 ;
1553 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1554 convFunc.resize( convSize ) ;
1555 boxFunc( convFunc, convSize ) ;
1556 }
1557 else if ( convType_ == "SF" ) {
1558 if ( convSupport_ < 0 )
1559 convSupport_ = 3 ;
1560 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1561 convFunc.resize( convSize ) ;
1562 spheroidalFunc( convFunc ) ;
1563 }
1564 else if ( convType_ == "GAUSS" ) {
1565 // to take into account Gaussian tail
1566 if ( convSupport_ < 0 )
1567 convSupport_ = 12 ; // 3 * 4
1568 else {
1569 convSupport_ = userSupport_ * 4 ;
1570 }
1571 Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
1572 convFunc.resize( convSize ) ;
1573 gaussFunc( convFunc ) ;
1574 }
1575 else if ( convType_ == "PB" ) {
1576 if ( convSupport_ < 0 )
1577 convSupport_ = 0 ;
1578 pbFunc( convFunc ) ;
1579 }
1580 else {
1581 throw AipsError( "Unsupported convolution function" ) ;
1582 }
1583}
1584
1585string STGrid::saveData( string outfile )
1586{
1587 LogIO os( LogOrigin("STGrid", "saveData", WHERE) ) ;
1588 double t0, t1 ;
1589 t0 = mathutil::gettimeofday_sec() ;
1590
1591 //Int polno = 0 ;
1592 String outfile_ ;
1593 if ( outfile.size() == 0 ) {
1594 if ( infileList_[0].lastchar() == '/' ) {
1595 outfile_ = infileList_[0].substr( 0, infileList_[0].size()-1 ) ;
1596 }
1597 else {
1598 outfile_ = infileList_[0] ;
1599 }
1600 outfile_ += ".grid" ;
1601 }
1602 else {
1603 outfile_ = outfile ;
1604 }
1605 Table tab ;
1606 prepareTable( tab, outfile_ ) ;
1607 IPosition dshape = data_.shape() ;
1608 Int nrow = nx_ * ny_ * npol_ ;
1609 tab.rwKeywordSet().define( "nPol", npol_ ) ;
1610 tab.addRow( nrow ) ;
1611 Vector<Double> cpix( 2 ) ;
1612 cpix(0) = Double( nx_ - 1 ) * 0.5 ;
1613 cpix(1) = Double( ny_ - 1 ) * 0.5 ;
1614 Vector<Double> dir( 2 ) ;
1615 ArrayColumn<Double> directionCol( tab, "DIRECTION" ) ;
1616 ArrayColumn<Float> spectraCol( tab, "SPECTRA" ) ;
1617 ScalarColumn<uInt> polnoCol( tab, "POLNO" ) ;
1618 Int irow = 0 ;
1619 Vector<Float> sp( nchan_ ) ;
1620 Bool bsp, bdata ;
1621 const Float *data_p = data_.getStorage( bdata ) ;
1622 Float *wsp_p, *sp_p ;
1623 const Float *wdata_p = data_p ;
1624 long step = nx_ * ny_ * npol_ ;
1625 long offset ;
1626 for ( Int iy = 0 ; iy < ny_ ; iy++ ) {
1627 dir(1) = center_(1) - ( cpix(1) - (Double)iy ) * celly_ ;
1628 for ( Int ix = 0 ; ix < nx_ ; ix++ ) {
1629 dir(0) = center_(0) - ( cpix(0) - (Double)ix ) * cellx_ ;
1630 for ( Int ipol = 0 ; ipol < npol_ ; ipol++ ) {
1631 offset = ix + iy * nx_ + ipol * nx_ * ny_ ;
1632 //os << "offset = " << offset << LogIO::POST ;
1633 sp_p = sp.getStorage( bsp ) ;
1634 wsp_p = sp_p ;
1635 wdata_p = data_p + offset ;
1636 for ( Int ichan = 0 ; ichan < nchan_ ; ichan++ ) {
1637 *wsp_p = *wdata_p ;
1638 wsp_p++ ;
1639 wdata_p += step ;
1640 }
1641 sp.putStorage( sp_p, bsp ) ;
1642 spectraCol.put( irow, sp ) ;
1643 directionCol.put( irow, dir ) ;
1644 polnoCol.put( irow, pollist_[ipol] ) ;
1645 irow++ ;
1646 }
1647 }
1648 }
1649 data_.freeStorage( data_p, bdata ) ;
1650
1651 t1 = mathutil::gettimeofday_sec() ;
1652 os << "saveData: elapsed time is " << t1-t0 << " sec." << LogIO::POST ;
1653
1654 return outfile_ ;
1655}
1656
1657void STGrid::prepareTable( Table &tab, String &name )
1658{
1659 Table t( infileList_[0], Table::Old ) ;
1660 t.deepCopy( name, Table::New, False, t.endianFormat(), True ) ;
1661 tab = Table( name, Table::Update ) ;
1662}
1663
1664}
Note: See TracBrowser for help on using the repository browser.