Interpolate percentiles between finite extremes without overflowing the gap

src/morphobricks_normalization/pooling.py

@@ -419,6 +419,22 @@ def _inverse_cdf_index(cumulative: np.ndarray, total: int, q: float) -> int:
     return int(np.searchsorted(cumulative, rank, side="left"))
 
 
+def _lerp(below: float, above: float, gamma: float) -> float:
+    """NumPy's ``_lerp``: start from the nearer endpoint; finite endpoints never overflow."""
+    with np.errstate(over="ignore"):  # an overflowing gap is handled just below
+        difference = above - below
+    if not math.isfinite(difference) and math.isfinite(below) and math.isfinite(above):
+        # The weighted mean of two finite values is representable; form it from the
+        # halved endpoints so the gap cannot overflow.
+        half = above * 0.5 - below * 0.5
+        if gamma >= 0.5:
+            return float(above - half * (2 * (1 - gamma)))
+        return float(below + half * (2 * gamma))
+    if gamma >= 0.5:
+        return float(above - difference * (1 - gamma))
+    return float(below + difference * gamma)
+
+
 def _linear_percentile(values: np.ndarray, cumulative: np.ndarray, total: float, q: float) -> float:
     # NumPy's "linear" method (Hyndman & Fan type 7), step by step as
     # numpy.lib's _quantile and _lerp compute it, so the result is bit-identical
@@ -436,10 +452,7 @@ def _linear_percentile(values: np.ndarray, cumulative: np.ndarray, total: float,
     below, above = (
         values[min(int(np.searchsorted(cumulative, i, side="right")), last)] for i in (k, k + 1)
     )
-    difference = above - below
-    if gamma >= 0.5:
-        return float(above - difference * (1 - gamma))
-    return float(below + difference * gamma)
+    return _lerp(below, above, gamma)
 
 
 def _weighted_index(cumulative: np.ndarray, total: float, q: float) -> int:

tests/test_pooling.py

@@ -8,7 +8,7 @@
     group_records,
     pooled_statistics,
 )
-from morphobricks_normalization.pooling import DEFAULT_PERCENTILES
+from morphobricks_normalization.pooling import DEFAULT_PERCENTILES, _linear_percentile
 from morphobricks_normalization.ranges import HistogramSettings
 from simulation import direct_stats, make_records, simulate_experiment
 
@@ -100,3 +100,19 @@ def test_opposing_infinite_tail_sums_are_a_named_refusal():
     with pytest.raises(ValueError, match="opposite sign") as caught:
         _ = PooledHistogram([histogram]).mean
     assert "tail 1 of 2" in str(caught.value) and "tail 2 of 2" in str(caught.value)
+
+
+@pytest.mark.parametrize("q", [0, 25, 50, 75, 100])
+@pytest.mark.parametrize("low, high", [(-1e308, 1e308), (-1.7e308, 1.7e308), (-1.7e308, 5e307)])
+def test_linear_percentile_between_finite_extremes_does_not_overflow(q, low, high):
+    # The gap between the endpoints exceeds float64 so above - below is infinite, yet
+    # the interpolation lies between them and is representable: it must be finite and
+    # equal the exact value (computed with fractions) to float64 precision.
+    from fractions import Fraction
+
+    values = np.array([low, high])
+    cumulative = np.array([1, 2])
+    got = _linear_percentile(values, cumulative, 2, q)
+    exact = Fraction(low) + (Fraction(high) - Fraction(low)) * Fraction(q, 100)
+    assert np.isfinite(got)
+    assert abs(Fraction(got) - exact) <= abs(exact) * Fraction(1, 10**15) + Fraction(1, 10**300)

tools/code_rules_baseline.json

@@ -34,7 +34,6 @@
     "src/morphobricks_normalization/histogram.py _integer_bins: private function of # lines without a docstring": 25,
     "src/morphobricks_normalization/histogram.py: # lines (max #)": 864,
     "src/morphobricks_normalization/pooling.py __init__: # Ruff findings missing-return-type-special-method": 1,
-    "src/morphobricks_normalization/pooling.py _linear_percentile: private function of # lines without a docstring": 13,
     "src/morphobricks_normalization/records.py _check_channels: private function of # lines without a docstring": 19,
     "src/morphobricks_normalization/scopes.py _incomplete_repeat_groups: private function of # lines without a docstring": 26,
     "src/morphobricks_normalization/settings.py _percentile_pair: # Ruff findings any-type": 1,