Line data Source code
1 : //# PolAverageTVI.h: This file contains the implementation of the PolAverageTVI class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 : #include <mstransform/TVI/PolAverageTVI.h>
23 :
24 : #include <casacore/casa/Arrays/Cube.h>
25 : #include <casacore/casa/Arrays/Matrix.h>
26 : #include <casacore/casa/Arrays/Vector.h>
27 : #include <casacore/casa/BasicSL/String.h>
28 : #include <casacore/casa/Logging/LogIO.h>
29 : #include <casacore/casa/Containers/Record.h>
30 : #include <casacore/casa/Exceptions/Error.h>
31 : #include <casacore/casa/Arrays/ArrayIter.h>
32 : #include <casacore/measures/Measures/Stokes.h>
33 : #include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
34 : #include <casacore/ms/MeasurementSets/MSPolColumns.h>
35 :
36 : #include <msvis/MSVis/VisBufferComponents2.h>
37 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
38 :
39 : #include <mstransform/TVI/UtilsTVI.h>
40 :
41 : using namespace casacore;
42 :
43 : namespace {
44 : template<class T>
45 23224320 : inline T replaceFlaggedDataWithZero(T v, Bool b) {
46 23224320 : return ((b) ? T(0) : v);
47 : }
48 :
49 : struct StokesTransformation {
50 : public:
51 : template<class T>
52 8496 : inline static void transformData(Cube<T> const &dataIn,
53 : Cube<Bool> const &flagIn, Int pid0, Int pid1, Cube<T> &dataOut) {
54 8496 : if (dataIn.empty()) {
55 0 : dataOut.resize();
56 0 : return;
57 : }
58 8496 : auto const cubeShape = dataIn.shape();
59 8496 : IPosition const newShape(3, 1, cubeShape[1], cubeShape[2]);
60 : // size_t const npol = cubeShape[0];
61 8496 : size_t const nchan = cubeShape[1];
62 8496 : size_t const nrow = cubeShape[2];
63 8496 : size_t const nelem = cubeShape[1] * cubeShape[2];
64 8496 : Cube<T> transformedData(newShape, T(0));
65 8496 : Cube<Float> weightSum(newShape, 0.0f);
66 :
67 8496 : Int pols[] = { pid0, pid1 };
68 25488 : for (size_t i = 0; i < 2; ++i) {
69 16992 : Int ipol = pols[i];
70 16992 : IPosition start(3, ipol, 0, 0);
71 16992 : IPosition end(3, ipol, nchan - 1, nrow - 1);
72 16992 : auto dslice = dataIn(start, end);
73 16992 : auto fslice = flagIn(start, end);
74 16992 : Array<Float> weight(dslice.shape());
75 16992 : Array<T> weightedData(dslice.shape());
76 16992 : arrayContTransform(dslice, fslice, weightedData,
77 : replaceFlaggedDataWithZero<T>);
78 11266272 : arrayContTransform(fslice, weight, [](Bool b) {
79 11249280 : return ((b) ? 0.0f : 1.0f);
80 : });
81 16992 : transformedData += weightedData;
82 16992 : weightSum += weight;
83 : }
84 :
85 : // transformedData, transformedFlag and nAccumulated should be contiguous array
86 8496 : auto p_tdata = transformedData.data();
87 8496 : auto p_wsum = weightSum.data();
88 :
89 5633136 : for (size_t i = 0; i < nelem; ++i) {
90 5624640 : if (p_wsum[i] > 0.0) {
91 5624640 : p_tdata[i] /= T(p_wsum[i]);
92 : }
93 : }
94 :
95 8496 : dataOut.reference(transformedData);
96 8496 : }
97 :
98 5624640 : static inline void AccumulateWeight(Float const wt, Double &wtsum) {
99 5624640 : wtsum += 1.0 / wt;
100 5624640 : }
101 :
102 2812320 : static inline void NormalizeWeight(Double const wtsum, Float &wt) {
103 2812320 : wt = 4.0 / wtsum;
104 2812320 : }
105 : };
106 :
107 : struct GeometricTransformation {
108 : template<class T>
109 1188 : inline static void transformData(Cube<T> const &dataIn,
110 : Cube<Bool> const &flagIn, Cube<Float> const &weightIn, Int pid0, Int pid1,
111 : Cube<T> &dataOut) {
112 1188 : if (dataIn.empty()) {
113 0 : dataOut.resize();
114 0 : return;
115 : }
116 1188 : auto const cubeShape = dataIn.shape();
117 1188 : IPosition const newShape(3, 1, cubeShape[1], cubeShape[2]);
118 1188 : Cube<T> transformedData(newShape, T(0));
119 1188 : Cube<Bool> transformedFlag(newShape, True);
120 1188 : Cube<Float> weightSum(newShape, 0.0f);
121 : // size_t const npol = cubeShape[0];
122 1188 : size_t const nchan = cubeShape[1];
123 1188 : size_t const nrow = cubeShape[2];
124 1188 : size_t const nelem = cubeShape[1] * cubeShape[2];
125 :
126 1188 : Int pols[] = { pid0, pid1 };
127 3564 : for (size_t i = 0; i < 2; ++i) {
128 2376 : Int ipol = pols[i];
129 2376 : IPosition start(3, ipol, 0, 0);
130 2376 : IPosition end(3, ipol, nchan - 1, nrow - 1);
131 2376 : auto const dslice = dataIn(start, end);
132 2376 : auto const wslice = weightIn(start, end);
133 2376 : auto const fslice = flagIn(start, end);
134 2376 : Array<Float> weight(dslice.shape());
135 2376 : Array<T> weightedData(dslice.shape());
136 2376 : arrayContTransform(dslice * wslice, fslice, weightedData,
137 : replaceFlaggedDataWithZero<T>);
138 2376 : arrayContTransform(wslice, fslice, weight,
139 : replaceFlaggedDataWithZero<Float>);
140 2376 : transformedData += weightedData;
141 2376 : weightSum += weight;
142 : }
143 :
144 : // transformedData, transformedFlag and nAccumulated should be contiguous array
145 1188 : T *p_tdata = transformedData.data();
146 1188 : Float *p_wsum = weightSum.data();
147 2994948 : for (size_t i = 0; i < nelem; ++i) {
148 2993760 : if (p_wsum[i] > 0.0) {
149 2993760 : p_tdata[i] /= T(p_wsum[i]);
150 : }
151 : }
152 :
153 1188 : dataOut.reference(transformedData);
154 1188 : }
155 :
156 : template<class T>
157 0 : inline static void transformData(Cube<T> const &dataIn,
158 : Cube<Bool> const &flagIn, Matrix<Float> const &weightIn, Int pid0,
159 : Int pid1, Cube<T> &dataOut) {
160 : // cout << "start " << __func__ << endl;
161 0 : if (dataIn.empty()) {
162 0 : dataOut.resize();
163 0 : return;
164 : }
165 0 : auto const cubeShape = dataIn.shape();
166 0 : IPosition const newShape(3, 1, cubeShape[1], cubeShape[2]);
167 0 : Cube<T> transformedData(newShape, T(0));
168 0 : Cube<Float> weightSum(newShape, 0.0f);
169 : // auto const npol = dataIn.shape()[0];
170 0 : auto const nchan = dataIn.shape()[1];
171 0 : auto const nrow = dataIn.shape()[2];
172 0 : auto const nelem = nchan * nrow;
173 :
174 0 : Int pols[] = { pid0, pid1 };
175 0 : for (ssize_t i = 0; i < 2; ++i) {
176 0 : Int ipol = pols[i];
177 0 : for (ssize_t j = 0; j < nrow; ++j) {
178 0 : IPosition start(3, ipol, 0, j);
179 0 : IPosition end(3, ipol, nchan - 1, j);
180 0 : auto dslice = dataIn(start, end);
181 0 : auto fslice = flagIn(start, end);
182 0 : auto w = weightIn(ipol, j);
183 0 : Array<Float> weight(dslice.shape());
184 0 : Array<T> weightedData(dslice.shape());
185 0 : arrayContTransform(dslice * w, fslice, weightedData,
186 : replaceFlaggedDataWithZero<T>);
187 0 : arrayContTransform(fslice, weight, [&w](Bool b) {
188 0 : return ((b) ? 0.0f: w);
189 : });
190 0 : IPosition tstart(3, 0, 0, j);
191 0 : IPosition tend(3, 0, nchan - 1, j);
192 0 : Array<T> tdSlice = transformedData(tstart, tend);
193 0 : Array<Float> twSlice = weightSum(tstart, tend);
194 0 : AlwaysAssert(tdSlice.conform(weightedData), AipsError);
195 0 : AlwaysAssert(twSlice.conform(weight), AipsError);
196 0 : tdSlice += weightedData;
197 0 : twSlice += weight;
198 : }
199 : }
200 :
201 : // transformedData, transformedFlag and nAccumulated should be contiguous array
202 0 : T *p_tdata = transformedData.data();
203 0 : Float *p_wsum = weightSum.data();
204 :
205 0 : for (ssize_t i = 0; i < nelem; ++i) {
206 0 : if (p_wsum[i] > 0.0) {
207 0 : p_tdata[i] /= T(p_wsum[i]);
208 : }
209 : }
210 :
211 0 : dataOut.reference(transformedData);
212 : // cout << "end " << __func__ << endl;
213 0 : }
214 :
215 2993760 : static inline void AccumulateWeight(Float const wt, Double &wtsum) {
216 2993760 : wtsum += wt;
217 2993760 : }
218 :
219 1496880 : static inline void NormalizeWeight(Double const wtsum, Float &wt) {
220 1496880 : wt = wtsum;
221 1496880 : }
222 : };
223 :
224 0 : inline Float weight2Sigma(Float x) {
225 0 : return 1.0 / sqrt(x);
226 : }
227 :
228 : template<class WeightHandler>
229 4842 : inline void transformWeight(Array<Float> const &weightIn, Int pid0, Int pid1,
230 : Array<Float> &weightOut) {
231 : // cout << "start " << __func__ << endl;
232 4842 : if (weightIn.empty()) {
233 : // cout << "input weight is empty" << endl;
234 0 : weightOut.resize();
235 0 : return;
236 : }
237 4842 : IPosition const shapeIn = weightIn.shape();
238 4842 : IPosition shapeOut(shapeIn);
239 : // set length of polarization axis to 1
240 4842 : shapeOut[0] = 1;
241 : // cout << "shapeIn = " << shapeIn << " shapeOut = " << shapeOut << endl;
242 :
243 : // initialization
244 4842 : weightOut.resize(shapeOut);
245 4842 : weightOut = 0.0f;
246 :
247 4842 : ssize_t numPol = shapeIn[0];
248 4842 : Int64 numElemPerPol = shapeOut.product();
249 : // cout << "numElemPerPol = " << numElemPerPol << endl;
250 : // cout << "numPol = " << numPol << endl;
251 :
252 : Bool b;
253 4842 : Float const *p_wIn = weightIn.getStorage(b);
254 4842 : Float *p_wOut = weightOut.data();
255 :
256 4842 : Int pols[] = { pid0, pid1 };
257 4314042 : for (Int64 i = 0; i < numElemPerPol; ++i) {
258 4309200 : ssize_t offsetIndex = i * numPol;
259 4309200 : Double sum = 0.0;
260 12927600 : for (ssize_t j = 0; j < 2; ++j) {
261 8618400 : Int ipol = pols[j];
262 8618400 : WeightHandler::AccumulateWeight(p_wIn[offsetIndex + ipol], sum);
263 : }
264 4309200 : WeightHandler::NormalizeWeight(sum, p_wOut[i]);
265 : }
266 :
267 4842 : weightIn.freeStorage(p_wIn, b);
268 4842 : }
269 : } // anonymous namespace
270 :
271 : namespace casa { //# NAMESPACE CASA - BEGIN
272 :
273 : namespace vi { //# NAMESPACE VI - BEGIN
274 : //////////
275 : // Base Class
276 : // PolAverageTVI
277 : /////////
278 5 : PolAverageTVI::PolAverageTVI(ViImplementation2 *inputVII) :
279 5 : TransformingVi2(inputVII) {
280 5 : configurePolAverage();
281 :
282 : // Initialize attached VisBuffer
283 5 : setVisBuffer(createAttachedVisBuffer(VbRekeyable));
284 5 : }
285 :
286 5 : PolAverageTVI::~PolAverageTVI() {
287 5 : }
288 :
289 54 : void PolAverageTVI::origin() {
290 54 : TransformingVi2::origin();
291 :
292 : // Configure the correlations per shape
293 54 : configureShapes();
294 :
295 : // Synchronize own VisBuffer
296 54 : configureNewSubchunk();
297 :
298 : // reconfigure if necessary
299 54 : reconfigurePolAverageIfNecessary();
300 :
301 : // warn if current dd is inappropriate for polarization averaging
302 54 : warnIfNoTransform();
303 54 : }
304 :
305 4842 : void PolAverageTVI::next() {
306 4842 : TransformingVi2::next();
307 :
308 : // Configure the correlations per shape
309 4842 : configureShapes();
310 :
311 : // Synchronize own VisBuffer
312 4842 : configureNewSubchunk();
313 :
314 : // reconfigure if necessary
315 4842 : reconfigurePolAverageIfNecessary();
316 :
317 : // warn if current dd is inappropriate for polarization averaging
318 4842 : warnIfNoTransform();
319 4842 : }
320 :
321 4896 : void PolAverageTVI::configureShapes() {
322 4896 : Vector <Int> corrs = getCorrelations ();
323 4896 : Int nCorrs = corrs.nelements();
324 :
325 4896 : nCorrelationsPerShape_ = casacore::Vector<casacore::Int> (1, nCorrs);
326 4896 : }
327 :
328 4896 : void PolAverageTVI::warnIfNoTransform() {
329 4896 : if (!doTransform_[dataDescriptionId()]) {
330 0 : auto const vb = getVii()->getVisBuffer();
331 0 : LogIO os(LogOrigin("PolAverageTVI", __func__, WHERE));
332 0 : String msg("Skip polarization average because");
333 0 : if (vb->nCorrelations() == 1) {
334 0 : msg += " number of polarizations is 1.";
335 0 : } else if (anyEQ(vb->correlationTypes(), (Int) Stokes::I)) {
336 0 : msg += " polarization type is Stokes.";
337 : } else {
338 0 : msg += " no valid polarization components are found.";
339 : }
340 0 : os << LogIO::WARN << msg << LogIO::POST;
341 0 : }
342 4896 : }
343 :
344 4842 : void PolAverageTVI::corrType(Vector<Int> & corrTypes) const {
345 4842 : if (doTransform_[dataDescriptionId()]) {
346 : // Always return (Stokes::I)
347 4842 : Vector<Int> myCorrTypes(1, (Int) Stokes::I);
348 4842 : corrTypes.reference(myCorrTypes);
349 4842 : } else {
350 0 : getVii()->corrType(corrTypes);
351 : }
352 4842 : }
353 :
354 4842 : void PolAverageTVI::flagRow(Vector<Bool> & rowflags) const {
355 4842 : Cube<Bool> const &flags = getVisBuffer()->flagCube();
356 4842 : accumulateFlagCube(flags, rowflags);
357 4842 : }
358 :
359 4842 : void PolAverageTVI::flag(Cube<Bool> & flags) const {
360 4842 : auto const vb = getVii()->getVisBuffer();
361 4842 : Cube<Bool> originalFlags = vb->flagCube();
362 4842 : Int ddid = dataDescriptionId();
363 :
364 4842 : if (doTransform_[ddid]) {
365 4842 : auto const cubeShape = originalFlags.shape();
366 4842 : IPosition const newShape(3, 1, cubeShape[1], cubeShape[2]);
367 4842 : Cube<Bool> transformedFlags(newShape, True);
368 : // accumulate first polarization component
369 4842 : IPosition start(3, polId0_[ddid], 0, 0);
370 4842 : IPosition end(3, polId0_[ddid], cubeShape[1] - 1, cubeShape[2] - 1);
371 4842 : transformedFlags = originalFlags(start, end);
372 :
373 : // accumulate second polarization component
374 4842 : start[0] = polId1_[ddid];
375 4842 : end[0] = polId1_[ddid];
376 4842 : transformedFlags &= originalFlags(start, end);
377 4842 : flags.reference(transformedFlags);
378 4842 : } else {
379 0 : flags.reference(originalFlags);
380 : }
381 4842 : }
382 :
383 0 : void PolAverageTVI::flag(Matrix<Bool> & flags) const {
384 0 : Cube<Bool> transformedFlags;
385 0 : flag(transformedFlags);
386 :
387 0 : flags.reference(transformedFlags.yzPlane(0));
388 0 : }
389 :
390 0 : void PolAverageTVI::jonesC(Vector<SquareMatrix<Complex, 2> > &cjones) const {
391 0 : if (doTransform_[dataDescriptionId()]) {
392 0 : throw AipsError("PolAverageTVI::jonesC should not be called.");
393 : } else {
394 0 : getVii()->jonesC(cjones);
395 : }
396 0 : }
397 :
398 0 : void PolAverageTVI::sigma(Matrix<Float> & sigmat) const {
399 0 : if (weightSpectrumExists()) {
400 0 : Cube<Float> const &sigmaSp = getVisBuffer()->sigmaSpectrum();
401 0 : Cube<Bool> const &flag = getVisBuffer()->flagCube();
402 0 : accumulateWeightCube(sigmaSp, flag, sigmat);
403 : } else {
404 0 : if (doTransform_[dataDescriptionId()]) {
405 0 : weight(sigmat);
406 0 : arrayTransformInPlace(sigmat, ::weight2Sigma);
407 : } else {
408 0 : getVii()->sigma(sigmat);
409 : }
410 : }
411 0 : }
412 :
413 4842 : void PolAverageTVI::visibilityCorrected(Cube<Complex> & vis) const {
414 4842 : if (getVii()->existsColumn(VisBufferComponent2::VisibilityCorrected)) {
415 4842 : Cube<Complex> dataCube;
416 4842 : getVii()->visibilityCorrected(dataCube);
417 4842 : if (doTransform_[dataDescriptionId()]) {
418 4842 : Cube<Bool> flagCube;
419 4842 : getVii()->flag(flagCube);
420 4842 : transformComplexData(dataCube, flagCube, vis);
421 4842 : } else {
422 0 : vis.reference(dataCube);
423 : }
424 4842 : } else {
425 0 : vis.resize();
426 : }
427 4842 : }
428 :
429 4842 : void PolAverageTVI::visibilityModel(Cube<Complex> & vis) const {
430 4842 : if (getVii()->existsColumn(VisBufferComponent2::VisibilityModel)) {
431 4842 : Cube<Complex> dataCube;
432 4842 : getVii()->visibilityModel(dataCube);
433 4842 : if (doTransform_[dataDescriptionId()]) {
434 4842 : Cube<Bool> flagCube;
435 4842 : getVii()->flag(flagCube);
436 4842 : transformComplexData(dataCube, flagCube, vis);
437 4842 : } else {
438 0 : vis.reference(dataCube);
439 : }
440 4842 : } else {
441 0 : vis.resize();
442 : }
443 4842 : }
444 :
445 0 : void PolAverageTVI::visibilityObserved(Cube<Complex> & vis) const {
446 0 : if (getVii()->existsColumn(VisBufferComponent2::VisibilityObserved)) {
447 0 : Cube<Complex> dataCube;
448 0 : getVii()->visibilityObserved(dataCube);
449 0 : if (doTransform_[dataDescriptionId()]) {
450 0 : Cube<Bool> flagCube;
451 0 : getVii()->flag(flagCube);
452 0 : transformComplexData(dataCube, flagCube, vis);
453 0 : } else {
454 0 : vis.reference(dataCube);
455 : }
456 0 : } else {
457 0 : vis.resize();
458 : }
459 0 : }
460 :
461 0 : void PolAverageTVI::floatData(casacore::Cube<casacore::Float> & fcube) const {
462 0 : if (getVii()->existsColumn(VisBufferComponent2::FloatData)) {
463 0 : Cube<Float> dataCube;
464 0 : getVii()->floatData(dataCube);
465 0 : if (doTransform_[dataDescriptionId()]) {
466 0 : Cube<Bool> flagCube;
467 0 : getVii()->flag(flagCube);
468 0 : transformFloatData(dataCube, flagCube, fcube);
469 0 : } else {
470 0 : fcube.reference(dataCube);
471 : }
472 0 : } else {
473 0 : fcube.resize();
474 : }
475 0 : }
476 :
477 0 : IPosition PolAverageTVI::visibilityShape() const {
478 0 : IPosition cubeShape = getVii()->visibilityShape();
479 0 : if (doTransform_[dataDescriptionId()]) {
480 : // Length of polarization (Stokes) axis is always 1 after polarizaton averaging
481 0 : cubeShape[0] = 1;
482 : }
483 0 : return cubeShape;
484 0 : }
485 :
486 0 : void PolAverageTVI::weight(Matrix<Float> & wtmat) const {
487 0 : if (weightSpectrumExists()) {
488 0 : Cube<Float> const &weightSp = getVisBuffer()->weightSpectrum();
489 0 : Cube<Bool> const &flag = getVisBuffer()->flagCube();
490 0 : accumulateWeightCube(weightSp, flag, wtmat);
491 : } else {
492 0 : Matrix<Float> wtmatOrg;
493 0 : getVii()->weight(wtmatOrg);
494 0 : if (doTransform_[dataDescriptionId()]) {
495 0 : transformWeight(wtmatOrg, wtmat);
496 : } else {
497 0 : wtmat.reference(wtmatOrg);
498 : }
499 0 : }
500 0 : }
501 :
502 4842 : void PolAverageTVI::weightSpectrum(Cube<Float> & wtsp) const {
503 4842 : if (weightSpectrumExists()) {
504 4842 : Cube<Float> wtspOrg;
505 4842 : getVii()->weightSpectrum(wtspOrg);
506 4842 : if (doTransform_[dataDescriptionId()]) {
507 4842 : transformWeight(wtspOrg, wtsp);
508 : } else {
509 0 : wtsp.reference(wtspOrg);
510 : }
511 4842 : } else {
512 0 : wtsp.resize();
513 : }
514 4842 : }
515 :
516 0 : void PolAverageTVI::sigmaSpectrum(Cube<Float> & wtsp) const {
517 0 : if (sigmaSpectrumExists()) {
518 0 : if (doTransform_[dataDescriptionId()]) {
519 : // sigma = (weight)^-1/2
520 0 : weightSpectrum(wtsp);
521 0 : arrayTransformInPlace(wtsp, ::weight2Sigma);
522 : } else {
523 0 : getVii()->sigmaSpectrum(wtsp);
524 : }
525 : } else {
526 0 : wtsp.resize();
527 : }
528 0 : }
529 :
530 0 : const VisImagingWeight & PolAverageTVI::getImagingWeightGenerator() const {
531 0 : if (doTransform_[dataDescriptionId()]) {
532 0 : throw AipsError(
533 0 : "PolAverageTVI::getImagingWeightGenerator should not be called.");
534 : }
535 :
536 0 : return getVii()->getImagingWeightGenerator();
537 : }
538 :
539 9792 : Vector<Int> PolAverageTVI::getCorrelations() const {
540 9792 : if (doTransform_[dataDescriptionId()]) {
541 : // Always return (Stokes::I)
542 19584 : return Vector<Int>(1, Stokes::I);
543 : } else {
544 0 : return getVii()->getCorrelations();
545 : }
546 : }
547 :
548 4896 : Vector<Stokes::StokesTypes> PolAverageTVI::getCorrelationTypesDefined() const {
549 4896 : if (doTransform_[dataDescriptionId()]) {
550 : // Always return (Stokes::I)
551 9792 : return Vector<Stokes::StokesTypes>(1, Stokes::I);
552 : } else {
553 0 : return getVii()->getCorrelationTypesDefined();
554 : }
555 : }
556 :
557 4896 : Vector<Stokes::StokesTypes> PolAverageTVI::getCorrelationTypesSelected() const {
558 4896 : if (doTransform_[dataDescriptionId()]) {
559 : // Always return (Stokes::I)
560 9792 : return Vector<Stokes::StokesTypes>(1, Stokes::I);
561 : } else {
562 0 : return getVii()->getCorrelationTypesSelected();
563 : }
564 : }
565 :
566 : const casacore::Vector<casacore::Int>&
567 4896 : PolAverageTVI::nCorrelationsPerShape() const {
568 4896 : return nCorrelationsPerShape_;
569 : }
570 :
571 493 : void PolAverageTVI::configurePolAverage() {
572 493 : MeasurementSet const &ms = getVii()->ms();
573 493 : auto const &msdd = ms.dataDescription();
574 493 : MSDataDescColumns msddcols(msdd);
575 493 : uInt ndd = msddcols.nrow();
576 493 : Vector<Int> polIds = msddcols.polarizationId().getColumn();
577 493 : doTransform_.resize(ndd);
578 493 : polId0_.resize(ndd);
579 493 : polId1_.resize(ndd);
580 493 : auto const &mspol = ms.polarization();
581 493 : MSPolarizationColumns mspolcols(mspol);
582 493 : doTransform_ = False;
583 2465 : for (uInt idd = 0; idd < ndd; ++idd) {
584 1972 : Vector<Int> corrType = mspolcols.corrType()(polIds[idd]);
585 1972 : polId0_[idd] = -1;
586 1972 : polId1_[idd] = -1;
587 9860 : for (size_t i = 0; i < corrType.size(); ++i) {
588 7888 : auto stokesType = Stokes::type(corrType[i]);
589 7888 : if (stokesType == Stokes::XX || stokesType == Stokes::RR) {
590 1972 : polId0_[idd] = i;
591 5916 : } else if (stokesType == Stokes::YY || stokesType == Stokes::LL) {
592 1972 : polId1_[idd] = i;
593 : }
594 : }
595 1972 : doTransform_[idd] = (polId0_[idd] >= 0 && polId1_[idd] >= 0);
596 1972 : }
597 493 : }
598 :
599 : //////////
600 : // GeometricPolAverageTVI
601 : /////////
602 2 : GeometricPolAverageTVI::GeometricPolAverageTVI(ViImplementation2 *inputVII) :
603 2 : PolAverageTVI(inputVII) {
604 2 : }
605 :
606 4 : GeometricPolAverageTVI::~GeometricPolAverageTVI() {
607 4 : }
608 :
609 1188 : void GeometricPolAverageTVI::transformComplexData(Cube<Complex> const &dataIn,
610 : Cube<Bool> const &flagIn, Cube<Complex> &dataOut) const {
611 1188 : Int ddid = dataDescriptionId();
612 1188 : Int pid0 = polId0_[ddid];
613 1188 : Int pid1 = polId1_[ddid];
614 1188 : transformData(dataIn, flagIn, pid0, pid1, dataOut);
615 1188 : }
616 :
617 0 : void GeometricPolAverageTVI::transformFloatData(Cube<Float> const &dataIn,
618 : Cube<Bool> const &flagIn, Cube<Float> &dataOut) const {
619 0 : Int ddid = dataDescriptionId();
620 0 : Int pid0 = polId0_[ddid];
621 0 : Int pid1 = polId1_[ddid];
622 0 : transformData(dataIn, flagIn, pid0, pid1, dataOut);
623 0 : }
624 :
625 594 : void GeometricPolAverageTVI::transformWeight(Array<Float> const &weightIn,
626 : Array<Float> &weightOut) const {
627 594 : Int ddid = dataDescriptionId();
628 594 : Int pid0 = polId0_[ddid];
629 594 : Int pid1 = polId1_[ddid];
630 594 : ::transformWeight<GeometricTransformation>(weightIn, pid0, pid1, weightOut);
631 594 : }
632 :
633 : template<class T>
634 1188 : void GeometricPolAverageTVI::transformData(Cube<T> const &dataIn,
635 : Cube<Bool> const &flagIn, Int pid0, Int pid1, Cube<T> &dataOut) const {
636 1188 : if (weightSpectrumExists()) {
637 1188 : Cube<Float> weightSp;
638 1188 : getVii()->weightSpectrum(weightSp);
639 1188 : ::GeometricTransformation::transformData<T>(dataIn, flagIn, weightSp, pid0,
640 : pid1, dataOut);
641 1188 : } else {
642 0 : Matrix<Float> weightMat;
643 0 : getVii()->weight(weightMat);
644 0 : ::GeometricTransformation::transformData<T>(dataIn, flagIn, weightMat, pid0,
645 : pid1, dataOut);
646 0 : }
647 1188 : }
648 :
649 : //////////
650 : // StokesPolAverageTVI
651 : /////////
652 3 : StokesPolAverageTVI::StokesPolAverageTVI(ViImplementation2 *inputVII) :
653 3 : PolAverageTVI(inputVII) {
654 3 : }
655 :
656 6 : StokesPolAverageTVI::~StokesPolAverageTVI() {
657 6 : }
658 :
659 8496 : void StokesPolAverageTVI::transformComplexData(Cube<Complex> const &dataIn,
660 : Cube<Bool> const &flagIn, Cube<Complex> &dataOut) const {
661 8496 : Int ddid = dataDescriptionId();
662 8496 : Int pid0 = polId0_[ddid];
663 8496 : Int pid1 = polId1_[ddid];
664 8496 : transformData(dataIn, flagIn, pid0, pid1, dataOut);
665 8496 : }
666 :
667 0 : void StokesPolAverageTVI::transformFloatData(Cube<Float> const &dataIn,
668 : Cube<Bool> const &flagIn, Cube<Float> &dataOut) const {
669 0 : Int ddid = dataDescriptionId();
670 0 : Int pid0 = polId0_[ddid];
671 0 : Int pid1 = polId1_[ddid];
672 0 : transformData(dataIn, flagIn, pid0, pid1, dataOut);
673 0 : }
674 :
675 4248 : void StokesPolAverageTVI::transformWeight(Array<Float> const &weightIn,
676 : Array<Float> &weightOut) const {
677 4248 : Int ddid = dataDescriptionId();
678 4248 : Int pid0 = polId0_[ddid];
679 4248 : Int pid1 = polId1_[ddid];
680 4248 : ::transformWeight<StokesTransformation>(weightIn, pid0, pid1, weightOut);
681 4248 : }
682 :
683 : template<class T>
684 8496 : void StokesPolAverageTVI::transformData(Cube<T> const &dataIn,
685 : Cube<Bool> const &flagIn, Int pid0, Int pid1, Cube<T> &dataOut) const {
686 8496 : ::StokesTransformation::transformData<T>(dataIn, flagIn, pid0, pid1, dataOut);
687 8496 : }
688 :
689 : //////////
690 : // PolAverageTVIFactory
691 : /////////
692 5 : PolAverageVi2Factory::PolAverageVi2Factory(Record const &configuration,
693 5 : ViImplementation2 *inputVII) :
694 5 : inputVII_p(inputVII), mode_(AveragingMode::DEFAULT) {
695 5 : inputVII_p = inputVII;
696 :
697 5 : mode_ = PolAverageVi2Factory::GetAverageModeFromConfig(configuration);
698 5 : }
699 :
700 0 : PolAverageVi2Factory::PolAverageVi2Factory(Record const &configuration,
701 : MeasurementSet const *ms, SortColumns const sortColumns,
702 0 : Double timeInterval, Bool isWritable) :
703 0 : inputVII_p(nullptr), mode_(AveragingMode::DEFAULT) {
704 0 : inputVII_p = new VisibilityIteratorImpl2(Block<MeasurementSet const *>(1, ms),
705 0 : sortColumns, timeInterval, isWritable);
706 :
707 0 : mode_ = PolAverageVi2Factory::GetAverageModeFromConfig(configuration);
708 0 : }
709 :
710 5 : PolAverageVi2Factory::~PolAverageVi2Factory() {
711 5 : }
712 :
713 5 : ViImplementation2 * PolAverageVi2Factory::createVi() const {
714 10 : LogIO os(LogOrigin("PolAverageVi2Factory", __func__, WHERE));
715 5 : if (mode_ == AveragingMode::GEOMETRIC) {
716 2 : String msg("Combining parallel-hand correlations directly.");
717 2 : os << LogIO::NORMAL << msg << LogIO::POST;
718 2 : return new GeometricPolAverageTVI(inputVII_p);
719 5 : } else if (mode_ == AveragingMode::STOKES) {
720 3 : String msg("Combining parallel-hand correlations to form Stokes I.");
721 3 : os << LogIO::NORMAL << msg << LogIO::POST;
722 3 : return new StokesPolAverageTVI(inputVII_p);
723 3 : }
724 :
725 0 : throw AipsError("Invalid Averaging Mode for PolAverageTVI.");
726 :
727 : return nullptr;
728 5 : }
729 :
730 5 : PolAverageTVILayerFactory::PolAverageTVILayerFactory(
731 5 : Record const &configuration) :
732 5 : ViiLayerFactory() {
733 5 : configuration_p = configuration;
734 5 : }
735 :
736 : ViImplementation2*
737 5 : PolAverageTVILayerFactory::createInstance(ViImplementation2* vii0) const {
738 : // Make the PolAverageTVI, using supplied ViImplementation2, and return it
739 5 : PolAverageVi2Factory factory(configuration_p, vii0);
740 5 : ViImplementation2 *vii = nullptr;
741 : try {
742 5 : vii = factory.createVi();
743 0 : } catch (...) {
744 0 : if (vii0) {
745 0 : delete vii0;
746 : }
747 0 : throw;
748 0 : }
749 5 : return vii;
750 5 : }
751 : } // # NAMESPACE VI - END
752 : } // #NAMESPACE CASA - END
|
