Line data Source code
1 : //# ChannelAverageTVI.h: This file contains the implementation of the ChannelAverageTVI 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 <mstransform/TVI/ChannelAverageTVI.h>
24 : #include <casacore/casa/Arrays/VectorIter.h>
25 :
26 : #ifdef _OPENMP
27 : #include <omp.h>
28 : #endif
29 :
30 :
31 : using namespace casacore;
32 : namespace casa { //# NAMESPACE CASA - BEGIN
33 :
34 : namespace vi { //# NAMESPACE VI - BEGIN
35 :
36 : //////////////////////////////////////////////////////////////////////////
37 : // ChannelAverageTVI class
38 : //////////////////////////////////////////////////////////////////////////
39 :
40 : // -----------------------------------------------------------------------
41 : //
42 : // -----------------------------------------------------------------------
43 21 : ChannelAverageTVI::ChannelAverageTVI( ViImplementation2 * inputVii,
44 21 : const Record &configuration):
45 21 : FreqAxisTVI (inputVii)
46 : {
47 : // Parse and check configuration parameters
48 : // Note: if a constructor finishes by throwing an exception, the memory
49 : // associated with the object itself is cleaned up — there is no memory leak.
50 21 : if (not parseConfiguration(configuration))
51 0 : throw AipsError("Error parsing ChannelAverageTVI configuration");
52 :
53 21 : if (inputVii == nullptr)
54 0 : throw AipsError("Input Vi is empty");
55 :
56 21 : initialize();
57 :
58 21 : return;
59 0 : }
60 :
61 : // -----------------------------------------------------------------------
62 : //
63 : // -----------------------------------------------------------------------
64 21 : Bool ChannelAverageTVI::parseConfiguration(const Record &configuration)
65 : {
66 21 : int exists = -1;
67 21 : Bool ret = true;
68 :
69 : // Parse chanbin parameter (mandatory)
70 21 : exists = -1;
71 21 : exists = configuration.fieldNumber ("chanbin");
72 21 : if (exists >= 0)
73 : {
74 21 : if ( configuration.type(exists) == casacore::TpInt )
75 : {
76 : Int freqbin;
77 21 : configuration.get (exists, freqbin);
78 21 : chanbin_p = Vector<Int>(spwInpChanIdxMap_p.size(),freqbin);
79 : }
80 0 : else if ( configuration.type(exists) == casacore::TpArrayInt)
81 : {
82 0 : configuration.get (exists, chanbin_p);
83 : }
84 : else
85 : {
86 0 : ret = false;
87 0 : logger_p << LogIO::SEVERE << LogOrigin("ChannelAverageTVI", __FUNCTION__)
88 : << "Wrong format for chanbin parameter (only Int and arrayInt are supported) "
89 0 : << LogIO::POST;
90 : }
91 :
92 42 : logger_p << LogIO::NORMAL << LogOrigin("ChannelAverageTVI", __FUNCTION__)
93 42 : << "Channel bin is " << chanbin_p << LogIO::POST;
94 21 : if (anyEQ(chanbin_p,0)) {
95 0 : logger_p << LogIO::NORMAL << LogOrigin("ChannelAverageTVI", __FUNCTION__)
96 0 : << " NB: Channel bin '0' means no channel averaging for the corresponding spw." << LogIO::POST;
97 : }
98 : }
99 : else
100 : {
101 0 : ret = false;
102 0 : logger_p << LogIO::SEVERE << LogOrigin("ChannelAverageTVI", __FUNCTION__)
103 : << "chanbin parameter not found in configuration "
104 0 : << LogIO::POST;
105 : }
106 :
107 21 : exists = configuration.fieldNumber ("flagdataFlagPropagation");
108 21 : flagdataFlagPropagation_p = exists >= 0;
109 :
110 : // Check consistency between chanbin vector and selected SPW/Chan map
111 21 : if (chanbin_p.size() != spwInpChanIdxMap_p.size())
112 : {
113 0 : ret = false;
114 0 : logger_p << LogIO::SEVERE << LogOrigin("ChannelAverageTVI", __FUNCTION__)
115 : << "Number of elements in chanbin vector does not match number of selected SPWs"
116 0 : << LogIO::POST;
117 : }
118 :
119 21 : return ret;
120 : }
121 :
122 : // -----------------------------------------------------------------------
123 : //
124 : // -----------------------------------------------------------------------
125 21 : void ChannelAverageTVI::initialize()
126 : {
127 : // Populate nchan input-output maps
128 : Int spw;
129 21 : uInt spw_idx = 0;
130 21 : map<Int,vector<Int> >::iterator iter;
131 :
132 89 : for(iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
133 : {
134 68 : spw = iter->first;
135 :
136 : // No averaging when user specifies 0
137 68 : if (chanbin_p(spw_idx)==0) {
138 0 : logger_p << LogIO::DEBUG1 << LogOrigin("ChannelAverageTVI", __FUNCTION__)
139 : << "Specified chanbin for spw " << spw
140 : << " of 0 means no averaging."
141 0 : << LogIO::POST;
142 :
143 0 : spwChanbinMap_p[spw] = 1;
144 : }
145 : // Make sure that chanbin is greater than 1
146 68 : else if ((uInt)chanbin_p(spw_idx) <= 1)
147 : {
148 46 : logger_p << LogIO::DEBUG1 << LogOrigin("ChannelAverageTVI", __FUNCTION__)
149 : << "Specified chanbin for spw " << spw
150 : << " less than 1 falls back to the default number of"
151 23 : << " existing/selected channels: " << iter->second.size()
152 46 : << LogIO::POST;
153 :
154 23 : spwChanbinMap_p[spw] = iter->second.size();
155 : }
156 : // Make sure that chanbin does not exceed number of selected channels
157 45 : else if ((uInt)chanbin_p(spw_idx) > iter->second.size())
158 : {
159 0 : logger_p << LogIO::DEBUG1 << LogOrigin("ChannelAverageTVI", __FUNCTION__)
160 0 : << "Number of selected channels " << iter->second.size()
161 : << " for SPW " << spw
162 0 : << " is smaller than specified chanbin " << (uInt)chanbin_p(spw_idx) << endl
163 0 : << "Setting chanbin to " << iter->second.size()
164 : << " for SPW " << spw
165 0 : << LogIO::POST;
166 0 : spwChanbinMap_p[spw] = iter->second.size();
167 : }
168 : else
169 : {
170 45 : spwChanbinMap_p[spw] = chanbin_p(spw_idx);
171 : }
172 :
173 : // Calculate number of output channels per spw
174 68 : spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size() / spwChanbinMap_p[spw];
175 68 : if (spwInpChanIdxMap_p[spw].size() % spwChanbinMap_p[spw] > 0) spwOutChanNumMap_p[spw] += 1;
176 :
177 68 : spw_idx++;
178 : }
179 :
180 42 : return;
181 : }
182 :
183 : #define DOJUSTO false
184 : // -----------------------------------------------------------------------
185 : //
186 : // -----------------------------------------------------------------------
187 192 : void ChannelAverageTVI::flag(Cube<Bool>& flagCube) const
188 : {
189 :
190 : // Pass-thru for single-channel case
191 192 : if (getVii()->visibilityShape()[1]==1) {
192 2 : getVii()->flag(flagCube);
193 2 : return;
194 : }
195 :
196 : // Get input VisBuffer and SPW
197 190 : VisBuffer2 *vb = getVii()->getVisBuffer();
198 190 : Int inputSPW = vb->spectralWindows()(0);
199 :
200 : // Pass-thru for chanbin=1 case:
201 190 : if (spwChanbinMap_p[inputSPW]==1) {
202 0 : getVii()->flag(flagCube);
203 0 : return;
204 : }
205 :
206 : #ifdef _OPENMP
207 : // Pre-load relevant input info and start clock
208 190 : vb->flagCube();
209 190 : Double time0=omp_get_wtime();
210 : #endif
211 :
212 :
213 : // Reshape output data before passing it to the DataCubeHolder
214 190 : flagCube.resize(getVisBuffer()->getShape(),false);
215 :
216 : // Gather input data
217 190 : DataCubeMap inputData;
218 190 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
219 190 : inputData.add(MS::FLAG,inputFlagCubeHolder);
220 :
221 : // Gather output data
222 190 : DataCubeMap outputData;
223 190 : DataCubeHolder<Bool> outputFlagCubeHolder(flagCube);
224 190 : outputData.add(MS::FLAG,outputFlagCubeHolder);
225 :
226 : // Configure Transformation Engine
227 190 : LogicalANDKernel<Bool> kernel;
228 190 : uInt width = spwChanbinMap_p[inputSPW];
229 190 : ChannelAverageTransformEngine<Bool> transformer(&kernel,&inputData,&outputData,width);
230 :
231 : // Transform data
232 : if (DOJUSTO) {
233 : transformFreqAxis2(vb->getShape(),transformer);
234 : } else {
235 190 : transformer.transformAll();
236 : }
237 :
238 : #ifdef _OPENMP
239 : // Accumulate elapsed time
240 190 : Tfl_+=omp_get_wtime()-time0;
241 : #endif
242 :
243 190 : return;
244 190 : }
245 :
246 : // -----------------------------------------------------------------------
247 : //
248 : // -----------------------------------------------------------------------
249 0 : void ChannelAverageTVI::floatData (Cube<Float> & vis) const
250 : {
251 :
252 : // Pass-thru for single-channel case
253 0 : if (getVii()->visibilityShape()[1]==1) {
254 0 : getVii()->floatData(vis);
255 0 : return;
256 : }
257 :
258 : // Get input VisBuffer and SPW
259 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
260 0 : Int inputSPW = vb->spectralWindows()(0);
261 :
262 : // Pass-thru for chanbin=1 case:
263 0 : if (spwChanbinMap_p[inputSPW]==1) {
264 0 : getVii()->floatData(vis);
265 0 : return;
266 : }
267 :
268 : // Reshape output data before passing it to the DataCubeHolder
269 0 : vis.resize(getVisBuffer()->getShape(),false);
270 :
271 : // Gather input data
272 0 : DataCubeMap inputData;
273 0 : DataCubeHolder<Float> inputVisCubeHolder(vb->visCubeFloat());
274 0 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
275 0 : DataCubeHolder<Float> weightCubeHolder(vb->weightSpectrum());
276 0 : inputData.add(MS::DATA,inputVisCubeHolder);
277 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
278 0 : inputData.add(MS::WEIGHT_SPECTRUM,weightCubeHolder);
279 :
280 : // Gather output data
281 0 : DataCubeMap outputData;
282 0 : DataCubeHolder<Float> outputVisCubeHolder(vis);
283 0 : outputData.add(MS::DATA,outputVisCubeHolder);
284 :
285 : // Configure Transformation Engine
286 0 : uInt width = spwChanbinMap_p[inputSPW];
287 0 : WeightedChannelAverageKernel<Float> kernel;
288 0 : ChannelAverageTransformEngine<Float> transformer(&kernel,&inputData,&outputData,width);
289 :
290 : // Transform data
291 : if (DOJUSTO) {
292 : transformFreqAxis2(vb->getShape(),transformer);
293 : } else {
294 0 : transformer.transformAll();
295 : }
296 :
297 0 : return;
298 0 : }
299 :
300 : // -----------------------------------------------------------------------
301 : //
302 : // -----------------------------------------------------------------------
303 180 : void ChannelAverageTVI::visibilityObserved (Cube<Complex> & vis) const
304 : {
305 :
306 : // Pass-thru for single-channel case
307 180 : if (getVii()->visibilityShape()[1]==1) {
308 2 : getVii()->visibilityObserved(vis);
309 2 : return;
310 : }
311 :
312 : // Get input VisBuffer and SPW
313 178 : VisBuffer2 *vb = getVii()->getVisBuffer();
314 178 : Int inputSPW = vb->spectralWindows()(0);
315 :
316 : // Reshape output data before passing it to the DataCubeHolder
317 178 : vis.resize(getVisBuffer()->getShape(),false);
318 :
319 : // Get weightSpectrum from sigmaSpectrum
320 178 : Cube<Float> weightSpFromSigmaSp;
321 178 : weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),false);
322 178 : weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
323 178 : arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
324 :
325 : // Gather input data
326 178 : DataCubeMap inputData;
327 178 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
328 178 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
329 178 : DataCubeHolder<Float> weightCubeHolder(weightSpFromSigmaSp);
330 178 : inputData.add(MS::DATA,inputVisCubeHolder);
331 178 : inputData.add(MS::FLAG,inputFlagCubeHolder);
332 178 : inputData.add(MS::WEIGHT_SPECTRUM,weightCubeHolder);
333 :
334 : // Gather output data
335 178 : DataCubeMap outputData;
336 178 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
337 178 : outputData.add(MS::DATA,outputVisCubeHolder);
338 :
339 : // Configure Transformation Engine
340 178 : uInt width = spwChanbinMap_p[inputSPW];
341 178 : WeightedChannelAverageKernel<Complex> kernel;
342 178 : ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
343 :
344 : // Transform data
345 : if (DOJUSTO) {
346 : transformFreqAxis2(vb->getShape(),transformer);
347 : }
348 : else {
349 178 : transformer.transformAll();
350 : }
351 :
352 178 : return;
353 178 : }
354 :
355 : // -----------------------------------------------------------------------
356 : //
357 : // -----------------------------------------------------------------------
358 6 : void ChannelAverageTVI::visibilityCorrected (Cube<Complex> & vis) const
359 : {
360 :
361 : // Pass-thru for single-channel case
362 6 : if (getVii()->visibilityShape()[1]==1) {
363 0 : getVii()->visibilityCorrected(vis);
364 0 : return;
365 : }
366 :
367 : // Get input VisBuffer and SPW
368 6 : VisBuffer2 *vb = getVii()->getVisBuffer();
369 6 : Int inputSPW = vb->spectralWindows()(0);
370 :
371 : // Pass-thru for chanbin=1 case:
372 6 : if (spwChanbinMap_p[inputSPW]==1) {
373 0 : getVii()->visibilityCorrected(vis);
374 0 : return;
375 : }
376 :
377 :
378 : #ifdef _OPENMP
379 : // Pre-load relevant input info and start clock
380 6 : vb->visCubeCorrected();
381 6 : vb->flagCube();
382 6 : vb->weightSpectrum();
383 6 : Double time0=omp_get_wtime();
384 : #endif
385 :
386 :
387 : // Reshape output data before passing it to the DataCubeHolder
388 6 : vis.resize(getVisBuffer()->getShape(),false);
389 :
390 : // Gather input data
391 6 : DataCubeMap inputData;
392 6 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
393 6 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
394 6 : DataCubeHolder<Float> weightCubeHolder(vb->weightSpectrum());
395 6 : inputData.add(MS::DATA,inputVisCubeHolder);
396 6 : inputData.add(MS::FLAG,inputFlagCubeHolder);
397 6 : inputData.add(MS::WEIGHT_SPECTRUM,weightCubeHolder);
398 :
399 : // Gather output data
400 6 : DataCubeMap outputData;
401 6 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
402 6 : outputData.add(MS::DATA,outputVisCubeHolder);
403 :
404 : // Configure Transformation Engine
405 6 : uInt width = spwChanbinMap_p[inputSPW];
406 6 : WeightedChannelAverageKernel<Complex> kernel;
407 6 : ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
408 :
409 : // Transform data
410 : if (DOJUSTO) {
411 : transformFreqAxis2(vb->getShape(),transformer);
412 : } else {
413 6 : transformer.transformAll();
414 :
415 : /*
416 : // Demo version upon which upgrades to DataCubeHolder/Map (in UtilsTVI.h)
417 : // and ChannelAverageTransformEngine::transformAll() are based
418 : // As written here, it averages _all_ channels (no
419 : // partial binning is supported).
420 : // NB: This is a bit faster than transformAll
421 : // (due mainly to fewer function calls? E.g. kernel.kernel()? )
422 : vis.set(0.0f); // initialize the output cube
423 : Cube<Complex> ivis(vb->visCubeCorrected());
424 : Cube<Float> iwtsp(vb->weightSpectrum());
425 : Cube<Bool> ifl(vb->flagCube());
426 : VectorIterator<Complex> vi(ivis,1);
427 : VectorIterator<Float> wi(iwtsp,1);
428 : VectorIterator<Bool> fi(ifl,1);
429 : VectorIterator<Complex> vo(vis,1);
430 :
431 : Vector<Complex>& viv = vi.vector();
432 : Vector<Float>& wiv = wi.vector();
433 : Vector<Bool>& fiv = fi.vector();
434 : Vector<Complex>& vov = vo.vector();
435 :
436 : Int nchan=viv.nelements();
437 :
438 : while (!vi.pastEnd()) {
439 :
440 : Float swt(0.0f);
441 : for (Int ich=0;ich<nchan;++ich) {
442 : if (!fiv(ich)) {
443 : vov(0)+=(viv(ich)*wiv(ich));
444 : swt+=wiv(ich);
445 : }
446 : }
447 : if (swt>0.0f)
448 : vov(0)/=swt;
449 : else
450 : vov(0)=0.0;
451 :
452 : vi.next();
453 : wi.next();
454 : fi.next();
455 : vo.next();
456 : }
457 : */
458 : }
459 :
460 : #ifdef _OPENMP
461 : // Accumulate elapsed time
462 6 : Tcd_+=omp_get_wtime()-time0;
463 : #endif
464 :
465 6 : return;
466 6 : }
467 :
468 : // -----------------------------------------------------------------------
469 : //
470 : // -----------------------------------------------------------------------
471 6 : void ChannelAverageTVI::visibilityModel (Cube<Complex> & vis) const
472 : {
473 :
474 : // Pass-thru for single-channel case
475 6 : if (getVii()->visibilityShape()[1]==1) {
476 0 : getVii()->visibilityModel(vis);
477 0 : return;
478 : }
479 :
480 : // Get input VisBuffer and SPW
481 6 : VisBuffer2 *vb = getVii()->getVisBuffer();
482 6 : Int inputSPW = vb->spectralWindows()(0);
483 :
484 : // Pass-thru for chanbin=1 case:
485 6 : if (spwChanbinMap_p[inputSPW]==1) {
486 0 : getVii()->visibilityModel(vis);
487 0 : return;
488 : }
489 :
490 : #ifdef _OPENMP
491 : // Pre-load relevant input info and start clock
492 6 : vb->visCubeModel();
493 6 : vb->flagCube();
494 6 : vb->weightSpectrum();
495 6 : Double time0=omp_get_wtime();
496 : #endif
497 :
498 :
499 : // Reshape output data before passing it to the DataCubeHolder
500 6 : vis.resize(getVisBuffer()->getShape(),false);
501 :
502 : // Gather input data
503 6 : DataCubeMap inputData;
504 6 : DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
505 6 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
506 6 : inputData.add(MS::DATA,inputVisCubeHolder);
507 6 : inputData.add(MS::FLAG,inputFlagCubeHolder);
508 :
509 : // Gather output data
510 6 : DataCubeMap outputData;
511 6 : DataCubeHolder<Complex> outputVisCubeHolder(vis);
512 6 : outputData.add(MS::DATA,outputVisCubeHolder);
513 :
514 : // Configure Transformation Engine
515 6 : uInt width = spwChanbinMap_p[inputSPW];
516 6 : FlaggedChannelAverageKernel<Complex> kernel;
517 6 : ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
518 :
519 : // Transform data
520 : if (DOJUSTO) {
521 : transformFreqAxis2(vb->getShape(),transformer);
522 : } else {
523 6 : transformer.transformAll();
524 : }
525 :
526 : #ifdef _OPENMP
527 : // Accumulate elapsed time
528 6 : Tmd_+=omp_get_wtime()-time0;
529 : #endif
530 :
531 6 : return;
532 6 : }
533 :
534 : // -----------------------------------------------------------------------
535 : //
536 : // -----------------------------------------------------------------------
537 6 : void ChannelAverageTVI::weightSpectrum(Cube<Float> &weightSp) const
538 : {
539 :
540 : // Pass-thru for single-channel case
541 6 : if (getVii()->visibilityShape()[1]==1) {
542 0 : getVii()->weightSpectrum(weightSp);
543 0 : return;
544 : }
545 :
546 : // Get input VisBuffer and SPW
547 6 : VisBuffer2 *vb = getVii()->getVisBuffer();
548 6 : Int inputSPW = vb->spectralWindows()(0);
549 :
550 : // Pass-thru for chanbin=1 case:
551 6 : if (spwChanbinMap_p[inputSPW]==1) {
552 0 : getVii()->weightSpectrum(weightSp);;
553 0 : return;
554 : }
555 :
556 : #ifdef _OPENMP
557 : // Pre-load relevant input info and start clock
558 6 : vb->weightSpectrum();
559 6 : vb->flagCube();
560 6 : Double time0=omp_get_wtime();
561 : #endif
562 :
563 :
564 : // Reshape output data before passing it to the DataCubeHolder
565 6 : weightSp.resize(getVisBuffer()->getShape(),false);
566 :
567 : // Gather input data
568 6 : DataCubeMap inputData;
569 6 : DataCubeHolder<Float> inputWeightCubeHolder(vb->weightSpectrum());
570 6 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
571 6 : inputData.add(MS::DATA,inputWeightCubeHolder);
572 6 : inputData.add(MS::FLAG,inputFlagCubeHolder);
573 :
574 : // Gather output data
575 6 : DataCubeMap outputData;
576 6 : DataCubeHolder<Float> outputWeightCubeHolder(weightSp);
577 6 : outputData.add(MS::DATA,outputWeightCubeHolder);
578 :
579 : // Configure Transformation Engine
580 6 : uInt width = spwChanbinMap_p[inputSPW];
581 6 : ChannelAccumulationKernel<Float> kernel;
582 6 : ChannelAverageTransformEngine<Float> transformer(&kernel,&inputData,&outputData,width);
583 :
584 : // Transform data
585 : if (DOJUSTO) {
586 : transformFreqAxis2(vb->getShape(),transformer);
587 : } else {
588 6 : transformer.transformAll();
589 : }
590 :
591 : #ifdef _OPENMP
592 : // Accumulate elapsed time
593 6 : Tws_+=omp_get_wtime()-time0;
594 : #endif
595 :
596 6 : return;
597 6 : }
598 :
599 : // -----------------------------------------------------------------------
600 : //
601 : // -----------------------------------------------------------------------
602 0 : void ChannelAverageTVI::sigmaSpectrum(Cube<Float> &sigmaSp) const
603 : {
604 : // Pass-thru for single-channel case
605 0 : if (getVii()->visibilityShape()[1]==1) {
606 0 : getVii()->sigmaSpectrum(sigmaSp);
607 0 : return;
608 : }
609 :
610 : // Get input VisBuffer and SPW
611 0 : VisBuffer2 *vb = getVii()->getVisBuffer();
612 0 : Int inputSPW = vb->spectralWindows()(0);
613 :
614 : // Pass-thru for chanbin=1 case:
615 0 : if (spwChanbinMap_p[inputSPW]==1) {
616 0 : getVii()->sigmaSpectrum(sigmaSp);;
617 0 : return;
618 : }
619 :
620 : #ifdef _OPENMP
621 : // Pre-load relevant input info and start clock
622 0 : vb->sigmaSpectrum();
623 0 : vb->flagCube();
624 0 : Double time0=omp_get_wtime();
625 : #endif
626 :
627 : // Reshape output data before passing it to the DataCubeHolder
628 0 : sigmaSp.resize(getVisBuffer()->getShape(),false);
629 :
630 : // Get weightSpectrum from sigmaSpectrum
631 0 : Cube<Float> weightSpFromSigmaSp;
632 0 : weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),false);
633 0 : weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
634 0 : arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
635 :
636 : // Gather input data
637 0 : DataCubeMap inputData;
638 0 : DataCubeHolder<Float> inputWeightCubeHolder(weightSpFromSigmaSp);
639 0 : DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
640 0 : inputData.add(MS::DATA,inputWeightCubeHolder);
641 0 : inputData.add(MS::FLAG,inputFlagCubeHolder);
642 :
643 : // Gather output data
644 0 : DataCubeMap outputData;
645 0 : DataCubeHolder<Float> outputWeightCubeHolder(sigmaSp);
646 0 : outputData.add(MS::DATA,outputWeightCubeHolder);
647 :
648 : // Configure Transformation Engine
649 0 : uInt width = spwChanbinMap_p[inputSPW];
650 0 : ChannelAccumulationKernel<Float> kernel;
651 0 : ChannelAverageTransformEngine<Float> transformer(&kernel,&inputData,&outputData,width);
652 :
653 : // Transform data
654 : if (DOJUSTO) {
655 : transformFreqAxis2(vb->getShape(),transformer);
656 : } else {
657 0 : transformer.transformAll();
658 : }
659 :
660 : // Transform back from weight format to sigma format
661 0 : arrayTransformInPlace (sigmaSp,weightToSigma);
662 :
663 : #ifdef _OPENMP
664 : // Accumulate elapsed time
665 0 : Tss_+=omp_get_wtime()-time0;
666 : #endif
667 :
668 0 : return;
669 0 : }
670 :
671 : // -----------------------------------------------------------------------
672 : //
673 : // -----------------------------------------------------------------------
674 39 : Vector<Double> ChannelAverageTVI::getFrequencies ( Double time,
675 : Int frameOfReference,
676 : Int spectralWindowId,
677 : Int msId) const
678 : {
679 :
680 : // Pass-thru for chanbin=1 case
681 39 : if (spwChanbinMap_p[spectralWindowId]==1) {
682 0 : return getVii()->getFrequencies(time,frameOfReference,spectralWindowId,msId);
683 : }
684 :
685 : // Get frequencies from input VI, which we will process
686 39 : Vector<Double> inputFrequencies = getVii()->getFrequencies(time,frameOfReference,
687 39 : spectralWindowId,msId);
688 :
689 : // Pass-thru for single-channel case
690 : // NB: This does NOT depend on current iteration having same spw as specified spectralWindowId
691 : // NB: Nor does it rely on sensing the apparent shape of the visibility data object (so that need not be pre-filled)
692 39 : if (inputFrequencies.nelements()==1)
693 0 : return inputFrequencies;
694 :
695 : // Produce output (transformed) frequencies
696 39 : Vector<Double> outputFrecuencies(spwOutChanNumMap_p[spectralWindowId],0);
697 :
698 : // Gather input data
699 39 : DataCubeMap inputData;
700 39 : DataCubeHolder<Double> inputFrecuenciesHolder(inputFrequencies);
701 39 : inputData.add(MS::DATA,inputFrecuenciesHolder);
702 :
703 : // Gather output data
704 39 : DataCubeMap outputData;
705 39 : DataCubeHolder<Double> outputFrecuenciesHolder(outputFrecuencies);
706 39 : outputData.add(MS::DATA,outputFrecuenciesHolder);
707 :
708 : // Configure Transformation Engine
709 39 : PlainChannelAverageKernel<Double> kernel;
710 39 : uInt width = spwChanbinMap_p[spectralWindowId];
711 39 : ChannelAverageTransformEngine<Double> transformer(&kernel,&inputData,&outputData,width);
712 :
713 : // Transform data
714 39 : transformer.transform();
715 :
716 39 : return outputFrecuencies;
717 39 : }
718 :
719 0 : Vector<Double> ChannelAverageTVI::getChanWidths ( Double time,
720 : Int frameOfReference,
721 : Int spectralWindowId,
722 : Int msId) const
723 : {
724 :
725 : // Pass-thru for chanbin=1 case
726 0 : if (spwChanbinMap_p[spectralWindowId]==1) {
727 0 : return getVii()->getChanWidths(time,frameOfReference,spectralWindowId,msId);
728 : }
729 :
730 : // Get widths from input VI
731 0 : Vector<Double> inputWidths = getVii()->getChanWidths(time,frameOfReference,
732 0 : spectralWindowId,msId);
733 : // Pass-thru for single-channel case
734 : // NB: This does NOT depend on current iteration having same spw as specified spectralWindowId
735 : // NB: Nor does it rely on sensing the apparent shape of the visibility data object (so that need not be pre-filled)
736 0 : if (inputWidths.nelements()==1)
737 0 : return inputWidths;
738 :
739 : // Produce output (summed) widths
740 0 : Vector<Double> outputWidths(spwOutChanNumMap_p[spectralWindowId],0);
741 :
742 : // Gather input data
743 0 : DataCubeMap inputData;
744 0 : DataCubeHolder<Double> inputWidthsHolder(inputWidths);
745 0 : inputData.add(MS::DATA,inputWidthsHolder);
746 : // We need flag vector (all false==unflagged) for ChannelAccumulationKernal (which just adds within bins)
747 0 : Vector<Bool> nonflagged(inputWidths.nelements(),false);
748 0 : DataCubeHolder<Bool> inputFlagHolder(nonflagged);
749 0 : inputData.add(MS::FLAG,inputFlagHolder);
750 :
751 : // Gather output data
752 0 : DataCubeMap outputData;
753 0 : DataCubeHolder<Double> outputWidthsHolder(outputWidths);
754 0 : outputData.add(MS::DATA,outputWidthsHolder);
755 :
756 : // Configure Transformation Engine
757 : //PlainChannelAverageKernel<Double> kernel;
758 0 : ChannelAccumulationKernel<Double> kernel;
759 0 : uInt width = spwChanbinMap_p[spectralWindowId];
760 0 : ChannelAverageTransformEngine<Double> transformer(&kernel,&inputData,&outputData,width);
761 :
762 : // Transform data
763 0 : transformer.transform();
764 :
765 0 : return outputWidths;
766 0 : }
767 :
768 : // -----------------------------------------------------------------------
769 : //
770 : // -----------------------------------------------------------------------
771 180 : void ChannelAverageTVI::writeFlag (const Cube<Bool> & flag)
772 : {
773 : // Create a flag cube with the input VI shape
774 180 : Cube<Bool> propagatedFlagCube;
775 180 : propagatedFlagCube = getVii()->getVisBuffer()->flagCube();
776 :
777 : // Propagate flags from the input cube to the propagated flag cube
778 180 : propagateChanAvgFlags(flag,propagatedFlagCube);
779 :
780 : // Pass propagated flag cube downstream for further propagation and/or writting
781 180 : getVii()->writeFlag(propagatedFlagCube);
782 :
783 360 : return;
784 180 : }
785 :
786 : /**
787 : * Strategy to support different ways of propagating flags from the 'transformed' cube to
788 : * the original ('propagated') cube. Iterates through rows, channels, correlations.
789 : *
790 : * This is meant to be used from propagateChanAvgFlags with at least two alternative
791 : * functors. One to propagate flags as required by flagdata (preserve pre-existing flags
792 : * in the original data cube), and a second one to propagate flags as required by plotms.
793 : * CAS-12737, CAS-9985, CAS-12205.
794 : *
795 : * @param transformedFlagCube Cube of flags after averaging
796 : * @param propagatedFlagClube Original cube of flags
797 : * @param inputOutputChan input->output channel mapping
798 : * @param propagate functor to implement the (back)propagation of flags (flagdata/plotms).
799 : */
800 : template <typename Functor>
801 180 : void cubePropagateFlags(const Cube<Bool> &transformedFlagCube,
802 : Cube<Bool> &propagatedFlagCube,
803 : const Vector<uInt> &inputOutputChan, Functor propagate)
804 : {
805 : // Get propagated (input) shape
806 180 : const auto inputShape = propagatedFlagCube.shape();
807 180 : const uInt nCorr = inputShape(0);
808 180 : const uInt nChan = inputShape(1);
809 180 : const uInt nRows = inputShape(2);
810 :
811 : // Get transformed (output) shape
812 180 : const auto transformedShape = transformedFlagCube.shape();
813 180 : const auto nTransChan = transformedShape(1);
814 :
815 : uInt outChan;
816 53098 : for (size_t row_i =0;row_i<nRows;row_i++)
817 : {
818 2216109 : for (size_t chan_i =0;chan_i<nChan;chan_i++)
819 : {
820 2163191 : outChan = inputOutputChan(chan_i);
821 2163191 : if (outChan < nTransChan) // outChan >= nChan may happen when channels are dropped
822 : {
823 7462982 : for (size_t corr_i =0;corr_i<nCorr;corr_i++)
824 5299791 : propagate(row_i, chan_i, corr_i, outChan);
825 : }
826 : }
827 : }
828 180 : }
829 :
830 : // -----------------------------------------------------------------------
831 : //
832 : // -----------------------------------------------------------------------
833 180 : void ChannelAverageTVI::propagateChanAvgFlags (const Cube<Bool> &transformedFlagCube,
834 : Cube<Bool> &propagatedFlagCube)
835 : {
836 : // Get current SPW and chanbin
837 180 : const VisBuffer2 *inputVB = getVii()->getVisBuffer();
838 180 : const Int inputSPW = inputVB->spectralWindows()(0);
839 180 : const uInt width = spwChanbinMap_p[inputSPW];
840 :
841 : // Map input-output channel
842 180 : const uInt nChan = propagatedFlagCube.shape()(1);
843 180 : uInt binCounts = 0;
844 180 : uInt transformedIndex = 0;
845 180 : Vector<uInt> inputOutputChan(nChan);
846 8161 : for (size_t chan_i =0;chan_i<nChan;chan_i++)
847 : {
848 7981 : binCounts += 1;
849 :
850 7981 : if (binCounts > width)
851 : {
852 763 : binCounts = 1;
853 763 : transformedIndex += 1;
854 : }
855 :
856 7981 : inputOutputChan(chan_i) = transformedIndex;
857 : }
858 :
859 : // Propagate chan-avg flags
860 : // Keeping two separate alternatives for 'flagdataFlagPropagation_p' (CAS-12737,
861 : // CAS-9985) until this issue is better settled.
862 180 : if (flagdataFlagPropagation_p)
863 : {
864 180 : cubePropagateFlags(transformedFlagCube, propagatedFlagCube, inputOutputChan,
865 5299791 : [&transformedFlagCube, &propagatedFlagCube]
866 7177470 : (size_t row_i, size_t chan_i, size_t corr_i, uInt outChan) {
867 5299791 : if (transformedFlagCube(corr_i, outChan, row_i))
868 1877679 : propagatedFlagCube(corr_i, chan_i, row_i) = true;
869 5299791 : });
870 : }
871 : else
872 : {
873 0 : cubePropagateFlags(transformedFlagCube, propagatedFlagCube, inputOutputChan,
874 0 : [&transformedFlagCube, &propagatedFlagCube]
875 0 : (size_t row_i, size_t chan_i, size_t corr_i, uInt outChan) {
876 0 : propagatedFlagCube(corr_i, chan_i, row_i) =
877 0 : transformedFlagCube(corr_i, outChan, row_i);
878 0 : });
879 : }
880 180 : }
881 :
882 : //////////////////////////////////////////////////////////////////////////
883 : // ChannelAverageTVIFactory class
884 : //////////////////////////////////////////////////////////////////////////
885 :
886 : // -----------------------------------------------------------------------
887 : //
888 : // -----------------------------------------------------------------------
889 0 : ChannelAverageTVIFactory::ChannelAverageTVIFactory (Record &configuration,
890 0 : ViImplementation2 *inputVii)
891 : {
892 0 : inputVii_p = inputVii;
893 0 : configuration_p = configuration;
894 0 : }
895 :
896 : // -----------------------------------------------------------------------
897 : //
898 : // -----------------------------------------------------------------------
899 0 : vi::ViImplementation2 * ChannelAverageTVIFactory::createVi() const
900 : {
901 0 : return new ChannelAverageTVI(inputVii_p,configuration_p);
902 : }
903 :
904 : //////////////////////////////////////////////////////////////////////////
905 : // ChannelAverageTVILayerFactory class
906 : //////////////////////////////////////////////////////////////////////////
907 :
908 21 : ChannelAverageTVILayerFactory::ChannelAverageTVILayerFactory(Record &configuration) :
909 : ViiLayerFactory(),
910 21 : configuration_p(configuration)
911 21 : {}
912 :
913 : ViImplementation2*
914 21 : ChannelAverageTVILayerFactory::createInstance(ViImplementation2* vii0) const
915 : {
916 : // Make the ChannelAverageTVi2, using supplied ViImplementation2, and return it
917 21 : ViImplementation2 *vii = new ChannelAverageTVI(vii0,configuration_p);
918 21 : return vii;
919 : }
920 :
921 :
922 : //////////////////////////////////////////////////////////////////////////
923 : // ChannelAverageTransformEngine class
924 : //////////////////////////////////////////////////////////////////////////
925 :
926 : // -----------------------------------------------------------------------
927 : //
928 : // -----------------------------------------------------------------------
929 425 : template<class T> ChannelAverageTransformEngine<T>::ChannelAverageTransformEngine( ChannelAverageKernel<T> *kernel,
930 : DataCubeMap *inputData,
931 : DataCubeMap *outputData,
932 : uInt width):
933 : FreqAxisTransformEngine2<T>(inputData,
934 425 : outputData)
935 : {
936 425 : width_p = width;
937 425 : chanAvgKernel_p = kernel;
938 425 : }
939 :
940 : // -----------------------------------------------------------------------
941 : //
942 : // -----------------------------------------------------------------------
943 386 : template<class T> void ChannelAverageTransformEngine<T>::transformAll()
944 : {
945 : // NB: Does NOT implement "parallelCorrAxis" option
946 : // (see, e.g., FreqAxisTVI::transformFreqAxis2(...))
947 :
948 : // Set up the VectorIterators inside the DataCubeMap/Holders
949 386 : inputData_p->setupVecIter();
950 386 : outputData_p->setupVecIter();
951 :
952 : // Iterate implicitly over row and correlation
953 289782 : while (!inputData_p->pastEnd()) {
954 289396 : this->transform(); // processes the current channel axis
955 289396 : inputData_p->next();
956 289396 : outputData_p->next();
957 : }
958 386 : }
959 :
960 :
961 : // -----------------------------------------------------------------------
962 : //
963 : // -----------------------------------------------------------------------
964 289435 : template<class T> void ChannelAverageTransformEngine<T>::transform()
965 : {
966 289435 : uInt startChan = 0;
967 289435 : uInt outChanIndex = 0;
968 289435 : uInt inputSize = inputData_p->getVectorShape()(0);
969 289435 : uInt outputSize = outputData_p->getVectorShape()(0);
970 289435 : uInt tail = inputSize % width_p;
971 289435 : uInt limit = inputSize - tail;
972 1109313 : while (startChan < limit)
973 : {
974 819878 : chanAvgKernel_p->kernel(inputData_p,outputData_p,
975 : startChan,outChanIndex,width_p);
976 819878 : startChan += width_p;
977 819878 : outChanIndex += 1;
978 : }
979 :
980 289435 : if (tail and (outChanIndex <= outputSize - 1) )
981 : {
982 65 : chanAvgKernel_p->kernel(inputData_p,outputData_p,
983 : startChan,outChanIndex,tail);
984 : }
985 :
986 289435 : return;
987 : }
988 :
989 : //////////////////////////////////////////////////////////////////////////
990 : // PlainChannelAverageKernel class
991 : // (numerical averaging, ignoring flags)
992 : //////////////////////////////////////////////////////////////////////////
993 :
994 : // -----------------------------------------------------------------------
995 : //
996 : // -----------------------------------------------------------------------
997 347 : template<class T> void PlainChannelAverageKernel<T>::kernel( DataCubeMap *inputData,
998 : DataCubeMap *outputData,
999 : uInt startInputPos,
1000 : uInt outputPos,
1001 : uInt width)
1002 : {
1003 347 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1004 347 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1005 :
1006 347 : uInt pos = startInputPos + 1;
1007 347 : uInt counts = 1;
1008 347 : T avg = inputVector(startInputPos);
1009 2496 : while (counts < width)
1010 : {
1011 2149 : avg += inputVector(pos);
1012 2149 : counts += 1;
1013 2149 : pos += 1;
1014 : }
1015 :
1016 347 : if (counts > 0)
1017 : {
1018 347 : avg /= counts;
1019 : }
1020 :
1021 347 : outputVector(outputPos) = avg;
1022 :
1023 347 : return;
1024 : }
1025 :
1026 : //////////////////////////////////////////////////////////////////////////
1027 : // FlaggedChannelAverageKernel class
1028 : // (numerical averaging, respecting flags)
1029 : //////////////////////////////////////////////////////////////////////////
1030 :
1031 : // -----------------------------------------------------------------------
1032 : //
1033 : // -----------------------------------------------------------------------
1034 25776 : template<class T> void FlaggedChannelAverageKernel<T>::kernel(DataCubeMap *inputData,
1035 : DataCubeMap *outputData,
1036 : uInt startInputPos,
1037 : uInt outputPos,
1038 : uInt width)
1039 : {
1040 25776 : T avg = 0;
1041 25776 : T normalization = 0;
1042 25776 : uInt inputPos = 0;
1043 25776 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1044 25776 : Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
1045 25776 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1046 25776 : Bool accumulatorFlag = inputFlagVector(startInputPos);
1047 :
1048 850608 : for (uInt sample_i=0;sample_i<width;sample_i++)
1049 : {
1050 : // Get input index
1051 824832 : inputPos = startInputPos + sample_i;
1052 :
1053 : // true/true or false/false
1054 824832 : if (accumulatorFlag == inputFlagVector(inputPos))
1055 : {
1056 824832 : normalization += 1.0f;
1057 824832 : avg += inputVector(inputPos);
1058 : }
1059 : // true/false: Reset accumulation when accumulator switches from flagged to unflag
1060 0 : else if ( (accumulatorFlag == true) and (inputFlagVector(inputPos) == false) )
1061 : {
1062 0 : accumulatorFlag = false;
1063 0 : normalization = 1.0f;
1064 0 : avg = inputVector(inputPos);
1065 : }
1066 :
1067 : }
1068 :
1069 :
1070 : // Apply normalization factor
1071 25776 : if (normalization > 0)
1072 : {
1073 25776 : avg /= normalization;
1074 25776 : outputVector(outputPos) = avg;
1075 : }
1076 : // If all weights are zero set accumulatorFlag to true
1077 : else
1078 : {
1079 0 : accumulatorFlag = true;
1080 0 : outputVector(outputPos) = 0; // If all weights are zero then the avg is 0 too
1081 : }
1082 :
1083 51552 : return;
1084 : }
1085 :
1086 : //////////////////////////////////////////////////////////////////////////
1087 : // WeightedChannelAverageKernel class
1088 : // (weighted averaging, respecting flags)
1089 : //////////////////////////////////////////////////////////////////////////
1090 :
1091 : // -----------------------------------------------------------------------
1092 : //
1093 : // -----------------------------------------------------------------------
1094 396900 : template<class T> void WeightedChannelAverageKernel<T>::kernel( DataCubeMap *inputData,
1095 : DataCubeMap *outputData,
1096 : uInt startInputPos,
1097 : uInt outputPos,
1098 : uInt width)
1099 : {
1100 396900 : T avg = 0;
1101 396900 : T normalization = 0;
1102 396900 : uInt inputPos = 0;
1103 396900 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1104 396900 : Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
1105 396900 : Vector<Float> &inputWeightVector = inputData->getVector<Float>(MS::WEIGHT_SPECTRUM);
1106 396900 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1107 396900 : Bool accumulatorFlag = inputFlagVector(startInputPos);
1108 :
1109 6797992 : for (uInt sample_i=0;sample_i<width;sample_i++)
1110 : {
1111 : // Get input index
1112 6401092 : inputPos = startInputPos + sample_i;
1113 :
1114 6401092 : Float& wt=inputWeightVector(inputPos);
1115 :
1116 : // true/true or false/false
1117 6401092 : if (accumulatorFlag == inputFlagVector(inputPos))
1118 : {
1119 6399764 : normalization += wt;
1120 6399764 : avg += inputVector(inputPos)*wt;
1121 : }
1122 : // true/false: Reset accumulation when accumulator switches from flagged to unflag
1123 1328 : else if ( (accumulatorFlag == true) and (inputFlagVector(inputPos) == false) )
1124 : {
1125 1328 : accumulatorFlag = false;
1126 1328 : normalization = inputWeightVector(inputPos);
1127 1328 : avg = inputVector(inputPos)*inputWeightVector(inputPos);
1128 : }
1129 : }
1130 :
1131 :
1132 : // Apply normalization factor
1133 396900 : if (normalization > 0)
1134 : {
1135 396900 : avg /= normalization;
1136 396900 : outputVector(outputPos) = avg;
1137 : }
1138 : // If all weights are zero set accumulatorFlag to true
1139 : else
1140 : {
1141 0 : accumulatorFlag = true;
1142 0 : outputVector(outputPos) = 0; // If all weights are zero then the avg is 0 too
1143 : }
1144 :
1145 793800 : return;
1146 : }
1147 :
1148 : //////////////////////////////////////////////////////////////////////////
1149 : // LogicalANDKernel class
1150 : //////////////////////////////////////////////////////////////////////////
1151 :
1152 : // -----------------------------------------------------------------------
1153 : //
1154 : // -----------------------------------------------------------------------
1155 396910 : template<class T> void LogicalANDKernel<T>::kernel( DataCubeMap *inputData,
1156 : DataCubeMap *outputData,
1157 : uInt startInputPos,
1158 : uInt outputPos,
1159 : uInt width)
1160 : {
1161 396910 : Vector<Bool> &inputVector = inputData->getVector<Bool>(MS::FLAG);
1162 396910 : Vector<Bool> &outputVector = outputData->getVector<Bool>(MS::FLAG);
1163 :
1164 396910 : Bool outputFlag = true;
1165 398254 : for (uInt sample_i=0;sample_i<width;sample_i++)
1166 : {
1167 398238 : if (not inputVector(startInputPos+sample_i))
1168 : {
1169 396894 : outputFlag = false;
1170 396894 : break;
1171 : }
1172 : }
1173 396910 : outputVector(outputPos) = outputFlag;
1174 :
1175 396910 : return;
1176 : }
1177 :
1178 : //////////////////////////////////////////////////////////////////////////
1179 : // ChannelAccumulationKernel class
1180 : // (numerical accumulation, respecting flags)
1181 : //////////////////////////////////////////////////////////////////////////
1182 :
1183 : // -----------------------------------------------------------------------
1184 : //
1185 : // -----------------------------------------------------------------------
1186 10 : template<class T> void ChannelAccumulationKernel<T>::kernel( DataCubeMap *inputData,
1187 : DataCubeMap *outputData,
1188 : uInt startInputPos,
1189 : uInt outputPos,
1190 : uInt width)
1191 : {
1192 10 : T acc = 0;
1193 10 : uInt inputPos = 0;
1194 10 : Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
1195 10 : Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
1196 10 : Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
1197 10 : Bool accumulatorFlag = inputFlagVector(startInputPos);
1198 :
1199 304 : for (uInt sample_i=0;sample_i<width;sample_i++)
1200 : {
1201 : // Get input index
1202 294 : inputPos = startInputPos + sample_i;
1203 :
1204 : // true/true or false/false
1205 294 : if (accumulatorFlag == inputFlagVector(inputPos))
1206 : {
1207 294 : acc += inputVector(inputPos);
1208 : }
1209 : // true/false: Reset accumulation when accumulator switches from flagged to unflag
1210 0 : else if ( (accumulatorFlag == true) and (inputFlagVector(inputPos) == false) )
1211 : {
1212 0 : accumulatorFlag = false;
1213 0 : acc = inputVector(inputPos);
1214 : }
1215 : }
1216 :
1217 :
1218 10 : outputVector(outputPos) = acc;
1219 :
1220 10 : return;
1221 : }
1222 :
1223 : } //# NAMESPACE VI - END
1224 :
1225 : } //# NAMESPACE CASA - END
1226 :
1227 :
|
