diff --git a/hist/histv7/doc/CodeArchitecture.md b/hist/histv7/doc/CodeArchitecture.md index 15bd648dc7291..0cb28176b5150 100644 --- a/hist/histv7/doc/CodeArchitecture.md +++ b/hist/histv7/doc/CodeArchitecture.md @@ -79,7 +79,7 @@ Each instance has a local `RHistStats` object to avoid contention on the global A single bin index, which is just an integer for normal bins. `Underflow()` and `Overflow()` are special values and not ordered with respect to others. -Objects of this type are passed by value; most notably to `GetBinContent`. +Objects of this type are passed by value; most notably to `GetBinContent` and `SetBinContent`. ### `RBinIndexRange` diff --git a/hist/histv7/doc/DesignImplementation.md b/hist/histv7/doc/DesignImplementation.md index e58b776d89ac2..836abf60c9b90 100644 --- a/hist/histv7/doc/DesignImplementation.md +++ b/hist/histv7/doc/DesignImplementation.md @@ -58,10 +58,12 @@ Special arguments are passed last. For example ```cpp template void Fill(const std::tuple &args, RWeight w); +template void SetBinContent(const std::array &indices, const V &value); ``` -The same works for the variadic function templates that will check the type of the last argument. +Note that we accept mandatory arguments with a template type as well to allow automatic conversion. -For profiles, we accept the value with a template type as well to allow automatic conversion to `double`, for example from `int`. +Variadic function templates receive all arguments in a single function parameter pack. +For optional arguments, the function will check the type of the last argument to determine if it was passed. ## Miscellaneous diff --git a/hist/histv7/inc/ROOT/RHist.hxx b/hist/histv7/inc/ROOT/RHist.hxx index 9054e00b91e6e..cff0c8fc0ca81 100644 --- a/hist/histv7/inc/ROOT/RHist.hxx +++ b/hist/histv7/inc/ROOT/RHist.hxx @@ -238,6 +238,64 @@ public: /// \return the multidimensional range RBinIndexMultiDimRange GetFullMultiDimRange() const { return fEngine.GetFullMultiDimRange(); } + /// Set the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RHist hist({/* two dimensions */}); + /// std::array indices = {3, 5}; + /// int value = /* ... */; + /// hist.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) const "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::RHist hist({/* two dimensions */}); + /// int value = /* ... */; + /// hist.SetBinContent(ROOT::Experimental::RBinIndex(3), ROOT::Experimental::RBinIndex(5), value); + /// // ... or construct the RBinIndex arguments implicitly from integers: + /// hist.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) const "function overload" + /// accepting `std::array` + template + void SetBinContent(const A &...args) + { + fEngine.SetBinContent(args...); + fStats.Taint(); + } + /// Add all bin contents and statistics of another histogram. /// /// Throws an exception if the axes configurations are not identical. diff --git a/hist/histv7/inc/ROOT/RHistEngine.hxx b/hist/histv7/inc/ROOT/RHistEngine.hxx index abd37de3274cc..d2c1b476a177f 100644 --- a/hist/histv7/inc/ROOT/RHistEngine.hxx +++ b/hist/histv7/inc/ROOT/RHistEngine.hxx @@ -29,14 +29,6 @@ class TBuffer; namespace ROOT { namespace Experimental { -// forward declarations for friend declaration -template -class RHistEngine; -namespace Internal { -template -static void SetBinContent(RHistEngine &hist, const std::array &indices, const T &value); -} // namespace Internal - /** A histogram data structure to bin data along multiple dimensions. @@ -73,9 +65,6 @@ class RHistEngine final { template friend class RHistEngine; - template - friend void Internal::SetBinContent(RHistEngine &, const std::array &, const T &); - /// The axis configuration for this histogram. Relevant methods are forwarded from the public interface. Internal::RAxes fAxes; /// The bin contents for this histogram @@ -247,6 +236,78 @@ public: /// \return the multidimensional range RBinIndexMultiDimRange GetFullMultiDimRange() const { return fAxes.GetFullMultiDimRange(); } + /// Set the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RHistEngine hist({/* two dimensions */}); + /// std::array indices = {3, 5}; + /// int value = /* ... */; + /// hist.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. + /// + /// \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) const "variadic function template overload" accepting arguments + /// directly + template + void SetBinContent(const std::array &indices, const V &value) + { + // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might + // be confusing for users. + if (N != GetNDimensions()) { + throw std::invalid_argument("invalid number of indices passed to SetBinContent"); + } + RLinearizedIndex index = fAxes.ComputeGlobalIndex(indices); + if (!index.fValid) { + throw std::invalid_argument("bin not found in SetBinContent"); + } + assert(index.fIndex < fBinContents.size()); + // To allow conversion, we have to accept value with a template type V to capture any argument. Otherwise it would + // select the variadic function template... + fBinContents[index.fIndex] = value; + } + +private: + template + void SetBinContentImpl(const std::tuple &args, std::index_sequence) + { + std::array indices{std::get(args)...}; + SetBinContent(indices, std::get(args)); + } + +public: + /// Set the content of a single bin. + /// + /// \code + /// ROOT::Experimental::RHistEngine hist({/* two dimensions */}); + /// int value = /* ... */; + /// hist.SetBinContent(ROOT::Experimental::RBinIndex(3), ROOT::Experimental::RBinIndex(5), value); + /// // ... or construct the RBinIndex arguments implicitly from integers: + /// hist.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. + /// + /// \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) const "function overload" + /// accepting `std::array` + template + void SetBinContent(const A &...args) + { + auto t = std::forward_as_tuple(args...); + SetBinContentImpl(t, std::make_index_sequence()); + } + /// Add all bin contents of another histogram. /// /// Throws an exception if the axes configurations are not identical. @@ -582,20 +643,6 @@ public: void Streamer(TBuffer &) { throw std::runtime_error("unable to store RHistEngine"); } }; -namespace Internal { -/// %Internal function to set the content of a single bin. -template -static void SetBinContent(RHistEngine &hist, const std::array &indices, const T &value) -{ - RLinearizedIndex index = hist.fAxes.ComputeGlobalIndex(indices); - if (!index.fValid) { - throw std::invalid_argument("bin not found in SetBinContent"); - } - assert(index.fIndex < hist.fBinContents.size()); - hist.fBinContents[index.fIndex] = value; -} -} // namespace Internal - } // namespace Experimental } // namespace ROOT diff --git a/hist/histv7/inc/ROOT/RHistStats.hxx b/hist/histv7/inc/ROOT/RHistStats.hxx index 462de88082aeb..33ee2d6acbdc1 100644 --- a/hist/histv7/inc/ROOT/RHistStats.hxx +++ b/hist/histv7/inc/ROOT/RHistStats.hxx @@ -114,6 +114,16 @@ private: double fSumW2 = 0.0; /// The sums per dimension std::vector fDimensionStats; + /// Whether this object is tainted, in which case read access will throw. This is used if a user modifies bin + /// contents explicitly or slices histograms without preserving all entries, for example. + bool fTainted = false; + + void ThrowIfTainted() const + { + if (fTainted) { + throw std::logic_error("statistics are tainted, for example after setting bin contents or slicing"); + } + } public: /// Construct a statistics object. @@ -129,9 +139,21 @@ public: std::size_t GetNDimensions() const { return fDimensionStats.size(); } - std::uint64_t GetNEntries() const { return fNEntries; } - double GetSumW() const { return fSumW; } - double GetSumW2() const { return fSumW2; } + std::uint64_t GetNEntries() const + { + ThrowIfTainted(); + return fNEntries; + } + double GetSumW() const + { + ThrowIfTainted(); + return fSumW; + } + double GetSumW2() const + { + ThrowIfTainted(); + return fSumW2; + } /// Get the statistics object for one dimension. /// @@ -141,6 +163,8 @@ public: /// \return the statistics object const RDimensionStats &GetDimensionStats(std::size_t dim = 0) const { + ThrowIfTainted(); + const RDimensionStats &stats = fDimensionStats.at(dim); if (!stats.fEnabled) { throw std::invalid_argument("dimension is disabled"); @@ -157,6 +181,13 @@ public: bool IsEnabled(std::size_t dim) const { return fDimensionStats.at(dim).fEnabled; } + /// Taint this statistics object. + /// + /// It can still be filled, but any read access will throw until Clear() is called. + void Taint() { fTainted = true; } + + bool IsTainted() const { return fTainted; } + /// Add all entries from another statistics object. /// /// Throws an exception if the number of dimensions are not identical. @@ -164,6 +195,8 @@ public: /// \param[in] other another statistics object void Add(const RHistStats &other) { + // NB: this method does *not* call ThrowIfTainted() to allow adding RHist which may contain a tainted statistics + // object. if (fDimensionStats.size() != other.fDimensionStats.size()) { throw std::invalid_argument("number of dimensions not identical in Add"); } @@ -178,6 +211,7 @@ public: fDimensionStats[i].Add(other.fDimensionStats[i]); } } + fTainted |= other.fTainted; } /// Add all entries from another statistics object using atomic instructions. @@ -187,6 +221,8 @@ public: /// \param[in] other another statistics object that must not be modified during the operation void AddAtomic(const RHistStats &other) { + // NB: this method does *not* call ThrowIfTainted() to allow adding RHist which may contain a tainted statistics + // object. if (fDimensionStats.size() != other.fDimensionStats.size()) { throw std::invalid_argument("number of dimensions not identical in Add"); } @@ -201,6 +237,7 @@ public: fDimensionStats[i].AddAtomic(other.fDimensionStats[i]); } } + fTainted |= other.fTainted; } /// Clear this statistics object. @@ -212,6 +249,7 @@ public: for (std::size_t i = 0; i < fDimensionStats.size(); i++) { fDimensionStats[i].Clear(); } + fTainted = false; } /// Compute the number of effective entries. @@ -223,6 +261,7 @@ public: /// \return the number of effective entries double ComputeNEffectiveEntries() const { + ThrowIfTainted(); if (fSumW2 == 0) { return std::numeric_limits::signaling_NaN(); } @@ -240,7 +279,7 @@ public: double ComputeMean(std::size_t dim = 0) const { // First get the statistics, which includes checking the argument. - auto &stats = fDimensionStats.at(dim); + auto &stats = GetDimensionStats(dim); if (fSumW == 0) { return std::numeric_limits::signaling_NaN(); } @@ -266,7 +305,7 @@ public: double ComputeVariance(std::size_t dim = 0) const { // First get the statistics, which includes checking the argument. - auto &stats = fDimensionStats.at(dim); + auto &stats = GetDimensionStats(dim); if (fSumW == 0) { return std::numeric_limits::signaling_NaN(); } @@ -310,7 +349,7 @@ public: double ComputeSkewness(std::size_t dim = 0) const { // First get the statistics, which includes checking the argument. - auto &stats = fDimensionStats.at(dim); + auto &stats = GetDimensionStats(dim); if (fSumW == 0) { return std::numeric_limits::signaling_NaN(); } @@ -347,7 +386,7 @@ public: double ComputeKurtosis(std::size_t dim = 0) const { // First get the statistics, which includes checking the argument. - auto &stats = fDimensionStats.at(dim); + auto &stats = GetDimensionStats(dim); if (fSumW == 0) { return std::numeric_limits::signaling_NaN(); } @@ -495,6 +534,8 @@ public: /// \param[in] factor the scale factor void Scale(double factor) { + // NB: this method does *not* call ThrowIfTainted() to allow scaling RHist which may contain a tainted statistics + // object. fSumW *= factor; fSumW2 *= factor * factor; for (std::size_t i = 0; i < fDimensionStats.size(); i++) { diff --git a/hist/histv7/test/hist_engine.cxx b/hist/histv7/test/hist_engine.cxx index bc2ae778f3a7c..a330eb148d910 100644 --- a/hist/histv7/test/hist_engine.cxx +++ b/hist/histv7/test/hist_engine.cxx @@ -132,19 +132,51 @@ TEST(RHistEngine, GetFullMultiDimRange) TEST(RHistEngine, SetBinContent) { - using ROOT::Experimental::Internal::SetBinContent; - static constexpr std::size_t Bins = 20; const RRegularAxis axis(Bins, {0, Bins}); RHistEngine engine({axis}); - std::array indices = {7}; - SetBinContent(engine, indices, 42); - EXPECT_EQ(engine.GetBinContent(indices), 42); + const RBinIndex index(7); + engine.SetBinContent(index, 42); + EXPECT_EQ(engine.GetBinContent(index), 42); + + const std::array indices = {index}; + engine.SetBinContent(indices, 43); + EXPECT_EQ(engine.GetBinContent(indices), 43); + + // This also works if the value must be converted to the bin content type. + RHistEngine engineF({axis}); + engineF.SetBinContent(index, 42); + EXPECT_EQ(engineF.GetBinContent(index), 42); + + engineF.SetBinContent(indices, 43); + EXPECT_EQ(engineF.GetBinContent(indices), 43); +} + +TEST(RHistEngine, SetBinContentInvalidNumberOfArguments) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHistEngine engine1({axis}); + ASSERT_EQ(engine1.GetNDimensions(), 1); + RHistEngine engine2({axis, axis}); + ASSERT_EQ(engine2.GetNDimensions(), 2); + + EXPECT_NO_THROW(engine1.SetBinContent(1, 0)); + EXPECT_THROW(engine1.SetBinContent(1, 2, 0), std::invalid_argument); + + EXPECT_THROW(engine2.SetBinContent(1, 0), std::invalid_argument); + EXPECT_NO_THROW(engine2.SetBinContent(1, 2, 0)); + EXPECT_THROW(engine2.SetBinContent(1, 2, 3, 0), std::invalid_argument); +} + +TEST(RHistEngine, SetBinContentNotFound) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHistEngine engine({axis}); - // "bin not found" - indices = {Bins}; - EXPECT_THROW(SetBinContent(engine, indices, 43), std::invalid_argument); + EXPECT_THROW(engine.SetBinContent(Bins, 0), std::invalid_argument); } TEST(RHistEngine, Add) diff --git a/hist/histv7/test/hist_hist.cxx b/hist/histv7/test/hist_hist.cxx index 2b03b4f7fa4db..90fba91b36aa4 100644 --- a/hist/histv7/test/hist_hist.cxx +++ b/hist/histv7/test/hist_hist.cxx @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -83,6 +84,31 @@ TEST(RHist, GetFullMultiDimRange) EXPECT_FLOAT_EQ(hist.GetStats().GetSumW(), sumW); } +TEST(RHist, SetBinContent) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + + { + RHist hist(axis); + ASSERT_FALSE(hist.GetStats().IsTainted()); + hist.SetBinContent(RBinIndex(1), 42); + EXPECT_EQ(hist.GetBinContent(RBinIndex(1)), 42); + EXPECT_TRUE(hist.GetStats().IsTainted()); + EXPECT_THROW(hist.GetNEntries(), std::logic_error); + } + + { + RHist hist(axis); + ASSERT_FALSE(hist.GetStats().IsTainted()); + const std::array indices = {2}; + hist.SetBinContent(indices, 42); + EXPECT_EQ(hist.GetBinContent(indices), 42); + EXPECT_TRUE(hist.GetStats().IsTainted()); + EXPECT_THROW(hist.GetNEntries(), std::logic_error); + } +} + TEST(RHist, Add) { static constexpr std::size_t Bins = 20; diff --git a/hist/histv7/test/hist_stats.cxx b/hist/histv7/test/hist_stats.cxx index 1b1f9542fb92a..3cfe5e78b716a 100644 --- a/hist/histv7/test/hist_stats.cxx +++ b/hist/histv7/test/hist_stats.cxx @@ -98,6 +98,36 @@ TEST(RHistStats, DisableDimension) stats.Fill(4, 5, 6, RWeight(1)); EXPECT_EQ(stats.GetNEntries(), 4); + EXPECT_THROW(stats.ComputeMean(1), std::invalid_argument); + EXPECT_THROW(stats.ComputeVariance(1), std::invalid_argument); + EXPECT_THROW(stats.ComputeSkewness(1), std::invalid_argument); + EXPECT_THROW(stats.ComputeKurtosis(1), std::invalid_argument); +} + +TEST(RHistStats, Taint) +{ + RHistStats stats(1); + stats.Taint(); + EXPECT_TRUE(stats.IsTainted()); + + // Modifications are still possible. + stats.Add(stats); + stats.Fill(1); + stats.Scale(2.0); + + // Any read access will throw. + EXPECT_THROW(stats.GetNEntries(), std::logic_error); + EXPECT_THROW(stats.GetSumW(), std::logic_error); + EXPECT_THROW(stats.GetSumW2(), std::logic_error); + EXPECT_THROW(stats.GetDimensionStats(), std::logic_error); + EXPECT_THROW(stats.ComputeMean(), std::logic_error); + EXPECT_THROW(stats.ComputeVariance(), std::logic_error); + EXPECT_THROW(stats.ComputeSkewness(), std::logic_error); + EXPECT_THROW(stats.ComputeKurtosis(), std::logic_error); + + // Clear resets the object, including its taint status. + stats.Clear(); + EXPECT_NO_THROW(stats.GetNEntries()); } TEST(RHistStats, Add) diff --git a/hist/histv7/test/hist_user.cxx b/hist/histv7/test/hist_user.cxx index c7791fa770da3..b5d348da0cdd9 100644 --- a/hist/histv7/test/hist_user.cxx +++ b/hist/histv7/test/hist_user.cxx @@ -11,6 +11,12 @@ struct UserWeight { struct User { double fValue = 0; + User &operator=(double value) + { + fValue = value; + return *this; + } + User &operator++() { fValue++; @@ -252,3 +258,19 @@ TEST(RHistEngineUser, Scale) EXPECT_EQ(engine.GetBinContent(8).fValue, Factor * 0.8); EXPECT_EQ(engine.GetBinContent(9).fValue, Factor * 0.9); } + +TEST(RHistEngineUser, SetBinContent) +{ + // Setting uses operator= + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHistEngine engine({axis}); + + const RBinIndex index(7); + engine.SetBinContent(index, 42); + EXPECT_EQ(engine.GetBinContent(index).fValue, 42); + + const std::array indices = {index}; + engine.SetBinContent(indices, 43); + EXPECT_EQ(engine.GetBinContent(indices).fValue, 43); +} diff --git a/hist/histv7util/inc/ROOT/Hist/ConvertToTH1.hxx b/hist/histv7util/inc/ROOT/Hist/ConvertToTH1.hxx index 197f9e62cd1a7..c52de19ad7c30 100644 --- a/hist/histv7util/inc/ROOT/Hist/ConvertToTH1.hxx +++ b/hist/histv7util/inc/ROOT/Hist/ConvertToTH1.hxx @@ -79,6 +79,9 @@ std::unique_ptr ConvertToTH1S(const RHist &hist); /// Convert a one-dimensional histogram to TH1I. /// +/// If the RHistStats are tainted, for example after setting bin contents, the number of entries and the total sum of +/// weights will be unset. +/// /// Throws an exception if the histogram has more than one dimension. /// /// \param[in] hist the RHist to convert diff --git a/hist/histv7util/src/ConvertToTH1.cxx b/hist/histv7util/src/ConvertToTH1.cxx index 86f7ee2cd6e8f..76bad7d8f06bd 100644 --- a/hist/histv7util/src/ConvertToTH1.cxx +++ b/hist/histv7util/src/ConvertToTH1.cxx @@ -69,6 +69,10 @@ std::unique_ptr ConvertToTH1Impl(const RHistEngine &engine) template void ConvertGlobalStatistics(Hist &h, const RHistStats &stats) { + if (stats.IsTainted()) { + return; + } + h.SetEntries(stats.GetNEntries()); Double_t hStats[4] = { diff --git a/hist/histv7util/test/hist_convert_TH1.cxx b/hist/histv7util/test/hist_convert_TH1.cxx index 1ab428ec1e5a0..fc4d87bf4e4f7 100644 --- a/hist/histv7util/test/hist_convert_TH1.cxx +++ b/hist/histv7util/test/hist_convert_TH1.cxx @@ -100,6 +100,27 @@ TEST(ConvertToTH1I, RHist) EXPECT_EQ(stats[3], 2470); } +TEST(ConvertToTH1I, RHistSetBinContentTainted) +{ + static constexpr std::size_t Bins = 20; + RHist hist(Bins, {0, Bins}); + hist.SetBinContent(RBinIndex(1), 42); + ASSERT_TRUE(hist.GetStats().IsTainted()); + + auto th1i = ConvertToTH1I(hist); + ASSERT_TRUE(th1i); + + EXPECT_EQ(th1i->GetBinContent(2), 42); + + EXPECT_EQ(th1i->GetEntries(), 0); + Double_t stats[4]; + th1i->GetStats(stats); + EXPECT_EQ(stats[0], 0); + EXPECT_EQ(stats[1], 0); + EXPECT_EQ(stats[2], 0); + EXPECT_EQ(stats[3], 0); +} + TEST(ConvertToTH1C, RHistEngine) { static constexpr std::size_t Bins = 20; @@ -149,7 +170,7 @@ TEST(ConvertToTH1L, RHistEngine) // Set one 64-bit long long value larger than what double can exactly represent. static constexpr long long Large = (1LL << 60) - 1; const std::array indices = {1}; - ROOT::Experimental::Internal::SetBinContent(engineLL, indices, Large); + engineLL.SetBinContent(indices, Large); th1l = ConvertToTH1L(engineLL); ASSERT_TRUE(th1l);