From cb13561454ef7c18db2917cb53a3f62f05ccbecb Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 7 Nov 2026 18:47:39 +0000 Subject: [PATCH] Bin narrow float64 ranges with a bounded fraction Halving rounded a subnormal endpoint (half range zero for distinct endpoints) and n / half_range overflowed for very narrow ranges, mis-binning finite pixels. Use (value - lo) / (hi - lo) * n there; ordinary ranges keep their exact arithmetic. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Scma6W3LC11s3qjXHZuGDm --- CHANGELOG.md | 1 + src/morphobricks_normalization/bin_centres.py | 41 ++++++++++++++++++- src/morphobricks_normalization/histogram.py | 14 +++---- tests/test_histogram.py | 34 +++++++++++++++ tools/code_rules_baseline.json | 4 +- 5 files changed, 82 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c421c03..bcc3eea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ new refusal or a numerical-range fix (see the pull requests for the exact messag pooled mean (#72). - 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). - 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 diff --git a/src/morphobricks_normalization/bin_centres.py b/src/morphobricks_normalization/bin_centres.py index b9b642a..f39731f 100644 --- a/src/morphobricks_normalization/bin_centres.py +++ b/src/morphobricks_normalization/bin_centres.py @@ -1,9 +1,41 @@ -"""Centres of equal-width float bins, without overflow near the float64 limits.""" +"""Positions and centres of equal-width float bins, safe at both ends of the float64 range.""" from __future__ import annotations +import math + import numpy as np +_SMALLEST_NORMAL = float(np.finfo(np.float64).tiny) +_SAFE_HALF_RANGE = 2.0**-1000 # n / half_range stays far below the float64 maximum + + +def float_bin_position(values: np.ndarray | float, lo: float, hi: float, n: int) -> np.ndarray: + """Position of ``values`` in bin units, ``n`` bins over ``[lo, hi]`` (float64, not floored). + + Where the halved arithmetic ``(value / 2 - lo / 2) * (n / (hi / 2 - lo / 2))`` is exact + and finite it is kept, so existing histograms do not change. It is not for a narrow + range: halving rounds a subnormal endpoint (a half range of zero for distinct + endpoints) and ``n / half_range`` overflows to infinity. There the position is a + bounded fraction times ``n``, ``(value - lo) / (hi - lo) * n``, exact for subnormals; + if ``hi - lo`` itself overflows the halves are used for the fraction. ``lo`` must be + below ``hi``. + """ + values = np.asarray(values, dtype=np.float64) + half_range = hi * 0.5 - lo * 0.5 # halving keeps the range finite + if _halving_is_exact(lo, hi, half_range, n): + return (values * 0.5 - lo * 0.5) * (n / half_range) + span = hi - lo + if math.isfinite(span): + return (values - lo) / span * n + return (values * 0.5 - lo * 0.5) / half_range * n + + +def _halving_is_exact(lo: float, hi: float, half_range: float, n: int) -> bool: + """Whether halving loses nothing and ``n / half_range`` is finite (far from both limits).""" + exact = all(x == 0 or abs(x) >= 2 * _SMALLEST_NORMAL for x in (lo, hi)) + return exact and half_range >= n * _SAFE_HALF_RANGE + 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). @@ -13,8 +45,13 @@ def float_bin_centres(lo: float, hi: float, bin_count: int) -> np.ndarray: 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 + span = hi - lo + if math.isfinite(span): + # Half the fraction keeps the product below the range; no halving of the range, + # which would round to zero between distinct subnormal endpoints. + return lo + span * (fraction * 0.5) + half_range = hi * 0.5 - lo * 0.5 # halving keeps the range finite with np.errstate(over="ignore"): # an overflowing step is redone just below centres = lo + half_range * fraction over = ~np.isfinite(centres) diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py index 6510bec..da6c977 100644 --- a/src/morphobricks_normalization/histogram.py +++ b/src/morphobricks_normalization/histogram.py @@ -50,7 +50,7 @@ import numpy as np -from .bin_centres import float_bin_centres +from .bin_centres import float_bin_centres, float_bin_position 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 @@ -661,8 +661,6 @@ def _float_bins( if lo == hi and extrema is None: counts[0] = valid return counts - # Halving keeps hi - lo finite even for ranges near the float64 limits. - half_range = hi * 0.5 - lo * 0.5 for values in checked_blocks(blocks(), dtype): finite = values[np.isfinite(values)] if finite.size == 0: @@ -672,11 +670,11 @@ def _float_bins( finite = finite.astype(np.float64, copy=False) if extrema is not None: finite = _split_tails(finite, lo, hi, under, over) - if half_range == 0: + if lo == hi: # Robust limits that coincide: every pixel inside them has that value. counts[0] += finite.size continue - index = np.floor((finite * 0.5 - lo * 0.5) * (n / half_range)).astype(np.int64) + index = np.floor(float_bin_position(finite, lo, hi, n)).astype(np.int64) np.clip(index, 0, n - 1, out=index) counts += np.bincount(index, minlength=n) return counts @@ -858,7 +856,7 @@ def _bin_index( elif lo == hi: index = 0 else: - # The arithmetic of _float_bins, in Python floats (IEEE float64 as NumPy's). - half_range = hi * 0.5 - lo * 0.5 - index = math.floor((value * 0.5 - lo * 0.5) * (n / half_range)) + # The arithmetic of _float_bins; the value is clipped first, as its bin is. + clipped = min(max(float(value), float(lo)), float(hi)) + index = math.floor(float(float_bin_position(clipped, lo, hi, n))) return min(max(index, 0), n - 1) diff --git a/tests/test_histogram.py b/tests/test_histogram.py index 5c311f0..90be97e 100644 --- a/tests/test_histogram.py +++ b/tests/test_histogram.py @@ -137,6 +137,40 @@ def test_float_extreme_range_does_not_overflow(): assert np.all(np.isfinite(hist.representatives)) +SMALLEST_SUBNORMAL = 5e-324 + + +def test_float_range_of_the_smallest_subnormal_keeps_distinct_bins(): + plane = np.array([0.0, 0.0, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL, SMALLEST_SUBNORMAL]) + hist = compute_plane_histogram(plane, 4) + np.testing.assert_array_equal(hist.counts, [2, 0, 0, 3]) + assert hist.representatives[0] == 0.0 and hist.representatives[-1] == SMALLEST_SUBNORMAL + block = compute_block_histogram(lambda: iter((plane[:2], plane[2:])), plane.dtype, 4) + np.testing.assert_array_equal(block.counts, hist.counts) + assert histogram_module._bin_index(0.0, plane.dtype, 4, 0.0, SMALLEST_SUBNORMAL) == 0 + assert ( + histogram_module._bin_index(SMALLEST_SUBNORMAL, plane.dtype, 4, 0.0, SMALLEST_SUBNORMAL) + == 3 + ) + + +def test_float_narrow_range_does_not_overflow_the_bin_scale(): + # n / half_range is infinite for 1024 bins over 1e-306: the position must stay bounded. + plane = np.array([0.0, 0.5e-306, 1e-306]) + hist = compute_plane_histogram(plane, 1024) + assert np.flatnonzero(hist.counts).tolist() == [0, 512, 1023] + assert hist.counts.sum() == 3 + assert np.all(np.diff(hist.representatives) > 0) + assert histogram_module._bin_index(0.5e-306, plane.dtype, 1024, 0.0, 1e-306) == 512 + + +def test_float_range_between_subnormals_keeps_every_pixel_in_its_bin(): + plane = np.array([1e-320, 1.5e-320, 2e-320, 1e-320]) + hist = compute_plane_histogram(plane, 4) + assert hist.counts.sum() == 4 + np.testing.assert_array_equal(hist.counts, [2, 0, 1, 1]) + + def test_float32_input(): plane = rng.normal(0, 1, (32, 32)).astype(np.float32) hist = compute_plane_histogram(plane, 64) diff --git a/tools/code_rules_baseline.json b/tools/code_rules_baseline.json index 997d7e3..c02f3c2 100644 --- a/tools/code_rules_baseline.json +++ b/tools/code_rules_baseline.json @@ -28,11 +28,11 @@ "src/morphobricks_normalization/histogram.py _bin_pass: # parameters (max #)": 9, "src/morphobricks_normalization/histogram.py _float_bins: # Ruff findings too-many-arguments": 1, "src/morphobricks_normalization/histogram.py _float_bins: # parameters (max #)": 10, - "src/morphobricks_normalization/histogram.py _float_bins: private function of # lines without a docstring": 21, + "src/morphobricks_normalization/histogram.py _float_bins: private function of # lines without a docstring": 20, "src/morphobricks_normalization/histogram.py _integer_bins: # Ruff findings too-many-arguments": 1, "src/morphobricks_normalization/histogram.py _integer_bins: # parameters (max #)": 9, "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/histogram.py: # lines (max #)": 862, "src/morphobricks_normalization/pooling.py __init__: # Ruff findings missing-return-type-special-method": 1, "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,