Line data    Source code

       1             : //#  CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
       2             : //#  Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All
       3             : //#  rights reserved.
       4             : //#  Copyright (C) European Southern Observatory, 2011, All rights reserved.
       5             : //#
       6             : //#  This library is free software; you can redistribute it and/or
       7             : //#  modify it under the terms of the GNU Lesser General Public
       8             : //#  License as published by the Free software Foundation; either
       9             : //#  version 2.1 of the License, or (at your option) any later version.
      10             : //#
      11             : //#  This library is distributed in the hope that it will be useful,
      12             : //#  but WITHOUT ANY WARRANTY, without even the implied warranty of
      13             : //#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      14             : //#  Lesser General Public License for more details.
      15             : //#
      16             : //#  You should have received a copy of the GNU Lesser General Public
      17             : //#  License along with this library; if not, write to the Free Software
      18             : //#  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
      19             : //#  MA 02111-1307  USA
      20             : 
      21             : #include <mstransform/TVI/StatWtVarianceAndWeightCalculator.h>
      22             : 
      23             : #include <casacore/casa/Arrays/ArrayMath.h>
      24             : #include <casacore/casa/Arrays/Cube.h>
      25             : 
      26             : // DEBUG
      27             : #include <casacore/casa/IO/ArrayIO.h>
      28             : 
      29             : #ifdef _OPENMP
      30             : #include <omp.h>
      31             : #endif
      32             : 
      33             : using namespace casacore;
      34             : using namespace std;
      35             : 
      36             : namespace casa {
      37             : 
      38          46 : StatWtVarianceAndWeightCalculator::StatWtVarianceAndWeightCalculator(
      39             :     const shared_ptr<
      40             :         StatisticsAlgorithm<
      41             :             Double, Array<Float>::const_iterator, Array<Bool>::const_iterator,
      42             :             Array<Double>::const_iterator
      43             :         >
      44             :     > statAlg,
      45             :     shared_ptr<map<uInt, pair<uInt, uInt>>> samples, Int minSamp
      46          46 : ) : _statAlg(statAlg->clone()), _samples(samples), _minSamp(minSamp) {}
      47             : 
      48          46 : StatWtVarianceAndWeightCalculator::~StatWtVarianceAndWeightCalculator() {}
      49             : 
      50      433691 : Double StatWtVarianceAndWeightCalculator::computeVariance(
      51             :     const Cube<Complex>& data,
      52             :     const Cube<Bool>& flags, const Vector<Double>& exposures,
      53             :     casacore::uInt spw
      54             : ) const {
      55      433691 :     const auto npts = data.size();
      56      433691 :     if ((Int)npts < _minSamp || (Int)nfalse(flags) < _minSamp) {
      57             :         // not enough points, trivial
      58       77287 :         return 0;
      59             :     }
      60             :     // called in multi-threaded mode, so need to clone this for each thread
      61             :     std::unique_ptr<
      62             :         StatisticsAlgorithm<
      63             :             Double, Array<Float>::const_iterator,
      64             :             Array<Bool>::const_iterator, Array<Double>::const_iterator
      65             :         >
      66      356404 :     > statAlg(_statAlg->clone());
      67             :     // some data not flagged
      68      356404 :     const auto realPart = real(data);
      69      356404 :     const auto imagPart = imag(data);
      70      356404 :     const auto mask = ! flags;
      71      356404 :     Cube<Double> exposureCube(data.shape());
      72      356404 :     const auto nPlanes = data.nplane();
      73      990612 :     for (size_t i=0; i<nPlanes; ++i) {
      74      634208 :         exposureCube.xyPlane(i) = exposures[i];
      75             :     }
      76      356404 :     auto riter = realPart.begin();
      77      356404 :     auto iiter = imagPart.begin();
      78      356404 :     auto miter = mask.begin();
      79      356404 :     auto eiter = exposureCube.begin();
      80      356404 :     statAlg->setData(riter, eiter, miter, npts);
      81      356404 :     auto realStats = statAlg->getStatistics();
      82      356404 :     auto realVar = realStats.nvariance/realStats.npts;
      83             :     // reset data to imaginary parts
      84      356404 :     statAlg->setData(iiter, eiter, miter, npts);
      85      356404 :     auto imagStats = statAlg->getStatistics();
      86      356404 :     auto imagVar = imagStats.nvariance/imagStats.npts;
      87      356404 :     auto varSum = realVar + imagVar;
      88             :     // _samples.second can be updated in two different places, so use
      89             :     // a local (per thread) variable and update the object's private field in one
      90             :     // place
      91      356404 :     uInt updateSecond = False;
      92      356404 :     if (varSum > 0) {
      93             : #ifdef _OPENMP
      94      356404 : #pragma omp atomic
      95             : #endif
      96      356404 :         ++((*_samples)[spw].first);
      97      356404 :         if (imagVar == 0 || realVar == 0) {
      98           0 :             updateSecond = True;
      99             :         }
     100             :         else {
     101      356404 :             auto ratio = imagVar/realVar;
     102      356404 :             auto inverse = 1/ratio;
     103      356404 :             updateSecond = ratio > 1.5 || inverse > 1.5;
     104             :         }
     105      356404 :         if (updateSecond) {
     106             : #ifdef _OPENMP
     107      217625 : #pragma omp atomic
     108             : #endif
     109      217625 :             ++((*_samples)[spw].second);
     110             :         }
     111             :     }
     112      356404 :     return varSum/2;
     113      356404 : }
     114             : 
     115      427529 : Double StatWtVarianceAndWeightCalculator::computeWeight(
     116             :     const Cube<Complex>& data, const Cube<Bool>& flags,
     117             :     const Vector<Double>& exposures, uInt spw, Double targetExposure
     118             : ) const {
     119      427529 :     auto varEq = computeVariance(data, flags, exposures, spw);
     120      427529 :     return varEq == 0 ? 0 : targetExposure/varEq;
     121             : }
     122             : 
     123           0 : Vector<Double> StatWtVarianceAndWeightCalculator::computeWeights(
     124             :     const Cube<Complex>& data, const Cube<Bool>& flags,
     125             :     const Vector<Double>& exposures, uInt spw
     126             : ) const {
     127           0 :     auto varEq = computeVariance(data, flags, exposures, spw);
     128           0 :     return varEq == 0 ? Vector<Double>(exposures.size(), 0) : exposures/varEq;
     129             : }
     130             : 
     131             : }