StatCollector.h

Go to the documentation of this file.

#ifndef LARDATA_UTILITIES_STATCOLLECTOR_H
#define LARDATA_UTILITIES_STATCOLLECTOR_H

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

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