interpolation code2 is deliberately discontinuous

[kratsg]

Summary

1. Discontinuity at $\alpha = \pm 1$: both fast and slow paths implement the documented I_quad|lin formula whose extrapolation term is zero at the boundary, value drops from 2.0 to 1.00125 crossing $\alpha = 1$. ROOT's flexibleInterpSingle case 2 (in RooFit/Detail/MathFuncs.h) adds the + $(I^+ − I^0) / + (I^- − I^0)$ offset and is continuous. This affects any minimizer exploring $|\alpha| > 1$. Present in src/pyhf/interpolators/code2.py on current main (the docstring's LaTeX formula matches the code, so it may be "as designed".
2. Fast/slow disagreement for $\alpha &lt; −1$ in v0.7.6 (last row of the table): fixed in fix: Correct mask in code2 fast interpolation for alpha < -1 #2668

---

(below is AI-assisted with Claude Fable)

pyhf's reference "slow" implementation,
src/pyhf/interpolators/code2.py::_slow_code2.summand:

class _slow_code2:
    def summand(self, down, nom, up, alpha):
        a = 0.5 * (up + down) - nom
        b = 0.5 * (up - down)
        if alpha > 1:
            delta = (b + 2 * a) * (alpha - 1)
        elif -1 <= alpha <= 1:
            delta = a * alpha * alpha + b * alpha
        else:
            delta = (b - 2 * a) * (alpha + 1)
        return delta

computes exactly the same unshifted extrapolation as pyhs3's interpolate_code2
— no + (up - nom) / + (down - nom) offset. pyhf's tensorized "fast"
code2.__call__ uses the identical unshifted formula
(value_gt1 = (alphasets - 1) * self.b_plus_2a, etc.).

Numerically (nom=1.0, hi=2.0, lo=0.5), pyhf's own summand gives:

alpha	value
0.999	1.99875025
1.0	2.0
1.001	1.0012499999999998
-0.999	0.50025025
-1.0	0.5
-1.001	0.99975

— a genuine jump crossing the boundary, not a documentation artifact.

By contrast, ROOT's FlexibleInterpVar/PiecewiseInterpolation code 2
(RooFit::Detail::MathFuncs::flexibleInterpSingle, case 2, in
roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h)
does include the offset:

} else if (code == 2) {
   // parabolic with linear
   double a = 0.5 * (high + low) - nominal;
   double b = 0.5 * (high - low);
   double c = 0;
   if (paramVal > 1) {
      return (2 * a + b) * (paramVal - 1) + high - nominal;
   } else if (paramVal < -1) {
      return -1 * (2 * a - b) * (paramVal + 1) + low - nominal;
   } else {
      return a * paramVal * paramVal + b * paramVal + c;
   }

Reference: scipp-atlas/pyhs3#259

OS / Environment

n/a

Steps to Reproduce

  """Reproduce code2 interpolator behavior at alpha = +-1 in pyhf.

  One systematic, one sample, one bin: (down, nominal, up) = (0.5, 1.0, 2.0).
  code2 is additive, so the interpolated value is nominal + delta.

  Finding 1 (present in 0.7.6 AND current main): the quadratic region gives
  delta(1.0) = a + b, but the extrapolation branches use
  delta = (b +- 2a)(alpha -+ 1), which is 0 at the boundary instead of a + b,
  so the interpolated value jumps as alpha crosses +-1. ROOT's
  flexibleInterpSingle (code 2) adds the offset and is continuous.

  Finding 2 (0.7.6 only, already fixed on main): the tensorized fast path
  builds the alpha < -1 term with `alphasets + self.mask_off` (zeros), i.e.
  (b - 2a) * alpha instead of (b - 2a) * (alpha + 1), so fast and slow
  disagree below alpha = -1. On main this line uses `+ self.mask_on`.
  """

  import numpy as np
  import pyhf
  from pyhf.interpolators.code2 import _slow_code2, code2

  histogramssets = [[[[0.5], [1.0], [2.0]]]]  # [sets][histos][down, nom, up][bins]
  nominal, down, up = 1.0, 0.5, 2.0

  a = 0.5 * (up + down) - nominal
  b = 0.5 * (up - down)

  fast = code2(histogramssets, subscribe=False)
  slow = _slow_code2(histogramssets, subscribe=False)


  def root_code2(alpha: float) -> float:
      """ROOT flexibleInterpSingle code 2: continuous at |alpha| = 1."""
      if alpha > 1:
          return nominal + (b + 2 * a) * (alpha - 1) + (a + b)
      if alpha < -1:
          return nominal + (b - 2 * a) * (alpha + 1) + (a - b)
      return nominal + a * alpha**2 + b * alpha


  print(f"pyhf {pyhf.__version__}")
  print(f"{'alpha':>8} {'pyhf fast':>22} {'pyhf slow':>22} {'ROOT code2':>22}")
  for alpha in [0.999, 1.0, 1.001, -0.999, -1.0, -1.001]:
      alphasets = np.array([[alpha]])
      delta_fast = np.asarray(fast(alphasets))[0][0][0][0]
      delta_slow = np.asarray(slow(alphasets))[0][0][0][0]
      print(
          f"{alpha:8.3f} {nominal + delta_fast:22.16f}"
          f" {nominal + delta_slow:22.16f} {root_code2(alpha):22.16f}"
      )

Output on pyhf 0.7.6:

     alpha              pyhf fast              pyhf slow             ROOT code2
     0.999     1.9987502500000001     1.9987502500000001     1.9987502500000001
     1.000     2.0000000000000000     2.0000000000000000     2.0000000000000000
     1.001     1.0012499999999998     1.0012499999999998     2.0012499999999998
    -0.999     0.5002502500000000     0.5002502500000000     0.5002502500000001
    -1.000     0.5000000000000000     0.5000000000000000     0.5000000000000000
    -1.001     0.7497500000000000     0.9997500000000000     0.4997500000000000

File Upload (optional)

No response

Expected Results

I expected pyhf to be continuous (matching ROOT) when we "extrapolate" beyond $|\alpha| &gt; 1$

Actual Results

pyhf is not continuous at the boundaries

pyhf Version

pyhf, version 0.7.6

Code of Conduct

• I agree to follow the Code of Conduct