Line data Source code
1 : //# SetJyGridFT.cc: Implementation of GridFT class
2 : //# Copyright (C) 2012-2015
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : /*
27 : * SetJyGridFT.cc
28 : *
29 : * Created on: Jun 11, 2012
30 : * Author: kgolap
31 : */
32 : #include <msvis/MSVis/VisibilityIterator2.h>
33 : #include <casacore/casa/Quanta/UnitMap.h>
34 : #include <casacore/casa/Quanta/UnitVal.h>
35 : #include <casacore/measures/Measures/Stokes.h>
36 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
37 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
38 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
39 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
40 : #include <casacore/coordinates/Coordinates/Projection.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <casacore/casa/BasicSL/Constants.h>
43 : #include <casacore/scimath/Mathematics/FFTServer.h>
44 : #include <synthesis/TransformMachines2/SetJyGridFT.h>
45 : #include <casacore/scimath/Mathematics/RigidVector.h>
46 : #include <msvis/MSVis/StokesVector.h>
47 : #include <synthesis/TransformMachines/StokesImageUtil.h>
48 : #include <msvis/MSVis/VisBuffer2.h>
49 : #include <casacore/images/Images/ImageInterface.h>
50 : #include <casacore/images/Images/PagedImage.h>
51 : #include <casacore/casa/Containers/Block.h>
52 : #include <casacore/casa/Containers/Record.h>
53 : #include <casacore/casa/Arrays/ArrayLogical.h>
54 : #include <casacore/casa/Arrays/ArrayMath.h>
55 : #include <casacore/casa/Arrays/Array.h>
56 : #include <casacore/casa/Arrays/MaskedArray.h>
57 : #include <casacore/casa/Arrays/Vector.h>
58 : #include <casacore/casa/Arrays/Matrix.h>
59 : #include <casacore/casa/Arrays/Cube.h>
60 : #include <casacore/casa/Arrays/MatrixIter.h>
61 : #include <casacore/casa/BasicSL/String.h>
62 : #include <casacore/casa/Utilities/Assert.h>
63 : #include <casacore/casa/Exceptions/Error.h>
64 : #include <casacore/lattices/Lattices/ArrayLattice.h>
65 : #include <casacore/measures/Measures/UVWMachine.h>
66 : #include <casacore/lattices/Lattices/SubLattice.h>
67 : #include <casacore/lattices/LRegions/LCBox.h>
68 : #include <casacore/lattices/Lattices/LatticeCache.h>
69 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
70 : #include <casacore/lattices/Lattices/LatticeIterator.h>
71 : #include <casacore/lattices/Lattices/LatticeStepper.h>
72 : #include <casacore/scimath/Mathematics/ConvolveGridder.h>
73 : #include <casacore/casa/Utilities/CompositeNumber.h>
74 : #include <casacore/casa/OS/Timer.h>
75 : #include <sstream>
76 : #ifdef _OPENMP
77 : #include <omp.h>
78 : #endif
79 :
80 : using namespace casacore;
81 : namespace casa { //# NAMESPACE CASA - BEGIN
82 : using namespace casacore;
83 : // using namespace casa::async;
84 : namespace refim {//# namespace for imaging refactor
85 :
86 : using namespace casacore;
87 : using namespace casa;
88 : using namespace casacore;
89 : using namespace casa::refim;
90 :
91 : #define NEED_UNDERSCORES
92 : #if defined(NEED_UNDERSCORES)
93 : #define sectdgridjy sectdgridjy_
94 : #endif
95 :
96 : extern "C" {
97 : void sectdgridjy(Complex*,
98 : const Int*,
99 : const Int*,
100 : const Int*,
101 : const Int*,
102 : const Int*,
103 : const Complex*,
104 : const Int*,
105 : const Int*,
106 : const Int *,
107 : const Int *,
108 : //support
109 : const Int*,
110 : const Int*,
111 : const Double*,
112 : const Int*,
113 : const Int*,
114 : //rbeg, rend, loc, off, phasor
115 : const Int*,
116 : const Int*,
117 : const Int*,
118 : const Int*,
119 : const Complex*,
120 : const Double*);
121 : }
122 0 : SetJyGridFT::SetJyGridFT(Long icachesize, Int itilesize, String iconvType,
123 : MPosition mLocation, MDirection mTangent, Float padding,
124 0 : Bool usezero, Bool useDoublePrec,const Vector<Double>& frequencyscale, const Vector<Double>& scale):GridFT(icachesize, itilesize, iconvType,
125 : mLocation, mTangent, padding, usezero, useDoublePrec), freqscale_p(frequencyscale),
126 0 : scale_p(scale) {
127 :
128 0 : machineName_p="SetJyGridFT";
129 :
130 :
131 0 : }
132 :
133 0 : SetJyGridFT::~SetJyGridFT(){
134 :
135 0 : }
136 :
137 0 : SetJyGridFT::SetJyGridFT(const RecordInterface& stateRec)
138 0 : : GridFT()
139 : {
140 : // Construct from the input state record
141 0 : String error;
142 0 : if (!fromRecord(error, stateRec)) {
143 0 : throw (AipsError("Failed to create gridder: " + error));
144 : };
145 0 : }
146 0 : SetJyGridFT::SetJyGridFT(const SetJyGridFT& other) : GridFT()
147 : {
148 0 : machineName_p=("SetJyGridFT");
149 0 : operator=(other);
150 0 : }
151 :
152 0 : SetJyGridFT& SetJyGridFT::operator=(const SetJyGridFT& other)
153 : {
154 0 : if(this!=&other) {
155 : //Do the base parameters
156 0 : GridFT::operator=(other);
157 0 : freqscale_p.resize();
158 0 : freqscale_p=other.freqscale_p;
159 0 : scale_p.resize();
160 0 : scale_p=other.scale_p;
161 : }
162 0 : return *this;
163 : }
164 :
165 0 : FTMachine* SetJyGridFT::cloneFTM(){
166 0 : return new SetJyGridFT(*this);
167 : }
168 0 : void SetJyGridFT::setScale(const Vector<Double>& freq, const Vector<Double>& scale){
169 0 : freqscale_p.resize();
170 0 : freqscale_p=freq;
171 0 : scale_p.resize();
172 0 : scale_p=scale;
173 :
174 0 : }
175 0 : String SetJyGridFT::name() const{
176 :
177 0 : return String("SetJyGridFT");
178 : }
179 0 : void SetJyGridFT::initializeToVis(ImageInterface<Complex>& image,
180 : const vi::VisBuffer2& vb){
181 0 : setFreqInterpolation(String("nearest"));
182 0 : GridFT::initializeToVis(image, vb);
183 0 : }
184 :
185 0 : Bool SetJyGridFT::toRecord(String& error,
186 : RecordInterface& outRec, Bool withImage, const String diskimage)
187 : {
188 0 : if(!GridFT::toRecord(error, outRec, withImage, diskimage))
189 0 : return false;
190 0 : outRec.define("freqscale", freqscale_p);
191 0 : outRec.define("scaleamp", scale_p);
192 0 : return true;
193 : }
194 0 : Bool SetJyGridFT::fromRecord(String& error,
195 : const RecordInterface& inRec)
196 : {
197 :
198 0 : if(!GridFT::fromRecord(error, inRec))
199 0 : return false;
200 0 : freqscale_p.resize();
201 0 : inRec.get("freqscale", freqscale_p);
202 0 : scale_p.resize();
203 0 : inRec.get("scaleamp", scale_p);
204 0 : machineName_p="SetJyGridFT";
205 0 : return true;
206 : }
207 0 : void SetJyGridFT::get(vi::VisBuffer2& vb, Int row){
208 : //Did somebody really want this version.
209 0 : if(scale_p.nelements()==0)
210 0 : return GridFT::get(vb,row);
211 :
212 0 : gridOk(gridder->cSupport()(0));
213 : // If row is -1 then we pass through all rows
214 : Int startRow, endRow, nRow;
215 0 : if (row < 0) {
216 0 : nRow=vb.nRows();
217 0 : startRow=0;
218 0 : endRow=nRow-1;
219 : //unnecessary zeroing
220 : //vb.modelVisCube()=Complex(0.0,0.0);
221 : } else {
222 0 : nRow=1;
223 0 : startRow=row;
224 0 : endRow=row;
225 : //unnecessary zeroing
226 : //vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
227 : }
228 :
229 : // Get the uvws in a form that Fortran can use
230 0 : Matrix<Double> uvw(negateUV(vb));
231 0 : Vector<Double> dphase(vb.nRows());
232 0 : dphase=0.0;
233 0 : rotateUVW(uvw, dphase, vb);
234 0 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
235 :
236 :
237 :
238 : //Check if ms has changed then cache new spw and chan selection
239 : //if(vb.newMS())
240 : // matchAllSpwChans(vb);
241 :
242 :
243 : //Here we redo the match or use previous match
244 :
245 : //Channel matching for the actual spectral window of buffer
246 : // if(doConversion_p[vb.spectralWindow()]){
247 0 : matchChannel(vb);
248 : // }
249 : // else{
250 : // chanMap.resize();
251 : // chanMap=multiChanMap_p[vb.spectralWindow()];
252 : // }
253 :
254 : //cerr << "chanMap " << chanMap << endl;
255 : //No point in reading data if its not matching in frequency
256 0 : if(max(chanMap)==-1)
257 0 : return;
258 :
259 0 : Cube<Complex> data;
260 0 : Cube<Int> flags;
261 0 : getInterpolateArrays(vb, data, flags);
262 :
263 0 : Vector<Double> lsrfreq;
264 : //Bool convert;
265 :
266 0 : InterpolateArray1D<Double, Double>::InterpolationMethod meth= InterpolateArray1D<Double, Double>::nearestNeighbour;
267 0 : if(freqscale_p.nelements() > 2)
268 0 : meth= InterpolateArray1D<Double, Double>::linear;
269 :
270 0 : lsrfreq=vb.getFrequencies(0,MFrequency::LSRK);
271 : //vb.lsrFrequency(vb.spectralWindow(), lsrfreq, convert);
272 0 : interpscale_p.resize(lsrfreq.nelements());
273 : InterpolateArray1D<Double,Double>::
274 0 : interpolate(interpscale_p,lsrfreq, freqscale_p, scale_p,meth);
275 :
276 : // IPosition s(data.shape());
277 0 : Int nvp=data.shape()(0);
278 0 : Int nvc=data.shape()(1);
279 0 : Int nvisrow=data.shape()(2);
280 :
281 :
282 : //cerr << "get flags " << min(flags) << " " << max(flags) << endl;
283 : Complex *datStorage;
284 : Bool isCopy;
285 0 : datStorage=data.getStorage(isCopy);
286 :
287 : ///
288 0 : Cube<Int> loc(2, nvc, nvisrow);
289 0 : Cube<Int> off(2, nvc, nRow);
290 0 : Matrix<Complex> phasor(nvc, nRow);
291 0 : Int csamp=gridder->cSampling();
292 : Bool delphase;
293 : Bool del;
294 0 : Complex * phasorstor=phasor.getStorage(delphase);
295 0 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
296 0 : const Double * scalestor=uvScale.getStorage(del);
297 0 : const Double * offsetstor=uvOffset.getStorage(del);
298 0 : const Double* uvstor= uvw.getStorage(del);
299 0 : Int * locstor=loc.getStorage(del);
300 0 : Int * offstor=off.getStorage(del);
301 0 : const Double *dpstor=dphase.getStorage(del);
302 0 : const Double *freqscalestor=interpscale_p.getStorage(del);
303 : //Vector<Double> f1=interpVisFreq_p.copy();
304 0 : Int nvchan=nvc;
305 : Int irow;
306 0 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvchan, scalestor, offsetstor, csamp, phasorstor, uvstor, locstor, offstor, dpstor) shared(startRow, endRow) num_threads(4)
307 : {
308 : #pragma omp for
309 : for (irow=startRow; irow<=endRow; ++irow){
310 : locateuvw(uvstor,dpstor, visfreqstor, nvchan, scalestor, offsetstor, csamp,
311 : locstor,
312 : offstor, phasorstor, irow);
313 : }
314 :
315 : }//end pragma parallel
316 :
317 0 : Int rbeg=startRow+1;
318 0 : Int rend=endRow+1;
319 :
320 :
321 0 : Vector<Int> rowFlags(vb.nRows());
322 0 : rowFlags=0;
323 0 : rowFlags(vb.flagRow())=true;
324 0 : if(!usezero_p) {
325 0 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
326 0 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
327 : }
328 : }
329 :
330 :
331 : //cerr <<"offset " << min(off) << " " <<max(off) << " length " << gridder->cFunction().shape() << endl;
332 :
333 : {
334 : Bool delgrid;
335 0 : const Complex* gridstor=griddedData.getStorage(delgrid);
336 0 : const Double * convfuncstor=(gridder->cFunction()).getStorage(del);
337 :
338 0 : const Int* pmapstor=polMap.getStorage(del);
339 0 : const Int *cmapstor=chanMap.getStorage(del);
340 0 : Int nc=nchan;
341 0 : Int np=npol;
342 0 : Int nxp=nx;
343 0 : Int nyp=ny;
344 0 : Int csupp=gridder->cSupport()(0);
345 0 : const Int * flagstor=flags.getStorage(del);
346 0 : const Int * rowflagstor=rowFlags.getStorage(del);
347 :
348 :
349 0 : Int npart=1;
350 : #ifdef _OPENMP
351 0 : Int nthreads=omp_get_max_threads();
352 0 : if (nthreads >3){
353 0 : npart=4;
354 : }
355 0 : else if(nthreads >1){
356 0 : npart=2;
357 : }
358 : #endif
359 0 : Int ix=0;
360 0 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(datStorage, flagstor, rowflagstor, convfuncstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csupp, nvp, nvc, nvisrow, phasorstor, locstor, offstor, freqscalestor) shared(npart) num_threads(npart)
361 : {
362 : #pragma omp for
363 : for (ix=0; ix< npart; ++ix){
364 : rbeg=ix*(nvisrow/npart)+1;
365 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)+nvisrow%npart-1) ;
366 : //cerr << "rbeg " << rbeg << " rend " << rend << " " << nvisrow << endl;
367 :
368 : sectdgridjy(datStorage,
369 : &nvp,
370 : &nvc,
371 : flagstor,
372 : rowflagstor,
373 : &nvisrow,
374 : gridstor,
375 : &nxp,
376 : &nyp,
377 : &np,
378 : &nc,
379 : &csupp,
380 : &csamp,
381 : convfuncstor,
382 : cmapstor,
383 : pmapstor,
384 : &rbeg, &rend, locstor, offstor, phasorstor, freqscalestor);
385 : }//end pragma parallel
386 : }
387 0 : data.putStorage(datStorage, isCopy);
388 0 : griddedData.freeStorage(gridstor, delgrid);
389 : //cerr << "Get min max " << min(data) << " " << max(data) << endl;
390 : }
391 0 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
392 :
393 :
394 : }
395 :
396 :
397 :
398 : }//END refim namespace
399 :
400 : }//END casa namespace
|
