From 5551e2e8a0af9bdcde787c3d6d9c97220d35a85f Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 11 Nov 2026 02:31:12 +0000 Subject: [PATCH 1/4] Report 64-bit integer representatives float64 cannot hold (#28) An int64/uint64 bin's representative was binning_min + (a + b) / 2, rounded twice, and its bound was (widest bin - 1) / 2 whatever float64 made of it, so [2**53 + 1] was reported exact with error 0 although its representative is 2**53. The representative is now (first + last) / 2 divided once as Python integers, the float64 nearest the midpoint. For a plane with a value of magnitude 2**51 or more the bound is the exact largest distance between each bin's (and tail's) integers and the float64 standing in for them, rounded up, so a histogram is exact only when every representative equals its value. Below 2**51 both are exact float64 arithmetic and unchanged. Bin counts and limits never change. Outputs that move: int64/uint64 lines of tools/dump_outputs.py (186 lines) and --corpus (90 lines), all checked against an exact Fraction oracle; the 400 golden fits are bit-identical. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G3qp9AxjpGPw9xfWAUJNqc --- AGENTS.md | 2 + CHANGELOG.md | 15 ++++ docs/contract.md | 8 +- pyproject.toml | 1 + src/morphobricks_normalization/histogram.py | 11 ++- .../integer_representatives.py | 65 ++++++++++++++ tests/test_histogram.py | 88 +++++++++++++++++++ tools/code_rules_baseline.json | 2 +- 8 files changed, 184 insertions(+), 8 deletions(-) create mode 100644 src/morphobricks_normalization/integer_representatives.py diff --git a/AGENTS.md b/AGENTS.md index 7001c12..602900b 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -34,6 +34,8 @@ src/morphobricks_normalization/ 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) + integer_representatives.py integer bin representatives (nearest float64 midpoint) and + their exact bound for 64-bit values float64 cannot hold (#28) 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/CHANGELOG.md b/CHANGELOG.md index 1317811..d59fd71 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,21 @@ Release notes and migration notes of `morphobricks-normalization`. Numbers are p issues of [CellMigrationLab/morphobricks-normalization](https://github.com/CellMigrationLab/morphobricks-normalization). The decisions are in [docs/plan.md](docs/plan.md) §12. +## Unreleased + +Numerical-correctness fixes. Each says which outputs change; every other output is +bit-identical (`tools/dump_outputs.py` and `--corpus` lines, and the 400 golden fits). + +- **64-bit integers float64 cannot hold (#28).** A bin's representative is the float64 nearest + its midpoint, rounded once (it was rounded twice, `binning_min` first, which could miss a + representable value by 2 above `2**53`). `max_representative_error` of an integer plane + with a value of magnitude `2**51` or more is the exact largest distance between its + integers and the float64 standing in for them, rounded up, so a histogram is exact only + when every representative equals its value: `[2**53 + 1]` is no longer reported exact with + error 0 (its representative is `2**53`; error 1), and `normalization_exact` and + `Histogram_Exact_C` are false for it. Only int64 and uint64 planes with such values + change; bin counts and limits never do. + ## v0.4.0 (2026-11-10) Breaking release: `PlaneRecord.position_id` is required (NM-04) and the readers of old table layouts diff --git a/docs/contract.md b/docs/contract.md index 90ab53d..03f670c 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -120,8 +120,14 @@ Notation in this section: `n` = `bin_count`, `lo` = `binning_min`, `hi` = `binni - If `W > n`: bin `i` holds the values `lo + ceil(i·W/n)` to `lo + ceil((i+1)·W/n) − 1`. Limits are whole numbers, consecutive bins touch without gaps or overlaps, and widths are `floor(W/n)` or `ceil(W/n)` values. -- Representative value of a bin: `(first value + last value) / 2`. +- Representative value of a bin: `(first value + last value) / 2`, the float64 nearest it. - Bound: every pixel is within `(widest bin − 1) / 2` of its representative. +- Planes with a value of magnitude `2**51` or more (64-bit types, #28): float64 cannot hold + every value there (`2**53 + 1` is stored as `2**53`), so the bound is the exact largest + distance between a bin's (or tail's) integers and the float64 standing in for them, rounded + up. A bin is exact only when its representative equals its value, so `normalization_exact` + is false for a group with such a value. Counts and bin limits stay exact integers. Below + `2**51` representatives and bound are exact float64 arithmetic, unchanged. ##### Float planes diff --git a/pyproject.toml b/pyproject.toml index eec9489..d5f89c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -103,6 +103,7 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.fitting" = "pure" "morphobricks_normalization.histogram" = "pure" "morphobricks_normalization.integrity" = "pure" +"morphobricks_normalization.integer_representatives" = "pure" "morphobricks_normalization.issues" = "pure" "morphobricks_normalization.normalization" = "pure" "morphobricks_normalization.parameter_checks" = "pure" diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py index da6c977..36a3f9e 100644 --- a/src/morphobricks_normalization/histogram.py +++ b/src/morphobricks_normalization/histogram.py @@ -52,6 +52,7 @@ from .bin_centres import float_bin_centres, float_bin_position from .dtypes import check_supported, default_bin_count, is_integer, validate_bin_count +from .integer_representatives import beyond_float_exact, integer_representatives, wide_integer_error from .integrity import check_histogram, check_tail, checked_blocks from .ranges import CENTRAL_PERCENTILES, HistogramSettings, RangeCheck from .sums import weighted_sum @@ -172,11 +173,7 @@ def representatives(self) -> np.ndarray: if self.is_empty: return np.full(self.bin_count, np.nan) if is_integer(self.dtype): - first, last = self._integer_bin_limits - return np.array( - [self.binning_min + (a + b) / 2 for a, b in zip(first, last, strict=True)], - dtype=np.float64, - ) + return integer_representatives(self) lo, hi = float(self.binning_min), float(self.binning_max) return float_bin_centres(lo, hi, self.bin_count) @@ -188,10 +185,12 @@ def max_representative_error(self) -> float: Float planes: half the bin width. Pixels of a tail are represented by the tail's mean, so they add the larger distance from that mean to the tail's ends: the plane's extremum and the binning limit. A histogram - with tails is therefore never exact. + with tails is therefore never exact. Beyond ``2**51``: see ``integer_representatives``. """ if self.is_empty: return 0.0 + if is_integer(self.dtype) and beyond_float_exact(self): + return wide_integer_error(self) if is_integer(self.dtype): first, last = self._integer_bin_limits widest = max(b - a + 1 for a, b in zip(first, last, strict=True)) diff --git a/src/morphobricks_normalization/integer_representatives.py b/src/morphobricks_normalization/integer_representatives.py new file mode 100644 index 0000000..0c14f09 --- /dev/null +++ b/src/morphobricks_normalization/integer_representatives.py @@ -0,0 +1,65 @@ +"""Representative values of integer bins, and their distance from the bins' pixels (#28). + +A bin of the integers ``first .. last`` stands for its pixels by its midpoint +``(first + last) / 2``, divided once as Python integers, so it is the float64 nearest the +midpoint. Below ``2**51`` every integer, every half-integer midpoint and every distance +between them is exact in float64, and the bound is ``(widest bin - 1) / 2``. From ``2**51`` +on (64-bit pixel types) float64 cannot hold every value: ``2**53 + 1`` is stored as +``2**53``. There the bound is the exact largest distance from a bin's or tail's integers to +the float64 standing in for them, rounded up, so a histogram is exact only when every +representative equals the pixel values it stands for. +""" + +from __future__ import annotations + +import math +from fractions import Fraction +from typing import TYPE_CHECKING + +import numpy as np + +if TYPE_CHECKING: # the histogram module imports this one + from .histogram import PlaneHistogram + +FLOAT_EXACT_LIMIT = 2**51 +"""Below this magnitude integer midpoints and their distances are exact float64 arithmetic.""" + + +def integer_representatives(histogram: PlaneHistogram) -> np.ndarray: + """Return the float64 nearest the midpoint of every bin of a non-empty integer plane.""" + return np.array( + [(first + last) / 2 for first, last in histogram.integer_bin_limits()], dtype=np.float64 + ) + + +def beyond_float_exact(histogram: PlaneHistogram) -> bool: + """Whether a non-empty integer plane has a pixel of magnitude ``2**51`` or more.""" + return max(-int(histogram.pixel_min), int(histogram.pixel_max)) >= FLOAT_EXACT_LIMIT + + +def wide_integer_error(histogram: PlaneHistogram) -> float: + """Return the exact bound of an integer plane beyond ``2**51``, rounded up to a float64. + + Every bin that holds an integer, and every non-empty tail (its mean against the plane's + extremum and the binning limit), counts with the larger exact distance from the value + standing in for it to its two ends. + """ + ends = [ + (rep, first, last) + for rep, (first, last) in zip( + histogram.representatives.tolist(), histogram.integer_bin_limits(), strict=True + ) + if first <= last # a bin that holds no integer has no pixel to stand in for + ] + for tail, first, last in ( + (histogram.underflow, histogram.pixel_min, histogram.binning_min), + (histogram.overflow, histogram.binning_max, histogram.pixel_max), + ): + if not tail.is_empty: + ends.append((tail.mean, first, last)) + worst = max( + max(abs(Fraction(rep) - int(first)), abs(Fraction(rep) - int(last))) + for rep, first, last in ends + ) + nearest = float(worst) + return nearest if Fraction(nearest) >= worst else math.nextafter(nearest, math.inf) diff --git a/tests/test_histogram.py b/tests/test_histogram.py index 90be97e..27ef002 100644 --- a/tests/test_histogram.py +++ b/tests/test_histogram.py @@ -1,12 +1,19 @@ import bisect +from fractions import Fraction import numpy as np import pytest from morphobricks_normalization import ( + NormalizationSettings, + PlaneRecord, + PooledHistogram, UnsupportedDtypeError, compute_block_histogram, compute_plane_histogram, + fit, + read_plane_histograms, + write_plane_histograms, ) from morphobricks_normalization import histogram as histogram_module @@ -344,3 +351,84 @@ def blocks(): floats = plane.astype(np.float32) compute_block_histogram(lambda: (calls.append(1), iter([floats]))[1], floats.dtype, 4) assert len(calls) == 2 # floats are always read twice + + +# 64-bit integers that float64 cannot hold (#28) --------------------------------------- + +UNREPRESENTABLE = [ + ("uint64", 2**53 - 1), + ("uint64", 2**53), + ("uint64", 2**53 + 1), + ("uint64", 2**64 - 1), + ("int64", 2**63 - 1), + ("int64", -(2**63)), + ("int64", -(2**63) + 1), + ("int64", -(2**53) - 1), +] + + +@pytest.mark.parametrize(("dtype", "value"), UNREPRESENTABLE) +def test_singleton_64_bit_value_is_exact_only_if_float64_holds_it(dtype, value): + hist = compute_plane_histogram(np.array([[value]], dtype=dtype), 256) + (occupied,) = np.flatnonzero(hist.counts) + representative = hist.representatives[occupied] + error = abs(Fraction(float(representative)) - value) # exact, independent of the package + assert hist.counts[occupied] == 1 and hist.pixel_min == hist.pixel_max == value + assert hist.is_exact == (error == 0) + assert hist.max_representative_error == float(error) # 0 or 1: representable + + +def test_adjacent_64_bit_integers_sharing_a_float_are_not_exact(): + plane = np.array([[2**53, 2**53 + 1, 2**53 + 1]], dtype=np.uint64) + hist = compute_plane_histogram(plane, 256) + np.testing.assert_array_equal(np.flatnonzero(hist.counts), [0, 1]) # bins stay exact + assert hist.representatives[0] == hist.representatives[1] == 2.0**53 + assert not hist.is_exact and hist.max_representative_error == 1.0 + pooled = PooledHistogram([hist]) + assert pooled.counts.tolist() == [3] and pooled.max_representative_error == 1.0 + record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a") + settings = NormalizationSettings(method="minmax", scope="image") + with pytest.raises(ValueError, match="zero"): # 2**53 + 1 - 2**53 rounds to 0 in float64 + fit([record], settings, group="a.tif") + plane = np.array([[2**53 - 2, 2**53 + 1]], dtype=np.uint64) + parameters = fit( + [ + PlaneRecord( + "a.tif", "c", "r", "t", [compute_plane_histogram(plane, 256)], position_id="a" + ) + ], + settings, + group="a.tif", + ) + assert parameters.exact is False + + +@pytest.mark.parametrize("dtype", ["int64", "uint64"]) +@pytest.mark.parametrize("bins", [2, 7, 256]) +def test_wide_64_bit_bound_covers_every_pixel_exactly(dtype, bins): + info = np.iinfo(dtype) + plane = np.random.default_rng(bins).integers( + info.min, info.max, (24, 31), dtype=dtype, endpoint=True + ) + plane[0, :2] = (info.min, info.max) + hist = compute_plane_histogram(plane, bins) + block = compute_block_histogram(lambda: iter((plane[:5], plane[5:])), plane.dtype, bins) + assert block.max_representative_error == hist.max_representative_error + starts = [first for first, _ in hist.integer_bin_limits()] + representatives = [Fraction(float(r)) for r in hist.representatives] + worst = max( + abs(representatives[bisect.bisect_right(starts, int(v)) - 1] - int(v)) + for v in plane.ravel().tolist() + ) + assert Fraction(hist.max_representative_error) >= worst and not hist.is_exact + + +def test_64_bit_bound_survives_the_statistics_file(tmp_path): + plane = np.array([[2**53 + 1, 2**53 + 3, 2**53 + 4]], dtype=np.uint64) + hist = compute_plane_histogram(plane, 256) + record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a") + (back,) = read_plane_histograms(write_plane_histograms([record], tmp_path / "p.csv")) + read = back.histograms[0] + assert read.pixel_min == 2**53 + 1 and read.pixel_max == 2**53 + 4 + assert read.max_representative_error == hist.max_representative_error == 1.0 + assert not read.is_exact diff --git a/tools/code_rules_baseline.json b/tools/code_rules_baseline.json index 4848f84..ac9f037 100644 --- a/tools/code_rules_baseline.json +++ b/tools/code_rules_baseline.json @@ -32,7 +32,7 @@ "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 #)": 862, + "src/morphobricks_normalization/histogram.py: # lines (max #)": 861, "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, From 47e5339c5e4e45e0320e1f38f9d991cd2e2db97d Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 11 Nov 2026 02:34:58 +0000 Subject: [PATCH 2/4] Keep pooled standard deviations of large finite values finite (#33) PooledHistogram.std and WeightedHistogram.std summed squared deviations in float64, which overflow beyond about 1e154: [-1e200, 1e200] gave inf (so zscore was refused) although its standard deviation is 1e200, and a tail could raise a bare OverflowError from float ** 2. Both now go through _std_with_tails. The plain sum is tried first, so every finite result is bit-identical; when it is not finite but the mean and every value are, it is redone with all values divided by a power of two near the largest (exact), and multiplied back. A tail whose stored sum of squares is infinite (its pixels' squares overflowed when the file was written) is refused by name, since its spread is not known; mean, percentiles, median and MAD remain available. Tests compare with an exact Fraction/Decimal oracle: +-1e200, 1e200 +- 1e185, +-1e308, constant 1e200, weighted parts and robust tails near 1e153, plus a CSV round trip of an infinite tail. tools/dump_outputs.py and --corpus are unchanged by this commit; the 400 golden fits are bit-identical. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G3qp9AxjpGPw9xfWAUJNqc --- CHANGELOG.md | 10 +++ docs/contract.md | 6 +- src/morphobricks_normalization/pooling.py | 65 +++++++++++++-- tests/test_pooling.py | 97 ++++++++++++++++++++++- 4 files changed, 168 insertions(+), 10 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d59fd71..35c0375 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,16 @@ bit-identical (`tools/dump_outputs.py` and `--corpus` lines, and the 400 golden error 0 (its representative is `2**53`; error 1), and `normalization_exact` and `Histogram_Exact_C` are false for it. Only int64 and uint64 planes with such values change; bin counts and limits never do. +- **Standard deviations of large finite values (#33).** `PooledHistogram.std` and + `WeightedHistogram.std` summed squared deviations that overflow for values beyond about + `1e154` (`[-1e200, 1e200]` gave `inf`, so `zscore` was refused, although its standard + deviation is `1e200`; a tail could raise a bare `OverflowError`). When the plain sum is not + finite but the mean and every value are, it is redone with the values divided by a power of + two, which is exact; every finite result is unchanged bit for bit. A tail whose stored sum + of squares is infinite (schema version 4 keeps raw sums of squares) is refused by name + instead of giving `inf`: its spread is not known. No dump output changes. Not changed: the + raw sums of squares of the file lose precision by cancellation for tails of large values + with a small spread; a centred representation needs a new schema version. ## v0.4.0 (2026-11-10) diff --git a/docs/contract.md b/docs/contract.md index 03f670c..debcd36 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -198,7 +198,7 @@ tail's mean and its count. | Statistic | Definition | Accuracy | |---|---|---| | min, max | Smallest `pixel_min`, largest `pixel_max` | Exact | -| mean, std | Count-weighted mean and population std (÷N) of representatives; tails count with their exact sums and sums of squares. The weighted sums are correctly rounded sums of the products (`math.fsum`), so they never depend on the CPU thread count or BLAS library (#67) | Exact for exact histograms; otherwise very close (errors cancel) | +| mean, std | Count-weighted mean and population std (÷N) of representatives; tails count with their exact sums and sums of squares. The weighted sums are correctly rounded sums of the products (`math.fsum`), so they never depend on the CPU thread count or BLAS library (#67). When the squares overflow float64 although every value is finite (`[-1e200, 1e200]`), the std is computed with the values divided by a power of two, exactly, so a finite spread is never reported as infinite (#33) | Exact for exact histograms; otherwise very close (errors cancel) | | percentile q | Smallest representative with cumulative count ≥ ceil(q/100 · N) (numpy `inverted_cdf`) | Within the group's largest per-plane bound | | median | 50th percentile | As percentiles | | MAD | Median of \|representative − median\| | Within twice the bound | @@ -260,6 +260,10 @@ pixels keeps its extrema as limits; a 512 × 512 plane has 26 pixels in each tai - Percentiles, median and MAD: each non-empty tail stands as one value, its mean, with its count. - Mean and standard deviation use the tails' sums: as exact as without tails. +- A tail whose sum of squares is not finite (its pixels' squares overflowed float64 when the + file was written, values beyond about `1.3e154`) has no known spread: the standard + deviation, and so `zscore`, is refused with a message naming the tail (#33). Mean, + percentiles, median and MAD do not use it. - Min-max parameters use the exact extrema (`pixel_min`, `pixel_max`). - `max_representative_error` is the larger of the bins' bound and, for each non-empty tail, the larger distance from the tail's mean to its two ends (`max(mean − pixel_min, diff --git a/src/morphobricks_normalization/pooling.py b/src/morphobricks_normalization/pooling.py index 10abd90..eebd8a2 100644 --- a/src/morphobricks_normalization/pooling.py +++ b/src/morphobricks_normalization/pooling.py @@ -167,9 +167,7 @@ def std(self) -> float: return _std_with_tails( self.bin_values, self.bin_counts, self.tails, self.pixel_count, self.mean ) - deviations = self.values - self.mean - squares = weighted_sum(self.counts, deviations * deviations) - return float(np.sqrt(squares / self.pixel_count)) + return _std_with_tails(self.values, self.counts, [], self.pixel_count, self.mean) def percentile(self, q: float) -> float: """Inverse-CDF percentile, ``q`` in percent (0-100).""" @@ -307,9 +305,7 @@ def std(self) -> float: return _std_with_tails( self.bin_values, self.bin_weights, self.tails, self.total_weight, self.mean ) - deviations = self.values - self.mean - squares = weighted_sum(self.weights, deviations * deviations) - return float(np.sqrt(squares / self.total_weight)) + return _std_with_tails(self.values, self.weights, [], self.total_weight, self.mean) def percentile(self, q: float) -> float: """Weighted inverse-CDF percentile, ``q`` in percent (0-100).""" @@ -393,14 +389,67 @@ def _std_with_tails( total: float, mean: float, ) -> float: - deviations = values - mean + """Return the population standard deviation of the bins (and tails) around ``mean`` (#33). + + The squares are first summed as they are, which keeps every finite result bit for + bit. When that overflows although the mean and every value are finite (``[-1e200, + 1e200]`` has standard deviation ``1e200`` but squares of ``1e400``), the sum is redone + with every value divided by a power of two near the largest, which is exact, so the finite + spread is found. A tail whose stored sum of squares is not finite is refused by name. + """ + squares = [s2 for _, _, s2 in tails] + if not all(math.isfinite(s2) for s2 in squares): + raise ValueError(_non_finite_squares_message(squares)) + try: + with np.errstate(over="ignore", invalid="ignore"): # an overflow is redone below + std = _std_scaled(values, weights, tails, total, mean, 1.0) + except OverflowError: # float ** 2 raises where NumPy gives inf + std = math.inf + if math.isfinite(std) or not (math.isfinite(mean) and np.isfinite(values).all()): + return std + largest = max(float(np.max(np.abs(values), initial=0.0)), abs(mean)) + largest = max([largest] + [abs(s / count) for count, s, _ in tails]) + # Values divided by 2**(exponent - 1) are below 2 in magnitude, and the scale is finite. + return _std_scaled( + values, weights, tails, total, mean, math.ldexp(1.0, math.frexp(largest)[1] - 1) + ) + + +def _std_scaled( + values: np.ndarray, + weights: np.ndarray, + tails: list[tuple[float, float, float]], + total: float, + mean: float, + scale: float, +) -> float: + """Return the standard deviation computed in units of ``scale``, a power of two (1: as is).""" + mean = mean / scale + deviations = values / scale - mean squares = weighted_sum(weights, deviations * deviations) for count, s, s2 in tails: # Each tail adds its spread around its own mean and its mean's distance # from the overall mean; the spread never goes negative through rounding. + s, s2 = s / scale, s2 / scale / scale tail_mean = s / count squares += max(s2 - s * tail_mean, 0.0) + count * (tail_mean - mean) ** 2 - return float(np.sqrt(squares / total)) + return float(np.sqrt(squares / total)) * scale + + +def _non_finite_squares_message(squares: list[float]) -> str: + """Refusal text for tails whose sum of squares overflowed float64 when it was computed.""" + bad = ", ".join( + f"tail {number} of {len(squares)} has sum of squares {s2!r}" + for number, s2 in enumerate(squares, 1) + if not math.isfinite(s2) + ) + return ( + f"The standard deviation cannot be pooled: {bad}. The squares of pixels outside the " + "robust range overflowed float64 when the statistics file was written (pixel values " + "beyond about 1.3e154), so their spread is not known, although it may be finite. Use " + "a method that does not need the standard deviation (percentile or robust z-score), " + "or the min-max range." + ) def _channel_histogram(record: PlaneRecord, channel: int | None) -> PlaneHistogram: diff --git a/tests/test_pooling.py b/tests/test_pooling.py index 154481d..4333506 100644 --- a/tests/test_pooling.py +++ b/tests/test_pooling.py @@ -2,13 +2,23 @@ import pytest from morphobricks_normalization import ( + NormalizationSettings, NoValidPixelsError, + PlaneRecord, PooledHistogram, compute_plane_histogram, + fit, group_records, pooled_statistics, + read_plane_histograms, + write_plane_histograms, + z_score, +) +from morphobricks_normalization.pooling import ( + DEFAULT_PERCENTILES, + WeightedHistogram, + _linear_percentile, ) -from morphobricks_normalization.pooling import DEFAULT_PERCENTILES, _linear_percentile from morphobricks_normalization.ranges import HistogramSettings from simulation import direct_stats, make_records, simulate_experiment @@ -142,3 +152,88 @@ def test_percentiles_from_pixels_near_float64_limits_are_finite_and_exact(q, lim assert np.isfinite(result) assert abs(Fraction(result) - expected) <= abs(expected) * Fraction(1, 10**14) assert np.isfinite(pooled.percentile(q)) + + +# Standard deviations of large finite values (#33) ------------------------------------- + + +def exact_std(values, weights, tails=()): + """Population standard deviation of weighted values and (count, sum, squares) tails. + + Every float is taken at its exact value (Fraction) and the square root in 40 digits, so + nothing here can overflow or round like the package's float64 arithmetic. + """ + from decimal import Decimal, localcontext + from fractions import Fraction + + values = [Fraction(v) for v in np.asarray(values, dtype=float).tolist()] + weights = [Fraction(w) for w in np.asarray(weights, dtype=float).tolist()] + tails = [tuple(Fraction(x) for x in tail) for tail in tails] + total = sum(weights) + sum(c for c, _, _ in tails) + mean = ( + sum(w * v for w, v in zip(weights, values, strict=True)) + sum(s for _, s, _ in tails) + ) / total + squares = sum(w * (v - mean) ** 2 for w, v in zip(weights, values, strict=True)) + squares += sum(s2 - s * s / c + c * (s / c - mean) ** 2 for c, s, s2 in tails) + variance = squares / total + with localcontext() as context: + context.prec = 40 + return float((Decimal(variance.numerator) / Decimal(variance.denominator)).sqrt()) + + +LARGE_PLANES = { + "opposite 1e200": [-1e200, 1e200], + "1e200 +- 1e185": [1e200 - 1e185, 1e200, 1e200 + 1e185], + "opposite 1e308": [-1e308, 0.0, 1e308], + "constant 1e200": [1e200] * 5, +} + + +@pytest.mark.parametrize("name", LARGE_PLANES) +def test_std_of_large_finite_values_is_finite_and_exact(name): + hist = compute_plane_histogram(np.array(LARGE_PLANES[name]), 64) + pooled = PooledHistogram([hist]) + expected = exact_std(pooled.values, pooled.counts) + assert np.isfinite(expected) + assert pooled.std == pytest.approx(expected, rel=1e-14, abs=0) + if name != "constant 1e200": + parameters = z_score(pooled) + assert parameters.scale == pooled.std and np.isfinite(parameters.offset) + + +def test_weighted_std_of_large_finite_values_is_finite_and_exact(): + parts = [ + PooledHistogram([compute_plane_histogram(np.array(LARGE_PLANES[name]), 64)]) + for name in ("opposite 1e200", "1e200 +- 1e185") + ] + weighted = WeightedHistogram(parts, [0.5, 2.0]) + expected = exact_std(weighted.values, weighted.weights) + assert weighted.std == pytest.approx(expected, rel=1e-14, abs=0) + + +def test_std_with_large_tails_is_found_by_scaling(): + plane = np.random.default_rng(7).normal(0, 1, 10000) * 1e153 + settings = HistogramSettings(histogram_range="robust", robust_quantiles=(0.05, 99.95)) + hist = compute_plane_histogram(plane, 256, settings=settings) + assert hist.has_tails and np.isfinite(hist.overflow.sum_of_squares) + pooled = PooledHistogram([hist]) + expected = exact_std(pooled.bin_values, pooled.bin_counts, pooled.tails) + assert pooled.std == pytest.approx(expected, rel=1e-12, abs=0) + assert pooled.std == pytest.approx(np.std(plane / 1e153) * 1e153, rel=1e-4) + + +def test_tail_whose_squares_overflowed_is_a_named_refusal(tmp_path): + plane = np.random.default_rng(1).normal(0, 1, 10000) * 1e150 + plane[:5] = (1e200, 2e200, -1e200, 3e200, -2e200) + with np.errstate(over="ignore"): # the tail's squares overflow as the file is written + hist = compute_plane_histogram(plane, 256, settings=ROBUST) + assert hist.overflow.sum_of_squares == np.inf + record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a") + (back,) = read_plane_histograms(write_plane_histograms([record], tmp_path / "p.csv")) + pooled = PooledHistogram([back.histograms[0]]) + assert np.isfinite(pooled.mean) and np.isfinite(pooled.percentile(50)) + message = "pooled: tail 1 of 2 has sum of squares inf, tail 2 of 2 has sum of squares inf" + with pytest.raises(ValueError, match=message): + pooled.summary() # the standard deviation; mean and percentiles stay available + with pytest.raises(ValueError, match="percentile or robust z-score"): + fit([back], NormalizationSettings(method="zscore", scope="image"), group="a.tif") From c084a35979a82adee630fd49f28a5686f988a981 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 11 Nov 2026 02:37:07 +0000 Subject: [PATCH 3/4] Check missing conditions in each per-time-point repeat group (#41) With scope="biological_repeat" and per_time_point=True the fitted groups are (label, t), but repeat_missing_conditions gathered the conditions of each label over all time points. Two rows (control, r1, t=0) and (treated, r1, t=1) gave the groups r1::t=0 and r1::t=1, one condition each, and no warning. The check now runs on the groups scope_groups returns, so it uses the same keys as the fit, and names them with group_name (r1::t=0). The expected set stays every condition of the statistics (contract 4.2). Without per_time_point the keys, items and message are unchanged. Controls: complete same-time groups, a condition missing an entire time point in every repeat, renamed pairs, pooled time and the image scope. tools/dump_outputs.py and --corpus are unchanged by this commit; the 400 golden fits are bit-identical. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01G3qp9AxjpGPw9xfWAUJNqc --- CHANGELOG.md | 6 ++++ docs/contract.md | 4 ++- src/morphobricks_normalization/scopes.py | 25 +++++++-------- tests/test_scopes.py | 39 ++++++++++++++++++++++++ tools/code_rules_baseline.json | 1 - 5 files changed, 61 insertions(+), 14 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 35c0375..eca6a70 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,12 @@ bit-identical (`tools/dump_outputs.py` and `--corpus` lines, and the 400 golden instead of giving `inf`: its spread is not known. No dump output changes. Not changed: the raw sums of squares of the file lose precision by cancellation for tails of large values with a small spread; a centred representation needs a new schema version. +- **Missing conditions per time point (#41).** With `scope="biological_repeat"` and + `per_time_point=True`, `repeat_missing_conditions` checked each repeat group over all its + time points, not the groups that are fitted, so a group such as `r1::t=0` with control + images only gave no warning when treated images of `r1` existed at another time point. + The groups that are fitted are now checked, and the warning's items and message name them + (`r1::t=0`). Without `per_time_point` the warning is unchanged. No dump output changes. ## v0.4.0 (2026-11-10) diff --git a/docs/contract.md b/docs/contract.md index debcd36..47b8abd 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -466,7 +466,9 @@ y_c = (x_c − offset) / scale for c = R, G, B forms a group with the other conditions' `single` rows, or with the repeat it is paired with. - **Incomplete groups:** a repeat group without rows from every condition of the statistics - gives the warning `repeat_missing_conditions` (§9.2) and is still used. + gives the warning `repeat_missing_conditions` (§9.2) and is still used. With + `per_time_point=true` each time point's group is checked on its own and named in the + warning as `