[guijacquemet]

Fixes #67.

Problem. PooledHistogram/WeightedHistogram mean and std (and so every zscore fit, at every scope) summed with np.dot. OpenBLAS threads a dot product above 10,000 elements and adds the partial sums in a thread-dependent order, so the fitted offset/scale changed in their last bits with OMP_NUM_THREADS/OPENBLAS_NUM_THREADS (MorphoFeature #130 finding: MorphoEmbed's per-plane z-score fit gave …030 vs …531).

Fix. New pure module sums.py: weighted_sum(weights, values) multiplies element by element (each product rounded once, as np.dot did) and adds the products with math.fsum, the correctly rounded sum. Same bits for any thread count, BLAS, NumPy version or platform; never less accurate than before. No new dependency.

Every thread-dependent reduction (audit): all seven np.dot calls in src/ go through it: PooledHistogram.mean/std, WeightedHistogram.mean/std, _mean_with_tails, _std_with_tails, and histogram._table_tail (the per-value tail sums). The remaining reductions (np.sum pairwise sums in _TailSums/integrity, cumsum, bincount, add.at, unique) are single-threaded NumPy loops whose order does not depend on threads. No matmul/einsum/linalg/scipy use.

Scientific precision (bits change). python tools/dump_outputs.py: digest d82a9f97… on main → dfd5e740… here. 86 of 1,788 lines differ, all pooled mean/std or zscore offset/scale; 134 numbers in total, max relative difference 3.7e-14; ulp histogram: 1 ulp ×89, 2 ×12, 3 ×1, 4 ×3, 5 ×11, 9 ×2, 10 ×1, 16 ×2, 17 ×2, 41 ×1, 52 ×3, 76 ×6, 317 ×1 (the large ones are int64 planes whose values are near 2**53). Every other output (percentiles, median, MAD, min/max, histograms, CSV text, other methods) is bitwise unchanged. The new values are the correctly rounded sums of the same products, so they are at least as accurate as the old ones. README "Unreleased" and contract §2 table say so.

Tests (tests/test_sums.py):

• test_fits_are_bitwise_equal_with_one_and_many_threads: runs the pooled/weighted/tail mean and std and zscore fits (image and biological_repeat scopes, >10,000 pooled values) in two subprocesses with OMP/OPENBLAS/MKL_NUM_THREADS = 1 and 4; asserts identical hex. Fails on main (tails std …d07f vs …d07e).
• test_weighted_sum_is_exactly_rounded_and_order_free: equals an exact Fraction sum, rounded once, and is permutation-invariant (fails against an np.dot stand-in: …462402 vs …462654).
• test_pooled_sums_are_correctly_rounded.

Gate: python tools/check.py (Python 3.12): exit 0; 630 passed, 2 skipped (csbdeep not installed: tests/test_rgb.py:78). Code-rules baseline unchanged (new module sums declared pure).

Release: v0.3.4 is not tagged. Merge this into main before the release PR #68; #68 must then be brought up to date with main (its README turns "Unreleased" into the v0.3.4 section, so it must carry this entry and the migration notes must say the z-score/mean/std bits change, instead of "valid inputs give bit-identical results"). MorphoFeature's pin then moves to v0.3.4 with the finals.

---

Generated by Claude Code

[guijacquemet]

Independent review: merged as b3999de.

Finding fixed (753bcc1): math.fsum raises ValueError on inf + -inf and OverflowError on an overflowing sum, where np.dot returned nan/inf and fit refused the group by name. Example: a float64 plane of ±1.7e308 now ended in the bare error "-inf + inf in fsum". weighted_sum now falls back to np.add.reduce (single-threaded, fixed order) in those cases only. There are two regression tests, both failing on 11bf51a. The dump digest is unchanged (dfd5e740…).

Verified:

• dtypes: values are always float64 (representatives, tail means); weights are int64 counts or float64 weights. np.multiply(..., dtype=float64) casts before multiplying, so no integer overflow (the same products as np.dot's int64·float64 promotion). The result for empty inputs is 0.0.
• No other BLAS call: there is no dot, matmul, einsum, linalg or scipy use in src/.
• The 1-vs-4-thread test fails on main (tails std …d07f vs …d07e).
• tools/dump_outputs.py reproduced: 86 lines, 134 numbers, max relative 3.70e-14, the same ulp histogram.
• Timing (3.12, OpenBLAS): a 65,536-term sum takes 2.5 ms (np.dot 0.05 ms, np.add.reduce 0.08 ms). Per plane with all 65,536 uint16 values: mean+std 4.0 ms vs 0.27 ms, against 22 ms for the histogram and 47 ms for fit(pixels). A 12-bit plane costs +0.2 ms. Not slow enough to give up correct rounding, so fsum stays.
• Gate: tools/check.py exit 0 (632 passed, 2 skipped: csbdeep). CI is green on 3 OSes × 3.11/3.13.

Follow-up (pre-existing on main, not this PR): a non-finite z-score scale is refused with "its intensity spread is zero (all pixels have the same value)", which is wrong for an overflow.

---

Generated by Claude Code