source: trunk/src/STGrid.cpp @ 2398

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

Reorganization of procedure:

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