_box_reduce's uint32 accumulator wraps every signed integer frame - #65
_box_reduce's uint32 accumulator wraps every signed integer frame#65CSSFrancis wants to merge 2 commits into
Conversation
_box_reduce picked its accumulator with `np.uint32 if is_int else np.float64`.
Both halves of that are wrong, and the integer half is a correctness bug.
Any SIGNED integer frame wraps. A 2048² int16 image of -100 comes back from the
divisible path as +1073741696 — imshow displays it, no warning, no clue:
a = np.full((2048, 2048), -100, np.int16)
ax.imshow(a) # overview reads 1073741696.0
The ragged path does not even get that far; `acc[:sh, :sw] += sub` raises
_UFuncOutputCastingError ("cannot cast ufunc 'add' output from int64 to uint32")
for the same frame on a non-divisible grid. Anything wider than 16 bits
overflows regardless of sign — an int64 frame of 3e9 reads back 4.7e7. The
tests only ever covered uint16 and float32, which is why none of this showed.
_acc_dtype now picks the narrowest accumulator that cannot wrap over the block:
same signedness, widened by the block size, falling back to float64 for the
absurd cases. bool keeps the cheap uint32 path.
The float half was only slow, and less slow than #64 supposed. float64 there
does not materialise a float64 copy of the region — numpy buffers the cast — so
switching to float32 alone is worth ~1.1x, not the ~2x the issue predicts. The
real cost is the reduction ORDER: sum(axis=(1, 3)) collapses a strided axis and
a contiguous one in one pass, which defeats numpy's fast inner loops.
Collapsing whole rows first and then the contiguous column blocks walks memory
in order both times. The two changes compound — fix the access pattern and the
cast becomes the bottleneck, and vice versa:
8192² float32 -> 1024² overview 64.6ms -> 15.7ms 4.13x
8192² float64 57.0ms -> 26.3ms 2.17x
8192² uint16 66.3ms -> 41.6ms 1.59x
float32 accumulation costs 2.4e-7 relative error against an output quantised to
1/255, so the overview cannot express the difference. It does bound a float32
frame at ~5e36 before the sum overflows to inf; such a frame cannot survive the
8-bit quantisation downstream either, and float64 input keeps float64.
16 of the new dtype tests fail on the previous implementation.
Assisted-by: Claude Opus 5 (1M context)
Named for the PR since it has one, per the upcoming_changes README — the
orphan +{slug} form is for work batched on a branch with no PR number.
Verified with `towncrier build --draft`.
Assisted-by: Claude Opus 5 (1M context)
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #65 +/- ##
==========================================
+ Coverage 90.90% 90.92% +0.01%
==========================================
Files 40 40
Lines 4550 4560 +10
==========================================
+ Hits 4136 4146 +10
Misses 414 414 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
Benchmarked on our Windows box to close the loop on the absolute-numbers caveat — as Machine: Windows 10, Xeon E5-1650 v2 workstation (6c/3.5 GHz), numpy 2.5.2, Python 3.12.10.
Per-rep spread was tight — the float32 8192² case ran 162–182 ms before and 53–61 ms Signed-integer bug confirmed here too, both grid paths, (With numpy 2.5.2 the ragged crash surfaces as On the numbers being different from Apple silicon: our before for float32 8192² Output equivalence, same fixed-seed frames: uint8/uint16/float64 overviews are |
3 participants
Closes #64.
Why
_box_reducechose its accumulator with one line:#64 is about the
np.float64half. Thenp.uint32half turned out to be worse:it is a silent correctness bug on every signed integer frame, and it has been
there as long as the tiled path has.
A 2048²
int16image of-100— a difference map, a background-subtractedframe, an ordinary signed detector readout — displays as
+1073741696:The ragged path does not even get that far. On a non-divisible grid the same
frame raises out of
acc[:sh, :sw] += sub:And the width is wrong independently of the sign —
uint32overflows anythingpast 16 bits, so an
int64frame of 3e9 reads back as 4.7e7.test_tile_backend.pyonly ever exerciseduint16andfloat32. That is thewhole reason this survived.
What
_acc_dtype(dtype, n)picks the narrowest accumulator that cannot wrap over ann-pixel block: same signedness, widened by the block size, falling back tofloat64 for absurdly large blocks.
boolkeeps the cheapuint32path.uint8,uint16,booluint32(unchanged)int8,int16int32uint32/int64uint64/float64float16,float32float32float64float64The float half — the issue's premise is right, its mechanism is not
#64 attributes the cost to "~537 MB of accumulator traffic". That is not what
happens: numpy does not materialise a float64 copy of the region for
sum(dtype=np.float64), it buffers the cast in chunks. The accumulator is the1024² output grid, 8 MB. So the proposed fix on its own underdelivers — measured
on an 8192² float32 frame:
sum(axis=(1,3), dtype=float64)— beforesum(axis=(1,3), dtype=float32)— as proposed in #64The dominant cost is the reduction order.
sum(axis=(1, 3))collapses astrided axis and a contiguous one in a single pass, which defeats numpy's fast
inner loops. Collapsing whole rows first and then the contiguous column blocks
walks memory in order on both passes:
The two changes compound rather than add (1.1x and 1.5x alone, 4.1x together) —
once the access pattern is contiguous the cast is the bottleneck, and vice versa.
Measured
Apple silicon, numpy 2.4.4. @TheDrOnos measured ~170 ms for the same call on
two other machines, so the absolute numbers travel less well than the ratios.
The trade this makes
float32 accumulation costs 2.4e-7 relative error against an output quantised to
1/255 ≈ 4e-3, so the overview cannot express the difference — the issue's core
argument, and it holds.
It does bound what a float32 frame may contain:
nterms overflow to inf pastfloat32.max / n, ~5e36 for the default 64-pixel block. A frame carrying valuesthat large cannot survive the 8-bit quantisation downstream either, and float64
input — where that dynamic range actually turns up — keeps its own width. The
bound is written down in
_acc_dtype's docstring rather than left to berediscovered.
Tests
A
TestSampleDtypesclass covering every integer width and both grid paths(divisible and ragged), plus
TestAccumulatorDtypefor the choice itself.16 of them fail on the previous implementation; all 59 pass now, and the
full suite is green (2090 passed, 58 skipped).