[hist] Implement weighted filling of RHistStats

Author: hahnjo

hist/histv7/inc/ROOT/RHistStats.hxx

@@ -5,7 +5,9 @@
 #ifndef ROOT_RHistStats
 #define ROOT_RHistStats
 
+#include "RHistUtils.hxx"
 #include "RLinearizedIndex.hxx"
+#include "RWeight.hxx"
 
 #include <cmath>
 #include <cstdint>
@@ -49,6 +51,14 @@ public:
          fSumWX3 += x * x * x;
          fSumWX4 += x * x * x * x;
       }
+
+      void Add(double x, double w)
+      {
+         fSumWX += w * x;
+         fSumWX2 += w * x * x;
+         fSumWX3 += w * x * x * x;
+         fSumWX4 += w * x * x * x * x;
+      }
    };
 
 private:
@@ -232,6 +242,15 @@ 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);
+      if constexpr (I + 1 < N) {
+         FillImpl<I + 1, N>(args, w);
+      }
+   }
+
 public:
    /// Fill an entry into this statistics object.
    ///
@@ -245,7 +264,8 @@ public:
    ///
    /// \param[in] args the arguments for each dimension
    /// \par See also
-   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly
+   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the
+   /// \ref Fill(const std::tuple<A...> &args, RWeight weight) "overload for weighted filling"
    template <typename... A>
    void Fill(const std::tuple<A...> &args)
    {
@@ -258,22 +278,67 @@ public:
       FillImpl<0>(args);
    }
 
+   /// Fill an entry into this statistics object with a weight.
+   ///
+   /// \code
+   /// ROOT::Experimental::RHistStats stats(2);
+   /// auto args = std::make_tuple(8.5, 10.5);
+   /// stats.Fill(args, ROOT::Experimental::RWeight(0.8));
+   /// \endcode
+   ///
+   /// \param[in] args the arguments for each dimension
+   /// \param[in] weight the weight for this entry
+   /// \par See also
+   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the
+   /// \ref Fill(const std::tuple<A...> &args) "overload for unweighted filling"
+   template <typename... A>
+   void Fill(const std::tuple<A...> &args, RWeight weight)
+   {
+      if (sizeof...(A) != fDimensionStats.size()) {
+         throw std::invalid_argument("invalid number of arguments to Fill");
+      }
+      fNEntries++;
+      double w = weight.fValue;
+      fSumW += w;
+      fSumW2 += w * w;
+      FillImpl<0, sizeof...(A)>(args, w);
+   }
+
    /// Fill an entry into this statistics object.
    ///
    /// \code
    /// ROOT::Experimental::RHistStats stats(2);
    /// stats.Fill(8.5, 10.5);
    /// \endcode
+   /// For weighted filling, pass an RWeight as the last argument:
+   /// \code
+   /// ROOT::Experimental::RHistStats stats(2);
+   /// stats.Fill(8.5, 10.5, ROOT::Experimental::RWeight(0.8));
+   /// \endcode
    ///
    /// Throws an exception if the number of arguments does not match the number of dimensions.
    ///
    /// \param[in] args the arguments for each dimension
    /// \par See also
-   /// the \ref Fill(const std::tuple<A...> &args) "function overload" accepting `std::tuple`
+   /// the function overloads accepting `std::tuple` \ref Fill(const std::tuple<A...> &args) "for unweighted filling"
+   /// and \ref Fill(const std::tuple<A...> &args, RWeight) "for weighted filling"
    template <typename... A>
    void Fill(const A &...args)
    {
-      Fill(std::forward_as_tuple(args...));
+      auto t = std::forward_as_tuple(args...);
+      if constexpr (std::is_same_v<typename Internal::LastType<A...>::type, RWeight>) {
+         static constexpr std::size_t N = sizeof...(A) - 1;
+         if (N != fDimensionStats.size()) {
+            throw std::invalid_argument("invalid number of arguments to Fill");
+         }
+         fNEntries++;
+         double w = std::get<N>(t).fValue;
+         fSumW += w;
+         fSumW2 += w * w;
+         FillImpl<0, N>(t, w);
+      } else {
+         Fill(t);
+      }
    }
 
    /// %ROOT Streamer function to throw when trying to store an object of this class.

hist/histv7/test/hist_stats.cxx

@@ -48,6 +48,40 @@ TEST(RHistStats, GetDimensionStats)
    }
 }
 
+TEST(RHistStats, GetDimensionStatsWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   {
+      const auto &dimensionStats = stats.GetDimensionStats(/*=0*/);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX, 93.1);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX2, 1330.0);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX3, 20489.98);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX4, 330265.6);
+   }
+   {
+      const auto &dimensionStats = stats.GetDimensionStats(1);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX, 2 * 93.1);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX2, 4 * 1330.0);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX3, 8 * 20489.98);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX4, 16 * 330265.6);
+   }
+   {
+      const auto &dimensionStats = stats.GetDimensionStats(2);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX, 1330.0);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX2, 330265.6);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX3, 93164182.0);
+      EXPECT_DOUBLE_EQ(dimensionStats.fSumWX4, 28108731464.8);
+   }
+}
+
 TEST(RHistStats, ComputeNEffectiveEntries)
 {
    RHistStats stats(1);
@@ -65,6 +99,24 @@ TEST(RHistStats, ComputeNEffectiveEntries)
    EXPECT_DOUBLE_EQ(stats.ComputeNEffectiveEntries(), Entries);
 }
 
+TEST(RHistStats, ComputeNEffectiveEntriesWeighted)
+{
+   RHistStats stats(1);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeNEffectiveEntries(), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(1, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   EXPECT_DOUBLE_EQ(stats.GetSumW(), 7.7);
+   EXPECT_DOUBLE_EQ(stats.GetSumW2(), 3.563);
+   // Cross-checked with TH1
+   EXPECT_DOUBLE_EQ(stats.ComputeNEffectiveEntries(), 16.640471512770137);
+}
+
 TEST(RHistStats, ComputeMean)
 {
    RHistStats stats(3);
@@ -84,6 +136,26 @@ TEST(RHistStats, ComputeMean)
    EXPECT_DOUBLE_EQ(stats.ComputeMean(2), 123.5);
 }
 
+TEST(RHistStats, ComputeMeanWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeMean(/*=0*/), 0);
+   EXPECT_EQ(stats.ComputeMean(1), 0);
+   EXPECT_EQ(stats.ComputeMean(2), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   // Cross-checked with TH1
+   EXPECT_DOUBLE_EQ(stats.ComputeMean(/*=0*/), 12.090909090909090);
+   EXPECT_DOUBLE_EQ(stats.ComputeMean(1), 24.181818181818180);
+   EXPECT_DOUBLE_EQ(stats.ComputeMean(2), 172.72727272727272);
+}
+
 TEST(RHistStats, ComputeVariance)
 {
    RHistStats stats(3);
@@ -103,6 +175,26 @@ TEST(RHistStats, ComputeVariance)
    EXPECT_DOUBLE_EQ(stats.ComputeVariance(2), 12881.05);
 }
 
+TEST(RHistStats, ComputeVarianceWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeVariance(/*=0*/), 0);
+   EXPECT_EQ(stats.ComputeVariance(1), 0);
+   EXPECT_EQ(stats.ComputeVariance(2), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   // Cross-checked with TH1::GetStdDev squared, numerical differences with EXPECT_DOUBLE_EQ
+   EXPECT_FLOAT_EQ(stats.ComputeVariance(/*=0*/), 26.5371900);
+   EXPECT_FLOAT_EQ(stats.ComputeVariance(1), 106.148760);
+   EXPECT_FLOAT_EQ(stats.ComputeVariance(2), 13056.9256);
+}
+
 TEST(RHistStats, ComputeStdDev)
 {
    RHistStats stats(3);
@@ -122,6 +214,26 @@ TEST(RHistStats, ComputeStdDev)
    EXPECT_DOUBLE_EQ(stats.ComputeStdDev(2), std::sqrt(12881.05));
 }
 
+TEST(RHistStats, ComputeStdDevWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeStdDev(/*=0*/), 0);
+   EXPECT_EQ(stats.ComputeStdDev(1), 0);
+   EXPECT_EQ(stats.ComputeStdDev(2), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   // Cross-checked with TH1, numerical differences with EXPECT_DOUBLE_EQ
+   EXPECT_FLOAT_EQ(stats.ComputeStdDev(/*=0*/), 5.15142602);
+   EXPECT_FLOAT_EQ(stats.ComputeStdDev(1), 10.3028520);
+   EXPECT_FLOAT_EQ(stats.ComputeStdDev(2), 114.266905);
+}
+
 TEST(RHistStats, ComputeSkewness)
 {
    RHistStats stats(3);
@@ -142,6 +254,26 @@ TEST(RHistStats, ComputeSkewness)
    EXPECT_FLOAT_EQ(stats.ComputeSkewness(2), 0.66125456);
 }
 
+TEST(RHistStats, ComputeSkewnessWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeSkewness(/*=0*/), 0);
+   EXPECT_EQ(stats.ComputeSkewness(1), 0);
+   EXPECT_EQ(stats.ComputeSkewness(2), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   // Cross-checked with TH1, numerical differences with EXPECT_DOUBLE_EQ
+   EXPECT_FLOAT_EQ(stats.ComputeSkewness(/*=0*/), -0.50554999);
+   EXPECT_FLOAT_EQ(stats.ComputeSkewness(1), -0.50554999);
+   EXPECT_FLOAT_EQ(stats.ComputeSkewness(2), 0.12072240);
+}
+
 TEST(RHistStats, ComputeKurtosis)
 {
    RHistStats stats(3);
@@ -162,6 +294,26 @@ TEST(RHistStats, ComputeKurtosis)
    EXPECT_FLOAT_EQ(stats.ComputeKurtosis(2), -0.84198253);
 }
 
+TEST(RHistStats, ComputeKurtosisWeighted)
+{
+   RHistStats stats(3);
+   ASSERT_EQ(stats.GetNEntries(), 0);
+   EXPECT_EQ(stats.ComputeKurtosis(/*=0*/), 0);
+   EXPECT_EQ(stats.ComputeKurtosis(1), 0);
+   EXPECT_EQ(stats.ComputeKurtosis(2), 0);
+
+   static constexpr std::size_t Entries = 20;
+   for (std::size_t i = 0; i < Entries; i++) {
+      stats.Fill(i, 2 * i, i * i, RWeight(0.1 + 0.03 * i));
+   }
+
+   ASSERT_EQ(stats.GetNEntries(), Entries);
+   // Cross-checked with TH1, numerical differences with EXPECT_DOUBLE_EQ
+   EXPECT_FLOAT_EQ(stats.ComputeKurtosis(/*=0*/), -0.74828797);
+   EXPECT_FLOAT_EQ(stats.ComputeKurtosis(1), -0.74828797);
+   EXPECT_FLOAT_EQ(stats.ComputeKurtosis(2), -1.2483086);
+}
+
 TEST(RHistStats, FillInvalidNumberOfArguments)
 {
    RHistStats stats1(1);
@@ -174,3 +326,29 @@ TEST(RHistStats, FillInvalidNumberOfArguments)
    EXPECT_NO_THROW(stats2.Fill(1, 2));
    EXPECT_THROW(stats2.Fill(1, 2, 3), std::invalid_argument);
 }
+
+TEST(RHistStats, FillWeightInvalidNumberOfArguments)
+{
+   RHistStats stats1(1);
+   RHistStats stats2(2);
+
+   EXPECT_NO_THROW(stats1.Fill(1, RWeight(1)));
+   EXPECT_THROW(stats1.Fill(1, 2, RWeight(1)), std::invalid_argument);
+
+   EXPECT_THROW(stats2.Fill(1, RWeight(1)), std::invalid_argument);
+   EXPECT_NO_THROW(stats2.Fill(1, 2, RWeight(1)));
+   EXPECT_THROW(stats2.Fill(1, 2, 3, RWeight(1)), std::invalid_argument);
+}
+
+TEST(RHistStats, FillTupleWeightInvalidNumberOfArguments)
+{
+   RHistStats stats1(1);
+   RHistStats stats2(2);
+
+   EXPECT_NO_THROW(stats1.Fill(std::make_tuple(1), RWeight(1)));
+   EXPECT_THROW(stats1.Fill(std::make_tuple(1, 2), RWeight(1)), std::invalid_argument);
+
+   EXPECT_THROW(stats2.Fill(std::make_tuple(1), RWeight(1)), std::invalid_argument);
+   EXPECT_NO_THROW(stats2.Fill(std::make_tuple(1, 2), RWeight(1)));
+   EXPECT_THROW(stats2.Fill(std::make_tuple(1, 2, 3), RWeight(1)), std::invalid_argument);
+}