diff --git a/hist/histv7/benchmark/CMakeLists.txt b/hist/histv7/benchmark/CMakeLists.txt index 40da903c90aa0..31d2b07487bcf 100644 --- a/hist/histv7/benchmark/CMakeLists.txt +++ b/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) diff --git a/hist/histv7/benchmark/hist_benchmark_axes.cxx b/hist/histv7/benchmark/hist_benchmark_axes.cxx new file mode 100644 index 0000000000000..50bede22702cb --- /dev/null +++ b/hist/histv7/benchmark/hist_benchmark_axes.cxx @@ -0,0 +1,69 @@ +#include + +#include + +#include +#include +#include + +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 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 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(); diff --git a/hist/histv7/headers.cmake b/hist/histv7/headers.cmake index 7b9d2163ea46a..c394f0d26d47d 100644 --- a/hist/histv7/headers.cmake +++ b/hist/histv7/headers.cmake @@ -1,4 +1,5 @@ set(histv7_headers + ROOT/RAxes.hxx ROOT/RLinearizedIndex.hxx ROOT/RRegularAxis.hxx ROOT/RVariableBinAxis.hxx diff --git a/hist/histv7/inc/LinkDef.h b/hist/histv7/inc/LinkDef.h index 9456f4a501ffc..69adac7caf80e 100644 --- a/hist/histv7/inc/LinkDef.h +++ b/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-; diff --git a/hist/histv7/inc/ROOT/RAxes.hxx b/hist/histv7/inc/ROOT/RAxes.hxx new file mode 100644 index 0000000000000..d44d01f353667 --- /dev/null +++ b/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 +#include +#include +#include +#include + +class TBuffer; + +namespace ROOT { +namespace Experimental { +namespace Internal { + +/** +Bin configurations for all dimensions of a histogram. +*/ +class RAxes final { +public: + using AxisVariant = std::variant; + +private: + std::vector fAxes; + +public: + /// \param[in] axes the axis objects, must have size > 0 + explicit RAxes(std::vector 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 &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(&axis)) { + totalNumBins *= regular->GetTotalNumBins(); + } else if (auto *variable = std::get_if(&axis)) { + totalNumBins *= variable->GetTotalNumBins(); + } else { + throw std::logic_error("unimplemented axis type"); + } + } + return totalNumBins; + } + +private: + template + RLinearizedIndex ComputeGlobalIndex(std::size_t index, const std::tuple &args) const + { + const auto &axis = fAxes[I]; + RLinearizedIndex linIndex; + if (auto *regular = std::get_if(&axis)) { + index *= regular->GetTotalNumBins(); + linIndex = regular->ComputeLinearizedIndex(std::get(args)); + } else if (auto *variable = std::get_if(&axis)) { + index *= variable->GetTotalNumBins(); + linIndex = variable->ComputeLinearizedIndex(std::get(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(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 + RLinearizedIndex ComputeGlobalIndex(const std::tuple &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 diff --git a/hist/histv7/inc/ROOT/RLinearizedIndex.hxx b/hist/histv7/inc/ROOT/RLinearizedIndex.hxx index e47377f8406a8..d373c24416943 100644 --- a/hist/histv7/inc/ROOT/RLinearizedIndex.hxx +++ b/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 diff --git a/hist/histv7/inc/ROOT/RVariableBinAxis.hxx b/hist/histv7/inc/ROOT/RVariableBinAxis.hxx index 9a4eabc981a95..20788ba008968 100644 --- a/hist/histv7/inc/ROOT/RVariableBinAxis.hxx +++ b/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 binEdges, bool enableFlowBins = true) + explicit RVariableBinAxis(std::vector binEdges, bool enableFlowBins = true) : fBinEdges(std::move(binEdges)), fEnableFlowBins(enableFlowBins) { if (fBinEdges.size() < 2) { diff --git a/hist/histv7/test/CMakeLists.txt b/hist/histv7/test/CMakeLists.txt index 852356fb3f9a7..e450ec9a650ba 100644 --- a/hist/histv7/test/CMakeLists.txt +++ b/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) diff --git a/hist/histv7/test/hist_axes.cxx b/hist/histv7/test/hist_axes.cxx new file mode 100644 index 0000000000000..006bb8bc5d35d --- /dev/null +++ b/hist/histv7/test/hist_axes.cxx @@ -0,0 +1,169 @@ +#include "hist_test.hxx" + +#include +#include +#include +#include + +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 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(&v[0]) != nullptr); + EXPECT_TRUE(std::get_if(&v[1]) != nullptr); + + std::vector 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 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 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 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 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); +} diff --git a/hist/histv7/test/hist_io.cxx b/hist/histv7/test/hist_io.cxx index c228ae9ca2300..63faea579f41a 100644 --- a/hist/histv7/test/hist_io.cxx +++ b/hist/histv7/test/hist_io.cxx @@ -34,3 +34,19 @@ TEST(RVariableBinAxis, Streamer) const RVariableBinAxis axis(bins); ExpectThrowOnWriteObject(axis); } + +TEST(RAxes, Streamer) +{ + static constexpr std::size_t BinsX = 20; + const RRegularAxis regularAxis(BinsX, 0, BinsX); + static constexpr std::size_t BinsY = 30; + std::vector 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}); + ExpectThrowOnWriteObject(axes); +} diff --git a/hist/histv7/test/hist_test.hxx b/hist/histv7/test/hist_test.hxx index c36ba8a722432..7ed6a31788d00 100644 --- a/hist/histv7/test/hist_test.hxx +++ b/hist/histv7/test/hist_test.hxx @@ -1,6 +1,7 @@ #ifndef hist_test #define hist_test +#include #include #include @@ -8,5 +9,6 @@ using ROOT::Experimental::RRegularAxis; using ROOT::Experimental::RVariableBinAxis; +using ROOT::Experimental::Internal::RAxes; #endif