Accumulator fixes

* Bug-fix of weighted_sum: variance was not calculated correctly when normal histogram was added to weighted histogram with weighted_sum storage
* Several minor improvements to accumulator interfaces

.appveyor.yml

@@ -43,7 +43,9 @@ install:
   - sh: ./bootstrap.sh; ./b2 headers; cd libs/histogram
 
 test_script:
+  # on windows
   - cmd: ..\..\b2 %B2_OPTS% cxxstd=latest test//minimal test//serialization
+  # on linux
   - sh:
       B2="../../b2 ${B2_OPTS}";
       $B2 toolset=gcc-7 cxxstd=14 exception-handling=off rtti=off test//minimal &&

.travis.yml

@@ -25,7 +25,7 @@ matrix:
         - cd build
         - cmake ..
       script:
-        ctest -j2
+        ctest -j2 --output-on-failure
 
     - name: osx and superproject cmake and b2 minimal
       os: osx
@@ -50,7 +50,7 @@ matrix:
       script:
         mkdir build && cd build &&
         cmake -DBOOST_ENABLE_CMAKE=1 -DBoost_VERBOSE=1 .. &&
-        ctest --output-on-failure -R boost_histogram &&
+        ctest -j2 --output-on-failure -R boost_histogram &&
         cd .. &&
         ./bootstrap.sh &&
         ./b2 headers &&

include/boost/histogram/accumulators/mean.hpp

@@ -26,18 +26,21 @@ namespace accumulators {
 template <class RealType>
 class mean {
 public:
+  using value_type = RealType;
+  using const_reference = const value_type&;
+
   mean() = default;
-  mean(const RealType& n, const RealType& mean, const RealType& variance) noexcept
+  mean(const_reference n, const_reference mean, const_reference variance) noexcept
       : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {}
 
-  void operator()(const RealType& x) noexcept {
-    sum_ += static_cast<RealType>(1);
+  void operator()(const_reference x) noexcept {
+    sum_ += static_cast<value_type>(1);
     const auto delta = x - mean_;
     mean_ += delta / sum_;
     sum_of_deltas_squared_ += delta * (x - mean_);
   }
 
-  void operator()(const weight_type<RealType>& w, const RealType& x) noexcept {
+  void operator()(const weight_type<value_type>& w, const_reference x) noexcept {
     sum_ += w.value;
     const auto delta = x - mean_;
     mean_ += w.value * delta / sum_;
@@ -47,15 +50,15 @@ public:
   template <class T>
   mean& operator+=(const mean<T>& rhs) noexcept {
     if (sum_ != 0 || rhs.sum_ != 0) {
-      const auto tmp = mean_ * sum_ + static_cast<RealType>(rhs.mean_ * rhs.sum_);
+      const auto tmp = mean_ * sum_ + static_cast<value_type>(rhs.mean_ * rhs.sum_);
       sum_ += rhs.sum_;
       mean_ = tmp / sum_;
     }
-    sum_of_deltas_squared_ += static_cast<RealType>(rhs.sum_of_deltas_squared_);
+    sum_of_deltas_squared_ += static_cast<value_type>(rhs.sum_of_deltas_squared_);
     return *this;
   }
 
-  mean& operator*=(const RealType& s) noexcept {
+  mean& operator*=(const_reference s) noexcept {
     mean_ *= s;
     sum_of_deltas_squared_ *= s * s;
     return *this;
@@ -72,9 +75,9 @@ public:
     return !operator==(rhs);
   }
 
-  const RealType& count() const noexcept { return sum_; }
-  const RealType& value() const noexcept { return mean_; }
-  RealType variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
+  const_reference count() const noexcept { return sum_; }
+  const_reference value() const noexcept { return mean_; }
+  value_type variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
 
   template <class Archive>
   void serialize(Archive& ar, unsigned version) {
@@ -82,7 +85,7 @@ public:
       // read only
       std::size_t sum;
       ar& make_nvp("sum", sum);
-      sum_ = static_cast<RealType>(sum);
+      sum_ = static_cast<value_type>(sum);
     } else {
       ar& make_nvp("sum", sum_);
     }
@@ -91,7 +94,9 @@ public:
   }
 
 private:
-  RealType sum_ = 0, mean_ = 0, sum_of_deltas_squared_ = 0;
+  value_type sum_{};
+  value_type mean_{};
+  value_type sum_of_deltas_squared_{};
 };
 
 } // namespace accumulators

include/boost/histogram/accumulators/sum.hpp

@@ -9,7 +9,7 @@
 
 #include <boost/core/nvp.hpp>
 #include <boost/histogram/fwd.hpp> // for sum<>
-#include <cmath>
+#include <cmath>                   // std::abs
 #include <type_traits>
 
 namespace boost {
@@ -19,35 +19,44 @@ namespace accumulators {
 /**
   Uses Neumaier algorithm to compute accurate sums.
 
-  The algorithm uses memory for two floats and is three to
-  five times slower compared to a simple floating point
-  number used to accumulate a sum, but the relative error
-  of the sum is at the level of the machine precision,
-  independent of the number of samples.
+  The algorithm is an improved Kahan algorithm
+  (https://en.wikipedia.org/wiki/Kahan_summation_algorithm). The algorithm uses memory for
+  two numbers and is three to five times slower compared to using a single number to
+  accumulate a sum, but the relative error of the sum is at the level of the machine
+  precision, independent of the number of samples.
 
-  A. Neumaier, Zeitschrift fuer Angewandte Mathematik
-  und Mechanik 54 (1974) 39-51.
+  A. Neumaier, Zeitschrift fuer Angewandte Mathematik und Mechanik 54 (1974) 39-51.
 */
 template <class RealType>
 class sum {
 public:
+  using value_type = RealType;
+  using const_reference = const value_type&;
+
   sum() = default;
 
-  /// Initialize sum to value
-  explicit sum(const RealType& value) noexcept : large_(value) {}
+  /// Initialize sum to value and allow implicit conversion
+  sum(const_reference value) noexcept : sum(value, 0) {}
 
-  /// Set sum to value
-  sum& operator=(const RealType& value) noexcept {
-    large_ = value;
-    small_ = 0;
+  template <class T>
+  sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {}
+
+  /// Initialize sum explicitly with large and small parts
+  sum(const_reference large, const_reference small) noexcept
+      : large_(large), small_(small) {}
+
+  template <class T>
+  sum& operator=(const sum<T>& s) noexcept {
+    large_ = s.large_;
+    small_ = s.small_;
     return *this;
   }
 
   /// Increment sum by one
-  sum& operator++() { return operator+=(1); }
+  sum& operator++() noexcept { return operator+=(1); }
 
   /// Increment sum by value
-  sum& operator+=(const RealType& value) {
+  sum& operator+=(const_reference value) noexcept {
     // prevent compiler optimization from destroying the algorithm
     // when -ffast-math is enabled
     volatile RealType l;
@@ -66,31 +75,34 @@ public:
     return *this;
   }
 
+  /// Add another sum
+  sum& operator+=(const sum& s) noexcept {
+    operator+=(s.large_);
+    small_ += s.small_;
+    return *this;
+  }
+
   /// Scale by value
-  sum& operator*=(const RealType& value) {
+  sum& operator*=(const_reference value) noexcept {
     large_ *= value;
     small_ *= value;
     return *this;
   }
 
-  template <class T>
-  bool operator==(const sum<T>& rhs) const noexcept {
-    return large_ == rhs.large_ && small_ == rhs.small_;
+  bool operator==(const sum& rhs) const noexcept {
+    return large_ + small_ == rhs.large_ + rhs.small_;
   }
 
-  template <class T>
-  bool operator!=(const T& rhs) const noexcept {
-    return !operator==(rhs);
-  }
+  bool operator!=(const sum& rhs) const noexcept { return !operator==(rhs); }
 
   /// Return large part of the sum.
-  const RealType& large() const { return large_; }
+  const_reference large() const noexcept { return large_; }
 
   /// Return small part of the sum.
-  const RealType& small() const { return small_; }
+  const_reference small() const noexcept { return small_; }
 
-  // allow implicit conversion to RealType
-  operator RealType() const { return large_ + small_; }
+  // lossy conversion to value type must be explicit
+  explicit operator value_type() const noexcept { return large_ + small_; }
 
   template <class Archive>
   void serialize(Archive& ar, unsigned /* version */) {
@@ -98,9 +110,55 @@ public:
     ar& make_nvp("small", small_);
   }
 
+  // begin: extra operators to make sum behave like a regular number
+
+  sum& operator*=(const sum& rhs) noexcept {
+    const auto scale = static_cast<value_type>(rhs);
+    large_ *= scale;
+    small_ *= scale;
+    return *this;
+  }
+
+  sum operator*(const sum& rhs) const noexcept {
+    sum s = *this;
+    s *= rhs;
+    return s;
+  }
+
+  sum& operator/=(const sum& rhs) noexcept {
+    const auto scale = 1.0 / static_cast<value_type>(rhs);
+    large_ *= scale;
+    small_ *= scale;
+    return *this;
+  }
+
+  sum operator/(const sum& rhs) const noexcept {
+    sum s = *this;
+    s /= rhs;
+    return s;
+  }
+
+  bool operator<(const sum& rhs) const noexcept {
+    return operator value_type() < rhs.operator value_type();
+  }
+
+  bool operator>(const sum& rhs) const noexcept {
+    return operator value_type() > rhs.operator value_type();
+  }
+
+  bool operator<=(const sum& rhs) const noexcept {
+    return operator value_type() <= rhs.operator value_type();
+  }
+
+  bool operator>=(const sum& rhs) const noexcept {
+    return operator value_type() >= rhs.operator value_type();
+  }
+
+  // end: extra operators
+
 private:
-  RealType large_ = RealType();
-  RealType small_ = RealType();
+  value_type large_{};
+  value_type small_{};
 };
 
 } // namespace accumulators

include/boost/histogram/accumulators/weighted_mean.hpp

@@ -26,18 +26,21 @@ namespace accumulators {
 template <class RealType>
 class weighted_mean {
 public:
+  using value_type = RealType;
+  using const_reference = const value_type&;
+
   weighted_mean() = default;
-  weighted_mean(const RealType& wsum, const RealType& wsum2, const RealType& mean,
-                const RealType& variance)
+  weighted_mean(const_reference wsum, const_reference wsum2, const_reference mean,
+                const_reference variance)
       : sum_of_weights_(wsum)
       , sum_of_weights_squared_(wsum2)
       , weighted_mean_(mean)
       , sum_of_weighted_deltas_squared_(
             variance * (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_)) {}
 
-  void operator()(const RealType& x) { operator()(weight(1), x); }
+  void operator()(const_reference x) { operator()(weight(1), x); }
 
-  void operator()(const weight_type<RealType>& w, const RealType& x) {
+  void operator()(const weight_type<value_type>& w, const_reference x) {
     sum_of_weights_ += w.value;
     sum_of_weights_squared_ += w.value * w.value;
     const auto delta = x - weighted_mean_;
@@ -49,9 +52,9 @@ public:
   weighted_mean& operator+=(const weighted_mean<T>& rhs) {
     if (sum_of_weights_ != 0 || rhs.sum_of_weights_ != 0) {
       const auto tmp = weighted_mean_ * sum_of_weights_ +
-                       static_cast<RealType>(rhs.weighted_mean_ * rhs.sum_of_weights_);
-      sum_of_weights_ += static_cast<RealType>(rhs.sum_of_weights_);
-      sum_of_weights_squared_ += static_cast<RealType>(rhs.sum_of_weights_squared_);
+                       static_cast<value_type>(rhs.weighted_mean_ * rhs.sum_of_weights_);
+      sum_of_weights_ += static_cast<value_type>(rhs.sum_of_weights_);
+      sum_of_weights_squared_ += static_cast<value_type>(rhs.sum_of_weights_squared_);
       weighted_mean_ = tmp / sum_of_weights_;
     }
     sum_of_weighted_deltas_squared_ +=
@@ -59,7 +62,7 @@ public:
     return *this;
   }
 
-  weighted_mean& operator*=(const RealType& s) {
+  weighted_mean& operator*=(const_reference s) {
     weighted_mean_ *= s;
     sum_of_weighted_deltas_squared_ *= s * s;
     return *this;
@@ -78,12 +81,12 @@ public:
     return !operator==(rhs);
   }
 
-  const RealType& sum_of_weights() const noexcept { return sum_of_weights_; }
-  const RealType& sum_of_weights_squared() const noexcept {
+  const_reference sum_of_weights() const noexcept { return sum_of_weights_; }
+  const_reference sum_of_weights_squared() const noexcept {
     return sum_of_weights_squared_;
   }
-  const RealType& value() const noexcept { return weighted_mean_; }
-  RealType variance() const {
+  const_reference value() const noexcept { return weighted_mean_; }
+  value_type variance() const {
     return sum_of_weighted_deltas_squared_ /
            (sum_of_weights_ - sum_of_weights_squared_ / sum_of_weights_);
   }
@@ -97,8 +100,10 @@ public:
   }
 
 private:
-  RealType sum_of_weights_ = RealType(), sum_of_weights_squared_ = RealType(),
-           weighted_mean_ = RealType(), sum_of_weighted_deltas_squared_ = RealType();
+  value_type sum_of_weights_{};
+  value_type sum_of_weights_squared_{};
+  value_type weighted_mean_{};
+  value_type sum_of_weighted_deltas_squared_{};
 };
 
 } // namespace accumulators

include/boost/histogram/accumulators/weighted_sum.hpp

@@ -20,54 +20,63 @@ template <class RealType>
 class weighted_sum {
 public:
   using value_type = RealType;
-  using const_reference = const RealType&;
+  using const_reference = const value_type&;
 
   weighted_sum() = default;
-  explicit weighted_sum(const RealType& value) noexcept
-      : sum_of_weights_(value), sum_of_weights_squared_(value) {}
-  weighted_sum(const RealType& value, const RealType& variance) noexcept
+
+  /// Initialize sum to value and allow implicit conversion
+  weighted_sum(const_reference value) noexcept : weighted_sum(value, value) {}
+
+  template <class T>
+  weighted_sum(const weighted_sum<T>& s) noexcept
+      : weighted_sum(s.value(), s.variance()) {}
+
+  /// Initialize sum to value and variance
+  weighted_sum(const_reference value, const_reference variance) noexcept
       : sum_of_weights_(value), sum_of_weights_squared_(variance) {}
 
-  /// Increment by one.
-  weighted_sum& operator++() { return operator+=(1); }
-
-  /// Increment by value.
   template <class T>
-  weighted_sum& operator+=(const T& value) {
-    sum_of_weights_ += value;
-    sum_of_weights_squared_ += value * value;
+  weighted_sum& operator=(const weighted_sum<T>& s) noexcept {
+    sum_of_weights_ = s.sum_of_weights_;
+    sum_of_weights_squared_ = s.sum_of_weights_squared_;
+    return *this;
+  }
+
+  /// Increment by one.
+  weighted_sum& operator++() {
+    ++sum_of_weights_;
+    ++sum_of_weights_squared_;
+    return *this;
+  }
+
+  /// Increment by weight.
+  template <class T>
+  weighted_sum& operator+=(const weight_type<T>& w) {
+    sum_of_weights_ += w.value;
+    sum_of_weights_squared_ += w.value * w.value;
     return *this;
   }
 
   /// Added another weighted sum.
-  template <class T>
-  weighted_sum& operator+=(const weighted_sum<T>& rhs) {
-    sum_of_weights_ += static_cast<RealType>(rhs.sum_of_weights_);
-    sum_of_weights_squared_ += static_cast<RealType>(rhs.sum_of_weights_squared_);
+  weighted_sum& operator+=(const weighted_sum& rhs) {
+    sum_of_weights_ += rhs.sum_of_weights_;
+    sum_of_weights_squared_ += rhs.sum_of_weights_squared_;
     return *this;
   }
 
   /// Scale by value.
-  weighted_sum& operator*=(const RealType& x) {
+  weighted_sum& operator*=(const_reference x) {
     sum_of_weights_ *= x;
     sum_of_weights_squared_ *= x * x;
     return *this;
   }
 
-  bool operator==(const RealType& rhs) const noexcept {
-    return sum_of_weights_ == rhs && sum_of_weights_squared_ == rhs;
-  }
-
-  template <class T>
-  bool operator==(const weighted_sum<T>& rhs) const noexcept {
+  bool operator==(const weighted_sum& rhs) const noexcept {
     return sum_of_weights_ == rhs.sum_of_weights_ &&
            sum_of_weights_squared_ == rhs.sum_of_weights_squared_;
   }
 
-  template <class T>
-  bool operator!=(const T& rhs) const noexcept {
-    return !operator==(rhs);
-  }
+  bool operator!=(const weighted_sum& rhs) const noexcept { return !operator==(rhs); }
 
   /// Return value of the sum.
   const_reference value() const noexcept { return sum_of_weights_; }
@@ -85,8 +94,8 @@ public:
   }
 
 private:
-  value_type sum_of_weights_ = value_type();
-  value_type sum_of_weights_squared_ = value_type();
+  value_type sum_of_weights_{};
+  value_type sum_of_weights_squared_{};
 };
 
 } // namespace accumulators

include/boost/histogram/detail/fill.hpp

@@ -18,6 +18,7 @@
 #include <boost/histogram/detail/linearize.hpp>
 #include <boost/histogram/detail/make_default.hpp>
 #include <boost/histogram/detail/optional_index.hpp>
+#include <boost/histogram/detail/priority.hpp>
 #include <boost/histogram/detail/tuple_slice.hpp>
 #include <boost/histogram/fwd.hpp>
 #include <boost/mp11/algorithm.hpp>
@@ -134,50 +135,62 @@ struct storage_grower {
 };
 
 template <class T, class... Us>
-void fill_storage_4(mp11::mp_false, T&& t, Us&&... args) noexcept {
-  t(std::forward<Us>(args)...);
-}
-
-template <class T>
-void fill_storage_4(mp11::mp_true, T&& t) noexcept {
-  ++t;
+auto fill_storage_element_impl(priority<2>, T&& t, const Us&... args) noexcept
+    -> decltype(t(args...), void()) {
+  t(args...);
 }
 
 template <class T, class U>
-void fill_storage_4(mp11::mp_true, T&& t, U&& w) noexcept {
+auto fill_storage_element_impl(priority<1>, T&& t, const weight_type<U>& w) noexcept
+    -> decltype(t += w, void()) {
+  t += w;
+}
+
+// fallback for arithmetic types and accumulators that do not handle the weight
+template <class T, class U>
+auto fill_storage_element_impl(priority<0>, T&& t, const weight_type<U>& w) noexcept
+    -> decltype(t += w.value, void()) {
   t += w.value;
 }
 
-template <class T, class... Us>
-void fill_storage_3(T&& t, Us&&... args) noexcept {
-  fill_storage_4(has_operator_preincrement<std::decay_t<T>>{}, std::forward<T>(t),
-                 std::forward<Us>(args)...);
+template <class T>
+auto fill_storage_element_impl(priority<1>, T&& t) noexcept -> decltype(++t, void()) {
+  ++t;
 }
 
+template <class T, class... Us>
+auto fill_storage_element(T&& t, const Us&... args) noexcept {
+  fill_storage_element_impl(priority<2>{}, std::forward<T>(t), args...);
+}
+
+// t may be a proxy and then it is an rvalue reference, not an lvalue reference
 template <class IW, class IS, class T, class U>
 void fill_storage_2(IW, IS, T&& t, U&& u) noexcept {
   mp11::tuple_apply(
-      [&](auto&&... args) {
-        fill_storage_3(std::forward<T>(t), std::get<IW::value>(u), args...);
+      [&](const auto&... args) {
+        fill_storage_element(std::forward<T>(t), std::get<IW::value>(u), args...);
       },
       std::get<IS::value>(u).value);
 }
 
+// t may be a proxy and then it is an rvalue reference, not an lvalue reference
 template <class IS, class T, class U>
 void fill_storage_2(mp11::mp_int<-1>, IS, T&& t, const U& u) noexcept {
   mp11::tuple_apply(
-      [&](const auto&... args) { fill_storage_3(std::forward<T>(t), args...); },
+      [&](const auto&... args) { fill_storage_element(std::forward<T>(t), args...); },
       std::get<IS::value>(u).value);
 }
 
+// t may be a proxy and then it is an rvalue reference, not an lvalue reference
 template <class IW, class T, class U>
 void fill_storage_2(IW, mp11::mp_int<-1>, T&& t, const U& u) noexcept {
-  fill_storage_3(std::forward<T>(t), std::get<IW::value>(u));
+  fill_storage_element(std::forward<T>(t), std::get<IW::value>(u));
 }
 
+// t may be a proxy and then it is an rvalue reference, not an lvalue reference
 template <class T, class U>
 void fill_storage_2(mp11::mp_int<-1>, mp11::mp_int<-1>, T&& t, const U&) noexcept {
-  fill_storage_3(std::forward<T>(t));
+  fill_storage_element(std::forward<T>(t));
 }
 
 template <class IW, class IS, class Storage, class Index, class Args>

include/boost/histogram/detail/fill_n.hpp

@@ -159,7 +159,7 @@ template <class S, class Index, class... Ts>
 void fill_n_storage(S& s, const Index idx, Ts&&... p) noexcept {
   if (is_valid(idx)) {
     BOOST_ASSERT(idx < s.size());
-    fill_storage_3(s[idx], *p.first...);
+    fill_storage_element(s[idx], *p.first...);
   }
   fold((p.second ? ++p.first : 0)...);
 }
@@ -168,8 +168,7 @@ template <class S, class Index, class T, class... Ts>
 void fill_n_storage(S& s, const Index idx, weight_type<T>&& w, Ts&&... ps) noexcept {
   if (is_valid(idx)) {
     BOOST_ASSERT(idx < s.size());
-    fill_storage_3(s[idx], weight_type<decltype(*w.value.first)>{*w.value.first},
-                   *ps.first...);
+    fill_storage_element(s[idx], weight(*w.value.first), *ps.first...);
   }
   if (w.value.second) ++w.value.first;
   fold((ps.second ? ++ps.first : 0)...);

include/boost/histogram/detail/priority.hpp

@@ -0,0 +1,28 @@
+// Copyright 2019 Hans Dembinski
+//
+// Distributed under the Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_HISTOGRAM_DETAIL_PRIORITY_HPP
+#define BOOST_HISTOGRAM_DETAIL_PRIORITY_HPP
+
+#include <cstdint>
+
+namespace boost {
+namespace histogram {
+namespace detail {
+
+// priority is used to priorise ambiguous overloads
+
+template <std::size_t N>
+struct priority : priority<(N - 1)> {};
+
+template <>
+struct priority<0> {};
+
+} // namespace detail
+} // namespace histogram
+} // namespace boost
+
+#endif

include/boost/histogram/storage_adaptor.hpp

@@ -247,8 +247,8 @@ struct map_impl : T {
     }
 
     template <class... Ts>
-    decltype(auto) operator()(Ts&&... args) {
-      return map->operator[](idx)(std::forward<Ts>(args)...);
+    auto operator()(const Ts&... args) -> decltype(std::declval<value_type>()(args...)) {
+      return (*map)[idx](args...);
     }
 
     map_impl* map;

test/accumulators_serialization_test.cpp

@@ -77,8 +77,8 @@ int main(int argc, char** argv) {
     const auto filename =
         join(argv[1], "accumulators_serialization_test_weighted_sum.xml");
     accumulators::weighted_sum<> a;
-    a += 1;
-    a += 10;
+    a += weight(1);
+    a += weight(10);
     print_xml(filename, a);
 
     accumulators::weighted_sum<> b;

test/accumulators_test.cpp

@@ -48,7 +48,7 @@ int main() {
     BOOST_TEST_EQ(w, 1);
     BOOST_TEST_NE(w, 2);
 
-    w += 2;
+    w += weight(2);
     BOOST_TEST_EQ(w.value(), 3);
     BOOST_TEST_EQ(w.variance(), 5);
     BOOST_TEST_EQ(w, w_t(3, 5));
@@ -66,7 +66,7 @@ int main() {
     BOOST_TEST_EQ(u, w_t(2, 4));
 
     w_t v(0);
-    v += 2;
+    v += weight(2);
     BOOST_TEST_EQ(u, v);
 
     // conversion to RealType
@@ -185,6 +185,19 @@ int main() {
     BOOST_TEST_EQ(sum.large(), 0);
     BOOST_TEST_EQ(sum.small(), 2);
 
+    sum = s_t{1e100, 1};
+    sum += s_t{1e100, 1};
+    BOOST_TEST_EQ(sum, (s_t{2e100, 2}));
+    sum = s_t{1e100, 1};
+    sum += s_t{1, 0};
+    BOOST_TEST_EQ(sum, (s_t{1e100, 2}));
+    sum = s_t{1, 0};
+    sum += s_t{1e100, 1};
+    BOOST_TEST_EQ(sum, (s_t{1e100, 2}));
+    sum = s_t{0, 1};
+    sum += s_t{1, 0};
+    BOOST_TEST_EQ(sum, (s_t{1, 1}));
+
     accumulators::sum<double> a(3), b(2), c(3);
     BOOST_TEST_LT(b, c);
     BOOST_TEST_LE(b, c);
@@ -201,9 +214,9 @@ int main() {
     s_t w;
 
     ++w;
-    w += 1e100;
+    w += weight(1e100);
     ++w;
-    w += -1e100;
+    w += weight(-1e100);
 
     BOOST_TEST_EQ(w.value(), 2);
     BOOST_TEST_EQ(w.variance(), 2e200);

test/histogram_operators_test.cpp

@@ -153,6 +153,16 @@ void run_tests() {
     BOOST_TEST_EQ(d.at(0).variance(), 1);
     BOOST_TEST_EQ(d.at(1).value(), 3);
     BOOST_TEST_EQ(d.at(1).variance(), 9);
+
+    // add unweighted histogram
+    auto e = make_s(Tag(), std::vector<int>(), ia);
+    std::fill(e.begin(), e.end(), 2);
+
+    d += e;
+    BOOST_TEST_EQ(d.at(0).value(), 3);
+    BOOST_TEST_EQ(d.at(0).variance(), 3);
+    BOOST_TEST_EQ(d.at(1).value(), 5);
+    BOOST_TEST_EQ(d.at(1).variance(), 11);
   }
 
   // bad operations

test/histogram_test.cpp

@@ -410,7 +410,7 @@ void run_tests() {
     BOOST_TEST_EQ(a[4], 3);
   }
 
-  // histogram_reset
+  // histogram reset
   {
     auto h = make(Tag(), axis::integer<int, axis::null_type, axis::option::none_t>(0, 2));
     h(0);
@@ -424,10 +424,9 @@ void run_tests() {
     BOOST_TEST_EQ(algorithm::sum(h), 0);
   }
 
-  // using containers or input and output
+  // using containers for input and output
   {
-    auto h = make_s(Tag(), weight_storage(), axis::integer<>(0, 2),
-                    axis::integer<double>(2, 4));
+    auto h = make(Tag(), axis::integer<>(0, 2), axis::integer<double>(2, 4));
     // tuple in
     h(std::make_tuple(0, 2.0));
     h(std::make_tuple(1, 3.0));
@@ -436,9 +435,9 @@ void run_tests() {
     auto i11 = std::make_tuple(1, 1);
 
     // tuple out
-    BOOST_TEST_EQ(h.at(i00).value(), 1);
-    BOOST_TEST_EQ(h[i00].value(), 1);
-    BOOST_TEST_EQ(h[i11].value(), 1);
+    BOOST_TEST_EQ(h.at(i00), 1);
+    BOOST_TEST_EQ(h[i00], 1);
+    BOOST_TEST_EQ(h[i11], 1);
 
     // iterable out
     int j11[] = {1, 1};
@@ -455,10 +454,8 @@ void run_tests() {
     h(std::make_tuple(weight(2), 0, 2.0));
     h(std::make_tuple(1, 3.0, weight(2)));
 
-    BOOST_TEST_EQ(h.at(i00).value(), 3);
-    BOOST_TEST_EQ(h[i00].value(), 3);
-    BOOST_TEST_EQ(h.at(i11).variance(), 5);
-    BOOST_TEST_EQ(h[i11].variance(), 5);
+    BOOST_TEST_EQ(h.at(i00), 3);
+    BOOST_TEST_EQ(h[i00], 3);
 
     // test special case of 1-dimensional histogram, which should unpack
     // 1-dimensional tuple normally, but forward larger tuples to the axis

test/storage_adaptor_test.cpp

@@ -218,7 +218,7 @@ int main() {
     ++a[0];
     a[0] += 1;
     a[0] += 2;
-    a[0] += accumulators::weighted_sum<double>(1, 0);
+    a[0] += accumulators::weighted_sum<double>(1, 2);
     BOOST_TEST_EQ(a[0].value(), 5);
     BOOST_TEST_EQ(a[0].variance(), 6);
     a[0] *= 2;