Division support for weighted_sum (#351)

include/boost/histogram/accumulators/weighted_sum.hpp

@@ -8,6 +8,7 @@
 #define BOOST_HISTOGRAM_ACCUMULATORS_WEIGHTED_SUM_HPP
 
 #include <boost/core/nvp.hpp>
+#include <boost/histogram/detail/square.hpp>
 #include <boost/histogram/fwd.hpp> // for weighted_sum<>
 #include <type_traits>
 
@@ -46,7 +47,7 @@ public:
   /// Increment by weight.
   weighted_sum& operator+=(const weight_type<value_type>& w) {
     sum_of_weights_ += w.value;
-    sum_of_weights_squared_ += w.value * w.value;
+    sum_of_weights_squared_ += detail::square(w.value);
     return *this;
   }
 
@@ -57,6 +58,20 @@ public:
     return *this;
   }
 
+  /// Divide by another weighted sum.
+  weighted_sum& operator/=(const weighted_sum& rhs) {
+    const auto v = sum_of_weights_;
+    const auto w = sum_of_weights_squared_;
+    const auto rv = rhs.sum_of_weights_;
+    const auto rw = rhs.sum_of_weights_squared_;
+    sum_of_weights_ /= rv;
+    // error propagation for independent a, b:
+    // c = a / b: var(c) = (var(a)/a^2 + var(b)/b^2) c^2
+    using detail::square;
+    sum_of_weights_squared_ = (w / square(v) + rw / square(rv)) * square(sum_of_weights_);
+    return *this;
+  }
+
   /// Scale by value.
   weighted_sum& operator*=(const_reference x) {
     sum_of_weights_ *= x;

test/accumulators_weighted_sum_test.cpp

@@ -8,6 +8,7 @@
 #include <boost/histogram/accumulators/ostream.hpp>
 #include <boost/histogram/accumulators/sum.hpp>
 #include <boost/histogram/accumulators/weighted_sum.hpp>
+#include <boost/histogram/detail/square.hpp>
 #include <boost/histogram/weight.hpp>
 #include <sstream>
 #include "throw_exception.hpp"
@@ -16,9 +17,10 @@
 using namespace boost::histogram;
 using namespace std::literals;
 
+using w_t = accumulators::weighted_sum<double>;
+
 int main() {
   {
-    using w_t = accumulators::weighted_sum<double>;
     w_t w;
     BOOST_TEST_EQ(str(w), "weighted_sum(0, 0)"s);
     BOOST_TEST_EQ(str(w, 20, false), "  weighted_sum(0, 0)"s);
@@ -78,5 +80,35 @@ int main() {
     BOOST_TEST_EQ(s_t() += s_t(), s_t());
   }
 
+  // division 1
+  {
+    w_t a{1, 2};
+    w_t b{2, 3};
+
+    auto c = a;
+    c /= b;
+
+    BOOST_TEST_EQ(c.value(), 0.5);
+    // error propagation for independent a and b:
+    // var(c)/c^2 = var(a)/a^2 + var(b)/b^2
+    BOOST_TEST_EQ(c.variance() / detail::square(c.value()),
+                  a.variance() / detail::square(a.value()) +
+                      b.variance() / detail::square(b.value()));
+  }
+
+  // division with implicit conversion
+  {
+    w_t a{1, 2};
+    double b = 2;
+
+    auto c = a;
+    c /= b;
+
+    BOOST_TEST_EQ(c.value(), 0.5);
+    // var(b) = b
+    BOOST_TEST_EQ(c.variance() / detail::square(c.value()),
+                  a.variance() / detail::square(a.value()) + 1.0 / b);
+  }
+
   return boost::report_errors();
 }

test/histogram_operators_test.cpp

@@ -160,8 +160,8 @@ void run_tests() {
     auto b = make_s(Tag(), std::vector<accumulators::weighted_sum<>>(), ia);
 
     a(0);
-    BOOST_TEST_EQ(a.at(0).variance(), 1);
     b(weight(3), 1);
+    BOOST_TEST_EQ(a.at(0).variance(), 1);
     BOOST_TEST_EQ(b.at(1).variance(), 9);
     auto c = a;
     c += b;
@@ -187,6 +187,50 @@ void run_tests() {
     BOOST_TEST_EQ(d.at(1).variance(), 11);
   }
 
+  // division operators with weighted storage
+  {
+    using w_t = accumulators::weighted_sum<>;
+    auto ia = axis::integer<int, axis::null_type, axis::option::none_t>(0, 2);
+    auto a = make_s(Tag(), std::vector<w_t>(), ia);
+    auto b = make_s(Tag(), std::vector<w_t>(), ia);
+
+    a(0);
+    a(1, weight(2));
+    b(weight(4), 0);
+    b(weight(3), 1);
+
+    w_t av[2] = {w_t{1, 1}, w_t{2, 4}};
+    w_t bv[2] = {w_t{4, 16}, w_t{3, 9}};
+    BOOST_TEST_EQ(a.at(0), av[0]);
+    BOOST_TEST_EQ(a.at(1), av[1]);
+    BOOST_TEST_EQ(b.at(0), bv[0]);
+    BOOST_TEST_EQ(b.at(1), bv[1]);
+
+    auto c = a;
+    c /= b;
+
+    w_t cv[2] = {av[0], av[1]};
+    cv[0] /= bv[0];
+    cv[1] /= bv[1];
+    BOOST_TEST_EQ(c.at(0), cv[0]);
+    BOOST_TEST_EQ(c.at(1), cv[1]);
+
+    // division by unweighted histogram
+    auto e = make_s(Tag(), std::vector<double>(), ia);
+
+    e(0);
+    e(0);
+    e(1);
+
+    auto f = a / e;
+    w_t fv[2] = {av[0], av[1]};
+    fv[0] /= e.at(0);
+    fv[1] /= e.at(1);
+
+    BOOST_TEST_EQ(f.at(0), fv[0]);
+    BOOST_TEST_EQ(f.at(1), fv[1]);
+  }
+
   // merging add
   {
     using C = axis::category<int, use_default, axis::option::growth_t>;