[hist] Implement RSliceSpec and Internal::RSliceBinIndexMapper

hist/histv7/doc/Terminology.md

@@ -47,6 +47,7 @@ During filling, it accepts an additional `double` value and accumulates its sum
 
 *Slicing* means to extract a subset of the normal bins in each dimension.
 Bin contents of excluded normal bins are added to the flow bins.
+During slicing, operations allow to *rebin* a dimension or *sum* the sliced bins.
 
 A *snapshot* is a consistent clone of the histogram during concurrent filling.
 

hist/histv7/headers.cmake

@@ -14,6 +14,8 @@ set(histv7_headers
     ROOT/RHistUtils.hxx
     ROOT/RLinearizedIndex.hxx
     ROOT/RRegularAxis.hxx
+    ROOT/RSliceBinIndexMapper.hxx
+    ROOT/RSliceSpec.hxx
     ROOT/RVariableBinAxis.hxx
     ROOT/RWeight.hxx
 )

hist/histv7/inc/ROOT/RBinIndexRange.hxx

@@ -58,6 +58,12 @@ public:
    RBinIndex GetEnd() const { return fEnd; }
    // fNNormalBins is not exposed because it might be confusing for partial ranges.
 
+   bool IsInvalid() const
+   {
+      // fEnd can legally be invalid for full ranges including the overflow bin.
+      return fBegin.IsInvalid();
+   }
+
    friend bool operator==(const RBinIndexRange &lhs, const RBinIndexRange &rhs)
    {
       return lhs.fBegin == rhs.fBegin && lhs.fEnd == rhs.fEnd && lhs.fNNormalBins == rhs.fNNormalBins;

hist/histv7/inc/ROOT/RSliceBinIndexMapper.hxx

@@ -0,0 +1,131 @@
+/// \file
+/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+/// Feedback is welcome!
+
+#ifndef ROOT_RSliceBinIndexMapper
+#define ROOT_RSliceBinIndexMapper
+
+#include "RBinIndex.hxx"
+#include "RBinIndexRange.hxx"
+#include "RSliceSpec.hxx"
+
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+namespace ROOT {
+namespace Experimental {
+namespace Internal {
+
+/**
+Mapper of bin indices for slice operations.
+*/
+class RSliceBinIndexMapper final {
+   /// The requested slice specifications
+   std::vector<RSliceSpec> fSliceSpecs;
+   /// The expected dimensionality of the mapped indices
+   std::size_t fMappedDimensionality;
+
+   static std::size_t ComputeMappedDimensionality(const std::vector<RSliceSpec> &sliceSpecs)
+   {
+      std::size_t dimensionality = 0;
+      for (auto &&spec : sliceSpecs) {
+         // A sum operation makes the dimension disappear.
+         if (spec.GetOperationSum() == nullptr) {
+            dimensionality++;
+         }
+      }
+      return dimensionality;
+   }
+
+public:
+   /// \param[in] sliceSpecs the slice specifications, must have size > 0
+   explicit RSliceBinIndexMapper(std::vector<RSliceSpec> sliceSpecs)
+      : fSliceSpecs(std::move(sliceSpecs)), fMappedDimensionality(ComputeMappedDimensionality(fSliceSpecs))
+   {
+      if (fSliceSpecs.empty()) {
+         throw std::invalid_argument("must have at least 1 slice specification");
+      }
+   }
+
+   const std::vector<RSliceSpec> &GetSliceSpecs() const { return fSliceSpecs; }
+   std::size_t GetMappedDimensionality() const { return fMappedDimensionality; }
+
+   /// Map a vector of RBinIndex according to the slice specifications.
+   ///
+   /// \param[in] original the original bin indices
+   /// \param[out] mapped the mapped bin indices
+   /// \return whether the mapping was successful or the bin content should be discarded
+   bool Map(const std::vector<RBinIndex> &original, std::vector<RBinIndex> &mapped) const
+   {
+      if (original.size() != fSliceSpecs.size()) {
+         throw std::invalid_argument("invalid number of original indices passed to RSliceBinIndexMapper::Map");
+      }
+      if (mapped.size() != fMappedDimensionality) {
+         throw std::invalid_argument("invalid size of mapped indices passed to RSliceBinIndexMapper::Map");
+      }
+
+      std::size_t mappedPos = 0;
+      for (std::size_t i = 0; i < original.size(); i++) {
+         RBinIndex index = original[i];
+         if (index.IsInvalid()) {
+            throw std::invalid_argument("invalid bin index passed to RSliceBinIndexMapper::Map");
+         }
+
+         const RSliceSpec &sliceSpec = fSliceSpecs[i];
+         const auto &range = sliceSpec.GetRange();
+         bool contained = true;
+         if (!range.IsInvalid()) {
+            // Underflow and overflow indices map to themselves, but they may not actually be contained in the range.
+            // This is important for the sum operation below.
+            if (index.IsUnderflow()) {
+               contained = range.GetBegin().IsUnderflow();
+            } else if (index.IsOverflow()) {
+               contained = range.GetEnd().IsInvalid();
+            } else if (index.IsNormal()) {
+               const auto &begin = range.GetBegin();
+               const auto &end = range.GetEnd();
+               if (begin.IsNormal() && index < begin) {
+                  index = RBinIndex::Underflow();
+                  contained = false;
+               } else if (end.IsNormal() && index >= end) {
+                  index = RBinIndex::Overflow();
+                  contained = false;
+               } else if (begin.IsNormal()) {
+                  // This normal bin is contained in the range. Its index must be shifted according to the begin of the
+                  // range.
+                  index -= begin.GetIndex();
+                  assert(!index.IsInvalid());
+               }
+            }
+         }
+
+         if (auto *opRebin = sliceSpec.GetOperationRebin()) {
+            if (index.IsNormal()) {
+               index = RBinIndex(index.GetIndex() / opRebin->GetNGroup());
+            }
+         } else if (sliceSpec.GetOperationSum() != nullptr) {
+            // This dimension disappears. If there is a range and the index is not contained, discard it.
+            if (!contained) {
+               return false;
+            }
+            // Otherwise got to the next dimension.
+            continue;
+         }
+
+         mapped[mappedPos] = index;
+         mappedPos++;
+      }
+
+      // If we got here, the loop should have filled all mapped indices.
+      assert(mappedPos == mapped.size());
+      return true;
+   }
+};
+
+} // namespace Internal
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/inc/ROOT/RSliceSpec.hxx

@@ -0,0 +1,99 @@
+/// \file
+/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+/// Feedback is welcome!
+
+#ifndef ROOT_RSliceSpec
+#define ROOT_RSliceSpec
+
+#include "RBinIndexRange.hxx"
+
+#include <cstdint>
+#include <stdexcept>
+#include <utility>
+#include <variant>
+
+namespace ROOT {
+namespace Experimental {
+
+/**
+Specification of a slice operation along one dimension.
+
+\code
+using ROOT::Experimental::RSliceSpec;
+// When not specifying a range, the slice will include all bins.
+RSliceSpec full;
+// In the following, assuming range is an RBinIndexRange.
+RSliceSpec slice(range);
+
+// Operations are specified with parameters.
+RSliceSpec rebin(RSliceSpec::ROperationRebin(2));
+RSliceSpec sum(RSliceSpec::ROperationSum{});
+
+// Finally, it is possible to combine a range and an operation.
+RSliceSpec sliceRebin(range, RSliceSpec::ROperationRebin(2));
+RSliceSpec sliceSum(range, RSliceSpec::ROperationSum{});
+\endcode
+
+\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+Feedback is welcome!
+*/
+class RSliceSpec final {
+public:
+   /// Rebin the dimension, grouping a number of original bins into a new one.
+   class ROperationRebin final {
+      std::uint64_t fNGroup = 1;
+
+   public:
+      /// \param[in] nGroup the number of bins to group, must be > 0
+      ROperationRebin(std::uint64_t nGroup) : fNGroup(nGroup)
+      {
+         if (nGroup == 0) {
+            throw std::invalid_argument("nGroup must be > 0");
+         }
+      }
+
+      std::uint64_t GetNGroup() const { return fNGroup; }
+   };
+
+   /// Sum bins along this dimension, effectively resulting in a projection.
+   class ROperationSum final {
+      // empty, no parameters at the moment
+   };
+
+private:
+   /// The range of the slice; can be invalid to signify the full range
+   RBinIndexRange fRange;
+   /// The operation to perform, if any
+   std::variant<std::monostate, ROperationRebin, ROperationSum> fOperation;
+
+public:
+   /// A default slice operation that keeps the dimension untouched.
+   RSliceSpec() = default;
+
+   /// A slice of a dimension.
+   ///
+   /// \param[in] range the range of the slice
+   RSliceSpec(RBinIndexRange range) : fRange(std::move(range)) {}
+
+   /// A rebin operation of a dimension.
+   RSliceSpec(ROperationRebin rebin) : fOperation(std::move(rebin)) {}
+
+   /// A sum operation of a dimension.
+   RSliceSpec(ROperationSum sum) : fOperation(std::move(sum)) {}
+
+   /// A rebin operation of a slice of the dimension.
+   RSliceSpec(RBinIndexRange range, ROperationRebin rebin) : fRange(std::move(range)), fOperation(std::move(rebin)) {}
+
+   /// A sum operation of a slice of the dimension.
+   RSliceSpec(RBinIndexRange range, ROperationSum sum) : fRange(std::move(range)), fOperation(std::move(sum)) {}
+
+   const RBinIndexRange &GetRange() const { return fRange; }
+   bool HasOperation() const { return !std::holds_alternative<std::monostate>(fOperation); }
+   const ROperationRebin *GetOperationRebin() const { return std::get_if<ROperationRebin>(&fOperation); }
+   const ROperationSum *GetOperationSum() const { return std::get_if<ROperationSum>(&fOperation); }
+};
+
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/test/CMakeLists.txt

@@ -9,6 +9,7 @@ HIST_ADD_GTEST(hist_engine_atomic hist_engine_atomic.cxx)
 HIST_ADD_GTEST(hist_hist hist_hist.cxx)
 HIST_ADD_GTEST(hist_index hist_index.cxx)
 HIST_ADD_GTEST(hist_regular hist_regular.cxx)
+HIST_ADD_GTEST(hist_slice hist_slice.cxx)
 HIST_ADD_GTEST(hist_stats hist_stats.cxx)
 HIST_ADD_GTEST(hist_user hist_user.cxx)
 HIST_ADD_GTEST(hist_variable hist_variable.cxx)

hist/histv7/test/hist_categorical.cxx

@@ -137,13 +137,15 @@ TEST(RCategoricalAxis, GetNormalRange)
       const auto normal = axis.GetNormalRange();
       EXPECT_EQ(normal.GetBegin(), index0);
       EXPECT_EQ(normal.GetEnd(), indexBins);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), categories.size());
    }
 
    {
       const auto normal = axis.GetNormalRange(index0, indexBins);
       EXPECT_EQ(normal.GetBegin(), index0);
       EXPECT_EQ(normal.GetEnd(), indexBins);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), categories.size());
    }
 
@@ -152,13 +154,15 @@ TEST(RCategoricalAxis, GetNormalRange)
       const auto normal = axis.GetNormalRange(index1, index2);
       EXPECT_EQ(normal.GetBegin(), index1);
       EXPECT_EQ(normal.GetEnd(), index2);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), 1);
    }
 
    {
       const auto empty = axis.GetNormalRange(index1, index1);
       EXPECT_EQ(empty.GetBegin(), index1);
       EXPECT_EQ(empty.GetEnd(), index1);
+      EXPECT_FALSE(empty.IsInvalid());
       EXPECT_EQ(empty.begin(), empty.end());
       EXPECT_EQ(std::distance(empty.begin(), empty.end()), 0);
    }
@@ -181,6 +185,7 @@ TEST(RCategoricalAxis, GetFullRange)
       const auto full = axis.GetFullRange();
       EXPECT_EQ(full.GetBegin(), RBinIndex(0));
       EXPECT_EQ(full.GetEnd(), RBinIndex());
+      EXPECT_FALSE(full.IsInvalid());
       EXPECT_EQ(std::distance(full.begin(), full.end()), categories.size() + 1);
    }
 
@@ -189,6 +194,7 @@ TEST(RCategoricalAxis, GetFullRange)
       const auto full = axisNoOverflow.GetFullRange();
       EXPECT_EQ(full.GetBegin(), RBinIndex(0));
       EXPECT_EQ(full.GetEnd(), RBinIndex(categories.size()));
+      EXPECT_FALSE(full.IsInvalid());
       EXPECT_EQ(std::distance(full.begin(), full.end()), categories.size());
    }
 }

hist/histv7/test/hist_index.cxx

@@ -177,15 +177,18 @@ TEST(RBinIndexRange, ConstructorCreate)
    const RBinIndexRange invalid;
    EXPECT_TRUE(invalid.GetBegin().IsInvalid());
    EXPECT_TRUE(invalid.GetEnd().IsInvalid());
+   EXPECT_TRUE(invalid.IsInvalid());
 
    const auto index0 = RBinIndex(0);
    const auto range0 = CreateBinIndexRange(index0, index0, 0);
    EXPECT_EQ(range0.GetBegin(), index0);
    EXPECT_EQ(range0.GetEnd(), index0);
+   EXPECT_FALSE(range0.IsInvalid());
 
    const auto range01 = CreateBinIndexRange(index0, RBinIndex(1), 1);
    EXPECT_EQ(range01.GetBegin(), index0);
    EXPECT_EQ(range01.GetEnd(), RBinIndex(1));
+   EXPECT_FALSE(range01.IsInvalid());
 }
 
 TEST(RBinIndexRange, Equality)

hist/histv7/test/hist_regular.cxx

@@ -197,13 +197,15 @@ TEST(RRegularAxis, GetNormalRange)
       const auto normal = axis.GetNormalRange();
       EXPECT_EQ(normal.GetBegin(), index0);
       EXPECT_EQ(normal.GetEnd(), indexBins);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), Bins);
    }
 
    {
       const auto normal = axis.GetNormalRange(index0, indexBins);
       EXPECT_EQ(normal.GetBegin(), index0);
       EXPECT_EQ(normal.GetEnd(), indexBins);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), Bins);
    }
 
@@ -212,13 +214,15 @@ TEST(RRegularAxis, GetNormalRange)
       const auto normal = axis.GetNormalRange(index1, index5);
       EXPECT_EQ(normal.GetBegin(), index1);
       EXPECT_EQ(normal.GetEnd(), index5);
+      EXPECT_FALSE(normal.IsInvalid());
       EXPECT_EQ(std::distance(normal.begin(), normal.end()), 4);
    }
 
    {
       const auto empty = axis.GetNormalRange(index1, index1);
       EXPECT_EQ(empty.GetBegin(), index1);
       EXPECT_EQ(empty.GetEnd(), index1);
+      EXPECT_FALSE(empty.IsInvalid());
       EXPECT_EQ(empty.begin(), empty.end());
       EXPECT_EQ(std::distance(empty.begin(), empty.end()), 0);
    }
@@ -241,6 +245,7 @@ TEST(RRegularAxis, GetFullRange)
       const auto full = axis.GetFullRange();
       EXPECT_EQ(full.GetBegin(), RBinIndex::Underflow());
       EXPECT_EQ(full.GetEnd(), RBinIndex());
+      EXPECT_FALSE(full.IsInvalid());
       EXPECT_EQ(std::distance(full.begin(), full.end()), Bins + 2);
    }
 
@@ -249,6 +254,7 @@ TEST(RRegularAxis, GetFullRange)
       const auto full = axisNoFlowBins.GetFullRange();
       EXPECT_EQ(full.GetBegin(), RBinIndex(0));
       EXPECT_EQ(full.GetEnd(), RBinIndex(Bins));
+      EXPECT_FALSE(full.IsInvalid());
       EXPECT_EQ(std::distance(full.begin(), full.end()), Bins);
    }
 }