casa
$Rev:20696$
|
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 <synthesis/MSVis/VisBuffer.h> 00038 #include <casa/Arrays/Array.h> 00039 #include <casa/Arrays/Vector.h> 00040 #include <synthesis/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 uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(), 00053 cfMap_p(), conjCFMap_p(), 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), 00054 timer_p() 00055 {}; 00056 // VisibilityResamplerBase(const CFStore& cfs): 00057 // uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), convFuncStore_p(), inc_p(), 00058 // cfMap_p(), conjCFMap_p() 00059 // {setConvFunc(cfs);}; 00060 00061 VisibilityResamplerBase(const VisibilityResamplerBase& other): 00062 uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(), 00063 cfMap_p(), conjCFMap_p() 00064 {copy(other);} 00065 00066 virtual ~VisibilityResamplerBase() {}; 00067 00068 VisibilityResamplerBase& operator=(const VisibilityResamplerBase& other); 00069 00070 virtual VisibilityResamplerBase* clone() = 0; 00071 00072 virtual void copy(const VisibilityResamplerBase& other); 00073 virtual void setParams(const Vector<Double>& uvwScale, 00074 const Vector<Double>& offset, 00075 const Vector<Double>& dphase) = 0; 00076 00077 virtual void setMaps(const Vector<Int>& chanMap, const Vector<Int>& polMap) = 0; 00078 virtual void setCFMaps(const Vector<Int>& cfMap, const Vector<Int>& conjCFMap)=0; 00079 virtual void setFreqMaps(const Matrix<Double>& spwChanFreqs, const Matrix<Double>& spwnChanConjFreqs) = 0; 00080 00081 virtual void setConvFunc(const CFStore& cfs) = 0; 00082 // 00083 //------------------------------------------------------------------------------ 00084 // 00085 // Re-sample the griddedData on the VisBuffer (a.k.a gridding). 00086 // 00087 // In this class, these just call the private templated version. 00088 // The first variant grids onto a double precision grid while the 00089 // second one does it on a single precision grid. 00090 // 00091 virtual void DataToGrid(Array<DComplex>& griddedData, VBStore& vbs, 00092 Matrix<Double>& sumwt, const Bool& dopsf, 00093 Bool useConjFreqCF=False) = 0; 00094 00095 virtual void DataToGrid(Array<Complex>& griddedData, VBStore& vbs, 00096 Matrix<Double>& sumwt, const Bool& dopsf, 00097 Bool useConjFreqCF=False) = 0; 00098 // 00099 //------------------------------------------------------------------------------ 00100 // 00101 // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding) 00102 // 00103 virtual void GridToData(VBStore& vbs,const Array<Complex>& griddedData) = 0; 00104 // virtual void GridToData(VBStore& vbs, Array<Complex>& griddedData); 00105 00106 virtual void ComputeResiduals(VBStore& vbs) = 0; 00107 00108 // Forward looking genealogical baggage -- required for the 00109 // MultiThreadedVisibilityResampler 00110 virtual void init(const Bool& doublePrecision) = 0; 00111 virtual void GatherGrids(Array<DComplex>& griddedData, Matrix<Double>& sumwt) = 0; 00112 virtual void GatherGrids(Array<Complex>& griddedData, Matrix<Double>& sumwt) = 0; 00113 virtual void initializePutBuffers(const Array<DComplex>& griddedData, 00114 const Matrix<Double>& sumwt) = 0; 00115 virtual void initializePutBuffers(const Array<Complex>& griddedData, 00116 const Matrix<Double>& sumwt) = 0; 00117 // 00118 // Aliases for more readable code at the FTMachine layer. 00119 // 00120 inline void finalizeToSky(Array<DComplex>& griddedData, Matrix<Double>& sumwt) 00121 {GatherGrids(griddedData, sumwt);}; 00122 inline void finalizeToSky(Array<Complex>& griddedData, Matrix<Double>& sumwt) 00123 {GatherGrids(griddedData, sumwt);}; 00124 inline void initializeToSky(const Array<DComplex>& griddedData,const Matrix<Double>& sumwt) 00125 {initializePutBuffers(griddedData, sumwt);}; 00126 inline void initializeToSky(const Array<Complex>& griddedData,const Matrix<Double>& sumwt) 00127 {initializePutBuffers(griddedData, sumwt);}; 00128 const Vector<Int> getCFMap() {return cfMap_p;}; 00129 const Vector<Int> getConjCFMap() {return conjCFMap_p;}; 00130 00131 00132 virtual void releaseBuffers() = 0; 00133 VBRow2CFMapType& getVBRow2CFMap() {return vbRow2CFMap_p;}; 00134 VBRow2CFBMapType& getVBRow2CFBMap() {return vbRow2CFBMap_p;}; 00135 virtual Int makeVBRow2CFMap(CFStore2& cfs, 00136 ConvolutionFunction& cf, 00137 const VisBuffer& vb, const Quantity& dPA, 00138 const Vector<Int>& dataChan2ImChanMap, 00139 const Vector<Int>& dataPol2ImPolMap, 00140 const Vector<Double>& pointingOffset); 00141 00142 Double runTimeG_p, runTimeDG_p, runTimeG1_p, runTimeG2_p, runTimeG3_p, runTimeG4_p, runTimeG5_p, runTimeG6_p, runTimeG7_p; 00143 Timer timer_p; 00144 // 00145 //------------------------------------------------------------------------------ 00146 //----------------------------Private parts------------------------------------- 00147 //------------------------------------------------------------------------------ 00148 // 00149 // private: 00150 protected: 00151 Vector<Double> uvwScale_p, offset_p, dphase_p; 00152 Vector<Int> chanMap_p, polMap_p; 00153 Matrix<Double> spwChanFreq_p, spwChanConjFreq_p; 00154 CFStore convFuncStore_p; 00155 // Int inc0_p, inc1_p, inc2_p, inc3_p; 00156 Vector<Int> inc_p; 00157 Int* __restrict__ incPtr_p; 00158 Vector<Int> cfMap_p, conjCFMap_p; 00159 VBRow2CFMapType vbRow2CFMap_p; 00160 VBRow2CFBMapType vbRow2CFBMap_p; 00161 00162 00163 void sgrid(Int& ndim, 00164 Double* __restrict__ pos, 00165 Int* __restrict__ loc, 00166 Int* __restrict__ off, 00167 Complex& phasor, const Int& irow, 00168 const Double* __restrict__ uvw, 00169 const Double& dphase, const Double& freq, 00170 const Double* __restrict__ scale, 00171 const Double* __restrict__ offset, 00172 const Float* __restrict__ sampling); 00173 00174 inline Bool onGrid (const Int& nx, const Int& ny, 00175 const Vector<Int>& __restrict__ loc, 00176 const Vector<Int>& __restrict__ support) __restrict__ 00177 { 00178 return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) && 00179 ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny)); 00180 }; 00181 inline Bool onGrid (const Int& nx, const Int& ny, 00182 const Int& loc0, const Int& loc1, 00183 const Int& support) __restrict__ 00184 { 00185 return (((loc0-support) >= 0 ) && ((loc0+support) < nx) && 00186 ((loc1-support) >= 0 ) && ((loc1+support) < ny)); 00187 }; 00188 00189 // Array assignment operator in CASACore requires lhs.nelements() 00190 // == 0 or lhs.nelements()=rhs.nelements() 00191 // template <class T> 00192 // inline void SETVEC(Vector<T>& lhs, const Vector<T>& rhs) 00193 // {lhs.resize(rhs.shape()); lhs = rhs;}; 00194 00195 00196 //=============================================================================== 00197 // CASACORE-LEVEL MATERIAL 00198 //=============================================================================== 00199 // Internal methods to address a 4D array. These should ulimately 00200 // moved to a Array4D class in CASACore 00201 // 00202 00203 // This is called less frequently. Currently once per VisBuffer 00204 inline void cacheAxisIncrements(const Int& n0, const Int& n1, const Int& n2, const Int& n3) 00205 { 00206 // inc0_p=1, inc1_p=inc0_p*n0, inc2_p=inc1_p*n1, inc3_p=inc2_p*n2;(void)n3; 00207 inc_p.resize(4); 00208 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; 00209 Bool D; 00210 incPtr_p = inc_p.getStorage(D); 00211 } 00212 inline void cacheAxisIncrements(const Vector<Int>& n) 00213 {cacheAxisIncrements(n[0],n[1],n[2],n[3]);} 00214 00215 inline void cacheAxisIncrements(const Vector<Int>& n, Vector<Int>& inc) 00216 {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];} 00217 00218 // Version that use internally cached inc_p 00219 // template <class T> 00220 inline void addTo4DArray(DComplex* __restrict__& store, Int* __restrict__& iPos, 00221 Complex& nvalue, Double& wt) __restrict__ 00222 {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);} 00223 00224 inline void addTo4DArray(Complex* __restrict__& store, Int* __restrict__& iPos, 00225 Complex& nvalue, Double& wt) __restrict__ 00226 {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);} 00227 00228 00229 // Version where inc_p is supplied from outside 00230 inline void addTo4DArray(DComplex* __restrict__& store, Int* __restrict__& iPos, 00231 Int* __restrict__ inc, Complex& nvalue, Double& wt) __restrict__ 00232 {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*Complex(wt));} 00233 00234 inline void addTo4DArray(Complex* __restrict__& store, Int* __restrict__& iPos, 00235 Int* __restrict__ inc, Complex& nvalue, Double& wt) __restrict__ 00236 {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*Complex(wt));} 00237 00238 00239 inline Complex getFrom4DArray(const Complex* __restrict__& store, 00240 const Int* __restrict__& iPos, 00241 const Vector<Int>& inc) 00242 // __restrict__ 00243 {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];}; 00244 00245 inline Complex getFrom4DArray(const Complex* __restrict__& store, 00246 const Vector<Int> iPos, const Vector<Int>& inc) 00247 // __restrict__ 00248 {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];}; 00249 00250 inline DComplex getFrom4DArray(const DComplex* __restrict__& store, 00251 const Int* __restrict__& iPos, 00252 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 Vector<Int> iPos, const Vector<Int>& inc) 00258 // __restrict__ 00259 {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];}; 00260 00261 00262 // The following two methods are called in the innermost loop. 00263 inline Complex getFrom4DArray(const Complex* __restrict__& store, const Int* __restrict__& iPos) 00264 // __restrict__ 00265 {return getFrom4DArray(store, iPos, inc_p);} 00266 00267 inline DComplex getFrom4DArray(const DComplex* __restrict__& store, const Int* __restrict__& iPos) 00268 // __restrict__ 00269 {return getFrom4DArray(store, iPos, inc_p);} 00270 00271 }; 00272 }; //# NAMESPACE CASA - END 00273 00274 #endif //