StatisticsUtilities.h

Go to the documentation of this file.

// # Copyright (C) 2000,2001
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library is distributed in the hope that it will be useful, but WITHOUT
// # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
// #
 
#ifndef SCIMATH_STATISTICSUTILITIES_H
#define SCIMATH_STATISTICSUTILITIES_H
 
#include <casacore/casa/Exceptions/Error.h>
#include <casacore/scimath/StatsFramework/StatisticsTypes.h>
#include <casacore/scimath/StatsFramework/StatsHistogram.h>
#include <casacore/casa/Utilities/DataType.h>
#include <casacore/casa/aips.h>
 
#include <iostream>
#include <casacore/casa/iosfwd.h>
 
namespace casacore {
 
CASA_STATD class StatsDataProvider;
 
// Various statistics related methods for the statistics framework.
 
template <class AccumType>
class StatisticsUtilities {
 public:
  StatisticsUtilities() = delete;
 
  ~StatisticsUtilities() {}
 
  // <group>
  // accumulate values. It is the responsibility of the caller to keep track
  // of the accumulated values after each call. This class does not since it
  // has no state. The accumulation derivation for mean and variance can be
  // found at
  // www.itl.nist.gov/div898/software/dataplot/refman2/ch2/weighvar.pdf
  // nvariance is an accumulated value. It is related to the variance via
  // variance = nvariance/npts or nvariance/(npts-1) depending on your
  // preferred definition in the non-weighted case and
  // wvariance = wnvariance/sumofweights or wnvariance/(sumofweights-1)
  // in the weighted case Its basic definition is
  // nvariance = sum((x_i - mean)**2),
  // wnvariance = sum((weight_i*(x_i - mean)**2)
  // npts is a Double rather than an Int64 because of compilation issues when
  // T is a Complex
  inline static void accumulate(Double& npts, AccumType& sum, AccumType& mean,
                                const AccumType& datum);
 
  // in order to optimize performance, no checking is done for the
  // weight == 0 case callers should ensure that the weigth is not zero before
  // calling this method, and shouldn't call this method if the weight is 0.
  // Expect a segfault because of division by zero if sumweights and weight
  // are both zero.
  inline static void waccumulate(Double& npts, AccumType& sumweights, AccumType& wsum,
                                 AccumType& wmean, const AccumType& datum, const AccumType& weight);
 
  inline static void accumulate(Double& npts, AccumType& sum, AccumType& mean, AccumType& nvariance,
                                AccumType& sumsq, const AccumType& datum);
 
  // wsumsq is the weighted sum of squares, sum(w_i*x_i*x_i)
  inline static void waccumulate(Double& npts, AccumType& sumweights, AccumType& wsum,
                                 AccumType& wmean, AccumType& wnvariance, AccumType& wsumsq,
                                 const AccumType& datum, const AccumType& weight);
  // </group>
 
  // <group>
  // The assignment operator of class LocationType should use copy, not
  // reference, semantics.
  template <class LocationType>
  inline static void accumulate(Double& npts, AccumType& sum, AccumType& mean, AccumType& nvariance,
                                AccumType& sumsq, AccumType& datamin, AccumType& datamax,
                                LocationType& minpos, LocationType& maxpos, const AccumType& datum,
                                const LocationType& location);
 
  template <class LocationType, class DataType>
  inline static void accumulate(Double& npts, AccumType& sum, AccumType& mean, AccumType& nvariance,
                                AccumType& sumsq, DataType& datamin, DataType& datamax,
                                LocationType& minpos, LocationType& maxpos, const DataType& datum,
                                const LocationType& location);
 
  template <class LocationType>
  inline static void waccumulate(Double& npts, AccumType& sumofweights, AccumType& sum,
                                 AccumType& mean, AccumType& nvariance, AccumType& sumsq,
                                 AccumType& datamin, AccumType& datamax, LocationType& minpos,
                                 LocationType& maxpos, const AccumType& datum,
                                 const AccumType& weight, const LocationType& location);
  // </group>
 
  // <group>
  // return True if the max or min was updated, False otherwise.
  template <class LocationType>
  inline static Bool doMax(AccumType& datamax, LocationType& maxpos, Bool isFirst,
                           const AccumType& datum, const LocationType& location);
 
  template <class LocationType>
  inline static Bool doMin(AccumType& datamin, LocationType& minpos, Bool isFirst,
                           const AccumType& datum, const LocationType& location);
  // </group>
 
  // <group>
  // These versions are for symmetric accumulation about a specified center
  // point. The actual point is accumulated, as is a "virtual" point that is
  // symmetric about the specified center. Of course, the trivial relationship
  // that the mean is the specified center is used to simplify things
  inline static void accumulateSym(Double& npts, AccumType& nvariance, AccumType& sumsq,
                                   const AccumType& datum, const AccumType& center);
 
  // wsumsq is the weighted sum of squares, sum(w_i*x_i*x_i)
  inline static void waccumulateSym(Double& npts, AccumType& sumweights, AccumType& wnvariance,
                                    AccumType& wsumsq, const AccumType& datum,
                                    const AccumType& weight, const AccumType& center);
 
  // <src>maxpos</src> and <src>minpos</src> refer to actual, not
  // virtually created, data only.
  template <class LocationType>
  inline static void accumulateSym(Double& npts, AccumType& nvariance, AccumType& sumsq,
                                   AccumType& datamin, AccumType& datamax, LocationType& minpos,
                                   LocationType& maxpos, const AccumType& datum,
                                   const LocationType& location, const AccumType& center);
 
  template <class LocationType>
  inline static void waccumulateSym(Double& npts, AccumType& sumofweights, AccumType& nvariance,
                                    AccumType& sumsq, AccumType& datamin, AccumType& datamax,
                                    LocationType& minpos, LocationType& maxpos,
                                    const AccumType& datum, const AccumType& weight,
                                    const LocationType& location, const AccumType& center);
 
  // convert in place by taking the absolute value of the difference of the
  // std::vector and the median
  inline static void convertToAbsDevMedArray(DataArray& myArray, AccumType median);
  // </group>
 
  inline static Bool includeDatum(const AccumType& datum,
                                  typename DataRanges::const_iterator beginRange,
                                  typename DataRanges::const_iterator endRange, Bool isInclude);
 
  // The array can be changed by partially sorting it up to the largest index.
  // Return a map of index to value in the sorted array.
  static std::map<uInt64, AccumType> indicesToValues(std::vector<AccumType>& myArray,
                                                     const std::set<uInt64>& indices);
 
  static void mergeResults(
      std::vector<BinCountArray>& bins, std::vector<std::shared_ptr<AccumType>>& sameVal,
      std::vector<Bool>& allSame, const std::unique_ptr<std::vector<BinCountArray>[]>& tBins,
      const std::unique_ptr<std::vector<std::shared_ptr<AccumType>>[]>& tSameVal,
      const std::unique_ptr<std::vector<Bool>[]>& tAllSame, uInt nThreadsMax);
 
  // use two statistics sets to get the statistics set that would
  // result in combining the two data sets used to produce the
  // individual statistics sets. The quantile related stats are
  // not considered, since it is not in general possible to determine
  // the resultant quantiles from the information provided; only
  // the aggregate statistics make sense.
  static StatsData<AccumType> combine(const std::vector<StatsData<AccumType>>& stats);
 
  template <class DataIterator, class MaskIterator, class WeightsIterator>
  static uInt nThreadsMax(const StatsDataProvider<CASA_STATP>* const dataProvider);
 
  static uInt threadIdx();
 
 private:
  const static AccumType TWO;
};
 
}  // namespace casacore
 
#ifndef CASACORE_NO_AUTO_TEMPLATES
#include <casacore/scimath/StatsFramework/StatisticsUtilities.tcc>
#endif  // # CASACORE_NO_AUTO_TEMPLATES
 
#endif