source: trunk/src/STGrid.cpp @ 2422

Last change on this file since 2422 was 2422, checked in by Takeshi Nakazato, 12 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.