Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

New interpCode=6 for PiecewiseInterpolation: Multiplicative Polynominal Interpolation with Linear Extrapolation #16858

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

New interpCode=6 for PiecewiseInterpolation: Multiplicative Polynominal Interpolation with Linear Extrapolation #16858

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 25 additions & 7 deletions roofit/histfactory/src/PiecewiseInterpolation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,17 +5,35 @@
* more altered or distorted ones, it computes a new shape depending on the value of the nuisance
* parameters \f$ \alpha_i \f$:
* \f[
* A = \sum_i \mathrm{Interpolate}(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
* A = \mathrm{nominal} + \sum_i I(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
* \f]
* for additive interpolation modes (interpCodes 0, 2, 3, and 4), or:
* \f[
* A = \mathrm{nominal}\prod_i I(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
* \f]
* for multiplicative interpolation modes (interpCodes 1, 5, and 6). The interpCodes determine the function \f$ I \f$ (see table below).
*
* Note that a PiecewiseInterpolation with \f$ \mathrm{nominal}=1 \f$, N variations, and a multiplicative interpolation mode is equivalen to N
* PiecewiseInterpolations each with a single variation and the same interpolation code, all inside a RooProduct.
*
* If an \f$ \alpha_i \f$ is zero, the distribution is identical to the nominal distribution, at
* \f$ \pm 1 \f$ it is identical to the up/down distribution for that specific \f$ i \f$.
*
* PiecewiseInterpolation will behave identically (except for differences in the interpCode assignments) to a FlexibleInterpVar if both its nominal, and high and low variation sets
* are all RooRealVar.
*
* The class supports several interpolation methods, which can be selected for each parameter separately
* using setInterpCode(). The default interpolation code is 4. This performs
* - \f$ |\alpha | > 1 \f$: Linear extrapolation.
* - \f$ |\alpha | < 1 \f$: Polynomial interpolation. A sixth-order polynomial is used. Its coefficients
* are chosen such that function, first, and second derivative at \f$ \alpha \pm 1 \f$ match the values
* that the extrapolation procedure uses.
* using setInterpCode(). The default interpolation code is 0. The table below provides details of the interpCodes:

| **interpCode** | **Name** | **Description** |
|----------------|----------|-----------------|
| 0 (default) | Additive Piecewise Linear | \f$ I = \alpha_i(\mathrm{high}_i - \mathrm{nominal}) \f$ for \f$ \alpha_i>=0 \f$, otherwise \f$ I = \alpha_i(\mathrm{nominal}-\mathrm{low}_i) \f$. Not recommended except if using a symmetric variation, because of discontinuities in derivatives. |
| 1 | Multiplicative Piecewise Exponential | \f$ I = (\mathrm{high}_i/\mathrm{nominal})^{\alpha_i} \f$ for \f$ \alpha_i>=0 \f$, otherwise \f$ I = (\mathrm{low}_i/\mathrm{nominal})^{-\alpha_i} \f$. |
| 2 | Additive Quadratic Interp. + Linear Extrap. | Deprecated by interpCode 4. |
| 4 | Additive Poly Interp. + Linear Extrap. | interpCode 0 for \f$ |\alpha_i|>=1 \f$. 6th-order polynomial for \f$ |\alpha_i|<1 \f$ with matching 0th,1st,2nd derivatives at boundary. Recommended for normalization factors. |
| 5 | Multiplicative Poly Interp. + Exponential Extrap. | interpCode 1 for \f$ |\alpha_i|>=1 \f$. 6th-order polynomial for \f$ |\alpha_i|<1 \f$ with matching 0th,1st,2nd derivatives at boundary. Recommended for normalization factors. In FlexibleInterpVar this is interpCode=4. |
| 6 | Multiplicative Poly Interp. + Linear Extrap. | Multiplicative version of interpCode 4. Recommended for normalization factors that must not have roots (i.e. be equal to 0) outside of \f$ |\alpha_i|<1 \f$. |

*/

#include "RooStats/HistFactory/PiecewiseInterpolation.h"
Expand Down Expand Up @@ -576,7 +594,7 @@ void PiecewiseInterpolation::setAllInterpCodes(int code)
void PiecewiseInterpolation::setInterpCodeForParam(int iParam, int code)
{
RooAbsArg const &param = _paramSet[iParam];
if (code < 0 || code > 5) {
if (code < 0 || code > 6) {
coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName()
<< " with unknown interpolation code " << code << ", keeping current code "
<< _interpCode[iParam] << std::endl;
Expand Down
32 changes: 25 additions & 7 deletions roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,7 @@ inline unsigned int getUniformBinning(double low, double high, double val, unsig
return val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
}

inline double interpolate1d(double low, double high, double val, unsigned int numBins, double const* vals)
inline double interpolate1d(double low, double high, double val, unsigned int numBins, double const *vals)
{
double binWidth = (high - low) / numBins;
int idx = val >= high ? numBins - 1 : std::abs((val - low) / binWidth);
Expand All @@ -185,9 +185,9 @@ inline double interpolate1d(double low, double high, double val, unsigned int nu
if (val > low + 0.5 * binWidth && val < high - 0.5 * binWidth) {
double slope;
if (val < central) {
slope = vals[idx] - vals[idx - 1];
slope = vals[idx] - vals[idx - 1];
} else {
slope = vals[idx + 1] - vals[idx];
slope = vals[idx + 1] - vals[idx];
}
return vals[idx] + slope * (val - central) / binWidth;
}
Expand Down Expand Up @@ -239,10 +239,10 @@ inline double flexibleInterpSingle(unsigned int code, double low, double high, d
} else {
return a * std::pow(paramVal, 2) + b * paramVal + c;
}
// According to an old comment in the source code, code 3 was apparently
// meant to be a "parabolic version of log-normal", but it never got
// implemented. If someone would need it, it could be implemented as doing
// code 2 in log space.
// According to an old comment in the source code, code 3 was apparently
// meant to be a "parabolic version of log-normal", but it never got
// implemented. If someone would need it, it could be implemented as doing
// code 2 in log space.
} else if (code == 4) {
double x = paramVal;
if (x >= boundary) {
Expand Down Expand Up @@ -304,6 +304,24 @@ inline double flexibleInterpSingle(unsigned int code, double low, double high, d
mod = value;
}
return res * (mod - 1.0);
} else if (code == 6) {
// multiplicative version of code 4 (6th order poly interp with linear extrap)

double x = paramVal;
if (x >= boundary) {
return res * (x * (high - nominal) - 1.);
} else if (x <= -boundary) {
return res * (x * (nominal - low) - 1.);
}

// interpolate 6th degree
double t = x / boundary;
double eps_plus = high - nominal;
double eps_minus = nominal - low;
double S = 0.5 * (eps_plus + eps_minus);
double A = 0.0625 * (eps_plus - eps_minus);

return res * (x * (S + t * A * (15 + t * t * (-10 + t * t * 3))) - 1.);
}

return 0.0;
Expand Down