Compute float bin centres without overflowing near the float64 limits

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

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"

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

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:

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))