mean.hpp 3.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128
  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. #ifndef BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
  7. #define BOOST_HISTOGRAM_ACCUMULATORS_MEAN_HPP
  8. #include <boost/assert.hpp>
  9. #include <boost/core/nvp.hpp>
  10. #include <boost/histogram/fwd.hpp> // for mean<>
  11. #include <boost/throw_exception.hpp>
  12. #include <stdexcept>
  13. #include <type_traits>
  14. namespace boost {
  15. namespace histogram {
  16. namespace accumulators {
  17. /** Calculates mean and variance of sample.
  18. Uses Welfords's incremental algorithm to improve the numerical
  19. stability of mean and variance computation.
  20. */
  21. template <class RealType>
  22. class mean {
  23. public:
  24. mean() = default;
  25. mean(const RealType& n, const RealType& mean, const RealType& variance) noexcept
  26. : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {}
  27. void operator()(const RealType& x) noexcept {
  28. sum_ += static_cast<RealType>(1);
  29. const auto delta = x - mean_;
  30. mean_ += delta / sum_;
  31. sum_of_deltas_squared_ += delta * (x - mean_);
  32. }
  33. void operator()(const weight_type<RealType>& w, const RealType& x) noexcept {
  34. sum_ += w.value;
  35. const auto delta = x - mean_;
  36. mean_ += w.value * delta / sum_;
  37. sum_of_deltas_squared_ += w.value * delta * (x - mean_);
  38. }
  39. template <class T>
  40. mean& operator+=(const mean<T>& rhs) noexcept {
  41. if (sum_ != 0 || rhs.sum_ != 0) {
  42. const auto tmp = mean_ * sum_ + static_cast<RealType>(rhs.mean_ * rhs.sum_);
  43. sum_ += rhs.sum_;
  44. mean_ = tmp / sum_;
  45. }
  46. sum_of_deltas_squared_ += static_cast<RealType>(rhs.sum_of_deltas_squared_);
  47. return *this;
  48. }
  49. mean& operator*=(const RealType& s) noexcept {
  50. mean_ *= s;
  51. sum_of_deltas_squared_ *= s * s;
  52. return *this;
  53. }
  54. template <class T>
  55. bool operator==(const mean<T>& rhs) const noexcept {
  56. return sum_ == rhs.sum_ && mean_ == rhs.mean_ &&
  57. sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_;
  58. }
  59. template <class T>
  60. bool operator!=(const mean<T>& rhs) const noexcept {
  61. return !operator==(rhs);
  62. }
  63. const RealType& count() const noexcept { return sum_; }
  64. const RealType& value() const noexcept { return mean_; }
  65. RealType variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
  66. template <class Archive>
  67. void serialize(Archive& ar, unsigned version) {
  68. if (version == 0) {
  69. // read only
  70. std::size_t sum;
  71. ar& make_nvp("sum", sum);
  72. sum_ = static_cast<RealType>(sum);
  73. } else {
  74. ar& make_nvp("sum", sum_);
  75. }
  76. ar& make_nvp("mean", mean_);
  77. ar& make_nvp("sum_of_deltas_squared", sum_of_deltas_squared_);
  78. }
  79. private:
  80. RealType sum_ = 0, mean_ = 0, sum_of_deltas_squared_ = 0;
  81. };
  82. } // namespace accumulators
  83. } // namespace histogram
  84. } // namespace boost
  85. #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
  86. namespace boost {
  87. namespace serialization {
  88. template <class T>
  89. struct version;
  90. // version 1 for boost::histogram::accumulators::mean<RealType>
  91. template <class RealType>
  92. struct version<boost::histogram::accumulators::mean<RealType>>
  93. : std::integral_constant<int, 1> {};
  94. } // namespace serialization
  95. } // namespace boost
  96. namespace std {
  97. template <class T, class U>
  98. /// Specialization for boost::histogram::accumulators::mean.
  99. struct common_type<boost::histogram::accumulators::mean<T>,
  100. boost::histogram::accumulators::mean<U>> {
  101. using type = boost::histogram::accumulators::mean<common_type_t<T, U>>;
  102. };
  103. } // namespace std
  104. #endif
  105. #endif