Line data Source code
1 : //# SPWCombinationTVI.h: This file contains the implementation of the SPWCombinationTVI 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 :
23 : #include <iomanip>
24 : #include <mstransform/TVI/SPWCombinationTVI.h>
25 : #include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
26 :
27 : namespace casa { //# NAMESPACE CASA - BEGIN
28 :
29 : namespace vi { //# NAMESPACE VI - BEGIN
30 :
31 : //////////////////////////////////////////////////////////////////////////
32 : // SPWCombinationTVI class
33 : //////////////////////////////////////////////////////////////////////////
34 :
35 : // -----------------------------------------------------------------------
36 : //
37 : // -----------------------------------------------------------------------
38 0 : SPWCombinationTVI::SPWCombinationTVI(ViImplementation2 * inputVii) :
39 0 : FreqAxisTVI (inputVii)
40 : {
41 0 : std::cout << std::setprecision(10) << std::endl;
42 0 : initialize();
43 :
44 0 : return;
45 0 : }
46 :
47 : // -----------------------------------------------------------------------
48 : //
49 : // -----------------------------------------------------------------------
50 0 : void SPWCombinationTVI::initialize()
51 : {
52 0 : freqWidthChan_p = checkEqualWidth();
53 :
54 0 : spwInpChanOutFreqMap_p.clear();
55 :
56 : // Where do the output SPWs indexes start. This assumes that
57 : // this TVI is not doing any reindexing in the output it produces.
58 0 : outSpwStartIdx_p = inputVii_p->nSpectralWindows();
59 :
60 : // TODO: If the input has not been "reindexed" there can be
61 : // spwIds, polIds, ddIds which are not referenced in the main table
62 : // The proper way to do it is to filter all the input SPWs to those
63 : // ones for which there is real data.
64 0 : auto nPolIds = inputVii_p->nPolarizationIds();
65 : //There are as many output SPWs as polarizations.
66 0 : std::cout << " npol " << nPolIds << std::endl;
67 0 : for(ssize_t polId = 0; polId < nPolIds; polId++)
68 : {
69 0 : int outSpwId = polId + outSpwStartIdx_p;
70 0 : std::vector<double> thisOutSpwFreqs;
71 0 : std::vector<int> thisOutSpwChann;
72 0 : spwInpChanOutFreqMap_p[outSpwId].clear();
73 : // TODO: Filter SPWs that do not appear together with this polarization
74 : // (they are not together in the same DDId)
75 0 : for(auto inpSpw : spwInpChanIdxMap_p)
76 : {
77 0 : for(auto channel : inpSpw.second)
78 : {
79 0 : double channFreq = getChannelNominalFreq(inpSpw.first, channel);
80 0 : spwInpChanOutFreqMap_p[outSpwId][inpSpw.first].push_back(channFreq);
81 0 : thisOutSpwFreqs.push_back(channFreq);
82 : }
83 0 : }
84 0 : spwOutFirstFreq_p[outSpwId] = *std::min_element(thisOutSpwFreqs.begin(), thisOutSpwFreqs.end());
85 0 : for(auto inpSpw : spwInpChanIdxMap_p)
86 : {
87 0 : for(auto freq : spwInpChanOutFreqMap_p[outSpwId][inpSpw.first])
88 : {
89 : //TODO: Check all the frequencies fall more or less in the same bins. With a bin tolerance (0.01% of the bin?)
90 : //TODO: Check that no frequencies overlap.
91 0 : spwInpChanOutMap_p[outSpwId][inpSpw.first].push_back(std::floor((freq - spwOutFirstFreq_p[outSpwId]) / freqWidthChan_p + 0.1));
92 0 : thisOutSpwChann.push_back(std::floor((freq - spwOutFirstFreq_p[outSpwId]) / freqWidthChan_p + 0.1));
93 : }
94 0 : }
95 0 : spwOutChanNumMap_p[outSpwId] = *std::max_element(thisOutSpwChann.begin(), thisOutSpwChann.end()) + 1;
96 0 : }
97 :
98 :
99 0 : return;
100 : }
101 :
102 0 : double SPWCombinationTVI::getChannelNominalFreq(size_t inpSpw, size_t channel) const
103 : {
104 : //TODO: refactor with next one
105 0 : auto& inputSPWSubtablecols = inputVii_p->spectralWindowSubtablecols();
106 :
107 0 : auto& channelNominalFreqs = inputSPWSubtablecols.chanFreq();
108 : //std::cout << " channelNominalFreqs " << channelNominalFreqs.getColumn() << std::endl;
109 :
110 0 : return channelNominalFreqs(inpSpw)(casacore::IPosition(1, channel));
111 : }
112 :
113 0 : double SPWCombinationTVI::checkEqualWidth() const
114 : {
115 0 : double channelWidth = -1;
116 0 : auto& inputSPWSubtablecols = inputVii_p->spectralWindowSubtablecols();
117 :
118 0 : auto& channelWidths = inputSPWSubtablecols.chanWidth();
119 :
120 0 : for(auto& inpSpw : spwInpChanIdxMap_p)
121 : {
122 0 : if (channelWidth == -1)
123 0 : channelWidth = channelWidths(inpSpw.first)(casacore::IPosition(1, 0));
124 0 : auto spwChannelWidths = channelWidths(inpSpw.first); //This is a copy
125 0 : if(!std::all_of(spwChannelWidths.begin(), spwChannelWidths.end(), [&](double width){return width == channelWidth;}))
126 0 : throw casacore::AipsError("SPWs to combine have different widths");
127 0 : }
128 0 : return channelWidth;
129 : }
130 :
131 0 : casacore::Vector<double> SPWCombinationTVI::getFrequencies(double time, int frameOfReference,
132 : int spectralWindowId, int msId) const
133 : {
134 0 : if(spectralWindowId < outSpwStartIdx_p ||
135 0 : spectralWindowId > int(outSpwStartIdx_p + spwInpChanOutFreqMap_p.size()))
136 0 : throw casacore::AipsError("SPWId out of valid range");
137 :
138 0 : auto inputSpwForThisOutputSpw = spwInpChanOutMap_p.at(spectralWindowId);
139 0 : std::vector<double> freqs(spwOutChanNumMap_p.at(spectralWindowId));
140 0 : for(auto inputSpw : inputSpwForThisOutputSpw)
141 : {
142 0 : auto innerFreqs = getVii()->getFrequencies(time, frameOfReference, inputSpw.first, msId);
143 0 : auto outChanThisInputSpw = spwInpChanOutMap_p.at(spectralWindowId).at(inputSpw.first)[0];
144 0 : std::copy(innerFreqs.begin(), innerFreqs.end(), freqs.begin() + outChanThisInputSpw);
145 0 : }
146 0 : return casacore::Vector<double>(freqs);
147 0 : }
148 :
149 0 : void SPWCombinationTVI::origin()
150 : {
151 : // Drive underlying ViImplementation2
152 0 : getVii()->origin();
153 :
154 : // Set structure parameters for this subchunk iteration
155 0 : setUpCurrentSubchunkStructure();
156 :
157 : // Synchronize own VisBuffer
158 0 : configureNewSubchunk();
159 :
160 0 : return;
161 : }
162 :
163 0 : void SPWCombinationTVI::next()
164 : {
165 : // Drive underlying ViImplementation2
166 0 : getVii()->next();
167 :
168 : // Set structure parameters for this subchunk iteration
169 0 : setUpCurrentSubchunkStructure();
170 :
171 : // Synchronize own VisBuffer
172 0 : configureNewSubchunk();
173 :
174 0 : return;
175 : }
176 :
177 0 : void SPWCombinationTVI::setUpCurrentSubchunkStructure()
178 : {
179 0 : auto& innerNRowsPerShape = getVii()->nRowsPerShape();
180 0 : getVii()->polarizationIds(currentSubchunkInnerPolIds_p);
181 0 : std::set<casacore::Int> uniquePolIDs;
182 0 : std::copy(currentSubchunkInnerPolIds_p.begin(),currentSubchunkInnerPolIds_p.end(),
183 : std::inserter(uniquePolIDs, uniquePolIDs.end()));
184 :
185 0 : ssize_t thisSubchunkNPolIds = uniquePolIDs.size();
186 : // This VisBuffer contains one single timestamp with "all" DDIds.
187 : // After SPW combination then number of rows is equal to the number of distinct DDIds,
188 : // i.e., the number of polarizations, which is also the number of distinct shapes
189 0 : nRowsPerShape_p.resize(thisSubchunkNPolIds);
190 0 : nChannelsPerShape_p.resize(thisSubchunkNPolIds);
191 :
192 :
193 0 : currentSubchunkSpwIds_p.resize(thisSubchunkNPolIds);
194 0 : std::iota(currentSubchunkSpwIds_p.begin(), currentSubchunkSpwIds_p.end(), outSpwStartIdx_p);
195 :
196 0 : getVii()->spectralWindows(currentSubchunkInnerSpwIds_p);
197 :
198 : // Set up the channels for shape
199 0 : size_t iShape=0;
200 0 : for(auto outSpw : spwOutChanNumMap_p)
201 : {
202 0 : nChannelsPerShape_p[iShape] = outSpw.second;
203 0 : iShape++;
204 : }
205 :
206 0 : nRowsPerShape_p = thisSubchunkNPolIds;
207 :
208 : //TODO: Check all the NRows are the same.
209 : //TODO: For several polarizations (nShapes)
210 0 : }
211 :
212 0 : casacore::rownr_t SPWCombinationTVI::nShapes() const
213 : {
214 0 : return nRowsPerShape_p.size();
215 : }
216 :
217 0 : const casacore::Vector<casacore::rownr_t>& SPWCombinationTVI::nRowsPerShape() const
218 : {
219 0 : return nRowsPerShape_p;
220 : }
221 :
222 0 : const casacore::Vector<casacore::Int>& SPWCombinationTVI::nChannelsPerShape() const
223 : {
224 0 : return nChannelsPerShape_p;
225 : }
226 :
227 0 : void SPWCombinationTVI::spectralWindows(casacore::Vector<casacore::Int>& spws) const
228 : {
229 0 : spws = currentSubchunkSpwIds_p;
230 0 : }
231 :
232 0 : void SPWCombinationTVI::antenna1(casacore::Vector<casacore::Int> & ant1) const
233 : {
234 : // Resize vector
235 0 : ant1.resize(nShapes());
236 :
237 : // It is assumed the values are constant in this VisBuffer
238 0 : ant1 = getVii()->getVisBuffer()->antenna1()(0);
239 0 : }
240 :
241 0 : void SPWCombinationTVI::antenna2(casacore::Vector<casacore::Int> & ant2) const
242 : {
243 : // Resize vector
244 0 : ant2.resize(nShapes());
245 :
246 : // It is assumed the values are constant in this VisBuffer
247 0 : ant2 = getVii()->getVisBuffer()->antenna2()(0);
248 0 : }
249 :
250 0 : void SPWCombinationTVI::arrayIds (casacore::Vector<casacore::Int>&arrayIds) const
251 : {
252 : // Resize vector
253 0 : arrayIds.resize(nShapes());
254 :
255 : // It is assumed the values are constant in this VisBuffer
256 0 : arrayIds = getVii()->getVisBuffer()->arrayId()(0);
257 0 : }
258 :
259 0 : void SPWCombinationTVI::exposure(casacore::Vector<double> & expo) const
260 : {
261 : // Resize vector
262 0 : expo.resize(nShapes());
263 :
264 : // It is assumed the values are constant in this VisBuffer
265 0 : expo = getVii()->getVisBuffer()->exposure()(0);
266 0 : }
267 :
268 0 : void SPWCombinationTVI::feed1(casacore::Vector<casacore::Int> & fd1) const
269 : {
270 : // Resize vector
271 0 : fd1.resize(nShapes());
272 :
273 : // It is assumed the values are constant in this VisBuffer
274 0 : fd1 = getVii()->getVisBuffer()->feed1()(0);
275 0 : }
276 :
277 0 : void SPWCombinationTVI::feed2(casacore::Vector<casacore::Int> & fd2) const
278 : {
279 : // Resize vector
280 0 : fd2.resize(nShapes());
281 :
282 : // It is assumed the values are constant in this VisBuffer
283 0 : fd2 = getVii()->getVisBuffer()->feed2()(0);
284 0 : }
285 :
286 0 : void SPWCombinationTVI::fieldIds (casacore::Vector<casacore::Int>& fieldId) const
287 : {
288 : // Resize vector
289 0 : fieldId.resize(nShapes());
290 :
291 : // It is assumed the values are constant in this VisBuffer
292 0 : fieldId = getVii()->getVisBuffer()->fieldId()(0);
293 0 : }
294 :
295 0 : void SPWCombinationTVI::observationId (casacore::Vector<casacore::Int> & obsids) const
296 : {
297 : // Resize vector
298 0 : obsids.resize(nShapes());
299 :
300 : // It is assumed the values are constant in this VisBuffer
301 0 : obsids = getVii()->getVisBuffer()->observationId()(0);
302 0 : }
303 :
304 0 : void SPWCombinationTVI::processorId (casacore::Vector<casacore::Int> & procids) const
305 : {
306 : // Resize vector
307 0 : procids.resize(nShapes());
308 :
309 : // It is assumed the values are constant in this VisBuffer
310 0 : procids = getVii()->getVisBuffer()->processorId()(0);
311 0 : }
312 :
313 0 : void SPWCombinationTVI::stateId (casacore::Vector<casacore::Int> & stateids) const
314 : {
315 : // Resize vector
316 0 : stateids.resize(nShapes());
317 :
318 : // It is assumed the values are constant in this VisBuffer
319 0 : stateids = getVii()->getVisBuffer()->stateId()(0);
320 0 : }
321 :
322 0 : void SPWCombinationTVI::time(casacore::Vector<double> & t) const
323 : {
324 : // Resize vector
325 0 : t.resize(nShapes());
326 :
327 : // It is assumed the values are constant in this VisBuffer
328 0 : t = getVii()->getVisBuffer()->time()(0);
329 0 : }
330 :
331 0 : void SPWCombinationTVI::timeCentroid(casacore::Vector<double> & t) const
332 : {
333 : // Resize vector
334 0 : t.resize(nShapes());
335 :
336 : // It is assumed the values are constant in this VisBuffer
337 0 : t = getVii()->getVisBuffer()->timeCentroid()(0);
338 0 : }
339 :
340 0 : void SPWCombinationTVI::timeInterval(casacore::Vector<double> & t) const
341 : {
342 : // Resize vector
343 0 : t.resize(nShapes());
344 :
345 : // It is assumed the values are constant in this VisBuffer
346 0 : t = getVii()->getVisBuffer()->timeInterval()(0);
347 0 : }
348 :
349 0 : void SPWCombinationTVI::flag(casacore::Cube<casacore::Bool>& flag) const
350 : {
351 0 : auto& flagCubes = getVisBuffer()->flagCubes();
352 0 : flag = flagCubes[0];
353 0 : }
354 :
355 0 : void SPWCombinationTVI::flag(casacore::Vector<casacore::Cube<casacore::Bool>>& flagCubes) const
356 : {
357 : // Get input VisBuffer and flag cubes
358 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
359 0 : auto& innerFlagCubes = vb->flagCubes();
360 0 : transformCubesVector(innerFlagCubes, flagCubes);
361 0 : }
362 :
363 0 : void SPWCombinationTVI::floatData(casacore::Cube<casacore::Float> & vis) const
364 : {
365 0 : auto& visCubes = getVisBuffer()->visCubesFloat();
366 0 : vis = visCubes[0];
367 0 : }
368 :
369 0 : void SPWCombinationTVI::floatData(casacore::Vector<casacore::Cube<casacore::Float>> & vis) const
370 : {
371 : // Get input VisBuffer and visibility float cubes
372 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
373 0 : auto& innerVisCubes = vb->visCubesFloat();
374 0 : transformCubesVector(innerVisCubes, vis);
375 0 : }
376 :
377 0 : void SPWCombinationTVI::visibilityObserved(casacore::Cube<casacore::Complex> & vis) const
378 : {
379 0 : auto& visCubes = getVisBuffer()->visCubes();
380 0 : vis = visCubes[0];
381 0 : }
382 :
383 0 : void SPWCombinationTVI::visibilityObserved(casacore::Vector<casacore::Cube<casacore::Complex>> & vis) const
384 : {
385 : // Get input VisBuffer and visibility observed cubes
386 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
387 0 : auto& innerVisCubes = vb->visCubes();
388 0 : transformCubesVector(innerVisCubes, vis);
389 0 : }
390 :
391 0 : void SPWCombinationTVI::visibilityCorrected(casacore::Cube<casacore::Complex> & vis) const
392 : {
393 0 : auto& visCubes = getVisBuffer()->visCubesCorrected();
394 0 : vis = visCubes[0];
395 0 : }
396 :
397 0 : void SPWCombinationTVI::visibilityCorrected(casacore::Vector<casacore::Cube<casacore::Complex>> & vis) const
398 : {
399 : // Get input VisBuffer and visibility corrected cubes
400 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
401 0 : auto& innerVisCubes = vb->visCubesCorrected();
402 0 : transformCubesVector(innerVisCubes, vis);
403 0 : }
404 :
405 0 : void SPWCombinationTVI::visibilityModel(casacore::Cube<casacore::Complex> & vis) const
406 : {
407 0 : auto& visCubes = getVisBuffer()->visCubesModel();
408 0 : vis = visCubes[0];
409 0 : }
410 :
411 0 : void SPWCombinationTVI::visibilityModel(casacore::Vector<casacore::Cube<casacore::Complex>> & vis) const
412 : {
413 : // Get input VisBuffer and visibility model cubes
414 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
415 0 : auto& innerVisCubes = vb->visCubesModel();
416 0 : transformCubesVector(innerVisCubes, vis);
417 0 : }
418 :
419 0 : void SPWCombinationTVI::weight(casacore::Matrix<casacore::Float> &weight) const
420 : {
421 0 : auto& weightCubes = getVisBuffer()->weights();
422 0 : weight = weightCubes[0];
423 0 : }
424 :
425 0 : void SPWCombinationTVI::weight(casacore::Vector<casacore::Matrix<casacore::Float>> &weight) const
426 : {
427 : // Get input VisBuffer and weight cubes
428 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
429 0 : auto& innerWeight = vb->weights();
430 0 : transformMatricesVector(innerWeight, weight);
431 0 : }
432 :
433 0 : void SPWCombinationTVI::weightSpectrum(casacore::Cube<casacore::Float> &weightSp) const
434 : {
435 0 : auto& weightSpCubes = getVisBuffer()->weightSpectra();
436 0 : weightSp = weightSpCubes[0];
437 0 : }
438 :
439 0 : void SPWCombinationTVI::weightSpectrum(casacore::Vector<casacore::Cube<casacore::Float>> &weightSp) const
440 : {
441 : // Get input VisBuffer and weight Sp cubes
442 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
443 0 : auto& innerWeightSp = vb->weightSpectra();
444 0 : transformCubesVector(innerWeightSp, weightSp);
445 0 : }
446 :
447 0 : void SPWCombinationTVI::sigma(casacore::Matrix<casacore::Float> &sigma) const
448 : {
449 0 : auto& sigmaCubes = getVisBuffer()->sigmas();
450 0 : sigma = sigmaCubes[0];
451 0 : }
452 :
453 0 : void SPWCombinationTVI::sigma(casacore::Vector<casacore::Matrix<casacore::Float>> &sigma) const
454 : {
455 : // Get input VisBuffer and sigma cubes
456 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
457 0 : auto& innerSigmaMat = vb->sigmas();
458 0 : transformMatricesVector(innerSigmaMat, sigma);
459 0 : }
460 :
461 0 : void SPWCombinationTVI::sigmaSpectrum(casacore::Cube<casacore::Float> &sigmaSp) const
462 : {
463 0 : auto& sigmaSpCubes = getVisBuffer()->sigmaSpectra();
464 0 : sigmaSp = sigmaSpCubes[0];
465 0 : }
466 :
467 0 : void SPWCombinationTVI::sigmaSpectrum(casacore::Vector<casacore::Cube<casacore::Float>> &sigmaSp) const
468 : {
469 : // Get input VisBuffer and sigma Sp cubes
470 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
471 0 : auto& innerSigmaSp = vb->sigmaSpectra();
472 0 : transformCubesVector(innerSigmaSp, sigmaSp);
473 0 : }
474 :
475 0 : SPWCombinationTVIFactory::SPWCombinationTVIFactory (ViImplementation2 *inputVii)
476 0 : : inputVii_p(inputVii)
477 : {
478 0 : }
479 :
480 0 : SPWCombinationTVIFactory::~SPWCombinationTVIFactory()
481 : {
482 0 : }
483 :
484 0 : ViImplementation2 * SPWCombinationTVIFactory::createVi() const
485 : {
486 0 : ViImplementation2* vii = new SPWCombinationTVI(inputVii_p);
487 0 : return vii;
488 : }
489 :
490 0 : SPWCombinationTVILayerFactory::SPWCombinationTVILayerFactory()
491 : {
492 0 : }
493 :
494 0 : SPWCombinationTVILayerFactory::~SPWCombinationTVILayerFactory()
495 : {
496 0 : }
497 :
498 : ViImplementation2*
499 0 : SPWCombinationTVILayerFactory::createInstance(ViImplementation2* vii0) const
500 : {
501 : // Make the SPWCombinationTVI, using supplied ViImplementation2, and return it
502 0 : ViImplementation2 *vii = new SPWCombinationTVI(vii0);
503 0 : return vii;
504 : }
505 :
506 : } //# NAMESPACE VI - END
507 :
508 : } //# NAMESPACE CASA - END
|
