[hist] Implement initial RAxes

hist/histv7/benchmark/CMakeLists.txt

@@ -1,2 +1,5 @@
+add_executable(hist_benchmark_axes hist_benchmark_axes.cxx)
+target_link_libraries(hist_benchmark_axes 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_axes.cxx

@@ -0,0 +1,69 @@
+#include <ROOT/RAxes.hxx>
+
+#include <benchmark/benchmark.h>
+
+#include <random>
+#include <tuple>
+#include <vector>
+
+struct RAxes_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::Internal::RAxes axes{{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(RAxes_Regular1, ComputeGlobalIndex)(benchmark::State &state)
+{
+   for (auto _ : state) {
+      for (double number : fNumbers) {
+         benchmark::DoNotOptimize(axes.ComputeGlobalIndex(std::make_tuple(number)));
+      }
+   }
+}
+BENCHMARK_REGISTER_F(RAxes_Regular1, ComputeGlobalIndex)->Range(0, 32768);
+
+struct RAxes_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::Internal::RAxes axes{{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(RAxes_Regular2, ComputeGlobalIndex)(benchmark::State &state)
+{
+   for (auto _ : state) {
+      for (std::size_t i = 0; i < fNumbers.size(); i += 2) {
+         benchmark::DoNotOptimize(axes.ComputeGlobalIndex(std::make_tuple(fNumbers[i], fNumbers[i + 1])));
+      }
+   }
+}
+BENCHMARK_REGISTER_F(RAxes_Regular2, ComputeGlobalIndex)->Range(0, 32768);
+
+BENCHMARK_MAIN();

hist/histv7/headers.cmake

@@ -1,4 +1,5 @@
 set(histv7_headers
+    ROOT/RAxes.hxx
     ROOT/RLinearizedIndex.hxx
     ROOT/RRegularAxis.hxx
     ROOT/RVariableBinAxis.hxx

hist/histv7/inc/LinkDef.h

@@ -1,2 +1,3 @@
 #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

@@ -0,0 +1,115 @@
+/// \file
+/// \warning This is part of the ROOT 7 prototype! It will change without notice. It might trigger earthquakes. Feedback
+/// is welcome!
+
+#ifndef ROOT_RAxes
+#define ROOT_RAxes
+
+#include "RLinearizedIndex.hxx"
+#include "RRegularAxis.hxx"
+#include "RVariableBinAxis.hxx"
+
+#include <stdexcept>
+#include <tuple>
+#include <utility>
+#include <variant>
+#include <vector>
+
+class TBuffer;
+
+namespace ROOT {
+namespace Experimental {
+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;
+
+public:
+   /// \param[in] axes the axis objects, must have size > 0
+   explicit RAxes(std::vector<AxisVariant> axes) : fAxes(std::move(axes))
+   {
+      if (fAxes.empty()) {
+         throw std::invalid_argument("must have at least 1 axis object");
+      }
+   }
+
+   std::size_t GetNumDimensions() const { return fAxes.size(); }
+   const std::vector<AxisVariant> &Get() const { return fAxes; }
+
+   friend bool operator==(const RAxes &lhs, const RAxes &rhs) { return lhs.fAxes == rhs.fAxes; }
+
+   /// Compute the total number of bins for all axes.
+   ///
+   /// It is the product of each dimension's total number of bins.
+   ///
+   /// \return the total number of bins
+   std::size_t ComputeTotalNumBins() const
+   {
+      std::size_t totalNumBins = 1;
+      for (auto &&axis : fAxes) {
+         if (auto *regular = std::get_if<RRegularAxis>(&axis)) {
+            totalNumBins *= regular->GetTotalNumBins();
+         } else if (auto *variable = std::get_if<RVariableBinAxis>(&axis)) {
+            totalNumBins *= variable->GetTotalNumBins();
+         } else {
+            throw std::logic_error("unimplemented axis type");
+         }
+      }
+      return totalNumBins;
+   }
+
+private:
+   template <std::size_t I, typename... A>
+   RLinearizedIndex ComputeGlobalIndex(std::size_t index, const std::tuple<A...> &args) const
+   {
+      const auto &axis = fAxes[I];
+      RLinearizedIndex linIndex;
+      if (auto *regular = std::get_if<RRegularAxis>(&axis)) {
+         index *= regular->GetTotalNumBins();
+         linIndex = regular->ComputeLinearizedIndex(std::get<I>(args));
+      } else if (auto *variable = std::get_if<RVariableBinAxis>(&axis)) {
+         index *= variable->GetTotalNumBins();
+         linIndex = variable->ComputeLinearizedIndex(std::get<I>(args));
+      } else {
+         throw std::logic_error("unimplemented axis type");
+      }
+      if (!linIndex.fValid) {
+         return {0, false};
+      }
+      index += linIndex.fIndex;
+      if constexpr (I + 1 < sizeof...(A)) {
+         return ComputeGlobalIndex<I + 1>(index, args);
+      }
+      return {index, true};
+   }
+
+public:
+   /// Compute the global index for all axes.
+   ///
+   /// \param[in] args the arguments
+   /// \return the global index that may be invalid
+   template <typename... A>
+   RLinearizedIndex ComputeGlobalIndex(const std::tuple<A...> &args) const
+   {
+      if (sizeof...(A) != fAxes.size()) {
+         throw std::invalid_argument("invalid number of arguments to ComputeGlobalIndex");
+      }
+      return ComputeGlobalIndex<0, A...>(0, 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 RAxes"); }
+};
+
+} // namespace Internal
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/inc/ROOT/RLinearizedIndex.hxx

@@ -19,8 +19,8 @@ For example, when an argument is outside the axis and underflow / overflow bins
 welcome!
 */
 struct RLinearizedIndex final {
-   std::size_t fIndex;
-   bool fValid;
+   std::size_t fIndex = 0;
+   bool fValid = false;
 };
 
 } // namespace Experimental

hist/histv7/inc/ROOT/RVariableBinAxis.hxx

@@ -44,7 +44,7 @@ public:
    ///
    /// \param[in] binEdges the (ordered) edges of the normal bins, must define at least one bin (i.e. size >= 2)
    /// \param[in] enableFlowBins whether to enable underflow and overflow bins
-   RVariableBinAxis(std::vector<double> binEdges, bool enableFlowBins = true)
+   explicit RVariableBinAxis(std::vector<double> binEdges, bool enableFlowBins = true)
       : fBinEdges(std::move(binEdges)), fEnableFlowBins(enableFlowBins)
    {
       if (fBinEdges.size() < 2) {

hist/histv7/test/CMakeLists.txt

@@ -1,3 +1,4 @@
+HIST_ADD_GTEST(hist_axes hist_axes.cxx)
 HIST_ADD_GTEST(hist_regular hist_regular.cxx)
 HIST_ADD_GTEST(hist_variable hist_variable.cxx)
 

hist/histv7/test/hist_axes.cxx

@@ -0,0 +1,169 @@
+#include "hist_test.hxx"
+
+#include <stdexcept>
+#include <tuple>
+#include <variant>
+#include <vector>
+
+TEST(RAxes, Constructor)
+{
+   static constexpr std::size_t BinsX = 20;
+   const RRegularAxis regularAxis(BinsX, 0, BinsX);
+   static constexpr std::size_t BinsY = 30;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   RAxes axes({regularAxis, variableBinAxis});
+   EXPECT_EQ(axes.GetNumDimensions(), 2);
+   const auto &v = axes.Get();
+   ASSERT_EQ(v.size(), 2);
+   EXPECT_EQ(v[0].index(), 0);
+   EXPECT_EQ(v[1].index(), 1);
+   EXPECT_TRUE(std::get_if<RRegularAxis>(&v[0]) != nullptr);
+   EXPECT_TRUE(std::get_if<RVariableBinAxis>(&v[1]) != nullptr);
+
+   std::vector<RAxes::AxisVariant> newAxes{variableBinAxis, regularAxis};
+   axes = RAxes(newAxes);
+   EXPECT_EQ(axes.GetNumDimensions(), 2);
+
+   EXPECT_THROW(RAxes({}), std::invalid_argument);
+}
+
+TEST(RAxes, Equality)
+{
+   static constexpr std::size_t BinsX = 20;
+   const RRegularAxis regularAxis(BinsX, 0, BinsX);
+   static constexpr std::size_t BinsY = 30;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   const RAxes axesA({regularAxis, variableBinAxis});
+   const RAxes axesA2({regularAxis, variableBinAxis});
+   const RAxes axesB({variableBinAxis, regularAxis});
+   const RAxes axesC({regularAxis});
+   const RAxes axesD({variableBinAxis});
+
+   EXPECT_TRUE(axesA == axesA);
+   EXPECT_TRUE(axesA == axesA2);
+   EXPECT_TRUE(axesA2 == axesA);
+
+   EXPECT_FALSE(axesA == axesB);
+   EXPECT_FALSE(axesA == axesC);
+   EXPECT_FALSE(axesA == axesD);
+
+   EXPECT_FALSE(axesB == axesC);
+   EXPECT_FALSE(axesB == axesD);
+
+   EXPECT_FALSE(axesC == axesD);
+}
+
+TEST(RAxes, ComputeTotalNumBins)
+{
+   static constexpr std::size_t BinsX = 20;
+   const RRegularAxis regularAxis(BinsX, 0, BinsX);
+   static constexpr std::size_t BinsY = 30;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+   const RAxes axes({regularAxis, variableBinAxis});
+
+   // Both axes include underflow and overflow bins.
+   EXPECT_EQ(axes.ComputeTotalNumBins(), (BinsX + 2) * (BinsY + 2));
+}
+
+TEST(RAxes, ComputeGlobalIndex)
+{
+   static constexpr std::size_t BinsX = 20;
+   const RRegularAxis regularAxis(BinsX, 0, BinsX);
+   static constexpr std::size_t BinsY = 30;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+   const RAxes axes({regularAxis, variableBinAxis});
+
+   {
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(1.5, 2.5));
+      EXPECT_EQ(globalIndex.fIndex, 1 * (BinsY + 2) + 2);
+      EXPECT_TRUE(globalIndex.fValid);
+   }
+
+   {
+      // Underflow bin of the first axis.
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(-1, 2.5));
+      EXPECT_EQ(globalIndex.fIndex, BinsX * (BinsY + 2) + 2);
+      EXPECT_TRUE(globalIndex.fValid);
+   }
+
+   {
+      // Overflow bin of the second axis.
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(1.5, 42));
+      EXPECT_EQ(globalIndex.fIndex, 1 * (BinsY + 2) + BinsY + 1);
+      EXPECT_TRUE(globalIndex.fValid);
+   }
+}
+
+TEST(RAxes, ComputeGlobalIndexNoFlowBins)
+{
+   static constexpr std::size_t BinsX = 20;
+   const RRegularAxis regularAxis(BinsX, 0, BinsX, /*enableFlowBins=*/false);
+   static constexpr std::size_t BinsY = 30;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins, /*enableFlowBins=*/false);
+   const RAxes axes({regularAxis, variableBinAxis});
+   ASSERT_EQ(axes.ComputeTotalNumBins(), BinsX * BinsY);
+
+   {
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(1.5, 2.5));
+      EXPECT_EQ(globalIndex.fIndex, 1 * BinsY + 2);
+      EXPECT_TRUE(globalIndex.fValid);
+   }
+
+   {
+      // Underflow bin of the first axis.
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(-1, 2.5));
+      EXPECT_EQ(globalIndex.fIndex, 0);
+      EXPECT_FALSE(globalIndex.fValid);
+   }
+
+   {
+      // Overflow bin of the second axis.
+      const auto globalIndex = axes.ComputeGlobalIndex(std::make_tuple(1.5, 42));
+      EXPECT_EQ(globalIndex.fIndex, 0);
+      EXPECT_FALSE(globalIndex.fValid);
+   }
+}
+
+TEST(RAxes, ComputeGlobalIndexInvalidNumberOfArguments)
+{
+   static constexpr std::size_t Bins = 20;
+   const RRegularAxis axis(Bins, 0, Bins);
+   const RAxes axes1({axis});
+   ASSERT_EQ(axes1.GetNumDimensions(), 1);
+   const RAxes axes2({axis, axis});
+   ASSERT_EQ(axes2.GetNumDimensions(), 2);
+
+   EXPECT_NO_THROW(axes1.ComputeGlobalIndex(std::make_tuple(1)));
+   EXPECT_THROW(axes1.ComputeGlobalIndex(std::make_tuple(1, 2)), std::invalid_argument);
+
+   EXPECT_THROW(axes2.ComputeGlobalIndex(std::make_tuple(1)), std::invalid_argument);
+   EXPECT_NO_THROW(axes2.ComputeGlobalIndex(std::make_tuple(1, 2)));
+   EXPECT_THROW(axes2.ComputeGlobalIndex(std::make_tuple(1, 2, 3)), std::invalid_argument);
+}