Make it possible to build combine without VDT via CMake

This makes it possible to build combine in environments where VDT is not
available.

Nothing is changed for the default code, e.g. the case where
`COMBINE_NO_VDT` is not defined.

The GBRMath header was also generalized to a meta-header for math
includes plus some preprocessor checks and defines dpeneding on VDT
is requested or not.

Author: guitargeek

CMakeLists.txt

@@ -4,6 +4,7 @@ project(HiggsAnalysisCombinedLimit VERSION 0.0.1)
 
 option( MODIFY_ROOTMAP "Modify generated Rootmap to take out classes already bundled in StatAnalysis" FALSE )
 option( INSTALL_PYTHON "Install the Python library and scripts" TRUE )
+option( USE_VDT "Use VDT (fast and vectorisable mathematical functions)" TRUE )
 
 # Can build with CMake after e.g. setting up StatAnalysis release like this:
 # export ATLAS_LOCAL_ROOT_BASE=/cvmfs/atlas.cern.ch/repo/ATLASLocalRootBase
@@ -41,7 +42,13 @@ ROOT_GENERATE_DICTIONARY(G__${LIBNAME} HiggsAnalysis/CombinedLimit/src/classes.h
         OPTIONS ${ROOTCLING_OPTIONS})
 add_library(${LIBNAME} SHARED ${SOURCES} G__${LIBNAME}.cxx)
 set_target_properties(${LIBNAME} PROPERTIES PUBLIC_HEADER "${HEADERS}")
-target_link_libraries (${LIBNAME} Eigen3::Eigen ${ROOT_LIBRARIES} ${Boost_LIBRARIES} VDT::VDT)
+target_link_libraries (${LIBNAME} Eigen3::Eigen ${ROOT_LIBRARIES} ${Boost_LIBRARIES})
+
+if(NOT USE_VDT)
+  target_compile_definitions(${LIBNAME} PUBLIC COMBINE_NO_VDT)
+else()
+  target_link_libraries (${LIBNAME} VDT::VDT)
+endif()
 
 add_executable(combine bin/combine.cpp)
 target_link_libraries(combine PUBLIC ${LIBNAME})

interface/GBRMath.h

@@ -0,0 +1,48 @@
+#ifndef MathHeaders_h
+#define MathHeaders_h
+
+#ifndef COMBINE_NO_VDT
+#include "vdt/vdtMath.h"
+#endif
+
+#include <cmath>
+
+#ifdef COMBINE_NO_VDT
+
+#define my_exp std::exp
+#define my_log std::log
+inline double my_inv(double x) { return 1. / x; }
+
+#else
+
+#define my_exp vdt::fast_exp
+#define my_log vdt::fast_log
+#define my_inv vdt::fast_inv
+
+#endif
+
+namespace gbrmath {
+
+  inline double fast_pow(double base, double exponent) {
+    if (base == 0. && exponent > 0.)
+      return 0.;
+#ifdef COMBINE_NO_VDT
+    else if (base > 0.)
+      return std::exp(exponent * std::log(base));
+#else
+    else if (base > 0.)
+      return vdt::fast_exp(exponent * vdt::fast_log(base));
+#endif
+    else
+      return std::nan("");
+  }
+
+  //   inline float fast_powf(float base, float exponent) {
+  //     if (base==0. && exponent>0.) return 0.;
+  //     else if (base>0.) return vdt::fast_expf(exponent*vdt::fast_logf(base));
+  //     else return std::nanf("");
+  //   }
+
+}  // namespace gbrmath
+
+#endif

src/GBRMath.cc

@@ -2,13 +2,10 @@
 #include <math.h>
 #include "TMath.h"
 
-//#include "../interface/RooDoubleCBFast.h"
-//#include "../interface/RooFermi.h"
-//#include "../interface/RooRelBW.h"
 #include "../interface/RooDoubleCBFast.h"
 #include "RooRealVar.h"
 #include "RooRealConstant.h"
-#include "../interface/GBRMath.h"
+#include "./MathHeaders.h"
 
 using namespace RooFit;
 
@@ -52,21 +49,21 @@ using namespace RooFit;
 
  double RooDoubleCBFast::evaluate() const 
  { 
-   double t = (x-mean)*vdt::fast_inv(width);
+   double t = (x-mean)*my_inv(width);
    double val = -99.;
    if(t>-alpha1 && t<alpha2){
-     val = vdt::fast_exp(-0.5*t*t);
+     val = my_exp(-0.5*t*t);
    }else if(t<=-alpha1){
-     double alpha1invn1 = alpha1*vdt::fast_inv(n1);
-     val = vdt::fast_exp(-0.5*alpha1*alpha1)*gbrmath::fast_pow(1. - alpha1invn1*(alpha1+t), -n1);
+     double alpha1invn1 = alpha1*my_inv(n1);
+     val = my_exp(-0.5*alpha1*alpha1)*gbrmath::fast_pow(1. - alpha1invn1*(alpha1+t), -n1);
 
 //      double n1invalpha1 = n1*vdt::fast_inv(fabs(alpha1));
 //      double A1 = gbrmath::fast_pow(n1invalpha1,n1)*vdt::fast_exp(-alpha1*alpha1/2.);
 //      double B1 = n1invalpha1-fabs(alpha1);
 //      val = A1*gbrmath::fast_pow(B1-t,-n1);
    }else if(t>=alpha2){
-     double alpha2invn2 = alpha2*vdt::fast_inv(n2);
-     val = vdt::fast_exp(-0.5*alpha2*alpha2)*gbrmath::fast_pow(1. - alpha2invn2*(alpha2-t), -n2);     
+     double alpha2invn2 = alpha2*my_inv(n2);
+     val = my_exp(-0.5*alpha2*alpha2)*gbrmath::fast_pow(1. - alpha2invn2*(alpha2-t), -n2);     
 
 //      double n2invalpha2 = n2*vdt::fast_inv(fabs(alpha2));
 //      double A2 = gbrmath::fast_pow(n2invalpha2,n2)*vdt::fast_exp(-alpha2*alpha2/2.);
@@ -111,7 +108,7 @@ using namespace RooFit;
    static const double rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
    static const double invRoot2 = 1.0/sqrt(2);   
 
-   double invwidth = vdt::fast_inv(width);
+   double invwidth = my_inv(width);
 
    double tmin = (xmin-mean)*invwidth;
    double tmax = (xmax-mean)*invwidth;
@@ -129,7 +126,7 @@ using namespace RooFit;
    }
    //compute left tail;
    if (isfullrange  && (n1-1.0)>1.e-5) {
-    left = width*vdt::fast_exp(-0.5*alpha1*alpha1)*n1*vdt::fast_inv(alpha1*(n1-1.)); 
+    left = width*my_exp(-0.5*alpha1*alpha1)*n1*my_inv(alpha1*(n1-1.)); 
    }
    else {
 
@@ -138,11 +135,11 @@ using namespace RooFit;
     double thigh = (left_high-mean)*invwidth;
 
     if(left_low < left_high){ //is the left tail in range?
-     double n1invalpha1 = n1*vdt::fast_inv(fabs(alpha1));
+     double n1invalpha1 = n1*my_inv(fabs(alpha1));
       if(fabs(n1-1.0)>1.e-5) {
-	double invn1m1 = vdt::fast_inv(n1-1.);
+	double invn1m1 = my_inv(n1-1.);
 	double leftpow = gbrmath::fast_pow(n1invalpha1,-n1*invn1m1);
-	double left0 = width*vdt::fast_exp(-0.5*alpha1*alpha1)*invn1m1;
+	double left0 = width*my_exp(-0.5*alpha1*alpha1)*invn1m1;
 	double left1, left2;
 
 	if (xmax>(mean-alpha1*width)) left1 = n1invalpha1;
@@ -158,28 +155,28 @@ using namespace RooFit;
 	//left = A1*vdt::fast_inv(-n1+1.0)*width*(gbrmath::fast_pow(B1-(left_low-mean)*invwidth,-n1+1.)-gbrmath::fast_pow(B1-(left_high-mean)*invwidth,-n1+1.));
       }
       else {
-	double A1 = gbrmath::fast_pow(n1invalpha1,n1)*vdt::fast_exp(-0.5*alpha1*alpha1);
+	double A1 = gbrmath::fast_pow(n1invalpha1,n1)*my_exp(-0.5*alpha1*alpha1);
 	double B1 = n1invalpha1-fabs(alpha1);	
-	left = A1*width*(vdt::fast_log(B1-(left_low-mean)*invwidth) - vdt::fast_log(B1-(left_high-mean)*invwidth) );
+	left = A1*width*(my_log(B1-(left_low-mean)*invwidth) - my_log(B1-(left_high-mean)*invwidth) );
       }
     }
    }
 
    //compute right tail;
    if (isfullrange && (n2-1.0)>1.e-5) {
-     right = width*vdt::fast_exp(-0.5*alpha2*alpha2)*n2*vdt::fast_inv(alpha2*(n2-1.));
+     right = width*my_exp(-0.5*alpha2*alpha2)*n2*my_inv(alpha2*(n2-1.));
    }
    else {    
     double right_low=std::max(xmin,mean + alpha2*width);
     double right_high=xmax;
     double tlow = (right_low - mean)*invwidth;
 
     if(right_low < right_high){ //is the right tail in range?
-      double n2invalpha2 = n2*vdt::fast_inv(fabs(alpha2)); 
+      double n2invalpha2 = n2*my_inv(fabs(alpha2)); 
       if(fabs(n2-1.0)>1.e-5) {
-	double invn2m2 = vdt::fast_inv(n2-1.);
+	double invn2m2 = my_inv(n2-1.);
 	double rightpow = gbrmath::fast_pow(n2invalpha2,-n2*invn2m2);
-	double right0 = width*vdt::fast_exp(-0.5*alpha2*alpha2)*invn2m2;
+	double right0 = width*my_exp(-0.5*alpha2*alpha2)*invn2m2;
 	double right1, right2;
 
 	if (xmin<(mean+alpha2*width)) right1 = n2invalpha2;
@@ -193,9 +190,9 @@ using namespace RooFit;
 	//right = A2*vdt::fast_inv(-n2+1.0)*width*(gbrmath::fast_pow(B2+(right_high-mean)*invwidth,-n2+1.)-gbrmath::fast_pow(B2+(right_low-mean)*invwidth,-n2+1.));
       }
       else {
-	double A2 = gbrmath::fast_pow(n2invalpha2,n2)*vdt::fast_exp(-0.5*alpha2*alpha2);
+	double A2 = gbrmath::fast_pow(n2invalpha2,n2)*my_exp(-0.5*alpha2*alpha2);
 	double B2 = n2invalpha2-fabs(alpha2);
-	right = A2*width*(vdt::fast_log(B2+(right_high-mean)*invwidth) - vdt::fast_log(B2+(right_low-mean)*invwidth) );
+	right = A2*width*(my_log(B2+(right_high-mean)*invwidth) - my_log(B2+(right_low-mean)*invwidth) );
       }
     }
    }

src/MathHeaders.h

@@ -4,6 +4,7 @@
 #include "../interface/ProfilingTools.h"
 #include <RooRealVar.h>
 #include <stdexcept>
+#include "./MathHeaders.h"
 
 VectorizedCBShape::VectorizedCBShape(const RooCBShape &gaus, const RooAbsData &data, bool includeZeroWeights)
 {
@@ -125,9 +126,9 @@ void VectorizedCBShape::fill(std::vector<Double_t> &out) const {
         } else {
             for (unsigned int i = 0; i < n; ++i) {
                 if(work1_[i]>-alpha1){
-                    out[i] = norm * vdt::fast_exp(-0.5*work1_[i]*work1_[i]);
+                    out[i] = norm * my_exp(-0.5*work1_[i]*work1_[i]);
                 } else {
-                    out[i] = norm2 * vdt::fast_exp(-n1 * vdt::fast_log(1. - alpha1invn1*(alpha1+work1_[i])));
+                    out[i] = norm2 * my_exp(-n1 * my_log(1. - alpha1invn1*(alpha1+work1_[i])));
                 }
             }
         }
@@ -171,17 +172,17 @@ double VectorizedCBShape::getIntegral() const {
 
     //compute left tail;
     if (isfullrange && (n1-1.0)>1.e-5) {
-        left = width*std::exp(-0.5*alpha1*alpha1)*n1*vdt::inv(alpha1*(n1-1.));
+        left = width*std::exp(-0.5*alpha1*alpha1)*n1*(1./(alpha1*(n1-1.)));
     } else {
 
         double left_low=xmin;
         double left_high=std::min(xmax,mean - alpha1*width);
         double thigh = (left_high-mean)*invwidth;
 
         if(left_low < left_high){ //is the left tail in range?
-            double n1invalpha1 = n1*vdt::inv(fabs(alpha1));
+            double n1invalpha1 = n1 / std::abs(alpha1);
             if(fabs(n1-1.0)>1.e-5) {
-                double invn1m1 = vdt::inv(n1-1.);
+                double invn1m1 = 1. / (n1-1.);
                 double leftpow = std::pow(n1invalpha1,-n1*invn1m1);
                 double left0 = width*std::exp(-0.5*alpha1*alpha1)*invn1m1;
                 double left1, left2;
@@ -210,42 +211,55 @@ double VectorizedCBShape::getIntegral() const {
 }
 
 void VectorizedCBShape::cbGauss(double* __restrict__ t, unsigned int n, double norm, double* __restrict__ out,  double* __restrict__ work2) const {
+// The fast code code path is only available when VDT is available
+#ifndef COMBINE_NO_VDT
+    if (!hasFast()) {
+#endif
+       for (unsigned int i = 0; i < n; ++i) {
+           out[i] = std::exp(-0.5 * t[i] * t[i]) * norm;
+       }
+#ifndef COMBINE_NO_VDT
+       return;
+    }
     for (unsigned int i = 0; i < n; ++i) {
         work2[i] = -0.5*t[i]*t[i];
     }
-    if (hasFast()) {
-        vdt::fast_expv(n, work2, t);
-    } else {
-        vdt::expv(n, work2, t);
-    }
+    vdt::fast_expv(n, work2, t);
     for (unsigned int i = 0; i < n; ++i) {
         out[i] = t[i]*norm;
     }
+#endif
 }
+
 void VectorizedCBShape::cbCB(double* __restrict__ t, unsigned int n, double norm, double* __restrict__ out,  double* __restrict__ work2) const {
     // val = norm * std::exp(-0.5*alpha1*alpha1) * pow(1. - alpha1/n1 * (alpha1+t), -n1);
     double alpha1 = alpha_->getVal();
-    double n1 = n_->getVal(), notn1 = -n1;
+    double n1 = n_->getVal();
+    double notn1 = -n1;
     double alpha1invn1 = alpha1/n1;
     double prefactor = norm*std::exp(-0.5*std::pow(alpha1, 2));
 
+// The fast code code path is only available when VDT is available
+#ifndef COMBINE_NO_VDT
+    if (!hasFast()) {
+#endif
+       for (unsigned int i = 0; i < n; ++i) {
+           out[i] = prefactor * std::exp(notn1 * std::log(1. - alpha1invn1*(alpha1+t[i])));
+       }
+#ifndef COMBINE_NO_VDT
+       return;
+    }
+
     for (unsigned int i = 0; i < n; ++i) {
         work2[i] = 1 - alpha1invn1*(alpha1+t[i]);
     }
-    if (hasFast()) {
-        vdt::fast_logv(n, work2, t);
-    } else {
-        vdt::logv(n, work2, t);
-    }
+    vdt::fast_logv(n, work2, t);
     for (unsigned int i = 0; i < n; ++i) {
          t[i] *= notn1;
     }
-    if (hasFast()) {
-        vdt::fast_expv(n, t, work2);
-    } else {
-        vdt::expv(n, t, work2);
-    }
+    vdt::fast_expv(n, t, work2);
     for (unsigned int i = 0; i < n; ++i) {
         out[i] = prefactor*work2[i];
     }
+#endif
 }