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/core/nvp.hpp> 11 #include <boost/histogram/fwd.hpp> // for mean<> 12 #include <boost/throw_exception.hpp> 13 #include <cassert> 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 ValueType> 27 class mean { 28 public: 29 using value_type = ValueType; 30 using const_reference = const value_type&; 31 32 mean() = default; 33 34 /// Allow implicit conversion from mean<T> 35 template <class T> mean(const mean<T> & o)36 mean(const mean<T>& o) noexcept 37 : sum_{o.sum_}, mean_{o.mean_}, sum_of_deltas_squared_{o.sum_of_deltas_squared_} {} 38 39 /// Initialize to external count, mean, and variance mean(const_reference n,const_reference mean,const_reference variance)40 mean(const_reference n, const_reference mean, const_reference variance) noexcept 41 : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {} 42 43 /// Insert sample x operator ()(const_reference x)44 void operator()(const_reference x) noexcept { 45 sum_ += static_cast<value_type>(1); 46 const auto delta = x - mean_; 47 mean_ += delta / sum_; 48 sum_of_deltas_squared_ += delta * (x - mean_); 49 } 50 51 /// Insert sample x with weight w operator ()(const weight_type<value_type> & w,const_reference x)52 void operator()(const weight_type<value_type>& w, const_reference x) noexcept { 53 sum_ += w.value; 54 const auto delta = x - mean_; 55 mean_ += w.value * delta / sum_; 56 sum_of_deltas_squared_ += w.value * delta * (x - mean_); 57 } 58 59 /// Add another mean accumulator operator +=(const mean & rhs)60 mean& operator+=(const mean& rhs) noexcept { 61 if (sum_ != 0 || rhs.sum_ != 0) { 62 const auto tmp = mean_ * sum_ + rhs.mean_ * rhs.sum_; 63 sum_ += rhs.sum_; 64 mean_ = tmp / sum_; 65 } 66 sum_of_deltas_squared_ += rhs.sum_of_deltas_squared_; 67 return *this; 68 } 69 70 /** Scale by value 71 72 This acts as if all samples were scaled by the value. 73 */ operator *=(const_reference s)74 mean& operator*=(const_reference s) noexcept { 75 mean_ *= s; 76 sum_of_deltas_squared_ *= s * s; 77 return *this; 78 } 79 operator ==(const mean & rhs) const80 bool operator==(const mean& rhs) const noexcept { 81 return sum_ == rhs.sum_ && mean_ == rhs.mean_ && 82 sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_; 83 } 84 operator !=(const mean & rhs) const85 bool operator!=(const mean& rhs) const noexcept { return !operator==(rhs); } 86 87 /// Return how many samples were accumulated count() const88 const_reference count() const noexcept { return sum_; } 89 90 /// Return mean value of accumulated samples value() const91 const_reference value() const noexcept { return mean_; } 92 93 /// Return variance of accumulated samples variance() const94 value_type variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); } 95 96 template <class Archive> serialize(Archive & ar,unsigned version)97 void serialize(Archive& ar, unsigned version) { 98 if (version == 0) { 99 // read only 100 std::size_t sum; 101 ar& make_nvp("sum", sum); 102 sum_ = static_cast<value_type>(sum); 103 } else { 104 ar& make_nvp("sum", sum_); 105 } 106 ar& make_nvp("mean", mean_); 107 ar& make_nvp("sum_of_deltas_squared", sum_of_deltas_squared_); 108 } 109 110 private: 111 value_type sum_{}; 112 value_type mean_{}; 113 value_type sum_of_deltas_squared_{}; 114 }; 115 116 } // namespace accumulators 117 } // namespace histogram 118 } // namespace boost 119 120 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED 121 122 namespace boost { 123 namespace serialization { 124 125 template <class T> 126 struct version; 127 128 // version 1 for boost::histogram::accumulators::mean<T> 129 template <class T> 130 struct version<boost::histogram::accumulators::mean<T>> : std::integral_constant<int, 1> { 131 }; 132 133 } // namespace serialization 134 } // namespace boost 135 136 namespace std { 137 template <class T, class U> 138 /// Specialization for boost::histogram::accumulators::mean. 139 struct common_type<boost::histogram::accumulators::mean<T>, 140 boost::histogram::accumulators::mean<U>> { 141 using type = boost::histogram::accumulators::mean<common_type_t<T, U>>; 142 }; 143 } // namespace std 144 145 #endif 146 147 #endif 148