boost/histogram/accumulators/mean.hpp

   1 // Copyright 2015-2018 Hans Dembinski
   2 //
   3 // Distributed under the Boost Software License, version 1.0.
   4 // (See accompanying file LICENSE_1_0.txt
   5 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
   8 #define BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
   9
  10 #include <boost/assert.hpp>
  11 #include <boost/core/nvp.hpp>
  12 #include <boost/histogram/fwd.hpp> // for mean<>
  13 #include <boost/throw_exception.hpp>
  14 #include <stdexcept>
  15 #include <type_traits>
  16
  17 namespace boost {
  18 namespace histogram {
  19 namespace accumulators {
  20
  21 /** Calculates mean and variance of sample.
  22
  23   Uses Welfords's incremental algorithm to improve the numerical
  24   stability of mean and variance computation.
  25 */
  26 template <class RealType>
  27 class mean {
  28 public:
  29   mean() = default;
  30   mean(const RealType& n, const RealType& mean, const RealType& variance) noexcept
  31       : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {}
  32
  33   void operator()(const RealType& x) noexcept {
  34     sum_ += static_cast<RealType>(1);
  35     const auto delta = x - mean_;
  36     mean_ += delta / sum_;
  37     sum_of_deltas_squared_ += delta * (x - mean_);
  38   }
  39
  40   void operator()(const weight_type<RealType>& w, const RealType& x) noexcept {
  41     sum_ += w.value;
  42     const auto delta = x - mean_;
  43     mean_ += w.value * delta / sum_;
  44     sum_of_deltas_squared_ += w.value * delta * (x - mean_);
  45   }
  46
  47   template <class T>
  48   mean& operator+=(const mean<T>& rhs) noexcept {
  49     if (sum_ != 0 || rhs.sum_ != 0) {
  50       const auto tmp = mean_ * sum_ + static_cast<RealType>(rhs.mean_ * rhs.sum_);
  51       sum_ += rhs.sum_;
  52       mean_ = tmp / sum_;
  53     }
  54     sum_of_deltas_squared_ += static_cast<RealType>(rhs.sum_of_deltas_squared_);
  55     return *this;
  56   }
  57
  58   mean& operator*=(const RealType& s) noexcept {
  59     mean_ *= s;
  60     sum_of_deltas_squared_ *= s * s;
  61     return *this;
  62   }
  63
  64   template <class T>
  65   bool operator==(const mean<T>& rhs) const noexcept {
  66     return sum_ == rhs.sum_ && mean_ == rhs.mean_ &&
  67            sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_;
  68   }
  69
  70   template <class T>
  71   bool operator!=(const mean<T>& rhs) const noexcept {
  72     return !operator==(rhs);
  73   }
  74
  75   const RealType& count() const noexcept { return sum_; }
  76   const RealType& value() const noexcept { return mean_; }
  77   RealType variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
  78
  79   template <class Archive>
  80   void serialize(Archive& ar, unsigned version) {
  81     if (version == 0) {
  82       // read only
  83       std::size_t sum;
  84       ar& make_nvp("sum", sum);
  85       sum_ = static_cast<RealType>(sum);
  86     } else {
  87       ar& make_nvp("sum", sum_);
  88     }
  89     ar& make_nvp("mean", mean_);
  90     ar& make_nvp("sum_of_deltas_squared", sum_of_deltas_squared_);
  91   }
  92
  93 private:
  94   RealType sum_ = 0, mean_ = 0, sum_of_deltas_squared_ = 0;
  95 };
  96
  97 } // namespace accumulators
  98 } // namespace histogram
  99 } // namespace boost
 100
 101 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
 102
 103 namespace boost {
 104 namespace serialization {
 105
 106 template <class T>
 107 struct version;
 108
 109 // version 1 for boost::histogram::accumulators::mean<RealType>
 110 template <class RealType>
 111 struct version<boost::histogram::accumulators::mean<RealType>>
 112     : std::integral_constant<int, 1> {};
 113
 114 } // namespace serialization
 115 } // namespace boost
 116
 117 namespace std {
 118 template <class T, class U>
 119 /// Specialization for boost::histogram::accumulators::mean.
 120 struct common_type<boost::histogram::accumulators::mean<T>,
 121                    boost::histogram::accumulators::mean<U>> {
 122   using type = boost::histogram::accumulators::mean<common_type_t<T, U>>;
 123 };
 124 } // namespace std
 125
 126 #endif
 127
 128 #endif