[hist] Allow disabling of dimensions in RHistStats

hist/histv7/inc/ROOT/RHist.hxx

@@ -6,6 +6,7 @@
 #define ROOT_RHist
 
 #include "RBinIndex.hxx"
+#include "RCategoricalAxis.hxx"
 #include "RHistEngine.hxx"
 #include "RHistStats.hxx"
 #include "RWeight.hxx"
@@ -17,6 +18,7 @@
 #include <stdexcept>
 #include <tuple>
 #include <utility>
+#include <variant>
 #include <vector>
 
 class TBuffer;
@@ -73,7 +75,16 @@ public:
    /// Construct a histogram.
    ///
    /// \param[in] axes the axis objects, must have size > 0
-   explicit RHist(std::vector<RAxisVariant> axes) : fEngine(std::move(axes)), fStats(fEngine.GetNDimensions()) {}
+   explicit RHist(std::vector<RAxisVariant> axes) : fEngine(std::move(axes)), fStats(fEngine.GetNDimensions())
+   {
+      // The axes parameter was moved, use from the engine.
+      const auto &engineAxes = fEngine.GetAxes();
+      for (std::size_t i = 0; i < engineAxes.size(); i++) {
+         if (std::get_if<RCategoricalAxis>(&engineAxes[i])) {
+            fStats.DisableDimension(i);
+         }
+      }
+   }
 
    /// Construct a one-dimensional histogram engine with a regular axis.
    ///

hist/histv7/inc/ROOT/RHistFillContext.hxx

@@ -36,7 +36,16 @@ private:
    RHistStats fStats;
 
    /// \sa RHistConcurrentFiller::CreateFillContent()
-   explicit RHistFillContext(RHist<BinContentType> &hist) : fHist(&hist), fStats(hist.GetNDimensions()) {}
+   explicit RHistFillContext(RHist<BinContentType> &hist) : fHist(&hist), fStats(hist.GetNDimensions())
+   {
+      // Propagate disabled dimensions to the local histogram statistics object.
+      const auto &histStats = hist.GetStats();
+      for (std::size_t i = 0; i < histStats.GetNDimensions(); i++) {
+         if (!histStats.IsEnabled(i)) {
+            fStats.DisableDimension(i);
+         }
+      }
+   }
    RHistFillContext(const RHistFillContext<BinContentType> &) = delete;
    RHistFillContext(RHistFillContext<BinContentType> &&) = default;
    RHistFillContext<BinContentType> &operator=(const RHistFillContext<BinContentType> &) = delete;

hist/histv7/inc/ROOT/RHistStats.hxx

@@ -9,11 +9,13 @@
 #include "RLinearizedIndex.hxx"
 #include "RWeight.hxx"
 
+#include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <cstdint>
 #include <stdexcept>
 #include <tuple>
+#include <type_traits>
 #include <vector>
 
 class TBuffer;
@@ -40,13 +42,15 @@ class RHistStats final {
 public:
    /// Statistics for one dimension.
    struct RDimensionStats final {
+      bool fEnabled = true;
       double fSumWX = 0.0;
       double fSumWX2 = 0.0;
       double fSumWX3 = 0.0;
       double fSumWX4 = 0.0;
 
       void Add(double x)
       {
+         assert(fEnabled);
          fSumWX += x;
          fSumWX2 += x * x;
          fSumWX3 += x * x * x;
@@ -55,6 +59,7 @@ public:
 
       void Add(double x, double w)
       {
+         assert(fEnabled);
          fSumWX += w * x;
          fSumWX2 += w * x * x;
          fSumWX3 += w * x * x * x;
@@ -63,6 +68,7 @@ public:
 
       void Add(const RDimensionStats &other)
       {
+         assert(fEnabled);
          fSumWX += other.fSumWX;
          fSumWX2 += other.fSumWX2;
          fSumWX3 += other.fSumWX3;
@@ -74,6 +80,7 @@ public:
       /// \param[in] other another statistics object that must not be modified during the operation
       void AddAtomic(const RDimensionStats &other)
       {
+         assert(fEnabled);
          Internal::AtomicAdd(&fSumWX, other.fSumWX);
          Internal::AtomicAdd(&fSumWX2, other.fSumWX2);
          Internal::AtomicAdd(&fSumWX3, other.fSumWX3);
@@ -90,6 +97,7 @@ public:
 
       void Scale(double factor)
       {
+         assert(fEnabled);
          fSumWX *= factor;
          fSumWX2 *= factor;
          fSumWX3 *= factor;
@@ -125,7 +133,29 @@ public:
    double GetSumW() const { return fSumW; }
    double GetSumW2() const { return fSumW2; }
 
-   const RDimensionStats &GetDimensionStats(std::size_t dim = 0) const { return fDimensionStats.at(dim); }
+   /// Get the statistics object for one dimension.
+   ///
+   /// Throws an exception if the dimension is disabled.
+   ///
+   /// \param[in] dim the dimension index, starting at 0
+   /// \return the statistics object
+   const RDimensionStats &GetDimensionStats(std::size_t dim = 0) const
+   {
+      const RDimensionStats &stats = fDimensionStats.at(dim);
+      if (!stats.fEnabled) {
+         throw std::runtime_error("dimension is disabled");
+      }
+      return stats;
+   }
+
+   /// Disable one dimension of this statistics object.
+   ///
+   /// All future calls to Fill will ignore the corresponding argument. Once disabled, a dimension cannot be reenabled.
+   ///
+   /// \param[in] dim the dimension index, starting at 0
+   void DisableDimension(std::size_t dim) { fDimensionStats.at(dim).fEnabled = false; }
+
+   bool IsEnabled(std::size_t dim) const { return fDimensionStats.at(dim).fEnabled; }
 
    /// Add all entries from another statistics object.
    ///
@@ -141,7 +171,12 @@ public:
       fSumW += other.fSumW;
       fSumW2 += other.fSumW2;
       for (std::size_t i = 0; i < fDimensionStats.size(); i++) {
-         fDimensionStats[i].Add(other.fDimensionStats[i]);
+         if (fDimensionStats[i].fEnabled != other.fDimensionStats[i].fEnabled) {
+            throw std::invalid_argument("dimension enablement not identical in Add");
+         }
+         if (fDimensionStats[i].fEnabled) {
+            fDimensionStats[i].Add(other.fDimensionStats[i]);
+         }
       }
    }
 
@@ -159,7 +194,12 @@ public:
       Internal::AtomicAdd(&fSumW, other.fSumW);
       Internal::AtomicAdd(&fSumW2, other.fSumW2);
       for (std::size_t i = 0; i < fDimensionStats.size(); i++) {
-         fDimensionStats[i].AddAtomic(other.fDimensionStats[i]);
+         if (fDimensionStats[i].fEnabled != other.fDimensionStats[i].fEnabled) {
+            throw std::invalid_argument("dimension enablement not identical in Add");
+         }
+         if (fDimensionStats[i].fEnabled) {
+            fDimensionStats[i].AddAtomic(other.fDimensionStats[i]);
+         }
       }
    }
 
@@ -326,7 +366,14 @@ private:
    template <std::size_t I, typename... A>
    void FillImpl(const std::tuple<A...> &args)
    {
-      fDimensionStats[I].Add(std::get<I>(args));
+      using ArgumentType = std::tuple_element_t<I, std::tuple<A...>>;
+      if (fDimensionStats[I].fEnabled) {
+         if constexpr (std::is_convertible_v<ArgumentType, double>) {
+            fDimensionStats[I].Add(std::get<I>(args));
+         } else {
+            throw std::invalid_argument("invalid type of argument");
+         }
+      }
       if constexpr (I + 1 < sizeof...(A)) {
          FillImpl<I + 1>(args);
       }
@@ -335,7 +382,16 @@ private:
    template <std::size_t I, std::size_t N, typename... A>
    void FillImpl(const std::tuple<A...> &args, double w)
    {
-      fDimensionStats[I].Add(std::get<I>(args), w);
+      using ArgumentType = std::tuple_element_t<I, std::tuple<A...>>;
+      if (fDimensionStats[I].fEnabled) {
+         if constexpr (std::is_convertible_v<ArgumentType, double>) {
+            fDimensionStats[I].Add(std::get<I>(args), w);
+         } else {
+            throw std::invalid_argument("invalid type of argument");
+         }
+      }
+      // Avoid compiler warning about unused parameter...
+      (void)w;
       if constexpr (I + 1 < N) {
          FillImpl<I + 1, N>(args, w);
       }
@@ -439,7 +495,9 @@ public:
       fSumW *= factor;
       fSumW2 *= factor * factor;
       for (std::size_t i = 0; i < fDimensionStats.size(); i++) {
-         fDimensionStats[i].Scale(factor);
+         if (fDimensionStats[i].fEnabled) {
+            fDimensionStats[i].Scale(factor);
+         }
       }
    }
 

hist/histv7/test/hist_concurrent.cxx

@@ -2,7 +2,9 @@
 
 #include <memory>
 #include <stdexcept>
+#include <string>
 #include <utility>
+#include <vector>
 
 TEST(RHistConcurrentFiller, Constructor)
 {
@@ -141,6 +143,53 @@ TEST(RHistFillContext, StressFillWeight)
    EXPECT_FLOAT_EQ(hist->ComputeMean(), 0.5);
 }
 
+TEST(RHistFillContext, FillCategorical)
+{
+   const std::vector<std::string> categories = {"a", "b", "c"};
+   const RCategoricalAxis axis(categories);
+   const std::vector<RAxisVariant> axes = {axis};
+   auto hist = std::make_shared<RHist<int>>(axes);
+
+   {
+      RHistConcurrentFiller filler(hist);
+      auto context = filler.CreateFillContext();
+      context->Fill("b");
+      context->Fill(std::make_tuple("c"));
+   }
+
+   EXPECT_EQ(hist->GetBinContent(RBinIndex(1)), 1);
+   std::array<RBinIndex, 1> indices = {2};
+   EXPECT_EQ(hist->GetBinContent(indices), 1);
+
+   EXPECT_EQ(hist->GetNEntries(), 2);
+   EXPECT_FLOAT_EQ(hist->ComputeNEffectiveEntries(), 2);
+}
+
+TEST(RHistFillContext, FillCategoricalWeight)
+{
+   const std::vector<std::string> categories = {"a", "b", "c"};
+   const RCategoricalAxis axis(categories);
+   const std::vector<RAxisVariant> axes = {axis};
+   auto hist = std::make_shared<RHist<float>>(axes);
+
+   {
+      RHistConcurrentFiller filler(hist);
+      auto context = filler.CreateFillContext();
+      context->Fill("b", RWeight(0.8));
+      context->Fill(std::make_tuple("c"), RWeight(0.9));
+   }
+
+   EXPECT_FLOAT_EQ(hist->GetBinContent(RBinIndex(1)), 0.8);
+   std::array<RBinIndex, 1> indices = {2};
+   EXPECT_FLOAT_EQ(hist->GetBinContent(indices), 0.9);
+
+   EXPECT_EQ(hist->GetNEntries(), 2);
+   EXPECT_FLOAT_EQ(hist->GetStats().GetSumW(), 1.7);
+   EXPECT_FLOAT_EQ(hist->GetStats().GetSumW2(), 1.45);
+   // Cross-checked with TH1
+   EXPECT_FLOAT_EQ(hist->ComputeNEffectiveEntries(), 1.9931034);
+}
+
 TEST(RHistFillContext, Flush)
 {
    static constexpr std::size_t Bins = 20;

hist/histv7/test/hist_hist.cxx

@@ -172,6 +172,43 @@ TEST(RHist, FillWeight)
    EXPECT_FLOAT_EQ(hist.ComputeStdDev(), 0.49913420);
 }
 
+TEST(RHist, FillCategorical)
+{
+   const std::vector<std::string> categories = {"a", "b", "c"};
+   const RCategoricalAxis axis(categories);
+   RHist<int> hist({axis});
+
+   hist.Fill("b");
+   hist.Fill(std::make_tuple("c"));
+
+   EXPECT_EQ(hist.GetBinContent(RBinIndex(1)), 1);
+   std::array<RBinIndex, 1> indices = {2};
+   EXPECT_EQ(hist.GetBinContent(indices), 1);
+
+   EXPECT_EQ(hist.GetNEntries(), 2);
+   EXPECT_FLOAT_EQ(hist.ComputeNEffectiveEntries(), 2);
+}
+
+TEST(RHist, FillCategoricalWeight)
+{
+   const std::vector<std::string> categories = {"a", "b", "c"};
+   const RCategoricalAxis axis(categories);
+   RHist<float> hist({axis});
+
+   hist.Fill("b", RWeight(0.8));
+   hist.Fill(std::make_tuple("c"), RWeight(0.9));
+
+   EXPECT_FLOAT_EQ(hist.GetBinContent(RBinIndex(1)), 0.8);
+   std::array<RBinIndex, 1> indices = {2};
+   EXPECT_FLOAT_EQ(hist.GetBinContent(indices), 0.9);
+
+   EXPECT_EQ(hist.GetNEntries(), 2);
+   EXPECT_FLOAT_EQ(hist.GetStats().GetSumW(), 1.7);
+   EXPECT_FLOAT_EQ(hist.GetStats().GetSumW2(), 1.45);
+   // Cross-checked with TH1
+   EXPECT_FLOAT_EQ(hist.ComputeNEffectiveEntries(), 1.9931034);
+}
+
 TEST(RHist, Scale)
 {
    static constexpr std::size_t Bins = 20;