[RDF] Adding support for HistoNSparseD in RDF

THnSparseD model is added; Interfaces for HistoNSparseD are added
Adding tests for HistoNSparseD including concurrency

Author: pjacka

README/ReleaseNotes/v638/index.md

@@ -32,6 +32,7 @@ The following people have contributed to this new version:
  Florian Uhlig, GSI,\
  Devajith Valaparambil Sreeramaswamy, CERN/EP-SFT,\
  Vassil Vassilev, Princeton
+ Petr Jacka, Czech Technical University in Prague
 
 ## Deprecation and Removal
 
@@ -112,6 +113,7 @@ If you want to keep using `TList*` return values, you can write a small adapter
   to numbers such as 8 would share one 3-d histogram among 8 threads, greatly reducing the memory consumption. This might slow down execution if the histograms
   are filled at very high rates. Use lower number in this case.
 - The Snapshot method has been refactored so that it does not need anymore compile-time information (i.e. either template arguments or JIT-ting) to know the input column types. This means that any Snapshot call that specifies the template arguments, e.g. `Snapshot<int, float>(..., {"intCol", "floatCol"})` is now redundant and the template arguments can safely be removed from the call. At the same time, Snapshot does not need to JIT compile the column types, practically giving huge speedups depending on the number of columns that need to be written to disk. In certain cases (e.g. when writing O(10000) columns) the speedup can be larger than an order of magnitude. The Snapshot template is now deprecated and it will issue a compile-time warning when called. The function overload is scheduled for removal in ROOT 6.40.
+- Add HistoNSparseD action that fills a sparse N-dimensional histogram.
 
 ## Python Interface
 

bindings/distrdf/python/DistRDF/Operation.py

@@ -51,7 +51,7 @@ def __init__(self, name: str, *args, **kwargs):
             # positional argument. In all Histo*D overload where it is present,
             # it is always the first argument.
             if not isinstance(self.args[0],
-                              (tuple, ROOT.RDF.TH1DModel, ROOT.RDF.TH2DModel, ROOT.RDF.TH3DModel, ROOT.RDF.THnDModel)):
+                              (tuple, ROOT.RDF.TH1DModel, ROOT.RDF.TH2DModel, ROOT.RDF.TH3DModel, ROOT.RDF.THnDModel, ROOT.RDF.THnSparseDModel)):
                 message = (
                     "Creating a histogram without a model is not supported in distributed mode. Please make sure to "
                     "specify the histogram model when rerunning the distributed RDataFrame application. For example:\n\n"
@@ -104,6 +104,7 @@ class Transformation(Operation):
     "Histo2D": Histo,
     "Histo3D": Histo,
     "HistoND": Histo,
+    "HistoNSparseD": Histo,
     "Max": Action,
     "Mean": Action,
     "Min": Action,

bindings/distrdf/test/test_operation.py

@@ -112,3 +112,14 @@ def test_histond_without_model(self):
         """Creating a histogram without model raises ValueError."""
         with self.assertRaises(ValueError):
             _ = Operation.create_op("HistoND", ["a", "b", "c", "d"])
+    
+    def test_histonsparsed_with_thnsparsedmodel(self):
+        """THnDModel"""
+        op = Operation.create_op("HistoNSparseD", ROOT.RDF.THnSparseDModel(), ["a", "b", "c", "d"])
+        self.assertIsInstance(op, Operation.Histo)
+        self.assertEqual(op.name, "HistoNSparseD")
+
+    def test_histonsparsed_without_model(self):
+        """Creating a histogram without model raises ValueError."""
+        with self.assertRaises(ValueError):
+            _ = Operation.create_op("HistoNSparseD", ["a", "b", "c", "d"])

roottest/python/distrdf/backends/check_reducer_merge.py

@@ -125,6 +125,28 @@ def test_histond_merge(self, payload):
         assert histond_distrdf.GetEntries() == histond_rdf.GetEntries()
         assert histond_distrdf.GetNbins() == histond_rdf.GetNbins()
 
+    def test_histonsparsed_merge(self, payload):
+        """Check the working of HistoND merge operation in the reducer."""
+        nbins = (10, 10, 10, 10)
+        xmin = (0., 0., 0., 0.)
+        xmax = (100., 100., 100., 100.)
+        modelTHNSparseD = ("name", "title", 4, nbins, xmin, xmax)
+        colnames = ("x0", "x1", "x2", "x3")
+
+        connection, _ = payload
+        distrdf = ROOT.RDataFrame(100, executor=connection)
+
+        rdf = ROOT.RDataFrame(100)
+
+        distrdf_withcols = self.define_four_columns(distrdf, colnames)
+        rdf_withcols = self.define_four_columns(rdf, colnames)
+
+        histond_distrdf = distrdf_withcols.HistoNSparseD(modelTHNSparseD, colnames)
+        histond_rdf = rdf_withcols.HistoNSparseD(modelTHNSparseD, colnames)
+
+        assert histond_distrdf.GetEntries() == histond_rdf.GetEntries()
+        assert histond_distrdf.GetNbins() == histond_rdf.GetNbins()
+
     def test_profile1d_merge(self, payload):
         """Check the working of Profile1D merge operation in the reducer."""
         # Operations with DistRDF

tree/dataframe/inc/ROOT/RDF/HistoModels.hxx

@@ -20,6 +20,10 @@ class TH3D;
 template <typename T>
 class THnT;
 using THnD = THnT<double>;
+template <typename T>
+class THnSparseT;
+class TArrayD;
+using THnSparseD = THnSparseT<TArrayD>;
 class TProfile;
 class TProfile2D;
 
@@ -123,6 +127,33 @@ struct THnDModel {
    std::shared_ptr<::THnD> GetHistogram() const;
 };
 
+struct THnSparseDModel {
+   TString fName;
+   TString fTitle;
+   int fDim;
+   std::vector<int> fNbins;
+   std::vector<double> fXmin;
+   std::vector<double> fXmax;
+   std::vector<std::vector<double>> fBinEdges;
+   Int_t fChunkSize;
+
+   THnSparseDModel() = default;
+   THnSparseDModel(const THnSparseDModel &) = default;
+   ~THnSparseDModel();
+   THnSparseDModel(const ::THnSparseD &h);
+   THnSparseDModel(const char *name, const char *title, int dim, const int *nbins, const double *xmin, const double *xmax, Int_t chunksize=1024 *16);
+   // alternate version with std::vector to allow more convenient initialization from PyRoot
+   THnSparseDModel(const char *name, const char *title, int dim, const std::vector<int> &nbins,
+             const std::vector<double> &xmin, const std::vector<double> &xmax, Int_t chunksize=1024 *16);
+   THnSparseDModel(const char *name, const char *title, int dim, const int *nbins,
+             const std::vector<std::vector<double>> &xbins, Int_t chunksize=1024 *16);
+   THnSparseDModel(const char *name, const char *title, int dim, const std::vector<int> &nbins,
+             const std::vector<std::vector<double>> &xbins, Int_t chunksize=1024 *16);
+   std::shared_ptr<::THnSparseD> GetHistogram() const;
+};
+
+
+
 struct TProfile1DModel {
    TString fName;
    TString fTitle;

tree/dataframe/inc/ROOT/RDF/InterfaceUtils.hxx

@@ -89,6 +89,7 @@ struct Histo1D{};
 struct Histo2D{};
 struct Histo3D{};
 struct HistoND{};
+struct HistoNSparseD{};
 struct Graph{};
 struct GraphAsymmErrors{};
 struct Profile1D{};
@@ -121,7 +122,7 @@ struct HistoUtils<T, false> {
    static bool HasAxisLimits(T &) { return true; }
 };
 
-// Generic filling (covers Histo2D, HistoND, Profile1D and Profile2D actions, with and without weights)
+// Generic filling (covers Histo2D, HistoND, HistoNSparseD, Profile1D and Profile2D actions, with and without weights)
 template <typename... ColTypes, typename ActionTag, typename ActionResultType, typename PrevNodeType>
 std::unique_ptr<RActionBase>
 BuildAction(const ColumnNames_t &bl, const std::shared_ptr<ActionResultType> &h, const unsigned int nSlots,

tree/dataframe/inc/ROOT/RDF/RInterface.hxx

@@ -40,6 +40,7 @@
 #include "TH2.h" // For Histo actions
 #include "TH3.h" // For Histo actions
 #include "THn.h"
+#include "THnSparse.h"
 #include "TProfile.h"
 #include "TProfile2D.h"
 #include "TStatistic.h"
@@ -2286,6 +2287,82 @@ public:
                                                                                       columnList.size());
    }
 
+
+   ////////////////////////////////////////////////////////////////////////////
+   /// \brief Fill and return a sparse N-dimensional histogram (*lazy action*).
+   /// \tparam FirstColumn The first type of the column the values of which are used to fill the object. Inferred if not
+   /// present.
+   /// \tparam OtherColumns A list of the other types of the columns the values of which are used to fill the
+   /// object.
+   /// \param[in] model The returned histogram will be constructed using this as a model.
+   /// \param[in] columnList
+   /// A list containing the names of the columns that will be passed when calling `Fill`.
+   ///  (N columns for unweighted filling, or N+1 columns for weighted filling)
+   /// \return the N-dimensional histogram wrapped in a RResultPtr.
+   ///
+   /// This action is *lazy*: upon invocation of this method the calculation is
+   /// booked but not executed. See RResultPtr documentation.
+   ///
+   /// ### Example usage:
+   /// ~~~{.cpp}
+   /// auto myFilledObj = myDf.HistoND<float, float, float, float>({"name","title", 4,
+   ///                                                {40,40,40,40}, {20.,20.,20.,20.}, {60.,60.,60.,60.}},
+   ///                                               {"col0", "col1", "col2", "col3"});
+   /// ~~~
+   ///
+   template <typename FirstColumn, typename... OtherColumns> // need FirstColumn to disambiguate overloads
+   RResultPtr<::THnSparseD> HistoNSparseD(const THnSparseDModel &model, const ColumnNames_t &columnList)
+   {
+      std::shared_ptr<::THnSparseD> h(nullptr);
+      {
+         ROOT::Internal::RDF::RIgnoreErrorLevelRAII iel(kError);
+         h = model.GetHistogram();
+
+         if (int(columnList.size()) == (h->GetNdimensions() + 1)) {
+            h->Sumw2();
+         } else if (int(columnList.size()) != h->GetNdimensions()) {
+            throw std::runtime_error("Wrong number of columns for the specified number of histogram axes.");
+         }
+      }
+      return CreateAction<RDFInternal::ActionTags::HistoNSparseD, FirstColumn, OtherColumns...>(columnList, h, h,
+                                                                                          fProxiedPtr);
+   }
+
+   ////////////////////////////////////////////////////////////////////////////
+   /// \brief Fill and return a sparse N-dimensional histogram (*lazy action*).
+   /// \param[in] model The returned histogram will be constructed using this as a model.
+   /// \param[in] columnList A list containing the names of the columns that will be passed when calling `Fill`
+   ///  (N columns for unweighted filling, or N+1 columns for weighted filling)
+   /// \return the N-dimensional histogram wrapped in a RResultPtr.
+   ///
+   /// This action is *lazy*: upon invocation of this method the calculation is
+   /// booked but not executed. Also see RResultPtr.
+   ///
+   /// ### Example usage:
+   /// ~~~{.cpp}
+   /// auto myFilledObj = myDf.HistoNSparseD({"name","title", 4,
+   ///                                                {40,40,40,40}, {20.,20.,20.,20.}, {60.,60.,60.,60.}},
+   ///                                               {"col0", "col1", "col2", "col3"});
+   /// ~~~
+   ///
+   RResultPtr<::THnSparseD> HistoNSparseD(const THnSparseDModel &model, const ColumnNames_t &columnList)
+   {
+      std::shared_ptr<::THnSparseD> h(nullptr);
+      {
+         ROOT::Internal::RDF::RIgnoreErrorLevelRAII iel(kError);
+         h = model.GetHistogram();
+
+         if (int(columnList.size()) == (h->GetNdimensions() + 1)) {
+            h->Sumw2();
+         } else if (int(columnList.size()) != h->GetNdimensions()) {
+            throw std::runtime_error("Wrong number of columns for the specified number of histogram axes.");
+         }
+      }
+      return CreateAction<RDFInternal::ActionTags::HistoNSparseD, RDFDetail::RInferredType>(columnList, h, h, fProxiedPtr,
+                                                                                      columnList.size());
+   }
+
+
    ////////////////////////////////////////////////////////////////////////////
    /// \brief Fill and return a TGraph object (*lazy action*).
    /// \tparam X The type of the column used to fill the x axis.

tree/dataframe/inc/ROOT/RDF/RMergeableValue.hxx

@@ -245,6 +245,7 @@ actions:
 - [Histo{1D,2D,3D}]
   (classROOT_1_1RDF_1_1RInterface.html#a247ca3aeb7ce5b95015b7fae72983055)
 - [HistoND](classROOT_1_1RDF_1_1RInterface.html#a0c9956a0f48c26f8e4294e17376c7fea)
+- [HistoNSparseD](classROOT_1_1RDF_1_1RInterface.html)
 - [Profile{1D,2D}]
   (classROOT_1_1RDF_1_1RInterface.html#a8ef7dc16b0e9f7bc9cfbe2d9e5de0cef)
 - [Stats](classROOT_1_1RDF_1_1RInterface.html#abc68922c464e472f5f856e8981955af6)

tree/dataframe/src/RDFHistoModels.cxx

@@ -20,6 +20,7 @@
 #include "TH2.h"
 #include "TH3.h"
 #include "THn.h"
+#include "THnSparse.h"
 
 /**
  * \class ROOT::RDF::TH1DModel
@@ -46,6 +47,12 @@
  * \note It stores only basic settings such as name, title, bins, bin edges,
  * but not others such as fSumw2.
  *
+ * \class ROOT::RDF::THnSparseDModel
+ * \ingroup dataframe
+ * \brief A struct which stores some basic parameters of a THnSparseD
+ * \note It stores only basic settings such as name, title, bins, bin edges,
+ * but not others such as fSumw2.
+ *
  * \class ROOT::RDF::TProfile1DModel
  * \ingroup dataframe
  * \brief A struct which stores some basic parameters of a TProfile
@@ -295,6 +302,78 @@ std::shared_ptr<::THnD> THnDModel::GetHistogram() const
 }
 THnDModel::~THnDModel() {}
 
+
+THnSparseDModel::THnSparseDModel(const ::THnSparseD &h)
+   : fName(h.GetName()), fTitle(h.GetTitle()), fDim(h.GetNdimensions()), fNbins(fDim), fXmin(fDim), fXmax(fDim),
+     fBinEdges(fDim), fChunkSize(h.GetChunkSize())
+{
+   for (int idim = 0; idim < fDim; ++idim) {
+      fNbins[idim] = h.GetAxis(idim)->GetNbins();
+      SetAxisProperties(h.GetAxis(idim), fXmin[idim], fXmax[idim], fBinEdges[idim]);
+   }
+}
+
+THnSparseDModel::THnSparseDModel(const char *name, const char *title, int dim, const int *nbins, const double *xmin,
+                     const double *xmax, Int_t chunksize)
+   : fName(name), fTitle(title), fDim(dim), fBinEdges(dim), fChunkSize(chunksize)
+{
+   fNbins.reserve(fDim);
+   fXmin.reserve(fDim);
+   fXmax.reserve(fDim);
+   for (int idim = 0; idim < fDim; ++idim) {
+      fNbins.push_back(nbins[idim]);
+      fXmin.push_back(xmin[idim]);
+      fXmax.push_back(xmax[idim]);
+   }
+}
+
+THnSparseDModel::THnSparseDModel(const char *name, const char *title, int dim, const std::vector<int> &nbins,
+                     const std::vector<double> &xmin, const std::vector<double> &xmax, Int_t chunksize)
+   : fName(name), fTitle(title), fDim(dim), fNbins(nbins), fXmin(xmin), fXmax(xmax), fBinEdges(dim), fChunkSize(chunksize)
+{
+}
+
+THnSparseDModel::THnSparseDModel(const char *name, const char *title, int dim, const int *nbins,
+                     const std::vector<std::vector<double>> &xbins, Int_t chunksize)
+   : fName(name), fTitle(title), fDim(dim), fXmin(dim, 0.), fXmax(dim, 64.), fBinEdges(xbins), fChunkSize(chunksize)
+{
+   fNbins.reserve(fDim);
+   for (int idim = 0; idim < fDim; ++idim) {
+      fNbins.push_back(nbins[idim]);
+   }
+}
+
+THnSparseDModel::THnSparseDModel(const char *name, const char *title, int dim, const std::vector<int> &nbins,
+                     const std::vector<std::vector<double>> &xbins, Int_t chunksize)
+   : fName(name), fTitle(title), fDim(dim), fNbins(nbins), fXmin(dim, 0.), fXmax(dim, 64.), fBinEdges(xbins), fChunkSize(chunksize)
+{
+}
+
+std::shared_ptr<::THnSparseD> THnSparseDModel::GetHistogram() const
+{
+   bool varbinning = false;
+   for (const auto &bins : fBinEdges) {
+      if (!bins.empty()) {
+         varbinning = true;
+         break;
+      }
+   }
+   std::shared_ptr<::THnSparseD> h;
+   if (varbinning) {
+      std::vector<TAxis> axes(fDim);
+      for (int idim = 0; idim < fDim; ++idim) {
+         axes[idim] = TAxis(fNbins[idim], fBinEdges[idim].data());
+      }
+      h = std::make_shared<::THnSparseD>(fName, fTitle, axes, fChunkSize);
+      // h = std::make_shared<::THnSparseD>(fName, fTitle, fDim, fNbins.data(), fBinEdges, fChunkSize);
+   } else {
+      h = std::make_shared<::THnSparseD>(fName, fTitle, fDim, fNbins.data(), fXmin.data(), fXmax.data(), fChunkSize);
+   }
+   return h;
+}
+THnSparseDModel::~THnSparseDModel() {}
+
+
 // Profiles
 
 TProfile1DModel::TProfile1DModel(const ::TProfile &h)

tree/dataframe/src/RDataFrame.cxx

@@ -135,6 +135,7 @@ produce many different results in one event loop. Instant actions trigger the ev
 | GraphAsymmErrors() | Fills a TGraphAsymmErrors. Should be used for any type of graph with errors, including cases with errors on one of the axes only. If multi-threading is enabled, the order of the points may not be the one expected, it is therefore suggested to sort if before drawing. |
 | Histo1D(), Histo2D(), Histo3D() | Fill a one-, two-, three-dimensional histogram with the processed column values. |
 | HistoND() | Fill an N-dimensional histogram with the processed column values. |
+| HistoNSparseD() | Fill an N-dimensional sparse histogram with the processed column values. Memory is allocated only for non-empty bins. |
 | Max() | Return the maximum of processed column values. If the type of the column is inferred, the return type is `double`, the type of the column otherwise.|
 | Mean() | Return the mean of processed column values.|
 | Min() | Return the minimum of processed column values. If the type of the column is inferred, the return type is `double`, the type of the column otherwise.|
@@ -737,7 +738,7 @@ parts of the RDataFrame API currently work with this package. The subset that is
 - FilterMissing
 - Graph
 - Histo[1,2,3]D
-- HistoND
+- HistoND, HistoNSparseD
 - Max
 - Mean
 - Min