_box_reduce's uint32 accumulator wraps every signed integer frame

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

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)

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.