diff --git a/AGENTS.md b/AGENTS.md index 47f6dc1..cde7114 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -30,6 +30,7 @@ src/morphobricks_normalization/ compute_block_histogram for planes read in blocks), memory estimate; robust range with HistogramTail underflow/overflow and histogram_with_range_check, quantiles by radix selection (Milestone 5) + bin_centres.py float_bin_centres: bin centres that do not overflow near the float64 limits ranges.py HistogramSettings, RangeCheck, range_check_issues (Milestone 5) integrity.py invariants every PlaneHistogram and HistogramTail checks when made (counts, limits, tail sums) and the pixel type of every block diff --git a/pyproject.toml b/pyproject.toml index a936580..8c02e20 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -113,4 +113,5 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.scopes" = "pure" "morphobricks_normalization.settings" = "pure" "morphobricks_normalization.summaries" = "pure" +"morphobricks_normalization.bin_centres" = "pure" "morphobricks_normalization.sums" = "pure" diff --git a/src/morphobricks_normalization/bin_centres.py b/src/morphobricks_normalization/bin_centres.py new file mode 100644 index 0000000..b9b642a --- /dev/null +++ b/src/morphobricks_normalization/bin_centres.py @@ -0,0 +1,24 @@ +"""Centres of equal-width float bins, without overflow near the float64 limits.""" + +from __future__ import annotations + +import numpy as np + + +def float_bin_centres(lo: float, hi: float, bin_count: int) -> np.ndarray: + """Centre of each of ``bin_count`` equal bins between ``lo`` and ``hi`` (float64). + + Valid outputs are the plain arithmetic ``lo + half_range * fraction``. That product + reaches twice the half range (up to the full range) for the last bins, which + overflows near the float64 limits although the centre is representable; only then are + the affected bins stepped from the midpoint instead. + """ + half_range = hi * 0.5 - lo * 0.5 # halving keeps the range finite + fraction = (2 * np.arange(bin_count) + 1) / bin_count + with np.errstate(over="ignore"): # an overflowing step is redone just below + centres = lo + half_range * fraction + over = ~np.isfinite(centres) + if over.any(): + midpoint = lo * 0.5 + hi * 0.5 + centres[over] = midpoint + half_range * (fraction[over] - 1) + return centres diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py index 830c2fd..6510bec 100644 --- a/src/morphobricks_normalization/histogram.py +++ b/src/morphobricks_normalization/histogram.py @@ -50,6 +50,7 @@ import numpy as np +from .bin_centres import float_bin_centres from .dtypes import check_supported, default_bin_count, is_integer, validate_bin_count from .integrity import check_histogram, check_tail, checked_blocks from .ranges import CENTRAL_PERCENTILES, HistogramSettings, RangeCheck @@ -177,8 +178,7 @@ def representatives(self) -> np.ndarray: dtype=np.float64, ) lo, hi = float(self.binning_min), float(self.binning_max) - half_range = hi * 0.5 - lo * 0.5 - return lo + half_range * ((2 * np.arange(self.bin_count) + 1) / self.bin_count) + return float_bin_centres(lo, hi, self.bin_count) @property def max_representative_error(self) -> float: diff --git a/tests/test_pooling.py b/tests/test_pooling.py index 10ea918..154481d 100644 --- a/tests/test_pooling.py +++ b/tests/test_pooling.py @@ -116,3 +116,29 @@ def test_linear_percentile_between_finite_extremes_does_not_overflow(q, low, hig 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) + + +@pytest.mark.parametrize("q", [0, 25, 50, 75, 100]) +@pytest.mark.parametrize("limit", [np.finfo(np.float64).max, 1e308]) +def test_percentiles_from_pixels_near_float64_limits_are_finite_and_exact(q, limit): + # End to end (#75): pixels at -limit and +limit -> histogram -> pooling -> percentile. + # Bin centres reach twice the half range, which overflowed before the guard. + from fractions import Fraction + + from morphobricks_normalization.histogram import compute_plane_histogram + + pixels = np.array([[-limit, limit]], dtype=np.float64) + n = 4 + histogram = compute_plane_histogram(pixels, bin_count=n) + low, high = Fraction(-limit), Fraction(limit) + centres = [low + (high - low) * Fraction(2 * i + 1, 2 * n) for i in range(n)] + assert np.isfinite(histogram.representatives).all() + for got, want in zip(histogram.representatives, centres, strict=True): + assert abs(Fraction(float(got)) - want) <= abs(want) * Fraction(1, 10**14) + pooled = PooledHistogram([histogram]) + gamma = Fraction(q, 100) + expected = centres[0] + (centres[-1] - centres[0]) * gamma + result = pooled.linear_percentile(q) + assert np.isfinite(result) + assert abs(Fraction(result) - expected) <= abs(expected) * Fraction(1, 10**14) + assert np.isfinite(pooled.percentile(q))