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 5000 pixels
82 0 : convSampling=(max(nx, ny) < 50) ? 100: Int(ceil(5000.0/max(nx, ny)));
83 0 : if(convSampling <4)
84 0 : convSampling=4;
85 : // TESTOO
86 : //convSampling = 1;
87 : // TESTOO
88 :
89 :
90 :
91 0 : if(awConvs_p.use_count()==0){
92 0 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
93 0 : awConvs_p=std::make_shared<AWConvFuncHolder>((*image).coordinates(), nx, ny,
94 0 : doSquint_p, paInc_p, observatory, convSampling);
95 0 : vi::VisibilityIterator2 *vi= const_cast<VisibilityIterator2 *>(vb.getVi());
96 :
97 0 : if(pbConvFunc_p.null())
98 0 : pbConvFunc_p=new HetArrayConvFunc();
99 0 : if(sj_p)
100 0 : pbConvFunc_p->setSkyJones(sj_p.get());
101 :
102 0 : if((pbConvFunc_p->getVBUtil()).null()){
103 0 : if(vbutil_p.null()){
104 0 : vbutil_p=new VisBufferUtil(vb);
105 : }
106 0 : pbConvFunc_p->setVBUtil(vbutil_p);
107 : }
108 0 : std::vector<Double> freqs;
109 0 : std::vector<Double> pAs;
110 0 : Double maxW=0.0;
111 0 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
112 0 : for (vi->origin(); vi->more(); vi->next()) {
113 0 : std::vector<Double> chunkfreq;
114 0 : pbConvFunc_p->findUsefulChannels(chunkfreq, vb);
115 : //cerr << "chunkfreq " << chunkfreq << endl;
116 0 : std::move(chunkfreq.begin(), chunkfreq.end(), std::back_inserter(freqs));
117 0 : if(doSquint_p)
118 0 : pAs.push_back(getPA(vb));
119 0 : if(nw_p > 1)
120 0 : maxW=max(maxW, max(abs(vb.uvw().row(2)*max(vb.getFrequencies(0))))/C::c);
121 0 : }
122 : }
123 :
124 : //return vi to origin
125 0 : vi->originChunks(); vi->origin();
126 :
127 0 : std::sort(freqs.begin(), freqs.end());
128 0 : auto last = std::unique(freqs.begin(), freqs.end());
129 0 : freqs.erase(last, freqs.end());
130 :
131 0 : Double paMax=0.0;
132 0 : if(pAs.size()>1){
133 0 : std::sort(pAs.begin(), pAs.end());
134 0 : last=std::unique(pAs.begin(), pAs.end());
135 0 : pAs.erase(last, pAs.end());
136 0 : if(pAs.size()==1)
137 0 : paMax=pAs[0];
138 0 : if(pAs.size() >1){
139 0 : auto [minpa, maxpa]=std::minmax_element(begin(pAs), end(pAs));
140 0 : paMax=max(abs(*minpa), *maxpa);
141 : }
142 : }
143 :
144 0 : cerr << "PAMax in data " << paMax << endl;
145 0 : if (nw_p == 0)
146 0 : nw_p = 1;
147 0 : Vector<Double> wVals(nw_p,0);
148 0 : if(nw_p >1){
149 0 : Double st=maxW/(Double(nw_p-1)*Double(nw_p-1));
150 0 : for (int k=0; k <nw_p; ++k)
151 0 : wVals[k]=Double(k*k)*st;
152 : }
153 0 : (*awConvs_p).addConvFunc(Vector<Double>(freqs), wVals, paMax);
154 :
155 0 : }
156 :
157 0 : }
158 :
159 0 : void AWPLPG::findConvFunction(const ImageInterface<Complex>& iimage, const vi::VisBuffer2& vb, const Matrix<Double>& rotuvw ){
160 : //
161 : // pbConvFunc_p.phasegradient
162 0 : convFunc.resize();
163 0 : convFunc.assign(awConvs_p->getConvFunc());
164 :
165 0 : weightConvFunc_p.resize();
166 0 : weightConvFunc_p.assign(awConvs_p->getWeightConvFunc());
167 : /*{
168 : ////TESTOO
169 : IPosition elshp = convFunc.shape().getFirst(4);
170 : IPosition elblc(5, 0);
171 :
172 : IPosition eltrc = convFunc.shape()-1;
173 : elblc[4] = eltrc[4];
174 : CoordinateSystem csysA = iimage.coordinates();
175 : Vector<Int> stoks(4);
176 : stoks(0) = Stokes::RR;
177 : stoks(1) = Stokes::RL;
178 : stoks(2) = Stokes::LR;
179 : stoks(3) = Stokes::LL;
180 : StokesCoordinate stokey(stoks);
181 : csysA.replaceCoordinate(stokey, 1);
182 : PagedImage<Complex> lastplane(elshp, csysA, "COOBOO");
183 : lastplane.put(convFunc(elblc, eltrc).nonDegenerate());
184 : PagedImage<Complex> lastplaneW(elshp, csysA, "WOOBOO");
185 : lastplaneW.put(weightConvFunc_p(elblc, eltrc).nonDegenerate());
186 : //////
187 : } */
188 0 : convSizePlanes_p.resize();
189 0 : convSizePlanes_p = awConvs_p->getConvSizes();
190 0 : convSupportPlanes_p.resize();
191 0 : convSupportPlanes_p = awConvs_p->getConvSupports();
192 0 : awConvs_p->getConvIndices(convPolMap_p, convChanMap_p, convRowMap_p, vb, rotuvw);
193 : //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;
194 0 : std::vector<Int> pmapused=convPolMap_p.tovector();
195 : {
196 0 : std::sort(pmapused.begin(), pmapused.end());
197 0 : auto last = std::unique(pmapused.begin(), pmapused.end());
198 0 : pmapused.erase(last, pmapused.end());
199 : }
200 0 : std::vector<Int> cmapused=convChanMap_p.tovector();
201 : {
202 0 : std::sort(cmapused.begin(), cmapused.end());
203 0 : auto last = std::unique(cmapused.begin(), cmapused.end());
204 0 : cmapused.erase(last, cmapused.end());
205 : }
206 0 : std::vector<Int> rmapused=abs(convRowMap_p).tovector();
207 : {
208 0 : std::sort(rmapused.begin(), rmapused.end());
209 0 : auto last = std::unique(rmapused.begin(), rmapused.end());
210 0 : rmapused.erase(last, rmapused.end());
211 : }
212 : //cerr << "pmap " << Vector<Int>(pmapused) << " cmp " << Vector<Int>(cmapused) << " rmap " << Vector<Int>(rmapused) << endl;
213 0 : pbConvFunc_p->rephaseConvFunc(iimage, vb, convSampling, convFunc, weightConvFunc_p, pmapused, cmapused, rmapused, MVDirection(-(movingDirShift_p.getAngle())), fixMovingSource_p);
214 0 : convSupport =max(convSupportPlanes_p);
215 0 : convSize = max(convSizePlanes_p);
216 :
217 :
218 0 : }
219 :
220 : /////==============================================
221 : //// some fortran defn
222 : #define NEED_UNDERSCORES
223 : #if defined(NEED_UNDERSCORES)
224 : #define sectgmosd3 sectgmosd3_
225 : #define sectdmos3 sectdmos3_
226 : #define gmoswgtd2 gmoswgtd2_
227 : #define locuvw locuvw_
228 : #endif
229 :
230 : extern "C" {
231 : void locuvw(const Double*, const Double*, const Double*, const Int*, const Double*, const Double*, const Int*,
232 : Int*, Int*, Complex*, const Int*, const Int*, const Double*);
233 : void gmoswgtd2(const Int*/*nvispol*/, const Int*/*nvischan*/,
234 : const Int*/*flag*/, const Int*/*rflag*/, const Float*/*weight*/, const Int*/*nrow*/,
235 : const Int*/*nx*/, const Int*/*ny*/, const Int*/*npol*/, const Int*/*nchan*/,
236 : const Int*/*support*/, const Int*/*convsize*/, const Int*/*sampling*/,
237 : const Int*/*chanmap*/, const Int*/*polmap*/,
238 : DComplex* /*weightgrid*/, Double* /*sumwt*/, const Complex*/*convweight*/, const Int*/*convplanemap*/,
239 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
240 : const Int*/*nconvplane*/, const Int*/*nconvchan*/, const Int*/*nconvpol*/, const Int*/*rbeg*/,
241 : const Int*/*rend*/, const Int*/*loc*/, const Int*/*off*/, const Complex*/*phasor*/);
242 :
243 :
244 : void sectgmosd3(const Complex* /*values*/,
245 : Int* /*nvispol*/, Int* /*nvischan*/,
246 : Int* /*dopsf*/, const Int* /*flag*/, const Int* /*rflag*/, const Float* /*weight*/,
247 : Int* /* nrow*/, DComplex* /*grid*/, Int* /*nx*/, Int* /*ny*/, Int * /*npol*/, Int * /*nchan */,
248 : const Int*/*support*/, Int*/*convsize*/, Int*/*sampling*/, const Complex*/*convfunc*/,
249 : const Int*/*chanmap*/, const Int*/*polmap*/,
250 : Double*/*sumwgt*/, const Int*/*convplanemap*/,
251 : const Int*/*convchanmap*/, const Int*/*convpolmap*/,
252 : Int*/*nconvplane*/, Int*/*nconvchan*/, Int* /*nconvpol*/,
253 : const Int*/*x0*/,const Int*/*y0*/, const Int*/*nxsub*/, const Int*/*nysub*/, const Int*/*rbeg*/,
254 : const Int* /*rend*/, const Int*/*loc*/, const Int* /*off*/, const Complex*/*phasor*/);
255 :
256 : void sectdmos3(Complex*,
257 : Int*,
258 : Int*,
259 : const Int*,
260 : const Int*,
261 : Int*,
262 : const Complex*,
263 : Int*,
264 : Int*,
265 : Int *,
266 : Int *,
267 : const Int*, //support
268 : Int*,
269 : Int*,
270 : const Complex*,
271 : const Int*,
272 : const Int*,
273 : const Int*,
274 : const Int*,
275 : const Int*,
276 : Int*, Int*, Int*,
277 : //rbeg
278 : const Int*,
279 : const Int*,
280 : const Int*,
281 : const Int*,
282 : const Complex*);
283 :
284 :
285 :
286 : }
287 :
288 :
289 :
290 :
291 :
292 :
293 :
294 :
295 :
296 :
297 :
298 : //===================================================
299 : /* void AWPLPG::put(const vi::VisBuffer2& vb, Int row, Bool dopsf,
300 : FTMachine::Type type)
301 : {
302 :
303 :
304 :
305 :
306 : Timer tim;
307 : tim.mark();
308 :
309 : matchChannel(vb);
310 :
311 :
312 : //cerr << "CHANMAP " << chanMap << endl;
313 : //No point in reading data if its not matching in frequency
314 : if(max(chanMap)==-1)
315 : return;
316 :
317 : //const Matrix<Float> *imagingweight;
318 : //imagingweight=&(vb.imagingWeight());
319 : Matrix<Float> imagingweight;
320 : getImagingWeight(imagingweight, vb);
321 :
322 : if(dopsf) type=FTMachine::PSF;
323 :
324 : Cube<Complex> data;
325 : //Fortran gridder need the flag as ints
326 : Cube<Int> flags;
327 : Matrix<Float> elWeight;
328 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, type);
329 :
330 :
331 :
332 : Bool iswgtCopy;
333 : const Float *wgtStorage;
334 : wgtStorage=elWeight.getStorage(iswgtCopy);
335 :
336 :
337 :
338 :
339 : Bool isCopy;
340 : const Complex *datStorage=0;
341 :
342 : // cerr << "dopsf " << dopsf << " isWeightCopy " << iswgtCopy << " " << wgtStorage<< endl;
343 : if(!dopsf)
344 : datStorage=data.getStorage(isCopy);
345 :
346 :
347 : // If row is -1 then we pass through all rows
348 : Int startRow, endRow, nRow;
349 : if (row==-1) {
350 : nRow=vb.nRows();
351 : startRow=0;
352 : endRow=nRow-1;
353 : } else {
354 : nRow=1;
355 : startRow=row;
356 : endRow=row;
357 : }
358 :
359 : // Get the uvws in a form that Fortran can use and do that
360 : // necessary phase rotation.
361 : Matrix<Double> uvw(negateUV(vb));
362 : Vector<Double> dphase(vb.nRows());
363 : dphase=0.0;
364 :
365 : doUVWRotation_p=true;
366 : girarUVW(uvw, dphase, vb);
367 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
368 : // This needs to be after the interp to get the interpolated channels
369 : //Also has to be after rotateuvw in case tracking is on
370 : findConvFunction(*image, vb);
371 : //cerr << "Put convsup " << convSupport << " max min convFunc " << max(convFunc) << " " << min(convFunc) << " " << max(weightConvFunc_p) << min(weightConvFunc_p) << "SHP " << convFunc.shape() << " " << weightConvFunc_p.shape() << endl;
372 : //cerr << "convRowMap " << convRowMap_p << " " << convChanMap_p << " " << convPolMap_p << endl;
373 : //nothing to grid here as the pointing resulted in a zero support convfunc
374 : if(convSupport <= 0)
375 : return;
376 :
377 : // Get the pointing positions. This can easily consume a lot
378 : // of time thus we are for now assuming a field per
379 : // vb chunk...need to change that accordingly if we start using
380 : // multiple pointings per vb.
381 : //Warning
382 :
383 : // Take care of translation of Bools to Integer
384 : Int idopsf=0;
385 : if(dopsf) idopsf=1;
386 :
387 :
388 : Vector<Int> rowFlags(vb.nRows());
389 : rowFlags=0;
390 : rowFlags(vb.flagRow())=true;
391 : if(!usezero_p) {
392 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
393 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
394 : }
395 : }
396 :
397 :
398 :
399 : //cerr << "convSamp " << convSampling << " convsupp " << convSupport << " consize " << convSize << " convFunc " << convFunc.shape() << endl;
400 : //TESTOO
401 :
402 : //TESTOO
403 :
404 : //Tell the gridder to grid the weights too ...need to do that once only
405 : //Int doWeightGridding=1;
406 : //if(doneWeightImage_p)
407 : // doWeightGridding=-1;
408 : Bool del;
409 : // IPosition s(flags.shape());
410 : const IPosition& fs=flags.shape();
411 : //cerr << "flags shape " << fs << endl;
412 : std::vector<Int>s(fs.begin(), fs.end());
413 : Int nvp=s[0];
414 : Int nvc=s[1];
415 : Int nvisrow=s[2];
416 : Int csamp=convSampling;
417 : Bool uvwcopy;
418 : const Double *uvwstor=uvw.getStorage(uvwcopy);
419 : Bool gridcopy;
420 : Bool convcopy;
421 : Bool wconvcopy;
422 : const Complex *convstor=convFunc.getStorage(convcopy);
423 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
424 : Int nPolConv=convFunc.shape()[2];
425 : Int nChanConv=convFunc.shape()[3];
426 : Int nConvFunc=convFunc.shape()(4);
427 : Bool weightcopy;
428 : ////////
429 : Cube<Int> loc(2, nvc, nRow);
430 : Cube<Int> off(2, nvc, nRow);
431 : Matrix<Complex> phasor(nvc, nRow);
432 : Bool delphase;
433 : Complex * phasorstor=phasor.getStorage(delphase);
434 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
435 : const Double * scalestor=uvScale.getStorage(del);
436 : const Double * offsetstor=uvOffset.getStorage(del);
437 : Int * locstor=loc.getStorage(del);
438 : Int * offstor=off.getStorage(del);
439 : const Double *dpstor=dphase.getStorage(del);
440 : Int irow;
441 : Int nth=1;
442 : #ifdef _OPENMP
443 : if(numthreads_p >0){
444 : nth=min(numthreads_p, omp_get_max_threads());
445 : }
446 : else{
447 : nth= omp_get_max_threads();
448 : }
449 : //nth=min(4,nth);
450 : #endif
451 : Double cinv=Double(1.0)/C::c;
452 :
453 : Int dow=0;
454 : #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)
455 : {
456 : #pragma omp for
457 : for (irow=startRow; irow<=endRow;irow++){
458 :
459 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
460 : }
461 :
462 : }//end pragma parallel
463 :
464 :
465 :
466 : timemass_p +=tim.real();
467 : Int ixsub, iysub, icounter;
468 : ixsub=1;
469 : iysub=1;
470 : //////@@@@@@@@@@@@@DEBUGGING
471 : //nth=1;
472 : ////////@@@@@@@@@@@@@
473 : if (nth >3){
474 : ixsub=8;
475 : iysub=8;
476 : }
477 : else if(nth >1){
478 : ixsub=2;
479 : iysub=2;
480 : }
481 : Int rbeg=startRow+1;
482 : Int rend=endRow+1;
483 : Block<Matrix<Double> > sumwgt(ixsub*iysub);
484 : Vector<Double *> swgtptr(ixsub*iysub);
485 : Vector<Bool> swgtdel(ixsub*iysub);
486 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
487 : sumwgt[icounter].resize(sumWeight.shape());
488 : sumwgt[icounter].set(0.0);
489 : swgtptr[icounter]=sumwgt[icounter].getStorage(swgtdel(icounter));
490 : }
491 : //cerr << "done thread " << doneThreadPartition_p << " " << ixsub*iysub << endl;
492 : if(doneThreadPartition_p < 0){
493 : xsect_p.resize(ixsub*iysub);
494 : ysect_p.resize(ixsub*iysub);
495 : nxsect_p.resize(ixsub*iysub);
496 : nysect_p.resize(ixsub*iysub);
497 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
498 : findGridSector(nx, ny, ixsub, iysub, 0, 0, icounter, xsect_p(icounter), ysect_p(icounter), nxsect_p(icounter), nysect_p(icounter), true);
499 : }
500 : }
501 : Vector<Int> xsect, ysect, nxsect, nysect;
502 : xsect=xsect_p; ysect=ysect_p; nxsect=nxsect_p; nysect=nysect_p;
503 : //cerr << xsect.shape() << " " << xsect << endl;
504 : const Int* pmapstor=polMap.getStorage(del);
505 : const Int* cmapstor=chanMap.getStorage(del);
506 : // Dummy sumwt for gridweight part
507 : Matrix<Double> dumSumWeight(npol, nchan);
508 : dumSumWeight=sumWeight;
509 : Bool isDSWC;
510 : Double *dsumwtstor=dumSumWeight.getStorage(isDSWC);
511 : Int nc=nchan;
512 : Int np=npol;
513 : Int nxp=nx;
514 : Int nyp=ny;
515 : Int csize=convSize;
516 : const Int * flagstor=flags.getStorage(del);
517 : const Int * rowflagstor=rowFlags.getStorage(del);
518 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
519 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
520 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
521 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
522 : ///
523 :
524 :
525 : ////////
526 : tim.mark();
527 :
528 : // if(useDoubleGrid_p) { //always using double prec here
529 : {
530 : DComplex *gridstor=griddedData2.getStorage(gridcopy);
531 :
532 : #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)
533 : {
534 : #pragma omp for schedule(dynamic)
535 : for(icounter=0; icounter < ixsub*iysub; ++icounter){
536 : Int x0=xsect(icounter);
537 : Int y0=ysect(icounter);
538 : Int nxsub=nxsect(icounter);
539 : Int nysub=nysect(icounter);
540 :
541 :
542 : sectgmosd3(datStorage,
543 : &nvp,
544 : &nvc,
545 : &idopsf,
546 : flagstor,
547 : rowflagstor,
548 : wgtStorage,
549 : &nvisrow,
550 : gridstor,
551 : &nxp,
552 : &nyp,
553 : &np,
554 : &nc,
555 : convsupportstor,
556 : &csize,
557 : &csamp,
558 : convstor,
559 : cmapstor,
560 : pmapstor,
561 : swgtptr[icounter],
562 : convrowmapstor,
563 : convchanmapstor,
564 : convpolmapstor,
565 : &nConvFunc, &nChanConv, &nPolConv,
566 : &x0, &y0, &nxsub, &nysub, &rbeg, &rend, locstor, offstor,
567 : phasorstor
568 : );
569 : }
570 : }//end pragma parallel
571 : for (icounter=0; icounter < ixsub*iysub; ++icounter){
572 : sumwgt[icounter].putStorage(swgtptr[icounter],swgtdel[icounter]);
573 : sumWeight=sumWeight+sumwgt[icounter];
574 : }
575 :
576 : //cerr << "SUMWEIG " << sumWeight << endl;
577 : griddedData2.putStorage(gridstor, gridcopy);
578 : if(dopsf && (nth >4))
579 : tweakGridSector(nx, ny, ixsub, iysub);
580 : timegrid_p+=tim.real();
581 : tim.mark();
582 : if(!doneWeightImage_p){
583 : //This can be parallelized by making copy of the central part of the griddedWeight
584 : //and adding it after dooing the gridding
585 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
586 : gmoswgtd2(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
587 : &nxp, &nyp, &np, &nc, convsupportstor, &csize, &csamp,
588 : cmapstor, pmapstor,
589 : gridwgtstor, dsumwtstor, wconvstor, convrowmapstor,
590 : convchanmapstor, convpolmapstor,
591 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
592 : &rend, locstor, offstor, phasorstor);
593 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
594 :
595 : }
596 : timemass_p+=tim.real();
597 : }
598 :
599 : convFunc.freeStorage(convstor, convcopy);
600 : weightConvFunc_p.freeStorage(wconvstor, wconvcopy);
601 : dumSumWeight.putStorage(dsumwtstor, isDSWC);
602 : //cerr << "dumSumwe " << dumSumWeight << endl;
603 : uvw.freeStorage(uvwstor, uvwcopy);
604 : if(!dopsf)
605 : data.freeStorage(datStorage, isCopy);
606 :
607 : elWeight.freeStorage(wgtStorage,iswgtCopy);
608 :
609 :
610 :
611 :
612 : }
613 : */
614 :
615 :
616 : /*
617 : void AWPLPG::gridImgWeights(const vi::VisBuffer2& vb){
618 :
619 : if(doneWeightImage_p)
620 : return;
621 : matchChannel(vb);
622 :
623 :
624 : //cerr << "CHANMAP " << chanMap << endl;
625 : //No point in reading data if its not matching in frequency
626 : if(max(chanMap)==-1)
627 : return;
628 :
629 : Int startRow, endRow, nRow;
630 : nRow=vb.nRows();
631 : startRow=0;
632 : endRow=nRow-1;
633 :
634 :
635 : //const Matrix<Float> *imagingweight;
636 : //imagingweight=&(vb.imagingWeight());
637 : Matrix<Float> imagingweight;
638 : getImagingWeight(imagingweight, vb);
639 :
640 :
641 : Cube<Complex> data;
642 : //Fortran gridder need the flag as ints
643 : Cube<Int> flags;
644 : Matrix<Float> elWeight;
645 : interpolateFrequencyTogrid(vb, imagingweight,data, flags, elWeight, FTMachine::PSF);
646 :
647 :
648 :
649 : Bool iswgtCopy;
650 : const Float *wgtStorage;
651 : wgtStorage=elWeight.getStorage(iswgtCopy);
652 : Bool issumWgtCopy;
653 : Double* sumwgtstor=sumWeight.getStorage(issumWgtCopy);
654 :
655 :
656 :
657 : // Get the uvws in a form that Fortran can use and do that
658 : // necessary phase rotation.
659 : Matrix<Double> uvw(negateUV(vb));
660 : Vector<Double> dphase(vb.nRows());
661 : dphase=0.0;
662 :
663 : doUVWRotation_p=true;
664 : girarUVW(uvw, dphase, vb);
665 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
666 : // This needs to be after the interp to get the interpolated channels
667 : //Also has to be after rotateuvw in case tracking is on
668 : findConvFunction(*image, vb);
669 : //nothing to grid here as the pointing resulted in a zero support convfunc
670 : if(convSupport <= 0)
671 : return;
672 :
673 : Bool del;
674 :
675 : const Int* pmapstor=polMap.getStorage(del);
676 : const Int* cmapstor=chanMap.getStorage(del);
677 :
678 : Vector<Int> rowFlags(vb.nRows());
679 : rowFlags=0;
680 : rowFlags(vb.flagRow())=true;
681 : if(!usezero_p) {
682 : for (uInt rownr=0; rownr< vb.nRows(); rownr++) {
683 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
684 : }
685 : }
686 :
687 : //Fortran indexing
688 :
689 : Int rbeg=1;
690 : Int rend=vb.nRows();
691 :
692 : const Int * flagstor=flags.getStorage(del);
693 : const Int * rowflagstor=rowFlags.getStorage(del);
694 :
695 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
696 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
697 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
698 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
699 : //Tell the gridder to grid the weights too ...need to do that once only
700 : //Int doWeightGridding=1;
701 : //if(doneWeightImage_p)
702 : // doWeightGridding=-1;
703 : // IPosition s(flags.shape());
704 : const IPosition& fs=flags.shape();
705 : //cerr << "flags shape " << fs << endl;
706 : std::vector<Int>s(fs.begin(), fs.end());
707 : Int nvp=s[0];
708 : Int nvc=s[1];
709 : Int nvisrow=s[2];
710 : Int csamp=convSampling;
711 : Bool uvwcopy;
712 : const Double *uvwstor=uvw.getStorage(uvwcopy);
713 : Bool gridcopy;
714 : Bool convcopy;
715 : Bool wconvcopy;
716 : const Complex *wconvstor=weightConvFunc_p.getStorage(wconvcopy);
717 : Int nPolConv=convFunc.shape()[2];
718 : Int nChanConv=convFunc.shape()[3];
719 : Int nConvFunc=convFunc.shape()(4);
720 : Bool weightcopy;
721 : ////////@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
722 : Cube<Int> loc(2, nvc, vb.nRows());
723 : Cube<Int> off(2, nvc, vb.nRows());
724 : Matrix<Complex> phasor(nvc, vb.nRows());
725 : Bool delphase;
726 : Complex * phasorstor=phasor.getStorage(delphase);
727 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
728 : const Double * scalestor=uvScale.getStorage(del);
729 : const Double * offsetstor=uvOffset.getStorage(del);
730 : Int * locstor=loc.getStorage(del);
731 : Int * offstor=off.getStorage(del);
732 : const Double *dpstor=dphase.getStorage(del);
733 :
734 : Int irow;
735 : Int nth=1;
736 : #ifdef _OPENMP
737 : if(numthreads_p >0){
738 : nth=min(numthreads_p, omp_get_max_threads());
739 : }
740 : else{
741 : nth= omp_get_max_threads();
742 : }
743 : //nth=min(4,nth);
744 : #endif
745 :
746 : Double cinv=Double(1.0)/C::c;
747 :
748 : Int dow=0;
749 :
750 : #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)
751 : {
752 : #pragma omp for
753 : for (irow=startRow; irow<=endRow;irow++){
754 :
755 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor, offstor, phasorstor, &irow, &dow, &cinv);
756 : }
757 :
758 : }//end pragma parallel
759 :
760 :
761 :
762 : //always using double prec in this gridder
763 : // if(useDoubleGrid_p) {
764 : {
765 : //This can be parallelized by making copy of the central part of the griddedWeight
766 : //and adding it after dooing the gridding
767 : DComplex *gridwgtstor=griddedWeight2.getStorage(weightcopy);
768 : gmoswgtd2(&nvp, &nvc,flagstor, rowflagstor, wgtStorage, &nvisrow,
769 : &nx, &ny, &npol, &nchan, convsupportstor, &convSize, &convSampling,
770 : cmapstor, pmapstor,
771 : gridwgtstor, sumwgtstor, wconvstor, convrowmapstor,
772 : convchanmapstor, convpolmapstor,
773 : &nConvFunc, &nChanConv, &nPolConv, &rbeg,
774 : &rend, locstor, offstor, phasorstor);
775 : griddedWeight2.putStorage(gridwgtstor, weightcopy);
776 :
777 :
778 :
779 :
780 :
781 : }
782 :
783 : sumWeight.putStorage(sumwgtstor, issumWgtCopy);
784 : elWeight.freeStorage(wgtStorage,iswgtCopy);
785 :
786 : }
787 : */
788 : /*
789 : void AWPLPG::get(vi::VisBuffer2& vb, Int row)
790 : {
791 :
792 :
793 :
794 : // If row is -1 then we pass through all rows
795 : Int startRow, endRow, nRow;
796 : if (row==-1) {
797 : nRow=vb.nRows();
798 : startRow=0;
799 : endRow=nRow-1;
800 : // vb.modelVisCube()=Complex(0.0,0.0);
801 : } else {
802 : nRow=1;
803 : startRow=row;
804 : endRow=row;
805 : // vb.modelVisCube().xyPlane(row)=Complex(0.0,0.0);
806 : }
807 :
808 :
809 :
810 :
811 : matchChannel(vb);
812 :
813 : //No point in reading data if its not matching in frequency
814 : if(max(chanMap)==-1)
815 : return;
816 :
817 : // Get the uvws in a form that Fortran can use
818 : Matrix<Double> uvw(negateUV(vb));
819 : Vector<Double> dphase(vb.nRows());
820 : dphase=0.0;
821 :
822 : doUVWRotation_p=true;
823 : girarUVW(uvw, dphase, vb);
824 : refocus(uvw, vb.antenna1(), vb.antenna2(), dphase, vb);
825 :
826 :
827 :
828 :
829 : Cube<Complex> data;
830 : Cube<Int> flags;
831 : getInterpolateArrays(vb, data, flags);
832 :
833 : //Need to get interpolated freqs
834 : findConvFunction(*image, vb);
835 :
836 : // no valid pointing in this buffer
837 : if(convSupport <= 0)
838 : return;
839 : Complex *datStorage;
840 : Bool isCopy;
841 : datStorage=data.getStorage(isCopy);
842 :
843 :
844 : Vector<Int> rowFlags(vb.nRows());
845 : rowFlags=0;
846 : rowFlags(vb.flagRow())=true;
847 : if(!usezero_p) {
848 : for (Int rownr=startRow; rownr<=endRow; rownr++) {
849 : if(vb.antenna1()(rownr)==vb.antenna2()(rownr)) rowFlags(rownr)=1;
850 : }
851 : }
852 : Int nvp=data.shape()[0];
853 : Int nvc=data.shape()[1];
854 : Int nvisrow=data.shape()[2];
855 : Int csamp=convSampling;
856 : Int csize=convSize;
857 : //Int csupp=convSupport;
858 : Int nc=nchan;
859 : Int np=npol;
860 : Int nxp=nx;
861 : Int nyp=ny;
862 : Bool uvwcopy;
863 : const Double *uvwstor=uvw.getStorage(uvwcopy);
864 : Int nPolConv=convFunc.shape()[2];
865 : Int nChanConv=convFunc.shape()[3];
866 : Int nConvFunc=convFunc.shape()(4);
867 : ////////@@@@@@@@@
868 : Cube<Int> loc(2, nvc, nRow);
869 : Cube<Int> off(2, nvc, nRow);
870 : Matrix<Complex> phasor(nvc, nRow);
871 : Bool delphase;
872 : Bool del;
873 : const Int* pmapstor=polMap.getStorage(del);
874 : const Int* cmapstor=chanMap.getStorage(del);
875 : Complex * phasorstor=phasor.getStorage(delphase);
876 : const Double * visfreqstor=interpVisFreq_p.getStorage(del);
877 : const Double * scalestor=uvScale.getStorage(del);
878 : const Double * offsetstor=uvOffset.getStorage(del);
879 : const Int * flagstor=flags.getStorage(del);
880 : const Int * rowflagstor=rowFlags.getStorage(del);
881 : Int * locstor=loc.getStorage(del);
882 : Int * offstor=off.getStorage(del);
883 : const Double *dpstor=dphase.getStorage(del);
884 : const Int *convrowmapstor=convRowMap_p.getStorage(del);
885 : const Int *convchanmapstor=convChanMap_p.getStorage(del);
886 : const Int *convpolmapstor=convPolMap_p.getStorage(del);
887 : const Int *convsupportstor=convSupportPlanes_p.getStorage(del);
888 : ////////@@@@@@@@@@@@@@@@@@@@@
889 :
890 : Int irow;
891 : Int nth=1;
892 : #ifdef _OPENMP
893 : if(numthreads_p >0){
894 : nth=min(numthreads_p, omp_get_max_threads());
895 : }
896 : else{
897 : nth= omp_get_max_threads();
898 : }
899 : //nth=min(4,nth);
900 : #endif
901 :
902 : Timer tim;
903 : tim.mark();
904 :
905 : Int dow=0;
906 : Double cinv=Double(1.0)/C::c;
907 : #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)
908 : {
909 : #pragma omp for
910 : for (irow=startRow; irow<=endRow;irow++){
911 : /////////////////locateuvw(uvwstor,dpstor, visfreqstor, nvc, scalestor, offsetstor, csamp,
912 : // locstor,
913 : /////////// offstor, phasorstor, irow, false);
914 : //using the fortran version which is significantly faster ...this can account for 10% less overall degridding time
915 : locuvw(uvwstor, dpstor, visfreqstor, &nvc, scalestor, offsetstor, &csamp, locstor,
916 : offstor, phasorstor, &irow, &dow, &cinv);
917 : }
918 :
919 : }//end pragma parallel
920 : Int rbeg=startRow+1;
921 : Int rend=endRow+1;
922 : Int npart=nth;
923 :
924 : Bool gridcopy;
925 : const Complex *gridstor=griddedData.getStorage(gridcopy);
926 : Bool convcopy;
927 : ////Degridding needs the conjugate ...doing it here
928 : Array<Complex> conjConvFunc=conj(convFunc);
929 : const Complex *convstor=conjConvFunc.getStorage(convcopy);
930 : Int ix=0;
931 : #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)
932 : {
933 : #pragma omp for schedule(dynamic)
934 : for (ix=0; ix< npart; ++ix){
935 : rbeg=ix*(nvisrow/npart)+1;
936 : rend=(ix != (npart-1)) ? (rbeg+(nvisrow/npart)-1) : (rbeg+(nvisrow/npart)-1+nvisrow%npart) ;
937 : //cerr << "maps " << convChanMap_p << " " << chanMap << endl;
938 : //cerr << "nchan " << nchan << " nchanconv " << nChanConv << " npolconv " << nPolConv << " nRowConv " << nConvFunc << endl;
939 : sectdmos3(
940 : datStorage,
941 : &nvp,
942 : &nvc,
943 : flagstor,
944 : rowflagstor,
945 : &nvisrow,
946 : gridstor,
947 : &nxp,
948 : &nyp,
949 : &np,
950 : &nc,
951 : convsupportstor,
952 : &csize,
953 : &csamp,
954 : convstor,
955 : cmapstor,
956 : pmapstor,
957 : convrowmapstor, convchanmapstor,
958 : convpolmapstor,
959 : &nConvFunc, &nChanConv, &nPolConv,
960 : &rbeg, &rend, locstor, offstor, phasorstor
961 : );
962 :
963 :
964 : }
965 : }//end pragma omp
966 :
967 :
968 : data.putStorage(datStorage, isCopy);
969 : griddedData.freeStorage(gridstor, gridcopy);
970 : convFunc.freeStorage(convstor, convcopy);
971 :
972 : timedegrid_p+=tim.real();
973 :
974 : interpolateFrequencyFromgrid(vb, data, FTMachine::MODEL);
975 : }
976 : */
977 :
978 : } // REFIM ends
979 : } //# NAMESPACE CASA - END
|
