Unsafe access to accumulator data

benchmark/histogram_parallel_filling.cpp

@@ -5,7 +5,7 @@
 // or copy at http://www.boost.org/LICENSE_1_0.txt)
 
 #include <benchmark/benchmark.h>
-#include <boost/histogram/accumulators/thread_safe.hpp>
+#include <boost/histogram/accumulators/count.hpp>
 #include <boost/histogram/axis/regular.hpp>
 #include <boost/histogram/histogram.hpp>
 #include <boost/histogram/make_histogram.hpp>
@@ -29,7 +29,7 @@ using namespace boost::histogram;
 using namespace std::chrono_literals;
 
 using DS = dense_storage<unsigned>;
-using DSTS = dense_storage<accumulators::thread_safe<unsigned>>;
+using DSTS = dense_storage<accumulators::count<unsigned, true>>;
 
 static void NoThreads(benchmark::State& state) {
   std::default_random_engine gen(1);

include/boost/histogram/accumulators/mean.hpp

@@ -30,36 +30,41 @@ class mean {
   using value_type = ValueType;
   using const_reference = const value_type&;
 
+  struct data_type {
+    value_type sum_;
+    value_type mean_;
+    value_type sum_of_deltas_squared_;
+  };
+
   mean() = default;
 
   /// Allow implicit conversion from mean<T>.
   template <class T>
-  mean(const mean<T>& o) noexcept
-      : sum_{o.sum_}, mean_{o.mean_}, sum_of_deltas_squared_{o.sum_of_deltas_squared_} {}
+  mean(const mean<T>& o) noexcept : data_{o.data_} {}
 
   /// Initialize to external count, mean, and variance.
   mean(const_reference n, const_reference mean, const_reference variance) noexcept
-      : sum_(n), mean_(mean), sum_of_deltas_squared_(variance * (n - 1)) {}
+      : data_{n, mean, variance * (n - 1)} {}
 
   /// Insert sample x.
   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_);
+    data_.sum_ += static_cast<value_type>(1);
+    const auto delta = x - data_.mean_;
+    data_.mean_ += delta / data_.sum_;
+    data_.sum_of_deltas_squared_ += delta * (x - data_.mean_);
   }
 
   /// Insert sample x with weight w.
   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_;
-    sum_of_deltas_squared_ += w.value * delta * (x - mean_);
+    data_.sum_ += w.value;
+    const auto delta = x - data_.mean_;
+    data_.mean_ += w.value * delta / data_.sum_;
+    data_.sum_of_deltas_squared_ += w.value * delta * (x - data_.mean_);
   }
 
   /// Add another mean accumulator.
   mean& operator+=(const mean& rhs) noexcept {
-    if (rhs.sum_ == 0) return *this;
+    if (rhs.data_.sum_ == 0) return *this;
 
     /*
       sum_of_deltas_squared
@@ -75,20 +80,20 @@ class mean {
 
       Putting it together:
       sum_of_deltas_squared
-        = sum_of_deltas_squared_1 + n1 (mu1 - mu))^2
-        + sum_of_deltas_squared_2 + n2 (mu2 - mu))^2
+        = sum_of_deltas_squared_1 + n1 (mu - mu1))^2
+        + sum_of_deltas_squared_2 + n2 (mu - mu2))^2
     */
 
-    const auto n1 = sum_;
-    const auto mu1 = mean_;
-    const auto n2 = rhs.sum_;
-    const auto mu2 = rhs.mean_;
+    const auto n1 = data_.sum_;
+    const auto mu1 = data_.mean_;
+    const auto n2 = rhs.data_.sum_;
+    const auto mu2 = rhs.data_.mean_;
 
-    sum_ += rhs.sum_;
-    mean_ = (n1 * mu1 + n2 * mu2) / sum_;
-    sum_of_deltas_squared_ += rhs.sum_of_deltas_squared_;
-    sum_of_deltas_squared_ += n1 * detail::square(mean_ - mu1);
-    sum_of_deltas_squared_ += n2 * detail::square(mean_ - mu2);
+    data_.sum_ += rhs.data_.sum_;
+    data_.mean_ = (n1 * mu1 + n2 * mu2) / data_.sum_;
+    data_.sum_of_deltas_squared_ += rhs.data_.sum_of_deltas_squared_;
+    data_.sum_of_deltas_squared_ += n1 * detail::square(data_.mean_ - mu1);
+    data_.sum_of_deltas_squared_ += n2 * detail::square(data_.mean_ - mu2);
 
     return *this;
   }
@@ -98,14 +103,14 @@ class mean {
    This acts as if all samples were scaled by the value.
   */
   mean& operator*=(const_reference s) noexcept {
-    mean_ *= s;
-    sum_of_deltas_squared_ *= s * s;
+    data_.mean_ *= s;
+    data_.sum_of_deltas_squared_ *= s * s;
     return *this;
   }
 
   bool operator==(const mean& rhs) const noexcept {
-    return sum_ == rhs.sum_ && mean_ == rhs.mean_ &&
-           sum_of_deltas_squared_ == rhs.sum_of_deltas_squared_;
+    return data_.sum_ == rhs.data_.sum_ && data_.mean_ == rhs.data_.mean_ &&
+           data_.sum_of_deltas_squared_ == rhs.data_.sum_of_deltas_squared_;
   }
 
   bool operator!=(const mean& rhs) const noexcept { return !operator==(rhs); }
@@ -116,38 +121,40 @@ class mean {
     see documentation of value() and variance(). count() can be used to compute
     the variance of the mean by dividing variance() by count().
   */
-  const_reference count() const noexcept { return sum_; }
+  const_reference count() const noexcept { return data_.sum_; }
 
   /** Return mean value of accumulated samples.
 
     The result is undefined, if `count() < 1`.
   */
-  const_reference value() const noexcept { return mean_; }
+  const_reference value() const noexcept { return data_.mean_; }
 
   /** Return variance of accumulated samples.
 
     The result is undefined, if `count() < 2`.
   */
-  value_type variance() const noexcept { return sum_of_deltas_squared_ / (sum_ - 1); }
+  value_type variance() const noexcept {
+    return data_.sum_of_deltas_squared_ / (data_.sum_ - 1);
+  }
 
   template <class Archive>
   void serialize(Archive& ar, unsigned version) {
     if (version == 0) {
       // read only
       std::size_t sum;
       ar& make_nvp("sum", sum);
-      sum_ = static_cast<value_type>(sum);
+      data_.sum_ = static_cast<value_type>(sum);
     } else {
-      ar& make_nvp("sum", sum_);
+      ar& make_nvp("sum", data_.sum_);
     }
-    ar& make_nvp("mean", mean_);
-    ar& make_nvp("sum_of_deltas_squared", sum_of_deltas_squared_);
+    ar& make_nvp("mean", data_.mean_);
+    ar& make_nvp("sum_of_deltas_squared", data_.sum_of_deltas_squared_);
   }
 
 private:
-  value_type sum_{};
-  value_type mean_{};
-  value_type sum_of_deltas_squared_{};
+  data_type data_{0, 0, 0};
+
+  friend struct ::boost::histogram::unsafe_access;
 };
 
 } // namespace accumulators

include/boost/histogram/accumulators/ostream.hpp

@@ -63,7 +63,8 @@ std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>&
 template <class CharT, class Traits, class U>
 std::basic_ostream<CharT, Traits>& operator<<(std::basic_ostream<CharT, Traits>& os,
                                               const sum<U>& x) {
-  if (os.width() == 0) return os << "sum(" << x.large() << " + " << x.small() << ")";
+  if (os.width() == 0)
+    return os << "sum(" << x.large_part() << " + " << x.small_part() << ")";
   return detail::handle_nonzero_width(os, x);
 }
 

include/boost/histogram/accumulators/sum.hpp

@@ -43,11 +43,11 @@ class sum {
 
   /// Allow implicit conversion from sum<T>
   template <class T>
-  sum(const sum<T>& s) noexcept : sum(s.large(), s.small()) {}
+  sum(const sum<T>& s) noexcept : sum(s.large_part(), s.small_part()) {}
 
   /// Initialize sum explicitly with large and small parts
-  sum(const_reference large, const_reference small) noexcept
-      : large_(large), small_(small) {}
+  sum(const_reference large_part, const_reference small_part) noexcept
+      : large_(large_part), small_(small_part) {}
 
   /// Increment sum by one
   sum& operator++() noexcept { return operator+=(1); }
@@ -96,10 +96,10 @@ class sum {
   value_type value() const noexcept { return large_ + small_; }
 
   /// Return large part of the sum.
-  const_reference large() const noexcept { return large_; }
+  const_reference large_part() const noexcept { return large_; }
 
   /// Return small part of the sum.
-  const_reference small() const noexcept { return small_; }
+  const_reference small_part() const noexcept { return small_; }
 
   // lossy conversion to value type must be explicit
   explicit operator value_type() const noexcept { return value(); }
@@ -156,6 +156,25 @@ class sum {
 
   // end: extra operators
 
+  // windows.h illegially uses `#define small char` which breaks this now deprecated API
+#if !defined(small)
+
+  /// Return large part of the sum.
+  [[deprecated("use large_part() instead; "
+               "large() will be removed in boost-1.80")]] const_reference
+  large() const noexcept {
+    return large_;
+  }
+
+  /// Return small part of the sum.
+  [[deprecated("use small_part() instead; "
+               "small() will be removed in boost-1.80")]] const_reference
+  small() const noexcept {
+    return small_;
+  }
+
+#endif
+
 private:
   value_type large_{};
   value_type small_{};

include/boost/histogram/algorithm/reduce.hpp

@@ -33,7 +33,8 @@ namespace algorithm {
 */
 using reduce_command = detail::reduce_command;
 
-using reduce_option [[deprecated("use reduce_command instead")]] =
+using reduce_option [[deprecated("use reduce_command instead; "
+                                 "reduce_option will be removed in boost-1.80")]] =
     reduce_command; ///< deprecated
 
 /** Shrink command to be used in `reduce`.

include/boost/histogram/axis/traits.hpp

@@ -194,7 +194,9 @@ template <class Axis>
 using get_options = decltype(detail::traits_options<Axis>(detail::priority<2>{}));
 
 template <class Axis>
-using static_options [[deprecated("use get_options instead")]] = get_options<Axis>;
+using static_options [[deprecated("use get_options instead; "
+                                  "static_options will be removed in boost-1.80")]] =
+    get_options<Axis>;
 
 #else
 struct get_options;
@@ -223,7 +225,10 @@ template <class Axis>
 using is_inclusive = decltype(detail::traits_is_inclusive<Axis>(detail::priority<1>{}));
 
 template <class Axis>
-using static_is_inclusive [[deprecated("use is_inclusive instead")]] = is_inclusive<Axis>;
+using static_is_inclusive
+    [[deprecated("use is_inclusive instead; "
+                 "static_is_inclusive will be removed in boost-1.80")]] =
+        is_inclusive<Axis>;
 
 #else
 struct is_inclusive;

include/boost/histogram/detail/fill_n.hpp

@@ -242,6 +242,7 @@ void fill_n_1(const std::size_t offset, S& storage, A& axes, const std::size_t v
   for_each_axis(axes,
                 [&](const auto& ax) { all_inclusive &= axis::traits::inclusive(ax); });
   if (axes_rank(axes) == 1) {
+    // Optimization: benchmark shows that this makes filling dynamic 1D histogram faster
     axis::visit(
         [&](auto& ax) {
           std::tuple<decltype(ax)> axes{ax};

include/boost/histogram/indexed.hpp

@@ -64,7 +64,9 @@ class BOOST_ATTRIBUTE_NODISCARD indexed_range {
   using value_type = typename std::iterator_traits<value_iterator>::value_type;
 
   class iterator;
-  using range_iterator [[deprecated("use iterator instead")]] = iterator; ///< deprecated
+  using range_iterator [[deprecated("use iterator instead; "
+                                    "range_iterator will be removed in boost-1.80")]] =
+      iterator; ///< deprecated
 
   /** Lightweight view to access value and index of current cell.
 
@@ -84,7 +86,8 @@ class BOOST_ATTRIBUTE_NODISCARD indexed_range {
 
     public:
       using const_reference = const axis::index_type&;
-      using reference [[deprecated("use const_reference instead")]] =
+      using reference [[deprecated("use const_reference instead; "
+                                   "reference will be removed in boost-1.80")]] =
           const_reference; ///< deprecated
 
       /// implementation detail

include/boost/histogram/unsafe_access.hpp

@@ -97,8 +97,8 @@ struct unsafe_access {
     Get buffer of unlimited_storage.
     @param storage instance of unlimited_storage.
   */
-  template <class Allocator>
-  static constexpr auto& unlimited_storage_buffer(unlimited_storage<Allocator>& storage) {
+  template <class T>
+  static constexpr auto& unlimited_storage_buffer(T& storage) {
     return storage.buffer_;
   }
 
@@ -107,8 +107,17 @@ struct unsafe_access {
     @param storage instance of storage_adaptor.
   */
   template <class T>
-  static constexpr auto& storage_adaptor_impl(storage_adaptor<T>& storage) {
-    return static_cast<typename storage_adaptor<T>::impl_type&>(storage);
+  static constexpr auto& storage_adaptor_impl(T& storage) {
+    return static_cast<typename T::impl_type&>(storage);
+  }
+
+  /**
+    Get internal data of accumulator.
+    @param obj instance of accumulator.
+  */
+  template <class T>
+  static constexpr auto& accumulator_data(T& m) {
+    return m.data_;
   }
 };
 

test/Jamfile

@@ -157,11 +157,14 @@ alias libserial :
 # for builds without optional boost dependencies
 alias minimal : cxx14 cxx17 failure threading ;
 
+# for builds with optional boost dependencies
+alias optional_boost : accumulators range units serialization ;
+
 # all tests
-alias all : minimal not_windows odr accumulators range units serialization ;
+alias all : minimal not_windows odr optional_boost ;
 
 # all except "failure", because it is distracting during development
-alias develop : cxx14 cxx17 threading not_windows odr accumulators range units serialization ;
+alias develop : odr cxx14 cxx17 threading not_windows optional_boost ;
 
 explicit minimal ;
 explicit all ;
@@ -176,3 +179,4 @@ explicit range ;
 explicit units ;
 explicit serialization ;
 explicit libserial ;
+explicit optional_boost ;

test/accumulators_mean_test.cpp

@@ -7,6 +7,7 @@
 #include <boost/core/lightweight_test.hpp>
 #include <boost/histogram/accumulators/mean.hpp>
 #include <boost/histogram/accumulators/ostream.hpp>
+#include <boost/histogram/unsafe_access.hpp>
 #include <boost/histogram/weight.hpp>
 #include <sstream>
 #include "is_close.hpp"
@@ -110,5 +111,16 @@ int main() {
     BOOST_TEST_IS_CLOSE(a.variance(), b.variance(), 1e-3);
   }
 
+  // unsafe_access
+  {
+    m_t a;
+    a(1);
+    a(2);
+
+    BOOST_TEST_EQ(a.count(), 2);
+    unsafe_access::accumulator_data(a).sum_ = 1;
+    BOOST_TEST_EQ(a.count(), 1);
+  }
+
   return boost::report_errors();
 }