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: aips2-request@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 : class VisibilityResamplerBase
54 : {
55 : public:
56 :
57 :
58 143 : VisibilityResamplerBase():
59 : 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),
60 : griddingTime(0.0),
61 : timer_p(),
62 : uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
63 143 : cfMap_p(), conjCFMap_p(), paTolerance_p(360.0), cached_phaseGrad_p(),vb2CFBMap_p()
64 :
65 143 : {};
66 : // VisibilityResamplerBase(const CFStore& cfs):
67 : // uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), convFuncStore_p(), inc_p(),
68 : // cfMap_p(), conjCFMap_p()
69 : // {setConvFunc(cfs);};
70 :
71 : VisibilityResamplerBase(const VisibilityResamplerBase& other):
72 : uvwScale_p(), offset_p(), chanMap_p(), polMap_p(), spwChanFreq_p(), spwChanConjFreq_p (), convFuncStore_p(), inc_p(),
73 : cfMap_p(), conjCFMap_p(), paTolerance_p(360.0),vb2CFBMap_p()
74 : {copy(other);}
75 :
76 200 : virtual String name() {return String("CPUResampler");}
77 : virtual Bool needCFPhaseScreen()=0;
78 143 : virtual ~VisibilityResamplerBase() {};
79 :
80 : VisibilityResamplerBase& operator=(const VisibilityResamplerBase& other);
81 :
82 : virtual VisibilityResamplerBase* clone() = 0;
83 :
84 : virtual void copy(const VisibilityResamplerBase& other);
85 : virtual void setParams(const casacore::Vector<casacore::Double>& uvwScale,
86 : const casacore::Vector<casacore::Double>& offset,
87 : const casacore::Vector<casacore::Double>& dphase) = 0;
88 :
89 : virtual void setMaps(const casacore::Vector<casacore::Int>& chanMap, const casacore::Vector<casacore::Int>& polMap) = 0;
90 : virtual void setCFMaps(const casacore::Vector<casacore::Int>& cfMap, const casacore::Vector<casacore::Int>& conjCFMap)=0;
91 : virtual void setFreqMaps(const casacore::Matrix<casacore::Double>& spwChanFreqs, const casacore::Matrix<casacore::Double>& spwnChanConjFreqs) = 0;
92 :
93 : virtual void setConvFunc(const CFStore& cfs) = 0;
94 : virtual void setPATolerance(const double& dPA) = 0;
95 : //
96 : //------------------------------------------------------------------------------
97 : //
98 : // Re-sample the griddedData on the VisBuffer (a.k.a gridding).
99 : //
100 : // In this class, these just call the private templated version.
101 : // The first variant grids onto a double precision grid while the
102 : // second one does it on a single precision grid.
103 : //
104 : virtual void DataToGrid(casacore::Array<casacore::DComplex>& griddedData, VBStore& vbs,
105 : casacore::Matrix<casacore::Double>& sumwt, const casacore::Bool& dopsf,
106 : casacore::Bool useConjFreqCF=false) = 0;
107 :
108 : virtual void DataToGrid(casacore::Array<casacore::Complex>& griddedData, VBStore& vbs,
109 : casacore::Matrix<casacore::Double>& sumwt, const casacore::Bool& dopsf,
110 : casacore::Bool useConjFreqCF=false) = 0;
111 : //
112 : //------------------------------------------------------------------------------
113 : //
114 : // Re-sample VisBuffer to a regular grid (griddedData) (a.k.a. de-gridding)
115 : //
116 : virtual void GridToData(VBStore& vbs,const casacore::Array<casacore::Complex>& griddedData) = 0;
117 : // virtual void GridToData(VBStore& vbs, casacore::Array<casacore::Complex>& griddedData);
118 :
119 : virtual void ComputeResiduals(VBStore& vbs) = 0;
120 :
121 : // Forward looking genealogical baggage -- required for the
122 : // MultiThreadedVisibilityResampler
123 : virtual void init(const casacore::Bool& doublePrecision) = 0;
124 : virtual void GatherGrids(casacore::Array<casacore::DComplex>& griddedData, casacore::Matrix<casacore::Double>& sumwt) = 0;
125 : virtual void GatherGrids(casacore::Array<casacore::Complex>& griddedData, casacore::Matrix<casacore::Double>& sumwt) = 0;
126 : virtual void initializePutBuffers(const casacore::Array<casacore::DComplex>& griddedData,
127 : const casacore::Matrix<casacore::Double>& sumwt) = 0;
128 : virtual void initializePutBuffers(const casacore::Array<casacore::Complex>& griddedData,
129 : const casacore::Matrix<casacore::Double>& sumwt) = 0;
130 : virtual void initializeDataBuffers(VBStore& vbs)=0;
131 0 : virtual void initGridder(const uInt& /*NProcs*/,
132 : const std::array<unsigned, 4>& /*grid_size*/,
133 0 : const std::array<float, 2>& /*grid_scale*/) {};
134 0 : virtual sumofweight_fp resetGridder(casacore::Array<casacore::DComplex>& /*griddedData2*/) {sumofweight_fp tt; return tt;};
135 : //
136 : // Aliases for more readable code at the FTMachine layer.
137 : //
138 214 : virtual inline void finalizeToSky(casacore::Array<casacore::DComplex>& griddedData, casacore::Matrix<casacore::Double>& sumwt)
139 214 : {GatherGrids(griddedData, sumwt);};
140 0 : virtual inline void finalizeToSky(casacore::Array<casacore::Complex>& griddedData, casacore::Matrix<casacore::Double>& sumwt)
141 0 : {GatherGrids(griddedData, sumwt);};
142 226 : inline void initializeToSky(const casacore::Array<casacore::DComplex>& griddedData,const casacore::Matrix<casacore::Double>& sumwt)
143 226 : {initializePutBuffers(griddedData, sumwt);};
144 0 : inline void initializeToSky(const casacore::Array<casacore::Complex>& griddedData,const casacore::Matrix<casacore::Double>& sumwt)
145 0 : {initializePutBuffers(griddedData, sumwt);};
146 : const casacore::Vector<casacore::Int> getCFMap() {return cfMap_p;};
147 : const casacore::Vector<casacore::Int> getConjCFMap() {return conjCFMap_p;};
148 :
149 : // get pointers to storage for gridded visibilities or weight;
150 : // size to number of elements in storage; may not be implemented
151 : // by all sub-classes, if not, will return empty shared_ptr and
152 : // size equal to 0
153 : virtual std::shared_ptr<std::complex<double>> getGridPtr(size_t& size) const;
154 : virtual std::shared_ptr<double> getSumWeightsPtr(size_t& size) const;
155 :
156 : virtual void releaseBuffers() = 0;
157 : // VBRow2CFMapType& getVBRow2CFMap() {return vbRow2CFMap_p;};
158 : VB2CFBMap& getVBRow2CFBMap() {return *vb2CFBMap_p;};
159 : // virtual casacore::Int makeVBRow2CFBMap(CFStore2& cfs,
160 : // ConvolutionFunction& cf,
161 : // const VisBuffer2& vb, const casacore::Quantity& dPA,
162 : // const casacore::Vector<casacore::Int>& dataChan2ImChanMap,
163 : // const casacore::Vector<casacore::Int>& dataPol2ImPolMap,
164 : // const casacore::Vector<casacore::Double>& pointingOffset);
165 :
166 3628308 : void setFieldPhaseGrad(const casacore::Matrix<casacore::Complex>& phaseGrad) {cached_phaseGrad_p.reference(phaseGrad);};
167 : void setVB2CFMap(const casacore::CountedPtr<refim::VB2CFBMap>& thisMap);
168 : // virtual void setModelImage(const std::unique_ptr<hpg::GridValueArray> /*HPGModelImage*/) {};
169 0 : virtual void setModelImage(const std::string& /*HPGModelImage*/) {};
170 0 : virtual void setModelImage(std::shared_ptr<casacore::ImageInterface<casacore::Complex> > /*HPGModelImage*/) {};
171 : casacore::Double runTimeG_p, runTimeDG_p, runTimeG1_p, runTimeG2_p, runTimeG3_p, runTimeG4_p, runTimeG5_p, runTimeG6_p, runTimeG7_p,griddingTime;
172 : casacore::Timer timer_p;
173 : //
174 : //------------------------------------------------------------------------------
175 : //----------------------------Private parts-------------------------------------
176 : //------------------------------------------------------------------------------
177 : //
178 : // private:
179 : protected:
180 : casacore::Vector<casacore::Double> uvwScale_p, offset_p, dphase_p;
181 : casacore::Vector<casacore::Int> chanMap_p, polMap_p;
182 : casacore::Matrix<casacore::Double> spwChanFreq_p, spwChanConjFreq_p;
183 : CFStore convFuncStore_p;
184 : // casacore::Int inc0_p, inc1_p, inc2_p, inc3_p;
185 : casacore::Vector<casacore::Int> inc_p;
186 : casacore::Int* __restrict__ incPtr_p;
187 : casacore::Vector<casacore::Int> cfMap_p, conjCFMap_p;
188 : // VBRow2CFMapType vbRow2CFMap_p;
189 : double paTolerance_p;
190 : casacore::Matrix<casacore::Complex> cached_phaseGrad_p;
191 : casacore::CountedPtr<refim::VB2CFBMap> vb2CFBMap_p;
192 :
193 : void sgrid(casacore::Int& ndim,
194 : casacore::Double* __restrict__ pos,
195 : casacore::Int* __restrict__ loc,
196 : casacore::Int* __restrict__ off,
197 : casacore::Complex& phasor, const casacore::Int& irow,
198 : const casacore::Double* __restrict__ uvw,
199 : const casacore::Double& dphase, const casacore::Double& freq,
200 : const casacore::Double* __restrict__ scale,
201 : const casacore::Double* __restrict__ offset,
202 : const casacore::Float* __restrict__ sampling);
203 :
204 0 : inline casacore::Bool onGrid (const casacore::Int& nx, const casacore::Int& ny,
205 : const casacore::Vector<casacore::Int>& __restrict__ loc,
206 : const casacore::Vector<casacore::Int>& __restrict__ support) __restrict__
207 : {
208 0 : return (((loc(0)-support[0]) >= 0 ) && ((loc(0)+support[0]) < nx) &&
209 0 : ((loc(1)-support[1]) >= 0 ) && ((loc(1)+support[1]) < ny));
210 : };
211 : inline casacore::Bool onGrid (const casacore::Int& nx, const casacore::Int& ny,
212 : const casacore::Int& loc0, const casacore::Int& loc1,
213 : const casacore::Int& support) __restrict__
214 : {
215 : return (((loc0-support) >= 0 ) && ((loc0+support) < nx) &&
216 : ((loc1-support) >= 0 ) && ((loc1+support) < ny));
217 : };
218 :
219 : // casacore::Array assignment operator in CASACore requires lhs.nelements()
220 : // == 0 or lhs.nelements()=rhs.nelements()
221 : // template <class T>
222 : // inline void SETVEC(casacore::Vector<T>& lhs, const casacore::Vector<T>& rhs)
223 : // {lhs.resize(rhs.shape()); lhs = rhs;};
224 :
225 :
226 : //===============================================================================
227 : // CASACORE-LEVEL MATERIAL
228 : //===============================================================================
229 : // Internal methods to address a 4D array. These should ulimately
230 : // moved to a Array4D class in CASACore
231 : //
232 :
233 : // This is called less frequently. Currently once per VisBuffer
234 0 : inline void cacheAxisIncrements(const casacore::Int& n0, const casacore::Int& n1, const casacore::Int& n2, const casacore::Int& n3)
235 : {
236 : // inc0_p=1, inc1_p=inc0_p*n0, inc2_p=inc1_p*n1, inc3_p=inc2_p*n2;(void)n3;
237 0 : inc_p.resize(4);
238 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;
239 : casacore::Bool D;
240 0 : incPtr_p = inc_p.getStorage(D);
241 0 : }
242 0 : inline void cacheAxisIncrements(const casacore::Vector<casacore::Int>& n)
243 0 : {cacheAxisIncrements(n[0],n[1],n[2],n[3]);}
244 :
245 7248095 : inline void cacheAxisIncrements(const casacore::Vector<casacore::Int>& n, casacore::Vector<casacore::Int>& inc)
246 7248095 : {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];}
247 :
248 : inline void cacheAxisIncrements(const casacore::Int n[4], casacore::Int inc[4])
249 : {inc[0]=1; inc[1]=inc[0]*n[0]; inc[2]=inc[1]*n[1]; inc[3]=inc[2]*n[2];(void)n[3];}
250 :
251 : // Version that use internally cached inc_p
252 : // template <class T>
253 0 : inline void addTo4DArray(casacore::DComplex* __restrict__& store, casacore::Int* __restrict__& iPos,
254 : casacore::Complex& nvalue, casacore::Double& wt) __restrict__
255 0 : {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
256 :
257 0 : inline void addTo4DArray(casacore::Complex* __restrict__& store, casacore::Int* __restrict__& iPos,
258 : casacore::Complex& nvalue, casacore::Double& wt) __restrict__
259 0 : {addTo4DArray(store, iPos, incPtr_p, nvalue, wt);}
260 :
261 :
262 : // Version where inc_p is supplied from outside
263 0 : inline void addTo4DArray(casacore::DComplex* __restrict__& store, casacore::Int* __restrict__& iPos,
264 : casacore::Int* __restrict__ inc, casacore::Complex& nvalue, casacore::Double& wt) __restrict__
265 0 : {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*casacore::Complex(wt));}
266 :
267 0 : inline void addTo4DArray(casacore::Complex* __restrict__& store, casacore::Int* __restrict__& iPos,
268 : casacore::Int* __restrict__ inc, casacore::Complex& nvalue, casacore::Double& wt) __restrict__
269 0 : {store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]] += (nvalue*casacore::Complex(wt));}
270 :
271 :
272 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store,
273 : const casacore::Int* __restrict__& iPos,
274 : const casacore::Vector<casacore::Int>& inc)
275 : // __restrict__
276 0 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
277 :
278 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store,
279 : const casacore::Vector<casacore::Int> iPos, const casacore::Vector<casacore::Int>& inc)
280 : // __restrict__
281 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
282 :
283 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex* __restrict__& store,
284 : const casacore::Int* __restrict__& iPos,
285 : 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::Vector<casacore::Int> iPos, const casacore::Vector<casacore::Int>& inc)
291 : // __restrict__
292 : {return store[iPos[0] + iPos[1]*inc[1] + iPos[2]*inc[2] +iPos[3]*inc[3]];};
293 :
294 :
295 : // The following two methods are called in the innermost loop.
296 0 : inline casacore::Complex getFrom4DArray(const casacore::Complex* __restrict__& store, const casacore::Int* __restrict__& iPos)
297 : // __restrict__
298 0 : {return getFrom4DArray(store, iPos, inc_p);}
299 :
300 : inline casacore::DComplex getFrom4DArray(const casacore::DComplex* __restrict__& store, const casacore::Int* __restrict__& iPos)
301 : // __restrict__
302 : {return getFrom4DArray(store, iPos, inc_p);}
303 :
304 : };
305 : }; //# NAMESPACE CASA - END
306 : };
307 : #endif //
|
