Fixes #28, fixes #33, fixes #41

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

CHANGELOG.md

@@ -4,6 +4,59 @@ 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<i>` are false for it. Only int64 and uint64 planes with such values
+  change; bin counts and limits never do.
+
+  **Who is affected, and by how much.** Only int64 and uint64 images with a pixel of
+  magnitude `2**51` (about `2.25e15`) or more. Every other pixel type (8-, 16- and 32-bit
+  integers, float32, float64), and 64-bit images whose values all stay below `2**51`, give
+  bit-identical numbers. For the affected images:
+  - representatives move by at most one float64 step (at most 1536 for int64 and 2048 for
+    uint64 values near `9e18` and `1.8e19`, about `1e-16` relative);
+  - `max_representative_error` is the true largest distance: for bins of one value float64
+    cannot hold it grows from 0 to at most half a float64 step (for example 1 for
+    `[2**53 + 1]`, 170 and 341 for the constant planes `3074457345618258602` and
+    `6148914691236517205`), so `is_exact`, `normalization_exact` and
+    `Histogram_Exact_C<i>` change from true to false; for wider bins it changes in its last
+    bits only;
+  - pooled means, standard deviations and percentiles, and fitted offsets and scales, move by
+    at most about `3e-14` relative in the `tools/dump_outputs.py` cases (a zscore offset
+    `9618469625711740` becomes `9618469625711942`, nearer the exact model mean
+    `9618469625711938`).
+
+  Statistics files need not be recomputed: they store counts and limits, which do not change,
+  and representatives and bounds are derived when the file is read, so an existing
+  `plane_histograms.csv` gives the new numbers. Parameters fitted with v0.4.0 from such
+  images differ from a refit by the amounts above.
+- **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.
+- **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)
 
 Breaking release: `PlaneRecord.position_id` is required (NM-04) and the readers of old table layouts

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
 
@@ -192,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 |
@@ -254,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,
@@ -456,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 `<label>::t=<t>` (#41).
 
 ### 4.3 Pooling weights
 

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"

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

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)

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: