[hist] Implement basic operations on RProfile

Add, AddAtomic, Clear, Clone, Scale; implemented similarly to RHist,
forwarding to RHistEngine and RHistStats.

Author: hahnjo

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:
@@ -271,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<ROOT::Experimental::RBinIndex, 2> 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 <std::size_t N, typename V>
+   void SetBinContent(const std::array<RBinIndex, N> &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<RBinIndex, N> &indices, const V &value) "function overload" accepting
+   /// `std::array`
+   template <typename... A>
+   void SetBinContent(const A &...args)
+   {
+      fEngine.SetBinContent(args...);
+      fStats.Taint();
+   }
+
    /// \}
    /// \name Filling
    /// \{
@@ -379,6 +469,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.

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<double> distX(0, Bins);
+   std::uniform_real_distribution<double> distV(Bins, 2 * Bins);
+   std::uniform_real_distribution<double> 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<RBinIndex, 1> 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<std::string> 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;
@@ -225,3 +404,26 @@ TEST(RProfile, FillExceptionSafety)
    EXPECT_EQ(profile.GetStats().GetDimensionStats(1).fSumWX, 2.5);
    EXPECT_EQ(profile.GetStats().GetDimensionStats(2).fSumWX, 3.5);
 }
+
+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);
+}