[RF] Vectorize weight evaluation with RooDataHist::weights()

This PR improves the speed for evaluating weights in `RooHistPdf`
and `RooHistFunc` for one dimensional histograms with no interpolation.
In the future, `RooDataHist::weights()` can be extended to cover
cases with higher dimensions and interpolation.

The function `RooDataHist::weights()` was implemented to enable
vectorized evaluations of bin weights. In `RooHistPdf` it is implemented
using the new function `RooHistPdf::computeBatch()`, which calls
`RooDataHist::weights()` in the case of no interpolation and 1D
histograms, and `RooAbsReal::computeBatch()` otherwise. In
`RooHistFunc::computeBatch`, `RooDataHist::weights()` is called in the
case of no interpolation and 1D histograms and is unchanged in the other cases.

To calculate the weight, bin indices are stored as a vector
using `RooAbsBinning::binNumbers`, which was implemented in root-project#11151.

Author: HannaOlvhammar

roofit/roofitcore/inc/RooDataHist.h

@@ -106,6 +106,7 @@ class RooDataHist : public RooAbsData, public RooDirItem {
                const std::map<const RooAbsArg*, std::pair<double, double> >& ranges,
                std::function<double(int)> getBinScale = [](int){ return 1.0; } );
 
+  void weights(double* output, RooSpan<double const> xVals, bool correctForBinSize);
   /// Return weight of i-th bin. \see getIndex()
   double weight(std::size_t i) const { return _wgt[i]; }
   double weightFast(const RooArgSet& bin, int intOrder, bool correctForBinSize, bool cdfBoundaries);

roofit/roofitcore/inc/RooHistPdf.h

@@ -102,6 +102,8 @@ class RooHistPdf : public RooAbsPdf {
   std::list<double>* binBoundaries(RooAbsRealLValue& /*obs*/, double /*xlo*/, double /*xhi*/) const override ;
   bool isBinnedDistribution(const RooArgSet&) const override { return _intOrder==0 ; }
 
+  void computeBatch(cudaStream_t*, double* output, size_t size, RooFit::Detail::DataMap const&) const override;
+
 
 protected:
 

roofit/roofitcore/src/RooDataHist.cxx

@@ -1142,6 +1142,33 @@ RooPlot *RooDataHist::plotOn(RooPlot *frame, PlotOpt o) const
   return RooAbsData::plotOn(frame,o) ;
 }
 
+
+////////////////////////////////////////////////////////////////////////////////
+/// A vectorized version of RooDataHist::weight for one dimensional histograms 
+/// with no interpolation. 
+/// \param[out] output An array of weights corresponding the values in xVals.
+/// \param[in] xVals An array of coordinates for which the weights should be
+///                  calculated.
+/// \param[in] correctForBinSize Enable the inverse bin volume correction factor.
+
+void RooDataHist::weights(double* output, RooSpan<double const> xVals, bool correctForBinSize) 
+{
+  auto const nEvents = xVals.size();
+  RooAbsBinning const& binning = *_lvbins[0];
+
+  // Reuse the output buffer for bin indices and zero-initialize it
+  auto binIndices = reinterpret_cast<int*>(output + nEvents) - nEvents; 
+  std::fill(binIndices, binIndices + nEvents, 0);
+
+  binning.binNumbers(xVals.data(), binIndices, nEvents);
+
+  for (std::size_t i=0; i < nEvents; ++i) {
+    auto binIdx = binIndices[i];
+    output[i] = correctForBinSize ? _wgt[binIdx] / _binv[binIdx] : _wgt[binIdx];
+  }
+}
+
+
 ////////////////////////////////////////////////////////////////////////////////
 /// A faster version of RooDataHist::weight that assumes the passed arguments
 /// are aligned with the histogram variables.
@@ -1243,7 +1270,6 @@ double RooDataHist::weightInterpolated(const RooArgSet& bin, int intOrder, bool
 
   double wInt{0} ;
   if (varInfo.nRealVars == 1) {
-
     // buffer needs to be 2 x (interpolation order + 1), with the factor 2 for x and y.
     _interpolationBuffer.resize(2 * intOrder + 2);
 

roofit/roofitcore/src/RooHistFunc.cxx

@@ -196,6 +196,12 @@ double RooHistFunc::evaluate() const
 
 
 void RooHistFunc::computeBatch(cudaStream_t*, double* output, size_t size, RooFit::Detail::DataMap const& dataMap) const {
+  if (_depList.size() == 1 && _intOrder == 0) {
+    auto xVals = dataMap.at(_depList[0]);
+    _dataHist->weights(output, xVals, false);
+    return;
+  }
+  
   std::vector<RooSpan<const double>> inputValues;
   for (const auto& obs : _depList) {
     auto realObs = dynamic_cast<const RooAbsReal*>(obs);

roofit/roofitcore/src/RooHistPdf.cxx

@@ -188,7 +188,17 @@ RooHistPdf::~RooHistPdf()
 }
 
 
+void RooHistPdf::computeBatch(cudaStream_t*, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const {
 
+  // For interpolation and histograms of higher dimension, use base function
+  if(_pdfObsList.size() > 1 || _intOrder > 0) {
+      RooAbsReal::computeBatch(nullptr, output, nEvents, dataMap);
+      return;
+  }
+
+  auto xVals = dataMap.at(_pdfObsList[0]);
+  _dataHist->weights(output, xVals, true);
+}
 
 
 ////////////////////////////////////////////////////////////////////////////////