From 11bf51a69df57db00f13b763ffcc5f944230ecfa Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 5 Nov 2026 07:48:43 +0000 Subject: [PATCH 1/2] Make pooled means and standard deviations independent of CPU threads (#67) The pooled mean and standard deviation (and so every zscore fit) summed with np.dot, which OpenBLAS splits across its threads above 10,000 elements, so the last bits of a fit changed with OMP_NUM_THREADS/OPENBLAS_NUM_THREADS. sums.weighted_sum multiplies element by element (each product rounded once, as np.dot did) and adds the products with math.fsum: the correctly rounded sum, identical for any thread count, BLAS or platform. All seven np.dot calls of the library go through it (PooledHistogram and WeightedHistogram mean/std, the tail-aware mean/std, and the per-value tail sums). The other reductions (np.sum, cumsum, bincount, add.at) are single-threaded NumPy loops. Scientific precision: tools/dump_outputs.py changes on 86 of 1,788 lines (134 numbers, pooled mean/std and zscore offset/scale only), by at most 3.7e-14 relative; 89 numbers by 1 ulp, the largest 317 ulps (int64 planes near 2**53). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WX6rdMV1r5jUPUR94kxV5a --- README.md | 8 ++ docs/contract.md | 2 +- pyproject.toml | 1 + src/morphobricks_normalization/histogram.py | 4 +- src/morphobricks_normalization/pooling.py | 15 +-- src/morphobricks_normalization/sums.py | 35 ++++++ tests/test_sums.py | 113 ++++++++++++++++++++ 7 files changed, 169 insertions(+), 9 deletions(-) create mode 100644 src/morphobricks_normalization/sums.py create mode 100644 tests/test_sums.py diff --git a/README.md b/README.md index edc56f5..3708d35 100644 --- a/README.md +++ b/README.md @@ -172,6 +172,14 @@ warnings = range_check_issues([("img.tif::z=0,t=0", 0, check)]) min-max histogram. A `plane_histograms.csv` row with one blank intensity limit and valid pixels is refused with a `CsvFormatError` naming the column. Valid inputs give bit-identical outputs (`python tools/dump_outputs.py` prints the digest to compare). +- **Means and standard deviations no longer depend on CPU threads (#67).** They summed with + `np.dot`, which the BLAS library splits across its threads, so a z-score fit changed in its + last bits with `OMP_NUM_THREADS` or `OPENBLAS_NUM_THREADS`. They now add the same products + with `math.fsum` (the correctly rounded sum): the same bits for any thread count, BLAS or + platform, never less accurate than before. **The last bits change:** pooled `mean` and + `std` and the `zscore` offset and scale differ from v0.3.3 by at most 3.7e-14 relative + (134 numbers on 86 of the 1,788 lines of `tools/dump_outputs.py`; 89 by 1 ulp, the largest + 317 ulps, for int64 planes near 2**53). Every other output is bitwise unchanged. ### Migrating to v0.3.3 diff --git a/docs/contract.md b/docs/contract.md index 9232cc2..d9bda27 100644 --- a/docs/contract.md +++ b/docs/contract.md @@ -185,7 +185,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 | 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) | 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 | diff --git a/pyproject.toml b/pyproject.toml index 9d9f916..afc4ff8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -113,3 +113,4 @@ baseline = "tools/code_rules_baseline.json" "morphobricks_normalization.scopes" = "pure" "morphobricks_normalization.settings" = "pure" "morphobricks_normalization.summaries" = "pure" +"morphobricks_normalization.sums" = "pure" diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py index 1a67eb6..830c2fd 100644 --- a/src/morphobricks_normalization/histogram.py +++ b/src/morphobricks_normalization/histogram.py @@ -53,6 +53,7 @@ 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 +from .sums import weighted_sum # Above this many distinct values, integer planes are binned with a search over # bin starts instead of a full per-value count. @@ -721,9 +722,8 @@ def _table_tail(per_value: np.ndarray, first_value: int) -> HistogramTail: if count == 0: return EMPTY_TAIL values = first_value + np.arange(per_value.size, dtype=np.float64) - weights = per_value.astype(np.float64) return HistogramTail( - count, float(np.dot(values, weights)), float(np.dot(values * values, weights)) + count, weighted_sum(per_value, values), weighted_sum(per_value, values * values) ) diff --git a/src/morphobricks_normalization/pooling.py b/src/morphobricks_normalization/pooling.py index abff086..5ee29a1 100644 --- a/src/morphobricks_normalization/pooling.py +++ b/src/morphobricks_normalization/pooling.py @@ -41,6 +41,7 @@ from .histogram import PlaneHistogram from .integrity import INT64_MAX from .records import PlaneRecord +from .sums import weighted_sum DEFAULT_PERCENTILES = (0.1, 1.0, 5.0, 25.0, 50.0, 75.0, 95.0, 99.0, 99.9) @@ -157,7 +158,7 @@ def mean(self) -> float: """Count-weighted mean; pixels outside the bins count with their exact values.""" if self.tails: return _mean_with_tails(self.bin_values, self.bin_counts, self.tails, self.pixel_count) - return float(np.dot(self.counts, self.values) / self.pixel_count) + return weighted_sum(self.counts, self.values) / self.pixel_count @property def std(self) -> float: @@ -167,7 +168,8 @@ def std(self) -> float: self.bin_values, self.bin_counts, self.tails, self.pixel_count, self.mean ) deviations = self.values - self.mean - return float(np.sqrt(np.dot(self.counts, deviations * deviations) / self.pixel_count)) + squares = weighted_sum(self.counts, deviations * deviations) + return float(np.sqrt(squares / self.pixel_count)) def percentile(self, q: float) -> float: """Inverse-CDF percentile, ``q`` in percent (0-100).""" @@ -296,7 +298,7 @@ def mean(self) -> float: return _mean_with_tails( self.bin_values, self.bin_weights, self.tails, self.total_weight ) - return float(np.dot(self.weights, self.values) / self.total_weight) + return weighted_sum(self.weights, self.values) / self.total_weight @property def std(self) -> float: @@ -306,7 +308,8 @@ def std(self) -> float: self.bin_values, self.bin_weights, self.tails, self.total_weight, self.mean ) deviations = self.values - self.mean - return float(np.sqrt(np.dot(self.weights, deviations * deviations) / self.total_weight)) + squares = weighted_sum(self.weights, deviations * deviations) + return float(np.sqrt(squares / self.total_weight)) def percentile(self, q: float) -> float: """Weighted inverse-CDF percentile, ``q`` in percent (0-100).""" @@ -358,7 +361,7 @@ def pooled_statistics( def _mean_with_tails( values: np.ndarray, weights: np.ndarray, tails: list[tuple[float, float, float]], total: float ) -> float: - return float((np.dot(weights, values) + math.fsum(s for _, s, _ in tails)) / total) + return (weighted_sum(weights, values) + math.fsum(s for _, s, _ in tails)) / total def _std_with_tails( @@ -369,7 +372,7 @@ def _std_with_tails( mean: float, ) -> float: deviations = values - mean - squares = float(np.dot(weights, deviations * deviations)) + 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. diff --git a/src/morphobricks_normalization/sums.py b/src/morphobricks_normalization/sums.py new file mode 100644 index 0000000..32c0140 --- /dev/null +++ b/src/morphobricks_normalization/sums.py @@ -0,0 +1,35 @@ +"""Sums whose result never depends on the CPU thread count (issue #67). + +``np.dot`` hands a long float64 dot product to the BLAS library, which splits it across its +threads and adds the partial sums in an order that depends on how many threads it has +(OpenBLAS threads a dot product of more than 10,000 elements). A pooled mean or standard +deviation, and so a z-score fit, then changed in its last bits with ``OMP_NUM_THREADS`` or +``OPENBLAS_NUM_THREADS``. :func:`weighted_sum` multiplies element by element (each product +rounded once, as in ``np.dot``) and adds the products with :func:`math.fsum`, which returns +their correctly rounded sum: the same bits for any thread count, BLAS library, NumPy version or +platform, and never less accurate than the dot product it replaces. +""" + +from __future__ import annotations + +import math + +import numpy as np + + +def weighted_sum(weights: np.ndarray, values: np.ndarray) -> float: + """Return ``sum(weights * values)``, correctly rounded and independent of threads. + + Parameters + ---------- + weights, values + One-dimensional arrays of the same length (counts or weights, and the values they + weigh); converted to float64 before multiplying, as ``np.dot`` does. + + Returns + ------- + total : float + The correctly rounded sum of the float64 products. + """ + products = np.multiply(weights, values, dtype=np.float64) + return math.fsum(products.tolist()) diff --git a/tests/test_sums.py b/tests/test_sums.py new file mode 100644 index 0000000..b3f910e --- /dev/null +++ b/tests/test_sums.py @@ -0,0 +1,113 @@ +"""Pooled means, standard deviations and z-score fits never depend on CPU threads (#67). + +``np.dot`` let the BLAS library split a long sum across its threads, so the last bits of a +z-score fit changed with ``OPENBLAS_NUM_THREADS``/``OMP_NUM_THREADS`` (OpenBLAS threads a dot +product of more than 10,000 elements). The sums are now correctly rounded sums of the products. +""" + +from __future__ import annotations + +import json +import os +import subprocess +import sys +from fractions import Fraction +from pathlib import Path + +import numpy as np + +from morphobricks_normalization import ( + HistogramSettings, + NormalizationSettings, + PlaneRecord, + PooledHistogram, + compute_plane_histogram, + fit, +) +from morphobricks_normalization.pooling import WeightedHistogram +from morphobricks_normalization.sums import weighted_sum + +THREAD_VARIABLES = ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS") +TIGHT = HistogramSettings(histogram_range="robust", robust_quantiles=(2.0, 97.0)) + + +def _records(settings: HistogramSettings | None = None) -> list[PlaneRecord]: + """Planes whose pooled list holds far more than 10,000 distinct values.""" + rng = np.random.default_rng(67) + records = [] + for i in range(60): + plane = rng.gamma(2.0, 300.0 + 7 * i, size=(64, 64)) + hist = compute_plane_histogram(plane, 256, settings=settings) + records.append(PlaneRecord(f"c/r/img{i}.tif", "c", f"r{i % 2}", "t", [hist])) + return records + + +def _numbers() -> dict[str, str]: + """Every sum that once went through ``np.dot``, as bit-exact text.""" + plain = PooledHistogram(_records()) + tails = PooledHistogram(_records(TIGHT)) + assert plain.values.size > 10_000 and tails.tails + weighted = WeightedHistogram([plain, tails], [1.0, 0.37]) + numbers = { + "mean": plain.mean, + "std": plain.std, + "tails mean": tails.mean, + "tails std": tails.std, + "weighted mean": weighted.mean, + "weighted std": weighted.std, + } + for scope in ("image", "biological_repeat"): + parameters = fit(_records(), NormalizationSettings(method="zscore", scope=scope)) + numbers[f"fit {scope}"] = repr(parameters.to_columns()) + return {k: v.hex() if isinstance(v, float) else v for k, v in numbers.items()} + + +def _exact(weights: np.ndarray, values: np.ndarray) -> float: + """The correctly rounded sum of the float64 products, with exact fractions.""" + products = np.multiply(weights, values, dtype=np.float64) + return float(sum(Fraction(p) for p in products.tolist())) + + +def test_pooled_sums_are_correctly_rounded() -> None: + """The mean and spread sums equal the exact sums of their products, rounded once.""" + pooled = PooledHistogram(_records()) + total = pooled.pixel_count + assert pooled.mean == _exact(pooled.counts, pooled.values) / total + deviations = pooled.values - pooled.mean + assert pooled.std == float(np.sqrt(_exact(pooled.counts, deviations * deviations) / total)) + + +def test_weighted_sum_is_exactly_rounded_and_order_free() -> None: + """``weighted_sum`` is the correctly rounded sum, whatever the order of the terms.""" + rng = np.random.default_rng(3) + counts = rng.integers(1, 10**6, 50_000) + values = rng.normal(0.0, 1e4, 50_000) + assert weighted_sum(counts, values) == _exact(counts, values) + order = rng.permutation(counts.size) + assert weighted_sum(counts[order], values[order]) == weighted_sum(counts, values) + + +def test_fits_are_bitwise_equal_with_one_and_many_threads(tmp_path: Path) -> None: + """Every pooled sum and z-score fit has the same bits with 1 and 4 BLAS/OpenMP threads.""" + script = tmp_path / "fit_bits.py" # not numbers.py, which would hide the standard module + script.write_text( + "import json, sys\n" + f"sys.path.insert(0, {str(Path(__file__).parent)!r})\n" + "from test_sums import _numbers\n" + "print(json.dumps(_numbers()))\n", + encoding="utf-8", + ) + results = [] + for threads in ("1", "4"): + environment = {**os.environ, **dict.fromkeys(THREAD_VARIABLES, threads)} + done = subprocess.run( + [sys.executable, str(script)], + env=environment, + capture_output=True, + text=True, + check=False, + ) + assert done.returncode == 0, done.stderr + results.append(json.loads(done.stdout)) + assert results[0] == results[1] + assert results[0] == _numbers() From 753bcc17a78927c3828312646c99e24a24cb1936 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 5 Nov 2026 08:07:39 +0000 Subject: [PATCH 2/2] Keep nan and inf for non-finite weighted sums instead of a bare fsum error (#67) math.fsum raises ValueError on inf + -inf and OverflowError on an overflowing sum, where np.dot returned nan or inf and fit refused the group by name. Fall back to NumPy's single-threaded sum in those cases. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WX6rdMV1r5jUPUR94kxV5a --- README.md | 3 ++- src/morphobricks_normalization/sums.py | 12 ++++++++++-- tests/test_sums.py | 25 +++++++++++++++++++++++++ 3 files changed, 37 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 3708d35..1121519 100644 --- a/README.md +++ b/README.md @@ -179,7 +179,8 @@ warnings = range_check_issues([("img.tif::z=0,t=0", 0, check)]) platform, never less accurate than before. **The last bits change:** pooled `mean` and `std` and the `zscore` offset and scale differ from v0.3.3 by at most 3.7e-14 relative (134 numbers on 86 of the 1,788 lines of `tools/dump_outputs.py`; 89 by 1 ulp, the largest - 317 ulps, for int64 planes near 2**53). Every other output is bitwise unchanged. + 317 ulps, for int64 planes near 2**53). Every other output is bitwise unchanged. A sum + that overflows float64 is still `nan` or `inf`, refused as before. ### Migrating to v0.3.3 diff --git a/src/morphobricks_normalization/sums.py b/src/morphobricks_normalization/sums.py index 32c0140..95efe28 100644 --- a/src/morphobricks_normalization/sums.py +++ b/src/morphobricks_normalization/sums.py @@ -29,7 +29,15 @@ def weighted_sum(weights: np.ndarray, values: np.ndarray) -> float: Returns ------- total : float - The correctly rounded sum of the float64 products. + The correctly rounded sum of the float64 products. When that sum is not finite (an + infinite product of each sign, or finite products whose sum overflows), the IEEE + result ``nan`` or ``inf`` that ``np.dot`` gave, never an exception. """ products = np.multiply(weights, values, dtype=np.float64) - return math.fsum(products.tolist()) + try: + return math.fsum(products.tolist()) + except (OverflowError, ValueError): + # fsum raises on inf + -inf and on an overflowing sum, where np.dot returned nan or + # inf for the callers to refuse by name. NumPy's sum of a contiguous float64 array is + # one thread in a fixed order, so the fallback is deterministic too. + return float(np.add.reduce(products)) diff --git a/tests/test_sums.py b/tests/test_sums.py index b3f910e..f32c2ec 100644 --- a/tests/test_sums.py +++ b/tests/test_sums.py @@ -111,3 +111,28 @@ def test_fits_are_bitwise_equal_with_one_and_many_threads(tmp_path: Path) -> Non results.append(json.loads(done.stdout)) assert results[0] == results[1] assert results[0] == _numbers() + + +def test_non_finite_sums_stay_nan_or_inf_as_with_np_dot() -> None: + """``math.fsum`` raises on ``inf - inf`` and overflow; the sums keep np.dot's nan and inf.""" + big = np.finfo(np.float64).max + with np.errstate(over="ignore", invalid="ignore"): + assert np.isnan(weighted_sum(np.array([1, 1]), np.array([np.inf, -np.inf]))) + assert np.isnan(weighted_sum(np.array([4, 4]), np.array([big, -big]))) + assert weighted_sum(np.array([1, 1]), np.array([big, big * 0.9])) == np.inf + assert weighted_sum(np.array([2, 1]), np.array([-np.inf, 1.0])) == -np.inf + assert np.isnan(weighted_sum(np.array([1, 1]), np.array([np.nan, 1.0]))) + assert weighted_sum(np.array([], dtype=np.int64), np.array([])) == 0.0 + + +def test_huge_planes_are_refused_by_name_not_by_fsum() -> None: + """A float64 plane whose sums overflow gets the z-score refusal, not a bare fsum error.""" + big = np.finfo(np.float64).max + for plane in (np.array([[big, -big], [big, -big]]), np.array([[big, big * 0.9], [1.0, 2.0]])): + with np.errstate(over="ignore", invalid="ignore"): + try: + fit(plane, NormalizationSettings(method="zscore")) + except ValueError as error: + assert "z-score normalization cannot be used" in str(error) + else: + raise AssertionError("expected the z-score refusal")