Line data    Source code

       1             : //# UtilsTVI.h: This file contains the interface definition of the MSTransformManager 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/UtilsTVI.h>
      24             : 
      25             : using namespace casacore;
      26             : namespace casa { //# NAMESPACE CASA - BEGIN
      27             : 
      28             : namespace vi { //# NAMESPACE VI - BEGIN
      29             : 
      30             : //////////////////////////////////////////////////////////////////////////
      31             : // DataCubeHolderBase class
      32             : //////////////////////////////////////////////////////////////////////////
      33             : 
      34             : // -----------------------------------------------------------------------
      35             : //
      36             : // -----------------------------------------------------------------------
      37           0 : uInt DataCubeHolderBase::getMatrixIndex()
      38             : {
      39           0 :         return matrixIndex_p;
      40             : }
      41             : 
      42             : // -----------------------------------------------------------------------
      43             : //
      44             : // -----------------------------------------------------------------------
      45           0 : uInt DataCubeHolderBase::getVectorIndex()
      46             : {
      47           0 :         return vectorIndex_p;
      48             : }
      49             : 
      50             : // -----------------------------------------------------------------------
      51             : //
      52             : // -----------------------------------------------------------------------
      53           0 : IPosition & DataCubeHolderBase::getCubeShape()
      54             : {
      55           0 :         return cubeShape_p;
      56             : }
      57             : 
      58             : // -----------------------------------------------------------------------
      59             : //
      60             : // -----------------------------------------------------------------------
      61           0 : IPosition & DataCubeHolderBase::getMatrixShape()
      62             : {
      63           0 :         return matrixShape_p;
      64             : }
      65             : 
      66             : // -----------------------------------------------------------------------
      67             : //
      68             : // -----------------------------------------------------------------------
      69           0 : IPosition & DataCubeHolderBase::getVectorShape()
      70             : {
      71           0 :         return vectorShape_p;
      72             : }
      73             : 
      74             : //////////////////////////////////////////////////////////////////////////
      75             : // DataCubeMap class
      76             : //////////////////////////////////////////////////////////////////////////
      77             : 
      78             : // -----------------------------------------------------------------------
      79             : //
      80             : // -----------------------------------------------------------------------
      81           0 : DataCubeMap::DataCubeMap()
      82             : {
      83           0 :         dataCubeMap_p.clear();
      84           0 : }
      85             : 
      86             : // -----------------------------------------------------------------------
      87             : //
      88             : // -----------------------------------------------------------------------
      89           0 : DataCubeMap::DataCubeMap(DataCubeMap& other)
      90             : {
      91           0 :         dataCubeMap_p.clear();
      92             : 
      93           0 :         for (dataCubeMapIter_p = other.dataCubeMap_p.begin();dataCubeMapIter_p!= other.dataCubeMap_p.end();dataCubeMapIter_p++)
      94             :         {
      95           0 :                 add(dataCubeMapIter_p->first,dataCubeMapIter_p->second->selfReference());
      96             :         }
      97             : 
      98           0 :         return;
      99             : }
     100             : 
     101             : // -----------------------------------------------------------------------
     102             : //
     103             : // -----------------------------------------------------------------------
     104           0 : DataCubeMap::~DataCubeMap()
     105             : {
     106           0 :         dataCubeMap_p.clear();
     107           0 : }
     108             : 
     109             : // -----------------------------------------------------------------------
     110             : //
     111             : // -----------------------------------------------------------------------
     112           0 : void DataCubeMap::add(MS::PredefinedColumns key,DataCubeHolderBase* dataCubeHolder)
     113             : {
     114           0 :         dataCubeMap_p[key] = dataCubeHolder;
     115           0 : }
     116             : 
     117             : // -----------------------------------------------------------------------
     118             : //
     119             : // -----------------------------------------------------------------------
     120           0 : void DataCubeMap::add(MS::PredefinedColumns key,DataCubeHolderBase &dataCubeHolder)
     121             : {
     122           0 :         dataCubeMap_p[key] = &dataCubeHolder;
     123           0 : }
     124             : 
     125             : // -----------------------------------------------------------------------
     126             : //
     127             : // -----------------------------------------------------------------------
     128           0 : Bool DataCubeMap::present(MS::PredefinedColumns key)
     129             : {
     130           0 :         if (dataCubeMap_p.find(key) != dataCubeMap_p.end())
     131             :         {
     132           0 :                 return true;
     133             :         }
     134             :         else
     135             :         {
     136           0 :                 return false;
     137             :         }
     138             : }
     139             : 
     140             : // -----------------------------------------------------------------------
     141             : //
     142             : // -----------------------------------------------------------------------
     143           0 : void DataCubeMap::setMatrixIndex(uInt rowIndex)
     144             : {
     145           0 :         for (dataCubeMapIter_p = dataCubeMap_p.begin();dataCubeMapIter_p!= dataCubeMap_p.end();dataCubeMapIter_p++)
     146             :         {
     147           0 :                 dataCubeMapIter_p->second->setMatrixIndex(rowIndex);
     148             :         }
     149             : 
     150           0 :         return;
     151             : }
     152             : 
     153             : // -----------------------------------------------------------------------
     154             : //
     155             : // -----------------------------------------------------------------------
     156           0 : void DataCubeMap::setVectorIndex(uInt vectorIndex)
     157             : {
     158           0 :         for (dataCubeMapIter_p = dataCubeMap_p.begin();dataCubeMapIter_p!= dataCubeMap_p.end();dataCubeMapIter_p++)
     159             :         {
     160           0 :                 dataCubeMapIter_p->second->setVectorIndex(vectorIndex);
     161             :         }
     162             : 
     163           0 :         return;
     164             : }
     165             : 
     166             : // -----------------------------------------------------------------------
     167             : //
     168             : // -----------------------------------------------------------------------
     169           0 : IPosition & DataCubeMap::getCubeShape()
     170             : {
     171           0 :         return dataCubeMap_p.begin()->second->getCubeShape();
     172             : }
     173             : 
     174             : // -----------------------------------------------------------------------
     175             : //
     176             : // -----------------------------------------------------------------------
     177           0 : IPosition & DataCubeMap::getMatrixShape()
     178             : {
     179           0 :         return dataCubeMap_p.begin()->second->getMatrixShape();
     180             : }
     181             : 
     182             : // -----------------------------------------------------------------------
     183             : //
     184             : // -----------------------------------------------------------------------
     185           0 : IPosition & DataCubeMap::getVectorShape()
     186             : {
     187           0 :         return dataCubeMap_p.begin()->second->getVectorShape();
     188             : }
     189             : 
     190             : // -----------------------------------------------------------------------
     191             : //
     192             : // -----------------------------------------------------------------------
     193           0 : size_t DataCubeMap::nelements()
     194             : {
     195           0 :         return dataCubeMap_p.size();
     196             : }
     197             : 
     198             : 
     199             : // Methods controlling iteration
     200             : // gmoellen (2017Mar06)
     201           0 : void DataCubeMap::setupVecIter() 
     202             : {
     203           0 :   for (dataCubeMapIter_p = dataCubeMap_p.begin();
     204           0 :        dataCubeMapIter_p!= dataCubeMap_p.end();
     205           0 :        dataCubeMapIter_p++)
     206             :     {
     207           0 :       dataCubeMapIter_p->second->setupVecIter();
     208             :     }
     209           0 :   return;
     210             : }
     211           0 : void DataCubeMap::reset() 
     212             : {
     213           0 :   for (dataCubeMapIter_p = dataCubeMap_p.begin();
     214           0 :        dataCubeMapIter_p!= dataCubeMap_p.end();
     215           0 :        dataCubeMapIter_p++)
     216             :     {
     217           0 :       dataCubeMapIter_p->second->reset();
     218             :     }
     219           0 :   return;
     220             : }
     221             : 
     222           0 : void DataCubeMap::next()
     223             : {
     224           0 :   for (dataCubeMapIter_p = dataCubeMap_p.begin();
     225           0 :        dataCubeMapIter_p!= dataCubeMap_p.end();
     226           0 :        dataCubeMapIter_p++)
     227             :     {
     228           0 :       dataCubeMapIter_p->second->next();
     229             :     }
     230           0 :   return;
     231             : }
     232             : 
     233           0 : casacore::Bool DataCubeMap::pastEnd()
     234             : {
     235             :   // All elements in the map are sync'd, so just ask the first one
     236           0 :   return dataCubeMap_p.begin()->second->pastEnd();
     237             : }
     238             : 
     239             : 
     240             : //////////////////////////////////////////////////////////////////////////
     241             : // Convenience methods
     242             : //////////////////////////////////////////////////////////////////////////
     243             : 
     244             : // -----------------------------------------------------------------------
     245             : //
     246             : // -----------------------------------------------------------------------
     247           0 : void accumulateWeightCube (const Cube<Float> &weightCube, const Cube<Bool> &flags, Matrix<Float> &result)
     248             : {
     249           0 :         IPosition shape = weightCube.shape();
     250           0 :         uInt nRows = shape(2);
     251           0 :         uInt nCorrelations = shape (0);
     252             : 
     253           0 :         result.resize(nCorrelations,nRows,false);
     254             : 
     255           0 :         Vector<Float> currentVector;
     256           0 :         for (uInt row=0; row < nRows; row++)
     257             :         {
     258           0 :                 currentVector.reference(result.column(row));
     259           0 :                 accumulateWeightMatrix (weightCube.xyPlane(row), flags.xyPlane(row),currentVector);
     260             :         }
     261             : 
     262           0 :     return;
     263             : }
     264             : 
     265             : // -----------------------------------------------------------------------
     266             : //
     267             : // -----------------------------------------------------------------------
     268           0 : void accumulateWeightMatrix (const Matrix<Float> &weightMatrix, const Matrix<Bool> &flags, Vector<Float> &result)
     269             : {
     270           0 :     IPosition shape = weightMatrix.shape();
     271           0 :     result.resize(shape(0),false);
     272           0 :     Vector<uInt> samples(shape(0),0);
     273           0 :     uInt nCorrelations = shape (0);
     274           0 :     uInt nChannels = shape (1);
     275             : 
     276           0 :     for (uInt correlation = 0; correlation < nCorrelations; correlation++)
     277             :     {
     278           0 :         int nSamples = 0;
     279           0 :         float sum = 0;
     280           0 :         bool accumulatorFlag = true;
     281             : 
     282           0 :         for (uInt channel=0; channel< nChannels; channel++)
     283             :         {
     284           0 :             Bool inputFlag = flags(correlation,channel);
     285             :             // true/true or false/false
     286           0 :             if (accumulatorFlag == inputFlag)
     287             :             {
     288           0 :                 nSamples ++;
     289           0 :                 sum += weightMatrix (correlation, channel);
     290             :             }
     291             :             // true/false: Reset accumulation when accumulator switches from flagged to unflagged
     292           0 :             else if ( accumulatorFlag and ! inputFlag )
     293             :             {
     294           0 :                 accumulatorFlag = false;
     295           0 :                 nSamples = 1;
     296           0 :                 sum = weightMatrix (correlation, channel);
     297             :             }
     298             :             // else ignore case where accumulator is valid and data is not
     299             : 
     300             :         }
     301             : 
     302           0 :         result (correlation) = sum / (nSamples != 0 ? nSamples : 1);
     303             :     }
     304             : 
     305           0 :     return;
     306             : }
     307             : 
     308             : // -----------------------------------------------------------------------
     309             : //
     310             : // -----------------------------------------------------------------------
     311           0 : void accumulateFlagCube (const Cube<Bool> &flagCube, Vector<Bool> &flagRow)
     312             : {
     313             :         // Get original shape
     314           0 :         IPosition shape = flagCube.shape();
     315           0 :         size_t nCorr = shape(0);
     316           0 :         size_t nChan = shape(1);
     317           0 :         size_t nRows = shape(2);
     318             : 
     319             :         // Reshape flag cube to match the input shape
     320           0 :         flagRow.resize(nRows,false);
     321           0 :         flagRow = false;
     322             : 
     323           0 :         Bool rowFlagValue = false;
     324           0 :         for (size_t row_i =0;row_i<nRows;row_i++)
     325             :         {
     326           0 :                 rowFlagValue = true;
     327           0 :                 for (size_t chan_i =0;chan_i<nChan;chan_i++)
     328             :                 {
     329           0 :                         if (rowFlagValue)
     330             :                         {
     331           0 :                                 for (size_t corr_i =0;corr_i<nCorr;corr_i++)
     332             :                                 {
     333           0 :                                         if (not flagCube(corr_i,chan_i,row_i))
     334             :                                         {
     335           0 :                                                 rowFlagValue = false;
     336           0 :                                                 break;
     337             :                                         }
     338             :                                 }
     339             :                         }
     340             :                         else
     341             :                         {
     342           0 :                                 break;
     343             :                         }
     344             :                 }
     345           0 :                 flagRow(row_i) = rowFlagValue;
     346             :         }
     347             : 
     348           0 :         return;
     349             : }
     350             : 
     351             : 
     352             : } //# NAMESPACE VI - END
     353             : 
     354             : } //# NAMESPACE CASA - END
     355             :