Line data Source code
1 : // -*- C++ -*-
2 : //# VisibilityResamplerBase.h: Definition of the VisibilityResamplerBase class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_TRANSFORM2_VISIBILITYRESAMPLERBASE_H
30 : #define SYNTHESIS_TRANSFORM2_VISIBILITYRESAMPLERBASE_H
31 :
32 : #include <synthesis/TransformMachines2/CFStore.h>
33 : #include <synthesis/TransformMachines2/CFStore2.h>
34 : #include <synthesis/TransformMachines2/ConvolutionFunction.h>
35 : #include <synthesis/TransformMachines2/Utils.h>
36 : #include <synthesis/TransformMachines2/VBStore.h>
37 : #include <msvis/MSVis/VisBuffer2.h>
38 :
39 : //#include <hpg.hpp>
40 : #include <casacore/casa/Arrays/Array.h>
41 : #include <casacore/casa/Arrays/Vector.h>
42 : #include <stdcasa/thread/AsynchronousTools.h>
43 :
44 : #include <casacore/casa/Logging/LogIO.h>
45 : #include <casacore/casa/Logging/LogSink.h>
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/OS/Timer.h>
48 : using sumofweight_fp = std::vector<std::vector<double>>;
49 :
50 : namespace casa { //# NAMESPACE CASA - BEGIN
51 : using namespace vi;
52 : namespace refim{
53 : //forward decl
54 : class AWConvFuncHolder;
55 :
56 : class VisibilityResamplerBase
57 : {
58 : public:
59 :
60 :
61 0 : VisibilityResamplerBase():
62 0 : 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),
63 0 : griddingTime(0.0),
64 0 : timer_p(),
65 0 : uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
66 0 : cfMap_p(), conjCFMap_p(), paTolerance_p(360.0), cached_phaseGrad_p(),vb2CFBMap_p()
67 :
68 0 : {};
69 : // VisibilityResamplerBase(const CFStore& cfs):
70 : // uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), convFuncStore_p(), inc_p(),
71 : // cfMap_p(), conjCFMap_p()
72 : // {setConvFunc(cfs);};
73 :
74 : VisibilityResamplerBase(const VisibilityResamplerBase& other):
75 : uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
76 : cfMap_p(), conjCFMap_p(), paTolerance_p(360.0),vb2CFBMap_p()
77 : {copy(other);}
78 :
79 0 : virtual String name() {return String("CPUResampler");}
80 : virtual Bool needCFPhaseScreen()=0;
81 0 : virtual ~VisibilityResamplerBase() {};
82 :
83 : VisibilityResamplerBase& operator=(const VisibilityResamplerBase& other);
84 :
85 : virtual VisibilityResamplerBase* clone() = 0;
86 :
87 : virtual void copy(const VisibilityResamplerBase& other);
88 : virtual void setParams(const casacore::Vector<casacore::Double>& uvwScale,
89 : const casacore::Vector<casacore::Double>& offset,
90 : const casacore::Vector<casacore::Double>& dphase) = 0;
91 :
92 : virtual void setMaps(const casacore::Vector<casacore::Int>& chanMap, const casacore::Vector<casacore::Int>& polMap) = 0;
93 : virtual void setCFMaps(const casacore::Vector<casacore::Int>& cfMap, const casacore::Vector<casacore::Int>& conjCFMap)=0;
94 : virtual void setFreqMaps(const casacore::Matrix<casacore::Double>& spwChanFreqs, const casacore::Matrix<casacore::Double>& spwnChanConjFreqs) = 0;
95 :
96 : virtual void setConvFunc(const CFStore& cfs) = 0;
97 : virtual void setPATolerance(const double& dPA) = 0;
98 :
99 : virtual void setConvFunc(std::shared_ptr<AWConvFuncHolder>& awh) = 0;
100 : //
101 : //------------------------------------------------------------------------------
102 : //
103 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding).
104 : //
105 : // In this class, these just call the private templated version.
106 : // The first variant grids onto a double precision grid while the
107 : // second one does it on a single precision grid.
108 : //
109 : virtual void DataToGrid(casacore::Array<casacore::DComplex>& griddedData, VBStore& vbs,
110 : casacore::Matrix<casacore::Double>& sumwt, const casacore::Bool& dopsf,
111 : casacore::Bool useConjFreqCF=false) = 0;
112 :
113 : virtual void DataToGrid(casacore::Array<casacore::Complex>& griddedData, VBStore& vbs,
114 : casacore::Matrix<casacore::Double>& sumwt, const casacore::Bool& dopsf,
115 : casacore::Bool useConjFreqCF=false) = 0;
116 : //
117 : //------------------------------------------------------------------------------
118 : //
119 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
120 : //
121 : virtual void GridToData(VBStore& vbs,const casacore::Array<casacore::Complex>& griddedData) = 0;
122 : // virtual void GridToData(VBStore& vbs, casacore::Array<casacore::Complex>& griddedData);
123 :
124 : virtual void ComputeResiduals(VBStore& vbs) = 0;
125 :
126 : // Forward looking genealogical baggage -- required for the
127 : // MultiThreadedVisibilityResampler
128 : virtual void init(const casacore::Bool& doublePrecision) = 0;
129 : virtual void GatherGrids(casacore::Array<casacore::DComplex>& griddedData, casacore::Matrix<casacore::Double>& sumwt) = 0;
130 : virtual void GatherGrids(casacore::Array<casacore::Complex>& griddedData, casacore::Matrix<casacore::Double>& sumwt) = 0;
131 : virtual void initializePutBuffers(const casacore::Array<casacore::DComplex>& griddedData,
132 : const casacore::Matrix<casacore::Double>& sumwt) = 0;
133 : virtual void initializePutBuffers(const casacore::Array<casacore::Complex>& griddedData,
134 : const casacore::Matrix<casacore::Double>& sumwt) = 0;
135 : virtual void initializeDataBuffers(VBStore& vbs)=0;
136 0 : virtual void initGridder(const uInt& /*NProcs*/,
137 : const std::array<unsigned, 4>& /*grid_size*/,
138 0 : const std::array<float, 2>& /*grid_scale*/) {};
139 0 : virtual sumofweight_fp resetGridder(casacore::Array<casacore::DComplex>& /*griddedData2*/) {sumofweight_fp tt; return tt;};
140 : //
141 : // Aliases for more readable code at the FTMachine layer.
142 : //
143 0 : virtual inline void finalizeToSky(casacore::Array<casacore::DComplex>& griddedData, casacore::Matrix<casacore::Double>& sumwt)
144 0 : {GatherGrids(griddedData, sumwt);};
145 0 : virtual inline void finalizeToSky(casacore::Array<casacore::Complex>& griddedData, casacore::Matrix<casacore::Double>& sumwt)
146 0 : {GatherGrids(griddedData, sumwt);};
147 0 : inline void initializeToSky(const casacore::Array<casacore::DComplex>& griddedData,const casacore::Matrix<casacore::Double>& sumwt)
148 0 : {initializePutBuffers(griddedData, sumwt);};
149 0 : inline void initializeToSky(const casacore::Array<casacore::Complex>& griddedData,const casacore::Matrix<casacore::Double>& sumwt)
150 0 : {initializePutBuffers(griddedData, sumwt);};
151 : const casacore::Vector<casacore::Int> getCFMap() {return cfMap_p;};
152 : const casacore::Vector<casacore::Int> getConjCFMap() {return conjCFMap_p;};
153 :
154 : // get pointers to storage for gridded visibilities or weight;
155 : // size to number of elements in storage; may not be implemented
156 : // by all sub-classes, if not, will return empty shared_ptr and
157 : // size equal to 0
158 : virtual std::shared_ptr<std::complex<double>> getGridPtr(size_t& size) const;
159 : virtual std::shared_ptr<double> getSumWeightsPtr(size_t& size) const;
160 :
161 : virtual void releaseBuffers() = 0;
162 : // VBRow2CFMapType& getVBRow2CFMap() {return vbRow2CFMap_p;};
163 : VB2CFBMap& getVBRow2CFBMap() {return *vb2CFBMap_p;};
164 : const VB2CFBMap& getVBRow2CFBMap()const {return *vb2CFBMap_p;};
165 : // virtual casacore::Int makeVBRow2CFBMap(CFStore2& cfs,
166 : // ConvolutionFunction& cf,
167 : // const VisBuffer2& vb, const casacore::Quantity& dPA,
168 : // const casacore::Vector<casacore::Int>& dataChan2ImChanMap,
169 : // const casacore::Vector<casacore::Int>& dataPol2ImPolMap,
170 : // const casacore::Vector<casacore::Double>& pointingOffset);
171 :
172 0 : void setFieldPhaseGrad(const casacore::Matrix<casacore::Complex>& phaseGrad) {cached_phaseGrad_p.reference(phaseGrad);};
173 : void setVB2CFMap(const casacore::CountedPtr<refim::VB2CFBMap>& thisMap);
174 : // virtual void setModelImage(const std::unique_ptr<hpg::GridValueArray> /*HPGModelImage*/) {};
175 0 : virtual void setModelImage(const std::string& /*HPGModelImage*/) {};
176 0 : virtual void setModelImage(std::shared_ptr<casacore::ImageInterface<casacore::Complex> > /*HPGModelImage*/) {};
177 : casacore::Double runTimeG_p, runTimeDG_p, runTimeG1_p, runTimeG2_p, runTimeG3_p, runTimeG4_p, runTimeG5_p, runTimeG6_p, runTimeG7_p,griddingTime;
178 : casacore::Timer timer_p;
179 : //
180 : //------------------------------------------------------------------------------
181 : //----------------------------Private parts-------------------------------------
182 : //------------------------------------------------------------------------------
183 : //
184 : // private:
185 : protected:
186 : casacore::Vector<casacore::Double> uvwScale_p, offset_p, dphase_p;
187 : casacore::Vector<casacore::Int> chanMap_p, polMap_p;
188 : casacore::Matrix<casacore::Double> spwChanFreq_p, spwChanConjFreq_p;
189 : CFStore convFuncStore_p;
190 : // casacore::Int inc0_p, inc1_p, inc2_p, inc3_p;
191 : casacore::Vector<casacore::Int> inc_p;
192 : casacore::Int* __restrict__ incPtr_p;
193 : casacore::Vector<casacore::Int> cfMap_p, conjCFMap_p;
194 : // VBRow2CFMapType vbRow2CFMap_p;
195 : double paTolerance_p;
196 : casacore::Matrix<casacore::Complex> cached_phaseGrad_p;
197 : casacore::CountedPtr<refim::VB2CFBMap> vb2CFBMap_p;
198 :
199 : void sgrid(casacore::Int& ndim,
200 : casacore::Double* __restrict__ pos,
201 : casacore::Int* __restrict__ loc,
202 : casacore::Int* __restrict__ off,
203 : casacore::Complex& phasor, const casacore::Int& irow,
204 : const casacore::Double* __restrict__ uvw,
205 : const casacore::Double& dphase, const casacore::Double& freq,
206 : const casacore::Double* __restrict__ scale,
207 : const casacore::Double* __restrict__ offset,
208 : const casacore::Float* __restrict__ sampling);
209 :
210 0 : inline casacore::Bool onGrid (const casacore::Int& nx, const casacore::Int& ny,
211 : const casacore::Vector<casacore::Int>& __restrict__ loc,
212 : const casacore::Vector<casacore::Int>& __restrict__ support) __restrict__
213 : {
214 0 : return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) &&
215 0 : ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny));
216 : };
217 : inline casacore::Bool onGrid (const casacore::Int& nx, const casacore::Int& ny,
218 : const casacore::Int& loc0, const casacore::Int& loc1,
219 : const casacore::Int& support) __restrict__
220 : {
221 : return (((loc0-support) >= 0 ) && ((loc0+support) < nx) &&
222 : ((loc1-support) >= 0 ) && ((loc1+support) < ny));
223 : };
224 :
225 : // casacore::Array assignment operator in CASACore requires lhs.nelements()
226 : // == 0 or lhs.nelements()=rhs.nelements()
227 : // template <class T>
228 : // inline void SETVEC(casacore::Vector<T>& lhs, const casacore::Vector<T>& rhs)
229 : // {lhs.resize(rhs.shape()); lhs = rhs;};
230 :
231 :
232 : //===============================================================================
233 : // CASACORE-LEVEL MATERIAL
234 : //===============================================================================
235 : // Internal methods to address a 4D array. These should ulimately
236 : // moved to a Array4D class in CASACore
237 : //
238 :
239 : // This is called less frequently. Currently once per VisBuffer
240 0 : inline void cacheAxisIncrements(const casacore::Int& n0, const casacore::Int& n1, const casacore::Int& n2, const casacore::Int& n3)
241 : {
242 : // inc0_p=1, inc1_p=inc0_p*n0, inc2_p=inc1_p*n1, inc3_p=inc2_p*n2;(void)n3;
243 0 : inc_p.resize(4);
244 0 : 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;
245 : casacore::Bool D;
246 0 : incPtr_p = inc_p.getStorage(D);
247 0 : }
248 0 : inline void cacheAxisIncrements(const casacore::Vector<casacore::Int>& n)
249 0 : {cacheAxisIncrements(n[0],n[1],n[2],n[3]);}
250 :
251 0 : inline void cacheAxisIncrements(const casacore::Vector<casacore::Int>& n, casacore::Vector<casacore::Int>& inc)
252 0 : {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];}
253 :
254 : inline void cacheAxisIncrements(const casacore::Int n[4], casacore::Int inc[4])
255 : {inc[0]=1; inc[1]=inc[0]*n[0]; inc[2]=inc[1]*n[1]; inc[3]=inc[2]*n[2];(void)n[3];}
256 :
257 : // Version that use internally cached inc_p
258 : // template <class T>
259 0 : inline void addTo4DArray(casacore::DComplex* __restrict__& store, casacore::Int* __restrict__& iPos,
260 : casacore::Complex& nvalue, casacore::Double& wt) __restrict__
261 0 : {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
262 :
263 0 : inline void addTo4DArray(casacore::Complex* __restrict__& store, casacore::Int* __restrict__& iPos,
264 : casacore::Complex& nvalue, casacore::Double& wt) __restrict__
265 0 : {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
266 :
267 :
268 : // Version where inc_p is supplied from outside
269 0 : inline void addTo4DArray(casacore::DComplex* __restrict__& store, casacore::Int* __restrict__& iPos,
270 : casacore::Int* __restrict__ inc, casacore::Complex& nvalue, casacore::Double& wt) __restrict__
271 0 : {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*casacore::Complex(wt));}
272 :
273 0 : inline void addTo4DArray(casacore::Complex* __restrict__& store, casacore::Int* __restrict__& iPos,
274 : casacore::Int* __restrict__ inc, casacore::Complex& nvalue, casacore::Double& wt) __restrict__
275 0 : {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*casacore::Complex(wt));}
276 :
277 :
278 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store,
279 : const casacore::Int* __restrict__& iPos,
280 : const casacore::Vector<casacore::Int>& inc)
281 : // __restrict__
282 0 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
283 :
284 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store,
285 : const casacore::Vector<casacore::Int> iPos, const casacore::Vector<casacore::Int>& inc)
286 : // __restrict__
287 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
288 :
289 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex* __restrict__& store,
290 : const casacore::Int* __restrict__& iPos,
291 : const casacore::Vector<casacore::Int>& inc)
292 : // __restrict__
293 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
294 :
295 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex* __restrict__& store,
296 : const casacore::Vector<casacore::Int> iPos, const casacore::Vector<casacore::Int>& inc)
297 : // __restrict__
298 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
299 :
300 :
301 : // The following two methods are called in the innermost loop.
302 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store, const casacore::Int* __restrict__& iPos)
303 : // __restrict__
304 0 : {return getFrom4DArray(store, iPos, inc_p);}
305 :
306 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex* __restrict__& store, const casacore::Int* __restrict__& iPos)
307 : // __restrict__
308 : {return getFrom4DArray(store, iPos, inc_p);}
309 :
310 : };
311 : }; //# NAMESPACE CASA - END
312 : };
313 : #endif //
|
