source: trunk/src/STGrid.cpp@ 2423

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

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