diff --git a/hist/histv7/inc/ROOT/RHist.hxx b/hist/histv7/inc/ROOT/RHist.hxx index 84de3c8a33e38..06b42a22db25c 100644 --- a/hist/histv7/inc/ROOT/RHist.hxx +++ b/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 #include #include +#include #include class TBuffer; @@ -73,7 +75,16 @@ public: /// Construct a histogram. /// /// \param[in] axes the axis objects, must have size > 0 - explicit RHist(std::vector axes) : fEngine(std::move(axes)), fStats(fEngine.GetNDimensions()) {} + explicit RHist(std::vector 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(&engineAxes[i])) { + fStats.DisableDimension(i); + } + } + } /// Construct a one-dimensional histogram engine with a regular axis. /// diff --git a/hist/histv7/inc/ROOT/RHistFillContext.hxx b/hist/histv7/inc/ROOT/RHistFillContext.hxx index 44e184d34854c..f4c0372ef3cfa 100644 --- a/hist/histv7/inc/ROOT/RHistFillContext.hxx +++ b/hist/histv7/inc/ROOT/RHistFillContext.hxx @@ -36,7 +36,16 @@ private: RHistStats fStats; /// \sa RHistConcurrentFiller::CreateFillContent() - explicit RHistFillContext(RHist &hist) : fHist(&hist), fStats(hist.GetNDimensions()) {} + explicit RHistFillContext(RHist &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 &) = delete; RHistFillContext(RHistFillContext &&) = default; RHistFillContext &operator=(const RHistFillContext &) = delete; diff --git a/hist/histv7/inc/ROOT/RHistStats.hxx b/hist/histv7/inc/ROOT/RHistStats.hxx index fcd302c0327ad..247a1372306f2 100644 --- a/hist/histv7/inc/ROOT/RHistStats.hxx +++ b/hist/histv7/inc/ROOT/RHistStats.hxx @@ -9,11 +9,13 @@ #include "RLinearizedIndex.hxx" #include "RWeight.hxx" +#include #include #include #include #include #include +#include #include class TBuffer; @@ -40,6 +42,7 @@ 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; @@ -47,6 +50,7 @@ public: 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 void FillImpl(const std::tuple &args) { - fDimensionStats[I].Add(std::get(args)); + using ArgumentType = std::tuple_element_t>; + if (fDimensionStats[I].fEnabled) { + if constexpr (std::is_convertible_v) { + fDimensionStats[I].Add(std::get(args)); + } else { + throw std::invalid_argument("invalid type of argument"); + } + } if constexpr (I + 1 < sizeof...(A)) { FillImpl(args); } @@ -335,7 +382,16 @@ private: template void FillImpl(const std::tuple &args, double w) { - fDimensionStats[I].Add(std::get(args), w); + using ArgumentType = std::tuple_element_t>; + if (fDimensionStats[I].fEnabled) { + if constexpr (std::is_convertible_v) { + fDimensionStats[I].Add(std::get(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(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); + } } } diff --git a/hist/histv7/test/hist_concurrent.cxx b/hist/histv7/test/hist_concurrent.cxx index 71e0995ddec68..73691cbe5762b 100644 --- a/hist/histv7/test/hist_concurrent.cxx +++ b/hist/histv7/test/hist_concurrent.cxx @@ -2,7 +2,9 @@ #include #include +#include #include +#include TEST(RHistConcurrentFiller, Constructor) { @@ -141,6 +143,53 @@ TEST(RHistFillContext, StressFillWeight) EXPECT_FLOAT_EQ(hist->ComputeMean(), 0.5); } +TEST(RHistFillContext, FillCategorical) +{ + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + const std::vector axes = {axis}; + auto hist = std::make_shared>(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 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 categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + const std::vector axes = {axis}; + auto hist = std::make_shared>(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 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; diff --git a/hist/histv7/test/hist_hist.cxx b/hist/histv7/test/hist_hist.cxx index 8bcbc12d48b24..a543c1653d5dc 100644 --- a/hist/histv7/test/hist_hist.cxx +++ b/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 categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + RHist hist({axis}); + + hist.Fill("b"); + hist.Fill(std::make_tuple("c")); + + EXPECT_EQ(hist.GetBinContent(RBinIndex(1)), 1); + std::array 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 categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + RHist 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 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; diff --git a/hist/histv7/test/hist_stats.cxx b/hist/histv7/test/hist_stats.cxx index 662d4d023bd60..eb9327b835ac5 100644 --- a/hist/histv7/test/hist_stats.cxx +++ b/hist/histv7/test/hist_stats.cxx @@ -82,15 +82,35 @@ TEST(RHistStats, GetDimensionStatsWeighted) } } +TEST(RHistStats, DisableDimension) +{ + RHistStats stats(3); + stats.DisableDimension(1); + EXPECT_TRUE(stats.IsEnabled(0)); + EXPECT_FALSE(stats.IsEnabled(1)); + EXPECT_TRUE(stats.IsEnabled(2)); + EXPECT_THROW(stats.GetDimensionStats(1), std::runtime_error); + + // The argument for the disabled dimension will be ignored. + stats.Fill(1, "b", 3); + stats.Fill(1, "b", 3, RWeight(1)); + stats.Fill(4, 5, 6); + stats.Fill(4, 5, 6, RWeight(1)); + + EXPECT_EQ(stats.GetNEntries(), 4); +} + TEST(RHistStats, Add) { - RHistStats statsA(2); - RHistStats statsB(2); + RHistStats statsA(3); + RHistStats statsB(3); + statsA.DisableDimension(1); + statsB.DisableDimension(1); static constexpr std::size_t Entries = 20; for (std::size_t i = 0; i < Entries; i++) { - statsA.Fill(i, 2 * i); - statsB.Fill(2 * i, 3 * i, RWeight(0.1 + 0.03 * i)); + statsA.Fill(i, 1, 2 * i); + statsB.Fill(2 * i, 1, 3 * i, RWeight(0.1 + 0.03 * i)); } statsA.Add(statsB); @@ -107,7 +127,7 @@ TEST(RHistStats, Add) EXPECT_FLOAT_EQ(dimensionStats.fSumWX4, 5846915.6); } { - const auto &dimensionStats = statsA.GetDimensionStats(1); + const auto &dimensionStats = statsA.GetDimensionStats(2); EXPECT_FLOAT_EQ(dimensionStats.fSumWX, 659.3); EXPECT_FLOAT_EQ(dimensionStats.fSumWX2, 21850); EXPECT_FLOAT_EQ(dimensionStats.fSumWX3, 842029.46); @@ -119,19 +139,25 @@ TEST(RHistStats, AddDifferent) { RHistStats statsA(2); RHistStats statsB(3); + RHistStats statsC(3); + statsC.DisableDimension(1); EXPECT_THROW(statsA.Add(statsB), std::invalid_argument); + EXPECT_THROW(statsB.Add(statsC), std::invalid_argument); + EXPECT_THROW(statsC.Add(statsB), std::invalid_argument); } TEST(RHistStats, AddAtomic) { - RHistStats statsA(2); - RHistStats statsB(2); + RHistStats statsA(3); + RHistStats statsB(3); + statsA.DisableDimension(1); + statsB.DisableDimension(1); static constexpr std::size_t Entries = 20; for (std::size_t i = 0; i < Entries; i++) { - statsA.Fill(i, 2 * i); - statsB.Fill(2 * i, 3 * i, RWeight(0.1 + 0.03 * i)); + statsA.Fill(i, 1, 2 * i); + statsB.Fill(2 * i, 1, 3 * i, RWeight(0.1 + 0.03 * i)); } statsA.AddAtomic(statsB); @@ -148,7 +174,7 @@ TEST(RHistStats, AddAtomic) EXPECT_FLOAT_EQ(dimensionStats.fSumWX4, 5846915.6); } { - const auto &dimensionStats = statsA.GetDimensionStats(1); + const auto &dimensionStats = statsA.GetDimensionStats(2); EXPECT_FLOAT_EQ(dimensionStats.fSumWX, 659.3); EXPECT_FLOAT_EQ(dimensionStats.fSumWX2, 21850); EXPECT_FLOAT_EQ(dimensionStats.fSumWX3, 842029.46); @@ -190,8 +216,12 @@ TEST(RHistStats, AddAtomicDifferent) { RHistStats statsA(2); RHistStats statsB(3); + RHistStats statsC(3); + statsC.DisableDimension(1); EXPECT_THROW(statsA.AddAtomic(statsB), std::invalid_argument); + EXPECT_THROW(statsB.AddAtomic(statsC), std::invalid_argument); + EXPECT_THROW(statsC.AddAtomic(statsB), std::invalid_argument); } TEST(RHistStats, Clear) @@ -501,6 +531,14 @@ TEST(RHistStats, FillInvalidNumberOfArguments) EXPECT_THROW(stats2.Fill(1, 2, 3), std::invalid_argument); } +TEST(RHistStats, FillInvalidArgumentType) +{ + RHistStats stats(1); + + EXPECT_NO_THROW(stats.Fill(1)); + EXPECT_THROW(stats.Fill("a"), std::invalid_argument); +} + TEST(RHistStats, FillWeightInvalidNumberOfArguments) { RHistStats stats1(1); @@ -514,6 +552,14 @@ TEST(RHistStats, FillWeightInvalidNumberOfArguments) EXPECT_THROW(stats2.Fill(1, 2, 3, RWeight(1)), std::invalid_argument); } +TEST(RHistStats, FillWeightInvalidArgumentType) +{ + RHistStats stats(1); + + EXPECT_NO_THROW(stats.Fill(1, RWeight(1))); + EXPECT_THROW(stats.Fill("a", RWeight(1)), std::invalid_argument); +} + TEST(RHistStats, FillTupleWeightInvalidNumberOfArguments) { RHistStats stats1(1);