fix(eval): the sign-flip interval of deltas near the float maximum scales instead of overflowing

dokimi/changelog.d/20261114_130000_sign-flip-interval-float-range.md

@@ -0,0 +1,3 @@
+### Fixed
+
+- **Sign-flip interval of deltas near the float maximum** — the interval search of `repo compare` and the arms track (`sign_flip_interval`) emitted `RuntimeWarning: invalid value encountered in matmul` for deltas such as `[1e308, 0.0]`, and a `nan` flowed through its comparisons so the ends could be silently wrong. A row whose largest delta is above `max float / (8 * n)` is now searched on the deltas divided by a power of two (exact) and each end is multiplied back, so it neither warns nor leaves the float range. An end beyond the float range is open (`None`), never `inf`. Rows below that threshold give the same ends as before.

dokimi/docs/benchmarking-and-comparison.md

@@ -403,7 +403,8 @@ count again, so both ends are open (`null`, printed `-inf` / `+inf`) at 5 or
 fewer pairs and bounded from 6 on. Six or more pairs that all tie get the
 interval [0, 0] with no p and a `—`: nothing speaks against 0, and every other
 shift moves all the pairs to one side, which the test rejects, as it does every
-shift but +20 for eight +20 pairs. `ci_method` is `sign_flip_inversion`; an
+shift but +20 for eight +20 pairs. An end the float range cannot hold (deltas
+near 1.8e308) is open too, never `inf`. `ci_method` is `sign_flip_inversion`; an
 `arms.json` written before this method holds `student_t`, `percentile_bootstrap`
 or `degenerate` there, and `p_ttest` (with `p_wilcoxon`, `n_wilcoxon`) instead of
 `p`, and the report still renders it.

dokimi/dokimi/eval/sign_flip.py

@@ -40,6 +40,7 @@
 import itertools
 import math
 import random
+import sys
 from typing import TYPE_CHECKING
 
 if TYPE_CHECKING:
@@ -70,12 +71,16 @@ def sign_flip_p(deltas: list[float], seed: int = 0, n_resamples: int = N_RESAMPL
     """Two-sided p of ``|sum(deltas)|`` when each delta's sign is a fair coin under the null.
 
     ``None`` when every delta is a tie: there is no evidence either way, and 1.0 would read as a measured null.
-    Rounded to 6 places as printed, unless ``exact``.
+    Rounded to 6 places as printed, unless ``exact``. Deltas whose largest magnitude is above ``max float / (8 * n)``
+    are divided by a power of two (exact) first, as in :func:`sign_flip_interval`, so no sum leaves the float range.
     """
     d = [x for x in deltas if abs(x) > EPS]
     if not d:
         return None
-    observed = abs(sum(d)) - EPS
+    eps, k = EPS, _scale_exponent(d)
+    if k:  # near the float maximum the sums go to inf/nan: the same p on the deltas divided by a power of two
+        d, eps = [math.ldexp(x, -k) for x in d], math.ldexp(EPS, -k)
+    observed = abs(sum(d)) - eps
     if len(d) <= EXACT_MAX:
         hits = sum(abs(sum(s * x for s, x in zip(signs, d, strict=True))) >= observed
                    for signs in itertools.product((1, -1), repeat=len(d)))
@@ -120,29 +125,31 @@ class ShiftTest:
     can round to one side of it in the product and to the other in ``sum()``, so ``p(0)`` can differ there: the
     row's own p, not ``p(0)``, places 0 (:func:`sign_flip_interval`)."""
 
-    def __init__(self, deltas: list[float], seed: int = 0, n_resamples: int = N_RESAMPLES) -> None:
+    def __init__(self, deltas: list[float], seed: int = 0, n_resamples: int = N_RESAMPLES, *,
+                 eps: float = EPS) -> None:
         import numpy as np
 
-        changed = [x for x in deltas if abs(x) > EPS]
-        self.d = np.array(changed + [x for x in deltas if abs(x) <= EPS], dtype=float)  # sign_flip_p's order first
-        self.changed, self.seed, self.n_resamples = len(changed), seed, n_resamples
+        changed = [x for x in deltas if abs(x) > eps]
+        self.d = np.array(changed + [x for x in deltas if abs(x) <= eps], dtype=float)  # sign_flip_p's order first
+        self.changed, self.seed, self.n_resamples, self.eps = len(changed), seed, n_resamples, eps
 
     def p(self, theta: float) -> float | None:
         """``None`` when every shifted delta ties, as :func:`sign_flip_p` has no p then."""
         import numpy as np
 
         x = self.d - theta
-        x[np.abs(x) <= EPS] = 0.0
+        eps = self.eps
+        x[np.abs(x) <= eps] = 0.0
         moved = x[x != 0.0]
         if not moved.size:
             return None
         observed = abs(math.fsum(moved.tolist()))  # exactly rounded: sum() rounds differently before 3.12
         if moved.size <= EXACT_MAX:
             sums = np.abs(_patterns(moved.size) @ moved)  # rows 0 and -1: the observed signs and their mirror
             least = min(observed, float(sums[0]), float(sums[-1]))  # the product rounds them otherwise from 1e6 on
-            return int(np.count_nonzero(sums >= least - EPS)) / 2 ** moved.size
+            return int(np.count_nonzero(sums >= least - eps)) / 2 ** moved.size
         signs = _draws(self.seed, self.changed, self.d.size, self.n_resamples)
-        return (int(np.count_nonzero(np.abs(signs @ x) >= observed - EPS)) + 1) / (self.n_resamples + 1)
+        return (int(np.count_nonzero(np.abs(signs @ x) >= observed - eps)) + 1) / (self.n_resamples + 1)
 
 
 def sign_flip_interval(deltas: list[float], p_zero: float | None, seed: int = 0, *,
@@ -155,23 +162,56 @@ def sign_flip_interval(deltas: list[float], p_zero: float | None, seed: int = 0,
 
     Every delta tied gives ``(0.0, 0.0)`` from 6 deltas on, and ``p_zero`` is ``None``: no evidence against 0, but
     any shift off 0 moves all n deltas to one side, which the test rejects at ``2 / 2**n``, as eight +20 deltas give
-    ``(20.0, 20.0)``. At 5 or fewer the ends are open, as at every row of so few."""
-    test, mean = ShiftTest(deltas, seed), math.fsum(deltas) / len(deltas)
+    ``(20.0, 20.0)``. At 5 or fewer the ends are open, as at every row of so few.
+
+    A largest ``|delta|`` above ``max float / (8 * n)`` would take the far bound, the shifted deltas or the pattern
+    sums out of the float range: the search then runs on the deltas divided by a power of two (exact), with
+    :data:`EPS` and :data:`TOL` divided the same, and each end is multiplied back. An end the float range cannot hold
+    is open (``None``), never ``inf``."""
+    k = _scale_exponent(deltas)
+    if k:
+        scaled = [math.ldexp(x, -k) for x in deltas]
+        ends = _ends(scaled, p_zero, seed, math.ldexp(EPS, -k), math.ldexp(TOL, -k))
+        ends = (_scaled_back(ends[0], k), _scaled_back(ends[1], k))
+    else:
+        ends = _ends(deltas, p_zero, seed, EPS, TOL)
+    return (_rounded(ends[0]), _rounded(ends[1])) if rounded else ends
+
+
+def _scale_exponent(deltas: list[float]) -> int:
+    """``k`` such that ``deltas / 2**k`` stays clear of the float maximum in the search; 0 when they already do."""
+    largest, limit = max(map(abs, deltas), default=0.0),sys.float_info.max / (8 * len(deltas))
+    if not math.isfinite(largest) or largest <= limit:
+        return 0
+    return math.frexp(largest)[1] - math.frexp(limit)[1] + 1
+
+
+def _scaled_back(end: float | None, k: int) -> float | None:
+    if end is None:
+        return None
+    try:
+        return math.ldexp(end, k)
+    except OverflowError:
+        return None
+
+
+def _ends(deltas: list[float], p_zero: float | None, seed: int, eps: float, tol: float
+          ) -> tuple[float | None, float | None]:
+    test, mean = ShiftTest(deltas, seed, eps=eps), math.fsum(deltas) / len(deltas)
     zero_inside, low, high = p_zero is None or p_zero >= ALPHA, min(*deltas, 0.0), max(*deltas, 0.0)
     # far is 1 past the extreme, or twice it from 1 on: 1e16 + 1.0 rounds back onto 1e16, all tied
-    ends = (_end(test, mean, low - max(1.0, -low), zero_inside),
-            _end(test, mean, high + max(1.0, high), zero_inside))
-    return (_rounded(ends[0]), _rounded(ends[1])) if rounded else ends
+    return (_end(test, mean, low - max(1.0, -low), zero_inside, tol),
+            _end(test, mean, high + max(1.0, high), zero_inside, tol))
 
 
-def _end(test: ShiftTest, mean: float, far: float, zero_inside: bool) -> float | None:
+def _end(test: ShiftTest, mean: float, far: float, zero_inside: bool, tol: float = TOL) -> float | None:
     """The end between ``mean`` and ``far``: the last shift bisection finds not rejected; ``None`` if ``far`` is not."""
     inside, outside = mean, far
     if min(mean, far) < 0 < max(mean, far):  # 0 is on this side: the row's p places it
         inside, outside = (0.0, far) if zero_inside else (mean, 0.0)
     if outside == far and _accepts(test.p(far)):
         return None
-    while abs(outside - inside) > TOL:
+    while abs(outside - inside) > tol:
         mid = (inside + outside) / 2
         if mid in (inside, outside):  # no float left between them
             break

dokimi/tests/test_repo_compare_overflow.py

@@ -71,9 +71,12 @@ def test_sign_flip_summary_refuses_finite_deltas_whose_sum_or_spread_leaves_the_
         sign_flip_summary([_task(d) for d in deltas])
 
 
-def test_sign_flip_summary_refuses_deltas_the_test_cannot_shift_within_the_range():
-    with pytest.raises(CompareError, match="too large to test"):
-        sign_flip_summary([_task(1.0e307)] * 6 + [_task(-1.0e307)] * 6)
+@pytest.mark.filterwarnings("ignore:overflow encountered:RuntimeWarning")  # the mde of a spread this large
+def test_sign_flip_summary_shifts_deltas_near_the_range_on_a_scaled_copy():
+    row = sign_flip_summary([_task(1.0e307), _task(-1.0e307)] * 6)  # cancels pairwise: sum() before 3.12 does not compensate
+    assert row["delta"]["mean"] == 0.0
+    assert math.isfinite(row["delta"]["lo"]) and math.isfinite(row["delta"]["hi"])
+    assert row["delta"]["lo"] < 0 < row["delta"]["hi"]
 
 
 def test_checked_mean_is_the_plain_mean_where_it_is_finite():

dokimi/tests/test_sign_flip_float_range.py

@@ -0,0 +1,87 @@
+# Copyright (C) 2026 Intel Corporation
+# SPDX-License-Identifier: Apache-2.0
+
+"""The sign-flip interval of deltas near the float maximum: no RuntimeWarning, no nan-driven end (#725)."""
+
+from __future__ import annotations
+
+import math
+import sys
+import warnings
+
+import pytest
+
+from dokimi.eval.sign_flip import sign_flip_interval, sign_flip_p, sign_flip_row
+from dokimi.eval.stats import _describe_deltas
+
+BIG = sys.float_info.max
+SHIFT = 25  # a power of two: scaling by it is exact, so a scaled row has the scaled row's ends
+MC_BASE = [1e299 * (0.5 + (i % 7) / 14) * (1 if i % 3 else -1) for i in range(20)]  # > EXACT_MAX deltas
+SIX_BASE = [1e300 * f for f in (0.6, 0.9, -0.3, 0.5, 0.7, 0.8)]
+ROWS = {
+    "1e308 and 0": [1e308, 0.0],
+    "1.7e308 and 0": [1.7e308, 0.0],
+    "max and 0": [BIG, 0.0],
+    "max and -max": [BIG, -BIG],
+    "six deltas": [math.ldexp(x, SHIFT) for x in SIX_BASE],
+    "monte carlo": [math.ldexp(x, SHIFT) for x in MC_BASE],
+}
+
+
+@pytest.mark.parametrize("name", ROWS)
+def test_a_row_near_the_maximum_returns_a_dict_and_warns_nothing(name):
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        row = sign_flip_row(ROWS[name], exact=False)
+    assert isinstance(row, dict)
+    for end in (row["lo"], row["hi"]):
+        assert end is None or math.isfinite(end)
+
+
+@pytest.mark.parametrize("name", [n for n in ROWS if n != "max and -max"])  # its magnitudes sum past the range
+def test_describe_deltas_of_a_row_near_the_maximum_warns_nothing(name):
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        described = _describe_deltas(ROWS[name], confidence=0.95, seed=0)
+    assert described["ci_method"] == "sign_flip_inversion"
+    assert math.isfinite(described["delta_mean"])
+    for end in (described["ci_lo"], described["ci_hi"]):
+        assert end is None or math.isfinite(end)
+
+
+@pytest.mark.parametrize("base", [SIX_BASE, MC_BASE], ids=["exact", "monte-carlo"])
+def test_the_ends_of_a_scaled_row_are_the_ends_of_the_row_scaled(base):
+    big = [math.ldexp(x, SHIFT) for x in base]
+    small = sign_flip_interval(base, sign_flip_p(base), rounded=False)
+    large = sign_flip_interval(big, sign_flip_p(big), rounded=False)
+    assert all(end is not None for end in small + large)
+    for got, want in zip(large, small, strict=True):
+        assert got == pytest.approx(math.ldexp(want, SHIFT), rel=1e-12)
+
+
+@pytest.mark.parametrize("row", [[BIG, 1.0, 1.0, 1.0, 1.0, 1.0], [-BIG, -1.0, -1.0, -1.0, -1.0, -1.0]])
+def test_an_end_beyond_the_float_range_is_open(row):
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        for rounded in (False, True):
+            ends = sign_flip_interval(row, sign_flip_p(row), rounded=rounded)
+            assert (ends[1] if row[0] > 0 else ends[0]) is None
+            assert all(end is None or math.isfinite(end) for end in ends)
+
+
+@pytest.mark.parametrize("row", [[1.7e308] * 3 + [1.5e308] * 3, [BIG, 1.0, 1.0, 1.0, 1.0, 1.0],
+                                 [math.ldexp(x, SHIFT) for x in MC_BASE]], ids=["six", "max and ones", "monte-carlo"])
+def test_the_p_and_ends_of_a_row_near_the_maximum_are_those_of_the_row_scaled(row):
+    small = [math.ldexp(x, -4) for x in row]
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        big_row, small_row = sign_flip_row(row, exact=False), sign_flip_row(small, exact=False)
+    assert big_row["p"] == small_row["p"]
+    for end in ("lo", "hi"):
+        try:
+            want = None if small_row[end] is None else math.ldexp(small_row[end], 4)
+        except OverflowError:
+            want = None  # the scaled end is past the float range: open
+        assert (big_row[end] is None) == (want is None)
+        if want is not None:
+            assert big_row[end] == pytest.approx(want, rel=1e-12)