Line data Source code
1 : //# AWPLPG.h: Implementation for a CPU based gridder for A,W (LPG= LowPerformanceGridder)
2 : //# Copyright (C) 2023
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 General Public License as published by
7 : //# the Free Software Foundation; either version 3 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 General Public
13 : //# License for more details.
14 : //#
15 : //# https://www.gnu.org/licenses/
16 : //#
17 : //# You should have received a copy of the GNU General Public License
18 : //# along with this library; if not, write to the Free Software Foundation,
19 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20 : //#
21 : //# Queries concerning CASA should be submitted at
22 : //# https://help.nrao.edu
23 : //#
24 : //# Postal address: CASA Project Manager
25 : //# National Radio Astronomy Observatory
26 : //# 520 Edgemont Road
27 : //# Charlottesville, VA 22903-2475 USA
28 : //#
29 : //#
30 :
31 : #include <synthesis/TransformMachines2/AWPLPG.h>
32 : #include <synthesis/TransformMachines2/AWConvFuncHolder.h>
33 : #include <synthesis/TransformMachines2/HetArrayConvFunc.h>
34 : #include <synthesis/TransformMachines2/SimplePBConvFunc.h>
35 : #include <msvis/MSVis/VisBuffer2.h>
36 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
37 : #include <casacore/casa/Arrays/ArrayMath.h>
38 : #include <casacore/casa/Arrays/Matrix.h>
39 : #include <casacore/casa/Arrays/Vector.h>
40 :
41 :
42 :
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 : namespace refim {//# namespace for wstatimaging refactor
45 : using namespace casacore;
46 : using namespace casa;
47 : using namespace casa::refim;
48 :
49 :
50 0 : AWPLPG::AWPLPG(SkyJones* sj, const Int nw, Bool dosquint, const Double painc, MPosition mloc, String stokes, const Bool usezero, const Bool useDoublePrec, const casacore::Bool usePointing): MosaicFTNew(sj, mloc, stokes, Long(1000000), 16, usezero, True, False, usePointing), doSquint_p(dosquint), paInc_p(painc), nw_p(nw) {
51 :
52 0 : useDoubleGrid_p=True; //We'll always use double prec for this grid
53 :
54 0 : }
55 0 : AWPLPG::AWPLPG(const AWPLPG& other) : MosaicFTNew(other)
56 : {
57 0 : operator=(other);
58 0 : }
59 :
60 0 : AWPLPG& AWPLPG::operator=(const AWPLPG& other) {
61 0 : if(this!=&other) {
62 :
63 : //Do the base parameters
64 0 : MosaicFTNew::operator=(other);
65 0 : awConvs_p = other.awConvs_p;
66 0 : doSquint_p = other.doSquint_p;
67 0 : paInc_p = other.paInc_p;
68 0 : nw_p = other.nw_p;
69 : }
70 0 : return *this;
71 :
72 : }
73 0 : refim::FTMachine* AWPLPG::cloneFTM(){
74 0 : return new AWPLPG(*this);
75 : }
76 0 : void AWPLPG::init(const vi::VisBuffer2& vb){
77 0 : MosaicFTNew::init(vb);
78 :
79 :
80 : //oversample if image is small
81 : //But not more than 10000 pixels
82 0 : convSampling=(max(nx, ny) < 100) ? 100: Int(ceil(10000.0/max(nx, ny)));
83 :
84 0 : if(convSampling <4)
85 0 : convSampling=4;
86 : //For multiple pa angle reduce mem consumed
87 : //if(doSquint_p)
88 : // convSampling = 2;
89 : // TESTOO
90 : // convSampling = 1;
91 : // TESTOO
92 :
93 0 : CoordinateSystem cs=image->coordinates();
94 :
95 0 : SpectralCoordinate spCS = cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL));
96 : double f1, f2;
97 : { //Lets get the frame to convert to
98 : MFrequency::Types fframe;
99 :
100 0 : Int spw = vb.spectralWindows()(0);
101 0 : MDirection d;
102 0 : cs.directionCoordinate(0).toWorld(d, Vector<Double>(2,0));
103 0 : MPosition p= cs.obsInfo().telescopePosition();
104 0 : MEpoch e = cs.obsInfo().obsDate();
105 0 : fframe=(MFrequency::Types)vb.subtableColumns().spectralWindow().measFreqRef()(spw);
106 0 : spCS.setReferenceConversion(fframe, e, p, d);
107 :
108 0 : }
109 :
110 0 : nchan = image->shape()(3);
111 0 : spCS.toWorld(f1, double(-0.5));
112 0 : spCS.toWorld(f2, double(nchan)-0.5);
113 0 : auto frange=std::make_pair(f1, f2);
114 :
115 0 : if (pbConvFunc_p.null())
116 0 : pbConvFunc_p=new HetArrayConvFunc();
117 0 : awConvs_p=pbConvFunc_p->getAWConvFuncHolder();
118 0 : if(awConvs_p.use_count()==0){
119 0 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
120 0 : awConvs_p=std::make_shared<AWConvFuncHolder>((*image).coordinates(), nx, ny,
121 0 : doSquint_p, paInc_p, observatory, convSampling);
122 0 : vi::VisibilityIterator2 *vi= const_cast<VisibilityIterator2 *>(vb.getVi());
123 :
124 0 : if(sj_p)
125 0 : pbConvFunc_p->setSkyJones(sj_p.get());
126 :
127 0 : if((pbConvFunc_p->getVBUtil()).null()){
128 0 : if(vbutil_p.null()){
129 0 : vbutil_p=new VisBufferUtil(vb);
130 : }
131 0 : pbConvFunc_p->setVBUtil(vbutil_p);
132 : }
133 0 : std::vector<Double> freqs;
134 0 : std::set<Int> fields;
135 0 : std::vector<Double> pAs;
136 0 : Double maxW=0.0;
137 0 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
138 0 : for (vi->origin(); vi->more(); vi->next()) {
139 0 : std::vector<Double> chunkfreq;
140 :
141 0 : fields.insert(vb.fieldId()(0));
142 0 : pbConvFunc_p->findUsefulChannels(chunkfreq, vb, frange);
143 : //cerr << "chunkfreq " << chunkfreq << endl;
144 0 : if (chunkfreq.size() > 0) {
145 : //cerr << "SPW " << vb.spectralWindows()(0) << " freqs " << Vector<Double>(chunkfreq) << endl;
146 0 : std::move(chunkfreq.begin(), chunkfreq.end(), std::back_inserter(freqs));
147 0 : double maxfreqused = *(std::max_element(chunkfreq.begin(), chunkfreq.end()));
148 0 : if (nw_p > 1) {
149 : // maxW=max(maxW, max(abs(vb.uvw().row(2)*max(vb.getFrequencies(0))))/C::c);
150 0 : maxW = max(maxW, max(abs(vb.uvw().row(2) * maxfreqused)) / C::c);
151 : }
152 : }
153 0 : if (doSquint_p)
154 0 : pAs.push_back(getPA(vb));
155 : //if(nw_p > 1)
156 : // maxW=max(maxW, max(abs(vb.uvw().row(2)*max(vb.getFrequencies(0))))/C::c);
157 0 : }
158 : }
159 : ///TESTOO
160 : //Double imMaxW = 0.25 / abs(cs.increment()(0));
161 : //cerr << " maxW " << maxW << " imMaxW " << imMaxW << endl;
162 :
163 : ////
164 : //return vi to origin
165 0 : vi->originChunks(); vi->origin();
166 : //cerr << "FREQS " << Vector<Double>(freqs) << endl;
167 0 : std::sort(freqs.begin(), freqs.end());
168 0 : auto last = std::unique(freqs.begin(), freqs.end());
169 0 : freqs.erase(last, freqs.end());
170 :
171 0 : if(freqs.size()==0){
172 0 : cerr << "No matching frequency in data in freq range of image " +
173 0 : String::toString(f1) + " to " + String::toString(f2)
174 0 : << endl;
175 : //Falling in a gap...channels in image does not match any data used
176 : //for now just calc pb for mid freq
177 0 : freqs.push_back((f1 + f2) / 2.0);
178 : }
179 : // tell holder it is a single field or not
180 0 : (*awConvs_p).setSingleField((fields.size()==1) && (nw_p==1));
181 :
182 :
183 0 : Double paMax=0.0;
184 0 : if(pAs.size()>1){
185 0 : std::sort(pAs.begin(), pAs.end());
186 0 : last=std::unique(pAs.begin(), pAs.end());
187 0 : pAs.erase(last, pAs.end());
188 0 : if(pAs.size()==1)
189 0 : paMax=pAs[0];
190 0 : if(pAs.size() >1){
191 0 : auto [minpa, maxpa]=std::minmax_element(begin(pAs), end(pAs));
192 0 : paMax=max(abs(*minpa), *maxpa);
193 : }
194 : }
195 :
196 :
197 0 : if (nw_p == 0)
198 0 : nw_p = 1;
199 0 : Vector<Double> wVals(nw_p,0);
200 0 : if(nw_p >1){
201 0 : Double st=maxW/(Double(nw_p-1)*Double(nw_p-1));
202 0 : for (int k=0; k <nw_p; ++k)
203 0 : wVals[k]=Double(k*k)*st;
204 : }
205 0 : (*awConvs_p).addConvFunc(Vector<Double>(freqs), wVals, paMax);
206 :
207 0 : pbConvFunc_p->setAWConvFuncHolder(awConvs_p);
208 :
209 0 : }
210 :
211 0 : }
212 :
213 0 : void AWPLPG::findConvFunction(const ImageInterface<Complex>& iimage, const vi::VisBuffer2& vb, const Matrix<Double>& rotuvw, const bool ispsf ){
214 : //
215 : // pbConvFunc_p.phasegradient=
216 : //double time0=omp_get_wtime();
217 : //Complex *oWgtPtr, *oConPtr;
218 : //Bool isCopy;
219 : //if(convFunc.size()==0 || (convFunc.shape() != awConvs_p->getConvFunc().shape())){
220 : // convFunc.resize(awConvs_p->getConvFunc().shape());
221 : // weightConvFunc_p.resize(awConvs_p->getWeightConvFunc().shape());
222 :
223 : //}
224 : //oWgtPtr=awConvs_p->getWeightConvFunc().getStorage(isCopy);
225 : //oConPtr=awConvs_p->getConvFunc().getStorage(isCopy);
226 : //convFunc.resize();
227 : //convFunc=(awConvs_p->getConvFunc());
228 : //Bool isCopy1, isCopy2;
229 : //cerr << "SIZEOF " << sizeof convFunc << " size elem wise " << convFunc.nelements() << endl;
230 : //Complex* convFuncPtr=convFunc.getStorage(isCopy1);
231 : // Complex* wgtFuncPtr=weightConvFunc_p.getStorage(isCopy2);
232 : //weightConvFunc_p.resize();
233 : //weightConvFunc_p=(awConvs_p->getWeightConvFunc());
234 : //std::memcpy(convFuncPtr, oConPtr, sizeof(Complex)*convFunc.nelements());
235 : //std::memcpy(wgtFuncPtr, oWgtPtr, sizeof(Complex)*weightConvFunc_p.nelements());
236 : //convFunc.putStorage(convFuncPtr, isCopy1);
237 : //weightConvFunc_p.putStorage(wgtFuncPtr, isCopy2);
238 : /*{
239 : ////TESTOO
240 : IPosition elshp = convFunc.shape().getFirst(4);
241 : IPosition elblc(5, 0);
242 :
243 : IPosition eltrc = convFunc.shape()-1;
244 : elblc[4] = eltrc[4];
245 : CoordinateSystem csysA = iimage.coordinates();
246 : Vector<Int> stoks(4);
247 : stoks(0) = Stokes::RR;
248 : stoks(1) = Stokes::RL;
249 : stoks(2) = Stokes::LR;
250 : stoks(3) = Stokes::LL;
251 : StokesCoordinate stokey(stoks);
252 : csysA.replaceCoordinate(stokey, 1);
253 : PagedImage<Complex> lastplane(elshp, csysA, "COOBOO");
254 : lastplane.put(convFunc(elblc, eltrc).nonDegenerate());
255 : PagedImage<Complex> lastplaneW(elshp, csysA, "WOOBOO");
256 : lastplaneW.put(weightConvFunc_p(elblc, eltrc).nonDegenerate());
257 : //////
258 : } */
259 :
260 :
261 0 : awConvs_p->getConvFuncs(convPolMap_p, convChanMap_p, convRowMap_p, convFunc,
262 0 : weightConvFunc_p, vb, rotuvw, interpVisFreq_p, toVis_p, ispsf);
263 :
264 : //double time1=omp_get_wtime();
265 : //cerr << " assign time " << time1-time0 << endl;
266 0 : convSizePlanes_p.resize();
267 0 : convSizePlanes_p = awConvs_p->getConvSizes();
268 0 : convSupportPlanes_p.resize();
269 0 : convSupportPlanes_p = awConvs_p->getConvSupports();
270 :
271 : //awConvs_p->getConvIndices(convPolMap_p, convChanMap_p, convRowMap_p, vb, rotuvw);
272 : //cerr << "min max convrowmap " << min(convRowMap_p) << " " << max(convRowMap_p) << " supp " << max(convSupportPlanes_p) << " csize " << max(convSizePlanes_p) << " convchanmap "<< min(convChanMap_p) << " " << max(convChanMap_p) << " convsamp " << convSampling << endl;
273 : //cerr << "LENGTHS bef" << convRowMap_p.size() << " " << convChanMap_p.size() << " " << convPolMap_p.size() << endl;
274 0 : std::vector<Int> pmapused = convPolMap_p.tovector();
275 : {
276 0 : std::sort(pmapused.begin(), pmapused.end());
277 0 : auto last = std::unique(pmapused.begin(), pmapused.end());
278 0 : pmapused.erase(last, pmapused.end());
279 : }
280 0 : std::vector<Int> cmapused=convChanMap_p.tovector();
281 : {
282 0 : std::sort(cmapused.begin(), cmapused.end());
283 0 : auto last = std::unique(cmapused.begin(), cmapused.end());
284 0 : cmapused.erase(last, cmapused.end());
285 : }
286 0 : std::vector<Int> rmapused=abs(convRowMap_p).tovector();
287 : {
288 0 : std::sort(rmapused.begin(), rmapused.end());
289 0 : auto last = std::unique(rmapused.begin(), rmapused.end());
290 0 : rmapused.erase(last, rmapused.end());
291 : }
292 :
293 0 : pbConvFunc_p->rephaseConvFunc(
294 0 : iimage, vb, convSampling, convFunc, weightConvFunc_p, pmapused,
295 0 : cmapused, rmapused, MVDirection(-(movingDirShift_p.getAngle())),
296 0 : fixMovingSource_p);
297 0 : convSupport =max(convSupportPlanes_p);
298 0 : convSize = max(convSizePlanes_p);
299 0 : if(convFunc.nelements()==0){
300 0 : convSupport = 0;
301 0 : convSize = 0;
302 : }
303 : //cerr << "csup " << convSupport << " csize "<< convSize << " csamp " << convSampling << endl;
304 :
305 :
306 0 : }
307 :
308 : /////==============================================
309 : //// some fortran defn
310 : #define NEED_UNDERSCORES
311 : #if defined(NEED_UNDERSCORES)
312 : #define sectgmosd3 sectgmosd3_
313 : #define sectdmos3 sectdmos3_
314 : #define gmoswgtd2 gmoswgtd2_
315 : #define locuvw locuvw_
316 : #endif
317 :
318 :
319 : /////==============================================
320 : //// some fortran defn
321 : #define NEED_UNDERSCORES
322 : #if defined(NEED_UNDERSCORES)
323 : #define sectgmosd3 sectgmosd3_
324 : #define sectdmos3 sectdmos3_
325 : #define gmoswgtd2 gmoswgtd2_
326 : #define locuvw locuvw_
327 : #endif
328 :
329 :
330 : extern "C" {
331 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
332 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
333 : void gmoswgtd2(const Int*/*nvispol*/, const Int*/*nvischan*/,
334 : const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/,
335 : const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/,
336 : const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/,
337 : const Int*/*chanmap*/, const Int*/*polmap*/,
338 : DComplex* /*weightgrid*/, Double* /*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/,
339 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
340 : const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/,
341 : const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
342 :
343 :
344 : void sectgmosd3(const Complex* /*values*/,
345 : Int* /*nvispol*/, Int* /*nvischan*/,
346 : Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
347 : Int* /* nrow*/, DComplex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan */,
348 : const Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
349 : const Int*/*chanmap*/, const Int*/*polmap*/,
350 : Double*/*sumwgt*/, const Int*/*convplanemap*/,
351 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
352 : Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
353 : const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/,
354 : const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);
355 :
356 : void sectdmos3(Complex*,
357 : Int*,
358 : Int*,
359 : const Int*,
360 : const Int*,
361 : Int*,
362 : const Complex*,
363 : Int*,
364 : Int*,
365 : Int *,
366 : Int *,
367 : const Int*, //support
368 : Int*,
369 : Int*,
370 : const Complex*,
371 : const Int*,
372 : const Int*,
373 : const Int*,
374 : const Int*,
375 : const Int*,
376 : Int*, Int*, Int*,
377 : //rbeg
378 : const Int*,
379 : const Int*,
380 : const Int*,
381 : const Int*,
382 : const Complex*);
383 :
384 :
385 :
386 : }
387 :
388 :
389 :
390 :
391 :
392 :
393 :
394 :
395 :
396 :
397 :
398 : //===================================================
399 : /* void AWPLPG::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
400 : FTMachine::Type type)
401 : {
402 :
403 :
404 :
405 :
406 : Timer tim;
407 : tim.mark();
408 :
409 : matchChannel(vb);
410 :
411 :
412 : //cerr << "CHANMAP " << chanMap << endl;
413 : //No point in reading data if its not matching in frequency
414 : if(max(chanMap)==-1)
415 : return;
416 :
417 : //const Matrix<Float> *imagingweight;
418 : //imagingweight=&(vb.imagingWeight());
419 : Matrix<Float> imagingweight;
420 : getImagingWeight(imagingweight, vb);
421 :
422 : if(dopsf) type=FTMachine::PSF;
423 :
424 : Cube<Complex> data;
425 : //Fortran gridder need the flag as ints
426 : Cube<Int> flags;
427 : Matrix<Float> elWeight;
428 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
429 :
430 :
431 :
432 : Bool iswgtCopy;
433 : const Float *wgtStorage;
434 : wgtStorage=elWeight.getStorage(iswgtCopy);
435 :
436 :
437 :
438 :
439 : Bool isCopy;
440 : const Complex *datStorage=0;
441 :
442 : // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << " " << wgtStorage<< endl;
443 : if(!dopsf)
444 : datStorage=data.getStorage(isCopy);
445 :
446 :
447 : // If row is -1 then we pass through all rows
448 : Int startRow, endRow, nRow;
449 : if (row==-1) {
450 : nRow=vb.nRows();
451 : startRow=0;
452 : endRow=nRow-1;
453 : } else {
454 : nRow=1;
455 : startRow=row;
456 : endRow=row;
457 : }
458 :
459 : // Get the uvws in a form that Fortran can use and do that
460 : // necessary phase rotation.
461 : Matrix<Double> uvw(negateUV(vb));
462 : Vector<Double> dphase(vb.nRows());
463 : dphase=0.0;
464 :
465 : doUVWRotation_p=true;
466 : girarUVW(uvw, dphase, vb);
467 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
468 : // This needs to be after the interp to get the interpolated channels
469 : //Also has to be after rotateuvw in case tracking is on
470 : findConvFunction(*image, vb);
471 : //cerr << "Put convsup " << convSupport << " max min convFunc " << max(convFunc) << " " << min(convFunc) << " " << max(weightConvFunc_p) << min(weightConvFunc_p) << "SHP " << convFunc.shape() << " " << weightConvFunc_p.shape() << endl;
472 : //cerr << "convRowMap " << convRowMap_p << " " << convChanMap_p << " " << convPolMap_p << endl;
473 : //nothing to grid here as the pointing resulted in a zero support convfunc
474 : if(convSupport <= 0)
475 : return;
476 :
477 : // Get the pointing positions. This can easily consume a lot
478 : // of time thus we are for now assuming a field per
479 : // vb chunk...need to change that accordingly if we start using
480 : // multiple pointings per vb.
481 : //Warning
482 :
483 : // Take care of translation of Bools to Integer
484 : Int idopsf=0;
485 : if(dopsf) idopsf=1;
486 :
487 :
488 : Vector<Int> rowFlags(vb.nRows());
489 : rowFlags=0;
490 : rowFlags(vb.flagRow())=true;
491 : if(!usezero_p) {
492 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
493 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
494 : }
495 : }
496 :
497 :
498 :
499 : //cerr << "convSamp " << convSampling << " convsupp " << convSupport << " consize " << convSize << " convFunc " << convFunc.shape() << endl;
500 : //TESTOO
501 :
502 : //TESTOO
503 :
504 : //Tell the gridder to grid the weights too ...need to do that once only
505 : //Int doWeightGridding=1;
506 : //if(doneWeightImage_p)
507 : // doWeightGridding=-1;
508 : Bool del;
509 : // IPosition s(flags.shape());
510 : const IPosition& fs=flags.shape();
511 : //cerr << "flags shape " << fs << endl;
512 : std::vector<Int>s(fs.begin(), fs.end());
513 : Int nvp=s[0];
514 : Int nvc=s[1];
515 : Int nvisrow=s[2];
516 : Int csamp=convSampling;
517 : Bool uvwcopy;
518 : const Double *uvwstor=uvw.getStorage(uvwcopy);
519 : Bool gridcopy;
520 : Bool convcopy;
521 : Bool wconvcopy;
522 : const Complex *convstor=convFunc.getStorage(convcopy);
523 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
524 : Int nPolConv=convFunc.shape()[2];
525 : Int nChanConv=convFunc.shape()[3];
526 : Int nConvFunc=convFunc.shape()(4);
527 : Bool weightcopy;
528 : ////////
529 : Cube<Int> loc(2, nvc, nRow);
530 : Cube<Int> off(2, nvc, nRow);
531 : Matrix<Complex> phasor(nvc, nRow);
532 : Bool delphase;
533 : Complex * phasorstor=phasor.getStorage(delphase);
534 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
535 : const Double * scalestor=uvScale.getStorage(del);
536 : const Double * offsetstor=uvOffset.getStorage(del);
537 : Int * locstor=loc.getStorage(del);
538 : Int * offstor=off.getStorage(del);
539 : const Double *dpstor=dphase.getStorage(del);
540 : Int irow;
541 : Int nth=1;
542 : #ifdef _OPENMP
543 : if(numthreads_p >0){
544 : nth=min(numthreads_p, omp_get_max_threads());
545 : }
546 : else{
547 : nth= omp_get_max_threads();
548 : }
549 : //nth=min(4,nth);
550 : #endif
551 : Double cinv=Double(1.0)/C::c;
552 :
553 : Int dow=0;
554 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)
555 : {
556 : #pragma omp for
557 : for (irow=startRow; irow<=endRow;irow++){
558 :
559 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
560 : }
561 :
562 : }//end pragma parallel
563 :
564 :
565 :
566 : timemass_p +=tim.real();
567 : Int ixsub, iysub, icounter;
568 : ixsub=1;
569 : iysub=1;
570 : //////@@@@@@@@@@@@@DEBUGGING
571 : //nth=1;
572 : ////////@@@@@@@@@@@@@
573 : if (nth >3){
574 : ixsub=8;
575 : iysub=8;
576 : }
577 : else if(nth >1){
578 : ixsub=2;
579 : iysub=2;
580 : }
581 : Int rbeg=startRow+1;
582 : Int rend=endRow+1;
583 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
584 : Vector<Double *> swgtptr(ixsub*iysub);
585 : Vector<Bool> swgtdel(ixsub*iysub);
586 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
587 : sumwgt[icounter].resize(sumWeight.shape());
588 : sumwgt[icounter].set(0.0);
589 : swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
590 : }
591 : //cerr << "done thread " << doneThreadPartition_p << " " << ixsub*iysub << endl;
592 : if(doneThreadPartition_p < 0){
593 : xsect_p.resize(ixsub*iysub);
594 : ysect_p.resize(ixsub*iysub);
595 : nxsect_p.resize(ixsub*iysub);
596 : nysect_p.resize(ixsub*iysub);
597 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
598 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
599 : }
600 : }
601 : Vector<Int> xsect, ysect, nxsect, nysect;
602 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
603 : //cerr << xsect.shape() << " " << xsect << endl;
604 : const Int* pmapstor=polMap.getStorage(del);
605 : const Int* cmapstor=chanMap.getStorage(del);
606 : // Dummy sumwt for gridweight part
607 : Matrix<Double> dumSumWeight(npol, nchan);
608 : dumSumWeight=sumWeight;
609 : Bool isDSWC;
610 : Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
611 : Int nc=nchan;
612 : Int np=npol;
613 : Int nxp=nx;
614 : Int nyp=ny;
615 : Int csize=convSize;
616 : const Int * flagstor=flags.getStorage(del);
617 : const Int * rowflagstor=rowFlags.getStorage(del);
618 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
619 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
620 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
621 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
622 : ///
623 :
624 :
625 : ////////
626 : tim.mark();
627 :
628 : // if(useDoubleGrid_p) { //always using double prec here
629 : {
630 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
631 :
632 : #pragma omp parallel default(none) private(icounter, del) firstprivate(idopsf, datStorage, wgtStorage, flagstor, rowflagstor, convstor, wconvstor, pmapstor, cmapstor, gridstor, convsupportstor, nxp, nyp, np, nc,ixsub, iysub, rend, rbeg, csamp, csize, nvp, nvc, nvisrow, phasorstor, locstor, offstor, convrowmapstor, convchanmapstor, convpolmapstor, nPolConv, nChanConv, nConvFunc,xsect, ysect, nxsect, nysect) shared(swgtptr)
633 : {
634 : #pragma omp for schedule(dynamic)
635 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
636 : Int x0=xsect(icounter);
637 : Int y0=ysect(icounter);
638 : Int nxsub=nxsect(icounter);
639 : Int nysub=nysect(icounter);
640 :
641 :
642 : sectgmosd3(datStorage,
643 : &nvp,
644 : &nvc,
645 : &idopsf,
646 : flagstor,
647 : rowflagstor,
648 : wgtStorage,
649 : &nvisrow,
650 : gridstor,
651 : &nxp,
652 : &nyp,
653 : &np,
654 : &nc,
655 : convsupportstor,
656 : &csize,
657 : &csamp,
658 : convstor,
659 : cmapstor,
660 : pmapstor,
661 : swgtptr[icounter],
662 : convrowmapstor,
663 : convchanmapstor,
664 : convpolmapstor,
665 : &nConvFunc, &nChanConv, &nPolConv,
666 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
667 : phasorstor
668 : );
669 : }
670 : }//end pragma parallel
671 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
672 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
673 : sumWeight=sumWeight+sumwgt[icounter];
674 : }
675 :
676 : //cerr << "SUMWEIG " << sumWeight << endl;
677 : griddedData2.putStorage(gridstor, gridcopy);
678 : if(dopsf && (nth >4))
679 : tweakGridSector(nx, ny, ixsub, iysub);
680 : timegrid_p+=tim.real();
681 : tim.mark();
682 : if(!doneWeightImage_p){
683 : //This can be parallelized by making copy of the central part of the griddedWeight
684 : //and adding it after dooing the gridding
685 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
686 : gmoswgtd2(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
687 : &nxp, &nyp, &np, &nc, convsupportstor, &csize, &csamp,
688 : cmapstor, pmapstor,
689 : gridwgtstor, dsumwtstor, wconvstor, convrowmapstor,
690 : convchanmapstor, convpolmapstor,
691 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
692 : &rend, locstor, offstor, phasorstor);
693 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
694 :
695 : }
696 : timemass_p+=tim.real();
697 : }
698 :
699 : convFunc.freeStorage(convstor, convcopy);
700 : weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
701 : dumSumWeight.putStorage(dsumwtstor, isDSWC);
702 : //cerr << "dumSumwe " << dumSumWeight << endl;
703 : uvw.freeStorage(uvwstor, uvwcopy);
704 : if(!dopsf)
705 : data.freeStorage(datStorage, isCopy);
706 :
707 : elWeight.freeStorage(wgtStorage,iswgtCopy);
708 :
709 :
710 :
711 :
712 : }
713 : */
714 :
715 :
716 : /*
717 : void AWPLPG::gridImgWeights(const vi::VisBuffer2& vb){
718 :
719 : if(doneWeightImage_p)
720 : return;
721 : matchChannel(vb);
722 :
723 :
724 : //cerr << "CHANMAP " << chanMap << endl;
725 : //No point in reading data if its not matching in frequency
726 : if(max(chanMap)==-1)
727 : return;
728 :
729 : Int startRow, endRow, nRow;
730 : nRow=vb.nRows();
731 : startRow=0;
732 : endRow=nRow-1;
733 :
734 :
735 : //const Matrix<Float> *imagingweight;
736 : //imagingweight=&(vb.imagingWeight());
737 : Matrix<Float> imagingweight;
738 : getImagingWeight(imagingweight, vb);
739 :
740 :
741 : Cube<Complex> data;
742 : //Fortran gridder need the flag as ints
743 : Cube<Int> flags;
744 : Matrix<Float> elWeight;
745 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, FTMachine::PSF);
746 :
747 :
748 :
749 : Bool iswgtCopy;
750 : const Float *wgtStorage;
751 : wgtStorage=elWeight.getStorage(iswgtCopy);
752 : Bool issumWgtCopy;
753 : Double* sumwgtstor=sumWeight.getStorage(issumWgtCopy);
754 :
755 :
756 :
757 : // Get the uvws in a form that Fortran can use and do that
758 : // necessary phase rotation.
759 : Matrix<Double> uvw(negateUV(vb));
760 : Vector<Double> dphase(vb.nRows());
761 : dphase=0.0;
762 :
763 : doUVWRotation_p=true;
764 : girarUVW(uvw, dphase, vb);
765 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
766 : // This needs to be after the interp to get the interpolated channels
767 : //Also has to be after rotateuvw in case tracking is on
768 : findConvFunction(*image, vb);
769 : //nothing to grid here as the pointing resulted in a zero support convfunc
770 : if(convSupport <= 0)
771 : return;
772 :
773 : Bool del;
774 :
775 : const Int* pmapstor=polMap.getStorage(del);
776 : const Int* cmapstor=chanMap.getStorage(del);
777 :
778 : Vector<Int> rowFlags(vb.nRows());
779 : rowFlags=0;
780 : rowFlags(vb.flagRow())=true;
781 : if(!usezero_p) {
782 : for (uInt rownr=0; rownr< vb.nRows(); rownr++) {
783 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
784 : }
785 : }
786 :
787 : //Fortran indexing
788 :
789 : Int rbeg=1;
790 : Int rend=vb.nRows();
791 :
792 : const Int * flagstor=flags.getStorage(del);
793 : const Int * rowflagstor=rowFlags.getStorage(del);
794 :
795 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
796 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
797 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
798 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
799 : //Tell the gridder to grid the weights too ...need to do that once only
800 : //Int doWeightGridding=1;
801 : //if(doneWeightImage_p)
802 : // doWeightGridding=-1;
803 : // IPosition s(flags.shape());
804 : const IPosition& fs=flags.shape();
805 : //cerr << "flags shape " << fs << endl;
806 : std::vector<Int>s(fs.begin(), fs.end());
807 : Int nvp=s[0];
808 : Int nvc=s[1];
809 : Int nvisrow=s[2];
810 : Int csamp=convSampling;
811 : Bool uvwcopy;
812 : const Double *uvwstor=uvw.getStorage(uvwcopy);
813 : Bool gridcopy;
814 : Bool convcopy;
815 : Bool wconvcopy;
816 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
817 : Int nPolConv=convFunc.shape()[2];
818 : Int nChanConv=convFunc.shape()[3];
819 : Int nConvFunc=convFunc.shape()(4);
820 : Bool weightcopy;
821 : ////////@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
822 : Cube<Int> loc(2, nvc, vb.nRows());
823 : Cube<Int> off(2, nvc, vb.nRows());
824 : Matrix<Complex> phasor(nvc, vb.nRows());
825 : Bool delphase;
826 : Complex * phasorstor=phasor.getStorage(delphase);
827 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
828 : const Double * scalestor=uvScale.getStorage(del);
829 : const Double * offsetstor=uvOffset.getStorage(del);
830 : Int * locstor=loc.getStorage(del);
831 : Int * offstor=off.getStorage(del);
832 : const Double *dpstor=dphase.getStorage(del);
833 :
834 : Int irow;
835 : Int nth=1;
836 : #ifdef _OPENMP
837 : if(numthreads_p >0){
838 : nth=min(numthreads_p, omp_get_max_threads());
839 : }
840 : else{
841 : nth= omp_get_max_threads();
842 : }
843 : //nth=min(4,nth);
844 : #endif
845 :
846 : Double cinv=Double(1.0)/C::c;
847 :
848 : Int dow=0;
849 :
850 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)
851 : {
852 : #pragma omp for
853 : for (irow=startRow; irow<=endRow;irow++){
854 :
855 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
856 : }
857 :
858 : }//end pragma parallel
859 :
860 :
861 :
862 : //always using double prec in this gridder
863 : // if(useDoubleGrid_p) {
864 : {
865 : //This can be parallelized by making copy of the central part of the griddedWeight
866 : //and adding it after dooing the gridding
867 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
868 : gmoswgtd2(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
869 : &nx, &ny, &npol, &nchan, convsupportstor, &convSize, &convSampling,
870 : cmapstor, pmapstor,
871 : gridwgtstor, sumwgtstor, wconvstor, convrowmapstor,
872 : convchanmapstor, convpolmapstor,
873 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
874 : &rend, locstor, offstor, phasorstor);
875 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
876 :
877 :
878 :
879 :
880 :
881 : }
882 :
883 : sumWeight.putStorage(sumwgtstor, issumWgtCopy);
884 : elWeight.freeStorage(wgtStorage,iswgtCopy);
885 :
886 : }
887 : */
888 : /*
889 : void AWPLPG::get(vi::VisBuffer2& vb, Int row)
890 : {
891 :
892 :
893 :
894 : // If row is -1 then we pass through all rows
895 : Int startRow, endRow, nRow;
896 : if (row==-1) {
897 : nRow=vb.nRows();
898 : startRow=0;
899 : endRow=nRow-1;
900 : // vb.modelVisCube()=Complex(0.0,0.0);
901 : } else {
902 : nRow=1;
903 : startRow=row;
904 : endRow=row;
905 : // vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
906 : }
907 :
908 :
909 :
910 :
911 : matchChannel(vb);
912 :
913 : //No point in reading data if its not matching in frequency
914 : if(max(chanMap)==-1)
915 : return;
916 :
917 : // Get the uvws in a form that Fortran can use
918 : Matrix<Double> uvw(negateUV(vb));
919 : Vector<Double> dphase(vb.nRows());
920 : dphase=0.0;
921 :
922 : doUVWRotation_p=true;
923 : girarUVW(uvw, dphase, vb);
924 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
925 :
926 :
927 :
928 :
929 : Cube<Complex> data;
930 : Cube<Int> flags;
931 : getInterpolateArrays(vb, data, flags);
932 :
933 : //Need to get interpolated freqs
934 : findConvFunction(*image, vb);
935 :
936 : // no valid pointing in this buffer
937 : if(convSupport <= 0)
938 : return;
939 : Complex *datStorage;
940 : Bool isCopy;
941 : datStorage=data.getStorage(isCopy);
942 :
943 :
944 : Vector<Int> rowFlags(vb.nRows());
945 : rowFlags=0;
946 : rowFlags(vb.flagRow())=true;
947 : if(!usezero_p) {
948 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
949 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
950 : }
951 : }
952 : Int nvp=data.shape()[0];
953 : Int nvc=data.shape()[1];
954 : Int nvisrow=data.shape()[2];
955 : Int csamp=convSampling;
956 : Int csize=convSize;
957 : //Int csupp=convSupport;
958 : Int nc=nchan;
959 : Int np=npol;
960 : Int nxp=nx;
961 : Int nyp=ny;
962 : Bool uvwcopy;
963 : const Double *uvwstor=uvw.getStorage(uvwcopy);
964 : Int nPolConv=convFunc.shape()[2];
965 : Int nChanConv=convFunc.shape()[3];
966 : Int nConvFunc=convFunc.shape()(4);
967 : ////////@@@@@@@@@
968 : Cube<Int> loc(2, nvc, nRow);
969 : Cube<Int> off(2, nvc, nRow);
970 : Matrix<Complex> phasor(nvc, nRow);
971 : Bool delphase;
972 : Bool del;
973 : const Int* pmapstor=polMap.getStorage(del);
974 : const Int* cmapstor=chanMap.getStorage(del);
975 : Complex * phasorstor=phasor.getStorage(delphase);
976 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
977 : const Double * scalestor=uvScale.getStorage(del);
978 : const Double * offsetstor=uvOffset.getStorage(del);
979 : const Int * flagstor=flags.getStorage(del);
980 : const Int * rowflagstor=rowFlags.getStorage(del);
981 : Int * locstor=loc.getStorage(del);
982 : Int * offstor=off.getStorage(del);
983 : const Double *dpstor=dphase.getStorage(del);
984 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
985 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
986 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
987 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
988 : ////////@@@@@@@@@@@@@@@@@@@@@
989 :
990 : Int irow;
991 : Int nth=1;
992 : #ifdef _OPENMP
993 : if(numthreads_p >0){
994 : nth=min(numthreads_p, omp_get_max_threads());
995 : }
996 : else{
997 : nth= omp_get_max_threads();
998 : }
999 : //nth=min(4,nth);
1000 : #endif
1001 :
1002 : Timer tim;
1003 : tim.mark();
1004 :
1005 : Int dow=0;
1006 : Double cinv=Double(1.0)/C::c;
1007 : #pragma omp parallel default(none) private(irow) firstprivate(visfreqstor, nvc, scalestor, offsetstor, csamp, phasorstor, uvwstor, locstor, offstor, dpstor, dow, cinv) shared(startRow, endRow) num_threads(nth)
1008 : {
1009 : #pragma omp for
1010 : for (irow=startRow; irow<=endRow;irow++){
1011 : /////////////////locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
1012 : // locstor,
1013 : /////////// offstor, phasorstor, irow, false);
1014 : //using the fortran version which is significantly faster ...this can account for 10% less overall degridding time
1015 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor,
1016 : offstor, phasorstor, &irow, &dow, &cinv);
1017 : }
1018 :
1019 : }//end pragma parallel
1020 : Int rbeg=startRow+1;
1021 : Int rend=endRow+1;
1022 : Int npart=nth;
1023 :
1024 : Bool gridcopy;
1025 : const Complex *gridstor=griddedData.getStorage(gridcopy);
1026 : Bool convcopy;
1027 : ////Degridding needs the conjugate ...doing it here
1028 : Array<Complex> conjConvFunc=conj(convFunc);
1029 : const Complex *convstor=conjConvFunc.getStorage(convcopy);
1030 : Int ix=0;
1031 : #pragma omp parallel default(none) private(ix, rbeg, rend) firstprivate(uvwstor, datStorage, flagstor, rowflagstor, convstor, pmapstor, cmapstor, gridstor, nxp, nyp, np, nc, csamp, csize, convsupportstor, nvp, nvc, nvisrow, phasorstor, locstor, offstor, nPolConv, nChanConv, nConvFunc, convrowmapstor, convpolmapstor, convchanmapstor, npart) num_threads(npart)
1032 : {
1033 : #pragma omp for schedule(dynamic)
1034 : for (ix=0; ix< npart; ++ix){
1035 : rbeg=ix*(nvisrow/npart)+1;
1036 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)-1+nvisrow%npart) ;
1037 : //cerr << "maps " << convChanMap_p << " " << chanMap << endl;
1038 : //cerr << "nchan " << nchan << " nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
1039 : sectdmos3(
1040 : datStorage,
1041 : &nvp,
1042 : &nvc,
1043 : flagstor,
1044 : rowflagstor,
1045 : &nvisrow,
1046 : gridstor,
1047 : &nxp,
1048 : &nyp,
1049 : &np,
1050 : &nc,
1051 : convsupportstor,
1052 : &csize,
1053 : &csamp,
1054 : convstor,
1055 : cmapstor,
1056 : pmapstor,
1057 : convrowmapstor, convchanmapstor,
1058 : convpolmapstor,
1059 : &nConvFunc, &nChanConv, &nPolConv,
1060 : &rbeg, &rend, locstor, offstor, phasorstor
1061 : );
1062 :
1063 :
1064 : }
1065 : }//end pragma omp
1066 :
1067 :
1068 : data.putStorage(datStorage, isCopy);
1069 : griddedData.freeStorage(gridstor, gridcopy);
1070 : convFunc.freeStorage(convstor, convcopy);
1071 :
1072 : timedegrid_p+=tim.real();
1073 :
1074 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
1075 : }
1076 : */
1077 :
1078 : } // REFIM ends
1079 : } //# NAMESPACE CASA - END
|
