[hist] Implement initial RHistEngine

hist/histv7/benchmark/CMakeLists.txt

@@ -1,5 +1,8 @@
 add_executable(hist_benchmark_axes hist_benchmark_axes.cxx)
 target_link_libraries(hist_benchmark_axes ROOTHist benchmark::benchmark)
 
+add_executable(hist_benchmark_engine hist_benchmark_engine.cxx)
+target_link_libraries(hist_benchmark_engine ROOTHist benchmark::benchmark)
+
 add_executable(hist_benchmark_regular hist_benchmark_regular.cxx)
 target_link_libraries(hist_benchmark_regular ROOTHist benchmark::benchmark)

hist/histv7/benchmark/hist_benchmark_engine.cxx

@@ -0,0 +1,69 @@
+#include <ROOT/RHistEngine.hxx>
+#include <ROOT/RRegularAxis.hxx>
+
+#include <benchmark/benchmark.h>
+
+#include <random>
+#include <vector>
+
+struct RHistEngine_int_Regular1 : public benchmark::Fixture {
+   // The objects are stored and constructed in the fixture to avoid compiler optimizations in the benchmark body taking
+   // advantage of the (constant) constructor parameters.
+   ROOT::Experimental::RRegularAxis axis{20, 0.0, 1.0};
+   ROOT::Experimental::RHistEngine<int> engine{{axis}};
+   std::vector<double> fNumbers;
+
+   // Avoid GCC warning
+   using benchmark::Fixture::SetUp;
+   void SetUp(benchmark::State &state) final
+   {
+      std::mt19937 gen;
+      std::uniform_real_distribution<> dis;
+      fNumbers.resize(state.range(0));
+      for (std::size_t i = 0; i < fNumbers.size(); i++) {
+         fNumbers[i] = dis(gen);
+      }
+   }
+};
+
+BENCHMARK_DEFINE_F(RHistEngine_int_Regular1, Fill)(benchmark::State &state)
+{
+   for (auto _ : state) {
+      for (double number : fNumbers) {
+         engine.Fill(number);
+      }
+   }
+}
+BENCHMARK_REGISTER_F(RHistEngine_int_Regular1, Fill)->Range(0, 32768);
+
+struct RHistEngine_int_Regular2 : public benchmark::Fixture {
+   // The objects are stored and constructed in the fixture to avoid compiler optimizations in the benchmark body taking
+   // advantage of the (constant) constructor parameters.
+   ROOT::Experimental::RRegularAxis axis{20, 0.0, 1.0};
+   ROOT::Experimental::RHistEngine<int> engine{{axis, axis}};
+   std::vector<double> fNumbers;
+
+   // Avoid GCC warning
+   using benchmark::Fixture::SetUp;
+   void SetUp(benchmark::State &state) final
+   {
+      std::mt19937 gen;
+      std::uniform_real_distribution<> dis;
+      fNumbers.resize(2 * state.range(0));
+      for (std::size_t i = 0; i < fNumbers.size(); i++) {
+         fNumbers[i] = dis(gen);
+      }
+   }
+};
+
+BENCHMARK_DEFINE_F(RHistEngine_int_Regular2, Fill)(benchmark::State &state)
+{
+   for (auto _ : state) {
+      for (std::size_t i = 0; i < fNumbers.size(); i += 2) {
+         engine.Fill(fNumbers[i], fNumbers[i + 1]);
+      }
+   }
+}
+BENCHMARK_REGISTER_F(RHistEngine_int_Regular2, Fill)->Range(0, 32768);
+
+BENCHMARK_MAIN();

hist/histv7/doc/DesignImplementation.md

@@ -4,7 +4,7 @@ This document describes key design decisions and implementation choices.
 
 ## Templating
 
-Classes are only templated if required for data members, in particular the bin content type `T`.
+Classes are only templated if required for data members, in particular the `BinContentType`.
 We use member function templates to accept variable number of arguments (see also below).
 Classes are **not** templated to improve performance, in particular not on the axis type(s).
 This avoids an explosion of types and simplifies serialization.
@@ -28,7 +28,7 @@ Many member functions have two overloads: one accepting a function parameter pac
 ### Arguments with Different Types
 
 Functions that take arguments with different types expect a `std::tuple`.
-An example is `template <typename A...> void Fill(const std::tuple<A...> &args)`.
+An example is `template <typename... A> void Fill(const std::tuple<A...> &args)`.
 
 For user-convenience, a variadic function template forwards to the `std::tuple` overload:
 ```cpp
@@ -41,13 +41,13 @@ This will forward the arguments as references, so no copy-constructors are calle
 ### Arguments with Same Type
 
 In this case, the function has a `std::size_t N` template argument and accepts a `std::array`.
-An example is `template <std::size_t N> const T &GetBinContent(const std::array<RBinIndex, N> &args)`
+An example is `template <std::size_t N> const BinContentType &GetBinContent(const std::array<RBinIndex, N> &indices)`
 
 For user-convenience, a variadic function template forwards to the `std::array` overload:
 ```cpp
-template <typename... A> const T &GetBinContent(const A &...args) {
-   std::array<RBinIndex, sizeof...(A)> a{args...};
-   return GetBinContent(a);
+template <typename... A> const BinContentType &GetBinContent(const A &...args) {
+   std::array<RBinIndex, sizeof...(A)> indices{args...};
+   return GetBinContent(indices);
 }
 ```
 This will copy the arguments, which is fine in this case because `RBinIndex` is small (see below).
@@ -58,7 +58,7 @@ Special arguments are passed last.
 Examples include
 ```cpp
 template <typename... A> void Fill(const std::tuple<A...> &args, RWeight w);
-template <std::size_t N> void SetBinContent(const std::array<RBinIndex, N> &args, const T &content);
+template <std::size_t N> void SetBinContent(const std::array<RBinIndex, N> &indices, const BinContentType &content);
 ```
 The same works for the variadic function templates that will check the type of the last argument.
 

hist/histv7/headers.cmake

@@ -2,6 +2,7 @@ set(histv7_headers
     ROOT/RAxes.hxx
     ROOT/RBinIndex.hxx
     ROOT/RBinIndexRange.hxx
+    ROOT/RHistEngine.hxx
     ROOT/RLinearizedIndex.hxx
     ROOT/RRegularAxis.hxx
     ROOT/RVariableBinAxis.hxx

hist/histv7/inc/LinkDef.h

@@ -1,3 +1,11 @@
+// For RHistEngine, we request dictionaries for the most commonly used bin content types. This results in proper error
+// messages when trying to stream. Other instantiations will be caught by the RAxes member.
+#pragma link C++ class ROOT::Experimental::RHistEngine<int>-;
+#pragma link C++ class ROOT::Experimental::RHistEngine<long>-;
+#pragma link C++ class ROOT::Experimental::RHistEngine<long long>-;
+#pragma link C++ class ROOT::Experimental::RHistEngine<float>-;
+#pragma link C++ class ROOT::Experimental::RHistEngine<double>-;
+
 #pragma link C++ class ROOT::Experimental::RRegularAxis-;
 #pragma link C++ class ROOT::Experimental::RVariableBinAxis-;
 #pragma link C++ class ROOT::Experimental::Internal::RAxes-;

hist/histv7/inc/ROOT/RAxes.hxx

@@ -22,29 +22,29 @@ class TBuffer;
 
 namespace ROOT {
 namespace Experimental {
+
+/// Variant of all supported axis types.
+using RAxisVariant = std::variant<RRegularAxis, RVariableBinAxis>;
+
 namespace Internal {
 
 /**
 Bin configurations for all dimensions of a histogram.
 */
 class RAxes final {
-public:
-   using AxisVariant = std::variant<RRegularAxis, RVariableBinAxis>;
-
-private:
-   std::vector<AxisVariant> fAxes;
+   std::vector<RAxisVariant> fAxes;
 
 public:
    /// \param[in] axes the axis objects, must have size > 0
-   explicit RAxes(std::vector<AxisVariant> axes) : fAxes(std::move(axes))
+   explicit RAxes(std::vector<RAxisVariant> axes) : fAxes(std::move(axes))
    {
       if (fAxes.empty()) {
          throw std::invalid_argument("must have at least 1 axis object");
       }
    }
 
    std::size_t GetNDimensions() const { return fAxes.size(); }
-   const std::vector<AxisVariant> &Get() const { return fAxes; }
+   const std::vector<RAxisVariant> &Get() const { return fAxes; }
 
    friend bool operator==(const RAxes &lhs, const RAxes &rhs) { return lhs.fAxes == rhs.fAxes; }
 

hist/histv7/inc/ROOT/RHistEngine.hxx

@@ -0,0 +1,212 @@
+/// \file
+/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+/// Feedback is welcome!
+
+#ifndef ROOT_RHistEngine
+#define ROOT_RHistEngine
+
+#include "RAxes.hxx"
+#include "RBinIndex.hxx"
+#include "RLinearizedIndex.hxx"
+#include "RRegularAxis.hxx"
+
+#include <array>
+#include <cassert>
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+class TBuffer;
+
+namespace ROOT {
+namespace Experimental {
+
+/**
+A histogram data structure to bin data along multiple dimensions.
+
+Every call to \ref Fill(const A &... args) "Fill" bins the data according to the axis configuration and increments the
+bin content:
+\code
+ROOT::Experimental::RHistEngine<int> hist(10, 5, 15);
+hist.Fill(8.5);
+// hist.GetBinContent(ROOT::Experimental::RBinIndex(3)) will return 1
+\endcode
+
+The class is templated on the bin content type. For counting, as in the example above, it may be an integer type such as
+`int` or `long`. Narrower types such as `unsigned char` or `short` are supported, but may overflow due to their limited
+range and must be used with care. For weighted filling, the bin content type must be a floating-point type such as
+`float` or `double`. Note that `float` has a limited significant precision of 24 bits.
+
+An object can have arbitrary dimensionality determined at run-time. The axis configuration is passed as a vector of
+RAxisVariant:
+\code
+std::vector<ROOT::Experimental::RAxisVariant> axes;
+axes.push_back(ROOT::Experimental::RRegularAxis(10, 5, 15));
+axes.push_back(ROOT::Experimental::RVariableBinAxis({1, 10, 100, 1000}));
+ROOT::Experimental::RHistEngine<int> hist(axes);
+// hist.GetNDimensions() will return 2
+\endcode
+
+\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+Feedback is welcome!
+*/
+template <typename BinContentType>
+class RHistEngine final {
+   /// The axis configuration for this histogram. Relevant methods are forwarded from the public interface.
+   Internal::RAxes fAxes;
+   /// The bin contents for this histogram
+   std::vector<BinContentType> fBinContents;
+
+public:
+   /// Construct a histogram engine.
+   ///
+   /// \param[in] axes the axis objects, must have size > 0
+   explicit RHistEngine(std::vector<RAxisVariant> axes) : fAxes(std::move(axes))
+   {
+      fBinContents.resize(fAxes.ComputeTotalNBins());
+   }
+
+   /// Construct a one-dimensional histogram engine with a regular axis.
+   ///
+   /// \param[in] nNormalBins the number of normal bins, must be > 0
+   /// \param[in] low the lower end of the axis interval (inclusive)
+   /// \param[in] high the upper end of the axis interval (exclusive), must be > low
+   /// \par See also
+   /// the \ref RRegularAxis::RRegularAxis(std::size_t nNormalBins, double low, double high, bool enableFlowBins)
+   /// "constructor of RRegularAxis"
+   RHistEngine(std::size_t nNormalBins, double low, double high) : RHistEngine({RRegularAxis(nNormalBins, low, high)})
+   {
+   }
+
+   // Copy constructor and assignment operator are deleted to avoid surprises.
+   RHistEngine(const RHistEngine<BinContentType> &) = delete;
+   RHistEngine(RHistEngine<BinContentType> &&) = default;
+   RHistEngine<BinContentType> &operator=(const RHistEngine<BinContentType> &) = delete;
+   RHistEngine<BinContentType> &operator=(RHistEngine<BinContentType> &&) = default;
+   ~RHistEngine() = default;
+
+   const std::vector<RAxisVariant> &GetAxes() const { return fAxes.Get(); }
+   std::size_t GetNDimensions() const { return fAxes.GetNDimensions(); }
+   std::size_t GetTotalNBins() const { return fBinContents.size(); }
+
+   /// Get the content of a single bin.
+   ///
+   /// \code
+   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
+   /// std::array<ROOT::Experimental::RBinIndex, 2> indices = {3, 5};
+   /// int content = hist.GetBinContent(indices);
+   /// \endcode
+   ///
+   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
+   /// values. See also the class documentation of RBinIndex.
+   ///
+   /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found.
+   ///
+   /// \param[in] indices the array of indices for each axis
+   /// \return the bin content
+   /// \par See also
+   /// the \ref GetBinContent(const A &... args) const "variadic function template overload" accepting arguments
+   /// directly
+   template <std::size_t N>
+   const BinContentType &GetBinContent(const std::array<RBinIndex, N> &indices) const
+   {
+      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
+      // be confusing for users.
+      if (N != GetNDimensions()) {
+         throw std::invalid_argument("invalid number of indices passed to GetBinContent");
+      }
+      RLinearizedIndex index = fAxes.ComputeGlobalIndex(indices);
+      if (!index.fValid) {
+         throw std::invalid_argument("bin not found in GetBinContent");
+      }
+      assert(index.fIndex < fBinContents.size());
+      return fBinContents[index.fIndex];
+   }
+
+   /// Get the content of a single bin.
+   ///
+   /// \code
+   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
+   /// int content = hist.GetBinContent(ROOT::Experimental::RBinIndex(3), ROOT::Experimental::RBinIndex(5));
+   /// // ... or construct the RBinIndex arguments implicitly from integers:
+   /// content = hist.GetBinContent(3, 5);
+   /// \endcode
+   ///
+   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
+   /// values. See also the class documentation of RBinIndex.
+   ///
+   /// Throws an exception if the number of arguments does not match the axis configuration or the bin is not found.
+   ///
+   /// \param[in] args the arguments for each axis
+   /// \return the bin content
+   /// \par See also
+   /// the \ref GetBinContent(const std::array<RBinIndex, N> &indices) const "function overload" accepting
+   /// `std::array`
+   template <typename... A>
+   const BinContentType &GetBinContent(const A &...args) const
+   {
+      std::array<RBinIndex, sizeof...(A)> indices{args...};
+      return GetBinContent(indices);
+   }
+
+   /// Fill an entry into the histogram.
+   ///
+   /// \code
+   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
+   /// auto args = std::make_tuple(8.5, 10.5);
+   /// hist.Fill(args);
+   /// \endcode
+   ///
+   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
+   /// discarded.
+   ///
+   /// Throws an exception if the number of arguments does not match the axis configuration.
+   ///
+   /// \param[in] args the arguments for each axis
+   /// \par See also
+   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly
+   template <typename... A>
+   void Fill(const std::tuple<A...> &args)
+   {
+      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
+      // be confusing for users.
+      if (sizeof...(A) != GetNDimensions()) {
+         throw std::invalid_argument("invalid number of arguments to Fill");
+      }
+      RLinearizedIndex index = fAxes.ComputeGlobalIndex(args);
+      if (index.fValid) {
+         assert(index.fIndex < fBinContents.size());
+         fBinContents[index.fIndex]++;
+      }
+   }
+
+   /// Fill an entry into the histogram.
+   ///
+   /// \code
+   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
+   /// hist.Fill(8.5, 10.5);
+   /// \endcode
+   ///
+   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
+   /// discarded.
+   ///
+   /// Throws an exception if the number of arguments does not match the axis configuration.
+   ///
+   /// \param[in] args the arguments for each axis
+   /// \par See also
+   /// the \ref Fill(const std::tuple<A...> &args) "function overload" accepting `std::tuple`
+   template <typename... A>
+   void Fill(const A &...args)
+   {
+      Fill(std::forward_as_tuple(args...));
+   }
+
+   /// %ROOT Streamer function to throw when trying to store an object of this class.
+   void Streamer(TBuffer &) { throw std::runtime_error("unable to store RHistEngine"); }
+};
+
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/test/CMakeLists.txt

@@ -1,4 +1,5 @@
 HIST_ADD_GTEST(hist_axes hist_axes.cxx)
+HIST_ADD_GTEST(hist_engine hist_engine.cxx)
 HIST_ADD_GTEST(hist_index hist_index.cxx)
 HIST_ADD_GTEST(hist_regular hist_regular.cxx)
 HIST_ADD_GTEST(hist_variable hist_variable.cxx)