diff --git a/src/morphobricks_normalization/pooling.py b/src/morphobricks_normalization/pooling.py index 32d3323..10abd90 100644 --- a/src/morphobricks_normalization/pooling.py +++ b/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: diff --git a/tests/test_pooling.py b/tests/test_pooling.py index b825989..10ea918 100644 --- a/tests/test_pooling.py +++ b/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) diff --git a/tools/code_rules_baseline.json b/tools/code_rules_baseline.json index 8c20b47..b19c9a6 100644 --- a/tools/code_rules_baseline.json +++ b/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,