From 778fddf6a1b0f85df5e4cffc7ba3b44568df0c67 Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 11:09:28 +0200 Subject: [PATCH 1/7] [hist] Fix static_assert in test for RHist --- hist/histv7/test/hist_hist.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hist/histv7/test/hist_hist.cxx b/hist/histv7/test/hist_hist.cxx index 2cc6d2190f2cf..dac4669a9ae57 100644 --- a/hist/histv7/test/hist_hist.cxx +++ b/hist/histv7/test/hist_hist.cxx @@ -12,8 +12,8 @@ // Mostly, RHist = RHistEngine + RHistStats which are tested individually. Here we mostly check that the forwarding // works correctly. -static_assert(std::is_nothrow_move_constructible_v>); -static_assert(std::is_nothrow_move_assignable_v>); +static_assert(std::is_nothrow_move_constructible_v>); +static_assert(std::is_nothrow_move_assignable_v>); TEST(RHist, Constructor) { From a2c599680fbe90e631eb581afd2070da50fa3979 Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 13:09:49 +0200 Subject: [PATCH 2/7] [hist] Simplify dimension check in variadic Fill --- hist/histv7/inc/ROOT/RHistEngine.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hist/histv7/inc/ROOT/RHistEngine.hxx b/hist/histv7/inc/ROOT/RHistEngine.hxx index f438e60b2c789..ff94af75eef41 100644 --- a/hist/histv7/inc/ROOT/RHistEngine.hxx +++ b/hist/histv7/inc/ROOT/RHistEngine.hxx @@ -470,7 +470,7 @@ public: if constexpr (std::is_same_v::type, RWeight>) { static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types"); static constexpr std::size_t N = sizeof...(A) - 1; - if (N != fAxes.GetNDimensions()) { + if (N != GetNDimensions()) { throw std::invalid_argument("invalid number of arguments to Fill"); } RWeight weight = std::get(t); @@ -566,7 +566,7 @@ public: if constexpr (std::is_same_v::type, RWeight>) { static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types"); static constexpr std::size_t N = sizeof...(A) - 1; - if (N != fAxes.GetNDimensions()) { + if (N != GetNDimensions()) { throw std::invalid_argument("invalid number of arguments to Fill"); } RWeight weight = std::get(t); From 0cbbfbe8da8f3c2952a41c8178e2a7d2d0e6918f Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 16:52:49 +0200 Subject: [PATCH 3/7] [hist] Remove unused include in RHist.hxx --- hist/histv7/inc/ROOT/RHist.hxx | 1 - 1 file changed, 1 deletion(-) diff --git a/hist/histv7/inc/ROOT/RHist.hxx b/hist/histv7/inc/ROOT/RHist.hxx index 207240353c2bb..a3e0c112a76ca 100644 --- a/hist/histv7/inc/ROOT/RHist.hxx +++ b/hist/histv7/inc/ROOT/RHist.hxx @@ -8,7 +8,6 @@ #include "RAxisVariant.hxx" #include "RBinIndex.hxx" #include "RBinIndexMultiDimRange.hxx" -#include "RCategoricalAxis.hxx" #include "RHistEngine.hxx" #include "RHistStats.hxx" #include "RRegularAxis.hxx" From 13e2cd3a9c2b3b12a3ec9a59b9ff3756477c8d8c Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 11:49:59 +0200 Subject: [PATCH 4/7] [hist] Implement initial RProfile It builds on RHistEngine with an additional variable. --- hist/histv7/doc/CodeArchitecture.md | 4 + hist/histv7/headers.cmake | 1 + hist/histv7/inc/LinkDef.h | 4 + hist/histv7/inc/ROOT/RProfile.hxx | 303 ++++++++++++++++++++++++++++ hist/histv7/test/CMakeLists.txt | 1 + hist/histv7/test/hist_io.cxx | 7 + hist/histv7/test/hist_profile.cxx | 154 ++++++++++++++ hist/histv7/test/hist_test.hxx | 2 + 8 files changed, 476 insertions(+) create mode 100644 hist/histv7/inc/ROOT/RProfile.hxx create mode 100644 hist/histv7/test/hist_profile.cxx diff --git a/hist/histv7/doc/CodeArchitecture.md b/hist/histv7/doc/CodeArchitecture.md index 0cb28176b5150..87850b4b8e602 100644 --- a/hist/histv7/doc/CodeArchitecture.md +++ b/hist/histv7/doc/CodeArchitecture.md @@ -50,6 +50,10 @@ It also keeps statistics of the unbinned values for each dimension. This class combines `RHistEngine`, with its axes and bin contents, and `RHistStats`. During `Fill`, it delegates to `RHistEngine::Fill` but also updates the histogram statistics. +### `RProfile` + +This class builds on `RHistEngine`, but requires an additional value during `Fill`. + ## Classes for Weighted Filling ### `RBinWithError` diff --git a/hist/histv7/headers.cmake b/hist/histv7/headers.cmake index 465b2bd51930a..a42cb5d002116 100644 --- a/hist/histv7/headers.cmake +++ b/hist/histv7/headers.cmake @@ -14,6 +14,7 @@ set(histv7_headers ROOT/RHistStats.hxx ROOT/RHistUtils.hxx ROOT/RLinearizedIndex.hxx + ROOT/RProfile.hxx ROOT/RRegularAxis.hxx ROOT/RSliceBinIndexMapper.hxx ROOT/RSliceSpec.hxx diff --git a/hist/histv7/inc/LinkDef.h b/hist/histv7/inc/LinkDef.h index 23a702fab6d2a..2cdd338872abc 100644 --- a/hist/histv7/inc/LinkDef.h +++ b/hist/histv7/inc/LinkDef.h @@ -20,6 +20,10 @@ #pragma link C++ class ROOT::Experimental::RHistEngine-; #pragma link C++ class ROOT::Experimental::RHistEngine-; +#pragma link C++ class ROOT::Experimental::RProfile::RProfileBin-; +#pragma link C++ class ROOT::Experimental::RHistEngine-; +#pragma link C++ class ROOT::Experimental::RProfile-; + #pragma link C++ class ROOT::Experimental::RAxisVariant-; #pragma link C++ class ROOT::Experimental::RCategoricalAxis-; #pragma link C++ class ROOT::Experimental::RHistStats-; diff --git a/hist/histv7/inc/ROOT/RProfile.hxx b/hist/histv7/inc/ROOT/RProfile.hxx new file mode 100644 index 0000000000000..71e5294c7709d --- /dev/null +++ b/hist/histv7/inc/ROOT/RProfile.hxx @@ -0,0 +1,303 @@ +/// \file +/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes. +/// Feedback is welcome! + +#ifndef ROOT_RProfile +#define ROOT_RProfile + +#include "RAxisVariant.hxx" +#include "RBinIndex.hxx" +#include "RHistEngine.hxx" +#include "RRegularAxis.hxx" +#include "RWeight.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +class TBuffer; + +namespace ROOT { +namespace Experimental { + +/** +A profile histogram, computing statistical quantities of an additional variable per bin. + +Calling \ref Fill(const std::tuple &args, const V &value) "Fill" requires an additional value: +\code +ROOT::Experimental::RProfile profile(10, {5, 15}); +profile.Fill(std::make_tuple(8.2), 23.0); +profile.Fill(std::make_tuple(8.7), 25.0); +// Bin 3 has a mean of 24.0 and a standard deviation of 1.0 +\endcode + +The class is not templated, the bin content is always of type RProfileBin. + +An object can have arbitrary dimensionality determined at run-time. The axis configuration is passed as a vector of +RAxisVariant: +\code +std::vector axes; +axes.push_back(ROOT::Experimental::RRegularAxis(10, {5, 15})); +axes.push_back(ROOT::Experimental::RVariableBinAxis({1, 10, 100, 1000})); +ROOT::Experimental::RProfile profile(axes); +// profile.GetNDimensions() will return 2 +\endcode + +\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes. +Feedback is welcome! +*/ +class RProfile final { + struct RValueWrapper { + double fValue; + + explicit RValueWrapper(double value) : fValue(value) {} + }; + + struct RValueWeightWrapper { + double fValue; + double fWeight; + + explicit RValueWeightWrapper(double value, double weight) : fValue(value), fWeight(weight) {} + }; + +public: + /// The bin content type of a profile histogram. + struct RProfileBin final { + double fSumValues = 0; + double fSumValues2 = 0; + double fSum = 0; + double fSum2 = 0; + + RProfileBin &operator+=(const RValueWrapper &rhs) + { + fSumValues += rhs.fValue; + fSumValues2 += rhs.fValue * rhs.fValue; + fSum++; + fSum2++; + return *this; + } + + RProfileBin &operator+=(const RValueWeightWrapper &rhs) + { + fSumValues += rhs.fWeight * rhs.fValue; + fSumValues2 += rhs.fWeight * rhs.fValue * rhs.fValue; + fSum += rhs.fWeight; + fSum2 += rhs.fWeight * rhs.fWeight; + return *this; + } + }; + +private: + /// The histogram engine including the bin contents. + RHistEngine fEngine; + +public: + /// Construct a profile histogram. + /// + /// \param[in] axes the axis objects, must have size > 0 + explicit RProfile(std::vector axes) : fEngine(std::move(axes)) {} + + /// Construct a profile histogram. + /// + /// Note that there is no perfect forwarding of the axis objects. If that is needed, use the + /// \ref RHist(std::vector axes) "overload accepting a std::vector". + /// + /// \param[in] axes the axis objects, must have size > 0 + explicit RProfile(std::initializer_list axes) : RProfile(std::vector(axes)) {} + + /// Construct a profile histogram. + /// + /// Note that there is no perfect forwarding of the axis objects. If that is needed, use the + /// \ref RHist(std::vector axes) "overload accepting a std::vector". + /// + /// \param[in] axis1 the first axis object + /// \param[in] axes the remaining axis objects + template + explicit RProfile(const RAxisVariant &axis1, const Axes &...axes) + : RProfile(std::vector{axis1, axes...}) + { + } + + /// Construct a one-dimensional profile histogram with a regular axis. + /// + /// \param[in] nNormalBins the number of normal bins, must be > 0 + /// \param[in] interval the axis interval (lower end inclusive, upper end exclusive) + /// \par See also + /// the \ref RRegularAxis::RRegularAxis(std::uint64_t nNormalBins, std::pair interval, bool + /// enableFlowBins) "constructor of RRegularAxis" + RProfile(std::uint64_t nNormalBins, std::pair interval) + : RProfile(std::vector{RRegularAxis(nNormalBins, interval)}) + { + } + + /// The copy constructor is deleted. + /// + /// Copying all bin contents can be an expensive operation, depending on the number of bins. If required, users can + /// explicitly call Clone(). + RProfile(const RProfile &) = delete; + /// Efficiently move construct a profile histogram. + /// + /// After this operation, the moved-from object is invalid. + RProfile(RProfile &&) = default; + + /// The copy assignment operator is deleted. + /// + /// Copying all bin contents can be an expensive operation, depending on the number of bins. If required, users can + /// explicitly call Clone(). + RProfile &operator=(const RProfile &) = delete; + /// Efficiently move a profile histogram. + /// + /// After this operation, the moved-from object is invalid. + RProfile &operator=(RProfile &&) = default; + + ~RProfile() = default; + + /// \name Accessors + /// \{ + + const RHistEngine &GetEngine() const { return fEngine; } + + const std::vector &GetAxes() const { return fEngine.GetAxes(); } + std::size_t GetNDimensions() const { return fEngine.GetNDimensions(); } + std::uint64_t GetTotalNBins() const { return fEngine.GetTotalNBins(); } + + /// Get the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// std::array indices = {3, 5}; + /// const auto &content = profile.GetBinContent(indices); + /// \endcode + /// + /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special + /// values. See also the class documentation of RBinIndex. + /// + /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found. + /// + /// \param[in] indices the array of indices for each axis + /// \return the bin content + /// \par See also + /// the \ref GetBinContent(const A &... args) const "variadic function template overload" accepting arguments + /// directly + template + const RProfileBin &GetBinContent(const std::array &indices) const + { + return fEngine.GetBinContent(indices); + } + + /// Get the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// std::vector indices = {3, 5}; + /// const auto &content = profile.GetBinContent(indices); + /// \endcode + /// + /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special + /// values. See also the class documentation of RBinIndex. + /// + /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found. + /// + /// \param[in] indices the array of indices for each axis + /// \return the bin content + /// \par See also + /// the \ref GetBinContent(const A &... args) const "variadic function template overload" accepting arguments + /// directly + const RProfileBin &GetBinContent(const std::vector &indices) const + { + return fEngine.GetBinContent(indices); + } + + /// Get the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// const auto &content = profile.GetBinContent(3, 5); + /// \endcode + /// + /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special + /// values. See also the class documentation of RBinIndex. + /// + /// Throws an exception if the number of arguments does not match the axis configuration or the bin is not found. + /// + /// \param[in] args the arguments for each axis + /// \return the bin content + /// \par See also + /// the function overloads accepting \ref GetBinContent(const std::array &indices) const "`std::array`" + /// or \ref GetBinContent(const std::vector &indices) const "`std::vector`" + template + const RProfileBin &GetBinContent(const A &...args) const + { + return fEngine.GetBinContent(args...); + } + + /// \} + /// \name Filling + /// \{ + + /// Fill an entry into the profile histogram. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// auto args = std::make_tuple(8.5, 10.5); + /// profile.Fill(args, 23.0); + /// \endcode + /// + /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently + /// discarded. + /// + /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be + /// converted for the axis type at run-time. + /// + /// \param[in] args the arguments for each axis + /// \param[in] v the additional argument + /// \par See also + /// the \ref Fill(const std::tuple &args, const V &value, RWeight weight) "overload for weighted filling" + template + void Fill(const std::tuple &args, const V &value) + { + RValueWrapper wrapper(value); + fEngine.Fill(args, wrapper); + } + + /// Fill an entry into the profile histogram with a weight. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// auto args = std::make_tuple(8.5, 10.5); + /// profile.Fill(args, 23.0, ROOT::Experimental::RWeight(0.8)); + /// \endcode + /// + /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently + /// discarded. + /// + /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be + /// converted for the axis type at run-time. + /// + /// \param[in] args the arguments for each axis + /// \param[in] v the additional argument + /// \param[in] weight the weight for this entry + /// \par See also + /// the \ref Fill(const std::tuple &args, const V &value) "overload for unweighted filling" + template + void Fill(const std::tuple &args, const V &value, RWeight weight) + { + RValueWeightWrapper wrapper(value, weight.fValue); + fEngine.Fill(args, wrapper); + } + + /// \} + + /// %ROOT Streamer function to throw when trying to store an object of this class. + void Streamer(TBuffer &) { throw std::runtime_error("unable to store RHist"); } +}; + +} // namespace Experimental +} // namespace ROOT + +#endif diff --git a/hist/histv7/test/CMakeLists.txt b/hist/histv7/test/CMakeLists.txt index ae8fb5ba6e0b6..eeb9a3463878e 100644 --- a/hist/histv7/test/CMakeLists.txt +++ b/hist/histv7/test/CMakeLists.txt @@ -8,6 +8,7 @@ HIST_ADD_GTEST(hist_engine hist_engine.cxx) 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_profile hist_profile.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) diff --git a/hist/histv7/test/hist_io.cxx b/hist/histv7/test/hist_io.cxx index 185a22287842c..c1e267429378a 100644 --- a/hist/histv7/test/hist_io.cxx +++ b/hist/histv7/test/hist_io.cxx @@ -144,3 +144,10 @@ TEST(RHist, Streamer) const RHist histE({axis}); ExpectThrowOnWriteObject(histE); } + +TEST(RProfile, Streamer) +{ + static constexpr std::size_t Bins = 20; + const RProfile profile(Bins, {0, Bins}); + ExpectThrowOnWriteObject(profile); +} diff --git a/hist/histv7/test/hist_profile.cxx b/hist/histv7/test/hist_profile.cxx new file mode 100644 index 0000000000000..ffba82ed55df8 --- /dev/null +++ b/hist/histv7/test/hist_profile.cxx @@ -0,0 +1,154 @@ +#include "hist_test.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include + +static_assert(std::is_nothrow_move_constructible_v); +static_assert(std::is_nothrow_move_assignable_v); + +TEST(RProfile, Constructor) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis regularAxis(Bins, {0, Bins}); + + // The most generic constructor takes a vector of axis objects. + const std::vector axes = {regularAxis, regularAxis}; + RProfile profile(axes); + EXPECT_EQ(profile.GetNDimensions(), 2); + const auto &engine = profile.GetEngine(); + EXPECT_EQ(engine.GetNDimensions(), 2); + EXPECT_EQ(profile.GetAxes().size(), 2); + // Both axes include underflow and overflow bins. + EXPECT_EQ(profile.GetTotalNBins(), (Bins + 2) * (Bins + 2)); + + // Test other constructors, including move-assignment. + profile = RProfile(Bins, {0, Bins}); + ASSERT_EQ(profile.GetNDimensions(), 1); + auto *regular = profile.GetAxes()[0].GetRegularAxis(); + ASSERT_TRUE(regular != nullptr); + EXPECT_EQ(regular->GetNNormalBins(), Bins); + EXPECT_EQ(regular->GetLow(), 0); + EXPECT_EQ(regular->GetHigh(), Bins); + // std::make_pair will take the types of the arguments, std::size_t in this case. + profile = RProfile(Bins, std::make_pair(0, Bins)); + EXPECT_EQ(profile.GetNDimensions(), 1); + + // Brace-enclosed initializer list + profile = RProfile({regularAxis}); + EXPECT_EQ(profile.GetNDimensions(), 1); + profile = RProfile({regularAxis, regularAxis}); + EXPECT_EQ(profile.GetNDimensions(), 2); + + // Templated constructors + profile = RProfile(regularAxis); + EXPECT_EQ(profile.GetNDimensions(), 1); + profile = RProfile(regularAxis, regularAxis); + EXPECT_EQ(profile.GetNDimensions(), 2); + profile = RProfile(regularAxis, regularAxis, regularAxis); + EXPECT_EQ(profile.GetNDimensions(), 3); +} + +TEST(RProfile, Fill) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + + profile.Fill(std::make_tuple(9.5), 25.0); + + std::array indices = {9}; + auto &bin9 = profile.GetBinContent(indices); + EXPECT_EQ(bin9.fSumValues, 25.0); + EXPECT_EQ(bin9.fSumValues2, 625.0); + EXPECT_EQ(bin9.fSum, 1.0); + EXPECT_EQ(bin9.fSum2, 1.0); +} + +TEST(RProfile, FillWeight) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + + profile.Fill(std::make_tuple(9.5), 25.0, RWeight(0.9)); + + std::array indices = {9}; + auto &bin9 = profile.GetBinContent(indices); + EXPECT_FLOAT_EQ(bin9.fSumValues, 22.5); + EXPECT_FLOAT_EQ(bin9.fSumValues2, 562.5); + EXPECT_FLOAT_EQ(bin9.fSum, 0.9); + EXPECT_FLOAT_EQ(bin9.fSum2, 0.81); +} + +TEST(RProfile, FillCategorical) +{ + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + RProfile profile({axis}); + + profile.Fill(std::make_tuple("c"), 25.0); + + std::array indices = {2}; + auto &bin2 = profile.GetBinContent(indices); + EXPECT_EQ(bin2.fSumValues, 25.0); + EXPECT_EQ(bin2.fSumValues2, 625.0); + EXPECT_EQ(bin2.fSum, 1.0); + EXPECT_EQ(bin2.fSum2, 1.0); +} + +TEST(RProfile, FillCategoricalWeight) +{ + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + RProfile profile({axis}); + + profile.Fill(std::make_tuple("c"), 25.0, RWeight(0.9)); + + std::array indices = {2}; + auto &bin2 = profile.GetBinContent(indices); + EXPECT_FLOAT_EQ(bin2.fSumValues, 22.5); + EXPECT_FLOAT_EQ(bin2.fSumValues2, 562.5); + EXPECT_FLOAT_EQ(bin2.fSum, 0.9); + EXPECT_FLOAT_EQ(bin2.fSum2, 0.81); +} + +class CopyArgument final { + static int gCopies; + + double fArgument; + +public: + CopyArgument(double argument) : fArgument(argument) {} + CopyArgument(const CopyArgument &) { gCopies++; } + CopyArgument(CopyArgument &&) = default; + CopyArgument &operator=(const CopyArgument &) + { + gCopies++; + return *this; + } + CopyArgument &operator=(CopyArgument &&) = default; + + operator double() const { return fArgument; } + + static bool HasBeenCopied() { return gCopies > 0; } +}; + +int CopyArgument::gCopies = 0; + +TEST(RProfile, FillForward) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + CopyArgument value(23.0); + + std::tuple args(1.5); + profile.Fill(args, value); + profile.Fill(args, value, RWeight(0.5)); + EXPECT_EQ(profile.GetBinContent(1).fSumValues, 34.5); + + ASSERT_FALSE(CopyArgument::HasBeenCopied()); +} diff --git a/hist/histv7/test/hist_test.hxx b/hist/histv7/test/hist_test.hxx index bec8784885387..c2ef159821a7f 100644 --- a/hist/histv7/test/hist_test.hxx +++ b/hist/histv7/test/hist_test.hxx @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -31,6 +32,7 @@ using ROOT::Experimental::RHistAutoAxisFiller; using ROOT::Experimental::RHistConcurrentFiller; using ROOT::Experimental::RHistEngine; using ROOT::Experimental::RHistStats; +using ROOT::Experimental::RProfile; using ROOT::Experimental::RRegularAxis; using ROOT::Experimental::RSliceSpec; using ROOT::Experimental::RVariableBinAxis; From 654d8f2ba30c766a716cbc4c892b9b44d3975fbb Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 13:35:41 +0200 Subject: [PATCH 5/7] [hist] Implement variadic RProfile::Fill --- hist/histv7/inc/ROOT/RHistEngine.hxx | 27 +++++++++--- hist/histv7/inc/ROOT/RProfile.hxx | 62 +++++++++++++++++++++++--- hist/histv7/test/hist_profile.cxx | 66 ++++++++++++++++++++++++++++ 3 files changed, 145 insertions(+), 10 deletions(-) diff --git a/hist/histv7/inc/ROOT/RHistEngine.hxx b/hist/histv7/inc/ROOT/RHistEngine.hxx index ff94af75eef41..ac7912c05f165 100644 --- a/hist/histv7/inc/ROOT/RHistEngine.hxx +++ b/hist/histv7/inc/ROOT/RHistEngine.hxx @@ -75,6 +75,8 @@ class RHistEngine final { // For slicing, RHist needs to call SliceImpl. friend class RHist; + friend class RProfile; + /// The axis configuration for this histogram. Relevant methods are forwarded from the public interface. Internal::RAxes fAxes; /// The bin contents for this histogram @@ -407,6 +409,25 @@ public: } } + /// \} + // End the group to ensure that all contained member functions are public. + +private: + // Also used by RProfile::Fill(const A &... args) + template + void FillImpl(const std::tuple &args, const W &weight) + { + RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl(args); + if (index.fValid) { + assert(index.fIndex < fBinContents.size()); + fBinContents[index.fIndex] += weight; + } + } + +public: + /// \name Filling + /// \{ + /// Fill an entry into the histogram with a user-defined weight. /// /// This overload is only available for user-defined bin content types. @@ -430,11 +451,7 @@ public: if (sizeof...(A) != GetNDimensions()) { throw std::invalid_argument("invalid number of arguments to Fill"); } - RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl(args); - if (index.fValid) { - assert(index.fIndex < fBinContents.size()); - fBinContents[index.fIndex] += weight; - } + FillImpl(args, weight); } /// Fill an entry into the histogram. diff --git a/hist/histv7/inc/ROOT/RProfile.hxx b/hist/histv7/inc/ROOT/RProfile.hxx index 71e5294c7709d..49377f8ee8b24 100644 --- a/hist/histv7/inc/ROOT/RProfile.hxx +++ b/hist/histv7/inc/ROOT/RProfile.hxx @@ -8,6 +8,7 @@ #include "RAxisVariant.hxx" #include "RBinIndex.hxx" #include "RHistEngine.hxx" +#include "RHistUtils.hxx" #include "RRegularAxis.hxx" #include "RWeight.hxx" @@ -17,6 +18,7 @@ #include #include #include +#include #include #include @@ -28,11 +30,11 @@ namespace Experimental { /** A profile histogram, computing statistical quantities of an additional variable per bin. -Calling \ref Fill(const std::tuple &args, const V &value) "Fill" requires an additional value: +Calling \ref Fill(const A &... args) "Fill" requires an additional value: \code ROOT::Experimental::RProfile profile(10, {5, 15}); -profile.Fill(std::make_tuple(8.2), 23.0); -profile.Fill(std::make_tuple(8.7), 25.0); +profile.Fill(8.2, 23.0); +profile.Fill(8.7, 25.0); // Bin 3 has a mean of 24.0 and a standard deviation of 1.0 \endcode @@ -257,7 +259,8 @@ public: /// \param[in] args the arguments for each axis /// \param[in] v the additional argument /// \par See also - /// the \ref Fill(const std::tuple &args, const V &value, RWeight weight) "overload for weighted filling" + /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the + /// \ref Fill(const std::tuple &args, const V &value, RWeight weight) "overload for weighted filling" template void Fill(const std::tuple &args, const V &value) { @@ -283,7 +286,8 @@ public: /// \param[in] v the additional argument /// \param[in] weight the weight for this entry /// \par See also - /// the \ref Fill(const std::tuple &args, const V &value) "overload for unweighted filling" + /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the + /// \ref Fill(const std::tuple &args, const V &value) "overload for unweighted filling" template void Fill(const std::tuple &args, const V &value, RWeight weight) { @@ -291,6 +295,54 @@ public: fEngine.Fill(args, wrapper); } + /// Fill an entry into the profile histogram. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// profile.Fill(8.5, 10.5, 23.0); + /// \endcode + /// + /// For weighted filling, pass an RWeight as the last argument: + /// \code + /// profile.Fill(8.5, 10.5, 23.0, ROOT::Experimental::RWeight(0.8)); + /// \endcode + /// + /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently + /// discarded. + /// + /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be + /// converted for the axis type at run-time. + /// + /// \param[in] args the arguments for each axis + /// \par See also + /// the function overloads accepting `std::tuple` + /// \ref Fill(const std::tuple &args, const V &value) "for unweighted filling" and + /// \ref Fill(const std::tuple &args, const V &value, RWeight weight) "for weighted filling" + template + void Fill(const A &...args) + { + static_assert(sizeof...(A) >= 2, "need at least two arguments to Fill"); + if constexpr (sizeof...(A) >= 2) { + auto t = std::forward_as_tuple(args...); + if constexpr (std::is_same_v::type, RWeight>) { + static constexpr std::size_t N = sizeof...(A) - 2; + if (N != GetNDimensions()) { + throw std::invalid_argument("invalid number of arguments to Fill"); + } + RWeight weight = std::get(t); + RValueWeightWrapper wrapper(std::get(t), weight.fValue); + fEngine.FillImpl(t, wrapper); + } else { + static constexpr std::size_t N = sizeof...(A) - 1; + if (N != GetNDimensions()) { + throw std::invalid_argument("invalid number of arguments to Fill"); + } + RValueWrapper wrapper(std::get(t)); + fEngine.FillImpl(t, wrapper); + } + } + } + /// \} /// %ROOT Streamer function to throw when trying to store an object of this class. diff --git a/hist/histv7/test/hist_profile.cxx b/hist/histv7/test/hist_profile.cxx index ffba82ed55df8..6a4cd2d583ac3 100644 --- a/hist/histv7/test/hist_profile.cxx +++ b/hist/histv7/test/hist_profile.cxx @@ -59,8 +59,14 @@ TEST(RProfile, Fill) static constexpr std::size_t Bins = 20; RProfile profile(Bins, {0, Bins}); + profile.Fill(8.5, 23.0); profile.Fill(std::make_tuple(9.5), 25.0); + auto &bin8 = profile.GetBinContent(RBinIndex(8)); + EXPECT_EQ(bin8.fSumValues, 23.0); + EXPECT_EQ(bin8.fSumValues2, 529.0); + EXPECT_EQ(bin8.fSum, 1.0); + EXPECT_EQ(bin8.fSum2, 1.0); std::array indices = {9}; auto &bin9 = profile.GetBinContent(indices); EXPECT_EQ(bin9.fSumValues, 25.0); @@ -69,13 +75,36 @@ TEST(RProfile, Fill) EXPECT_EQ(bin9.fSum2, 1.0); } +TEST(RProfile, FillInvalidNumberOfArguments) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profile1(axis); + ASSERT_EQ(profile1.GetNDimensions(), 1); + RProfile profile2(axis, axis); + ASSERT_EQ(profile2.GetNDimensions(), 2); + + EXPECT_NO_THROW(profile1.Fill(1, 2)); + EXPECT_THROW(profile1.Fill(1, 2, 3), std::invalid_argument); + + EXPECT_THROW(profile2.Fill(1, 2), std::invalid_argument); + EXPECT_NO_THROW(profile2.Fill(1, 2, 3)); + EXPECT_THROW(profile2.Fill(1, 2, 3, 4), std::invalid_argument); +} + TEST(RProfile, FillWeight) { static constexpr std::size_t Bins = 20; RProfile profile(Bins, {0, Bins}); + profile.Fill(8.5, 23.0, RWeight(0.8)); profile.Fill(std::make_tuple(9.5), 25.0, RWeight(0.9)); + auto &bin8 = profile.GetBinContent(RBinIndex(8)); + EXPECT_FLOAT_EQ(bin8.fSumValues, 18.4); + EXPECT_FLOAT_EQ(bin8.fSumValues2, 423.2); + EXPECT_FLOAT_EQ(bin8.fSum, 0.8); + EXPECT_FLOAT_EQ(bin8.fSum2, 0.64); std::array indices = {9}; auto &bin9 = profile.GetBinContent(indices); EXPECT_FLOAT_EQ(bin9.fSumValues, 22.5); @@ -84,14 +113,38 @@ TEST(RProfile, FillWeight) EXPECT_FLOAT_EQ(bin9.fSum2, 0.81); } +TEST(RProfile, FillWeightInvalidNumberOfArguments) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profile1(axis); + ASSERT_EQ(profile1.GetNDimensions(), 1); + RProfile profile2(axis, axis); + ASSERT_EQ(profile2.GetNDimensions(), 2); + + EXPECT_THROW(profile1.Fill(1, RWeight(1)), std::invalid_argument); + EXPECT_NO_THROW(profile1.Fill(1, 2, RWeight(1))); + EXPECT_THROW(profile1.Fill(1, 2, 3, RWeight(1)), std::invalid_argument); + + EXPECT_THROW(profile2.Fill(1, 2, RWeight(1)), std::invalid_argument); + EXPECT_NO_THROW(profile2.Fill(1, 2, 3, RWeight(1))); + EXPECT_THROW(profile2.Fill(1, 2, 3, 4, RWeight(1)), std::invalid_argument); +} + TEST(RProfile, FillCategorical) { const std::vector categories = {"a", "b", "c"}; const RCategoricalAxis axis(categories); RProfile profile({axis}); + profile.Fill("b", 23.0); profile.Fill(std::make_tuple("c"), 25.0); + auto &bin1 = profile.GetBinContent(RBinIndex(1)); + EXPECT_EQ(bin1.fSumValues, 23.0); + EXPECT_EQ(bin1.fSumValues2, 529.0); + EXPECT_EQ(bin1.fSum, 1.0); + EXPECT_EQ(bin1.fSum2, 1.0); std::array indices = {2}; auto &bin2 = profile.GetBinContent(indices); EXPECT_EQ(bin2.fSumValues, 25.0); @@ -106,8 +159,14 @@ TEST(RProfile, FillCategoricalWeight) const RCategoricalAxis axis(categories); RProfile profile({axis}); + profile.Fill("b", 23.0, RWeight(0.8)); profile.Fill(std::make_tuple("c"), 25.0, RWeight(0.9)); + auto &bin1 = profile.GetBinContent(RBinIndex(1)); + EXPECT_FLOAT_EQ(bin1.fSumValues, 18.4); + EXPECT_FLOAT_EQ(bin1.fSumValues2, 423.2); + EXPECT_FLOAT_EQ(bin1.fSum, 0.8); + EXPECT_FLOAT_EQ(bin1.fSum2, 0.64); std::array indices = {2}; auto &bin2 = profile.GetBinContent(indices); EXPECT_FLOAT_EQ(bin2.fSumValues, 22.5); @@ -151,4 +210,11 @@ TEST(RProfile, FillForward) EXPECT_EQ(profile.GetBinContent(1).fSumValues, 34.5); ASSERT_FALSE(CopyArgument::HasBeenCopied()); + + CopyArgument arg(2.5); + profile.Fill(arg, value); + profile.Fill(arg, value, RWeight(0.5)); + EXPECT_EQ(profile.GetBinContent(2).fSumValues, 34.5); + + ASSERT_FALSE(CopyArgument::HasBeenCopied()); } From 18c58211c564b3a91107c4dc947e129d2a8b7b3c Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Thu, 11 Jun 2026 14:01:59 +0200 Subject: [PATCH 6/7] [hist] Integrate RHistStats into RProfile --- hist/histv7/inc/ROOT/RHistUtils.hxx | 15 +++++++ hist/histv7/inc/ROOT/RProfile.hxx | 37 +++++++++++++++- hist/histv7/test/hist_profile.cxx | 68 ++++++++++++++++++++++++++--- 3 files changed, 112 insertions(+), 8 deletions(-) diff --git a/hist/histv7/inc/ROOT/RHistUtils.hxx b/hist/histv7/inc/ROOT/RHistUtils.hxx index 452819ab4739b..6f4d421bfd9c8 100644 --- a/hist/histv7/inc/ROOT/RHistUtils.hxx +++ b/hist/histv7/inc/ROOT/RHistUtils.hxx @@ -5,7 +5,9 @@ #ifndef ROOT_RHistUtils #define ROOT_RHistUtils +#include #include +#include #ifdef _MSC_VER #include @@ -25,6 +27,19 @@ struct LastType { using type = T; }; +template +std::tuple +AppendReferenceImpl(const std::tuple &t, std::index_sequence, const A &a) +{ + return std::tuple(std::get(t)..., a); +} + +template +std::tuple AppendReference(const std::tuple &t, const A &a) +{ + return AppendReferenceImpl(t, std::make_index_sequence(), a); +} + #ifdef _MSC_VER namespace MSVC { template diff --git a/hist/histv7/inc/ROOT/RProfile.hxx b/hist/histv7/inc/ROOT/RProfile.hxx index 49377f8ee8b24..bdc1ce7fc9188 100644 --- a/hist/histv7/inc/ROOT/RProfile.hxx +++ b/hist/histv7/inc/ROOT/RProfile.hxx @@ -8,6 +8,7 @@ #include "RAxisVariant.hxx" #include "RBinIndex.hxx" #include "RHistEngine.hxx" +#include "RHistStats.hxx" #include "RHistUtils.hxx" #include "RRegularAxis.hxx" #include "RWeight.hxx" @@ -97,12 +98,23 @@ public: private: /// The histogram engine including the bin contents. RHistEngine fEngine; + /// The global histogram statistics. + RHistStats fStats; public: /// Construct a profile histogram. /// /// \param[in] axes the axis objects, must have size > 0 - explicit RProfile(std::vector axes) : fEngine(std::move(axes)) {} + explicit RProfile(std::vector axes) : fEngine(std::move(axes)), fStats(fEngine.GetNDimensions() + 1) + { + // 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 (engineAxes[i].GetCategoricalAxis() != nullptr) { + fStats.DisableDimension(i); + } + } + } /// Construct a profile histogram. /// @@ -163,11 +175,29 @@ public: /// \{ const RHistEngine &GetEngine() const { return fEngine; } + const RHistStats &GetStats() const { return fStats; } const std::vector &GetAxes() const { return fEngine.GetAxes(); } std::size_t GetNDimensions() const { return fEngine.GetNDimensions(); } std::uint64_t GetTotalNBins() const { return fEngine.GetTotalNBins(); } + std::uint64_t GetNEntries() const { return fStats.GetNEntries(); } + + /// \} + /// \name Computations + /// \{ + + /// \copydoc RHistStats::ComputeNEffectiveEntries() + double ComputeNEffectiveEntries() const { return fStats.ComputeNEffectiveEntries(); } + /// \copydoc RHistStats::ComputeMean() + double ComputeMean(std::size_t dim = 0) const { return fStats.ComputeMean(dim); } + /// \copydoc RHistStats::ComputeStdDev() + double ComputeStdDev(std::size_t dim = 0) const { return fStats.ComputeStdDev(dim); } + + /// \} + /// \name Accessors + /// \{ + /// Get the content of a single bin. /// /// \code @@ -266,6 +296,8 @@ public: { RValueWrapper wrapper(value); fEngine.Fill(args, wrapper); + // Avoid a second conversion of value, which we already did in wrapper. + fStats.Fill(Internal::AppendReference(args, wrapper.fValue)); } /// Fill an entry into the profile histogram with a weight. @@ -293,6 +325,8 @@ public: { RValueWeightWrapper wrapper(value, weight.fValue); fEngine.Fill(args, wrapper); + // Avoid a second conversion of value, which we already did in wrapper. + fStats.Fill(Internal::AppendReference(args, wrapper.fValue), weight); } /// Fill an entry into the profile histogram. @@ -340,6 +374,7 @@ public: RValueWrapper wrapper(std::get(t)); fEngine.FillImpl(t, wrapper); } + fStats.Fill(args...); } } diff --git a/hist/histv7/test/hist_profile.cxx b/hist/histv7/test/hist_profile.cxx index 6a4cd2d583ac3..a38c0b64ab703 100644 --- a/hist/histv7/test/hist_profile.cxx +++ b/hist/histv7/test/hist_profile.cxx @@ -23,6 +23,8 @@ TEST(RProfile, Constructor) EXPECT_EQ(profile.GetNDimensions(), 2); const auto &engine = profile.GetEngine(); EXPECT_EQ(engine.GetNDimensions(), 2); + const auto &stats = profile.GetStats(); + EXPECT_EQ(stats.GetNDimensions(), 3); EXPECT_EQ(profile.GetAxes().size(), 2); // Both axes include underflow and overflow bins. EXPECT_EQ(profile.GetTotalNBins(), (Bins + 2) * (Bins + 2)); @@ -73,6 +75,13 @@ TEST(RProfile, Fill) EXPECT_EQ(bin9.fSumValues2, 625.0); EXPECT_EQ(bin9.fSum, 1.0); EXPECT_EQ(bin9.fSum2, 1.0); + + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 2); + EXPECT_FLOAT_EQ(profile.ComputeMean(0), 9); + EXPECT_FLOAT_EQ(profile.ComputeStdDev(0), 0.5); + EXPECT_FLOAT_EQ(profile.ComputeMean(1), 24.0); + EXPECT_FLOAT_EQ(profile.ComputeStdDev(1), 1.0); } TEST(RProfile, FillInvalidNumberOfArguments) @@ -111,6 +120,14 @@ TEST(RProfile, FillWeight) EXPECT_FLOAT_EQ(bin9.fSumValues2, 562.5); EXPECT_FLOAT_EQ(bin9.fSum, 0.9); EXPECT_FLOAT_EQ(bin9.fSum2, 0.81); + + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW(), 1.7); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW2(), 1.45); + // Cross-checked with TH1 + EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 1.9931034); + EXPECT_FLOAT_EQ(profile.ComputeMean(0), 9.0294118); + EXPECT_FLOAT_EQ(profile.ComputeStdDev(0), 0.49913420); } TEST(RProfile, FillWeightInvalidNumberOfArguments) @@ -151,6 +168,9 @@ TEST(RProfile, FillCategorical) EXPECT_EQ(bin2.fSumValues2, 625.0); EXPECT_EQ(bin2.fSum, 1.0); EXPECT_EQ(bin2.fSum2, 1.0); + + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 2); } TEST(RProfile, FillCategoricalWeight) @@ -173,6 +193,37 @@ TEST(RProfile, FillCategoricalWeight) EXPECT_FLOAT_EQ(bin2.fSumValues2, 562.5); EXPECT_FLOAT_EQ(bin2.fSum, 0.9); EXPECT_FLOAT_EQ(bin2.fSum2, 0.81); + + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW(), 1.7); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW2(), 1.45); + // Cross-checked with TH1 + EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 1.9931034); +} + +TEST(RProfile, FillExceptionSafety) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profile({axis, axis}); + + profile.Fill(1.5, 2.5, 3.5); + ASSERT_EQ(profile.GetNEntries(), 1); + ASSERT_EQ(profile.GetBinContent(RBinIndex(1), RBinIndex(2)).fSumValues, 3.5); + + EXPECT_THROW(profile.Fill(1.5, "b", 3.5), std::invalid_argument); + EXPECT_THROW(profile.Fill(std::make_tuple(1.5, "b"), 3.5), std::invalid_argument); + EXPECT_THROW(profile.Fill(1.5, "b", 3.5, RWeight(1)), std::invalid_argument); + EXPECT_THROW(profile.Fill(std::make_tuple(1.5, "b"), 3.5, RWeight(1)), std::invalid_argument); + + // Verify exception safety. Only the first entry should be there. + EXPECT_EQ(profile.GetNEntries(), 1); + EXPECT_EQ(profile.GetBinContent(RBinIndex(1), RBinIndex(2)).fSumValues, 3.5); + EXPECT_EQ(profile.GetStats().GetSumW(), 1); + EXPECT_EQ(profile.GetStats().GetSumW2(), 1); + EXPECT_EQ(profile.GetStats().GetDimensionStats(0).fSumWX, 1.5); + EXPECT_EQ(profile.GetStats().GetDimensionStats(1).fSumWX, 2.5); + EXPECT_EQ(profile.GetStats().GetDimensionStats(2).fSumWX, 3.5); } class CopyArgument final { @@ -201,20 +252,23 @@ int CopyArgument::gCopies = 0; TEST(RProfile, FillForward) { static constexpr std::size_t Bins = 20; - RProfile profile(Bins, {0, Bins}); + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profile(axis, axis); CopyArgument value(23.0); - std::tuple args(1.5); + std::tuple args(1.5, 2.5); profile.Fill(args, value); profile.Fill(args, value, RWeight(0.5)); - EXPECT_EQ(profile.GetBinContent(1).fSumValues, 34.5); + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_EQ(profile.GetBinContent(1, 2).fSumValues, 34.5); ASSERT_FALSE(CopyArgument::HasBeenCopied()); - CopyArgument arg(2.5); - profile.Fill(arg, value); - profile.Fill(arg, value, RWeight(0.5)); - EXPECT_EQ(profile.GetBinContent(2).fSumValues, 34.5); + CopyArgument arg1(3.5), arg2(4.5); + profile.Fill(arg1, arg2, value); + profile.Fill(arg1, arg2, value, RWeight(0.5)); + EXPECT_EQ(profile.GetNEntries(), 4); + EXPECT_EQ(profile.GetBinContent(3, 4).fSumValues, 34.5); ASSERT_FALSE(CopyArgument::HasBeenCopied()); } From ce6f93e6d80a05a58ddc63e2c96132c9611e9c53 Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Fri, 12 Jun 2026 09:47:29 +0200 Subject: [PATCH 7/7] [hist] Implement basic methods on RProfile SetBinContent, Add, AddAtomic, Clear, Clone, Scale; implemented similarly to RHist, forwarding to RHistEngine and RHistStats. --- hist/histv7/inc/ROOT/RProfile.hxx | 147 ++++++++++++++++++++++ hist/histv7/test/hist_profile.cxx | 202 ++++++++++++++++++++++++++++++ 2 files changed, 349 insertions(+) diff --git a/hist/histv7/inc/ROOT/RProfile.hxx b/hist/histv7/inc/ROOT/RProfile.hxx index bdc1ce7fc9188..5d982a6245d68 100644 --- a/hist/histv7/inc/ROOT/RProfile.hxx +++ b/hist/histv7/inc/ROOT/RProfile.hxx @@ -7,6 +7,7 @@ #include "RAxisVariant.hxx" #include "RBinIndex.hxx" +#include "RBinIndexMultiDimRange.hxx" #include "RHistEngine.hxx" #include "RHistStats.hxx" #include "RHistUtils.hxx" @@ -93,6 +94,35 @@ public: fSum2 += rhs.fWeight * rhs.fWeight; return *this; } + + RProfileBin &operator+=(const RProfileBin &rhs) + { + fSumValues += rhs.fSumValues; + fSumValues2 += rhs.fSumValues2; + fSum += rhs.fSum; + fSum2 += rhs.fSum2; + return *this; + } + + RProfileBin &operator*=(double factor) + { + fSumValues *= factor; + fSumValues2 *= factor; + fSum *= factor; + fSum2 *= factor * factor; + return *this; + } + + /// Add another bin content using atomic instructions. + /// + /// \param[in] rhs another bin content that must not be modified during the operation + void AtomicAdd(const RProfileBin &rhs) + { + Internal::AtomicAdd(&fSumValues, rhs.fSumValues); + Internal::AtomicAdd(&fSumValues2, rhs.fSumValues2); + Internal::AtomicAdd(&fSum, rhs.fSum); + Internal::AtomicAdd(&fSum2, rhs.fSum2); + } }; private: @@ -101,6 +131,9 @@ private: /// The global histogram statistics. RHistStats fStats; + /// Private constructor based off an engine. + RProfile(RHistEngine engine) : fEngine(std::move(engine)), fStats(fEngine.GetNDimensions() + 1) {} + public: /// Construct a profile histogram. /// @@ -268,6 +301,66 @@ public: return fEngine.GetBinContent(args...); } + /// Get the multidimensional range of all bins. + /// + /// \return the multidimensional range + RBinIndexMultiDimRange GetFullMultiDimRange() const { return fEngine.GetFullMultiDimRange(); } + + /// Set the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// std::array indices = {3, 5}; + /// ROOT::Experimental::RProfile::RProfileBin value = /* ... */; + /// profile.SetBinContent(indices, value); + /// \endcode + /// + /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special + /// values. See also the class documentation of RBinIndex. + /// + /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found. + /// + /// \warning Setting the bin content will taint the global histogram statistics. Attempting to access its values, for + /// example calling GetNEntries(), will throw exceptions. + /// + /// \param[in] indices the array of indices for each axis + /// \param[in] value the new value of the bin content + /// \par See also + /// the \ref SetBinContent(const A &... args) "variadic function template overload" accepting arguments directly + template + void SetBinContent(const std::array &indices, const V &value) + { + fEngine.SetBinContent(indices, value); + fStats.Taint(); + } + + /// Set the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RProfile profile({/* two dimensions */}); + /// ROOT::Experimental::RProfile::RProfileBin value = /* ... */; + /// profile.SetBinContent(3, 5, value); + /// \endcode + /// + /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special + /// values. See also the class documentation of RBinIndex. + /// + /// Throws an exception if the number of arguments does not match the axis configuration or the bin is not found. + /// + /// \warning Setting the bin content will taint the global histogram statistics. Attempting to access its values, for + /// example calling GetNEntries(), will throw exceptions. + /// + /// \param[in] args the arguments for each axis and the new value of the bin content + /// \par See also + /// the \ref SetBinContent(const std::array &indices, const V &value) "function overload" accepting + /// `std::array` + template + void SetBinContent(const A &...args) + { + fEngine.SetBinContent(args...); + fStats.Taint(); + } + /// \} /// \name Filling /// \{ @@ -378,6 +471,60 @@ public: } } + /// \} + /// \name Operations + /// \{ + + /// Add all bin contents and statistics of another profile histogram. + /// + /// Throws an exception if the axes configurations are not identical. + /// + /// \param[in] other another profile histogram + void Add(const RProfile &other) + { + fEngine.Add(other.fEngine); + fStats.Add(other.fStats); + } + + /// Add all bin contents and statistics of another profile histogram using atomic instructions. + /// + /// Throws an exception if the axes configurations are not identical. + /// + /// \param[in] other another profile histogram that must not be modified during the operation + void AddAtomic(const RProfile &other) + { + fEngine.AddAtomic(other.fEngine); + fStats.AddAtomic(other.fStats); + } + + /// Clear all bin contents and statistics. + void Clear() + { + fEngine.Clear(); + fStats.Clear(); + } + + /// Clone this profile histogram. + /// + /// Copying all bin contents can be an expensive operation, depending on the number of bins. + /// + /// \return the cloned object + RProfile Clone() const + { + RProfile profile(fEngine.Clone()); + profile.fStats = fStats; + return profile; + } + + /// Scale all bin contents and statistics. + /// + /// \param[in] factor the scale factor + void Scale(double factor) + { + fEngine.Scale(factor); + fStats.Scale(factor); + } + /// \} /// %ROOT Streamer function to throw when trying to store an object of this class. diff --git a/hist/histv7/test/hist_profile.cxx b/hist/histv7/test/hist_profile.cxx index a38c0b64ab703..6f9ea26dd6da0 100644 --- a/hist/histv7/test/hist_profile.cxx +++ b/hist/histv7/test/hist_profile.cxx @@ -56,6 +56,185 @@ TEST(RProfile, Constructor) EXPECT_EQ(profile.GetNDimensions(), 3); } +TEST(RProfile, GetFullMultiDimRange) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + + std::mt19937 gen; + std::uniform_real_distribution distX(0, Bins); + std::uniform_real_distribution distV(Bins, 2 * Bins); + std::uniform_real_distribution distW(0.8, 1.2); + static constexpr std::uint64_t Entries = 1000; + for (std::uint64_t i = 0; i < Entries; i++) { + profile.Fill(distX(gen), distV(gen), RWeight(distW(gen))); + } + ASSERT_EQ(profile.GetNEntries(), Entries); + + auto range = profile.GetFullMultiDimRange(); + EXPECT_EQ(std::distance(range.begin(), range.end()), Bins + 2); + + double sumW = 0; + double sumValues = 0; + for (auto &&indices : range) { + auto &bin = profile.GetBinContent(indices); + sumW += bin.fSum; + sumValues += bin.fSumValues; + } + // Numerical differences with EXPECT_DOUBLE_EQ + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW(), sumW); + EXPECT_FLOAT_EQ(profile.GetStats().GetDimensionStats(1).fSumWX, sumValues); +} + +TEST(RProfile, SetBinContent) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + + RProfile::RProfileBin bin; + bin.fSumValues = 42; + + { + RProfile profile(axis); + ASSERT_FALSE(profile.GetStats().IsTainted()); + profile.SetBinContent(RBinIndex(1), bin); + EXPECT_EQ(profile.GetBinContent(RBinIndex(1)).fSumValues, 42); + EXPECT_TRUE(profile.GetStats().IsTainted()); + EXPECT_THROW(profile.GetNEntries(), std::logic_error); + } + + { + RProfile profile(axis); + ASSERT_FALSE(profile.GetStats().IsTainted()); + const std::array indices = {2}; + profile.SetBinContent(indices, bin); + EXPECT_EQ(profile.GetBinContent(indices).fSumValues, 42); + EXPECT_TRUE(profile.GetStats().IsTainted()); + EXPECT_THROW(profile.GetNEntries(), std::logic_error); + } +} + +TEST(RProfile, Add) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profileA(axis); + RProfile profileB(axis); + + profileA.Fill(8.5, 23.0); + profileB.Fill(9.5, 25.0); + + profileA.Add(profileB); + + EXPECT_EQ(profileA.GetNEntries(), 2); + EXPECT_EQ(profileA.GetBinContent(RBinIndex(8)).fSumValues, 23.0); + EXPECT_EQ(profileA.GetBinContent(RBinIndex(9)).fSumValues, 25.0); +} + +TEST(RProfile, AddAtomic) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RProfile profileA(axis); + RProfile profileB(axis); + + profileA.Fill(8.5, 23.0); + profileB.Fill(9.5, 25.0); + + profileA.AddAtomic(profileB); + + EXPECT_EQ(profileA.GetNEntries(), 2); + EXPECT_EQ(profileA.GetBinContent(RBinIndex(8)).fSumValues, 23.0); + EXPECT_EQ(profileA.GetBinContent(RBinIndex(9)).fSumValues, 25.0); +} + +TEST(RProfile, StressAddAtomic) +{ + static constexpr std::size_t NThreads = 4; + static constexpr std::size_t NAddsPerThread = 10000; + static constexpr std::size_t NAdds = NThreads * NAddsPerThread; + + // Fill a single bin, to maximize contention. + const RRegularAxis axis(1, {0, 1}); + RProfile profileA(axis); + RProfile profileB(axis); + profileB.Fill(0.5, 23.0); + + StressInParallel(NThreads, [&] { + for (std::size_t i = 0; i < NAddsPerThread; i++) { + profileA.AddAtomic(profileB); + } + }); + + EXPECT_EQ(profileA.GetNEntries(), NAdds); + EXPECT_EQ(profileA.GetBinContent(0).fSumValues, 23.0 * NAdds); +} + +TEST(RProfile, AddExceptionSafety) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis regularAxis(Bins, {0, Bins}); + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis categoricalAxis(categories); + + RProfile profileA({regularAxis, regularAxis}); + RProfile profileB({regularAxis, categoricalAxis}); + + profileA.Fill(1.5, 2.5, 23.0); + ASSERT_EQ(profileA.GetNEntries(), 1); + ASSERT_EQ(profileA.GetBinContent(RBinIndex(1), RBinIndex(2)).fSumValues, 23.0); + profileB.Fill(1.5, "b", 25.0); + + EXPECT_THROW(profileA.Add(profileB), std::invalid_argument); + EXPECT_THROW(profileA.AddAtomic(profileB), std::invalid_argument); + + // Verify exception safety. Only the original entry should be there. + EXPECT_EQ(profileA.GetNEntries(), 1); + EXPECT_EQ(profileA.GetBinContent(RBinIndex(1), RBinIndex(2)).fSumValues, 23.0); + EXPECT_EQ(profileA.GetStats().GetSumW(), 1); + EXPECT_EQ(profileA.GetStats().GetSumW2(), 1); + EXPECT_EQ(profileA.GetStats().GetDimensionStats(0).fSumWX, 1.5); + EXPECT_EQ(profileA.GetStats().GetDimensionStats(1).fSumWX, 2.5); +} + +TEST(RProfile, Clear) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + + profile.Fill(8.5, 23.0); + profile.Fill(9.5, 25.0); + + profile.Clear(); + + EXPECT_EQ(profile.GetNEntries(), 0); + EXPECT_EQ(profile.GetBinContent(RBinIndex(8)).fSumValues, 0); + EXPECT_EQ(profile.GetBinContent(RBinIndex(9)).fSumValues, 0); +} + +TEST(RProfile, Clone) +{ + static constexpr std::size_t Bins = 20; + RProfile profileA(Bins, {0, Bins}); + + profileA.Fill(8.5, 23.0); + + RProfile profileB = profileA.Clone(); + ASSERT_EQ(profileB.GetNDimensions(), 1); + ASSERT_EQ(profileB.GetTotalNBins(), Bins + 2); + + EXPECT_EQ(profileB.GetNEntries(), 1); + EXPECT_EQ(profileB.GetBinContent(8).fSumValues, 23.0); + + // Check that we can continue filling the clone. + profileB.Fill(9.5, 25.0); + + EXPECT_EQ(profileA.GetNEntries(), 1); + EXPECT_EQ(profileB.GetNEntries(), 2); + EXPECT_EQ(profileA.GetBinContent(9).fSumValues, 0); + EXPECT_EQ(profileB.GetBinContent(9).fSumValues, 25.0); +} + TEST(RProfile, Fill) { static constexpr std::size_t Bins = 20; @@ -272,3 +451,26 @@ TEST(RProfile, FillForward) ASSERT_FALSE(CopyArgument::HasBeenCopied()); } + +TEST(RProfile, Scale) +{ + static constexpr std::size_t Bins = 20; + RProfile profile(Bins, {0, Bins}); + + profile.Fill(8.5, 23.0, RWeight(0.8)); + profile.Fill(9.5, 25.0, RWeight(0.9)); + + static constexpr double Factor = 0.8; + profile.Scale(Factor); + + EXPECT_FLOAT_EQ(profile.GetBinContent(8).fSumValues, Factor * 0.8 * 23.0); + EXPECT_FLOAT_EQ(profile.GetBinContent(9).fSumValues, Factor * 0.9 * 25.0); + + EXPECT_EQ(profile.GetNEntries(), 2); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW(), Factor * 1.7); + EXPECT_FLOAT_EQ(profile.GetStats().GetSumW2(), Factor * Factor * 1.45); + // Cross-checked with TH1 - unchanged compared to FillWeight because the factor cancels out. + EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 1.9931034); + EXPECT_FLOAT_EQ(profile.ComputeMean(), 9.0294118); + EXPECT_FLOAT_EQ(profile.ComputeStdDev(), 0.49913420); +}