From 08ffc0b507596da1e41072959be0a9993b520475 Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 31 Aug 2026 18:06:10 -0500 Subject: [PATCH 1/2] fix(plot2d): size the box-reduce accumulator to the data, not to uint32 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _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) --- anyplotlib/plot2d/_tile_backend.py | 57 +++++++--- .../tests/test_plot2d/test_tile_backend.py | 101 +++++++++++++++++- 2 files changed, 142 insertions(+), 16 deletions(-) diff --git a/anyplotlib/plot2d/_tile_backend.py b/anyplotlib/plot2d/_tile_backend.py index be68a0fb5..954a2918d 100644 --- a/anyplotlib/plot2d/_tile_backend.py +++ b/anyplotlib/plot2d/_tile_backend.py @@ -59,30 +59,57 @@ def sample(self, x0: int, x1: int, y0: int, y1: int, # ── Integer-friendly area reductions (shared by the numpy backend) ───────────── +def _acc_dtype(dtype: np.dtype, n: int) -> np.dtype: + """Accumulator dtype for a box sum of ``n`` elements of ``dtype``. + + Integers get the narrowest same-SIGNEDNESS integer that cannot wrap over ``n`` + terms. A blanket ``uint32`` is wrong twice over: it reads a negative ``int16`` + as ~4.3e9, and it silently overflows on anything wider than 16 bits. Floats + accumulate at their own width but never below float32 — the result is quantised + to 8-bit tile bytes, so a float64 accumulator over a float32 frame buys + precision the output cannot express and pays a cast of the whole region for it. + + The float32 sum does bound what a float32 frame may hold: ``n`` terms overflow + to inf past ``float32.max / n`` (~5e36 for the default 64-pixel block). A frame + that large cannot survive the 8-bit quantisation downstream either, and float64 + input — where that dynamic range actually turns up — keeps its own width.""" + if dtype.kind in "bui": + hi = 1 if dtype.kind == "b" else int(np.iinfo(dtype).max) + lo = -int(np.iinfo(dtype).min) if dtype.kind == "i" else 0 + need = int(n) * max(hi, lo) + for cand in (np.int32, np.int64) if dtype.kind == "i" else (np.uint32, np.uint64): + if need <= np.iinfo(cand).max: + return np.dtype(cand) + return np.dtype(np.float64) # absurdly large block; float64 still beats wrapping + return np.promote_types(dtype, np.float32) + + def _box_reduce(region: np.ndarray, out_h: int, out_w: int, op: str) -> np.ndarray: """Reduce ``region`` (2-D) to ``(out_h, out_w)`` by a block ``op`` ("mean"|"max"). - Fast path (region divisible by the stride): a single VECTORISED reshape-reduce in - a wide integer accumulator (no full float cast of the source — the cast of a 16 MP - uint16 frame alone is ~34 ms). Ragged path (non-divisible): a strided-accumulate - box filter with a per-cell count, so the last partial block is reduced over only - its valid pixels. Either way the grid is nearest-resized to the exact (out_h, - out_w) the caller asked for.""" + Fast path (region divisible by the stride): TWO vectorised reshape-reduces in an + accumulator sized by :func:`_acc_dtype` (no full float cast of the source — the + cast of a 16 MP uint16 frame alone is ~34 ms). Collapsing whole ROWS first and + then the contiguous column blocks walks memory in order on both passes; the one + ``sum(axis=(1, 3))`` it replaces reduced a strided axis and a contiguous one + together, which defeats numpy's fast inner loops and cost ~4x as long on an + 8192² float32 frame. Ragged path (non-divisible): a strided-accumulate box filter + with a per-cell count, so the last partial block is reduced over only its valid + pixels. Either way the grid is nearest-resized to the exact (out_h, out_w) the + caller asked for.""" h, w = region.shape sy = max(1, h // out_h) sx = max(1, w // out_w) - is_int = np.issubdtype(region.dtype, np.integer) if h % sy == 0 and w % sx == 0: - # Divisible → one reshape-reduce (fast, vectorised). Integer sum stays uint32 - # (uint16 × up-to-64 block fits) so there's no giant float cast. gh, gw = h // sy, w // sx - blk = region.reshape(gh, sy, gw, sx) if op == "max": - out = blk.max(axis=(1, 3)) - else: - out = (blk.sum(axis=(1, 3), dtype=np.uint32 if is_int else np.float64) - .astype(np.float32) / (sy * sx)) + return _nearest_resize(region.reshape(gh, sy, gw, sx).max(axis=(1, 3)), + out_h, out_w) + acc = _acc_dtype(region.dtype, sy * sx) + rows = region.reshape(gh, sy, w).sum(axis=1, dtype=acc) + out = (rows.reshape(gh, gw, sx).sum(axis=2, dtype=acc) + .astype(np.float32) / (sy * sx)) return _nearest_resize(out, out_h, out_w) # Ragged → strided accumulate with a per-cell count (handles the partial block). @@ -101,7 +128,7 @@ def _box_reduce(region: np.ndarray, out_h: int, out_w: int, op: str) -> np.ndarr np.maximum(acc[:sh, :sw], sub, out=acc[:sh, :sw]) out = acc else: - acc = np.zeros((gh, gw), np.uint32 if is_int else np.float64) + acc = np.zeros((gh, gw), _acc_dtype(region.dtype, sy * sx)) cnt = np.zeros((gh, gw), np.uint32) for dy in range(sy): for dx in range(sx): diff --git a/anyplotlib/tests/test_plot2d/test_tile_backend.py b/anyplotlib/tests/test_plot2d/test_tile_backend.py index e8e055e6e..7d2da7d44 100644 --- a/anyplotlib/tests/test_plot2d/test_tile_backend.py +++ b/anyplotlib/tests/test_plot2d/test_tile_backend.py @@ -3,7 +3,7 @@ import pytest from anyplotlib.plot2d._tile_backend import ( - NumpyTileBackend, TileBackend, as_tile_backend, + NumpyTileBackend, TileBackend, _acc_dtype, as_tile_backend, ) @@ -50,6 +50,105 @@ def test_subregion(self): np.testing.assert_allclose(out, ref, rtol=1e-4) +class TestSampleDtypes: + """A box mean must mean the same thing in every dtype the caller can hand us. + + The accumulator used to be a blanket ``uint32`` for anything integral, which + wrapped every negative value and overflowed anything wider than 16 bits.""" + + @pytest.mark.parametrize("dt", [np.int8, np.int16, np.int32, np.int64]) + @pytest.mark.parametrize("shape,out", [((64, 64), 16), ((63, 65), 16)]) + def test_negative_integers_survive_the_mean(self, dt, shape, out): + # Both grids on purpose: (64, 64) takes the divisible reshape-reduce, + # (63, 65) the ragged strided accumulate. uint32 got the first one silently + # wrong (-100 read as ~4.29e9) and made the second raise. + a = np.full(shape, -100, dt) + got = NumpyTileBackend(a).sample(0, shape[1], 0, shape[0], out, out, "mean") + np.testing.assert_allclose(got, -100.0) + + def test_mixed_sign_block_averages_to_zero(self): + a = np.tile(np.array([-50, 50], np.int16), (32, 32)) + out = NumpyTileBackend(a).sample(0, 64, 0, 64, 16, 16, "mean") + np.testing.assert_allclose(out, 0.0, atol=1e-5) + + @pytest.mark.parametrize("dt,val", [(np.int32, 2_000_000), + (np.uint32, 4_000_000_000), + (np.int64, 3_000_000_000), + (np.uint64, 5_000_000_000)]) + def test_wide_integers_do_not_overflow(self, dt, val): + # 64 terms of these all exceed uint32; the sum has to widen, not wrap. + a = np.full((64, 64), val, dt) + out = NumpyTileBackend(a).sample(0, 64, 0, 64, 8, 8, "mean") + np.testing.assert_allclose(out, float(val), rtol=1e-6) + + def test_bool_mean_is_the_fraction_set(self): + a = np.zeros((16, 16), bool) + a[::2] = True # every other row + out = NumpyTileBackend(a).sample(0, 16, 0, 16, 4, 4, "mean") + np.testing.assert_allclose(out, 0.5) + + @pytest.mark.parametrize("dt", [np.float16, np.float32, np.float64, + np.uint8, np.uint16, np.int16, np.int64, np.bool_]) + @pytest.mark.parametrize("shape,out", [((64, 64), 16), ((128, 128), 32)]) + def test_mean_matches_an_exact_block_mean(self, dt, shape, out): + rs = np.random.RandomState(0) + h, w = shape + if dt is np.bool_: + a = rs.rand(h, w) > 0.5 + elif np.issubdtype(dt, np.integer): + info = np.iinfo(dt) + a = rs.randint(max(info.min, -20000), min(info.max, 20000), shape).astype(dt) + else: + a = ((rs.rand(h, w) - 0.5) * 2000).astype(dt) + got = NumpyTileBackend(a).sample(0, w, 0, h, out, out, "mean") + sy, sx = h // out, w // out + ref = a.astype(np.float64).reshape(out, sy, out, sx).mean(axis=(1, 3)) + # A mixed-sign block can average to ~0, so the bound that matters is + # absolute and set by the SOURCE magnitude, not by the mean it lands on. + np.testing.assert_allclose(got, ref, rtol=2e-6, + atol=1e-5 * float(np.abs(a).max())) + + def test_float64_input_is_not_narrowed_before_the_sum(self): + # The mean is returned as float32 either way, so the extra mantissa is not + # observable — the RANGE is. 16 terms of 1e37 overflow a float32 + # accumulator to inf; the mean itself is an ordinary float32. + a = np.full((64, 64), 1e37, np.float64) + out = NumpyTileBackend(a).sample(0, 64, 0, 64, 16, 16, "mean") + assert np.isfinite(out).all() + np.testing.assert_allclose(out, 1e37, rtol=1e-6) + + @pytest.mark.parametrize("dt", [np.uint16, np.int16, np.float32]) + def test_ragged_grid_averages_over_valid_pixels_only(self, dt): + # 100 // 32 = 3, so the last block of each axis is partial: the count has to + # follow the block, not assume a full 3x3. + a = np.arange(100 * 100, dtype=np.float64).reshape(100, 100) % 1000 + a = a.astype(dt) + got = NumpyTileBackend(a).sample(0, 100, 0, 100, 32, 32, "mean") + ref = np.array([[a[i * 3:(i + 1) * 3, j * 3:(j + 1) * 3].astype(np.float64).mean() + for j in range(34)] for i in range(34)]) + yi = (np.arange(32) * 34 // 32).clip(0, 33) + np.testing.assert_allclose(got, ref[yi][:, yi], rtol=2e-6) + + +class TestAccumulatorDtype: + @pytest.mark.parametrize("dt,expect", [ + (np.uint8, np.uint32), (np.uint16, np.uint32), (np.bool_, np.uint32), + (np.int8, np.int32), (np.int16, np.int32), + (np.float16, np.float32), (np.float32, np.float32), (np.float64, np.float64), + ]) + def test_narrow_dtypes_keep_a_cheap_same_signedness_accumulator(self, dt, expect): + assert _acc_dtype(np.dtype(dt), 64) == np.dtype(expect) + + def test_accumulator_widens_with_the_block(self): + # uint16 x 64 fits uint32; uint16 x 2**24 does not. + assert _acc_dtype(np.dtype(np.uint16), 64) == np.dtype(np.uint32) + assert _acc_dtype(np.dtype(np.uint16), 2 ** 24) == np.dtype(np.uint64) + + def test_signedness_is_preserved(self): + for dt in (np.int8, np.int16, np.int32, np.int64): + assert _acc_dtype(np.dtype(dt), 64).kind in "if" + + class TestSampleSubsampleMax: def test_subsample_drops_between_grid(self): a = np.zeros((16, 16), np.float32) From 5d2fe3947d375b77040b80245a179fe993f18ddc Mon Sep 17 00:00:00 2001 From: Carter Francis Date: Mon, 31 Aug 2026 18:06:57 -0500 Subject: [PATCH 2/2] docs: changelog fragment for the box-reduce accumulator fix MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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) --- upcoming_changes/65.bugfix.rst | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 upcoming_changes/65.bugfix.rst diff --git a/upcoming_changes/65.bugfix.rst b/upcoming_changes/65.bugfix.rst new file mode 100644 index 000000000..b2867697b --- /dev/null +++ b/upcoming_changes/65.bugfix.rst @@ -0,0 +1,6 @@ +Fixed a tiled :meth:`~anyplotlib.Axes.imshow` of a **signed** integer frame +displaying wrapped values — a ``uint32`` accumulator in the overview box-mean read +an ``int16`` ``-100`` as ``+1073741696``, raised +``_UFuncOutputCastingError`` on a non-divisible grid, and overflowed on any dtype +wider than 16 bits; the accumulator is now sized to the data, which also makes the +overview up to 4x faster on large float frames.