Document and test ultra-narrow float bins (no output change)

AGENTS.md

@@ -32,6 +32,8 @@ src/morphobricks_normalization/
                     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
+  bin_reachability.py  float_bin_indices, reachable_float_bins, centres_map_back: which float
+                    bins a value can fall in and whether centres map back (#86)
   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

CHANGELOG.md

@@ -14,6 +14,7 @@ new refusal or a numerical-range fix (see the pull requests for the exact messag
 - Percentiles are interpolated between finite extremes without overflowing the gap (#74).
 - Float bin centres are computed without overflowing near the float64 limits (#76).
 - Float histograms over a very narrow range (down to the smallest subnormal) bin pixels and place bin centres correctly: the position is a bounded fraction times the bin count where the halved arithmetic would round to zero or overflow. Histograms of ordinary ranges are bit-identical (#84).
+- Ultra-narrow float ranges (fewer representable floats than bins, e.g. `[0, 5e-324]` with 4 bins): documented and tested, no output changes. Bins no float64 value can fall in may share a representative centre; the centre of every bin holding two or more representable values maps back to it, and a bin holding exactly one can have a centre that rounds into the adjacent bin. New `bin_reachability` helpers state this; `tools/dump_outputs.py --corpus` dumps a generated corpus (all pixel types, extremes, subnormals, one-ulp to full-range float ranges, 2 to 1024 bins) that is bit-identical against `main` (#86).
 - Image summaries whose channels disagree, pooled summaries with different percentile
   definitions, and images whose Z/T planes have different pixel types are refused (#77).
 - RGB: channel count and pixel type are checked before histograms are built, RGB groups with

docs/contract.md

@@ -136,6 +136,13 @@ and ±infinity are ignored.
 - No finite pixels: `lo`, `hi` blank; all counts 0.
 - Representative value: the bin midpoint `lo + (i + 0.5)·width`.
 - Bound: every pixel is within `width / 2` of its representative.
+- Ultra-narrow ranges (#86): when `[lo, hi]` holds fewer representable floats than bins
+  (`[0, 5e-324]` with 4 bins), bins that no value can fall in may share a representative
+  (`0, 0, 5e-324, 5e-324`); they stay empty and the bin count is not reduced. Every bin that
+  holds two or more representable values has a representative that falls in that bin. A bin
+  holding exactly one can have a representative that rounds to the neighbouring float, in the
+  adjacent bin, because its midpoint is not representable. Bins and representatives are
+  unchanged by this statement.
 
 ##### Planes read in blocks
 

pyproject.toml

@@ -114,4 +114,5 @@ baseline = "tools/code_rules_baseline.json"
 "morphobricks_normalization.settings" = "pure"
 "morphobricks_normalization.summaries" = "pure"
 "morphobricks_normalization.bin_centres" = "pure"
+"morphobricks_normalization.bin_reachability" = "pure"
 "morphobricks_normalization.sums" = "pure"

src/morphobricks_normalization/bin_centres.py

@@ -44,6 +44,13 @@ def float_bin_centres(lo: float, hi: float, bin_count: int) -> np.ndarray:
     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.
+
+    A range that holds fewer representable floats than ``bin_count`` has bins no value can
+    fall in (``[0, 5e-324]`` with 4 bins: ``0, 0, 5e-324, 5e-324``). Such unreachable bins may
+    share a centre. For every bin that some float64 value falls in, the centre maps back to
+    that bin, except that a bin holding exactly one float can have a centre that rounds to the
+    neighbouring float, in the adjacent bin (its midpoint is not representable). The binning
+    and the centres do not depend on this; ``bin_reachability`` states and tests it.
     """
     fraction = (2 * np.arange(bin_count) + 1) / bin_count
     span = hi - lo

src/morphobricks_normalization/bin_reachability.py

@@ -0,0 +1,113 @@
+"""Which float bins a representable value can fall in, and whether centres map back to them.
+
+When ``[lo, hi]`` holds fewer float64 values than there are bins (for example
+``[0, 5e-324]`` with 4 bins), some bins can hold no value at all. Those bins may share a
+representative centre; the invariant is for the others. This module states it as code:
+
+- :func:`float_bin_indices` is the binning of the histograms, as a function;
+- :func:`reachable_float_bins` marks the bins some float64 value in the range falls in;
+- :func:`centres_map_back` marks the bins whose centre falls in that same bin.
+
+Every reachable bin must be one whose centre maps back to it. Nothing here changes how
+pixels are binned or where centres are placed.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from .bin_centres import float_bin_centres, float_bin_position
+
+_SIGN_BIT = np.uint64(1) << np.uint64(63)
+_MAGNITUDE = ~_SIGN_BIT
+
+
+def float_bin_indices(values: np.ndarray | float, lo: float, hi: float, n: int) -> np.ndarray:
+    """Bin of each value, ``n`` bins over ``[lo, hi]``, as the histograms bin pixels.
+
+    Parameters
+    ----------
+    values
+        Finite values; values outside the range are clipped to the first or last bin.
+    lo, hi
+        The binning range, ``lo <= hi``. A constant range puts every value in bin 0.
+    n
+        Number of bins.
+
+    Returns
+    -------
+    numpy.ndarray
+        ``int64`` bin index of every value.
+    """
+    values = np.asarray(values, dtype=np.float64)
+    if lo == hi:
+        return np.zeros(values.shape, dtype=np.int64)
+    index = np.floor(float_bin_position(values, lo, hi, n)).astype(np.int64)
+    return np.clip(index, 0, n - 1)
+
+
+def reachable_float_bins(lo: float, hi: float, n: int) -> np.ndarray:
+    """Mark the bins that some float64 value in ``[lo, hi]`` falls in.
+
+    The bin position is non-decreasing in the value, so the first value that reaches bin
+    ``b`` is found by bisection over the ordered float64 values, and bin ``b`` is reachable
+    when that value falls in ``b`` itself and not beyond it. This is exact and takes about
+    64 vectorised steps for any range.
+
+    Parameters
+    ----------
+    lo, hi
+        The binning range, ``lo <= hi``.
+    n
+        Number of bins.
+
+    Returns
+    -------
+    numpy.ndarray
+        Boolean array of length ``n``, True for a bin that can hold a value.
+    """
+    wanted = np.arange(n, dtype=np.int64)
+    left = np.full(n, _ordered_key(lo), dtype=np.int64)
+    right = np.full(n, _ordered_key(hi), dtype=np.int64)
+    while True:
+        active = left < right
+        if not active.any():
+            break
+        middle = (left >> 1) + (right >> 1) + (left & right & 1)  # floor mean, no overflow
+        reaches = float_bin_indices(_from_ordered_key(middle), lo, hi, n) >= wanted
+        right = np.where(active & reaches, middle, right)
+        left = np.where(active & ~reaches, middle + 1, left)
+    return float_bin_indices(_from_ordered_key(left), lo, hi, n) == wanted
+
+
+def centres_map_back(lo: float, hi: float, n: int) -> np.ndarray:
+    """Mark the bins whose centre falls in that same bin.
+
+    Parameters
+    ----------
+    lo, hi
+        The binning range, ``lo <= hi``.
+    n
+        Number of bins.
+
+    Returns
+    -------
+    numpy.ndarray
+        Boolean array of length ``n``.
+    """
+    centres = float_bin_centres(lo, hi, n)
+    return float_bin_indices(centres, lo, hi, n) == np.arange(n)
+
+
+def _ordered_key(value: float) -> np.int64:
+    """Integer that grows with the float (``-0.0`` and ``0.0`` share key 0)."""
+    bits = np.array(value, dtype=np.float64).view(np.uint64)
+    magnitude = np.int64(bits & _MAGNITUDE)
+    return -magnitude if bits & _SIGN_BIT else magnitude
+
+
+def _from_ordered_key(keys: np.ndarray) -> np.ndarray:
+    """Float64 values of integer keys made by :func:`_ordered_key`."""
+    negative = keys < 0
+    bits = np.abs(keys).astype(np.uint64) | np.where(negative, _SIGN_BIT, np.uint64(0))
+    return bits.view(np.float64)

tests/test_bin_reachability.py

@@ -0,0 +1,172 @@
+"""Ultra-narrow float ranges (#86): bins no float64 value can fall in, and the centres' invariant.
+
+The invariant, checked against every float64 value of narrow ranges (enumerated with
+``nextafter``) and against the histograms of those values:
+
+- a bin that holds no representable value ("unreachable") may share its centre with another;
+- a bin that holds two or more representable values has a centre that maps back to it;
+- a bin that holds exactly one may have a centre that rounds to the neighbouring float, which
+  is in the adjacent bin (the midpoint of a one-float bin is not representable).
+"""
+
+import numpy as np
+import pytest
+
+from morphobricks_normalization import compute_plane_histogram
+from morphobricks_normalization.bin_centres import float_bin_centres, float_bin_position
+from morphobricks_normalization.bin_reachability import (
+    centres_map_back,
+    float_bin_indices,
+    reachable_float_bins,
+)
+
+SMALLEST_SUBNORMAL = 5e-324
+SMALLEST_NORMAL = float(np.finfo(np.float64).tiny)
+MAX = float(np.finfo(np.float64).max)
+BIN_COUNTS = (1, 2, 3, 4, 5, 7, 16, 17, 64, 256)
+STARTS = (
+    0.0,
+    -0.0,
+    SMALLEST_SUBNORMAL,
+    -SMALLEST_SUBNORMAL,
+    1e-320,
+    SMALLEST_NORMAL,
+    1e-300,
+    0.1,
+    1.0,
+    -2.0,
+    3.7,
+    1e10,
+    1e300,
+    -1e-310,
+    -1e300,
+    1e308,
+    -MAX,
+    float(np.nextafter(MAX, 0.0)),
+)
+STEPS = (1, 2, 3, 4, 5, 6, 7, 8, 11, 20, 33, 100)
+
+
+def representable(lo: float, steps: int) -> np.ndarray:
+    """``steps + 1`` consecutive float64 values from ``lo`` (every value of the range)."""
+    values = [lo]
+    for _ in range(steps):
+        with np.errstate(over="ignore"):  # stepping past the maximum gives infinity
+            values.append(float(np.nextafter(values[-1], np.inf)))
+    return np.array(values)
+
+
+def narrow_ranges() -> list[tuple[float, int]]:
+    cases = []
+    for start in STARTS:
+        for steps in STEPS:
+            if np.isfinite(representable(start, steps)).all():
+                cases.append((start, steps))
+    return cases
+
+
+def test_issue_example_smallest_subnormal_with_four_bins():
+    lo, hi = 0.0, SMALLEST_SUBNORMAL
+    assert reachable_float_bins(lo, hi, 4).tolist() == [True, False, False, True]
+    centres = float_bin_centres(lo, hi, 4)
+    assert centres.tolist() == [0.0, 0.0, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL]
+    # The unreachable inner bins share a centre with a reachable one; the reachable map back.
+    assert centres_map_back(lo, hi, 4).tolist() == [True, False, False, True]
+    assert float_bin_indices(np.array([lo, hi]), lo, hi, 4).tolist() == [0, 3]
+    plane = np.array([0.0, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL])
+    hist = compute_plane_histogram(plane, 4)
+    assert hist.counts.tolist() == [1, 0, 0, 2]
+    assert hist.representatives.tolist() == centres.tolist()
+
+
+@pytest.mark.parametrize(("lo", "steps"), narrow_ranges())
+def test_reachable_bins_are_exactly_those_holding_a_representable_value(lo, steps):
+    values = representable(lo, steps)
+    hi = float(values[-1])
+    for n in BIN_COUNTS:
+        counts = np.bincount(float_bin_indices(values, lo, hi, n), minlength=n)
+        np.testing.assert_array_equal(reachable_float_bins(lo, hi, n), counts > 0)
+
+
+@pytest.mark.parametrize(("lo", "steps"), narrow_ranges())
+def test_centre_of_every_bin_with_two_values_maps_back(lo, steps):
+    values = representable(lo, steps)
+    hi = float(values[-1])
+    for n in BIN_COUNTS:
+        counts = np.bincount(float_bin_indices(values, lo, hi, n), minlength=n)
+        back = centres_map_back(lo, hi, n)
+        assert back[counts >= 2].all(), (lo, hi, n)
+        # Every value lands in a bin whose centre maps back, unless it is that bin's only value.
+        landed = float_bin_indices(values, lo, hi, n)
+        assert (back[landed] | (counts[landed] == 1)).all()
+
+
+@pytest.mark.parametrize(("lo", "steps"), narrow_ranges())
+def test_centre_of_a_one_value_bin_falls_in_an_adjacent_bin(lo, steps):
+    values = representable(lo, steps)
+    hi = float(values[-1])
+    for n in BIN_COUNTS:
+        counts = np.bincount(float_bin_indices(values, lo, hi, n), minlength=n)
+        where = float_bin_indices(float_bin_centres(lo, hi, n), lo, hi, n)
+        single = counts == 1
+        assert (abs(where[single] - np.arange(n)[single]) <= 1).all(), (lo, hi, n)
+
+
+def test_a_bin_with_one_value_can_have_a_centre_in_the_next_bin():
+    # [0, 3 * smallest subnormal], 3 bins: bin 1 holds only 1 * 5e-324, its midpoint 1.5 *
+    # 5e-324 rounds to 2 * 5e-324, which is in bin 2. Documented, not changed (#86).
+    hi = 3 * SMALLEST_SUBNORMAL
+    assert reachable_float_bins(0.0, hi, 3).all()
+    assert centres_map_back(0.0, hi, 3).tolist() == [True, False, True]
+    assert float_bin_centres(0.0, hi, 3).tolist() == [0.0, 2 * SMALLEST_SUBNORMAL, hi]
+
+
+@pytest.mark.parametrize(("lo", "steps"), narrow_ranges())
+def test_histograms_of_every_value_leave_exactly_the_unreachable_bins_empty(lo, steps):
+    values = representable(lo, steps)
+    hi = float(values[-1])
+    for n in (2, 4, 7, 64):
+        hist = compute_plane_histogram(values, n)
+        assert (hist.binning_min, hist.binning_max) == (lo + 0.0, hi)
+        reachable = reachable_float_bins(hist.binning_min, hist.binning_max, n)
+        np.testing.assert_array_equal(hist.counts > 0, reachable)
+        assert hist.counts.sum() == values.size
+        centres_back = centres_map_back(hist.binning_min, hist.binning_max, n)
+        np.testing.assert_array_equal(hist.representatives, float_bin_centres(lo, hi, n))
+        assert centres_back[hist.counts >= 2].all()
+
+
+def test_wide_ranges_map_every_centre_back():
+    for lo, hi in ((-MAX, MAX), (0.0, MAX), (-1e308, 1e308), (-1.0, 1.0), (0.0, 1e-306)):
+        for n in (2, 4, 16, 256, 1024):
+            assert reachable_float_bins(lo, hi, n).all()
+            assert centres_map_back(lo, hi, n).all()
+
+
+def test_random_ranges_of_every_magnitude_map_every_reachable_centre_back():
+    rng = np.random.default_rng(86)
+    checked = 0
+    while checked < 3000:
+        exponents = rng.integers(-1074, 1024, 2).astype(float)
+        signs = rng.choice([-1.0, 1.0], 2)
+        with np.errstate(all="ignore"):
+            lo, hi = sorted(signs * rng.random(2) * np.exp2(exponents))
+        if not (np.isfinite([lo, hi]).all() and lo < hi):
+            continue
+        n = int(rng.choice([2, 3, 7, 16, 256, 1024]))
+        reachable = reachable_float_bins(lo, hi, n)
+        assert centres_map_back(lo, hi, n)[reachable].all(), (lo, hi, n)
+        checked += 1
+
+
+def test_bin_position_never_decreases_with_the_value():
+    rng = np.random.default_rng(87)
+    for lo, hi, n in ((0.0, SMALLEST_SUBNORMAL, 4), (-MAX, MAX, 256), (1.0, 1.0 + 1e-15, 64)):
+        fractions = np.sort(rng.random(2000))
+        values = np.clip(lo * (1 - fractions) + hi * fractions, lo, hi)
+        assert (np.diff(float_bin_position(values, lo, hi, n)) >= 0).all()
+
+
+def test_constant_range_has_one_reachable_bin():
+    assert reachable_float_bins(2.5, 2.5, 4).tolist() == [True, False, False, False]
+    assert float_bin_indices(np.array([2.5]), 2.5, 2.5, 4).tolist() == [0]