StatCollector.h

Go to the documentation of this file.

#ifndef LARDATA_UTILITIES_STATCOLLECTOR_H
#define LARDATA_UTILITIES_STATCOLLECTOR_H

// C/C++ standard libraries
#include <cmath> // std::sqrt()
#include <array>
#include <tuple>
#include <limits> // std::numeric_limits<>
#include <initializer_list>
#include <iterator> // std::begin(), std::end()
#include <utility> // std::forward()
#include <algorithm> // std::for_each()
#include <stdexcept> // std::range_error


namespace lar {
  namespace util {
    
    
    struct identity {
      template<typename T>
      constexpr auto operator()(T&& v) const noexcept
        -> decltype(std::forward<T>(v))
      {
        return std::forward<T>(v);
      } // operator()
    }; // class identity
    
    
    namespace details {
      
      template <typename W>
      class WeightTracker {
          public:
        using Weight_t = W; 
        
        void add(Weight_t weight) { ++n; w += weight; }
        
        void clear() { n= 0; w = Weight_t(0); }
        
        int N() const { return n; }
        
        Weight_t Weights() const { return w; }
        
        Weight_t AverageWeight() const;
        
        
        template <typename V>
        static constexpr V sqr(V const& v) { return v*v; }
        
          protected:
        int n = 0;                
        Weight_t w = Weight_t(0); 
        
      }; // class WeightTracker<>
      
      
      template <unsigned int PWR, typename T, typename W = T>
      class DataTracker {
          public:
        static constexpr unsigned int Power = PWR;
        static_assert(Power >= 1, "DataTracker must have at least power 1");
        
        using Data_t = T; 
        using Weight_t = T; 
        
        DataTracker() { clear(); }
        
        void add(Data_t v, Weight_t w)
          {
            Weight_t x = w;
            for (size_t i = 0; i < Power; ++i) sums[i] += (x *= v);
          }
        
        void clear() { sums.fill(Data_t(0)); }
        
        template <unsigned int N>
        Weight_t SumN() const
          {
            static_assert
              ((N > 0) && (N <= Power), "Invalid sum of power requested");
            return sums[N - 1];
          } // SumN()
        
        Weight_t Sum(unsigned int n) const { return sums[n - 1]; }
        
        Weight_t Sum() const { return SumN<1>(); }
        
          protected:
        std::array<Weight_t, Power> sums;
        
      }; // class DataTracker<>
      
      
      template <typename T, typename W = T, unsigned int PWR = 2>
      class DataTracker2: public DataTracker<PWR, T, W> {
        using Base_t = DataTracker<PWR, T, W>; 
          public:
        using Base_t::Power;
        static_assert(Power >= 2, "DataTracker2 must have Power >= 2");
        using Weight_t = typename Base_t::Weight_t;
        
        Weight_t SumSq() const
          { return Base_t::sums[1]; }
        //{ return Base_t::SumN<2>(); } // I can't make gcc 4.9.1 digest this...
        
      }; // class DataTracker2
      
      template <typename T, typename W = T, unsigned int PWR = 3>
      class DataTracker3: public DataTracker2<T, W, PWR> {
        using Base_t = DataTracker2<T, W, PWR>; 
          public:
        using Base_t::Power;
        static_assert(Power >= 3, "DataTracker3 must have Power >= 3");
        using Weight_t = typename Base_t::Weight_t;
        
        Weight_t SumCube() const
          { return Base_t::sums[2]; }
        //{ return Base_t::SumN<3>(); } // I can't make gcc 4.9.1 digest this...
        
      }; // class DataTracker3
      
    } // namespace details
    
    
    
    template <typename T, typename W = T>
    class StatCollector: protected details::WeightTracker<W> {
      using Base_t = details::WeightTracker<W>;
      using Base_t::sqr;
        public:
      using This_t = StatCollector<T, W>; 
      using Data_t = T; 
      using Weight_t = typename Base_t::Weight_t; 
      
      // default constructor, destructor and all the rest
      
      
      void add(Data_t value, Weight_t weight = Weight_t(1.0));
      
      template <typename Iter>
      void add_unweighted(Iter begin, Iter end)
        {  add_unweighted(begin, end, identity()); }
      
      template <typename Iter, typename Pred>
      void add_unweighted(Iter begin, Iter end, Pred extractor);
      
      template <typename Cont, typename Pred>
      void add_unweighted(Cont cont, Pred extractor)
        { add_unweighted(std::begin(cont), std::end(cont), extractor); }
      
      template <typename Cont>
      void add_unweighted(Cont cont)
        { add_unweighted(std::begin(cont), std::end(cont)); }
      
      
      template <
        typename VIter, typename WIter,
        typename VPred, typename WPred = identity
        >
      void add_weighted(
        VIter begin_value, VIter end_value,
        WIter begin_weight,
        VPred value_extractor,
        WPred weight_extractor = WPred()
        );
      
      
      template <typename Iter>
      void add_weighted(Iter begin, Iter end);
      
      template <typename Cont>
      void add_weighted(Cont cont)
        { add_weighted(std::begin(cont), std::end(cont)); }
      
      
      void clear();
      
      
      int N() const { return Base_t::N(); }
      
      Weight_t Weights() const { return Base_t::Weights(); }
      
      Weight_t Sum() const { return x.Sum(); }
      
      Weight_t SumSq() const { return x.SumSq(); }
      
      Weight_t Average() const;
      
      Weight_t Variance() const;
      
      Weight_t RMS() const;
      
      
      Weight_t AverageWeight() const { return Base_t::AverageWeight(); }
      
      
      
        protected:
      using Variable_t = details::DataTracker2<Data_t, Weight_t>;
      
      Variable_t x; 
      
    }; // class StatCollector<>
    
    
    template <typename T, typename W = T>
    class StatCollector2D: public details::WeightTracker<W> {
        public:
      using Base_t = details::WeightTracker<W>;
      using Base_t::sqr;
        public:
      using This_t = StatCollector2D<T, W>; 
      using Data_t = T; 
      using Weight_t = typename Base_t::Weight_t; 

      using Pair_t = std::tuple<Data_t, Data_t>;
      using WeightedPair_t = std::tuple<Data_t, Data_t, Weight_t>;
      
      // default constructor, destructor and all the rest
      
      
      void add(Data_t x, Data_t y, Weight_t weight = Weight_t(1.0));
      
      void add(Pair_t value, Weight_t weight = Weight_t(1.0))
        { add(std::get<0>(value), std::get<1>(value), weight); }
      
      void add(WeightedPair_t value)
        { add(std::get<0>(value), std::get<1>(value), std::get<2>(value)); }
      
      template <typename Iter>
      void add_unweighted(Iter begin, Iter end)
        {  add_unweighted(begin, end, identity()); }
      
      template <typename Iter, typename Pred>
      void add_unweighted(Iter begin, Iter end, Pred extractor);
      
      template <typename Cont, typename Pred>
      void add_unweighted(Cont cont, Pred extractor)
        { add_unweighted(std::begin(cont), std::end(cont), extractor); }
      
      template <typename Cont>
      void add_unweighted(Cont cont)
        { add_unweighted(std::begin(cont), std::end(cont)); }
      
      
      template <
        typename VIter, typename WIter,
        typename VPred, typename WPred = identity
        >
      void add_weighted(
        VIter begin_value, VIter end_value,
        WIter begin_weight,
        VPred value_extractor,
        WPred weight_extractor = WPred()
        );
      
      
      template <typename Iter>
      void add_weighted(Iter begin, Iter end);
      
      template <typename Cont>
      void add_weighted(Cont cont)
        { add_weighted(std::begin(cont), std::end(cont)); }
      
      
      void clear();
      
      
      int N() const { return Base_t::N(); }
      
      Weight_t Weights() const { return Base_t::Weights(); }
      
      Weight_t SumX() const { return x.Sum(); }
      
      Weight_t SumY() const { return y.Sum(); }
      
      Weight_t SumSqX() const { return x.SumSq(); }
      
      Weight_t SumSqY() const { return y.SumSq(); }
      
      Weight_t SumXY() const { return sum_xy; }
      
      Weight_t AverageX() const;
      
      Weight_t AverageY() const;
      
      Weight_t VarianceX() const;
      
      Weight_t VarianceY() const;
      
      Weight_t Covariance() const;
      
      Weight_t RMSx() const;
      
      Weight_t RMSy() const;
      
      
      Weight_t LinearCorrelation() const;
      
      
      Weight_t AverageWeight() const { return Base_t::AverageWeight(); }
      
      
        protected:
      using Variable_t = details::DataTracker2<Data_t, Weight_t>;
      
      Variable_t x; 
      Variable_t y; 
      Weight_t sum_xy = Weight_t(0); 
      
    }; // class StatCollector2D<>
    
    
    
    template <typename T>
    class MinMaxCollector {
        public:
      using Data_t = T; 
      using This_t = MinMaxCollector<T>; 
      
      MinMaxCollector() = default;
      
      MinMaxCollector(std::initializer_list<Data_t> init)
        { add(init); }
      
      template <typename Iter>
      MinMaxCollector(Iter begin, Iter end)
        { add(begin, end); }
      
      
      // default copy and move constructor and assignment, and destructor
      

      This_t& add(Data_t value);
      
      This_t& add(std::initializer_list<Data_t> values);
      
      template <typename Iter>
      This_t& add(Iter begin, Iter end);
      
      
      bool has_data() const { return minimum <= maximum; }
      
      Data_t min() const { return minimum; }
      
      Data_t max() const { return maximum; }
      
      
      void clear();
      
        protected:
      Data_t minimum = std::numeric_limits<Data_t>::max();
      
      Data_t maximum = std::numeric_limits<Data_t>::min();
      
    }; // class MinMaxCollector<>
    
    
  } // namespace util
} // namespace lar


//==============================================================================
//=== template implementation
//==============================================================================


//******************************************************************************
//***  details::WeightTracker<>
//***

template <typename W>
typename lar::util::details::WeightTracker<W>::Weight_t
  lar::util::details::WeightTracker<W>::AverageWeight() const
{
  if (N() == 0)
    throw std::range_error("WeightTracker<>::AverageWeight(): divide by 0");
  return Weights() / N();
} // details::WeightTracker<W>::AverageWeight()


//******************************************************************************
//***  StatCollector<>
//***

template <typename T, typename W>
void lar::util::StatCollector<T, W>::add
  (Data_t value, Weight_t weight /* = Weight_t(1.0) */)
{
  Base_t::add(weight);
  x.add(value, weight);
} // StatCollector<T, W>::add()


template <typename T, typename W>
template <typename Iter, typename Pred>
void lar::util::StatCollector<T, W>::add_unweighted
  (Iter begin, Iter end, Pred extractor)
{
  std::for_each
    (begin, end, [this, extractor](auto item) { this->add(extractor(item)); });
} // StatCollector<T, W>::add_unweighted(Iter, Iter, Pred)


template <typename T, typename W>
template <typename VIter, typename WIter, typename VPred, typename WPred>
void lar::util::StatCollector<T, W>::add_weighted(
  VIter begin_value, VIter end_value,
  WIter begin_weight,
  VPred value_extractor,
  WPred weight_extractor /* = WPred() */
) {
  while (begin_value != end_value) {
    add(value_extractor(*begin_value), weight_extractor(*begin_weight));
    ++begin_value;
    ++begin_weight;
  } // while
} // StatCollector<T, W>::add_weighted(VIter, VIter, WIter, VPred, WPred)


template <typename T, typename W>
template <typename Iter>
void lar::util::StatCollector<T, W>::add_weighted(Iter begin, Iter end) {
  
  std::for_each(begin, end, [this](auto p) { this->add(p.first, p.second); });
  
} // StatCollector<T, W>::add_weighted(Iter, Iter)



template <typename T, typename W>
inline void lar::util::StatCollector<T, W>::clear() {
  Base_t::clear();
  x.clear();
} // StatCollector<T, W>::clear()


template <typename T, typename W>
typename lar::util::StatCollector<T, W>::Weight_t
  lar::util::StatCollector<T, W>::Average() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector<>::Average(): divide by 0");
  return Sum() / Weights();
} // StatCollector<T, W>::Average()


template <typename T, typename W>
typename lar::util::StatCollector<T, W>::Weight_t
  lar::util::StatCollector<T, W>::Variance() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector<>::Variance(): divide by 0");
  return std::max(Weight_t(0), (SumSq() - sqr(Sum()) / Weights()) / Weights());
} // StatCollector<T, W>::Variance()


template <typename T, typename W>
typename lar::util::StatCollector<T, W>::Weight_t
  lar::util::StatCollector<T, W>::RMS() const
{
  const Weight_t rms2 = Variance();
  if (rms2 < Weight_t(0)) return 0.;
//    throw std::range_error("StatCollector<>::RMS(): negative RMS^2");
  return std::sqrt(rms2);
} // StatCollector<T, W>::RMS()



//******************************************************************************
//***  StatCollector2D<>
//***

template <typename T, typename W>
void lar::util::StatCollector2D<T, W>::add
  (Data_t x_value, Data_t y_value, Weight_t weight /* = Weight_t(1.0) */)
{
  Base_t::add(weight);
  x.add(x_value, weight);
  y.add(y_value, weight);
  sum_xy += weight * x_value * y_value;
} // StatCollector2D<T, W>::add()


template <typename T, typename W>
template <typename Iter, typename Pred>
void lar::util::StatCollector2D<T, W>::add_unweighted
  (Iter begin, Iter end, Pred extractor)
{
  std::for_each
    (begin, end, [this, extractor](auto item) { this->add(extractor(item)); });
} // StatCollector2D<T, W>::add_unweighted(Iter, Iter, Pred)


template <typename T, typename W>
template <typename VIter, typename WIter, typename VPred, typename WPred>
void lar::util::StatCollector2D<T, W>::add_weighted(
  VIter begin_value, VIter end_value,
  WIter begin_weight,
  VPred value_extractor,
  WPred weight_extractor /* = WPred() */
) {
  while (begin_value != end_value) {
    add(value_extractor(*begin_value), weight_extractor(*begin_weight));
    ++begin_value;
    ++begin_weight;
  } // while
} // StatCollector2D<T, W>::add_weighted(VIter, VIter, WIter, VPred, WPred)


template <typename T, typename W>
template <typename Iter>
void lar::util::StatCollector2D<T, W>::add_weighted(Iter begin, Iter end) {
  
  std::for_each(begin, end, [this](auto p) { this->add(p); });
  
} // StatCollector2D<T, W>::add_weighted(Iter, Iter)



template <typename T, typename W>
inline void lar::util::StatCollector2D<T, W>::clear() {
  Base_t::clear();
  x.clear();
  y.clear();
  sum_xy = Weight_t(0);
} // StatCollector<T, W>::clear()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::AverageX() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector2D<>::AverageX(): divide by 0");
  return SumX() / Weights();
} // StatCollector2D<T, W>::AverageX()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::VarianceX() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector2D<>::VarianceX(): divide by 0");
  return (SumSqX() - sqr(SumX()) / Weights()) / Weights();
} // StatCollector2D<T, W>::VarianceX()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::RMSx() const
{
  const Weight_t rms2 = VarianceX();
  if (rms2 < Weight_t(0)) return 0.;
//    throw std::range_error("StatCollector2D<>::RMSx(): negative RMS^2");
  return std::sqrt(rms2);
} // StatCollector2D<T, W>::RMSx()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::AverageY() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector2D<>::AverageY(): divide by 0");
  return SumY() / Weights();
} // StatCollector2D<T, W>::AverageY()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::VarianceY() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector2D<>::VarianceY(): divide by 0");
  return (SumSqY() - sqr(SumY()) / Weights()) / Weights();
} // StatCollector2D<T, W>::VarianceY()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::Covariance() const
{
  if (Weights() == Weight_t(0))
    throw std::range_error("StatCollector2D<>::Covariance(): divide by 0");
  return (SumXY() - SumX() * SumY() / Weights()) / Weights();
} // StatCollector2D<T, W>::VarianceY()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::RMSy() const
{
  const Weight_t rms2 = VarianceY();
  if (rms2 < Weight_t(0)) return 0.;
//    throw std::range_error("StatCollector2D<>::RMSy(): negative RMS^2");
  return std::sqrt(rms2);
} // StatCollector2D<T, W>::RMSy()


template <typename T, typename W>
typename lar::util::StatCollector2D<T, W>::Weight_t
  lar::util::StatCollector2D<T, W>::LinearCorrelation() const
{
  if (Weights() == Data_t(0))
    throw std::range_error("StatCollector2D<>::LinearCorrelation(): divide by 0");
  
  const Weight_t var_prod = VarianceX() * VarianceY();
  if (var_prod <= Weight_t(0))
    throw std::range_error("StatCollector2D<>::LinearCorrelation(): variance is 0");
  
  return Covariance() / std::sqrt(var_prod);
} // StatCollector2D<T, W>::LinearCorrelation()



//******************************************************************************
//*** MinMaxCollector
//***

template <typename T>
lar::util::MinMaxCollector<T>& lar::util::MinMaxCollector<T>::add(Data_t value)
{
  if (value < minimum) minimum = value;
  if (value > maximum) maximum = value;
  return *this;
} // lar::util::MinMaxCollector<T>::add


template <typename T>
inline lar::util::MinMaxCollector<T>& lar::util::MinMaxCollector<T>::add
  (std::initializer_list<Data_t> values)
  { return add(values.begin(), values.end()); }


template <typename T> template <typename Iter>
inline lar::util::MinMaxCollector<T>& lar::util::MinMaxCollector<T>::add
  (Iter begin, Iter end)
{
  std::for_each(begin, end, [this](Data_t value) { this->add(value); });
  return *this;
} // lar::util::MinMaxCollector<T>::add(Iter)


template <typename T>
inline void lar::util::MinMaxCollector<T>::clear() {
  minimum = std::numeric_limits<Data_t>::max();
  maximum = std::numeric_limits<Data_t>::min();
} // lar::util::MinMaxCollector<T>::clear()

//******************************************************************************


#endif // LARDATA_UTILITIES_STATCOLLECTOR_H