Vectorized PDFs

math/mathmore/inc/Math/SpecFuncMathMore.h

@@ -121,7 +121,9 @@ namespace Math {
   double assoc_legendre(unsigned l, unsigned m, double x);
 
 
-
+  namespace internal{
+    double legendre(unsigned l, unsigned m, double x);
+  }
 
 
   /**

math/mathmore/src/SpecFuncMathMore.cxx

@@ -55,8 +55,11 @@ double assoc_legendre(unsigned l, unsigned m, double x) {
 
 }
 
-
-
+namespace internal{
+  double legendre(unsigned l, unsigned m, double x) {
+    return gsl_sf_legendre_Plm(l, m, x);
+  }
+}
 
 
 // [5.2.1.4] (complete) elliptic integral of the first kind

roofit/roofit/inc/RooArgusBG.h

@@ -43,6 +43,8 @@ class RooArgusBG : public RooAbsPdf {
   RooRealProxy p ;
 
   Double_t evaluate() const ;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 //   void initGenerator();
 
 private:

roofit/roofit/inc/RooBernstein.h

@@ -40,6 +40,8 @@ class RooBernstein : public RooAbsPdf {
   RooListProxy _coefList ;
 
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
   ClassDef(RooBernstein,1) // Bernstein polynomial PDF
 };

roofit/roofit/inc/RooBifurGauss.h

@@ -43,6 +43,8 @@ class RooBifurGauss : public RooAbsPdf {
   RooRealProxy sigmaR;
 
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
 private:
 

roofit/roofit/inc/RooBreitWigner.h

@@ -41,6 +41,7 @@ class RooBreitWigner : public RooAbsPdf {
   RooRealProxy width ;
 
   Double_t evaluate() const ;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
 
 //   void initGenerator();
 //   Int_t generateDependents();

roofit/roofit/inc/RooCBShape.h

@@ -51,6 +51,7 @@ class RooCBShape : public RooAbsPdf {
   RooRealProxy n;
 
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
 
 private:
 

roofit/roofit/inc/RooGamma.h

@@ -42,6 +42,8 @@ class RooGamma : public RooAbsPdf {
   RooRealProxy mu ;
 
   Double_t evaluate() const ;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
 private:
 

roofit/roofit/inc/RooLegendre.h

@@ -42,6 +42,7 @@ class RooLegendre : public RooAbsReal {
   int _l2,_m2;
 
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
 
   ClassDef(RooLegendre,1) // Legendre polynomial
 };

roofit/roofit/inc/RooPolynomial.h

@@ -50,6 +50,8 @@ class RooPolynomial : public RooAbsPdf {
 
   /// Evaluation
   Double_t evaluate() const;
+  RooSpan<double> evaluateBatch(std::size_t begin, std::size_t batchSize) const;
+
 
   ClassDef(RooPolynomial,1) // Polynomial PDF
 };

roofit/roofit/src/BatchHelpers.cxx

@@ -34,7 +34,7 @@ size_t findSize(std::vector< RooSpan<const double> > parameters)
 /* This function returns the minimum size of the non-zero-sized batches
  * as well as the number of parameters that are batches, wrapped in a
  * EvaluateInfo struct (see BatchHelpers.h). It will be used when the 
- * number of parameters is > 3 and the BracketAdapterWithBranch will be used.
+ * number of parameters is > 3 and the BracketAdapterWithMask will be used.
  */
 EvaluateInfo getInfo(std::vector<const RooRealProxy*> parameters, size_t begin, size_t batchSize)
 {
@@ -49,33 +49,4 @@ EvaluateInfo getInfo(std::vector<const RooRealProxy*> parameters, size_t begin,
   return ret;
 }
 
-/* This function returns the minimum size of the non-zero-sized batches
- * as well as the number of parameters that are batches, wrapped in a
- * EvaluateInfo struct (see BatchHelpers.h). The function assigns the pointers
- * the span pointer value (in case it's a batch), or fills an array with the 
- * corresponding parameter value and sets the pointer to the array (in case it's
- * a scalar). init is used when the number of parameters is > 3 and the two indices
- * trick is used.
- */
-EvaluateInfo init(std::vector< RooRealProxy > parameters, 
-                  std::vector< ArrayWrapper* > wrappers,
-                  std::vector< double* > arrays,
-                  size_t begin, size_t batchSize )
-{
-  EvaluateInfo ret = {SIZE_MAX, 0};
-  for (size_t i=0; i<parameters.size(); i++) {
-    RooSpan<const double> span = parameters[i].getValBatch(begin,batchSize);
-    wrappers[i]->_batch = !span.empty();
-    if ( span.empty() ) {
-      for (size_t j=0; j<block; j++) arrays[i][j] = parameters[i];
-      wrappers[i]->ptr = arrays[i];
-    }
-    else {
-      wrappers[i]->ptr = span.data();
-      ret.nBatches++;
-      if (span.size() < ret.size) ret.size = span.size();
-    }
-  }
-  return ret;
-}
-}
+};

roofit/roofit/src/BatchHelpers.h

@@ -45,6 +45,7 @@ struct EvaluateInfo {
 };
 
 size_t findSize(std::vector< RooSpan<const double> > parameters);
+
 EvaluateInfo getInfo(std::vector<const RooRealProxy*> parameters, size_t begin, size_t batchSize);
 EvaluateInfo init(std::vector< RooRealProxy > parameters, 
                   std::vector<  ArrayWrapper* > wrappers,
@@ -63,29 +64,39 @@ class BracketAdapter {
     constexpr double operator[](std::size_t) const {
       return _payload;
     }
+
+    constexpr bool isBatch() const noexcept {
+      return false;
+    }
 
   private:
     const T _payload;
 };
 
 
-class BracketAdapterWithBranch {
+class BracketAdapterWithMask {
   public:
-    explicit BracketAdapterWithBranch(double payload, const RooSpan<const double>& batch) noexcept :
+    explicit BracketAdapterWithMask(double payload, const RooSpan<const double>& batch) noexcept :
+    _isBatch(!batch.empty()),
     _payload(payload),
-    _pointer(batch.data()),
-    _batchEmpty(batch.empty())
+    _pointer(batch.empty() ? &_payload : batch.data()),
+    _mask(batch.empty() ? 0 : ~static_cast<size_t>(0))
     {
     }
 
-    inline double operator[](std::size_t i) const noexcept {
-      return _batchEmpty ? _payload : _pointer[i];
+    constexpr double operator[](std::size_t i) const noexcept {
+      return _pointer[ i & _mask];
+    }
+
+    constexpr bool isBatch() const noexcept {
+      return _isBatch;
     }
 
   private:
+    const bool _isBatch;
     const double _payload;
     const double* __restrict__ const _pointer;
-    const bool _batchEmpty;
+    const size_t _mask;
 };
 
 }

roofit/roofit/src/RooArgusBG.cxx

@@ -38,6 +38,9 @@ RooArgusBG is a RooAbsPdf implementation describing the ARGUS background shape.
 #include "TMath.h"
 
 #include "TError.h"
+#include "BatchHelpers.h"
+#include "vdt/exp.h"
+#include "vdt/log.h"
 
 using namespace std;
 
@@ -94,6 +97,54 @@ Double_t RooArgusBG::evaluate() const {
 
 ////////////////////////////////////////////////////////////////////////////////
 
+namespace ArgusBatchEvaluate {
+//Author: Emmanouil Michalainas, CERN 19 AUGUST 2019  
+
+template<class Tm, class Tm0, class Tc, class Tp>
+void compute(  size_t batchSize,
+               double * __restrict__ output,
+               Tm M, Tm0 M0, Tc C, Tp P)
+{
+  for (size_t i=0; i<batchSize; i++) {
+    const double t = M[i]/M0[i];
+    const double u = 1 - t*t;
+    output[i] = C[i]*u + P[i]*vdt::fast_log(u);
+  }
+  for (size_t i=0; i<batchSize; i++) {
+    if (M[i] >= M0[i]) output[i] = 0.0;
+    else output[i] = M[i]*vdt::fast_exp(output[i]);
+  }
+}
+};
+
+RooSpan<double> RooArgusBG::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  using namespace BatchHelpers;
+  using namespace ArgusBatchEvaluate;
+
+  EvaluateInfo info = getInfo( {&m,&m0,&c,&p}, begin, batchSize );
+  auto output = _batchData.makeWritableBatchUnInit(begin, info.size);
+  auto mData = m.getValBatch(begin, info.size);
+  if (info.nBatches == 0) {
+    throw std::logic_error("Requested a batch computation, but no batch data available.");
+  }
+  else if (info.nBatches==1 && !mData.empty()) {
+    compute(info.size, output.data(), mData.data(), 
+    BracketAdapter<double> (m0), 
+    BracketAdapter<double> (c), 
+    BracketAdapter<double> (p));
+  }
+  else {
+    compute(info.size, output.data(), 
+    BracketAdapterWithMask (m,m.getValBatch(begin,batchSize)), 
+    BracketAdapterWithMask (m0,m0.getValBatch(begin,batchSize)), 
+    BracketAdapterWithMask (c,c.getValBatch(begin,batchSize)), 
+    BracketAdapterWithMask (p,p.getValBatch(begin,batchSize)));
+  }
+  return output;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
 Int_t RooArgusBG::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
 {
   if (p.arg().isConstant()) {

roofit/roofit/src/RooBernstein.cxx

@@ -131,6 +131,85 @@ Double_t RooBernstein::evaluate() const
   return TMath::SignalingNaN();
 }
 
+////////////////////////////////////////////////////////////////////////////////
+
+namespace BernsteinEvaluate {
+//Author: Emmanouil Michalainas, CERN 16 AUGUST 2019  
+
+void compute(  size_t batchSize, double xmax, double xmin,
+               double * __restrict__ output,
+               const double * __restrict__ const xData,
+               const RooListProxy& coefList)
+{
+  constexpr size_t block = 128;
+  const int nCoef = coefList.size();
+  const int degree = nCoef-1;
+  double X[block], _1_X[block], powX[block], pow_1_X[block];
+  double *Binomial = new double[nCoef+5];
+  //Binomial stores values c(degree,i) for i in [0..degree]
+
+  Binomial[0] = 1.0;
+  for (int i=1; i<=degree; i++) {
+    Binomial[i] = Binomial[i-1]*(degree-i+1)/i;
+  }
+
+  for (size_t i=0; i<batchSize; i+=block) {
+    const size_t stop = (i+block > batchSize) ? batchSize-i : block;
+
+    //initialization
+    for (size_t j=0; j<stop; j++) {
+      powX[j] = pow_1_X[j] = 1.0;
+      X[j] = (xData[i+j]-xmin) / (xmax-xmin);
+      _1_X[j] = 1-X[j];
+      output[i+j] = 0.0;
+    }
+
+    //raising 1-x to the power of degree
+    for (int k=2; k<=degree; k+=2) 
+      for (size_t j=0; j<stop; j++) 
+        pow_1_X[j] *= _1_X[j]*_1_X[j];
+
+    if (degree%2 == 1)
+      for (size_t j=0; j<stop; j++) 
+        pow_1_X[j] *= _1_X[j];
+
+    //inverting 1-x ---> 1/(1-x)
+    for (size_t j=0; j<stop; j++) 
+      _1_X[j] = 1/_1_X[j];
+
+    for (int k=0; k<nCoef; k++) {
+      double coef = static_cast<RooAbsReal&>(coefList[k]).getVal();
+      for (size_t j=0; j<stop; j++) {
+        output[i+j] += coef*Binomial[k]*powX[j]*pow_1_X[j];
+
+        //calculating next power for x and 1-x
+        powX[j] *= X[j];
+        pow_1_X[j] *= _1_X[j];
+      }
+    }
+  }
+  delete[] Binomial;
+}
+};
+
+////////////////////////////////////////////////////////////////////////////////
+
+RooSpan<double> RooBernstein::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
+  auto xData = _x.getValBatch(begin, batchSize);
+  batchSize = xData.size();
+  auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
+
+  if (xData.empty()) {
+        throw std::logic_error("Requested a batch computation, but no batch data available.");
+  }
+  else {
+    const double xmax = _x.max();
+    const double xmin = _x.min();
+    BernsteinEvaluate::compute(batchSize, xmax, xmin, output.data(), xData.data(), _coefList);
+  }
+  return output;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// No analytical calculation available (yet) of integrals over subranges