Add codegen support for VerticalInterpPdf

.github/workflows/clang-tidy-review-post.yml

@@ -21,10 +21,10 @@ jobs:
     steps:
       - name: Post review comments
         id: post-review
-        uses: ZedThree/clang-tidy-review/post@v0.20.1
+        uses: ZedThree/clang-tidy-review/post@v0.21.0
         with:
           max_comments: 10
 
        # If there are any comments, fail the check
       - if: steps.post-review.outputs.total_comments > 0
-        run: exit 1
+        run: exit 1

.github/workflows/clang-tidy-review.yml

@@ -31,7 +31,7 @@ jobs:
         run: pip install lit
 
       - name: Run clang-tidy
-        uses: ZedThree/clang-tidy-review@v0.20.1
+        uses: ZedThree/clang-tidy-review@v0.21.0
         id: review
         with:
           build_dir: build
@@ -55,4 +55,4 @@ jobs:
             micromamba run -n clang-tidy-review cmake . -B build -DINSTALL_PYTHON=FALSE -DUSE_VDT=OFF -DCMAKE_EXPORT_COMPILE_COMMANDS=ON -DCMAKE_PREFIX_PATH=/github/home/micromamba/envs/clang-tidy-review
 
       - name: Upload artifacts
-        uses: ZedThree/clang-tidy-review/upload@v0.20.1
+        uses: ZedThree/clang-tidy-review/upload@v0.21.0

interface/CombineCodegenImpl.h

@@ -2,17 +2,28 @@
 #define HiggsAnalysis_CombinedLimit_CombineCodegenImpl_h
 
 #include <ROOT/RConfig.hxx> // for ROOT_VERSION
+#include <RooAbsReal.h>
+#include <string>
 
 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0)
 # define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME) \
     namespace RooFit { namespace Experimental { void codegenImpl(CLASS_NAME &arg, CodegenContext &ctx); }}
+# define COMBINE_DECLARE_CODEGEN_INTEGRAL_IMPL(CLASS_NAME) \
+    namespace RooFit { namespace Experimental { std::string codegenIntegralImpl(CLASS_NAME &arg, int code, const char *rangeName, CodegenContext &ctx); }}
 # define COMBINE_DECLARE_TRANSLATE
+# define COMBINE_DECLARE_ANALYTICAL_INTEGRAL
 #elif ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0)
 # define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME)
-# define COMBINE_DECLARE_TRANSLATE void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+# define COMBINE_DECLARE_CODEGEN_INTEGRAL_IMPL(CLASS_NAME)
+# define COMBINE_DECLARE_TRANSLATE \
+    void translate(RooFit::Detail::CodeSquashContext &ctx) const override;
+# define COMBINE_DECLARE_ANALYTICAL_INTEGRAL \
+    std::string buildCallToAnalyticIntegral(Int_t code, const char *rangeName, RooFit::Detail::CodeSquashContext &ctx) const override;
 #else
 # define COMBINE_DECLARE_CODEGEN_IMPL(_)
+# define COMBINE_DECLARE_CODEGEN_INTEGRAL_IMPL(_)
 # define COMBINE_DECLARE_TRANSLATE
+# define COMBINE_DECLARE_ANALYTICAL_INTEGRAL
 #endif
 
 #endif

interface/CombineMathFuncs.h

@@ -1,6 +1,11 @@
 #ifndef CombineMathFuncs_h
 #define CombineMathFuncs_h
 
+#include <RooAbsReal.h>
+#include <RooArgList.h>
+#include <RooArgSet.h>
+#include <RooConstVar.h>
+
 #include <cmath>
 
 namespace RooFit {
@@ -145,6 +150,162 @@ inline double processNormalization(double nominalValue, std::size_t nThetas, std
    return norm;
 }
 
+// Interpolation  (from VerticalInterpPdf)
+inline Double_t interpolate(Double_t const coeff, Double_t const central, Double_t const fUp, 
+                            Double_t const fDn, Double_t const quadraticRegion, Int_t const quadraticAlgo)
+{
+   if (quadraticAlgo == -1) {
+      Double_t kappa = (coeff > 0 ? fUp/central : central/fDn);
+      return pow(kappa, sqrt(pow(coeff, 2)));
+   }
+
+   if (fabs(coeff) >= quadraticRegion) {
+      return coeff * (coeff > 0 ? fUp - central : central - fDn);
+   } else {
+      // quadratic interpolation coefficients between the three
+      if (quadraticAlgo == 0) {
+         // quadratic interpolation null at zero and continuous at boundaries, but not differentiable at boundaries
+         // conditions:
+         //   c_up (+quadraticRegion) = +quadraticRegion
+         //   c_cen(+quadraticRegion) = -quadraticRegion
+         //   c_dn (+quadraticRegion) = 0
+         //   c_up (-quadraticRegion) = 0 
+         //   c_cen(-quadraticRegion) = -quadraticRegion
+         //   c_dn (-quadraticRegion) = +quadraticRegion
+         //   c_up(0) = c_dn(0) = c_cen(0) = 0
+         Double_t c_up  = + coeff * (quadraticRegion + coeff) / (2 * quadraticRegion);
+         Double_t c_dn  = - coeff * (quadraticRegion - coeff) / (2 * quadraticRegion);
+         Double_t c_cen = - coeff * coeff / quadraticRegion;
+         return c_up * fUp + c_dn * fDn + c_cen * central;
+      } else if (quadraticAlgo == 1) { 
+         // quadratic interpolation that is everywhere differentiable, but it's not null at zero
+         // conditions on the function
+         //   c_up (+quadraticRegion) = +quadraticRegion
+         //   c_cen(+quadraticRegion) = -quadraticRegion
+         //   c_dn (+quadraticRegion) = 0
+         //   c_up (-quadraticRegion) = 0 
+         //   c_cen(-quadraticRegion) = -quadraticRegion
+         //   c_dn (-quadraticRegion) = +quadraticRegion
+         // conditions on the derivatives
+         //   c_up '(+quadraticRegion) = +1
+         //   c_cen'(+quadraticRegion) = -1
+         //   c_dn '(+quadraticRegion) = 0
+         //   c_up '(-quadraticRegion) = 0
+         //   c_cen'(-quadraticRegion) = +1
+         //   c_dn '(-quadraticRegion) = -1
+         Double_t c_up  = (quadraticRegion + coeff) * (quadraticRegion + coeff) / (4 * quadraticRegion);
+         Double_t c_dn  = (quadraticRegion - coeff) * (quadraticRegion - coeff) / (4 * quadraticRegion);
+         Double_t c_cen = - c_up - c_dn;
+         return c_up * fUp + c_dn * fDn + c_cen * central;
+      } else/* if (quadraticAlgo == 1)*/ {
+         // P(6) interpolation that is everywhere differentiable and null at zero
+         /* === how the algorithm works, in theory ===
+         * let  dhi = h_hi - h_nominal
+         *      dlo = h_lo - h_nominal
+         * and x be the morphing parameter
+         * we define alpha = x * 0.5 * ((dhi-dlo) + (dhi+dlo)*smoothStepFunc(x));
+         * which satisfies:
+         *     alpha(0) = 0
+         *     alpha(+1) = dhi
+         *     alpha(-1) = dlo
+         *     alpha(x >= +1) = |x|*dhi
+         *     alpha(x <= -1) = |x|*dlo
+         *     alpha is continuous and has continuous first and second derivative, as smoothStepFunc has them
+         * === and in practice ===
+         * we already have computed the histogram for diff=(dhi-dlo) and sum=(dhi+dlo)
+         * so we just do template += (0.5 * x) * (diff + smoothStepFunc(x) * sum)
+         * ========================================== */
+         Double_t cnorm = coeff/quadraticRegion;
+         Double_t cnorm2 = pow(cnorm, 2);
+         Double_t hi = fUp - central;
+         Double_t lo = fDn - central;
+         Double_t sum = hi+lo;
+         Double_t diff = hi-lo;
+         Double_t a = coeff/2.; // cnorm*quadraticRegion
+         Double_t b = 0.125 * cnorm * (cnorm2 * (3.*cnorm2 - 10.) + 15.);
+         Double_t result = a*(diff + b*sum);
+         return result;
+      }
+   }
+}
+
+template <typename Operation>
+inline Double_t opInterpolate(RooArgList const& coefList, RooArgList const& funcList, Double_t const pdfFloorVal,
+                                    Double_t const quadraticRegion, Int_t const quadraticAlgo, const RooArgSet* normSet2=nullptr)
+{
+   // Do running sum of coef/func pairs, calculate lastCoef.
+   RooAbsReal* func = &(RooAbsReal&)funcList[0];
+   Double_t central = func->getVal();
+   Double_t value = central;
+
+   Operation op;
+
+   for (int iCoef = 0; iCoef < coefList.getSize(); ++iCoef) {
+      Double_t coefVal = static_cast<RooAbsReal&>(coefList[iCoef]).getVal(normSet2) ;
+      RooAbsReal* funcUp = &(RooAbsReal&)funcList[2 * iCoef + 1];
+      RooAbsReal* funcDn = &(RooAbsReal&)funcList[2 * iCoef + 2];
+      value = op(value, interpolate(coefVal, central, funcUp->getVal(), funcDn->getVal(), quadraticRegion, quadraticAlgo));
+   }
+   return ( value > 0. ? value : pdfFloorVal);
+}
+
+inline Double_t additiveInterpolate(double const* coefList, std::size_t nCoeffs, double const* funcList, std::size_t nFuncs,
+                                  Double_t const pdfFloorVal, Double_t const quadraticRegion, Int_t const quadraticAlgo)
+{
+   // Do running sum of coef/func pairs, calculate lastCoef.
+   Double_t central = funcList[0];
+   Double_t value = central;
+
+   for (std::size_t iCoef = 0; iCoef < nCoeffs; ++iCoef) {
+      double coefVal = coefList[iCoef];
+      double funcUp = funcList[2 * iCoef + 1];
+      double funcDn = funcList[2 * iCoef + 2];
+      value += interpolate(coefVal, central, funcUp, funcDn, quadraticRegion, quadraticAlgo);
+   }
+   return ( value > 0. ? value : pdfFloorVal);
+}
+
+inline Double_t multiplicativeInterpolate(double const* coefList, std::size_t nCoeffs, double const* funcList, std::size_t nFuncs,
+                                  Double_t const pdfFloorVal, Double_t const quadraticRegion, Int_t const quadraticAlgo)
+{
+   // Do running sum of coef/func pairs, calculate lastCoef.
+   Double_t central = funcList[0];
+   Double_t value = central;
+
+   for (std::size_t iCoef = 0; iCoef < nCoeffs; ++iCoef) {
+      double coefVal = coefList[iCoef];
+      double funcUp = funcList[2 * iCoef + 1];
+      double funcDn = funcList[2 * iCoef + 2];
+      value *= interpolate(coefVal, central, funcUp, funcDn, quadraticRegion, quadraticAlgo);
+   }
+   return ( value > 0. ? value : pdfFloorVal);
+}
+
+inline Double_t verticalInterpolate(double const* coefList, std::size_t nCoeffs, double const* funcList, std::size_t nFuncs,
+                                    double const pdfFloorVal, double const quadraticRegion, Int_t const quadraticAlgo)
+{
+   // Do running sum of coef/func pairs, calculate lastCoef.
+   Double_t value = pdfFloorVal;
+   if (quadraticAlgo >= 0) {
+      value = RooFit::Detail::MathFuncs::additiveInterpolate(coefList, nCoeffs, funcList, nFuncs, pdfFloorVal, quadraticRegion, quadraticAlgo);
+   } else {
+      value = RooFit::Detail::MathFuncs::multiplicativeInterpolate(coefList, nCoeffs, funcList, nFuncs, pdfFloorVal, quadraticRegion, quadraticAlgo);
+   }
+   return value;
+}
+
+inline Double_t verticalInterpPdfIntegral(double const* coefList, std::size_t nCoeffs, double const* funcIntList, std::size_t nFuncs,
+                                          double const pdfFloorVal, double const integralFloorVal,
+                                          double const quadraticRegion, Int_t const quadraticAlgo)
+{
+   double value = RooFit::Detail::MathFuncs::additiveInterpolate(coefList, nCoeffs, funcIntList, nFuncs, 
+                                                                 pdfFloorVal, quadraticRegion, quadraticAlgo);
+   double normVal(1);
+   double result = 0;
+   if(normVal>0.) result = value / normVal;
+   return result > 0. ? result : integralFloorVal;
+}
+
 } // namespace MathFuncs
 } // namespace Detail
 } // namespace RooFit

interface/VerticalInterpPdf.h

@@ -8,6 +8,8 @@
 #include "RooListProxy.h"
 #include "RooObjCacheManager.h"
 
+#include "CombineCodegenImpl.h"
+
 class VerticalInterpPdf : public RooAbsPdf {
 public:
 
@@ -18,17 +20,23 @@ class VerticalInterpPdf : public RooAbsPdf {
 
   Double_t evaluate() const override ;
   Bool_t checkObservables(const RooArgSet* nset) const override ;	
+  COMBINE_DECLARE_TRANSLATE
 
   Bool_t forceAnalyticalInt(const RooAbsArg&) const override { return kTRUE ; }
   Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& numVars, const RooArgSet* normSet, const char* rangeName=0) const override ;
   Double_t analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName=0) const override ;
+  COMBINE_DECLARE_ANALYTICAL_INTEGRAL
 
   const RooArgList& funcList() const { return _funcList ; }
+  const RooArgList& funcIntListFromCache() const;
   const RooArgList& coefList() const { return _coefList ; }
 
   const Double_t quadraticRegion() const { return _quadraticRegion; }
   const Int_t quadraticAlgo() const { return _quadraticAlgo; }
 
+  const Double_t pdfFloorVal() const { return _pdfFloorVal; }
+  const Double_t integralFloorVal() const { return _integralFloorVal; }
+
   void setFloorVals(Double_t const& pdf_val, Double_t const& integral_val);
 
 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,34,06)
@@ -62,4 +70,7 @@ class VerticalInterpPdf : public RooAbsPdf {
   ClassDefOverride(VerticalInterpPdf,3) // PDF constructed from a sum of (non-pdf) functions
 };
 
+COMBINE_DECLARE_CODEGEN_IMPL(VerticalInterpPdf);
+COMBINE_DECLARE_CODEGEN_INTEGRAL_IMPL(VerticalInterpPdf);
+
 #endif

src/CombineCodegenImpl.cxx

@@ -5,16 +5,20 @@
 #include "../interface/AsymPow.h"
 #include "../interface/ProcessNormalization.h"
 #include "../interface/VerticalInterpHistPdf.h"
+#include "../interface/VerticalInterpPdf.h"
+#include "../interface/CombineMathFuncs.h"
 
 #include <RooUniformBinning.h>
 
 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0)
 namespace RooFit {
 namespace Experimental {
 # define CODEGEN_IMPL(CLASS_NAME) void codegenImpl(CLASS_NAME &arg0, CodegenContext &ctx)
+# define CODEGEN_INTEGRAL_IMPL(CLASS_NAME) std::string codegenIntegralImpl(CLASS_NAME &arg0, int code, const char *rangeName, CodegenContext &ctx)
 # define ARG_VAR auto &arg = arg0;
 #else
 # define CODEGEN_IMPL(CLASS_NAME) void CLASS_NAME::translate(RooFit::Detail::CodeSquashContext &ctx) const
+# define CODEGEN_INTEGRAL_IMPL(CLASS_NAME) std::string CLASS_NAME::buildCallToAnalyticIntegral(Int_t code, const char *rangeName, RooFit::Detail::CodeSquashContext &ctx) const
 # define ARG_VAR auto &arg = *this;
 #endif
 
@@ -217,6 +221,33 @@ CODEGEN_IMPL(FastVerticalInterpHistPdf2D2) {
   ctx.addResult(&arg, arrName + "[" + binIdx.str() + "]");
 }
 
+CODEGEN_IMPL(VerticalInterpPdf) {
+  ARG_VAR;
+  ctx.addResult(&arg,
+                ctx.buildCall("RooFit::Detail::MathFuncs::verticalInterpolate",
+                              arg.coefList(),
+                              arg.coefList().size(),
+                              arg.funcList(),
+                              arg.funcList().size(),
+                              arg.pdfFloorVal(),
+                              arg.quadraticRegion(),
+                              arg.quadraticAlgo()));
+
+}
+
+CODEGEN_INTEGRAL_IMPL(VerticalInterpPdf) {
+  ARG_VAR;
+  return ctx.buildCall("RooFit::Detail::MathFuncs::verticalInterpPdfIntegral",
+                       arg.coefList(),
+                       arg.coefList().size(),
+                       arg.funcIntListFromCache(),
+                       arg.funcIntListFromCache().size(),
+                       arg.pdfFloorVal(),
+                       arg.integralFloorVal(),
+                       arg.quadraticRegion(),
+                       arg.quadraticAlgo());
+}
+
 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0)
 } // namespace RooFit
 } // namespace Experimental