VisibilityResamplerBase.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //# VisibilityResamplerBase.h: Definition of the VisibilityResamplerBase class
00003 //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
00004 //# Associated Universities, Inc. Washington DC, USA.
00005 //#
00006 //# This library is free software; you can redistribute it and/or modify it
00007 //# under the terms of the GNU Library General Public License as published by
00008 //# the Free Software Foundation; either version 2 of the License, or (at your
00009 //# option) any later version.
00010 //#
00011 //# This library is distributed in the hope that it will be useful, but WITHOUT
00012 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00013 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00014 //# License for more details.
00015 //#
00016 //# You should have received a copy of the GNU Library General Public License
00017 //# along with this library; if not, write to the Free Software Foundation,
00018 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00019 //#
00020 //# Correspondence concerning AIPS++ should be addressed as follows:
00021 //#        Internet email: aips2-request@nrao.edu.
00022 //#        Postal address: AIPS++ Project Office
00023 //#                        National Radio Astronomy Observatory
00024 //#                        520 Edgemont Road
00025 //#                        Charlottesville, VA 22903-2475 USA
00026 //#
00027 //# $Id$
00028 
00029 #ifndef SYNTHESIS_VISIBILITYRESAMPLERBASE_H
00030 #define SYNTHESIS_VISIBILITYRESAMPLERBASE_H
00031 
00032 #include <synthesis/TransformMachines/CFStore.h>
00033 #include <synthesis/TransformMachines/CFStore2.h>
00034 #include <synthesis/TransformMachines/ConvolutionFunction.h>
00035 #include <synthesis/TransformMachines/Utils.h>
00036 #include <synthesis/TransformMachines/VBStore.h>
00037 #include <msvis/MSVis/VisBuffer.h>
00038 #include <casa/Arrays/Array.h>
00039 #include <casa/Arrays/Vector.h>
00040 #include <msvis/MSVis/AsynchronousTools.h>
00041 
00042 #include <casa/Logging/LogIO.h>
00043 #include <casa/Logging/LogSink.h>
00044 #include <casa/Logging/LogMessage.h>
00045 #include <casa/OS/Timer.h>
00046 
00047 namespace casa { //# NAMESPACE CASA - BEGIN
00048   class VisibilityResamplerBase
00049   {
00050   public: 
00051     VisibilityResamplerBase(): 
00052       runTimeG_p(0.0), runTimeDG_p(0.0),runTimeG1_p(0.0), runTimeG2_p(0.0), runTimeG3_p(0.0), runTimeG4_p(0.0), runTimeG5_p(0.0), runTimeG6_p(0.0), runTimeG7_p(0.0),
00053       timer_p(),
00054       uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
00055       cfMap_p(), conjCFMap_p() 
00056 
00057     {};
00058     // VisibilityResamplerBase(const CFStore& cfs): 
00059     //   uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), convFuncStore_p(), inc_p(),
00060     //   cfMap_p(), conjCFMap_p()
00061     // {setConvFunc(cfs);};
00062 
00063     VisibilityResamplerBase(const VisibilityResamplerBase& other):
00064       uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
00065       cfMap_p(), conjCFMap_p()
00066     {copy(other);}
00067 
00068     virtual ~VisibilityResamplerBase() {};
00069 
00070     VisibilityResamplerBase& operator=(const VisibilityResamplerBase& other);
00071 
00072     virtual VisibilityResamplerBase* clone() = 0;
00073 
00074     virtual void copy(const VisibilityResamplerBase& other);
00075     virtual void setParams(const Vector<Double>& uvwScale, 
00076                            const Vector<Double>& offset,
00077                            const Vector<Double>& dphase) = 0;
00078 
00079     virtual void setMaps(const Vector<Int>& chanMap, const Vector<Int>& polMap) = 0;
00080     virtual void setCFMaps(const Vector<Int>& cfMap, const Vector<Int>& conjCFMap)=0;
00081     virtual void setFreqMaps(const Matrix<Double>& spwChanFreqs, const Matrix<Double>& spwnChanConjFreqs) = 0;
00082 
00083     virtual void setConvFunc(const CFStore& cfs) = 0;
00084     //
00085     //------------------------------------------------------------------------------
00086     //
00087     // Re-sample the griddedData on the VisBuffer (a.k.a gridding).
00088     //
00089     // In this class, these just call the private templated version.
00090     // The first variant grids onto a double precision grid while the
00091     // second one does it on a single precision grid.
00092     //
00093     virtual void DataToGrid(Array<DComplex>& griddedData, VBStore& vbs, 
00094                             Matrix<Double>& sumwt, const Bool& dopsf,
00095                             Bool useConjFreqCF=False) = 0;
00096 
00097     virtual void DataToGrid(Array<Complex>& griddedData, VBStore& vbs, 
00098                             Matrix<Double>& sumwt, const Bool& dopsf,
00099                             Bool useConjFreqCF=False) = 0;
00100     //
00101     //------------------------------------------------------------------------------
00102     //
00103     // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
00104     //
00105     virtual void GridToData(VBStore& vbs,const Array<Complex>& griddedData) = 0; 
00106     //    virtual void GridToData(VBStore& vbs, Array<Complex>& griddedData); 
00107 
00108     virtual void ComputeResiduals(VBStore& vbs) = 0;
00109     
00110     // Forward looking genealogical baggage -- required for the
00111     // MultiThreadedVisibilityResampler
00112     virtual void init(const Bool& doublePrecision) = 0;
00113     virtual void GatherGrids(Array<DComplex>& griddedData, Matrix<Double>& sumwt) = 0;
00114     virtual void GatherGrids(Array<Complex>& griddedData, Matrix<Double>& sumwt) = 0;
00115     virtual void initializePutBuffers(const Array<DComplex>& griddedData,
00116                                       const Matrix<Double>& sumwt) = 0;
00117     virtual void initializePutBuffers(const Array<Complex>& griddedData,
00118                                       const Matrix<Double>& sumwt) = 0;
00119     virtual void initializeDataBuffers(VBStore& vbs)=0;
00120     //
00121     // Aliases for more readable code at the FTMachine layer.
00122     //
00123     inline void finalizeToSky(Array<DComplex>& griddedData, Matrix<Double>& sumwt) 
00124     {GatherGrids(griddedData, sumwt);};
00125     inline void finalizeToSky(Array<Complex>& griddedData, Matrix<Double>& sumwt) 
00126     {GatherGrids(griddedData, sumwt);};
00127     inline void initializeToSky(const Array<DComplex>& griddedData,const Matrix<Double>& sumwt) 
00128     {initializePutBuffers(griddedData, sumwt);};
00129     inline void initializeToSky(const Array<Complex>& griddedData,const Matrix<Double>& sumwt)
00130     {initializePutBuffers(griddedData, sumwt);};
00131     const Vector<Int> getCFMap() {return cfMap_p;};
00132     const Vector<Int> getConjCFMap() {return conjCFMap_p;};
00133 
00134     
00135     virtual void releaseBuffers() = 0;
00136     VBRow2CFMapType& getVBRow2CFMap() {return vbRow2CFMap_p;};
00137     VBRow2CFBMapType& getVBRow2CFBMap() {return vbRow2CFBMap_p;};
00138     virtual Int makeVBRow2CFMap(CFStore2& cfs,
00139                                 ConvolutionFunction& cf,
00140                                 const VisBuffer& vb, const Quantity& dPA,
00141                                 const Vector<Int>& dataChan2ImChanMap,
00142                                 const Vector<Int>& dataPol2ImPolMap,
00143                                 const Vector<Double>& pointingOffset);
00144 
00145     Double runTimeG_p, runTimeDG_p, runTimeG1_p, runTimeG2_p, runTimeG3_p, runTimeG4_p, runTimeG5_p, runTimeG6_p, runTimeG7_p;
00146     Timer timer_p;
00147     //
00148     //------------------------------------------------------------------------------
00149     //----------------------------Private parts-------------------------------------
00150     //------------------------------------------------------------------------------
00151     //
00152     //  private:
00153   protected:
00154     Vector<Double> uvwScale_p, offset_p, dphase_p;
00155     Vector<Int> chanMap_p, polMap_p;
00156     Matrix<Double> spwChanFreq_p, spwChanConjFreq_p;
00157     CFStore convFuncStore_p;
00158     //    Int inc0_p, inc1_p, inc2_p, inc3_p;
00159     Vector<Int> inc_p;
00160     Int* __restrict__ incPtr_p;
00161     Vector<Int> cfMap_p, conjCFMap_p;
00162     VBRow2CFMapType vbRow2CFMap_p;
00163     VBRow2CFBMapType vbRow2CFBMap_p;
00164     
00165 
00166     void sgrid(Int& ndim, 
00167                Double* __restrict__  pos, 
00168                Int* __restrict__  loc, 
00169                Int* __restrict__  off, 
00170                Complex& phasor, const Int& irow, 
00171                const Double* __restrict__  uvw, 
00172                const Double& dphase, const Double& freq, 
00173                const Double* __restrict__  scale, 
00174                const Double* __restrict__  offset,
00175                const Float* __restrict__  sampling);
00176 
00177     inline Bool onGrid (const Int& nx, const Int& ny, 
00178                         const Vector<Int>& __restrict__ loc, 
00179                         const Vector<Int>& __restrict__ support) __restrict__ 
00180     {
00181       return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) &&
00182               ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny));
00183     };
00184     inline Bool onGrid (const Int& nx, const Int& ny, 
00185                         const Int& loc0, const Int& loc1, 
00186                         const Int& support) __restrict__ 
00187     {
00188       return (((loc0-support) >= 0 ) && ((loc0+support) < nx) &&
00189               ((loc1-support) >= 0 ) && ((loc1+support) < ny));
00190     };
00191 
00192     // Array assignment operator in CASACore requires lhs.nelements()
00193     // == 0 or lhs.nelements()=rhs.nelements()
00194     // template <class T>
00195     // inline void SETVEC(Vector<T>& lhs, const Vector<T>& rhs)
00196     // {lhs.resize(rhs.shape()); lhs = rhs;};
00197 
00198 
00199     //===============================================================================
00200     // CASACORE-LEVEL MATERIAL
00201     //===============================================================================
00202     // Internal methods to address a 4D array.  These should ulimately
00203     // moved to a Array4D class in CASACore
00204     //
00205 
00206     // This is called less frequently.  Currently once per VisBuffer
00207     inline void cacheAxisIncrements(const Int& n0, const Int& n1, const Int& n2, const Int& n3)
00208     {
00209       //      inc0_p=1, inc1_p=inc0_p*n0, inc2_p=inc1_p*n1, inc3_p=inc2_p*n2;(void)n3;
00210       inc_p.resize(4);
00211       inc_p[0]=1; inc_p[1]=inc_p[0]*n0; inc_p[2]=inc_p[1]*n1; inc_p[3]=inc_p[2]*n2;(void)n3;
00212       Bool D;
00213       incPtr_p = inc_p.getStorage(D);
00214     }
00215     inline void cacheAxisIncrements(const Vector<Int>& n)
00216     {cacheAxisIncrements(n[0],n[1],n[2],n[3]);}
00217 
00218     inline void cacheAxisIncrements(const Vector<Int>& n, Vector<Int>& inc)
00219     {inc.resize(4);inc[0]=1; inc[1]=inc[0]*n[0]; inc[2]=inc[1]*n[1]; inc[3]=inc[2]*n[2];(void)n[3];}
00220 
00221     inline void cacheAxisIncrements(const Int n[4], Int inc[4])
00222     {inc[0]=1; inc[1]=inc[0]*n[0]; inc[2]=inc[1]*n[1]; inc[3]=inc[2]*n[2];(void)n[3];}
00223 
00224     // Version that use internally cached inc_p
00225     //    template <class T>
00226     inline void addTo4DArray(DComplex* __restrict__& store, Int* __restrict__& iPos, 
00227                              Complex& nvalue, Double& wt) __restrict__ 
00228     {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
00229 
00230     inline void addTo4DArray(Complex* __restrict__& store, Int* __restrict__& iPos, 
00231                              Complex& nvalue, Double& wt) __restrict__ 
00232     {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
00233 
00234     
00235     // Version where inc_p is supplied from outside
00236     inline void addTo4DArray(DComplex* __restrict__& store, Int* __restrict__& iPos, 
00237                              Int* __restrict__ inc, Complex& nvalue, Double& wt) __restrict__ 
00238     {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*Complex(wt));}
00239 
00240     inline void addTo4DArray(Complex* __restrict__& store, Int* __restrict__& iPos, 
00241                              Int* __restrict__ inc, Complex& nvalue, Double& wt) __restrict__ 
00242     {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*Complex(wt));}
00243 
00244 
00245     inline Complex getFrom4DArray(const Complex* __restrict__& store, 
00246                                   const Int* __restrict__& iPos, 
00247                                   const Vector<Int>& inc) 
00248     //  __restrict__ 
00249     {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
00250 
00251     inline Complex getFrom4DArray(const Complex* __restrict__& store, 
00252                                   const Vector<Int> iPos, const Vector<Int>& inc) 
00253     //  __restrict__ 
00254     {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
00255 
00256     inline DComplex getFrom4DArray(const DComplex* __restrict__& store, 
00257                                    const Int* __restrict__& iPos, 
00258                                    const Vector<Int>& inc) 
00259     //  __restrict__ 
00260     {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
00261 
00262     inline DComplex getFrom4DArray(const DComplex* __restrict__& store, 
00263                                   const Vector<Int> iPos, const Vector<Int>& inc) 
00264     //  __restrict__ 
00265     {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
00266 
00267 
00268     // The following two methods are called in the innermost loop.
00269     inline Complex getFrom4DArray(const Complex* __restrict__& store, const Int* __restrict__& iPos) 
00270     //  __restrict__ 
00271     {return getFrom4DArray(store, iPos, inc_p);}
00272 
00273     inline DComplex getFrom4DArray(const DComplex* __restrict__& store, const Int* __restrict__& iPos) 
00274     //  __restrict__ 
00275     {return getFrom4DArray(store, iPos, inc_p);}
00276 
00277   };
00278 }; //# NAMESPACE CASA - END
00279 
00280 #endif //