[RF] Better ParamHistFunc::computeBatch with RooAbsBinning::binNumbers

So far, the `ParamHistFunc` BatchMode implementation was still
sub-uptimal because it used the non-vectorized interface of the
RooDataHist.

Using the new `RooAbsBinning::binNumbers()` function too look up
multiple bin indices at once, the implementaiton can be improved and
sped.

In the three-dimensional many-bin case, the new implementation is a bit
more than three times faster than the old one. This should benefit
HistFactory fits with many bins.

A new `testParamHistFunc` unit test was also introduced to validate the
results of a ParamHistFunc both with and without the batch mode,
comparing to manually computed reference results.

Author: guitargeek

roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h

@@ -27,7 +27,6 @@ class ParamHistFunc : public RooAbsReal {
   ParamHistFunc() ;
   ParamHistFunc(const char *name, const char *title, const RooArgList& vars, const RooArgList& paramSet );
   ParamHistFunc(const char *name, const char *title, const RooArgList& vars, const RooArgList& paramSet, const TH1* hist );
-  ~ParamHistFunc() override ;
 
   ParamHistFunc(const ParamHistFunc& other, const char* name = 0);
   TObject* clone(const char* newname) const override { return new ParamHistFunc(*this, newname); }
@@ -67,8 +66,6 @@ class ParamHistFunc : public RooAbsReal {
 
   class CacheElem : public RooAbsCacheElement {
   public:
-    CacheElem()  {} ;
-    ~CacheElem() override {} ;
     RooArgList containedArgs(Action) override {
       RooArgList ret(_funcIntList) ;
       ret.add(_lowIntList);
@@ -85,7 +82,7 @@ class ParamHistFunc : public RooAbsReal {
   RooListProxy _dataVars;             ///< The RooRealVars
   RooListProxy _paramSet ;            ///< interpolation parameters
 
-  Int_t _numBins;
+  Int_t _numBins = 0;
   struct NumBins {
     NumBins() {}
     NumBins(int nx, int ny, int nz) : x{nx}, y{ny}, z{nz}, xy{x*y}, xz{x*z}, yz{y*z}, xyz{xy*z} {}

roofit/histfactory/src/ParamHistFunc.cxx

@@ -60,7 +60,7 @@ ClassImp(ParamHistFunc);
 ////////////////////////////////////////////////////////////////////////////////
 
 ParamHistFunc::ParamHistFunc()
-  : _normIntMgr(this), _numBins(0)
+  : _normIntMgr(this)
 {
   _dataSet.removeSelfFromDir(); // files must not delete _dataSet.
 }
@@ -84,7 +84,6 @@ ParamHistFunc::ParamHistFunc(const char* name, const char* title,
   _normIntMgr(this),
   _dataVars("!dataVars","data Vars",       this),
   _paramSet("!paramSet","bin parameters",  this),
-  _numBins(0),
   _dataSet( (std::string(name)+"_dataSet").c_str(), "", vars)
 {
 
@@ -130,7 +129,6 @@ ParamHistFunc::ParamHistFunc(const char* name, const char* title,
   //  _dataVar("!dataVar","data Var", this, (RooRealVar&) var),
   _dataVars("!dataVars","data Vars",       this),
   _paramSet("!paramSet","bin parameters",  this),
-  _numBins(0),
   _dataSet( (std::string(name)+"_dataSet").c_str(), "", vars, Hist)
 {
 
@@ -187,14 +185,6 @@ ParamHistFunc::ParamHistFunc(const ParamHistFunc& other, const char* name) :
 }
 
 
-////////////////////////////////////////////////////////////////////////////////
-
-ParamHistFunc::~ParamHistFunc()
-{
-  ;
-}
-
-
 ////////////////////////////////////////////////////////////////////////////////
 /// Get the parameter associated with the index.
 /// The index follows RooDataHist indexing conventions.
@@ -603,35 +593,34 @@ double ParamHistFunc::evaluate() const
 /// \param[in,out] evalData Input/output data for evaluating the ParamHistFunc.
 /// \param[in] normSet Normalisation set passed on to objects that are serving values to us.
 void ParamHistFunc::computeBatch(cudaStream_t*, double* output, size_t size, RooFit::Detail::DataMap const& dataMap) const {
-  std::vector<double> oldValues;
-  std::vector<RooSpan<const double>> data;
-  oldValues.reserve(_dataVars.size());
-  data.reserve(_dataVars.size());
-
-  // Retrieve data for all variables
-  for (auto arg : _dataVars) {
-    const auto* var = static_cast<RooRealVar*>(arg);
-    oldValues.push_back(var->getVal());
-    data.push_back(dataMap.at(var));
-  }
 
-  // Run computation for each entry in the dataset
-  for (std::size_t i = 0; i < size; ++i) {
-    for (unsigned int j = 0; j < _dataVars.size(); ++j) {
-      assert(i < data[j].size());
-      auto& var = static_cast<RooRealVar&>(_dataVars[j]);
-      var.setCachedValue(data[j][i], /*notifyClients=*/false);
-    }
+  auto const& n = _numBinsPerDim;
+  // check if _numBins needs to be filled
+  if(n.x == 0) {
+    _numBinsPerDim = getNumBinsPerDim(_dataVars);
+  }
 
-    const auto index = _dataSet.getIndex(_dataVars, /*fast=*/true);
-    const RooAbsReal& param = getParameter(index);
-    output[i] = param.getVal();
+  // Different from the evaluate() funnction that first retrieves the indices
+  // corresponding to the RooDataHist and then transforms them, we can use the
+  // right bin multiplicators to begin with.
+  std::array<int, 3> idxMult{{1, n.x, n.xy}};
+
+  // As a working buffer for the bin indices, we use the tail of the output
+  // buffer. We can't use the same starting pointer, otherwise we would
+  // overwrite the later bin indices as we fill the output.
+  auto indexBuffer = reinterpret_cast<int*>(output + size) - size;
+  std::fill(indexBuffer, indexBuffer + size, 0); // output buffer for bin indices needs to be zero-initialized
+
+  // Use the vectorized RooAbsBinning::binNumbers() to update the total bin
+  // index for each dimension, using the `coef` parameter to multiply with the
+  // right index multiplication factor for each dimension.
+  for (std::size_t iVar = 0; iVar < _dataVars.size(); ++iVar) {
+    _dataSet.getBinnings()[iVar]->binNumbers(dataMap.at(&_dataVars[iVar]).data(), indexBuffer, size, idxMult[iVar]);
   }
 
-  // Restore old values
-  for (unsigned int j = 0; j < _dataVars.size(); ++j) {
-    auto& var = static_cast<RooRealVar&>(_dataVars[j]);
-    var.setCachedValue(oldValues[j], /*notifyClients=*/false);
+  // Finally, look up the parameters and get their values to fill the output buffer
+  for (std::size_t i = 0; i < size; ++i) {
+    output[i] = static_cast<RooAbsReal const&>(_paramSet[indexBuffer[i]]).getVal();
   }
 }
 

roofit/histfactory/test/CMakeLists.txt

@@ -7,5 +7,7 @@
 # @author Stephan Hageboeck CERN, 2019
 
 ROOT_ADD_GTEST(testHistFactory testHistFactory.cxx
-  LIBRARIES RooFitCommon RooFitCore RooFit RooStats HistFactory RooBatchCompute
+  LIBRARIES RooFitCore RooFit RooStats HistFactory
   COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/ref_6.16_example_UsingC_channel1_meas_model.root ${CMAKE_CURRENT_SOURCE_DIR}/ref_6.16_example_UsingC_combined_meas_model.root)
+
+ROOT_ADD_GTEST(testParamHistFunc testParamHistFunc.cxx LIBRARIES RooFitCore HistFactory)

roofit/histfactory/test/testParamHistFunc.cxx

@@ -0,0 +1,76 @@
+// Tests for the ParamHistFunc
+// Authors: Jonas Rembser, CERN  08/2022
+
+#include <RooArgSet.h>
+#include <RooArgSet.h>
+#include <RooConstVar.h>
+#include <RooDataSet.h>
+#include <RooRandom.h>
+#include <RooRealVar.h>
+#include <RooStats/HistFactory/ParamHistFunc.h>
+
+#include <gtest/gtest.h>
+
+/// Validate the ParamHistFunc in the n-dimensional case, comparing both the
+/// BatchMode and the old implementation results to a manually-compute
+/// reference result.
+TEST(ParamHistFunc, ValidateND)
+{
+
+   // Define the number of bins in each dimension
+   std::array<std::size_t, 3> nbins{{11, 32, 8}};
+   std::size_t nbinstot = nbins[0] * nbins[1] * nbins[2];
+
+   // The bin mltiplication factors to look up the right parameters
+   std::array<std::size_t, 3> binMult{{1, nbins[0], nbins[0] * nbins[1]}};
+
+   // Create the variables and set their bin numbers. The range is tweaked
+   // such that each integer values falls in a different bin, starting from
+   // zero.
+   RooRealVar x{"x", "x", -0.5, nbins[0] - 0.5};
+   RooRealVar y{"y", "y", -0.5, nbins[1] - 0.5};
+   RooRealVar z{"z", "z", -0.5, nbins[2] - 0.5};
+   RooArgSet vars{x, y, z};
+   for (std::size_t i = 0; i < nbins.size(); ++i) {
+      static_cast<RooRealVar &>(*vars[i]).setBins(nbins[i]);
+   }
+
+   // Simple set of parameters that just return their index in the parameter
+   // list
+   RooArgList params;
+   for (std::size_t i = 0; i < nbinstot; ++i) {
+      params.add(RooFit::RooConst(i));
+   }
+
+   ParamHistFunc paramHistFunc{"phf", "phf", vars, params};
+
+   std::size_t nEntries = 100;
+
+   RooDataSet data{"data", "data", vars};
+   std::vector<double> resultsRef(nEntries);
+   std::vector<double> resultsScalar(nEntries);
+
+   // Do some things in one go:
+   //   * assing random integer values to each variable in each iteration
+   //   * fill the dataset used for batched evaluation
+   //   * compute the reference result manually
+   //   * compute the result with the ParamHistFunc without BatchMode
+   for (std::size_t i = 0; i < nEntries; ++i) {
+      for (std::size_t iVar = 0; iVar < vars.size(); ++iVar) {
+         auto var = static_cast<RooRealVar *>(vars[iVar]);
+         var->setVal(int(RooRandom::uniform() * nbins[iVar]));
+      }
+      data.add(vars);
+      resultsRef[i] = binMult[0] * x.getVal() + binMult[1] * y.getVal() + binMult[2] * z.getVal();
+      resultsScalar[i] = paramHistFunc.getVal();
+   }
+
+   // Get the results in BatchMode using the dataset
+   auto resultsBatch = paramHistFunc.getValues(data);
+
+   // Validate the results
+   for (std::size_t i = 0; i < nEntries; ++i) {
+      EXPECT_EQ(int(resultsScalar[i]), int(resultsRef[i])) << "Scalar result is not correct!";
+      EXPECT_EQ(int(resultsBatch[i]), int(resultsRef[i])) << "BatchMode result is not correct!";
+   }
+}