LCOV - code coverage report
Current view: top level - mstransform/TVI - ChannelAverageTVI.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 0 484 0.0 %
Date: 2024-10-29 13:38:20 Functions: 0 40 0.0 %

          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           0 : ChannelAverageTVI::ChannelAverageTVI(   ViImplementation2 * inputVii,
      44           0 :                                                                                 const Record &configuration):
      45           0 :                                                                                 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           0 :     if (not parseConfiguration(configuration))
      51           0 :         throw AipsError("Error parsing ChannelAverageTVI configuration");
      52             : 
      53           0 :     if (inputVii == nullptr)
      54           0 :         throw AipsError("Input Vi is empty");
      55             : 
      56           0 :     initialize();
      57             : 
      58           0 :     return;
      59           0 : }
      60             : 
      61             : // -----------------------------------------------------------------------
      62             : //
      63             : // -----------------------------------------------------------------------
      64           0 : Bool ChannelAverageTVI::parseConfiguration(const Record &configuration)
      65             : {
      66           0 :         int exists = -1;
      67           0 :         Bool ret = true;
      68             : 
      69             :         // Parse chanbin parameter (mandatory)
      70           0 :         exists = -1;
      71           0 :         exists = configuration.fieldNumber ("chanbin");
      72           0 :         if (exists >= 0)
      73             :         {
      74           0 :                 if ( configuration.type(exists) == casacore::TpInt )
      75             :                 {
      76             :                         Int freqbin;
      77           0 :                         configuration.get (exists, freqbin);
      78           0 :                         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           0 :                 logger_p << LogIO::NORMAL << LogOrigin("ChannelAverageTVI", __FUNCTION__)
      93           0 :                                 << "Channel bin is " << chanbin_p << LogIO::POST;
      94           0 :                 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           0 :         exists = configuration.fieldNumber ("flagdataFlagPropagation");
     108           0 :         flagdataFlagPropagation_p = exists >= 0;
     109             : 
     110             :         // Check consistency between chanbin vector and selected SPW/Chan map
     111           0 :         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           0 :         return ret;
     120             : }
     121             : 
     122             : // -----------------------------------------------------------------------
     123             : //
     124             : // -----------------------------------------------------------------------
     125           0 : void ChannelAverageTVI::initialize()
     126             : {
     127             :         // Populate nchan input-output maps
     128             :         Int spw;
     129           0 :         uInt spw_idx = 0;
     130           0 :         map<Int,vector<Int> >::iterator iter;
     131             : 
     132           0 :         for(iter=spwInpChanIdxMap_p.begin();iter!=spwInpChanIdxMap_p.end();iter++)
     133             :         {
     134           0 :                 spw = iter->first;
     135             : 
     136             :                 // No averaging when user specifies 0
     137           0 :                 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           0 :                 else if ((uInt)chanbin_p(spw_idx) <= 1)
     147             :                 {
     148           0 :                         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           0 :                                         << " existing/selected channels: " << iter->second.size()
     152           0 :                                         << LogIO::POST;
     153             : 
     154           0 :                         spwChanbinMap_p[spw] = iter->second.size();
     155             :                 }
     156             :                 // Make sure that chanbin does not exceed number of selected channels
     157           0 :                 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           0 :                         spwChanbinMap_p[spw] = chanbin_p(spw_idx);
     171             :                 }
     172             : 
     173             :                 // Calculate number of output channels per spw
     174           0 :                 spwOutChanNumMap_p[spw] = spwInpChanIdxMap_p[spw].size() / spwChanbinMap_p[spw];
     175           0 :                 if (spwInpChanIdxMap_p[spw].size() % spwChanbinMap_p[spw] > 0) spwOutChanNumMap_p[spw] += 1;
     176             : 
     177           0 :                 spw_idx++;
     178             :         }
     179             : 
     180           0 :         return;
     181             : }
     182             : 
     183             : #define DOJUSTO false
     184             : // -----------------------------------------------------------------------
     185             : //
     186             : // -----------------------------------------------------------------------
     187           0 : void ChannelAverageTVI::flag(Cube<Bool>& flagCube) const
     188             : {
     189             : 
     190             :         // Pass-thru for single-channel case
     191           0 :         if (getVii()->visibilityShape()[1]==1) {
     192           0 :           getVii()->flag(flagCube);
     193           0 :           return;
     194             :         }
     195             :     
     196             :         // Get input VisBuffer and SPW
     197           0 :         VisBuffer2 *vb = getVii()->getVisBuffer();
     198           0 :         Int inputSPW = vb->spectralWindows()(0);
     199             : 
     200             :         // Pass-thru for chanbin=1 case:
     201           0 :         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           0 :         vb->flagCube();
     209           0 :         Double time0=omp_get_wtime();
     210             : #endif
     211             : 
     212             : 
     213             :         // Reshape output data before passing it to the DataCubeHolder
     214           0 :         flagCube.resize(getVisBuffer()->getShape(),false);
     215             : 
     216             :         // Gather input data
     217           0 :         DataCubeMap inputData;
     218           0 :         DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
     219           0 :         inputData.add(MS::FLAG,inputFlagCubeHolder);
     220             : 
     221             :         // Gather output data
     222           0 :         DataCubeMap outputData;
     223           0 :         DataCubeHolder<Bool> outputFlagCubeHolder(flagCube);
     224           0 :         outputData.add(MS::FLAG,outputFlagCubeHolder);
     225             : 
     226             :         // Configure Transformation Engine
     227           0 :         LogicalANDKernel<Bool> kernel;
     228           0 :         uInt width = spwChanbinMap_p[inputSPW];
     229           0 :         ChannelAverageTransformEngine<Bool> transformer(&kernel,&inputData,&outputData,width);
     230             : 
     231             :         // Transform data
     232             :         if (DOJUSTO) {
     233             :           transformFreqAxis2(vb->getShape(),transformer);
     234             :         } else {
     235           0 :           transformer.transformAll();
     236             :         }
     237             : 
     238             : #ifdef _OPENMP
     239             :         // Accumulate elapsed time
     240           0 :         Tfl_+=omp_get_wtime()-time0;
     241             : #endif
     242             :           
     243           0 :           return;
     244           0 : }
     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           0 : void ChannelAverageTVI::visibilityObserved (Cube<Complex> & vis) const
     304             : {
     305             : 
     306             :         // Pass-thru for single-channel case
     307           0 :         if (getVii()->visibilityShape()[1]==1) {
     308           0 :           getVii()->visibilityObserved(vis);
     309           0 :           return;
     310             :         }
     311             : 
     312             :         // Get input VisBuffer and SPW
     313           0 :         VisBuffer2 *vb = getVii()->getVisBuffer();
     314           0 :         Int inputSPW = vb->spectralWindows()(0);
     315             : 
     316             :         // Reshape output data before passing it to the DataCubeHolder
     317           0 :         vis.resize(getVisBuffer()->getShape(),false);
     318             : 
     319             :         // Get weightSpectrum from sigmaSpectrum
     320           0 :         Cube<Float> weightSpFromSigmaSp;
     321           0 :         weightSpFromSigmaSp.resize(vb->sigmaSpectrum().shape(),false);
     322           0 :         weightSpFromSigmaSp = vb->sigmaSpectrum(); // = Operator makes a copy
     323           0 :         arrayTransformInPlace (weightSpFromSigmaSp,sigmaToWeight);
     324             : 
     325             :         // Gather input data
     326           0 :         DataCubeMap inputData;
     327           0 :         DataCubeHolder<Complex> inputVisCubeHolder(vb->visCube());
     328           0 :         DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
     329           0 :         DataCubeHolder<Float> weightCubeHolder(weightSpFromSigmaSp);
     330           0 :         inputData.add(MS::DATA,inputVisCubeHolder);
     331           0 :         inputData.add(MS::FLAG,inputFlagCubeHolder);
     332           0 :         inputData.add(MS::WEIGHT_SPECTRUM,weightCubeHolder);
     333             : 
     334             :         // Gather output data
     335           0 :         DataCubeMap outputData;
     336           0 :         DataCubeHolder<Complex> outputVisCubeHolder(vis);
     337           0 :         outputData.add(MS::DATA,outputVisCubeHolder);
     338             : 
     339             :         // Configure Transformation Engine
     340           0 :         uInt width = spwChanbinMap_p[inputSPW];
     341           0 :         WeightedChannelAverageKernel<Complex> kernel;
     342           0 :         ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
     343             : 
     344             :         // Transform data
     345             :         if (DOJUSTO) {
     346             :           transformFreqAxis2(vb->getShape(),transformer);
     347             :         }
     348             :         else {
     349           0 :           transformer.transformAll();
     350             :         }
     351             : 
     352           0 :         return;
     353           0 : }
     354             : 
     355             : // -----------------------------------------------------------------------
     356             : //
     357             : // -----------------------------------------------------------------------
     358           0 : void ChannelAverageTVI::visibilityCorrected (Cube<Complex> & vis) const
     359             : {
     360             : 
     361             :         // Pass-thru for single-channel case
     362           0 :         if (getVii()->visibilityShape()[1]==1) {
     363           0 :           getVii()->visibilityCorrected(vis);
     364           0 :           return;
     365             :         }
     366             : 
     367             :         // Get input VisBuffer and SPW
     368           0 :         VisBuffer2 *vb = getVii()->getVisBuffer();
     369           0 :         Int inputSPW = vb->spectralWindows()(0);
     370             : 
     371             :         // Pass-thru for chanbin=1 case:
     372           0 :         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           0 :         vb->visCubeCorrected();
     381           0 :         vb->flagCube();
     382           0 :         vb->weightSpectrum();
     383           0 :         Double time0=omp_get_wtime();
     384             : #endif
     385             : 
     386             : 
     387             :         // Reshape output data before passing it to the DataCubeHolder
     388           0 :         vis.resize(getVisBuffer()->getShape(),false);
     389             : 
     390             :         // Gather input data
     391           0 :         DataCubeMap inputData;
     392           0 :         DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeCorrected());
     393           0 :         DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
     394           0 :         DataCubeHolder<Float> weightCubeHolder(vb->weightSpectrum());
     395           0 :         inputData.add(MS::DATA,inputVisCubeHolder);
     396           0 :         inputData.add(MS::FLAG,inputFlagCubeHolder);
     397           0 :         inputData.add(MS::WEIGHT_SPECTRUM,weightCubeHolder);
     398             : 
     399             :         // Gather output data
     400           0 :         DataCubeMap outputData;
     401           0 :         DataCubeHolder<Complex> outputVisCubeHolder(vis);
     402           0 :         outputData.add(MS::DATA,outputVisCubeHolder);
     403             : 
     404             :         // Configure Transformation Engine
     405           0 :         uInt width = spwChanbinMap_p[inputSPW];
     406           0 :         WeightedChannelAverageKernel<Complex> kernel;
     407           0 :         ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
     408             : 
     409             :         // Transform data
     410             :         if (DOJUSTO) {
     411             :           transformFreqAxis2(vb->getShape(),transformer);
     412             :         } else {
     413           0 :           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           0 :         Tcd_+=omp_get_wtime()-time0;
     463             : #endif
     464             : 
     465           0 :         return;
     466           0 : }
     467             : 
     468             : // -----------------------------------------------------------------------
     469             : //
     470             : // -----------------------------------------------------------------------
     471           0 : void ChannelAverageTVI::visibilityModel (Cube<Complex> & vis) const
     472             : {
     473             : 
     474             :         // Pass-thru for single-channel case
     475           0 :         if (getVii()->visibilityShape()[1]==1) {
     476           0 :           getVii()->visibilityModel(vis);
     477           0 :           return;
     478             :         }
     479             : 
     480             :         // Get input VisBuffer and SPW
     481           0 :         VisBuffer2 *vb = getVii()->getVisBuffer();
     482           0 :         Int inputSPW = vb->spectralWindows()(0);
     483             : 
     484             :         // Pass-thru for chanbin=1 case:
     485           0 :         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           0 :         vb->visCubeModel();
     493           0 :         vb->flagCube();
     494           0 :         vb->weightSpectrum();
     495           0 :         Double time0=omp_get_wtime();
     496             : #endif
     497             : 
     498             : 
     499             :         // Reshape output data before passing it to the DataCubeHolder
     500           0 :         vis.resize(getVisBuffer()->getShape(),false);
     501             : 
     502             :         // Gather input data
     503           0 :         DataCubeMap inputData;
     504           0 :         DataCubeHolder<Complex> inputVisCubeHolder(vb->visCubeModel());
     505           0 :         DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
     506           0 :         inputData.add(MS::DATA,inputVisCubeHolder);
     507           0 :         inputData.add(MS::FLAG,inputFlagCubeHolder);
     508             : 
     509             :         // Gather output data
     510           0 :         DataCubeMap outputData;
     511           0 :         DataCubeHolder<Complex> outputVisCubeHolder(vis);
     512           0 :         outputData.add(MS::DATA,outputVisCubeHolder);
     513             : 
     514             :         // Configure Transformation Engine
     515           0 :         uInt width = spwChanbinMap_p[inputSPW];
     516           0 :         FlaggedChannelAverageKernel<Complex> kernel;
     517           0 :         ChannelAverageTransformEngine<Complex> transformer(&kernel,&inputData,&outputData,width);
     518             : 
     519             :         // Transform data
     520             :         if (DOJUSTO) {
     521             :           transformFreqAxis2(vb->getShape(),transformer);
     522             :         } else {
     523           0 :           transformer.transformAll();
     524             :         }
     525             : 
     526             : #ifdef _OPENMP
     527             :         // Accumulate elapsed time
     528           0 :         Tmd_+=omp_get_wtime()-time0;
     529             : #endif
     530             : 
     531           0 :         return;
     532           0 : }
     533             : 
     534             : // -----------------------------------------------------------------------
     535             : //
     536             : // -----------------------------------------------------------------------
     537           0 : void ChannelAverageTVI::weightSpectrum(Cube<Float> &weightSp) const
     538             : {
     539             : 
     540             :         // Pass-thru for single-channel case
     541           0 :         if (getVii()->visibilityShape()[1]==1) {
     542           0 :           getVii()->weightSpectrum(weightSp);
     543           0 :           return;
     544             :         }
     545             : 
     546             :         // Get input VisBuffer and SPW
     547           0 :         VisBuffer2 *vb = getVii()->getVisBuffer();
     548           0 :         Int inputSPW = vb->spectralWindows()(0);
     549             : 
     550             :         // Pass-thru for chanbin=1 case:
     551           0 :         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           0 :         vb->weightSpectrum();
     559           0 :         vb->flagCube();
     560           0 :         Double time0=omp_get_wtime();
     561             : #endif
     562             : 
     563             : 
     564             :         // Reshape output data before passing it to the DataCubeHolder
     565           0 :         weightSp.resize(getVisBuffer()->getShape(),false);
     566             : 
     567             :         // Gather input data
     568           0 :         DataCubeMap inputData;
     569           0 :         DataCubeHolder<Float> inputWeightCubeHolder(vb->weightSpectrum());
     570           0 :         DataCubeHolder<Bool> inputFlagCubeHolder(vb->flagCube());
     571           0 :         inputData.add(MS::DATA,inputWeightCubeHolder);
     572           0 :         inputData.add(MS::FLAG,inputFlagCubeHolder);
     573             : 
     574             :         // Gather output data
     575           0 :         DataCubeMap outputData;
     576           0 :         DataCubeHolder<Float> outputWeightCubeHolder(weightSp);
     577           0 :         outputData.add(MS::DATA,outputWeightCubeHolder);
     578             : 
     579             :         // Configure Transformation Engine
     580           0 :         uInt width = spwChanbinMap_p[inputSPW];
     581           0 :         ChannelAccumulationKernel<Float> kernel;
     582           0 :         ChannelAverageTransformEngine<Float> transformer(&kernel,&inputData,&outputData,width);
     583             : 
     584             :         // Transform data
     585             :         if (DOJUSTO) {
     586             :           transformFreqAxis2(vb->getShape(),transformer);
     587             :         } else {
     588           0 :           transformer.transformAll();
     589             :         }
     590             : 
     591             : #ifdef _OPENMP
     592             :         // Accumulate elapsed time
     593           0 :         Tws_+=omp_get_wtime()-time0;
     594             : #endif
     595             : 
     596           0 :         return;
     597           0 : }
     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           0 : 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           0 :         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           0 :         Vector<Double> inputFrequencies = getVii()->getFrequencies(time,frameOfReference,
     687           0 :                                                                    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           0 :         if (inputFrequencies.nelements()==1)
     693           0 :           return inputFrequencies;
     694             : 
     695             :         // Produce output (transformed) frequencies
     696           0 :         Vector<Double> outputFrecuencies(spwOutChanNumMap_p[spectralWindowId],0);
     697             : 
     698             :         // Gather input data
     699           0 :         DataCubeMap inputData;
     700           0 :         DataCubeHolder<Double> inputFrecuenciesHolder(inputFrequencies);
     701           0 :         inputData.add(MS::DATA,inputFrecuenciesHolder);
     702             : 
     703             :         // Gather output data
     704           0 :         DataCubeMap outputData;
     705           0 :         DataCubeHolder<Double> outputFrecuenciesHolder(outputFrecuencies);
     706           0 :         outputData.add(MS::DATA,outputFrecuenciesHolder);
     707             : 
     708             :         // Configure Transformation Engine
     709           0 :         PlainChannelAverageKernel<Double> kernel;
     710           0 :         uInt width = spwChanbinMap_p[spectralWindowId];
     711           0 :         ChannelAverageTransformEngine<Double> transformer(&kernel,&inputData,&outputData,width);
     712             : 
     713             :         // Transform data
     714           0 :         transformer.transform();
     715             : 
     716           0 :         return outputFrecuencies;
     717           0 : }
     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           0 : void ChannelAverageTVI::writeFlag (const Cube<Bool> & flag)
     772             : {
     773             :         // Create a flag cube with the input VI shape
     774           0 :         Cube<Bool> propagatedFlagCube;
     775           0 :         propagatedFlagCube = getVii()->getVisBuffer()->flagCube();
     776             : 
     777             :         // Propagate flags from the input cube to the propagated flag cube
     778           0 :         propagateChanAvgFlags(flag,propagatedFlagCube);
     779             : 
     780             :         // Pass propagated flag cube downstream for further propagation and/or writting
     781           0 :         getVii()->writeFlag(propagatedFlagCube);
     782             : 
     783           0 :         return;
     784           0 : }
     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           0 : void cubePropagateFlags(const Cube<Bool> &transformedFlagCube,
     802             :                         Cube<Bool> &propagatedFlagCube,
     803             :                         const Vector<uInt> &inputOutputChan, Functor propagate)
     804             : {
     805             :     // Get propagated (input) shape
     806           0 :     const auto inputShape = propagatedFlagCube.shape();
     807           0 :     const uInt nCorr = inputShape(0);
     808           0 :     const uInt nChan = inputShape(1);
     809           0 :     const uInt nRows = inputShape(2);
     810             : 
     811             :     // Get transformed (output) shape
     812           0 :     const auto transformedShape = transformedFlagCube.shape();
     813           0 :     const auto nTransChan = transformedShape(1);
     814             : 
     815             :     uInt outChan;
     816           0 :     for (size_t row_i =0;row_i<nRows;row_i++)
     817             :     {
     818           0 :         for (size_t chan_i =0;chan_i<nChan;chan_i++)
     819             :         {
     820           0 :             outChan = inputOutputChan(chan_i);
     821           0 :             if (outChan < nTransChan)  // outChan >= nChan  may happen when channels are dropped
     822             :             {
     823           0 :                 for (size_t corr_i =0;corr_i<nCorr;corr_i++)
     824           0 :                     propagate(row_i, chan_i, corr_i, outChan);
     825             :             }
     826             :         }
     827             :     }
     828           0 : }
     829             : 
     830             : // -----------------------------------------------------------------------
     831             : //
     832             : // -----------------------------------------------------------------------
     833           0 : void ChannelAverageTVI::propagateChanAvgFlags (const Cube<Bool> &transformedFlagCube,
     834             :                                                Cube<Bool> &propagatedFlagCube)
     835             : {
     836             :     // Get current SPW and chanbin
     837           0 :     const VisBuffer2 *inputVB = getVii()->getVisBuffer();
     838           0 :     const Int inputSPW = inputVB->spectralWindows()(0);
     839           0 :     const uInt width = spwChanbinMap_p[inputSPW];
     840             : 
     841             :     // Map input-output channel
     842           0 :     const uInt nChan = propagatedFlagCube.shape()(1);
     843           0 :     uInt binCounts = 0;
     844           0 :     uInt transformedIndex = 0;
     845           0 :     Vector<uInt> inputOutputChan(nChan);
     846           0 :     for (size_t chan_i =0;chan_i<nChan;chan_i++)
     847             :     {
     848           0 :         binCounts += 1;
     849             : 
     850           0 :         if (binCounts > width)
     851             :         {
     852           0 :             binCounts = 1;
     853           0 :             transformedIndex += 1;
     854             :         }
     855             : 
     856           0 :         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           0 :     if (flagdataFlagPropagation_p)
     863             :     {
     864           0 :         cubePropagateFlags(transformedFlagCube, propagatedFlagCube, inputOutputChan,
     865           0 :                            [&transformedFlagCube, &propagatedFlagCube]
     866           0 :                            (size_t row_i, size_t chan_i, size_t corr_i, uInt outChan) {
     867           0 :                                if (transformedFlagCube(corr_i, outChan, row_i))
     868           0 :                                    propagatedFlagCube(corr_i, chan_i, row_i) = true;
     869           0 :                            });
     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           0 : }
     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           0 : ChannelAverageTVILayerFactory::ChannelAverageTVILayerFactory(Record &configuration) :
     909             :   ViiLayerFactory(),
     910           0 :   configuration_p(configuration)
     911           0 : {}
     912             : 
     913             : ViImplementation2* 
     914           0 : ChannelAverageTVILayerFactory::createInstance(ViImplementation2* vii0) const 
     915             : {
     916             :   // Make the ChannelAverageTVi2, using supplied ViImplementation2, and return it
     917           0 :   ViImplementation2 *vii = new ChannelAverageTVI(vii0,configuration_p);
     918           0 :   return vii;
     919             : }
     920             : 
     921             : 
     922             : //////////////////////////////////////////////////////////////////////////
     923             : // ChannelAverageTransformEngine class
     924             : //////////////////////////////////////////////////////////////////////////
     925             : 
     926             : // -----------------------------------------------------------------------
     927             : //
     928             : // -----------------------------------------------------------------------
     929           0 : template<class T> ChannelAverageTransformEngine<T>::ChannelAverageTransformEngine(  ChannelAverageKernel<T> *kernel,
     930             :                                                                                                                                                                         DataCubeMap *inputData,
     931             :                                                                                                                                                                         DataCubeMap *outputData,
     932             :                                                                                                                                                                         uInt width):
     933             :                                                                                                                                                                         FreqAxisTransformEngine2<T>(inputData,
     934           0 :                                                                                                                                                                                                                                 outputData)
     935             : {
     936           0 :         width_p = width;
     937           0 :         chanAvgKernel_p = kernel;
     938           0 : }
     939             : 
     940             : // -----------------------------------------------------------------------
     941             : //
     942             : // -----------------------------------------------------------------------
     943           0 : 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           0 :   inputData_p->setupVecIter();
     950           0 :   outputData_p->setupVecIter();
     951             : 
     952             :   // Iterate implicitly over row and correlation
     953           0 :   while (!inputData_p->pastEnd()) {
     954           0 :     this->transform();   // processes the current channel axis
     955           0 :     inputData_p->next();
     956           0 :     outputData_p->next();
     957             :   }
     958           0 : }
     959             : 
     960             : 
     961             : // -----------------------------------------------------------------------
     962             : //
     963             : // -----------------------------------------------------------------------
     964           0 : template<class T> void ChannelAverageTransformEngine<T>::transform()
     965             : {
     966           0 :         uInt startChan = 0;
     967           0 :         uInt outChanIndex = 0;
     968           0 :         uInt inputSize = inputData_p->getVectorShape()(0);
     969           0 :         uInt outputSize = outputData_p->getVectorShape()(0);
     970           0 :         uInt tail = inputSize % width_p;
     971           0 :         uInt limit = inputSize - tail;
     972           0 :         while (startChan < limit)
     973             :         {
     974           0 :                 chanAvgKernel_p->kernel(inputData_p,outputData_p,
     975             :                                                                 startChan,outChanIndex,width_p);
     976           0 :                 startChan += width_p;
     977           0 :                 outChanIndex += 1;
     978             :         }
     979             : 
     980           0 :         if (tail and (outChanIndex <= outputSize - 1) )
     981             :         {
     982           0 :                 chanAvgKernel_p->kernel(inputData_p,outputData_p,
     983             :                                                                 startChan,outChanIndex,tail);
     984             :         }
     985             : 
     986           0 :         return;
     987             : }
     988             : 
     989             : //////////////////////////////////////////////////////////////////////////
     990             : // PlainChannelAverageKernel class
     991             : //   (numerical averaging, ignoring flags)
     992             : //////////////////////////////////////////////////////////////////////////
     993             : 
     994             : // -----------------------------------------------------------------------
     995             : //
     996             : // -----------------------------------------------------------------------
     997           0 : template<class T> void PlainChannelAverageKernel<T>::kernel(        DataCubeMap *inputData,
     998             :                                                                                                                                 DataCubeMap *outputData,
     999             :                                                                                                                                 uInt startInputPos,
    1000             :                                                                                                                                 uInt outputPos,
    1001             :                                                                                                                                 uInt width)
    1002             : {
    1003           0 :         Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
    1004           0 :         Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
    1005             : 
    1006           0 :         uInt pos = startInputPos + 1;
    1007           0 :         uInt counts = 1;
    1008           0 :         T avg = inputVector(startInputPos);
    1009           0 :         while (counts < width)
    1010             :         {
    1011           0 :                 avg += inputVector(pos);
    1012           0 :                 counts += 1;
    1013           0 :                 pos += 1;
    1014             :         }
    1015             : 
    1016           0 :         if (counts > 0)
    1017             :         {
    1018           0 :                 avg /= counts;
    1019             :         }
    1020             : 
    1021           0 :         outputVector(outputPos) = avg;
    1022             : 
    1023           0 :         return;
    1024             : }
    1025             : 
    1026             : //////////////////////////////////////////////////////////////////////////
    1027             : // FlaggedChannelAverageKernel class
    1028             : //   (numerical averaging, respecting flags)
    1029             : //////////////////////////////////////////////////////////////////////////
    1030             : 
    1031             : // -----------------------------------------------------------------------
    1032             : //
    1033             : // -----------------------------------------------------------------------
    1034           0 : template<class T> void FlaggedChannelAverageKernel<T>::kernel(DataCubeMap *inputData,
    1035             :                                                                                                                                 DataCubeMap *outputData,
    1036             :                                                                                                                                 uInt startInputPos,
    1037             :                                                                                                                                 uInt outputPos,
    1038             :                                                                                                                                 uInt width)
    1039             : {
    1040           0 :         T avg = 0;
    1041           0 :         T normalization = 0;
    1042           0 :         uInt inputPos = 0;
    1043           0 :         Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
    1044           0 :         Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
    1045           0 :         Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
    1046           0 :         Bool accumulatorFlag = inputFlagVector(startInputPos);
    1047             : 
    1048           0 :         for (uInt sample_i=0;sample_i<width;sample_i++)
    1049             :         {
    1050             :                 // Get input index
    1051           0 :                 inputPos = startInputPos + sample_i;
    1052             : 
    1053             :                 // true/true or false/false
    1054           0 :                 if (accumulatorFlag == inputFlagVector(inputPos))
    1055             :                 {
    1056           0 :                         normalization += 1.0f;
    1057           0 :                         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           0 :         if (normalization > 0)
    1072             :         {
    1073           0 :                 avg /= normalization;
    1074           0 :                 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           0 :         return;
    1084             : }
    1085             : 
    1086             : //////////////////////////////////////////////////////////////////////////
    1087             : // WeightedChannelAverageKernel class
    1088             : //   (weighted averaging, respecting flags)
    1089             : //////////////////////////////////////////////////////////////////////////
    1090             : 
    1091             : // -----------------------------------------------------------------------
    1092             : //
    1093             : // -----------------------------------------------------------------------
    1094           0 : template<class T> void WeightedChannelAverageKernel<T>::kernel(     DataCubeMap *inputData,
    1095             :                                                                                                                                 DataCubeMap *outputData,
    1096             :                                                                                                                                 uInt startInputPos,
    1097             :                                                                                                                                 uInt outputPos,
    1098             :                                                                                                                                 uInt width)
    1099             : {
    1100           0 :         T avg = 0;
    1101           0 :         T normalization = 0;
    1102           0 :         uInt inputPos = 0;
    1103           0 :         Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
    1104           0 :         Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
    1105           0 :         Vector<Float> &inputWeightVector = inputData->getVector<Float>(MS::WEIGHT_SPECTRUM);
    1106           0 :         Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
    1107           0 :         Bool accumulatorFlag = inputFlagVector(startInputPos);
    1108             : 
    1109           0 :         for (uInt sample_i=0;sample_i<width;sample_i++)
    1110             :         {
    1111             :                 // Get input index
    1112           0 :                 inputPos = startInputPos + sample_i;
    1113             : 
    1114           0 :                 Float& wt=inputWeightVector(inputPos);
    1115             : 
    1116             :                 // true/true or false/false
    1117           0 :                 if (accumulatorFlag == inputFlagVector(inputPos))
    1118             :                 {
    1119           0 :                         normalization += wt;
    1120           0 :                         avg += inputVector(inputPos)*wt;
    1121             :                 }
    1122             :                 // true/false: Reset accumulation when accumulator switches from flagged to unflag
    1123           0 :                 else if ( (accumulatorFlag == true) and (inputFlagVector(inputPos) == false) )
    1124             :                 {
    1125           0 :                         accumulatorFlag = false;
    1126           0 :                         normalization = inputWeightVector(inputPos);
    1127           0 :                         avg = inputVector(inputPos)*inputWeightVector(inputPos);
    1128             :                 }
    1129             :         }
    1130             : 
    1131             : 
    1132             :         // Apply normalization factor
    1133           0 :         if (normalization > 0)
    1134             :         {
    1135           0 :                 avg /= normalization;
    1136           0 :                 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           0 :         return;
    1146             : }
    1147             : 
    1148             : //////////////////////////////////////////////////////////////////////////
    1149             : // LogicalANDKernel class
    1150             : //////////////////////////////////////////////////////////////////////////
    1151             : 
    1152             : // -----------------------------------------------------------------------
    1153             : //
    1154             : // -----------------------------------------------------------------------
    1155           0 : template<class T> void LogicalANDKernel<T>::kernel( DataCubeMap *inputData,
    1156             :                                                                                                         DataCubeMap *outputData,
    1157             :                                                                                                         uInt startInputPos,
    1158             :                                                                                                         uInt outputPos,
    1159             :                                                                                                         uInt width)
    1160             : {
    1161           0 :         Vector<Bool> &inputVector = inputData->getVector<Bool>(MS::FLAG);
    1162           0 :         Vector<Bool> &outputVector = outputData->getVector<Bool>(MS::FLAG);
    1163             : 
    1164           0 :         Bool outputFlag = true;
    1165           0 :         for (uInt sample_i=0;sample_i<width;sample_i++)
    1166             :         {
    1167           0 :                 if (not inputVector(startInputPos+sample_i))
    1168             :                 {
    1169           0 :                         outputFlag = false;
    1170           0 :                         break;
    1171             :                 }
    1172             :                 }
    1173           0 :         outputVector(outputPos) = outputFlag;
    1174             : 
    1175           0 :         return;
    1176             : }
    1177             : 
    1178             : //////////////////////////////////////////////////////////////////////////
    1179             : // ChannelAccumulationKernel class
    1180             : //   (numerical accumulation, respecting flags)
    1181             : //////////////////////////////////////////////////////////////////////////
    1182             : 
    1183             : // -----------------------------------------------------------------------
    1184             : //
    1185             : // -----------------------------------------------------------------------
    1186           0 : template<class T> void ChannelAccumulationKernel<T>::kernel(        DataCubeMap *inputData,
    1187             :                                                                                                                                 DataCubeMap *outputData,
    1188             :                                                                                                                                 uInt startInputPos,
    1189             :                                                                                                                                 uInt outputPos,
    1190             :                                                                                                                                 uInt width)
    1191             : {
    1192           0 :         T acc = 0;
    1193           0 :         uInt inputPos = 0;
    1194           0 :         Vector<T> &inputVector = inputData->getVector<T>(MS::DATA);
    1195           0 :         Vector<Bool> &inputFlagVector = inputData->getVector<Bool>(MS::FLAG);
    1196           0 :         Vector<T> &outputVector = outputData->getVector<T>(MS::DATA);
    1197           0 :         Bool accumulatorFlag = inputFlagVector(startInputPos);
    1198             : 
    1199           0 :         for (uInt sample_i=0;sample_i<width;sample_i++)
    1200             :         {
    1201             :                 // Get input index
    1202           0 :                 inputPos = startInputPos + sample_i;
    1203             : 
    1204             :                 // true/true or false/false
    1205           0 :                 if (accumulatorFlag == inputFlagVector(inputPos))
    1206             :                 {
    1207           0 :                         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           0 :         outputVector(outputPos) = acc;
    1219             : 
    1220           0 :         return;
    1221             : }
    1222             : 
    1223             : } //# NAMESPACE VI - END
    1224             : 
    1225             : } //# NAMESPACE CASA - END
    1226             : 
    1227             : 

Generated by: LCOV version 1.16