Line data Source code
1 : //# AWConvFuncHolser.cc: Implementation for Holding class for AW convolution functions
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 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/Matrix.h>
32 : #include <casacore/casa/Arrays/Vector.h>
33 : #include <casacore/measures/Measures/MeasTable.h>
34 : #include <msvis/MSVis/VisBuffer2.h>
35 : #include <msvis/MSVis/VisibilityIterator2.h>
36 : #include <synthesis/TransformMachines2/AWConvFunc.h>
37 : #include <synthesis/TransformMachines2/EVLAAperture.h>
38 : #include <synthesis/TransformMachines2/AWConvFuncHolder.h>
39 :
40 : #include <casacore/coordinates/Coordinates/TabularCoordinate.h>
41 : #include <iomanip>
42 : #include <casacore/casa/OS/Timer.h>
43 : namespace casa { //# CASA namespace
44 : namespace refim { //# namespace refactor imaging
45 :
46 : using namespace casacore;
47 : using namespace casa;
48 : using namespace casa::refim;
49 : using namespace std;
50 :
51 0 : AWConvFuncHolder::AWConvFuncHolder(const CoordinateSystem& csys, const int nx, const int ny, const bool dosquint, const double paInc, const String& obs, const int& oversamp): painc_p(paInc), dosquint_p(dosquint), outcsys_p(csys), nx_p(nx), ny_p(ny), oversamp_p(oversamp), isSingleField_p(false) {
52 :
53 0 : convFunc_p.resize();
54 0 : wgtConvFunc_p.resize();
55 0 : polVals_p.resize();
56 0 : freqVals_p.resize();
57 0 : wVals_p.resize();
58 0 : paVals_p.resize();
59 0 : antpairVals_p.resize();
60 0 : rowAxisWVals_p.resize();
61 0 : rowAxisPAVals_p.resize();
62 0 : rowAxisAntennaPair_p.resize();
63 0 : convSizes_p.resize();
64 0 : convSupport_p.resize();
65 0 : convSizesHPG_p.resize();
66 0 : convSupportHPG_p.resize();
67 0 : aterm_p = std::make_shared<refim::EVLAAperture>();
68 0 : aterm_p->cacheVBInfo(obs, 25.0);
69 :
70 0 : }
71 0 : AWConvFuncHolder::AWConvFuncHolder(const AWConvFuncHolder& other) {
72 0 : operator = (other);
73 0 : }
74 0 : AWConvFuncHolder& AWConvFuncHolder::operator=(const AWConvFuncHolder& other) {
75 0 : if (this != &other) {
76 0 : nx_p = other.nx_p;
77 0 : ny_p = other.ny_p;
78 0 : oversamp_p = other.oversamp_p;
79 0 : outcsys_p = other.outcsys_p;
80 0 : calcNpix_p = other.calcNpix_p;
81 0 : calcCsys_p = other.calcCsys_p;
82 0 : dosquint_p = other.dosquint_p;
83 0 : painc_p = other.painc_p;
84 0 : convFunc_p.resize();
85 0 : convFunc_p = other.convFunc_p;
86 0 : wgtConvFunc_p.resize();
87 0 : wgtConvFunc_p = other.wgtConvFunc_p;
88 0 : polVals_p.resize();
89 0 : polVals_p = other.polVals_p;
90 0 : freqVals_p.resize();
91 0 : freqVals_p = other.freqVals_p;
92 0 : wVals_p.resize();
93 0 : wVals_p = other.wVals_p;
94 0 : paVals_p.resize();
95 0 : paVals_p = other.paVals_p;
96 0 : antpairVals_p.resize();
97 0 : antpairVals_p = other.antpairVals_p;
98 0 : rowAxisWVals_p.resize();
99 0 : rowAxisWVals_p = other.rowAxisWVals_p;
100 0 : rowAxisPAVals_p.resize();
101 0 : rowAxisPAVals_p = other.rowAxisPAVals_p;
102 0 : rowAxisAntennaPair_p.resize();
103 0 : rowAxisAntennaPair_p = other.rowAxisAntennaPair_p;
104 0 : convSizes_p.resize();
105 0 : convSizesHPG_p.resize();
106 0 : convSizes_p = other.convSizes_p;
107 0 : convSizesHPG_p = other.convSizesHPG_p;
108 :
109 0 : convSupport_p.resize();
110 0 : convSupport_p.resize();
111 0 : convSupportHPG_p.resize();
112 0 : convSupport_p = other.convSupport_p;
113 0 : convSupportHPG_p = other.convSupportHPG_p;
114 0 : aterm_p = other.aterm_p;
115 0 : isSingleField_p = other.isSingleField_p;
116 : }
117 0 : return *this;
118 : }
119 0 : bool AWConvFuncHolder::addConvFunc(const casacore::Vector<casacore::Double>& freqs, const casacore::Vector<Double>& wVals, const casacore::Double& paMax) {
120 :
121 0 : Vector<Double> freqsToCalc;
122 0 : if (freqVals_p.nelements() == 0) {
123 0 : freqVals_p = freqs;
124 0 : freqsToCalc = freqs;
125 :
126 : }
127 : else{
128 : //append new freqs freqVals_p;
129 : // assign missing freqs to freqsToCalc
130 0 : freqsToCalc =freqs;
131 : }
132 0 : if (wVals_p.nelements() == 0) {
133 : // make sure first wval is 0;
134 0 : wVals_p = wVals;
135 :
136 : }
137 0 : else if (wVals_p.nelements() != wVals.nelements())
138 0 : throw(AipsError("Cannot change W terms length right now"));
139 0 : if (!dosquint_p) {
140 0 : paVals_p.resize(1);
141 0 : paVals_p[0] = 0.0;
142 : }
143 : else{
144 0 : cerr << "paMax " << paMax << " painc " << painc_p << endl;
145 0 : Vector<Double> pavals(int(std::ceil(2 * paMax / painc_p)));
146 : // setting pavals from -paMax to paMax
147 0 : for (uint k = 0; k < pavals.nelements(); ++k)
148 0 : pavals[k] = double(k) * painc_p - paMax;
149 0 : if (paVals_p.nelements() == 0)
150 0 : paVals_p = pavals;
151 0 : else if ((paVals_p.nelements()) != pavals.nelements())
152 0 : throw(AipsError("Cannot change number of PA's in between"));
153 :
154 0 : }
155 0 : std::shared_ptr<refim::WPConvFunc>wptr;
156 0 : AWConvFunc a(aterm_p, wptr);
157 : //Array<Complex> aWConv;
158 : //Array<Complex> aWwtconv;
159 : //Matrix<Int> awSupport;
160 0 : calcCsys_p = outcsys_p;
161 0 : calcNpix_p = min(nx_p, ny_p);
162 : //cerr << "PAVALS " << paVals_p << " dosquint " << dosquint_p << endl;
163 : //cerr << "FREQS " << freqsToCalc << endl;
164 :
165 0 : for (uint k = 0; k < paVals_p.nelements(); ++k) {
166 0 : calcNpix_p = min(nx_p, ny_p);
167 0 : calcCsys_p = outcsys_p;
168 0 : Array<Complex> aWConv;
169 0 : Array<Complex> aWwtconv;
170 0 : Matrix<Int> awSupport;
171 0 : a.makeAWConvFunc(aWConv, aWwtconv, calcCsys_p, awSupport, calcNpix_p,
172 0 : freqsToCalc, wVals_p, dosquint_p, paVals_p[k],
173 0 : isSingleField_p);
174 : //cerr << "######MAX awsupp " << max(awSupport) << endl;
175 0 : int startrow = k * wVals_p.nelements();
176 0 : appendConvFuncs(aWConv, aWwtconv, awSupport, freqsToCalc, paVals_p[k], startrow);
177 : // Let's resize for all paVals as resizing is costly
178 0 : if (k == 0 && paVals_p.nelements() > 1) {
179 0 : IPosition shp = convFunc_p.shape();
180 0 : shp[4] = wVals.nelements() * paVals_p.nelements();
181 0 : Double memAmt = Double(shp.product()) * 16.0;
182 0 : Double memAvl = Double(HostInfo::memoryFree()) * 1024.0;
183 : //cerr << "Memory needed " << memAmt << " available " << memAvl << endl;
184 0 : if (memAmt > 0.8 * memAvl) {
185 0 : throw(AipsError(
186 0 : "Not enough memory to hold all AW with squint correction convolution functions "));
187 : }
188 0 : convFunc_p.resize(shp, true);
189 0 : wgtConvFunc_p.resize(shp, true);
190 0 : }
191 :
192 0 : }
193 :
194 0 : return true;;
195 0 : }
196 0 : void AWConvFuncHolder::appendConvFuncs(const Array<Complex>& awConv, const Array<Complex>& aWwtConv, const Matrix<Int>& awsupport, const Vector<Double>& newFreqs, const Double paVal, const int startrow) {
197 0 : Int nfreqs = awConv.shape()[3];
198 : // Make sure the polVals are in the stokes used in making convfun
199 0 : polVals_p.resize(4);
200 : // for now antpairVals are for all pair of possible antenna combination
201 0 : antpairVals_p.resize(1);
202 0 : antpairVals_p[0] = std::pair<int, int>(-1, -1);
203 :
204 0 : convertArray(polVals_p, calcCsys_p.stokesCoordinate(1).stokes());
205 0 : IPosition blc = freqVals_p.shape();
206 0 : IPosition trc(1, blc[0]+nfreqs-1);
207 0 : if (freqVals_p.nelements() != 0) {
208 0 : if (!allEQ(freqVals_p, newFreqs))
209 0 : throw(AipsError("Not implemented appending different freqs yet to conv"));
210 : }
211 : else{
212 0 : freqVals_p.resize(freqVals_p.nelements()+nfreqs, True);
213 0 : freqVals_p(blc, trc) = newFreqs;
214 : }
215 0 : Int indPA = -1;
216 0 : for (uint k = 0; k < paVals_p.nelements(); ++k) {
217 0 : if (fabs(paVal-paVals_p[k]) < 1e-4)
218 0 : indPA = k;
219 : }
220 0 : if (indPA <0) {
221 0 : paVals_p.resize(paVals_p.nelements()+1, True);
222 0 : indPA = paVals_p.nelements()-1;
223 0 : paVals_p[indPA] = paVal;
224 :
225 : }
226 0 : blc[0] = rowAxisPAVals_p.nelements();
227 0 : trc[0] = blc[0]+wVals_p.nelements()-1;
228 0 : rowAxisPAVals_p.resize(trc[0]+1, True);
229 0 : rowAxisPAVals_p(blc, trc).set(indPA);
230 0 : rowAxisWVals_p.resize(trc[0]+1, True);
231 0 : Vector<Int> waxis(wVals_p.nelements());
232 0 : indgen(waxis);
233 0 : rowAxisWVals_p(blc, trc) = waxis;
234 : //cerr << "rowAxisWVals " << rowAxisWVals_p << endl;
235 0 : rowAxisAntennaPair_p.resize(trc[0] + 1, True);
236 0 : rowAxisAntennaPair_p(blc, trc).set(0);
237 : /// Let us rescale convolution function to match nx, ny and incr of image that makes
238 : /// uvgrid
239 0 : Float factorX=fabs(calcCsys_p.increment()(0)/outcsys_p.increment()(0));
240 0 : Float factorY=fabs(calcCsys_p.increment()(1)/outcsys_p.increment()(1));
241 : //cerr << "####Factor " << factorX << " " << factorY << endl;
242 0 : factorX = Float(nx_p) *Float(oversamp_p)/Float(calcNpix_p)/factorX;
243 0 : factorY = Float(ny_p) *Float(oversamp_p)/Float(calcNpix_p)/factorY;
244 :
245 : //cerr << "factors " << factorX << " " << factorY << "nx, ny" << nx_p << " " << ny_p << " calcNpix " << calcNpix_p << " oversamp " << oversamp_p << endl;
246 0 : MathUtils m;
247 0 : Array<Complex> newAWConv;
248 0 : Array<Complex> newWtConv;
249 : // For small images or factor less than 1.0 use linear interpolation
250 : //cerr << "Shape before " << awConv.shape() << endl;
251 : // factor/oversamp is ratio of im fov to pb fov
252 0 : if ((factorX / oversamp_p) < 1.0 || (factorY / oversamp_p) < 1.0 ||
253 0 : nx_p < 200 || ny_p < 200) {
254 0 : newAWConv = m.resample(awConv, factorX, factorY);
255 0 : newWtConv = m.resample(aWwtConv, factorX, factorY);
256 :
257 :
258 : } else {
259 0 : newAWConv = m.resampleViaFFT(awConv, factorX, factorY);
260 0 : newWtConv = m.resampleViaFFT(aWwtConv, factorX, factorY);
261 : }
262 : //cerr << "Shape after " << newAWConv.shape() << endl;
263 : Float correcfac =
264 0 : float(awConv.shape()(0) * awConv.shape()(1) * oversamp_p * oversamp_p) /
265 0 : float(newAWConv.shape()(0) * newAWConv.shape()(1));
266 :
267 0 : newAWConv *= correcfac;
268 0 : newWtConv *= correcfac;
269 :
270 : /*{
271 : ////TESTOO
272 : IPosition elshp = newAWConv.shape().getFirst(4);
273 : IPosition elblc(5, 0);
274 :
275 : IPosition eltrc = newAWConv.shape()-1;
276 : elblc[4] = eltrc[4];
277 : PagedImage<Complex> lastplane(elshp, calcCsys_p, "MOOBOO");
278 : lastplane.put(newAWConv(elblc, eltrc).nonDegenerate());
279 :
280 : //////
281 : }
282 : */
283 : // have to slice if not zero
284 : //cerr << "convFunc nelements " << convFunc_p.nelements() << endl;
285 0 : if (convFunc_p.nelements() == 0) {
286 :
287 0 : Int npix = min(newAWConv.shape()[0], newAWConv.shape()[1]);
288 : //cerr << "####npix " << npix << " " << 2*max(awsupport)*oversamp_p << " oversamp " << oversamp_p << endl;
289 0 : if(npix < (2*max(awsupport+1)*oversamp_p)){
290 0 : npix=2*(max(awsupport)+1)*oversamp_p;
291 : //cerr << "aft npix " << npix << " shape " << newAWConv.shape() << endl;
292 :
293 0 : IPosition elshp=newAWConv.shape();
294 0 : elshp[0]=npix;
295 0 : elshp[1]=npix;
296 0 : convFunc_p=Array<Complex>(elshp,Complex(0.0));
297 0 : wgtConvFunc_p=Array<Complex>(elshp,Complex(0.0));
298 0 : MathUtils::putMiddle(convFunc_p, newAWConv);
299 0 : MathUtils::putMiddle(wgtConvFunc_p, newWtConv);
300 :
301 :
302 0 : }
303 : else{
304 0 : npix=2*(max(awsupport)+1)*oversamp_p;
305 :
306 0 : convFunc_p = MathUtils::getMiddle(newAWConv, npix, npix);
307 0 : wgtConvFunc_p = MathUtils::getMiddle(newWtConv, npix, npix);
308 :
309 : }
310 0 : convSizes_p.resize(trc[0] + 1);
311 0 : convSizes_p.set(npix);
312 0 : convSupport_p.resize(trc[0]+1);
313 0 : convSupport_p = awsupport.row(nfreqs-1);
314 : } else {
315 : // Appending PA changing only
316 : //in case support is bigger for this PA.
317 :
318 0 : IPosition newshp = convFunc_p.shape();
319 :
320 0 : if (newshp[0] < (2 * (max(awsupport) + 1) * oversamp_p)){
321 : //cerr << "@@@@RESHAPING " << endl;
322 0 : IPosition elshp = newshp;
323 0 : elshp[0] = (2 * (max(awsupport) + 1) * oversamp_p);
324 0 : elshp[1] = (2 * (max(awsupport) + 1) * oversamp_p);
325 0 : Array<Complex> tempC(elshp, Complex(0.0));
326 0 : Array<Complex> tempW(elshp,Complex(0.0));
327 0 : MathUtils::putMiddle(tempC, convFunc_p);
328 0 : MathUtils::putMiddle(tempW,wgtConvFunc_p);
329 0 : convFunc_p.reference(tempC);
330 0 : wgtConvFunc_p.reference(tempW);
331 0 : }
332 : // asumming same freqs for now
333 0 : newshp(4) = startrow + awConv.shape()(4);
334 0 : Int npix = min(newAWConv.shape()[0], newAWConv.shape()[1]);
335 0 : if (npix <= newshp[0]) {
336 0 : IPosition blcadded(5, 0, 0, 0, 0, startrow);
337 0 : IPosition trcadded = newshp - 1;
338 :
339 0 : if(convFunc_p.shape()(4)< newshp[4]){
340 0 : convFunc_p.resize(newshp, True);
341 0 : wgtConvFunc_p.resize(newshp, True);
342 : }
343 0 : convFunc_p(blcadded, trcadded).set(0.0);
344 0 : wgtConvFunc_p(blcadded, trcadded).set(0.0);
345 0 : Array<Complex> c = convFunc_p(blcadded, trcadded);
346 0 : MathUtils::putMiddle(c, newAWConv);
347 0 : Array<Complex> d=wgtConvFunc_p(blcadded, trcadded);
348 0 : MathUtils::putMiddle(d, newWtConv);
349 :
350 0 : convSizes_p.resize(trc[0]+1, true);
351 0 : convSizes_p(blc, trc).set(newshp[0]);
352 0 : convSupport_p.resize(trc[0]+1, true);
353 0 : convSupport_p(blc, trc)= awsupport.row(nfreqs-1);
354 0 : } else {
355 0 : IPosition blcadded(5, 0, 0, 0, 0, startrow);
356 0 : if(newshp.product()>0 && newshp[0] < npix)
357 0 : npix=newshp[0];
358 0 : IPosition trcadded = newshp-1;
359 0 : if (convFunc_p.shape()(4) < newshp[4]) {
360 0 : convFunc_p.resize(newshp, True);
361 0 : wgtConvFunc_p.resize(newshp, True);
362 : }
363 0 : convFunc_p(blcadded, trcadded) =
364 0 : MathUtils::getMiddle(newAWConv, npix, npix);
365 0 : wgtConvFunc_p(blcadded, trcadded) =
366 0 : MathUtils::getMiddle(newWtConv, npix, npix);
367 0 : convSizes_p.resize(trc[0] + 1, true);
368 0 : convSizes_p(blc, trc).set(npix);
369 0 : convSupport_p.resize(trc[0] + 1, true);
370 0 : convSupport_p(blc, trc) = awsupport.row(nfreqs - 1);
371 0 : }
372 0 : }
373 : //cerr << "Shape at the end " << convFunc_p.shape() << endl;
374 0 : }
375 0 : Array<Complex>& AWConvFuncHolder::getConvFunc() {
376 0 : return convFunc_p;
377 :
378 : }
379 0 : Array<Complex>& AWConvFuncHolder::getWeightConvFunc() {
380 0 : return wgtConvFunc_p;
381 :
382 : }
383 0 : Vector<Int> AWConvFuncHolder::getConvSizes() {
384 0 : return convSizes_p;
385 : }
386 0 : Vector<Int> AWConvFuncHolder::getConvSupports() {
387 :
388 0 : return convSupport_p;
389 : }
390 :
391 0 : Array<Complex> &AWConvFuncHolder::getConvFuncHPG() { return convFuncHPG_p; }
392 0 : Array<Complex> &AWConvFuncHolder::getWeightConvFuncHPG() { return wgtConvFuncHPG_p; }
393 0 : void AWConvFuncHolder::resetHPGConvFuncs(const vi::VisBuffer2 &vb){
394 : // A given set for hph will hold all w's
395 0 : wValsHPG_p.resize();
396 0 : wValsHPG_p = wVals_p;
397 0 : freqValsHPG_p.resize(freqVals_p.nelements(), False);
398 0 : Double fmin = min(vb.getFrequencies(0));
399 0 : Double fmax = max(vb.getFrequencies(0));
400 0 : uint indx = 0;
401 0 : vector<bool> usedfreq(freqVals_p.nelements());
402 0 : std::fill(usedfreq.begin(), usedfreq.end(), false);
403 : //cerr << "fmin " << fmin << " fmax " << fmax << " freqs " << freqVals_p << " " << (freqVals_p[0] >= fmin) << " " <<(freqVals_p[0] <= fmax) << endl;
404 0 : if(freqVals_p.nelements() >1){
405 0 : for (uint k = 0; k < freqVals_p.nelements(); ++k) {
406 0 : if(freqVals_p[k] >= fmin && freqVals_p[k] <= fmax){
407 0 : freqValsHPG_p[indx] = freqVals_p[k];
408 0 : usedfreq[k] = true;
409 0 : ++indx;
410 : }
411 :
412 : }
413 0 : if(indx==0){ //some single channel spw will do this
414 0 : Double diffFreq=1e40;
415 0 : for (uint k = 0; k < freqVals_p.nelements(); ++k) {
416 0 : if(abs(fmax-freqVals_p[k]) < diffFreq){
417 0 : diffFreq=abs(fmax-freqVals_p[k]);
418 0 : freqValsHPG_p[0]=freqVals_p[k];
419 0 : std::fill(usedfreq.begin(),usedfreq.end(), false);
420 0 : usedfreq[k]=true;
421 0 : indx=1;
422 : }
423 : }
424 :
425 : }
426 : }
427 : else{
428 : //only one freq so it has to match
429 0 : usedfreq[0]=true;
430 0 : indx=1;
431 :
432 : }
433 0 : freqValsHPG_p.resize(indx, True);
434 0 : IPosition cshap = convFunc_p.shape();
435 0 : cshap[3] = indx;
436 0 : wgtConvFuncHPG_p.resize(cshap);
437 0 : convFuncHPG_p.resize(cshap);
438 : //cerr << "spw " << vb.spectralWindows()(0) << " CSHAP " << cshap << " orig " << convFunc_p.shape() << endl;
439 0 : IPosition blcin(5, 0);
440 0 : IPosition trcin = convFunc_p.shape() - 1;
441 0 : IPosition blcout(5, 0);
442 0 : IPosition trcout = cshap - 1;
443 0 : indx = 0;
444 0 : for (uint k = 0; k < freqVals_p.nelements(); ++k) {
445 0 : if(usedfreq[k]){
446 0 : blcin[3] = k;
447 0 : trcin[3] = k;
448 0 : blcout[3] = indx;
449 0 : trcout[3] = indx;
450 0 : for (uint j = 0; j < wVals_p.nelements(); ++j) {
451 0 : blcin[4] = j;
452 0 : trcin[4] = j;
453 0 : blcout[4] = j;
454 0 : trcout[4] = j;
455 0 : wgtConvFuncHPG_p(blcout, trcout) = wgtConvFunc_p(blcin, trcin);
456 0 : convFuncHPG_p(blcout, trcout) = convFunc_p(blcin, trcin);
457 : }
458 0 : ++indx;
459 : }
460 : }
461 0 : }
462 :
463 : /////////////////////
464 0 : void AWConvFuncHolder::getConvFuncs(Vector<Int> &polMap, Vector<Int> &chanMap,
465 : Vector<Int> &rowMap,
466 : Array<Complex> &convFunc,
467 : Array<Complex> &wgtConvFunc,
468 : const vi::VisBuffer2 &vb,
469 : const Matrix<Double> &rotuvw,
470 : const Vector<Double> & interpFreqs,
471 : const Bool predictMode,
472 : const bool ispsf) {
473 :
474 0 : Vector<Int> cmap;
475 0 : Vector<Int> pmap;
476 0 : Vector<Int> rmap;
477 :
478 0 : getConvIndices(pmap, cmap, rmap, vb, rotuvw, interpFreqs, predictMode, ispsf);
479 : //cerr << "pmap "<< pmap << endl;
480 : //cerr << "MIN Max rmap" << min(rmap) << " " << max(rmap) << endl;
481 0 : std::vector<Int> pmapused = pmap.tovector();
482 : {
483 0 : std::sort(pmapused.begin(), pmapused.end());
484 0 : auto last = std::unique(pmapused.begin(), pmapused.end());
485 0 : pmapused.erase(last, pmapused.end());
486 : }
487 0 : std::vector<Int> cmapused = cmap.tovector();
488 : {
489 0 : std::sort(cmapused.begin(), cmapused.end());
490 0 : auto last = std::unique(cmapused.begin(), cmapused.end());
491 0 : cmapused.erase(last, cmapused.end());
492 : }
493 0 : std::vector<Int> rmapused = rmap.tovector();
494 : {
495 0 : std::sort(rmapused.begin(), rmapused.end());
496 0 : auto last = std::unique(rmapused.begin(), rmapused.end());
497 0 : rmapused.erase(last, rmapused.end());
498 : }
499 : {
500 0 : vector<Int> cpRmapUsed = rmapused;
501 :
502 : // lets move the -ve values to the end -ve means -w which means we have to
503 : // conjugate the plane
504 0 : vector<int>::iterator it = remove_if(rmapused.begin(), rmapused.end(),
505 0 : [](const int i) { return i < 0; });
506 0 : rmapused.erase(it, rmapused.end());
507 0 : for (auto cit = cpRmapUsed.rbegin(); cit != cpRmapUsed.rend(); ++cit){
508 0 : if(*cit <0)
509 0 : rmapused.push_back(*cit);
510 : }
511 0 : }
512 : //cerr << "#####rmapused " << rmapused << endl;
513 0 : IPosition shp(5, convFunc_p.shape()[0], convFunc_p.shape()[1],
514 0 : pmapused.size(), cmapused.size(), rmapused.size());
515 0 : polMap.resize(pmap.shape());
516 0 : for (uint j = 0; j < polMap.nelements(); ++j) {
517 0 : for (uint k = 0; k < pmapused.size(); ++k) {
518 0 : if (pmap[j] == pmapused[k]) {
519 0 : polMap[j] = k;
520 : }
521 : }
522 : }
523 0 : chanMap.resize(cmap.shape());
524 : //std::vector<int>cindex(cmapused.size());
525 0 : for (uint j = 0; j < chanMap.nelements(); ++j) {
526 0 : for (uint k = 0; k < cmapused.size(); ++k) {
527 0 : if (cmap[j] == cmapused[k]){
528 0 : chanMap[j] = k;
529 : }
530 : }
531 : }
532 0 : rowMap.resize(rmap.shape());
533 0 : for (uint j = 0; j < rowMap.nelements(); ++j) {
534 0 : for (uint k = 0; k < rmapused.size(); ++k) {
535 0 : if (rmap[j] == rmapused[k]){
536 : //rowmap is -ve for -ve w
537 0 : rowMap[j] = k;
538 : }
539 : }
540 : }
541 : //cerr << "old rmapused" << Vector<int>(rmapused) << " cmap " << Vector<int>(cmapused) << " pmap " << Vector<int>(pmapused) << endl;
542 0 : convFunc.resize(shp);
543 0 : wgtConvFunc.resize(shp);
544 0 : IPosition inblc(5, 0, 0, 0, 0, 0);
545 0 : IPosition intrc(5, shp[0] - 1, shp[1] - 1, 0, 0, 0);
546 0 : IPosition outblc(5, 0, 0, 0, 0, 0);
547 0 : IPosition outtrc(5, shp[0] - 1, shp[1] - 1, 0, 0, 0);
548 :
549 0 : for (uint r = 0; r < rmapused.size(); ++r){
550 0 : inblc[4] = abs(rmapused[r]);
551 0 : intrc[4] = abs(rmapused[r]);
552 0 : outblc[4] = r;
553 0 : outtrc[4] = r;
554 0 : for (uint c = 0; c < cmapused.size(); ++c) {
555 0 : inblc[3] = cmapused[c];
556 0 : intrc[3] = cmapused[c];
557 0 : outblc[3] = c;
558 0 : outtrc[3] = c;
559 0 : for (uint p = 0; p < pmapused.size(); ++p) {
560 0 : inblc[2] = pmapused[p];
561 0 : intrc[2] = pmapused[p];
562 0 : outblc[2] = p;
563 0 : outtrc[2] = p;
564 : // rowmap is -ve for -ve w
565 0 : convFunc(outblc, outtrc) = rmapused[r] > 0
566 0 : ? convFunc_p(inblc, intrc)
567 0 : : conj(convFunc_p(inblc, intrc));
568 0 : wgtConvFunc(outblc, outtrc) = wgtConvFunc_p(inblc, intrc);
569 : }
570 : }
571 : }
572 0 : }
573 : //////////////////////
574 :
575 :
576 : //////////////////////
577 :
578 0 : void AWConvFuncHolder::getConvIndices(Vector<Int>& polMap, Vector<Int>& chanMap, Vector<Int>& rowMap, const vi::VisBuffer2& vb, const Matrix<Double>& rotuvw, const Vector<Double>& interpFreqs,
579 : const Bool predictMode, const bool ispsf) {
580 : // Lets do the polmap
581 0 : Vector<Stokes::StokesTypes> visPolMap(vb.getCorrelationTypesSelected());
582 0 : polMap.resize(visPolMap.nelements());
583 0 : polMap.set(-1);
584 0 : if (!dosquint_p)
585 0 : polMap.set(0);
586 : else{
587 0 : for (uint k = 0; k < polMap.nelements(); ++k) {
588 0 : for (uint j = 0; j < polVals_p.nelements(); ++j) {
589 0 : if (visPolMap[k] == polVals_p[j])
590 0 : polMap[k] = j;
591 : }
592 : }
593 : }
594 : // Lets do chanMap
595 0 : chanMap.resize(interpFreqs.nelements());
596 0 : chanMap.set(-1);
597 0 : Vector<Double>visFreq = interpFreqs;
598 0 : for (uint k = 0; k < chanMap.nelements(); ++k) {
599 0 : Double minDiff = 1e40;
600 0 : Int indexF = -1;
601 0 : for (uint j = 0; j < freqVals_p.nelements(); ++j) {
602 0 : if (fabs(freqVals_p[j] -visFreq[k]) < minDiff) {
603 0 : minDiff = fabs(freqVals_p[j] -visFreq[k]);
604 0 : indexF = j;
605 : }
606 : }
607 0 : chanMap[k] = indexF;
608 :
609 : }
610 : // Now to the complicated rowMap
611 0 : Vector<Int> antPairIndex(vb.nRows(), -1);
612 0 : Vector<Int> wIndex(vb.nRows(), -1);
613 0 : Vector<Int> paIndex(vb.nRows(), -1);
614 : //Assuming pa is the same for this vb which is usually associated with time.
615 : // using utils.cc global function
616 0 : Double paval = refim::getPA(vb);
617 0 : if(predictMode){
618 0 : paval -= M_PI;
619 0 : paval = atan2(sin(paval), cos(paval));
620 : }
621 0 : Int tmpPAInd = -1;
622 0 : Double minDiff = 1e40;
623 0 : for (uint k = 0; k <paVals_p.nelements(); ++k) {
624 0 : if (fabs(paval-paVals_p[k]) < minDiff) {
625 0 : tmpPAInd = k;
626 0 : minDiff = fabs(paval-paVals_p[k]);
627 : }
628 : }
629 0 : paIndex.set(tmpPAInd);
630 : //cerr << "paVal " << paval << " predictMode " << predictMode << endl;
631 : // For antenna pairs ..for homogenous arrays only one pair is necessary
632 0 : if ( (antpairVals_p.nelements() == 1) && antpairVals_p[0] == std::pair<int, int>(-1, -1)) {
633 0 : antPairIndex.set(0);
634 : }
635 : else{
636 0 : for (uint k = 0; k < vb.nRows(); ++k) {
637 0 : std::pair<int, int> antpair = std::make_pair(vb.antenna1()[k], vb.antenna2()[k]);
638 0 : for (uint j = 0; j < antpairVals_p.nelements();++j ) {
639 0 : if (antpairVals_p[j] == antpair)
640 0 : antPairIndex[k] = j;
641 : }
642 :
643 : }
644 : }
645 0 : Double invlamda = mean(vb.getFrequencies(0))/C::c;
646 0 : for (uint k = 0; k < vb.nRows();++k) {
647 0 : minDiff = 1e40;
648 0 : Int tmpWInd = -1;
649 0 : Double w = ispsf ? 0 : rotuvw.row(2)[k] *invlamda;
650 0 : for (uint j =0; j < wVals_p.nelements();++j ) {
651 0 : if (fabs(fabs(w)-wVals_p[j]) < minDiff) {
652 0 : minDiff = fabs(fabs(w)-wVals_p[j]);
653 0 : tmpWInd = j;
654 : }
655 : }
656 0 : wIndex[k] = (w > 0)? -tmpWInd : tmpWInd;
657 : }
658 : // Now lets search for combination of all 3
659 0 : rowMap.resize(vb.nRows());
660 0 : for (uint k = 0; k < vb.nRows();++k) {
661 0 : for (uint j = 0; j < rowAxisWVals_p.nelements(); ++j) {
662 0 : if ( (abs(wIndex[k]) == rowAxisWVals_p[j]) && (paIndex[k] == rowAxisPAVals_p[j]) && (antPairIndex[k] == rowAxisAntennaPair_p[j]) ) {
663 0 : rowMap[k] = wIndex[k] >= 0 ? j : -j;
664 : }
665 : }
666 :
667 : }
668 : // A little dab will d'ya
669 :
670 0 : }
671 :
672 :
673 0 : void AWConvFuncHolder::getConvIndicesHPG(Vector<Int> &polMap, Vector<Int> &chanMap,
674 : Vector<Int> &rowMap,
675 : const vi::VisBuffer2 &vb,
676 : const Matrix<Double> &rotuvw) {
677 0 : Vector<Stokes::StokesTypes> visPolMap(vb.getCorrelationTypesSelected());
678 0 : polMap.resize(visPolMap.nelements());
679 : //HPG is doing I single plane gridding
680 0 : polMap.set(0);
681 : // Lets do chanMap
682 0 : chanMap.resize(vb.nChannels());
683 0 : chanMap.set(-1);
684 0 : Vector<Double> visFreq = vb.getFrequencies(0);
685 0 : for (uint k = 0; k < chanMap.nelements(); ++k) {
686 0 : Double minDiff = 1e40;
687 0 : Int indexF = -1;
688 0 : for (uint j = 0; j < freqValsHPG_p.nelements(); ++j) {
689 0 : if (fabs(freqValsHPG_p[j] - visFreq[k]) < minDiff) {
690 0 : minDiff = fabs(freqValsHPG_p[j] - visFreq[k]);
691 0 : indexF = j;
692 : }
693 : }
694 0 : chanMap[k] = indexF;
695 : }
696 : //As there is no PA or antenna pair to deal with HPG...windex should be rowMap
697 0 : Vector<Int> wIndex(vb.nRows(), 0);
698 0 : Double invlamda = mean(vb.getFrequencies(0)) / C::c;
699 0 : for (uint k = 0; k < vb.nRows(); ++k) {
700 0 : Double minDiff = 1e40;
701 0 : Int tmpWInd = -1;
702 0 : Double w = rotuvw.row(2)[k] * invlamda;
703 0 : for (uint j = 0; j < wValsHPG_p.nelements(); ++j) {
704 0 : if (fabs(fabs(w) - wValsHPG_p[j]) < minDiff) {
705 0 : minDiff = fabs(fabs(w) - wValsHPG_p[j]);
706 0 : tmpWInd = j;
707 : }
708 : }
709 0 : wIndex[k] = tmpWInd;
710 : }
711 0 : rowMap.resize();
712 0 : rowMap = wIndex;
713 : //cerr << "FID " << vb.fieldId()(0) << " SPID " << vb.spectralWindows()(0) << " winDex " << wIndex << endl;
714 0 : }
715 :
716 0 : Vector<Double> AWConvFuncHolder::getPointingPhaseShift(
717 : const vi::VisBuffer2 &vb, bool usePointingTable) {
718 0 : Bool hasValidPointing = False;
719 0 : if (vbutil_p.use_count() == 0)
720 0 : vbutil_p = std::make_shared<VisBufferUtil>(vb);
721 0 : MDirection ant1PointVal;
722 0 : if(Table::isReadable(vb.ms().pointingTableName())){
723 0 : hasValidPointing=usePointingTable && (vb.ms().pointing().nrow() >0);
724 : }
725 0 : DirectionCoordinate dc=outcsys_p.directionCoordinate(0);
726 :
727 0 : if(hasValidPointing){
728 : //ant1PointingCache_p[val]=vb.direction1()[0];
729 0 : ant1PointVal=vbutil_p->getPointingDir(vb, vb.antenna1()(0), 0, dc.directionType());
730 : }
731 : else
732 0 : ant1PointVal=vbutil_p->getPhaseCenter(vb);
733 0 : MSColumns mscol(vb.ms());
734 0 : String tel;
735 0 : if (vb.subtableColumns().observation().nrow() > 0) {
736 0 : tel =vb.subtableColumns().observation().telescopeName()(mscol.observationId()(0));
737 : }
738 : MEpoch::Types timeMType;
739 0 : casacore::Unit timeUnit;
740 0 : timeMType=MEpoch::castType(mscol.timeMeas()(0).getRef().getType());
741 0 : timeUnit=Unit(mscol.timeMeas().measDesc().getUnits()(0).getName());
742 0 : MPosition pos;
743 0 : MDirection dirOnImage;
744 0 : MeasTable::Observatory(pos,tel);
745 : //need to conver antpoint frame to image frame
746 0 : if(dc.directionType() != MDirection::castType(ant1PointVal.getRef().getType())){
747 :
748 0 : MEpoch timenow(Quantity(vb.time()(0), timeUnit), timeMType);
749 0 : MeasFrame pointFrame(timenow, pos);
750 0 : MDirection::Ref elRef(dc.directionType(), pointFrame);
751 0 : dirOnImage=MDirection::Convert(ant1PointVal, elRef)();
752 :
753 0 : }
754 : else{
755 0 : dirOnImage=ant1PointVal;
756 :
757 : }
758 :
759 0 : Vector<Double> thePix(2);
760 0 : dc.toPixel(thePix, dirOnImage);
761 : //shift from center
762 0 : thePix(0) = thePix(0) - Double(nx_p / 2);
763 0 : thePix(1) = thePix(1) - Double(ny_p / 2);
764 :
765 : //phase gradient per pixel to apply
766 0 : thePix(0) = -thePix(0)*2.0*M_PI/Double(nx_p)/Double(oversamp_p);
767 0 : thePix(1) = -thePix(1) * 2.0 * M_PI / Double(ny_p) / Double(oversamp_p);
768 : //cerr << std::setprecision(12) << "fid " << vb.fieldId()(0) << " POINT shift " << thePix << endl;
769 :
770 0 : return thePix;
771 :
772 :
773 :
774 :
775 :
776 0 : }
777 :
778 : } // # namespace refim ends
779 : }//namespace CASA ends
780 :
781 :
|
