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 <msvis/MSVis/VisBuffer2.h>
34 : #include <msvis/MSVis/VisibilityIterator2.h>
35 : #include <synthesis/TransformMachines2/AWConvFunc.h>
36 : #include <synthesis/TransformMachines2/EVLAAperture.h>
37 : #include <synthesis/TransformMachines2/AWConvFuncHolder.h>
38 :
39 : namespace casa { //# CASA namespace
40 : namespace refim { //# namespace refactor imaging
41 :
42 : using namespace casacore;
43 : using namespace casa;
44 : using namespace casa::refim;
45 : using namespace std;
46 :
47 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) {
48 :
49 0 : convFunc_p.resize();
50 0 : wgtConvFunc_p.resize();
51 0 : polVals_p.resize();
52 0 : freqVals_p.resize();
53 0 : wVals_p.resize();
54 0 : paVals_p.resize();
55 0 : antpairVals_p.resize();
56 0 : rowAxisWVals_p.resize();
57 0 : rowAxisPAVals_p.resize();
58 0 : rowAxisAntennaPair_p.resize();
59 0 : convSizes_p.resize();
60 0 : convSupport_p.resize();
61 0 : aterm_p = std::make_shared<refim::EVLAAperture>();
62 0 : aterm_p->cacheVBInfo(obs, 25.0);
63 :
64 0 : }
65 0 : AWConvFuncHolder::AWConvFuncHolder(const AWConvFuncHolder& other) {
66 0 : operator = (other);
67 0 : }
68 0 : AWConvFuncHolder& AWConvFuncHolder::operator=(const AWConvFuncHolder& other) {
69 0 : if (this != &other) {
70 0 : nx_p = other.nx_p;
71 0 : ny_p = other.ny_p;
72 0 : oversamp_p = other.oversamp_p;
73 0 : outcsys_p = other.outcsys_p;
74 0 : calcNpix_p = other.calcNpix_p;
75 0 : calcCsys_p = other.calcCsys_p;
76 0 : dosquint_p = other.dosquint_p;
77 0 : painc_p = other.painc_p;
78 0 : convFunc_p.resize();
79 0 : convFunc_p = other.convFunc_p;
80 0 : wgtConvFunc_p.resize();
81 0 : wgtConvFunc_p = other.wgtConvFunc_p;
82 0 : polVals_p.resize();
83 0 : polVals_p = other.polVals_p;
84 0 : freqVals_p.resize();
85 0 : freqVals_p = other.freqVals_p;
86 0 : wVals_p.resize();
87 0 : wVals_p = other.wVals_p;
88 0 : paVals_p.resize();
89 0 : paVals_p = other.paVals_p;
90 0 : antpairVals_p.resize();
91 0 : antpairVals_p = other.antpairVals_p;
92 0 : rowAxisWVals_p.resize();
93 0 : rowAxisWVals_p = other.rowAxisWVals_p;
94 0 : rowAxisPAVals_p.resize();
95 0 : rowAxisPAVals_p = other.rowAxisPAVals_p;
96 0 : rowAxisAntennaPair_p.resize();
97 0 : rowAxisAntennaPair_p = other.rowAxisAntennaPair_p;
98 0 : convSizes_p.resize();
99 0 : convSizes_p = other.convSizes_p;
100 0 : convSupport_p.resize();
101 0 : convSupport_p = other.convSupport_p;
102 0 : aterm_p = other.aterm_p;
103 :
104 :
105 :
106 :
107 : }
108 0 : return *this;
109 : }
110 0 : bool AWConvFuncHolder::addConvFunc(const casacore::Vector<casacore::Double>& freqs, const casacore::Vector<Double>& wVals, const casacore::Double& paMax) {
111 :
112 0 : Vector<Double> freqsToCalc;
113 0 : if (freqVals_p.nelements() == 0) {
114 0 : freqVals_p = freqs;
115 0 : freqsToCalc = freqs;
116 :
117 : }
118 : else{
119 : //append new freqs freqVals_p;
120 : // assign missing freqs to freqsToCalc
121 0 : freqsToCalc =freqs;
122 : }
123 0 : if (wVals_p.nelements() == 0) {
124 : // make sure first wval is 0;
125 0 : wVals_p = wVals;
126 :
127 : }
128 0 : else if (wVals_p.nelements() != wVals.nelements())
129 0 : throw(AipsError("Cannot change W terms length right now"));
130 0 : if (!dosquint_p) {
131 0 : paVals_p.resize(1);
132 0 : paVals_p[0] = 0.0;
133 : }
134 : else{
135 0 : cerr << "paMax " << paMax << " painc " << painc_p << endl;
136 0 : Vector<Double> pavals(int(std::ceil(2*paMax/painc_p)));
137 : //setting pavals from -paMax to paMax
138 0 : for (uint k = 0; k < pavals.nelements(); ++k )
139 0 : pavals[k] = double(k) *painc_p-paMax;
140 0 : if (paVals_p.nelements() == 0)
141 0 : paVals_p = pavals;
142 0 : else if ((paVals_p.nelements()) != pavals.nelements())
143 0 : throw(AipsError("Cannot change number of PA's in between"));
144 0 : }
145 0 : std::shared_ptr<refim::WPConvFunc>wptr;
146 0 : AWConvFunc a(aterm_p, wptr);
147 0 : Array<Complex> aWConv;
148 0 : Array<Complex> aWwtconv;
149 0 : Matrix<Int> awSupport;
150 0 : calcCsys_p = outcsys_p;
151 0 : calcNpix_p = min(nx_p, ny_p);
152 0 : cerr << "PAVALS " << paVals_p << " dosquint " << dosquint_p << endl;
153 0 : for (uint k=0; k<paVals_p.nelements(); ++k){
154 0 : a.makeAWConvFunc(aWConv, aWwtconv,calcCsys_p,awSupport, calcNpix_p, freqsToCalc, wVals_p, dosquint_p, paVals_p[k]);
155 :
156 :
157 : //append arrays and indices
158 0 : appendConvFuncs(aWConv, aWwtconv, awSupport, freqsToCalc, paVals_p[k]);
159 : }
160 :
161 :
162 :
163 :
164 0 : return true;;
165 0 : }
166 0 : void AWConvFuncHolder::appendConvFuncs(const Array<Complex>& awConv, const Array<Complex>& aWwtConv, const Matrix<Int>& awsupport, const Vector<Double>& newFreqs, const Double paVal) {
167 0 : Int nfreqs = awConv.shape()[3];
168 : // Make sure the polVals are in the stokes used in making convfun
169 0 : polVals_p.resize(4);
170 : // for now antpairVals are for all pair of possible antenna combination
171 0 : antpairVals_p.resize(1);
172 0 : antpairVals_p[0] = std::pair<int, int>(-1, -1);
173 :
174 0 : convertArray(polVals_p, calcCsys_p.stokesCoordinate(1).stokes());
175 0 : IPosition blc = freqVals_p.shape();
176 0 : IPosition trc(1, blc[0]+nfreqs-1);
177 0 : if (freqVals_p.nelements() != 0) {
178 0 : if (!allEQ(freqVals_p, newFreqs))
179 0 : throw(AipsError("Not implemented appending different freqs yet to conv"));
180 : }
181 : else{
182 0 : freqVals_p.resize(freqVals_p.nelements()+nfreqs, True);
183 0 : freqVals_p(blc, trc) = newFreqs;
184 : }
185 0 : Int indPA = -1;
186 0 : for (uint k = 0; k < paVals_p.nelements(); ++k) {
187 0 : if (fabs(paVal-paVals_p[k]) < 1e-4)
188 0 : indPA = k;
189 : }
190 0 : if (indPA <0) {
191 0 : paVals_p.resize(paVals_p.nelements()+1, True);
192 0 : indPA = paVals_p.nelements()-1;
193 0 : paVals_p[indPA] = paVal;
194 :
195 : }
196 0 : blc[0] = rowAxisPAVals_p.nelements();
197 0 : trc[0] = blc[0]+wVals_p.nelements()-1;
198 0 : rowAxisPAVals_p.resize(trc[0]+1, True);
199 0 : rowAxisPAVals_p(blc, trc).set(indPA);
200 0 : rowAxisWVals_p.resize(trc[0]+1, True);
201 0 : Vector<Int> waxis(wVals_p.nelements());
202 0 : indgen(waxis);
203 0 : rowAxisWVals_p(blc, trc) = waxis;
204 0 : rowAxisAntennaPair_p.resize(trc[0]+1, True);
205 0 : rowAxisAntennaPair_p(blc, trc).set(0);
206 : /// Let us rescale convolution function to match nx, ny and incr of image that makes
207 : /// uvgrid
208 0 : Float factorX=fabs(calcCsys_p.increment()(0)/outcsys_p.increment()(0));
209 0 : Float factorY=fabs(calcCsys_p.increment()(1)/outcsys_p.increment()(1));
210 : // cerr << "####Factor " << factorX << " " << factorY << endl;
211 0 : factorX = Float(nx_p) *Float(oversamp_p)/Float(calcNpix_p)/factorX;
212 0 : factorY = Float(ny_p) *Float(oversamp_p)/Float(calcNpix_p)/factorY;
213 : // cerr << "factors " << factorX << " " << factorY << "nx, ny" << nx_p << " " << ny_p << " calcNpix " << calcNpix_p << " oversamp " << oversamp_p << endl;
214 0 : MathUtils m;
215 0 : Array<Complex>newAWConv = m.resampleViaFFT(awConv, factorX, factorY);
216 0 : Array<Complex> newWtConv = m.resampleViaFFT(aWwtConv, factorX, factorY);
217 0 : Float correcfac = float(awConv.shape()(0) *awConv.shape()(1) *oversamp_p*oversamp_p)/float(newAWConv.shape()(0) *newAWConv.shape()(1));
218 : //cerr << "correcfac " << correcfac << " " << 1.0/correcfac << endl;
219 0 : newAWConv *= correcfac;
220 0 : newWtConv *= correcfac;
221 : /*{
222 : ////TESTOO
223 : IPosition elshp = newAWConv.shape().getFirst(4);
224 : IPosition elblc(5, 0);
225 :
226 : IPosition eltrc = newAWConv.shape()-1;
227 : elblc[4] = eltrc[4];
228 : PagedImage<Complex> lastplane(elshp, calcCsys_p, "MOOBOO");
229 : lastplane.put(newAWConv(elblc, eltrc).nonDegenerate());
230 :
231 : //////
232 : }*/
233 : // have to slice if not zero
234 0 : if (convFunc_p.nelements() == 0) {
235 0 : Int npix = min(newAWConv.shape()[0], newAWConv.shape()[1]);
236 : //cerr << "npix " << npix << " " << 2*max(awsupport)*oversamp_p << endl;
237 0 : if(npix <= 2*max(awsupport)*oversamp_p){
238 0 : npix=2*(max(awsupport)+1)*oversamp_p;
239 0 : cerr << "aft npix " << npix << endl;
240 0 : IPosition elshp=newAWConv.shape();
241 0 : elshp[0]=npix;
242 0 : elshp[1]=npix;
243 0 : convFunc_p=Array<Complex>(elshp,Complex(0.0));
244 0 : wgtConvFunc_p=Array<Complex>(elshp,Complex(0.0));
245 0 : MathUtils::putMiddle(convFunc_p, newAWConv);
246 0 : MathUtils::putMiddle(wgtConvFunc_p, newWtConv);
247 :
248 :
249 0 : }
250 : else{
251 0 : convFunc_p = MathUtils::getMiddle(newAWConv, npix, npix);
252 0 : wgtConvFunc_p = MathUtils::getMiddle(newWtConv, npix, npix);
253 : }
254 0 : convSizes_p.resize(trc[0]+1);
255 0 : convSizes_p.set(npix);
256 0 : convSupport_p.resize(trc[0]+1);
257 0 : convSupport_p = awsupport.row(nfreqs-1);
258 : }
259 : else{
260 : // Appending PA changing only
261 0 : IPosition newshp = convFunc_p.shape();
262 : // asumming same freqs for now
263 0 : newshp(4) = newshp(4)+awConv.shape()(4);
264 0 : Int npix = min(newAWConv.shape()[0], newAWConv.shape()[1]);
265 0 : if (npix != newshp[0])
266 : {
267 :
268 0 : cerr << "npix is not the same for a different PA" << endl;
269 : }
270 : else{
271 0 : IPosition blcadded(5, 0, 0, 0, 0, convFunc_p.shape()[4]);
272 0 : IPosition trcadded = newshp-1;
273 0 : convFunc_p.resize(newshp, True);
274 0 : wgtConvFunc_p.resize(newshp, True);
275 0 : convFunc_p(blcadded, trcadded) = MathUtils::getMiddle(newAWConv, npix, npix);
276 0 : wgtConvFunc_p(blcadded, trcadded) = MathUtils::getMiddle(newWtConv, npix, npix);
277 0 : convSizes_p.resize(trc[0]+1);
278 0 : convSizes_p(blc, trc).set(npix);
279 0 : convSupport_p.resize(trc[0]+1);
280 0 : convSupport_p(blc, trc)= awsupport.row(nfreqs-1);
281 :
282 :
283 0 : }
284 0 : }
285 :
286 :
287 0 : }
288 0 : Array<Complex>& AWConvFuncHolder::getConvFunc() {
289 0 : return convFunc_p;
290 :
291 : }
292 0 : Array<Complex>& AWConvFuncHolder::getWeightConvFunc() {
293 0 : return wgtConvFunc_p;
294 :
295 : }
296 0 : Vector<Int> AWConvFuncHolder::getConvSizes() {
297 0 : return convSizes_p;
298 : }
299 0 : Vector<Int> AWConvFuncHolder::getConvSupports() {
300 :
301 0 : return convSupport_p;
302 : }
303 0 : void AWConvFuncHolder::getConvIndices(Vector<Int>& polMap, Vector<Int>& chanMap, Vector<Int>& rowMap, const vi::VisBuffer2& vb, const Matrix<Double>& rotuvw) {
304 : // Lets do the polmap
305 0 : Vector<Stokes::StokesTypes> visPolMap(vb.getCorrelationTypesSelected());
306 0 : polMap.resize(visPolMap.nelements());
307 0 : polMap.set(-1);
308 0 : if (!dosquint_p)
309 0 : polMap.set(0);
310 : else{
311 0 : for (uint k = 0; k < polMap.nelements(); ++k) {
312 0 : for (uint j = 0; j < polVals_p.nelements(); ++j) {
313 0 : if (visPolMap[k] == polVals_p[j])
314 0 : polMap[k] = j;
315 : }
316 : }
317 : }
318 : // Lets do chanMap
319 0 : chanMap.resize(vb.nChannels());
320 0 : chanMap.set(-1);
321 0 : Vector<Double>visFreq = vb.getFrequencies(0);
322 0 : for (uint k = 0; k < chanMap.nelements(); ++k) {
323 0 : Double minDiff = 1e40;
324 0 : Int indexF = -1;
325 0 : for (uint j = 0; j < freqVals_p.nelements(); ++j) {
326 0 : if (fabs(freqVals_p[j] -visFreq[k]) < minDiff) {
327 0 : minDiff = fabs(freqVals_p[j] -visFreq[k]);
328 0 : indexF = j;
329 : }
330 : }
331 0 : chanMap[k] = indexF;
332 :
333 : }
334 : // Now to the complicated rowMap
335 0 : Vector<Int> antPairIndex(vb.nRows(), -1);
336 0 : Vector<Int> wIndex(vb.nRows(), -1);
337 0 : Vector<Int> paIndex(vb.nRows(), -1);
338 : //Assuming pa is the same for this vb which is usually associated with time.
339 : // using utils.cc global function
340 0 : Double paval = refim::getPA(vb);
341 0 : Int tmpPAInd = -1;
342 0 : Double minDiff = 1e40;
343 0 : for (uint k = 0; k <paVals_p.nelements(); ++k) {
344 0 : if (fabs(paval-paVals_p[k]) < minDiff) {
345 0 : tmpPAInd = k;
346 0 : minDiff = fabs(paval-paVals_p[k]);
347 : }
348 : }
349 0 : paIndex.set(tmpPAInd);
350 : // For antenna pairs ..for homogenous arrays only one pair is necessary
351 0 : if ( (antpairVals_p.nelements() == 1) && antpairVals_p[0] == std::pair<int, int>(-1, -1)) {
352 0 : antPairIndex.set(0);
353 : }
354 : else{
355 0 : for (uint k = 0; k < vb.nRows(); ++k) {
356 0 : std::pair<int, int> antpair = std::make_pair(vb.antenna1()[k], vb.antenna2()[k]);
357 0 : for (uint j = 0; j < antpairVals_p.nelements();++j ) {
358 0 : if (antpairVals_p[j] == antpair)
359 0 : antPairIndex[k] = j;
360 : }
361 :
362 : }
363 : }
364 0 : Double invlamda = mean(vb.getFrequencies(0))/C::c;
365 0 : for (uint k = 0; k < vb.nRows();++k) {
366 0 : minDiff = 1e40;
367 0 : Int tmpWInd = -1;
368 0 : Double w = rotuvw.row(2)[k] *invlamda;
369 0 : for (uint j =0; j < wVals_p.nelements();++j ) {
370 0 : if (fabs(fabs(w)-wVals_p[j]) < minDiff) {
371 0 : minDiff = fabs(fabs(w)-wVals_p[j]);
372 0 : tmpWInd = j;
373 : }
374 : }
375 0 : wIndex[k] = (w > 0)? -tmpWInd : tmpWInd;
376 : }
377 : // Now lets search for combination of all 3
378 0 : rowMap.resize(vb.nRows());
379 0 : for (uint k = 0; k < vb.nRows();++k) {
380 0 : for (uint j = 0; j < rowAxisWVals_p.nelements(); ++j) {
381 0 : if ( (abs(wIndex[k]) == rowAxisWVals_p[j]) && (paIndex[k] == rowAxisPAVals_p[j]) && (antPairIndex[k] == rowAxisAntennaPair_p[j]) ) {
382 0 : rowMap[k] = wIndex[k] > 0 ? j : -j;
383 : }
384 : }
385 :
386 : }
387 : // A little dab will d'ya
388 :
389 0 : }
390 : }//# namespace refim ends
391 : }//namespace CASA ends
392 :
393 :
|
