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 112 : 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 112 : useDoubleGrid_p=True; //We'll always use double prec for this grid
53 :
54 112 : }
55 3 : AWPLPG::AWPLPG(const AWPLPG& other) : MosaicFTNew(other)
56 : {
57 3 : operator=(other);
58 3 : }
59 :
60 3 : AWPLPG& AWPLPG::operator=(const AWPLPG& other) {
61 3 : if(this!=&other) {
62 :
63 : //Do the base parameters
64 3 : MosaicFTNew::operator=(other);
65 3 : awConvs_p = other.awConvs_p;
66 3 : doSquint_p = other.doSquint_p;
67 3 : paInc_p = other.paInc_p;
68 3 : nw_p = other.nw_p;
69 : }
70 3 : return *this;
71 :
72 : }
73 3 : refim::FTMachine* AWPLPG::cloneFTM(){
74 3 : return new AWPLPG(*this);
75 : }
76 95 : void AWPLPG::init(const vi::VisBuffer2& vb){
77 95 : MosaicFTNew::init(vb);
78 :
79 :
80 : //oversample if image is small
81 : //But not more than 10000 pixels
82 95 : convSampling=(max(nx, ny) < 100) ? 100: Int(ceil(10000.0/max(nx, ny)));
83 :
84 95 : 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 95 : CoordinateSystem cs=image->coordinates();
94 :
95 95 : 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 95 : Int spw = vb.spectralWindows()(0);
101 95 : MDirection d;
102 95 : cs.directionCoordinate(0).toWorld(d, Vector<Double>(2,0));
103 95 : MPosition p= cs.obsInfo().telescopePosition();
104 95 : MEpoch e = cs.obsInfo().obsDate();
105 95 : fframe=(MFrequency::Types)vb.subtableColumns().spectralWindow().measFreqRef()(spw);
106 95 : spCS.setReferenceConversion(fframe, e, p, d);
107 :
108 95 : }
109 :
110 95 : nchan = image->shape()(3);
111 95 : spCS.toWorld(f1, double(-0.5));
112 95 : spCS.toWorld(f2, double(nchan)-0.5);
113 95 : auto frange=std::make_pair(f1, f2);
114 :
115 95 : if (pbConvFunc_p.null())
116 0 : pbConvFunc_p=new HetArrayConvFunc();
117 95 : awConvs_p=pbConvFunc_p->getAWConvFuncHolder();
118 95 : if(awConvs_p.use_count()==0){
119 45 : String observatory=(vb.subtableColumns().observation()).telescopeName()(0);
120 90 : awConvs_p=std::make_shared<AWConvFuncHolder>((*image).coordinates(), nx, ny,
121 90 : doSquint_p, paInc_p, observatory, convSampling);
122 45 : vi::VisibilityIterator2 *vi= const_cast<VisibilityIterator2 *>(vb.getVi());
123 :
124 45 : if(sj_p)
125 0 : pbConvFunc_p->setSkyJones(sj_p.get());
126 :
127 45 : if((pbConvFunc_p->getVBUtil()).null()){
128 45 : if(vbutil_p.null()){
129 45 : vbutil_p=new VisBufferUtil(vb);
130 : }
131 45 : pbConvFunc_p->setVBUtil(vbutil_p);
132 : }
133 45 : std::vector<Double> freqs;
134 45 : std::set<Int> fields;
135 45 : std::vector<Double> pAs;
136 45 : Double maxW=0.0;
137 124 : for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
138 3647 : for (vi->origin(); vi->more(); vi->next()) {
139 3568 : std::vector<Double> chunkfreq;
140 :
141 3568 : fields.insert(vb.fieldId()(0));
142 3568 : pbConvFunc_p->findUsefulChannels(chunkfreq, vb, frange);
143 : //cerr << "chunkfreq " << chunkfreq << endl;
144 3568 : if (chunkfreq.size() > 0) {
145 : //cerr << "SPW " << vb.spectralWindows()(0) << " freqs " << Vector<Double>(chunkfreq) << endl;
146 3568 : std::move(chunkfreq.begin(), chunkfreq.end(), std::back_inserter(freqs));
147 3568 : double maxfreqused = *(std::max_element(chunkfreq.begin(), chunkfreq.end()));
148 3568 : if (nw_p > 1) {
149 : // maxW=max(maxW, max(abs(vb.uvw().row(2)*max(vb.getFrequencies(0))))/C::c);
150 160 : maxW = max(maxW, max(abs(vb.uvw().row(2) * maxfreqused)) / C::c);
151 : }
152 : }
153 3568 : if (doSquint_p)
154 140 : 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 3568 : }
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 45 : vi->originChunks(); vi->origin();
166 : //cerr << "FREQS " << Vector<Double>(freqs) << endl;
167 45 : std::sort(freqs.begin(), freqs.end());
168 45 : auto last = std::unique(freqs.begin(), freqs.end());
169 45 : freqs.erase(last, freqs.end());
170 :
171 45 : 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 45 : (*awConvs_p).setSingleField((fields.size()==1) && (nw_p==1));
181 :
182 :
183 45 : Double paMax=0.0;
184 45 : if(pAs.size()>1){
185 4 : std::sort(pAs.begin(), pAs.end());
186 4 : last=std::unique(pAs.begin(), pAs.end());
187 4 : pAs.erase(last, pAs.end());
188 4 : if(pAs.size()==1)
189 0 : paMax=pAs[0];
190 4 : if(pAs.size() >1){
191 4 : auto [minpa, maxpa]=std::minmax_element(begin(pAs), end(pAs));
192 4 : paMax=max(abs(*minpa), *maxpa);
193 : }
194 : }
195 :
196 :
197 45 : if (nw_p == 0)
198 0 : nw_p = 1;
199 45 : Vector<Double> wVals(nw_p,0);
200 45 : if(nw_p >1){
201 5 : Double st=maxW/(Double(nw_p-1)*Double(nw_p-1));
202 25 : for (int k=0; k <nw_p; ++k)
203 20 : wVals[k]=Double(k*k)*st;
204 : }
205 45 : (*awConvs_p).addConvFunc(Vector<Double>(freqs), wVals, paMax);
206 :
207 45 : pbConvFunc_p->setAWConvFuncHolder(awConvs_p);
208 :
209 45 : }
210 :
211 95 : }
212 :
213 11137 : 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 11137 : awConvs_p->getConvFuncs(convPolMap_p, convChanMap_p, convRowMap_p, convFunc,
262 11137 : weightConvFunc_p, vb, rotuvw, interpVisFreq_p, toVis_p, ispsf);
263 :
264 : //double time1=omp_get_wtime();
265 : //cerr << " assign time " << time1-time0 << endl;
266 11137 : convSizePlanes_p.resize();
267 11137 : convSizePlanes_p = awConvs_p->getConvSizes();
268 11137 : convSupportPlanes_p.resize();
269 11137 : 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 11137 : std::vector<Int> pmapused = convPolMap_p.tovector();
275 : {
276 11137 : std::sort(pmapused.begin(), pmapused.end());
277 11137 : auto last = std::unique(pmapused.begin(), pmapused.end());
278 11137 : pmapused.erase(last, pmapused.end());
279 : }
280 11137 : std::vector<Int> cmapused=convChanMap_p.tovector();
281 : {
282 11137 : std::sort(cmapused.begin(), cmapused.end());
283 11137 : auto last = std::unique(cmapused.begin(), cmapused.end());
284 11137 : cmapused.erase(last, cmapused.end());
285 : }
286 11137 : std::vector<Int> rmapused=abs(convRowMap_p).tovector();
287 : {
288 11137 : std::sort(rmapused.begin(), rmapused.end());
289 11137 : auto last = std::unique(rmapused.begin(), rmapused.end());
290 11137 : rmapused.erase(last, rmapused.end());
291 : }
292 :
293 11137 : pbConvFunc_p->rephaseConvFunc(
294 11137 : iimage, vb, convSampling, convFunc, weightConvFunc_p, pmapused,
295 22274 : cmapused, rmapused, MVDirection(-(movingDirShift_p.getAngle())),
296 11137 : fixMovingSource_p);
297 11137 : convSupport =max(convSupportPlanes_p);
298 11137 : convSize = max(convSizePlanes_p);
299 11137 : 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 11137 : }
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
|
