Skip to content

refactor: deduplicate interpolator precompute and code4 matrix construction #2713

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

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

refactor: deduplicate interpolator precompute and code4 matrix construction #2713

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
24 changes: 10 additions & 14 deletions src/pyhf/interpolators/code1.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,9 @@ def __init__(self, histogramssets, subscribe=True):
if subscribe:
events.subscribe("tensorlib_changed")(self._precompute)

def _precompute(self):
def _precompute_alphaset_dependent(self):
# terms that depend on the alphasets_shape, recomputed when it changes
tensorlib, _ = get_backend()
self.deltas_up = tensorlib.astensor(self._deltas_up)
self.deltas_dn = tensorlib.astensor(self._deltas_dn)
self.broadcast_helper = tensorlib.astensor(self._broadcast_helper)
self.bases_up = tensorlib.einsum(
"sa,shb->shab", tensorlib.ones(self.alphasets_shape), self.deltas_up
)
Expand All @@ -63,20 +61,18 @@ def _precompute(self):
self.mask_on = tensorlib.ones(self.alphasets_shape)
self.mask_off = tensorlib.zeros(self.alphasets_shape)

def _precompute(self):
tensorlib, _ = get_backend()
self.deltas_up = tensorlib.astensor(self._deltas_up)
self.deltas_dn = tensorlib.astensor(self._deltas_dn)
self.broadcast_helper = tensorlib.astensor(self._broadcast_helper)
self._precompute_alphaset_dependent()

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.mask_on = tensorlib.ones(self.alphasets_shape)
self.mask_off = tensorlib.zeros(self.alphasets_shape)
return
self._precompute_alphaset_dependent()

def __call__(self, alphasets):
"""Compute Interpolated Values."""
Expand Down
198 changes: 75 additions & 123 deletions src/pyhf/interpolators/code4.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,64 @@
log = logging.getLogger(__name__)


def _a_inverse(alpha0):
"""Return the inverse of the rank-6 boundary-condition matrix as a nested list.

The matrix inverts ``A x = b`` for the polynomial coefficients (see
:class:`_slow_code4` for the derivation). It depends only on ``alpha0``.
"""
return [
[
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)),
],
]


class code4:
r"""
The polynomial interpolation and exponential extrapolation strategy.
Expand Down Expand Up @@ -60,58 +118,7 @@ def __init__(self, histogramssets, subscribe=True, alpha0=1):
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)),
],
]
)
A_inverse = default_backend.astensor(_a_inverse(alpha0))
b = default_backend.stack(
[
deltas_up_alpha0 - self._broadcast_helper,
Expand All @@ -132,13 +139,9 @@ def __init__(self, histogramssets, subscribe=True, alpha0=1):
if subscribe:
events.subscribe("tensorlib_changed")(self._precompute)

def _precompute(self):
def _precompute_alphaset_dependent(self):
# terms that depend on the alphasets_shape, recomputed when it changes
tensorlib, _ = get_backend()
self.deltas_up = tensorlib.astensor(self._deltas_up)
self.deltas_dn = tensorlib.astensor(self._deltas_dn)
self.broadcast_helper = tensorlib.astensor(self._broadcast_helper)
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
)
Expand All @@ -151,23 +154,20 @@ def _precompute(self):
"sa,shb->shab", self.mask_on, self.broadcast_helper
)

def _precompute(self):
tensorlib, _ = get_backend()
self.deltas_up = tensorlib.astensor(self._deltas_up)
self.deltas_dn = tensorlib.astensor(self._deltas_dn)
self.broadcast_helper = tensorlib.astensor(self._broadcast_helper)
self.alpha0 = tensorlib.astensor(self._alpha0)
self.coefficients = tensorlib.astensor(self._coefficients)
self._precompute_alphaset_dependent()

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.mask_on = tensorlib.ones(self.alphasets_shape)
self.mask_off = tensorlib.zeros(self.alphasets_shape)
self.ones = tensorlib.einsum(
"sa,shb->shab", self.mask_on, self.broadcast_helper
)
return
self._precompute_alphaset_dependent()

def __call__(self, alphasets):
"""Compute Interpolated Values."""
Expand Down Expand Up @@ -296,59 +296,9 @@ def product(self, down, nom, up, alpha):
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
sum(A_i * b_j for A_i, b_j in zip(A_row, b))
for A_row in self._A_inverse
]
delta = 1
for i in range(1, 7):
Expand All @@ -360,6 +310,8 @@ def product(self, down, nom, up, alpha):
def __init__(self, histogramssets, alpha0=1, subscribe=True): # noqa: ARG002
self._histogramssets = histogramssets
self.alpha0 = alpha0
# alpha0 is fixed, so the inverse matrix can be built once
self._A_inverse = _a_inverse(alpha0)

def __call__(self, alphasets):
tensorlib, _ = get_backend()
Expand Down
Loading