# Source code for pyhf.interpolators.code4

"""Polynomial Interpolation (Code 4)."""

import logging
import math
import pyhf
from pyhf.tensor.manager import get_backend
from pyhf import events
from pyhf.interpolators import _slow_interpolator_looper

log = logging.getLogger(__name__)

[docs]
class code4:
r"""
The polynomial interpolation and exponential extrapolation strategy.

.. math::
\sigma_{sb} (\vec{\alpha}) = \sigma_{sb}^0(\vec{\alpha}) \underbrace{\prod_{p \in \text{Syst}} I_\text{poly|exp.} (\alpha_p; \sigma_{sb}^0, \sigma_{psb}^+, \sigma_{psb}^-, \alpha_0)}_\text{factors to calculate}

with

.. math::
I_\text{poly|exp.}(\alpha; I^0, I^+, I^-, \alpha_0) = \begin{cases} \left(\frac{I^+}{I^0}\right)^{\alpha} \qquad \alpha \geq \alpha_0\\ 1 + \sum_{i=1}^6 a_i \alpha^i \qquad |\alpha| < \alpha_0 \\ \left(\frac{I^-}{I^0}\right)^{-\alpha} \qquad \alpha < -\alpha_0 \end{cases}

and the :math:a_i are fixed by the boundary conditions

.. math::
\sigma_{sb}(\alpha=\pm\alpha_0), \left.\frac{\mathrm{d}\sigma_{sb}}{\mathrm{d}\alpha}\right|_{\alpha=\pm\alpha_0}, \mathrm{ and } \left.\frac{\mathrm{d}^2\sigma_{sb}}{\mathrm{d}\alpha^2}\right|_{\alpha=\pm\alpha_0}.

Namely that :math:\sigma_{sb}(\vec{\alpha}) is continuous, and its first- and second-order derivatives are continuous as well.

"""

[docs]
def __init__(self, histogramssets, subscribe=True, alpha0=1):
"""Polynomial Interpolation."""
default_backend = pyhf.default_backend

# alpha0 is assumed to be positive and non-zero. If alpha0 == 0, then
# we cannot calculate the coefficients (e.g. determinant == 0)
assert alpha0 > 0
self.__alpha0 = alpha0

self._histogramssets = default_backend.astensor(histogramssets)
# initial shape will be (nsysts, 1)
self.alphasets_shape = (self._histogramssets.shape[0], 1)
# precompute terms that only depend on the histogramssets
self._deltas_up = default_backend.divide(
self._histogramssets[:, :, 2], self._histogramssets[:, :, 1]
)
self._deltas_dn = default_backend.divide(
self._histogramssets[:, :, 0], self._histogramssets[:, :, 1]
)
default_backend.shape(self._deltas_up)
)

deltas_up_alpha0 = default_backend.power(self._deltas_up, self._alpha0)
deltas_dn_alpha0 = default_backend.power(self._deltas_dn, self._alpha0)
# x = A^{-1} b
A_inverse = default_backend.astensor(
[
[
15.0 / (16 * alpha0),
-15.0 / (16 * alpha0),
-7.0 / 16.0,
-7.0 / 16.0,
1.0 / 16 * alpha0,
-1.0 / 16.0 * alpha0,
],
[
3.0 / (2 * math.pow(alpha0, 2)),
3.0 / (2 * math.pow(alpha0, 2)),
-9.0 / (16 * alpha0),
9.0 / (16 * alpha0),
1.0 / 16,
1.0 / 16,
],
[
-5.0 / (8 * math.pow(alpha0, 3)),
5.0 / (8 * math.pow(alpha0, 3)),
5.0 / (8 * math.pow(alpha0, 2)),
5.0 / (8 * math.pow(alpha0, 2)),
-1.0 / (8 * alpha0),
1.0 / (8 * alpha0),
],
[
3.0 / (-2 * math.pow(alpha0, 4)),
3.0 / (-2 * math.pow(alpha0, 4)),
-7.0 / (-8 * math.pow(alpha0, 3)),
7.0 / (-8 * math.pow(alpha0, 3)),
-1.0 / (8 * math.pow(alpha0, 2)),
-1.0 / (8 * math.pow(alpha0, 2)),
],
[
3.0 / (16 * math.pow(alpha0, 5)),
-3.0 / (16 * math.pow(alpha0, 5)),
-3.0 / (16 * math.pow(alpha0, 4)),
-3.0 / (16 * math.pow(alpha0, 4)),
1.0 / (16 * math.pow(alpha0, 3)),
-1.0 / (16 * math.pow(alpha0, 3)),
],
[
1.0 / (2 * math.pow(alpha0, 6)),
1.0 / (2 * math.pow(alpha0, 6)),
-5.0 / (16 * math.pow(alpha0, 5)),
5.0 / (16 * math.pow(alpha0, 5)),
1.0 / (16 * math.pow(alpha0, 4)),
1.0 / (16 * math.pow(alpha0, 4)),
],
]
)
b = default_backend.stack(
[
default_backend.log(self._deltas_up) * deltas_up_alpha0,
-default_backend.log(self._deltas_dn) * deltas_dn_alpha0,
default_backend.power(default_backend.log(self._deltas_up), 2)
* deltas_up_alpha0,
default_backend.power(default_backend.log(self._deltas_dn), 2)
* deltas_dn_alpha0,
]
)
self._coefficients = default_backend.einsum(
)

self._precompute()
if subscribe:
events.subscribe('tensorlib_changed')(self._precompute)

[docs]
def _precompute(self):
tensorlib, _ = get_backend()
self.deltas_up = tensorlib.astensor(self._deltas_up)
self.deltas_dn = tensorlib.astensor(self._deltas_dn)
self.alpha0 = tensorlib.astensor(self._alpha0)
self.coefficients = tensorlib.astensor(self._coefficients)
self.bases_up = tensorlib.einsum(
'sa,shb->shab', tensorlib.ones(self.alphasets_shape), self.deltas_up
)
self.bases_dn = tensorlib.einsum(
'sa,shb->shab', tensorlib.ones(self.alphasets_shape), self.deltas_dn
)
self.ones = tensorlib.einsum(
)

[docs]
def _precompute_alphasets(self, alphasets_shape):
if alphasets_shape == self.alphasets_shape:
return
tensorlib, _ = get_backend()
self.alphasets_shape = alphasets_shape
self.bases_up = tensorlib.einsum(
'sa,shb->shab', tensorlib.ones(self.alphasets_shape), self.deltas_up
)
self.bases_dn = tensorlib.einsum(
'sa,shb->shab', tensorlib.ones(self.alphasets_shape), self.deltas_dn
)
self.ones = tensorlib.einsum(
)
return

def __call__(self, alphasets):
"""Compute Interpolated Values."""
tensorlib, _ = get_backend()
self._precompute_alphasets(tensorlib.shape(alphasets))

# select where alpha >= alpha0 and produce the mask
where_alphasets_gtalpha0 = tensorlib.where(
)
tensorlib.einsum(
),
dtype="bool",
)

# select where alpha > -alpha0 ["not(alpha <= -alpha0)"] and produce the mask
where_alphasets_not_ltalpha0 = tensorlib.where(
)
tensorlib.einsum(
),
dtype="bool",
)

# s: set under consideration (i.e. the modifier)
# a: alpha variation
# h: histogram affected by modifier
# b: bin of histogram
exponents = tensorlib.einsum(
)
# for |alpha| >= alpha0, we want to raise the bases to the exponent=alpha
# and for |alpha| < alpha0, we want to raise the bases to the exponent=1
exponents >= self.__alpha0, exponents, self.ones
)
# we need to produce the terms of alpha^i for summing up
alphasets_powers = tensorlib.stack(
[
alphasets,
tensorlib.power(alphasets, 2),
tensorlib.power(alphasets, 3),
tensorlib.power(alphasets, 4),
tensorlib.power(alphasets, 5),
tensorlib.power(alphasets, 6),
]
)
# this is the 1 + sum_i a_i alpha^i
value_btwn = tensorlib.ones(exponents.shape) + tensorlib.einsum(
'rshb,rsa->shab', self.coefficients, alphasets_powers
)

# first, build a result where:
#       alpha > alpha0   : fill with bases_up
#   not(alpha > alpha0)  : fill with 1 + sum(a_i alpha^i)
results_gtalpha0_btwn = tensorlib.where(
)
# then, build a result where:
#      alpha >= -alpha0  : do nothing (fill with previous result)
#   not(alpha >= -alpha0): fill with bases_dn
bases = tensorlib.where(
)

class _slow_code4:
"""
Reference Implementation of Code 4.

delta_up^alpha0 = 1 + a1 alpha0 + a2 alpha0^2 + a3 alpha0^3 + a4 alpha0^4 + a5 alpha0^5 + a6 alpha0^6
delta_down^alpha0 = 1 - a1 alpha0 + a2 alpha0^2 - a3 alpha0^3 + a4 alpha0^4 - a5 alpha0^5 + a6 alpha0^6

f[alpha_] := 1 + a1 * alpha + a2 * alpha^2 + a3 * alpha^3 + a4 * alpha^4 + a5 * alpha^5 + a6 * alpha^6
up[alpha_] := delta_up^alpha
down[alpha_] := delta_down^(-alpha)

We want to find the coefficients a1, a2, a3, a4, a5, a6 by solving:
f[alpha0] == up[alpha0]
f[-alpha0] == down[-alpha0]
f'[alpha0] == up'[alpha0]
f'[-alpha0] == down'[-alpha0]
f''[alpha0] == up''[alpha0]
f''[-alpha0] == down''[-alpha0]

Treating this as multiplication with a rank-6 matrix A: A*x = b, where x = [a1, a2, a3, a4, a5, a6]

[alpha0,  alpha0^2, alpha0^3,   alpha0^4,   alpha0^5,     alpha0^6  ] [a1] = [ delta_up^(alpha0) - 1                ]
[-alpha0, alpha0^2, -alpha0^3,  alpha0^4,   -alpha0^5,    alpha0^6  ] [a2] = [ delta_down^(alpha0) - 1              ]
[1,       2alpha0,  3alpha0^2,  4alpha0^3,  5alpha0^4,    6alpha0^5 ] [a3] = [ ln(delta_up) delta_up^(alpha0)       ]
[1,       -2alpha0, 3alpha0^2,  -4alpha0^3, 5alpha0^4,    -6alpha0^5] [a4] = [ - ln(delta_down) delta_down^(alpha0) ]
[0,       2,        6alpha0,    12alpha0^2, 20alpha0^3,   30alpha0^4] [a5] = [ ln(delta_up)^2 delta_up^(alpha0)     ]
[0,       2,        -6alpha0,   12alpha0^2, -20alpha0^3,  30alpha0^4] [a6] = [ ln(delta_down)^2 delta_down^(alpha0) ]

The determinant of this matrix is -2048*alpha0^15. The trace is 30*alpha0^4+16*alpha0^3+4*alpha0^2+alpha0. Therefore, this matrix is invertible if and only if alpha0 != 0.

[15/(16*alpha0),  -15/(16*alpha0),  -7/16,            -7/16,            1/16*alpha0,      -1/16*alpha0    ]
[3/(2*alpha0^2),  3/(2*alpha0^2),   -9/(16*alpha0),   9/(16*alpha0),    1/16,             1/16            ]
[-5/(8*alpha0^3), 5/(8*alpha0^3),   5/(8*alpha0^2),   5/(8*alpha0^2),   -1/(8*alpha0),    1/(8*alpha0)    ]
[3/(-2*alpha0^4), 3/(-2*alpha0^4),  -7/(-8*alpha0^3), 7/(-8*alpha0^3),  -1/(8*alpha0^2),  -1/(8*alpha0^2) ]
[3/(16*alpha0^5), -3/(16*alpha0^5), -3/(16*alpha0^4), -3/(16*alpha0^4), 1/(16*alpha0^3),  -1/(16*alpha0^3)]
[1/(2*alpha0^6),  1/(2*alpha0^6),   -5/(16*alpha0^5), 5/(16*alpha0^5),  1/(16*alpha0^4),  1/(16*alpha0^4) ]

"""

def product(self, down, nom, up, alpha):
delta_up = up / nom
delta_down = down / nom
if alpha >= self.alpha0:
delta = math.pow(delta_up, alpha)
elif -self.alpha0 < alpha < self.alpha0:
delta_up_alpha0 = math.pow(delta_up, self.alpha0)
delta_down_alpha0 = math.pow(delta_down, self.alpha0)
b = [
delta_up_alpha0 - 1,
delta_down_alpha0 - 1,
math.log(delta_up) * delta_up_alpha0,
-math.log(delta_down) * delta_down_alpha0,
math.pow(math.log(delta_up), 2) * delta_up_alpha0,
math.pow(math.log(delta_down), 2) * delta_down_alpha0,
]
A_inverse = [
[
15.0 / (16 * self.alpha0),
-15.0 / (16 * self.alpha0),
-7.0 / 16.0,
-7.0 / 16.0,
1.0 / 16 * self.alpha0,
-1.0 / 16.0 * self.alpha0,
],
[
3.0 / (2 * math.pow(self.alpha0, 2)),
3.0 / (2 * math.pow(self.alpha0, 2)),
-9.0 / (16 * self.alpha0),
9.0 / (16 * self.alpha0),
1.0 / 16,
1.0 / 16,
],
[
-5.0 / (8 * math.pow(self.alpha0, 3)),
5.0 / (8 * math.pow(self.alpha0, 3)),
5.0 / (8 * math.pow(self.alpha0, 2)),
5.0 / (8 * math.pow(self.alpha0, 2)),
-1.0 / (8 * self.alpha0),
1.0 / (8 * self.alpha0),
],
[
3.0 / (-2 * math.pow(self.alpha0, 4)),
3.0 / (-2 * math.pow(self.alpha0, 4)),
-7.0 / (-8 * math.pow(self.alpha0, 3)),
7.0 / (-8 * math.pow(self.alpha0, 3)),
-1.0 / (8 * math.pow(self.alpha0, 2)),
-1.0 / (8 * math.pow(self.alpha0, 2)),
],
[
3.0 / (16 * math.pow(self.alpha0, 5)),
-3.0 / (16 * math.pow(self.alpha0, 5)),
-3.0 / (16 * math.pow(self.alpha0, 4)),
-3.0 / (16 * math.pow(self.alpha0, 4)),
1.0 / (16 * math.pow(self.alpha0, 3)),
-1.0 / (16 * math.pow(self.alpha0, 3)),
],
[
1.0 / (2 * math.pow(self.alpha0, 6)),
1.0 / (2 * math.pow(self.alpha0, 6)),
-5.0 / (16 * math.pow(self.alpha0, 5)),
5.0 / (16 * math.pow(self.alpha0, 5)),
1.0 / (16 * math.pow(self.alpha0, 4)),
1.0 / (16 * math.pow(self.alpha0, 4)),
],
]

coefficients = [
sum(A_i * b_j for A_i, b_j in zip(A_row, b)) for A_row in A_inverse
]
delta = 1
for i in range(1, 7):
delta += coefficients[i - 1] * math.pow(alpha, i)
else:
delta = math.pow(delta_down, (-alpha))
return delta

def __init__(self, histogramssets, subscribe=True, alpha0=1):
self._histogramssets = histogramssets
self.alpha0 = alpha0

def __call__(self, alphasets):
tensorlib, _ = get_backend()
return tensorlib.astensor(
_slow_interpolator_looper(
self._histogramssets, tensorlib.tolist(alphasets), self.product
)
)