source: trunk/src/STGrid.cpp @ 2405

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

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


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