Make pooled means and standard deviations independent of CPU threads (#67)

README.md

@@ -172,6 +172,15 @@ 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. A sum
+  that overflows float64 is still `nan` or `inf`, refused as before.
 
 ### Migrating to v0.3.3
 

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 |

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"

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

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.

src/morphobricks_normalization/sums.py

@@ -0,0 +1,43 @@
+"""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. 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)
+    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))

tests/test_sums.py

@@ -0,0 +1,138 @@
+"""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()
+
+
+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")