source: trunk/src/STGrid.cpp@ 2407

Last change on this file since 2407 was 2405, 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...

  • Added copyright information.
  • Got rid of unnecessary header files.
  • copy subtable rows explicitly using TableCopy::copySubTables since Table::deepCopy with noRows=True doesn't copy all rows including subtables.


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