Add codegen support for VerticalInterpPdf
#1060
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,6 +1,12 @@ | ||
| #ifndef CombineMathFuncs_h | ||
| #define CombineMathFuncs_h | ||
|
|
||
| #include <RooAbsReal.h> | ||
| #include <RooArgList.h> | ||
| #include <RooArgSet.h> | ||
| #include <RooConstVar.h> | ||
| #include <RtypesCore.h> | ||
|
|
||
| #include <cmath> | ||
|
|
||
| namespace RooFit { | ||
|
|
@@ -145,6 +151,161 @@ 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) | ||
|
There was a problem hiding this comment. warning: no header providing "Int_t" is directly included [misc-include-cleaner] 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); | ||
| } | ||
| // 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); | ||
| } | ||
| 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); | ||
| } | ||
| // 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 = &dynamic_cast<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 = dynamic_cast<RooAbsReal&>(coefList[iCoef]).getVal(normSet2); | ||
| RooAbsReal* funcUp = &dynamic_cast<RooAbsReal&>(funcList[(2 * iCoef) + 1]); | ||
| RooAbsReal* funcDn = &dynamic_cast<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]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] Double_t central = funcList[0];
^ |
||
| Double_t value = central; | ||
|
|
||
| for (std::size_t iCoef = 0; iCoef < nCoeffs; ++iCoef) { | ||
| double coefVal = coefList[iCoef]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] double coefVal = coefList[iCoef];
^ |
||
| double funcUp = funcList[(2 * iCoef) + 1]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] double funcUp = funcList[(2 * iCoef) + 1];
^ |
||
| double funcDn = funcList[(2 * iCoef) + 2]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] 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]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] Double_t central = funcList[0];
^ |
||
| Double_t value = central; | ||
|
|
||
| for (std::size_t iCoef = 0; iCoef < nCoeffs; ++iCoef) { | ||
| double coefVal = coefList[iCoef]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] double coefVal = coefList[iCoef];
^ |
||
| double funcUp = funcList[(2 * iCoef) + 1]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] double funcUp = funcList[(2 * iCoef) + 1];
^ |
||
| double funcDn = funcList[(2 * iCoef) + 2]; | ||
|
There was a problem hiding this comment. warning: do not use pointer arithmetic [cppcoreguidelines-pro-bounds-pointer-arithmetic] 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 | ||
|
|
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
warning: no header providing "Double_t" is directly included [misc-include-cleaner]
interface/CombineMathFuncs.h:7:
+ #include <RtypesCore.h>