feat: sparse storage add/sub

include/boost/histogram/detail/detect.hpp

@@ -85,6 +85,10 @@ BOOST_HISTOGRAM_DETAIL_DETECT(is_map_like, ((typename T::key_type*)nullptr,
                                             (typename T::mapped_type*)nullptr,
                                             std::begin(t), std::end(t)));
 
+// storage exposes the container of non-empty cells (sparse storage)
+BOOST_HISTOGRAM_DETAIL_DETECT(has_node_access,
+                              (std::begin(t.node_access()), std::end(t.node_access())));
+
 // ok: is_axis is false for axis::variant, because T::index is templated
 BOOST_HISTOGRAM_DETAIL_DETECT(is_axis, (t.size(), &T::index));
 

include/boost/histogram/histogram.hpp

@@ -450,8 +450,7 @@ class histogram : detail::mutex_base<Axes, Storage> {
   operator+=(const histogram<A, S>& rhs) {
     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
-    auto rit = unsafe_access::storage(rhs).begin();
-    for (auto&& x : storage_) x += *rit++;
+    add_storage(unsafe_access::storage(rhs));
     return *this;
   }
 
@@ -467,8 +466,7 @@ class histogram : detail::mutex_base<Axes, Storage> {
   operator+=(const histogram<axes_type, S>& rhs) {
     const auto& raxes = unsafe_access::axes(rhs);
     if (detail::axes_equal(axes_, unsafe_access::axes(rhs))) {
-      auto rit = unsafe_access::storage(rhs).begin();
-      for (auto&& x : storage_) x += *rit++;
+      add_storage(unsafe_access::storage(rhs));
       return *this;
     }
 
@@ -501,8 +499,7 @@ class histogram : detail::mutex_base<Axes, Storage> {
   operator-=(const histogram<A, S>& rhs) {
     if (!detail::axes_equal(axes_, unsafe_access::axes(rhs)))
       BOOST_THROW_EXCEPTION(std::invalid_argument("axes of histograms differ"));
-    auto rit = unsafe_access::storage(rhs).begin();
-    for (auto&& x : storage_) x -= *rit++;
+    sub_storage(unsafe_access::storage(rhs));
     return *this;
   }
 
@@ -611,6 +608,36 @@ class histogram : detail::mutex_base<Axes, Storage> {
   }
 
 private:
+  template <class OtherStorage>
+  void add_storage(const OtherStorage& rs) {
+    detail::static_if<detail::has_node_access<OtherStorage>>(
+        [this](const auto& s) {
+          // sparse rhs: visit only the stored cells, absent cells add zero;
+          // the bounds check guards against stray keys from a corrupt archive
+          for (auto&& kv : s.node_access())
+            if (kv.first < storage_.size()) storage_[kv.first] += kv.second;
+        },
+        [this](const auto& s) {
+          auto rit = s.begin();
+          for (auto&& x : storage_) x += *rit++;
+        },
+        rs);
+  }
+
+  template <class OtherStorage>
+  void sub_storage(const OtherStorage& rs) {
+    detail::static_if<detail::has_node_access<OtherStorage>>(
+        [this](const auto& s) {
+          for (auto&& kv : s.node_access())
+            if (kv.first < storage_.size()) storage_[kv.first] -= kv.second;
+        },
+        [this](const auto& s) {
+          auto rit = s.begin();
+          for (auto&& x : storage_) x -= *rit++;
+        },
+        rs);
+  }
+
   axes_type axes_;
   storage_type storage_;
   std::size_t offset_ = 0;

include/boost/histogram/storage_adaptor.hpp

@@ -332,6 +332,10 @@ struct map_impl : T {
 
   std::size_t size() const noexcept { return size_; }
 
+  // access to the underlying map with only the non-empty cells,
+  // allows histogram operations to skip empty cells
+  const T& node_access() const noexcept { return *this; }
+
   template <class Archive>
   void serialize(Archive& ar, unsigned /* version */) {
     ar& make_nvp("size", size_);

test/detail_detect_test.cpp

@@ -12,6 +12,7 @@
 #include <boost/histogram/axis/variable.hpp>
 #include <boost/histogram/axis/variant.hpp>
 #include <boost/histogram/detail/detect.hpp>
+#include <boost/histogram/storage_adaptor.hpp>
 #include <boost/histogram/unlimited_storage.hpp>
 #include <deque>
 #include <initializer_list>
@@ -106,6 +107,18 @@ int main() {
     BOOST_TEST_TRAIT_TRUE((is_map_like<E>));
   }
 
+  // has_node_access
+  {
+    using A = storage_adaptor<std::map<std::size_t, double>>;
+    using B = storage_adaptor<std::unordered_map<std::size_t, double>>;
+    using C = storage_adaptor<std::vector<double>>;
+    using D = unlimited_storage<>;
+    BOOST_TEST_TRAIT_TRUE((has_node_access<A>));
+    BOOST_TEST_TRAIT_TRUE((has_node_access<B>));
+    BOOST_TEST_TRAIT_FALSE((has_node_access<C>));
+    BOOST_TEST_TRAIT_FALSE((has_node_access<D>));
+  }
+
   // is_axis
   {
     struct A {};

test/histogram_operators_test.cpp

@@ -12,9 +12,14 @@
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/ostream.hpp>
 #include <boost/throw_exception.hpp>
+#include <cmath>
+#include <cstddef>
+#include <limits>
+#include <map>
 #include <string>
 #include <type_traits>
 #include <vector>
+#include "allocator.hpp"
 #include "dummy_storage.hpp"
 #include "histogram.hpp"
 #include "ostream.hpp"
@@ -335,6 +340,126 @@ void run_tests() {
     }
   }
 
+  // arithmetic with sparse storage: map-backed histograms must not densify
+  {
+    using map_t = std::map<std::size_t, double, std::less<std::size_t>,
+                           tracing_allocator<std::pair<const std::size_t, double>>>;
+    tracing_allocator_db db;
+    tracing_allocator<char> alloc(db);
+    auto ia = axis::integer<int, axis::null_type, axis::option::none_t>(0, 100);
+
+    auto a = make_s(Tag(), map_t(alloc), ia);
+    auto b = make_s(Tag(), map_t(alloc), ia);
+    const auto baseline = db.first;
+    a(5); // causes one allocation
+    const auto node = db.first - baseline;
+    BOOST_TEST_GT(node, 0);
+    b(7); // causes one allocation
+    BOOST_TEST_EQ(db.first, baseline + 2 * node);
+
+    // dense references, all checks below compare against dense results
+    auto ra = make_s(Tag(), std::vector<double>(), ia);
+    auto rb = make_s(Tag(), std::vector<double>(), ia);
+    ra(5);
+    rb(7);
+
+    // map += map allocates one node, not one per cell
+    auto t0 = db.first;
+    a += b;
+    BOOST_TEST_EQ(db.first, t0 + node);
+    auto rab = ra;
+    rab += rb;
+    BOOST_TEST_EQ(a, rab);
+
+    // map -= map allocates one node holding a negative value
+    auto c = make_s(Tag(), map_t(alloc), ia);
+    t0 = db.first;
+    c -= b;
+    BOOST_TEST_EQ(db.first, t0 + node);
+    auto rc = make_s(Tag(), std::vector<double>(), ia);
+    rc -= rb;
+    BOOST_TEST_EQ(c, rc);
+
+    // self-addition doubles values without new allocations
+    auto d = a;
+    t0 = db.first;
+    d += d;
+    BOOST_TEST_EQ(db.first, t0);
+    auto rd = rab;
+    rd += rd;
+    BOOST_TEST_EQ(d, rd);
+
+    // map += dense allocates only for non-zero cells of rhs
+    auto e = make_s(Tag(), map_t(alloc), ia);
+    e(3);
+    t0 = db.first;
+    e += ra;
+    BOOST_TEST_EQ(db.first, t0 + node);
+    auto re = make_s(Tag(), std::vector<double>(), ia);
+    re(3);
+    re += ra;
+    BOOST_TEST_EQ(e, re);
+
+    // dense += map
+    auto f = make_s(Tag(), std::vector<double>(), ia);
+    f(2);
+    f += b;
+    auto rf = make_s(Tag(), std::vector<double>(), ia);
+    rf(2);
+    rf += rb;
+    BOOST_TEST_EQ(f, rf);
+
+    // default storage += map
+    auto g = make(Tag(), ia);
+    g(2);
+    g += b;
+    BOOST_TEST_EQ(g, rf);
+
+    // mixed value types: map of double += vector of int
+    auto h = make_s(Tag(), map_t(alloc), ia);
+    auto vi = make_s(Tag(), std::vector<int>(), ia);
+    vi(4);
+    t0 = db.first;
+    h += vi;
+    BOOST_TEST_EQ(db.first, t0 + node);
+    BOOST_TEST_EQ(h.at(4), 1);
+
+    // map *= map does not allocate
+    auto m = a; // a has entries at 5 and 7
+    t0 = db.first;
+    m *= b;
+    BOOST_TEST_EQ(db.first, t0);
+    BOOST_TEST_EQ(m.at(5), 0);
+    BOOST_TEST_EQ(m.at(7), 1);
+
+    // map /= map keeps dense semantics: 0/0 is NaN, 1/0 is inf
+    auto n = a;
+    n /= b;
+    BOOST_TEST_EQ(n.at(7), 1);
+    BOOST_TEST_EQ(n.at(5), std::numeric_limits<double>::infinity());
+    BOOST_TEST(std::isnan(static_cast<double>(n.at(0))));
+
+    // merging add with growing axes and map storage on both sides
+    using C = axis::category<int, use_default, axis::option::growth_t>;
+    auto p = make_s(Tag(), map_t(alloc), C{1, 2});
+    auto q = make_s(Tag(), map_t(alloc), C{2, 3});
+    p(1);
+    p(2);
+    q(2);
+    q(3);
+    q(3);
+    p += q;
+    auto rp = make_s(Tag(), std::vector<double>(), C{1, 2});
+    auto rq = make_s(Tag(), std::vector<double>(), C{2, 3});
+    rp(1);
+    rp(2);
+    rq(2);
+    rq(3);
+    rq(3);
+    rp += rq;
+    BOOST_TEST_EQ(p, rp);
+  }
+
   // bad operations
   {
     auto a = make(Tag(), axis::regular<>(2, 0, 4));
@@ -375,5 +500,20 @@ int main() {
   run_tests<static_tag>();
   run_tests<dynamic_tag>();
 
+  // histograms with different axes types and map-backed rhs
+  {
+    auto ia = axis::integer<int, axis::null_type, axis::option::none_t>(0, 3);
+    auto a = make_s(static_tag(), std::vector<double>(), ia);
+    auto b = make_s(dynamic_tag(), std::map<std::size_t, double>(), ia);
+    a(0);
+    b(1);
+    a += b;
+    BOOST_TEST_EQ(a.at(0), 1);
+    BOOST_TEST_EQ(a.at(1), 1);
+    a -= b;
+    BOOST_TEST_EQ(a.at(0), 1);
+    BOOST_TEST_EQ(a.at(1), 0);
+  }
+
   return boost::report_errors();
 }

test/storage_adaptor_test.cpp

@@ -305,6 +305,14 @@ int main() {
     BOOST_TEST_EQ(db.first, baseline + 3 * node);
     a[4] = 0; // causes one deallocation
     BOOST_TEST_EQ(db.first, baseline + 2 * node);
+    a[6] += 0.0; // no allocation, result is indistinguishable from empty cell
+    a[6] -= 0.0; // no allocation
+    BOOST_TEST_EQ(db.first, baseline + 2 * node);
+    BOOST_TEST_EQ(a[6], 0);
+    a[6] += 1.0; // causes one allocation
+    a[6] -= 1.0; // node stays, in-place zeroing does not erase
+    BOOST_TEST_EQ(db.first, baseline + 3 * node);
+    BOOST_TEST_EQ(a[6], 0);
 
     auto b = storage_adaptor<std::vector<int>>();
     b.reset(5);
@@ -314,5 +322,23 @@ int main() {
     BOOST_TEST_EQ(db.first, baseline + node);
   }
 
+  // adding a zero accumulator to an empty cell does not allocate
+  {
+    using ws_t = accumulators::weighted_sum<double>;
+    tracing_allocator_db db;
+    tracing_allocator<char> alloc(db);
+    using map_t = std::map<std::size_t, ws_t, std::less<std::size_t>,
+                           tracing_allocator<std::pair<const std::size_t, ws_t>>>;
+    auto a = storage_adaptor<map_t>(alloc);
+    const auto baseline = db.first;
+    a.reset(10);
+    a[1] += ws_t();
+    BOOST_TEST_EQ(db.first, baseline);
+    BOOST_TEST(a[1] == ws_t());
+    a[1] += ws_t(1, 1); // causes one allocation
+    BOOST_TEST_GT(db.first, baseline);
+    BOOST_TEST(a[1] == ws_t(1, 1));
+  }
+
   return boost::report_errors();
 }