Unsafe access to accumulator data

[HDembinski]

• Test this with downstream boost-histogram

[HDembinski]

@henryiii Before I change all accumulators in this way, we should try first whether this change works the way we expect in boost-histogram. I can try to make the changes to test this in boost-histogram.

[henryiii]

Can this be rebased on top of current develop?

[HDembinski]

@henryiii Rebased on develop as you requested

[henryiii]

Thanks, I'll work on this further soon - might be after Python 3.10 / pybind11 2.8 / Windows 11 though, all to be released next Monday. (I don't think Windows 11 affects me, but recently noticed it had the same release date as Python 3.10) First initial try showed it will take some work, but didn't get far enough for good feedback yet.

[henryiii]

FYI, I did attempt this a little while ago, but it will take a bit of work to properly implement. Boost-histogram has to know the memory layout exactly, since it can cast between a block of Python memory to C++ without copy. Convincing pybind11 of this without access to the actual structure just requires manual work for things normally covered by macros. I'll try to get back to this in ~1 week.

[henryiii]

I tried again, with AI this time. I can make it work. I have to use the internal PYBIND11_FIELD_DESCRIPTOR_EX since PYBIND11_NUMPY_DTYPE requires public members, but that's not too risky as we don't change this part of pybind11 much.

🤖 AI text below 🤖

It works. Full boost-histogram suite passes (1134 tests) using the real boost::histogram::accumulators::mean through unsafe_access::accumulator_data.

boostorg/histogram (unsafe_access_accumulators, PR #323): rebased onto develop (one trivial include conflict in test/accumulators_mean_test.cpp); all 89 CMake tests pass. Not pushed — it's Hans's PR branch, so force-pushing is your call.

boost-histogram (new branch henryiii/try/unsafe-access-mean, commit 336c406, extern/histogram pointed at the rebased commit):

• bh_python/accumulators/mean.hpp is now an alias accumulators::mean<T> = bh::accumulators::mean<T> plus a mean_data() helper, a NaN-guarded variance(), and save/load overloads that keep the old positional pickle layout (Boost's own serialize would read a size_t for version-0 pickles).
• register_accumulators.cpp: dtype registered by hand with npy_format_descriptor<mean>::register_dtype + PYBIND11_FIELD_DESCRIPTOR_EX(mean::data_type, sum_, "count") etc. — PYBIND11_NUMPY_DTYPE can't be used because it needs public members. A static_assert checks sizeof(mean) == sizeof(data_type). Properties/__getitem__/__setitem__/_make go through mean_data().
• ostream.hpp: shift_to_string overload for the Boost mean so repr stays Mean(count=…, value=…, variance=…); the histogram repr also unchanged.

Things worth noting for the upstream design:

1. The manual dtype is ~10 lines and would be the same for every accumulator, so the approach scales fine to weighted_sum/weighted_mean once Hans converts them (he said he'd wait for this trial).
2. PYBIND11_FIELD_DESCRIPTOR_EX is an internal pybind11 macro; if that's a concern, it's trivial to write the field_descriptor initializer by hand.
3. Boost's mean::variance() returns -0/nan unguarded for count ≤ 1; the binding keeps the NaN behavior locally. Could move the guard upstream.
4. extern/histogram points at a local-only commit, so the boost-histogram branch can't build in CI until Unsafe access to accumulator data #323 is pushed/merged.

boost-histogram diff:

diff --git a/extern/histogram b/extern/histogram
index 5e9b7b1..dcf0c7a 160000
--- a/extern/histogram
+++ b/extern/histogram
@@ -1 +1 @@
-Subproject commit 5e9b7b162e8a2801b77ccf199dedea88fa44325c
+Subproject commit dcf0c7a613ab98f85d3ff9c8e2958ea5767acbde
diff --git a/include/bh_python/accumulators/mean.hpp b/include/bh_python/accumulators/mean.hpp
index bded524..d93d94f 100644
--- a/include/bh_python/accumulators/mean.hpp
+++ b/include/bh_python/accumulators/mean.hpp
@@ -4,109 +4,52 @@
 // (See accompanying file LICENSE_1_0.txt
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 //
-// Based on boost/histogram/accumulators/mean.hpp
-// Changes:
-//  * Internal values are public for access from Python
-//  * A special constructor added for construction from Python
+// Uses boost::histogram::accumulators::mean directly; the internal data is
+// reached through boost::histogram::unsafe_access::accumulator_data.
 
 #pragma once
 
 #include <boost/core/nvp.hpp>
-#include <boost/histogram/weight.hpp>
+#include <boost/histogram/accumulators/mean.hpp>
+#include <boost/histogram/unsafe_access.hpp>
 
 #include <limits>
+#include <type_traits>
 
 namespace accumulators {
 
-/** Calculates mean and variance of sample.
-
-  Uses Welford's incremental algorithm to improve the numerical
-  stability of mean and variance computation.
-*/
 template <class ValueType>
-struct mean {
-    using value_type      = ValueType;
-    using const_reference = const value_type&;
-
-    mean() = default;
-
-    mean(const value_type& n,
-         const value_type& mean,
-         const value_type& variance) noexcept
-        : count(n)
-        , value(mean)
-        , _sum_of_deltas_squared(variance * (n - 1)) {}
-
-    mean(const value_type& sum,
-         const value_type& mean,
-         const value_type& _sum_of_deltas_squared,
-         bool /* Tag to trigger python internal constructor */)
-        : count(sum)
-        , value(mean)
-        , _sum_of_deltas_squared(_sum_of_deltas_squared) {}
-
-    void operator()(const value_type& x) noexcept {
-        count += static_cast<value_type>(1);
-        const auto delta = x - value;
-        value += delta / count;
-        _sum_of_deltas_squared += delta * (x - value);
-    }
-
-    void operator()(const boost::histogram::weight_type<value_type>& w,
-                    const value_type& x) noexcept {
-        count += w.value;
-        const auto delta = x - value;
-        value += w.value * delta / count;
-        _sum_of_deltas_squared += w.value * delta * (x - value);
-    }
-
-    mean& operator+=(const mean& rhs) noexcept {
-        if(rhs.count == 0)
-            return *this;
-
-        const auto mu1 = value;
-        const auto mu2 = rhs.value;
-        const auto n1  = count;
-        const auto n2  = rhs.count;
+using mean = boost::histogram::accumulators::mean<ValueType>;
 
-        count += rhs.count;
-        value = (n1 * mu1 + n2 * mu2) / count;
-        _sum_of_deltas_squared += rhs._sum_of_deltas_squared;
-        _sum_of_deltas_squared
-            += n1 * (value - mu1) * (value - mu1) + n2 * (value - mu2) * (value - mu2);
+template <class T>
+decltype(auto) mean_data(T& m) {
+    return boost::histogram::unsafe_access::accumulator_data(m);
+}
 
-        return *this;
-    }
-
-    mean& operator*=(const value_type& s) noexcept {
-        value *= s;
-        _sum_of_deltas_squared *= s * s;
-        return *this;
-    }
-
-    bool operator==(const mean& rhs) const noexcept {
-        return count == rhs.count && value == rhs.value
-               && _sum_of_deltas_squared == rhs._sum_of_deltas_squared;
-    }
-
-    bool operator!=(const mean& rhs) const noexcept { return !operator==(rhs); }
-
-    value_type variance() const noexcept {
-        if(count <= 1)
-            return std::numeric_limits<value_type>::quiet_NaN();
-        return _sum_of_deltas_squared / (count - 1);
-    }
-
-    template <class Archive>
-    void serialize(Archive& ar, unsigned) {
-        ar& boost::make_nvp("count", count);
-        ar& boost::make_nvp("value", value);
-        ar& boost::make_nvp("_sum_of_deltas_squared", _sum_of_deltas_squared);
-    }
-
-    value_type count{};
-    value_type value{};
-    value_type _sum_of_deltas_squared{};
-};
+/// Boost's variance() divides by count-1 without a guard; keep the NaN result
+template <class T>
+T variance(const mean<T>& m) {
+    if(m.count() <= 1)
+        return std::numeric_limits<T>::quiet_NaN();
+    return m.variance();
+}
 
 } // namespace accumulators
+
+// Positional pickle layout matches the old bh_python mean (count, value, sods).
+// Boost's own serialize reads a size_t count for version 0, which old pickles lack.
+template <class Archive, class T>
+void save(Archive& ar, const accumulators::mean<T>& m, unsigned /* version */) {
+    const auto& d = accumulators::mean_data(m);
+    ar& boost::make_nvp("count", d.sum_);
+    ar& boost::make_nvp("value", d.mean_);
+    ar& boost::make_nvp("_sum_of_deltas_squared", d.sum_of_deltas_squared_);
+}
+
+template <class Archive, class T>
+void load(Archive& ar, accumulators::mean<T>& m, unsigned /* version */) {
+    auto& d = accumulators::mean_data(m);
+    ar& boost::make_nvp("count", d.sum_);
+    ar& boost::make_nvp("value", d.mean_);
+    ar& boost::make_nvp("_sum_of_deltas_squared", d.sum_of_deltas_squared_);
+}
diff --git a/include/bh_python/accumulators/ostream.hpp b/include/bh_python/accumulators/ostream.hpp
index c182adf..6458fab 100644
--- a/include/bh_python/accumulators/ostream.hpp
+++ b/include/bh_python/accumulators/ostream.hpp
@@ -59,15 +59,6 @@ std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>&
     return handle_nonzero_width(os, x);
 }
 
-template <class CharT, class Traits, class W>
-std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
-                                              const mean<W>& x) {
-    if(os.width() == 0)
-        return os << "count=" << x.count << ", value=" << x.value
-                  << ", variance=" << x.variance();
-    return handle_nonzero_width(os, x);
-}
-
 template <class CharT, class Traits, class W>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
                                               const weighted_mean<W>& x) {
@@ -95,3 +86,12 @@ shift_to_string(const ::boost::histogram::accumulators::sum<double>& x) {
     out << x.large_part() << " + " << x.small_part();
     return out.str();
 }
+
+// boost::histogram's operator<< for mean prints "mean(...)"; Python wants fields.
+inline std::string
+shift_to_string(const ::boost::histogram::accumulators::mean<double>& x) {
+    std::ostringstream out;
+    out << "count=" << x.count() << ", value=" << x.value()
+        << ", variance=" << accumulators::variance(x);
+    return out.str();
+}
diff --git a/src/register_accumulators.cpp b/src/register_accumulators.cpp
index 1c634e3..c4678c2 100644
--- a/src/register_accumulators.cpp
+++ b/src/register_accumulators.cpp
@@ -14,6 +14,9 @@
 #include <boost/histogram/accumulators/sum.hpp>
 #include <pybind11/operators.h>
 
+#include <type_traits>
+#include <vector>
+
 namespace {
 /// The mean fill can be implemented once. (sum fill varies slightly)
 template <class T>
@@ -267,8 +270,29 @@ void register_accumulators(py::module& accumulators) {
 
         ;
 
-    using mean = accumulators::mean<double>;
-    PYBIND11_NUMPY_DTYPE(mean, count, value, _sum_of_deltas_squared);
+    using mean      = accumulators::mean<double>;
+    using mean_data = mean::data_type;
+    static_assert(std::is_standard_layout<mean>::value
+                      && sizeof(mean) == sizeof(mean_data),
+                  "mean must wrap data_type with no padding");
+
+    // Boost's members are private, so PYBIND11_NUMPY_DTYPE cannot be used; the dtype
+    // is built from data_type instead, with the public Python field names.
+    py::detail::npy_format_descriptor<mean>::register_dtype(
+        std::vector<py::detail::field_descriptor>{
+            PYBIND11_FIELD_DESCRIPTOR_EX(mean_data, sum_, "count"),
+            PYBIND11_FIELD_DESCRIPTOR_EX(mean_data, mean_, "value"),
+            PYBIND11_FIELD_DESCRIPTOR_EX(
+                mean_data, sum_of_deltas_squared_, "_sum_of_deltas_squared")});
+
+    auto make_mean_raw = [](double count, double value, double sods) {
+        mean m;
+        auto& d                  = accumulators::mean_data(m);
+        d.sum_                   = count;
+        d.mean_                  = value;
+        d.sum_of_deltas_squared_ = sods;
+        return m;
+    };
 
     register_accumulator<mean>(accumulators, "Mean", py::buffer_protocol())
         .def_buffer(make_buffer<mean>())
@@ -278,11 +302,16 @@ void register_accumulators(py::module& accumulators) {
              "value"_a,
              "variance"_a)
 
-        .def_readonly("count", &mean::count)
-        .def_readonly("value", &mean::value)
-        .def_readonly("_sum_of_deltas_squared", &mean::_sum_of_deltas_squared)
+        .def_property_readonly("count", [](const mean& self) { return self.count(); })
+        .def_property_readonly("value", [](const mean& self) { return self.value(); })
+        .def_property_readonly(
+            "_sum_of_deltas_squared",
+            [](const mean& self) {
+                return accumulators::mean_data(self).sum_of_deltas_squared_;
+            })
 
-        .def_property_readonly("variance", &mean::variance)
+        .def_property_readonly(
+            "variance", [](const mean& self) { return accumulators::variance(self); })
 
         .def("__call__",
              make_mean_call<mean>(),
@@ -298,11 +327,7 @@ void register_accumulators(py::module& accumulators) {
              "weight"_a = py::none(),
              "Fill the accumulator with values. Optional weight parameter.")
 
-        .def_static(
-            "_make",
-            py::vectorize([](const double& a, const double& b, const double& c) {
-                return mean(a, b, c, true);
-            }))
+        .def_static("_make", py::vectorize(make_mean_raw))
 
         .def_static(
             "_array",
@@ -312,24 +337,26 @@ void register_accumulators(py::module& accumulators) {
 
         .def("__getitem__",
              [](const mean& self, const py::str& key) {
+                 const auto& d = accumulators::mean_data(self);
                  if(key.equal(py::str("count")))
-                     return self.count;
+                     return d.sum_;
                  if(key.equal(py::str("value")))
-                     return self.value;
+                     return d.mean_;
                  if(key.equal(py::str("_sum_of_deltas_squared")))
-                     return self._sum_of_deltas_squared;
+                     return d.sum_of_deltas_squared_;
                  throw py::key_error(
                      py::str("{0} not one of count, value, _sum_of_deltas_squared")
                          .format(key));
              })
         .def("__setitem__",
              [](mean& self, const py::str& key, double value) {
+                 auto& d = accumulators::mean_data(self);
                  if(key.equal(py::str("count")))
-                     self.count = value;
+                     d.sum_ = value;
                  else if(key.equal(py::str("value")))
-                     self.value = value;
+                     d.mean_ = value;
                  else if(key.equal(py::str("_sum_of_deltas_squared")))
-                     self._sum_of_deltas_squared = value;
+                     d.sum_of_deltas_squared_ = value;
                  else
                      throw py::key_error(
                          py::str("{0} not one of count, value, _sum_of_deltas_squared")