diff --git a/dokimi/changelog.d/20261114_130000_sign-flip-interval-float-range.md b/dokimi/changelog.d/20261114_130000_sign-flip-interval-float-range.md new file mode 100644 index 000000000..1c44ccd0e --- /dev/null +++ b/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. diff --git a/dokimi/docs/benchmarking-and-comparison.md b/dokimi/docs/benchmarking-and-comparison.md index 9fa66e983..e0d023640 100644 --- a/dokimi/docs/benchmarking-and-comparison.md +++ b/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. diff --git a/dokimi/dokimi/eval/sign_flip.py b/dokimi/dokimi/eval/sign_flip.py index 0b00de55f..80601dd59 100644 --- a/dokimi/dokimi/eval/sign_flip.py +++ b/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,19 +125,21 @@ 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 @@ -140,9 +147,9 @@ def p(self, theta: float) -> float | None: 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 diff --git a/dokimi/tests/test_repo_compare_overflow.py b/dokimi/tests/test_repo_compare_overflow.py index c181b6252..646756fcf 100644 --- a/dokimi/tests/test_repo_compare_overflow.py +++ b/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(): diff --git a/dokimi/tests/test_sign_flip_float_range.py b/dokimi/tests/test_sign_flip_float_range.py new file mode 100644 index 000000000..5903c12ae --- /dev/null +++ b/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)