diff --git a/hist/histv7/inc/ROOT/RAxisVariant.hxx b/hist/histv7/inc/ROOT/RAxisVariant.hxx index e912bd35dbaed..9884003460c6e 100644 --- a/hist/histv7/inc/ROOT/RAxisVariant.hxx +++ b/hist/histv7/inc/ROOT/RAxisVariant.hxx @@ -127,6 +127,26 @@ public: } } + /// Slice this axis according to the specification. + /// + /// Axes throw exceptions if the slicing cannot be performed. For example, the rebin operation must divide the + /// number of normal bins for RRegularAxis and RVariableBinAxis, while RCategoricalAxis cannot be sliced at all. + /// + /// \param[in] sliceSpec the slice specification + /// \return the sliced axis + RAxisVariant Slice(const RSliceSpec &sliceSpec) const + { + if (auto *regular = GetRegularAxis()) { + return regular->Slice(sliceSpec); + } else if (auto *variable = GetVariableBinAxis()) { + return variable->Slice(sliceSpec); + } else if (auto *categorical = GetCategoricalAxis()) { + return categorical->Slice(sliceSpec); + } else { + throw std::logic_error("unimplemented axis type"); // GCOVR_EXCL_LINE + } + } + friend bool operator==(const RAxisVariant &lhs, const RAxisVariant &rhs) { return lhs.fVariant == rhs.fVariant; } /// %ROOT Streamer function to throw when trying to store an object of this class. diff --git a/hist/histv7/inc/ROOT/RCategoricalAxis.hxx b/hist/histv7/inc/ROOT/RCategoricalAxis.hxx index 60bb4fe69c4f7..5cb161562f6bb 100644 --- a/hist/histv7/inc/ROOT/RCategoricalAxis.hxx +++ b/hist/histv7/inc/ROOT/RCategoricalAxis.hxx @@ -8,6 +8,7 @@ #include "RBinIndex.hxx" #include "RBinIndexRange.hxx" #include "RLinearizedIndex.hxx" +#include "RSliceSpec.hxx" #include #include @@ -167,6 +168,31 @@ public: : GetNormalRange(); } + /// Slice this axis according to the specification. + /// + /// A categorical axis cannot be sliced. The method will throw if a specification other than the default slice + /// operation is passed. + /// + /// \param[in] sliceSpec the slice specification + /// \return the sliced / copied axis + RCategoricalAxis Slice(const RSliceSpec &sliceSpec) const + { + if (sliceSpec.GetOperationSum() != nullptr) { + throw std::runtime_error("sum operation makes dimension disappear"); + } + + if (!sliceSpec.GetRange().IsInvalid()) { + throw std::runtime_error("slicing of RCategoricalAxis not implemented"); + } + if (sliceSpec.GetOperationRebin() != nullptr) { + throw std::runtime_error("cannot rebin RCategoricalAxis"); + } + + // The sliced axis always has flow bins enabled, for symmetry with other axis types. + bool enableOverflowBin = true; + return RCategoricalAxis(fCategories, enableOverflowBin); + } + /// %ROOT Streamer function to throw when trying to store an object of this class. void Streamer(TBuffer &) { throw std::runtime_error("unable to store RCategoricalAxis"); } }; diff --git a/hist/histv7/inc/ROOT/RRegularAxis.hxx b/hist/histv7/inc/ROOT/RRegularAxis.hxx index 67c27346063f5..d6078ec4560b8 100644 --- a/hist/histv7/inc/ROOT/RRegularAxis.hxx +++ b/hist/histv7/inc/ROOT/RRegularAxis.hxx @@ -8,6 +8,7 @@ #include "RBinIndex.hxx" #include "RBinIndexRange.hxx" #include "RLinearizedIndex.hxx" +#include "RSliceSpec.hxx" #include #include @@ -83,6 +84,24 @@ public: double GetHigh() const { return fHigh; } bool HasFlowBins() const { return fEnableFlowBins; } + /// Compute the low edge of a bin. + /// + /// \param[in] bin the index, must be \f$< fNNormalBins\f$ + double ComputeLowEdge(std::uint64_t bin) const + { + assert(bin < fNNormalBins); + return fLow + (fHigh - fLow) * bin / fNNormalBins; + } + + /// Compute the high edge of a bin. + /// + /// \param[in] bin the index, must be \f$< fNNormalBins\f$ + double ComputeHighEdge(std::uint64_t bin) const + { + assert(bin < fNNormalBins); + return fLow + (fHigh - fLow) * (bin + 1) / fNNormalBins; + } + friend bool operator==(const RRegularAxis &lhs, const RRegularAxis &rhs) { return lhs.fNNormalBins == rhs.fNNormalBins && lhs.fLow == rhs.fLow && lhs.fHigh == rhs.fHigh && @@ -194,6 +213,59 @@ public: : GetNormalRange(); } + /// Slice this axis according to the specification. + /// + /// Throws an exception if the axis cannot be sliced: + /// * A sum operation makes the dimension disappear. + /// * The rebin operation must divide the number of normal bins. + /// + /// \param[in] sliceSpec the slice specification + /// \return the sliced axis, with enabled underflow and overflow bins + RRegularAxis Slice(const RSliceSpec &sliceSpec) const + { + if (sliceSpec.GetOperationSum() != nullptr) { + throw std::runtime_error("sum operation makes dimension disappear"); + } + + // Figure out the properties of the sliced axis. + std::uint64_t nNormalBins = fNNormalBins; + double low = fLow; + double high = fHigh; + + const auto &range = sliceSpec.GetRange(); + if (!range.IsInvalid()) { + std::uint64_t begin = 0; + std::uint64_t end = nNormalBins; + if (range.GetBegin().IsNormal()) { + begin = range.GetBegin().GetIndex(); + // Only compute a new lower end of the axis interval if needed. + if (begin > 0) { + low = ComputeLowEdge(begin); + } + } + if (range.GetEnd().IsNormal()) { + end = range.GetEnd().GetIndex(); + // Only compute a new upper end of the axis interval if needed, to avoid floating-point inaccuracies. + if (end < nNormalBins) { + high = ComputeHighEdge(end - 1); + } + } + nNormalBins = end - begin; + } + + if (auto *opRebin = sliceSpec.GetOperationRebin()) { + if (nNormalBins % opRebin->GetNGroup() != 0) { + throw std::runtime_error("rebin operation does not divide number of normal bins"); + } + nNormalBins /= opRebin->GetNGroup(); + } + + // The sliced axis always has flow bins enabled to preserve all entries. This is the least confusing for users, + // even if not always strictly necessary. + bool enableFlowBins = true; + return RRegularAxis(nNormalBins, {low, high}, enableFlowBins); + } + /// %ROOT Streamer function to throw when trying to store an object of this class. void Streamer(TBuffer &) { throw std::runtime_error("unable to store RRegularAxis"); } }; diff --git a/hist/histv7/inc/ROOT/RVariableBinAxis.hxx b/hist/histv7/inc/ROOT/RVariableBinAxis.hxx index 47e5f8970b2a2..58faa12ada28b 100644 --- a/hist/histv7/inc/ROOT/RVariableBinAxis.hxx +++ b/hist/histv7/inc/ROOT/RVariableBinAxis.hxx @@ -8,6 +8,7 @@ #include "RBinIndex.hxx" #include "RBinIndexRange.hxx" #include "RLinearizedIndex.hxx" +#include "RSliceSpec.hxx" #include #include @@ -196,6 +197,57 @@ public: : GetNormalRange(); } + /// Slice this axis according to the specification. + /// + /// Throws an exception if the axis cannot be sliced: + /// * A sum operation makes the dimension disappear. + /// * The rebin operation must divide the number of normal bins. + /// + /// \param[in] sliceSpec the slice specification + /// \return the sliced axis, with enabled underflow and overflow bins + RVariableBinAxis Slice(const RSliceSpec &sliceSpec) const + { + if (sliceSpec.GetOperationSum() != nullptr) { + throw std::runtime_error("sum operation makes dimension disappear"); + } + + // Figure out the properties of the sliced axis. + std::size_t nNormalBins = fBinEdges.size() - 1; + std::size_t begin = 0; + std::size_t end = nNormalBins; + + const auto &range = sliceSpec.GetRange(); + if (!range.IsInvalid()) { + if (range.GetBegin().IsNormal()) { + begin = range.GetBegin().GetIndex(); + } + if (range.GetEnd().IsNormal()) { + end = range.GetEnd().GetIndex(); + } + nNormalBins = end - begin; + } + + std::uint64_t nGroup = 1; + if (auto *opRebin = sliceSpec.GetOperationRebin()) { + nGroup = opRebin->GetNGroup(); + if (nNormalBins % nGroup != 0) { + throw std::runtime_error("rebin operation does not divide number of normal bins"); + } + } + + // Prepare the bin edges. + std::vector binEdges; + for (std::size_t bin = begin; bin < end; bin += nGroup) { + binEdges.push_back(fBinEdges[bin]); + } + binEdges.push_back(fBinEdges[end]); + + // The sliced axis always has flow bins enabled to preserve all entries. This is the least confusing for users, + // even if not always strictly necessary. + bool enableFlowBins = true; + return RVariableBinAxis(std::move(binEdges), enableFlowBins); + } + /// %ROOT Streamer function to throw when trying to store an object of this class. void Streamer(TBuffer &) { throw std::runtime_error("unable to store RVariableBinAxis"); } }; diff --git a/hist/histv7/test/hist_axis_variant.cxx b/hist/histv7/test/hist_axis_variant.cxx index 09c84923182eb..76f89accecae9 100644 --- a/hist/histv7/test/hist_axis_variant.cxx +++ b/hist/histv7/test/hist_axis_variant.cxx @@ -23,6 +23,11 @@ TEST(RAxisVariant, RegularAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), Bins + 2); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 0); + EXPECT_TRUE(slicedAxis.GetRegularAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), Bins); } { @@ -36,6 +41,11 @@ TEST(RAxisVariant, RegularAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), Bins); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 0); + EXPECT_TRUE(slicedAxis.GetRegularAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), Bins); } } @@ -64,6 +74,11 @@ TEST(RAxisVariant, VariableBinAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), Bins + 2); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 1); + EXPECT_TRUE(slicedAxis.GetVariableBinAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), Bins); } { @@ -77,6 +92,11 @@ TEST(RAxisVariant, VariableBinAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), Bins); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 1); + EXPECT_TRUE(slicedAxis.GetVariableBinAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), Bins); } } @@ -100,6 +120,11 @@ TEST(RAxisVariant, CategoricalAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), 4); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 2); + EXPECT_TRUE(slicedAxis.GetCategoricalAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), 3); } { @@ -113,5 +138,10 @@ TEST(RAxisVariant, CategoricalAxis) EXPECT_EQ(std::distance(normal12.begin(), normal12.end()), 1); const auto full = axis.GetFullRange(); EXPECT_EQ(std::distance(full.begin(), full.end()), 3); + + const auto slicedAxis = axis.Slice(RSliceSpec{}); + EXPECT_EQ(slicedAxis.GetVariant().index(), 2); + EXPECT_TRUE(slicedAxis.GetCategoricalAxis() != nullptr); + EXPECT_EQ(slicedAxis.GetNNormalBins(), 3); } } diff --git a/hist/histv7/test/hist_categorical.cxx b/hist/histv7/test/hist_categorical.cxx index ef9f5269932b9..a5ecd4d2646d2 100644 --- a/hist/histv7/test/hist_categorical.cxx +++ b/hist/histv7/test/hist_categorical.cxx @@ -198,3 +198,35 @@ TEST(RCategoricalAxis, GetFullRange) EXPECT_EQ(std::distance(full.begin(), full.end()), categories.size()); } } + +static void Test_RCategoricalAxis_Slice(bool enableOverflowBin) +{ + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis origAxis(categories, enableOverflowBin); + ASSERT_EQ(origAxis.HasOverflowBin(), enableOverflowBin); + + { + // The only allowed "slicing" is to keep the entire axis. + const RSliceSpec full; + const auto axis = origAxis.Slice(full); + EXPECT_EQ(axis.GetNNormalBins(), 3); + EXPECT_EQ(axis.GetCategories(), categories); + EXPECT_TRUE(axis.HasOverflowBin()); + } + + // Slicing and other operations are not allowed / implemented. + EXPECT_THROW(origAxis.Slice(origAxis.GetFullRange()), std::runtime_error); + EXPECT_THROW(origAxis.Slice(origAxis.GetNormalRange()), std::runtime_error); + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationRebin(2)), std::runtime_error); + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationSum{}), std::runtime_error); +} + +TEST(RCategoricalAxis, Slice) +{ + Test_RCategoricalAxis_Slice(true); +} + +TEST(RCategoricalAxis, SliceNoOverflowBin) +{ + Test_RCategoricalAxis_Slice(false); +} diff --git a/hist/histv7/test/hist_regular.cxx b/hist/histv7/test/hist_regular.cxx index 8b4793af912a9..f010f738688fd 100644 --- a/hist/histv7/test/hist_regular.cxx +++ b/hist/histv7/test/hist_regular.cxx @@ -32,6 +32,17 @@ TEST(RRegularAxis, Constructor) EXPECT_THROW(RRegularAxis(Bins, {0, NaN}), std::invalid_argument); } +TEST(RRegularAxis, ComputeLowHighEdge) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + + for (std::size_t i = 0; i < Bins; i++) { + EXPECT_EQ(axis.ComputeLowEdge(i), i); + EXPECT_EQ(axis.ComputeHighEdge(i), i + 1); + } +} + TEST(RRegularAxis, Equality) { static constexpr std::size_t Bins = 20; @@ -258,3 +269,62 @@ TEST(RRegularAxis, GetFullRange) EXPECT_EQ(std::distance(full.begin(), full.end()), Bins); } } + +static void Test_RegularAxis_Slice(bool enableFlowBins) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis origAxis(Bins, {0, Bins}, enableFlowBins); + ASSERT_EQ(origAxis.HasFlowBins(), enableFlowBins); + + // Three different ways of "slicing" which will keep the entire axis. + for (auto sliceSpec : {RSliceSpec{}, RSliceSpec(origAxis.GetFullRange()), RSliceSpec(origAxis.GetNormalRange())}) { + const auto axis = origAxis.Slice(sliceSpec); + EXPECT_EQ(axis.GetNNormalBins(), Bins); + EXPECT_EQ(axis.GetLow(), 0); + EXPECT_EQ(axis.GetHigh(), Bins); + EXPECT_TRUE(axis.HasFlowBins()); + } + + { + const RSliceSpec slice(origAxis.GetNormalRange(1, 5)); + const auto axis = origAxis.Slice(slice); + EXPECT_EQ(axis.GetNNormalBins(), 4); + EXPECT_EQ(axis.GetLow(), 1); + EXPECT_EQ(axis.GetHigh(), 5); + EXPECT_TRUE(axis.HasFlowBins()); + } + + { + const RSliceSpec rebin(RSliceSpec::ROperationRebin(2)); + const auto axis = origAxis.Slice(rebin); + EXPECT_EQ(axis.GetNNormalBins(), Bins / 2); + EXPECT_EQ(axis.GetLow(), 0); + EXPECT_EQ(axis.GetHigh(), Bins); + EXPECT_TRUE(axis.HasFlowBins()); + } + + // Rebin grouping must divide the number of normal bins. + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationRebin(3)), std::runtime_error); + + // Sum operation makes dimension disappear. + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationSum{}), std::runtime_error); + + { + const RSliceSpec sliceRebin(origAxis.GetNormalRange(1, 5), RSliceSpec::ROperationRebin(2)); + const auto axis = origAxis.Slice(sliceRebin); + EXPECT_EQ(axis.GetNNormalBins(), 2); + EXPECT_EQ(axis.GetLow(), 1); + EXPECT_EQ(axis.GetHigh(), 5); + EXPECT_TRUE(axis.HasFlowBins()); + } +} + +TEST(RRegularAxis, Slice) +{ + Test_RegularAxis_Slice(true); +} + +TEST(RRegularAxis, SliceNoFlowBins) +{ + Test_RegularAxis_Slice(false); +} diff --git a/hist/histv7/test/hist_variable.cxx b/hist/histv7/test/hist_variable.cxx index 284958d703100..fb0d74012a9b0 100644 --- a/hist/histv7/test/hist_variable.cxx +++ b/hist/histv7/test/hist_variable.cxx @@ -332,3 +332,69 @@ TEST(RVariableBinAxis, GetFullRange) EXPECT_EQ(std::distance(full.begin(), full.end()), Bins); } } + +static void Test_RVariableBinAxis_Slice(bool enableFlowBins) +{ + static constexpr std::size_t Bins = 20; + std::vector bins; + for (std::size_t i = 0; i < Bins; i++) { + bins.push_back(i); + } + bins.push_back(Bins); + const RVariableBinAxis origAxis(bins, enableFlowBins); + ASSERT_EQ(origAxis.HasFlowBins(), enableFlowBins); + + // Three different ways of "slicing" which will keep the entire axis. + for (auto sliceSpec : {RSliceSpec{}, RSliceSpec(origAxis.GetFullRange()), RSliceSpec(origAxis.GetNormalRange())}) { + const auto axis = origAxis.Slice(sliceSpec); + EXPECT_EQ(axis.GetNNormalBins(), Bins); + EXPECT_EQ(axis.GetBinEdges(), bins); + EXPECT_TRUE(axis.HasFlowBins()); + } + + { + const RSliceSpec slice(origAxis.GetNormalRange(1, 5)); + const auto axis = origAxis.Slice(slice); + EXPECT_EQ(axis.GetNNormalBins(), 4); + EXPECT_EQ(axis.GetBinEdges().front(), 1); + EXPECT_EQ(axis.GetBinEdges().back(), 5); + EXPECT_TRUE(axis.HasFlowBins()); + } + + { + const RSliceSpec rebin(RSliceSpec::ROperationRebin(2)); + const auto axis = origAxis.Slice(rebin); + EXPECT_EQ(axis.GetNNormalBins(), Bins / 2); + EXPECT_EQ(axis.GetBinEdges().front(), 0); + EXPECT_EQ(axis.GetBinEdges().back(), Bins); + EXPECT_TRUE(axis.HasFlowBins()); + } + + // Rebin grouping must divide the number of normal bins. + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationRebin(3)), std::runtime_error); + + // Sum operation makes dimension disappear. + EXPECT_THROW(origAxis.Slice(RSliceSpec::ROperationSum{}), std::runtime_error); + + { + const RSliceSpec sliceRebin(origAxis.GetNormalRange(1, 5), RSliceSpec::ROperationRebin(2)); + const auto axis = origAxis.Slice(sliceRebin); + EXPECT_EQ(axis.GetNNormalBins(), 2); + const auto &binEdges = axis.GetBinEdges(); + ASSERT_EQ(binEdges.size(), 3); + EXPECT_EQ(binEdges[0], 1); + EXPECT_EQ(binEdges[1], 3); + EXPECT_EQ(binEdges[2], 5); + EXPECT_TRUE(axis.HasFlowBins()); + } +} + +TEST(RVariableBinAxis, Slice) +{ + Test_RVariableBinAxis_Slice(true); +} + +TEST(RVariableBinAxis, SliceNoFlowBins) +{ + Test_RVariableBinAxis_Slice(false); +}