histogram/accumulators/weighted_sum.hpp

0001 // Copyright 2015-2018 Hans Dembinski
0002 //
0003 // Distributed under the Boost Software License, version 1.0.
0004 // (See accompanying file LICENSE_1_0.txt
0005 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_SUM_HPP
0008 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_SUM_HPP
0009
0010 #include <boost/core/nvp.hpp>
0011 #include <boost/histogram/detail/square.hpp>
0012 #include <boost/histogram/fwd.hpp> // for weighted_sum<>
0013 #include <type_traits>
0014
0015 namespace boost {
0016 namespace histogram {
0017 namespace accumulators {
0018
0019 /// Holds sum of weights and its variance estimate
0020 template <class ValueType>
0021 class weighted_sum {
0022 public:
0023   using value_type = ValueType;
0024   using const_reference = const value_type&;
0025
0026   weighted_sum() = default;
0027
0028   /// Initialize sum to value and allow implicit conversion
0029   weighted_sum(const_reference value) noexcept : weighted_sum(value, value) {}
0030
0031   /// Allow implicit conversion from sum<T>
0032   template <class T>
0033   weighted_sum(const weighted_sum<T>& s) noexcept
0034       : weighted_sum(s.value(), s.variance()) {}
0035
0036   /// Initialize sum to value and variance
0037   weighted_sum(const_reference value, const_reference variance) noexcept
0038       : sum_of_weights_(value), sum_of_weights_squared_(variance) {}
0039
0040   /// Increment by one.
0041   weighted_sum& operator++() {
0042     ++sum_of_weights_;
0043     ++sum_of_weights_squared_;
0044     return *this;
0045   }
0046
0047   /// Increment by weight.
0048   weighted_sum& operator+=(const weight_type<value_type>& w) {
0049     sum_of_weights_ += w.value;
0050     sum_of_weights_squared_ += detail::square(w.value);
0051     return *this;
0052   }
0053
0054   /// Added another weighted sum.
0055   weighted_sum& operator+=(const weighted_sum& rhs) {
0056     sum_of_weights_ += rhs.sum_of_weights_;
0057     sum_of_weights_squared_ += rhs.sum_of_weights_squared_;
0058     return *this;
0059   }
0060
0061   /// Divide by another weighted sum.
0062   weighted_sum& operator/=(const weighted_sum& rhs) {
0063     const auto v = sum_of_weights_;
0064     const auto w = sum_of_weights_squared_;
0065     const auto rv = rhs.sum_of_weights_;
0066     const auto rw = rhs.sum_of_weights_squared_;
0067     sum_of_weights_ /= rv;
0068     // error propagation for independent a, b:
0069     // c = a / b: var(c) = (var(a)/a^2 + var(b)/b^2) c^2
0070     using detail::square;
0071     sum_of_weights_squared_ = (w / square(v) + rw / square(rv)) * square(sum_of_weights_);
0072     return *this;
0073   }
0074
0075   /// Scale by value.
0076   weighted_sum& operator*=(const_reference x) {
0077     sum_of_weights_ *= x;
0078     sum_of_weights_squared_ *= x * x;
0079     return *this;
0080   }
0081
0082   bool operator==(const weighted_sum& rhs) const noexcept {
0083     return sum_of_weights_ == rhs.sum_of_weights_ &&
0084            sum_of_weights_squared_ == rhs.sum_of_weights_squared_;
0085   }
0086
0087   bool operator!=(const weighted_sum& rhs) const noexcept { return !operator==(rhs); }
0088
0089   /// Return value of the sum.
0090   const_reference value() const noexcept { return sum_of_weights_; }
0091
0092   /// Return estimated variance of the sum.
0093   const_reference variance() const noexcept { return sum_of_weights_squared_; }
0094
0095   // lossy conversion must be explicit
0096   explicit operator const_reference() const { return sum_of_weights_; }
0097
0098   template <class Archive>
0099   void serialize(Archive& ar, unsigned /* version */) {
0100     ar& make_nvp("sum_of_weights", sum_of_weights_);
0101     ar& make_nvp("sum_of_weights_squared", sum_of_weights_squared_);
0102   }
0103
0104 private:
0105   value_type sum_of_weights_{};
0106   value_type sum_of_weights_squared_{};
0107 };
0108
0109 } // namespace accumulators
0110 } // namespace histogram
0111 } // namespace boost
0112
0113 #ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
0114 namespace std {
0115 template <class T, class U>
0116 struct common_type<boost::histogram::accumulators::weighted_sum<T>,
0117                    boost::histogram::accumulators::weighted_sum<U>> {
0118   using type = boost::histogram::accumulators::weighted_sum<common_type_t<T, U>>;
0119 };
0120
0121 template <class T, class U>
0122 struct common_type<boost::histogram::accumulators::weighted_sum<T>, U> {
0123   using type = boost::histogram::accumulators::weighted_sum<common_type_t<T, U>>;
0124 };
0125
0126 template <class T, class U>
0127 struct common_type<T, boost::histogram::accumulators::weighted_sum<U>> {
0128   using type = boost::histogram::accumulators::weighted_sum<common_type_t<T, U>>;
0129 };
0130 } // namespace std
0131 #endif
0132
0133 #endif