[hist] Implement RProfileConcurrentFiller

hist/histv7/headers.cmake

@@ -15,6 +15,8 @@ set(histv7_headers
     ROOT/RHistUtils.hxx
     ROOT/RLinearizedIndex.hxx
     ROOT/RProfile.hxx
+    ROOT/RProfileConcurrentFiller.hxx
+    ROOT/RProfileFillContext.hxx
     ROOT/RRegularAxis.hxx
     ROOT/RSliceBinIndexMapper.hxx
     ROOT/RSliceSpec.hxx

hist/histv7/inc/ROOT/RHistEngine.hxx

@@ -78,6 +78,7 @@ class RHistEngine final {
    friend class RHist<BinContentType>;
 
    friend class RProfile;
+   friend class RProfileFillContext;
 
    /// The axis configuration for this histogram. Relevant methods are forwarded from the public interface.
    Internal::RAxes fAxes;
@@ -565,6 +566,27 @@ public:
       }
    }
 
+   /// \}
+   // End the group to ensure that all contained member functions are public.
+
+private:
+   // Also used by RProfileFillContext::Fill(const A &...args) - similar to the variadic
+   // RHistEngine::FillAtomic(const A &...args) below, is has all arguments in the forwarded std::tuple and needs to
+   // explicitly specify how many of them should be used by RAxes::ComputeGlobalIndexImpl<N>(args).
+   template <std::size_t N, typename... A, typename W>
+   void FillAtomicImpl(const std::tuple<A...> &args, const W &weight)
+   {
+      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<N>(args);
+      if (index.fValid) {
+         assert(index.fIndex < fBinContents.size());
+         Internal::AtomicAddRelease(&fBinContents[index.fIndex], weight);
+      }
+   }
+
+public:
+   /// \name Filling
+   /// \{
+
    /// Fill an entry into the histogram with a user-defined weight using atomic instructions.
    ///
    /// This overload is only available for user-defined bin content types.
@@ -587,11 +609,7 @@ public:
       if (sizeof...(A) != GetNDimensions()) {
          throw std::invalid_argument("invalid number of arguments to Fill");
       }
-      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<sizeof...(A)>(args);
-      if (index.fValid) {
-         assert(index.fIndex < fBinContents.size());
-         Internal::AtomicAddRelease(&fBinContents[index.fIndex], weight);
-      }
+      FillAtomicImpl<sizeof...(A)>(args, weight);
    }
 
    /// Fill an entry into the histogram using atomic instructions.

hist/histv7/inc/ROOT/RProfile.hxx

@@ -56,6 +56,8 @@ ROOT::Experimental::RProfile profile(axes);
 Feedback is welcome!
 */
 class RProfile final {
+   friend class RProfileFillContext;
+
    struct RValueWrapper {
       double fValue;
 
@@ -113,6 +115,29 @@ public:
          return *this;
       }
 
+   private:
+      void AtomicAddRelease(double v, double v2, double w, double w2)
+      {
+         Internal::AtomicAddRelease(&fSumValues, v);
+         Internal::AtomicAddRelease(&fSumValues2, v2);
+         Internal::AtomicAddRelease(&fSum, w);
+         Internal::AtomicAddRelease(&fSum2, w2);
+      }
+
+   public:
+      void AtomicAddRelease(const RValueWrapper &rhs)
+      {
+         double v = rhs.fValue;
+         AtomicAddRelease(v, v * v, 1.0, 1.0);
+      }
+
+      void AtomicAddRelease(const RValueWeightWrapper &rhs)
+      {
+         double v = rhs.fValue;
+         double w = rhs.fWeight;
+         AtomicAddRelease(w * v, w * v * v, w, w * w);
+      }
+
       /// Add another bin content using atomic instructions.
       ///
       /// \param[in] rhs another bin content that must not be modified during the operation

hist/histv7/inc/ROOT/RProfileConcurrentFiller.hxx

@@ -0,0 +1,93 @@
+/// \file
+/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+/// Feedback is welcome!
+
+#ifndef ROOT_RProfileConcurrentFiller
+#define ROOT_RProfileConcurrentFiller
+
+#include "RProfile.hxx"
+#include "RProfileFillContext.hxx"
+
+#include <exception>
+#include <memory>
+#include <mutex>
+#include <stdexcept>
+#include <vector>
+
+namespace ROOT {
+namespace Experimental {
+
+/**
+A histogram filler to concurrently fill an RProfile.
+
+\code
+auto profile = std::make_shared<ROOT::Experimental::RProfile>(10, std::make_pair(5, 15));
+{
+   ROOT::Experimental::RProfileConcurrentFiller filler(profile);
+   auto context = filler.CreateFillContext();
+   context->Fill(8.5, 23.0);
+}
+// profile->GetBinContent(ROOT::Experimental::RBinIndex(3)) will return the filled bin content
+\endcode
+
+\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+Feedback is welcome!
+*/
+class RProfileConcurrentFiller final {
+   /// A pointer to the filled profile histogram
+   std::shared_ptr<RProfile> fProfile;
+
+   /// Mutex to protect access to the list of fill contexts (not for filling itself!)
+   std::mutex fMutex;
+   /// The list of fill contexts, for checks during destruction
+   std::vector<std::weak_ptr<RProfileFillContext>> fFillContexts;
+
+public:
+   /// Create a filler object.
+   ///
+   /// \param[in] hist a pointer to the histogram
+   explicit RProfileConcurrentFiller(std::shared_ptr<RProfile> profile) : fProfile(profile)
+   {
+      if (!profile) {
+         throw std::invalid_argument("profile must not be nullptr");
+      }
+   }
+
+   RProfileConcurrentFiller(const RProfileConcurrentFiller &) = delete;
+   RProfileConcurrentFiller(RProfileConcurrentFiller &&) = delete;
+   RProfileConcurrentFiller &operator=(const RProfileConcurrentFiller &) = delete;
+   RProfileConcurrentFiller &operator=(RProfileConcurrentFiller &&) = delete;
+
+   ~RProfileConcurrentFiller()
+   {
+      for (const auto &context : fFillContexts) {
+         if (!context.expired()) {
+            // According to C++ Core Guideline C.36 "A destructor must not fail" and (C.37) "If a destructor tries to
+            // exit with an exception, it’s a bad design error and the program had better terminate".
+            std::terminate(); // GCOVR_EXCL_LINE
+         }
+      }
+   }
+
+   const std::shared_ptr<RProfile> &GetProfile() const { return fProfile; }
+
+   /// Create a new context for concurrent filling.
+   std::shared_ptr<RProfileFillContext> CreateFillContext()
+   {
+      // Cannot use std::make_shared because the constructor of RProfileFillContext is private. Also it would mean that
+      // the (direct) memory of all contexts stays around until the vector of weak_ptr's is cleared.
+      std::shared_ptr<RProfileFillContext> context(new RProfileFillContext(*fProfile));
+
+      {
+         std::lock_guard g(fMutex);
+         fFillContexts.push_back(context);
+      }
+
+      return context;
+   }
+};
+
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/inc/ROOT/RProfileFillContext.hxx

@@ -0,0 +1,147 @@
+/// \file
+/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+/// Feedback is welcome!
+
+#ifndef ROOT_RProfileFillContext
+#define ROOT_RProfileFillContext
+
+#include "RProfile.hxx"
+#include "RHistEngine.hxx"
+#include "RHistStats.hxx"
+#include "RHistUtils.hxx"
+#include "RWeight.hxx"
+
+#include <cstddef>
+#include <tuple>
+
+namespace ROOT {
+namespace Experimental {
+
+/**
+A context to concurrently fill an RProfile.
+
+\sa RProfileConcurrentFiller
+
+\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
+Feedback is welcome!
+*/
+class RProfileFillContext final {
+   friend class RProfileConcurrentFiller;
+
+private:
+   /// A pointer to the filled profile histogram
+   RProfile *fProfile = nullptr;
+
+   /// Local histogram statistics
+   RHistStats fStats;
+
+   /// \sa RProfileConcurrentFiller::CreateFillContent()
+   explicit RProfileFillContext(RProfile &profile) : fProfile(&profile), fStats(profile.GetNDimensions() + 1)
+   {
+      // Propagate disabled dimensions to the local histogram statistics object.
+      const auto &stats = profile.GetStats();
+      for (std::size_t i = 0; i < stats.GetNDimensions(); i++) {
+         if (!stats.IsEnabled(i)) {
+            fStats.DisableDimension(i);
+         }
+      }
+   }
+   RProfileFillContext(const RProfileFillContext &) = delete;
+   RProfileFillContext(RProfileFillContext &&) = default;
+   RProfileFillContext &operator=(const RProfileFillContext &) = delete;
+   RProfileFillContext &operator=(RProfileFillContext &&) = default;
+
+public:
+   ~RProfileFillContext() { Flush(); }
+
+   /// Fill an entry into the profile histogram.
+   ///
+   /// 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
+   /// \sa RProfile::Fill(const std::tuple<A...> &args, const V &value)
+   template <typename... A, typename V>
+   void Fill(const std::tuple<A...> &args, const V &value)
+   {
+      RProfile::RValueWrapper wrapper(value);
+      fProfile->fEngine.FillAtomic(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.
+   ///
+   /// 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
+   /// \sa RProfile::Fill(const std::tuple<A...> &args, const V &value, RWeight weight)
+   template <typename... A, typename V>
+   void Fill(const std::tuple<A...> &args, const V &value, RWeight weight)
+   {
+      RProfile::RValueWeightWrapper wrapper(value, weight.fValue);
+      fProfile->fEngine.FillAtomic(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.
+   ///
+   /// For weighted filling, pass an RWeight as the last argument.
+   ///
+   /// 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
+   /// \sa RProfile::Fill(const A &...args)
+   template <typename... A>
+   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<typename Internal::LastType<A...>::type, RWeight>) {
+            static constexpr std::size_t N = sizeof...(A) - 2;
+            if (N != fProfile->GetNDimensions()) {
+               throw std::invalid_argument("invalid number of arguments to Fill");
+            }
+            RWeight weight = std::get<N + 1>(t);
+            RProfile::RValueWeightWrapper wrapper(std::get<N>(t), weight.fValue);
+            fProfile->fEngine.FillAtomicImpl<N>(t, wrapper);
+         } else {
+            static constexpr std::size_t N = sizeof...(A) - 1;
+            if (N != fProfile->GetNDimensions()) {
+               throw std::invalid_argument("invalid number of arguments to Fill");
+            }
+            RProfile::RValueWrapper wrapper(std::get<N>(t));
+            fProfile->fEngine.FillAtomicImpl<N>(t, wrapper);
+         }
+         fStats.Fill(args...);
+      }
+   }
+
+   /// Flush locally accumulated entries to the profile histogram.
+   void Flush()
+   {
+      fProfile->fStats.AddAtomic(fStats);
+      fStats.Clear();
+   }
+};
+
+} // namespace Experimental
+} // namespace ROOT
+
+#endif

hist/histv7/test/CMakeLists.txt

@@ -9,6 +9,7 @@ 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_profile_concurrent hist_profile_concurrent.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)

hist/histv7/test/hist_profile.cxx

@@ -303,10 +303,12 @@ TEST(RProfile, FillWeight)
    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
+   // Cross-checked with TProfile
    EXPECT_FLOAT_EQ(profile.ComputeNEffectiveEntries(), 1.9931034);
    EXPECT_FLOAT_EQ(profile.ComputeMean(0), 9.0294118);
    EXPECT_FLOAT_EQ(profile.ComputeStdDev(0), 0.49913420);
+   EXPECT_FLOAT_EQ(profile.ComputeMean(1), 24.058824);
+   EXPECT_FLOAT_EQ(profile.ComputeStdDev(1), 0.99826840);
 }
 
 TEST(RProfile, FillWeightInvalidNumberOfArguments)
@@ -408,23 +410,22 @@ TEST(RProfile, FillExceptionSafety)
 TEST(RProfile, FillForward)
 {
    static constexpr std::size_t Bins = 20;
-   const RRegularAxis axis(Bins, {0, Bins});
-   RProfile profile(axis, axis);
+   RProfile profile(Bins, {0, Bins});
    CopyArgument value(23.0);
 
-   std::tuple<CopyArgument, CopyArgument> args(1.5, 2.5);
+   std::tuple<CopyArgument> args(1.5);
    profile.Fill(args, value);
    profile.Fill(args, value, RWeight(0.5));
    EXPECT_EQ(profile.GetNEntries(), 2);
-   EXPECT_EQ(profile.GetBinContent(1, 2).fSumValues, 34.5);
+   EXPECT_EQ(profile.GetBinContent(1).fSumValues, 34.5);
 
    ASSERT_FALSE(CopyArgument::HasBeenCopied());
 
-   CopyArgument arg1(3.5), arg2(4.5);
-   profile.Fill(arg1, arg2, value);
-   profile.Fill(arg1, arg2, value, RWeight(0.5));
+   CopyArgument arg(2.5);
+   profile.Fill(arg, value);
+   profile.Fill(arg, value, RWeight(0.5));
    EXPECT_EQ(profile.GetNEntries(), 4);
-   EXPECT_EQ(profile.GetBinContent(3, 4).fSumValues, 34.5);
+   EXPECT_EQ(profile.GetBinContent(2).fSumValues, 34.5);
 
    ASSERT_FALSE(CopyArgument::HasBeenCopied());
 }