Scan M4 decimation buckets with SIMD argmin/argmax

Cargo.toml

@@ -8,12 +8,15 @@ license = "Apache-2.0"
 # The compute core remains exposed through a plain C ABI, so one cdylib per
 # platform covers every CPython version. Static chart export uses the focused
 # `png` crate's fdeflate fast path; compression belongs beside rasterization so
-# the hot path never round-trips a full framebuffer through Python.
+# the hot path never round-trips a full framebuffer through Python. M4
+# decimation scans buckets of 32+ rows in dense windows with `argminmax`'s
+# runtime-dispatched SIMD argmin/argmax (spec/design/rust-engine.md §3.4).
 [lib]
 name = "xy_core"
 crate-type = ["cdylib", "rlib"]
 
 [dependencies]
+argminmax = { version = "0.7.1", default-features = false, features = ["float"] }
 png = "0.18.1"
 
 # Whole-program optimization; `strip` keeps .dynsym, which ctypes binds

news/535.performance.md

@@ -0,0 +1 @@
+Long line and area charts decimate faster. When pixel columns average 8 or more rows, M4 now scans one column's rows at a time instead of updating a per-row state machine, with SIMD (the `argminmax` crate) for columns of 32 rows or more. Deep zooms keep the per-row loop. On one core, the decimation kernel is 3–9× faster at 1M points, the first payload of a 1M-point line chart is prepared 1.2–1.6× faster, and zooming into 10% of it re-decimates about 2× faster. The decimated rows do not change.

spec/design-dossier.md

@@ -1228,7 +1228,8 @@ complexity lives almost entirely here, not in rendering. Miss any piece and the
 hits a source build requiring a Rust toolchain — an instant adoption cliff.
 
 1. **The native core as prebuilt wheels.** Built as a Rust **`cdylib` with a
-   plain C ABI** and a focused PNG-encoding dependency, compiled by the Hatchling build hook and loaded
+   plain C ABI** and two focused dependencies (PNG encoding, and `argminmax` for M4's SIMD
+   argmin/argmax), compiled by the Hatchling build hook and loaded
    from Python with `ctypes`. There is no CPython extension ABI at all, so one
    `py3-none-<platform>` wheel covers every supported Python version on that platform
    without PyO3 or `abi3`. Wheel matrix in CI, release-blocking: manylinux

spec/design/rust-engine.md

@@ -1,8 +1,9 @@
 # Rust Data Engine — Module Boundaries & FFI Protocol
 
 **Status:** design. Decides what lives in Rust vs Python and how the FFI seam
-evolves without rewrites. Grounded in the shipped engine: single-dependency
-cdylib (`src/lib.rs`, ABI v44; `png` is the one crate, for static export),
+evolves without rewrites. Grounded in the shipped engine: two-dependency
+cdylib (`src/lib.rs`, ABI v44; `png` for static export, `argminmax` for M4's
+SIMD argmin/argmax),
 ctypes binding (`_native.py`), and
 dispatch in `kernels.py`. The native core is required — `kernels.py` raises a
 clear ImportError when it can't load, with no pure-Python fallback.
@@ -88,7 +89,7 @@ src/                                          # 15,423 lines shipped, 8 modules
                         #     · row-banded multithreaded paint (std::thread::scope,
                         #       byte-identical to the serial path) and the fused
                         #       raster→PNG encode via `png`/fdeflate.
-                        #   No unsafe; owns the crate's one third-party dep.
+                        #   No unsafe; owns the `png` dependency.
   font.rs       (2503)  # generated by scripts/gen_font.py — do not hand-edit.
                         #   Baked DejaVu Sans grayscale coverage atlas for
                         #   raster.rs: 424 glyph records (advance, w, h, left,
@@ -106,7 +107,8 @@ src/                                          # 15,423 lines shipped, 8 modules
                         #   input by returning None. Deliberately narrow: it
                         #   exists to kill one Python string per point, and the
                         #   rest of SVG scene construction stays in Python.
-  simd.rs        (448)  # AVX2 twins of eligible kernels, runtime-dispatched (§3.4).
+  simd.rs        (448)  # AVX2 twins of eligible kernels, runtime-dispatched (§3.4),
+                        #   and `try_argminmax` (the `argminmax` crate, for M4).
                         #   The one module besides lib.rs allowed `unsafe`.
   tiles.rs      (1050)  # pyramid build/compose/incremental append, plus the
                         #   mean-color planes for channel-bearing traces
@@ -130,6 +132,11 @@ the C-ABI/one-cdylib-per-platform property is preserved. Note the dev sandbox
 cannot reach crates.io, so required crates must be vendored or the sandbox
 loses local build/test; prefer feature-gated optional deps (e.g. SIMD
 argminmax, tsdownsample-class speed) with the lean build as default.
+`argminmax` (M4, §3.4) is a plain dependency, not a gated one: the release
+wheels build xy-core with its default Cargo features (`hatch_build.py` runs a
+plain `cargo build --release`), so a gated crate would not reach users. It is
+used with `default-features = false, features = ["float"]`; its tree is
+`argminmax` + `num-traits` (+ `autocfg` at build).
 
 Release profile (`Cargo.toml`, 2026-07-27): `lto = "fat"`, `codegen-units =
 1`, `strip = true`. Fat LTO measured no runtime change over thin on the §12
@@ -306,7 +313,9 @@ registers unless we ask. `simd.rs` holds branch-free clones of selected
 kernels compiled under `#[target_feature(enable = "avx2")]` (LLVM
 autovectorizes them — no hand-written intrinsics unless a loop demonstrably
 fails to vectorize) with runtime `is_x86_feature_detected!` dispatch and a
-`XY_SIMD=0` kill switch.
+`XY_SIMD=0` kill switch. It also wraps the `argminmax` crate
+(`try_argminmax`), which dispatches at run time on its own (SSE4.1/AVX2/
+AVX-512 on x86, NEON on aarch64); the kill switch covers it too.
 
 Rules, in priority order:
 
@@ -317,14 +326,21 @@ Rules, in priority order:
    `simd.rs` comparing SIMD vs scalar on hostile data.
 2. **Only measured wins ship.** Every dispatch is justified by a before/after
    number (via the kill switch); a kernel where the two-phase restructure
-   loses (M4's sequential bucket state machine, histogram's scatter-dominated
-   loop) is documented at its scalar definition and NOT dispatched.
+   loses (histogram's scatter-dominated loop) is documented at its scalar
+   definition and NOT dispatched. It lost on M4's per-row bucket state machine
+   too, so M4 windows that average at least 8 rows per bucket
+   (`M4_BUCKET_WALK_MIN_ROWS`) walk whole buckets instead and hand the y of
+   each bucket of 32+ rows to `try_argminmax` (3–9× at 1M points, single
+   thread); shorter buckets take a scalar scan. Sparser
+   windows (deep zooms) keep the per-row loop, which is faster there. A fuzz
+   test in `kernels.rs` checks that both M4 paths emit the same rows.
 3. **`unsafe` containment.** This module is the one exception to "unsafe only
    in `lib.rs`": the `#[target_feature]` wrappers are unsafe to call, and
    every call site sits behind a safe `try_*` fn that checks detection first.
    Kernels code never writes `unsafe` — it calls `simd::try_*` and falls back.
 4. **aarch64 needs no twin.** NEON is part of the aarch64 baseline, so the
-   scalar kernels already autovectorize at full width there.
+   scalar kernels already autovectorize at full width there (argminmax picks
+   NEON itself).
 
 ## 4. What Python keeps forever
 

src/kernels.rs

@@ -2523,10 +2523,9 @@ pub(crate) fn m4_flush(first: u32, min_i: u32, max_i: u32, last: u32, out: &mut
 
 /// Serial M4 over rows `[lo, hi)` — the shared building block: the serial path
 /// runs it once over the whole window; the parallel path runs it per bucket-
-/// aligned segment and concatenates. Stays scalar deliberately: a two-phase
-/// SIMD restructure (precomputed bucket ids) measured *slower* (472 vs
-/// 915 Mpt/s at 1M) — the sequential bucket state machine dominates, and the
-/// extra block store/load pass outweighs vectorizing the cheap float math.
+/// aligned segment and concatenates. Windows whose buckets average at least
+/// `M4_BUCKET_WALK_MIN_ROWS` rows walk bucket by bucket; sparser ones (deep
+/// zooms) keep the per-row loop, which is faster there. Both emit the same rows.
 fn m4_range(
     x: &[f64],
     y: &[f64],
@@ -2535,6 +2534,38 @@ fn m4_range(
     n_buckets: usize,
     lo: usize,
     hi: usize,
+) -> Vec<u32> {
+    let bucket_of = |v: f64| (((v - x0) * inv_bucket_w) as usize).min(n_buckets - 1);
+    let spanned = bucket_of(x[hi - 1]) - bucket_of(x[lo]) + 1;
+    if hi - lo < M4_BUCKET_WALK_MIN_ROWS * spanned {
+        m4_range_rows(x, y, x0, inv_bucket_w, n_buckets, lo, hi)
+    } else {
+        m4_range_buckets(x, y, x0, inv_bucket_w, n_buckets, lo, hi)
+    }
+}
+
+/// Average rows per bucket from which the bucket walk beats the per-row loop.
+/// At 4–5 rows its per-bucket search and call made it up to 1.3× slower on
+/// Tiger Lake (Zen 3 gained); at 8 rows it is 1.1–1.4× faster on both (1M
+/// rows, 2048 buckets, one core).
+const M4_BUCKET_WALK_MIN_ROWS: usize = 8;
+
+/// Per-row M4 state machine over rows `[lo, hi)`: one pass, one bucket id per
+/// row. Stays scalar deliberately: a two-phase SIMD restructure (precomputed
+/// bucket ids) measured *slower* (472 vs 915 Mpt/s at 1M) — the sequential
+/// bucket state machine dominates, and the extra block store/load pass
+/// outweighs vectorizing the cheap float math. Kept out of line: inlined next
+/// to the bucket walk, its loop lost registers to stack reloads (+8% on
+/// Tiger Lake for deep zooms).
+#[inline(never)]
+fn m4_range_rows(
+    x: &[f64],
+    y: &[f64],
+    x0: f64,
+    inv_bucket_w: f64,
+    n_buckets: usize,
+    lo: usize,
+    hi: usize,
 ) -> Vec<u32> {
     let mut out: Vec<u32> = Vec::with_capacity((n_buckets * 4).min((hi - lo) * 4));
 
@@ -2583,7 +2614,112 @@ fn m4_range(
     out
 }
 
-fn m4_bucket_end<F>(values: &[f64], bucket: usize, bucket_of: &F) -> usize
+/// Bucket-at-a-time M4 over rows `[lo, hi)`. x is ascending, so every
+/// bucket's rows are contiguous: `m4_bucket_end` probes x for where a bucket
+/// stops and `m4_bucket` scans its y, so x is only probed and y is read once.
+fn m4_range_buckets(
+    x: &[f64],
+    y: &[f64],
+    x0: f64,
+    inv_bucket_w: f64,
+    n_buckets: usize,
+    lo: usize,
+    hi: usize,
+) -> Vec<u32> {
+    let bucket_of = |v: f64| (((v - x0) * inv_bucket_w) as usize).min(n_buckets - 1);
+    let mut out: Vec<u32> = Vec::with_capacity((n_buckets * 4).min((hi - lo) * 4));
+    let (mut start, mut len) = (lo, 1);
+    while start < hi {
+        let b = bucket_of(x[start]);
+        let end = m4_bucket_end(x, start, hi, len, |v| bucket_of(v) <= b);
+        m4_bucket(y, start, end, &mut out);
+        (start, len) = (end, end - start);
+    }
+    out
+}
+
+/// End of the bucket that starts at row `start`: the first row in
+/// `(start, hi]` that fails `in_bucket` (monotone, as x is ascending). Probes
+/// around `start + guess` (the previous bucket's length) first, so
+/// equal-width buckets cost two probes; otherwise gallops, then bisects.
+fn m4_bucket_end(
+    x: &[f64],
+    start: usize,
+    hi: usize,
+    guess: usize,
+    in_bucket: impl Fn(f64) -> bool,
+) -> usize {
+    let g = start + guess;
+    // Rows [start, lo) are in the bucket; the end lies in [lo, end].
+    let (mut lo, mut end) = (start + 1, hi);
+    if g < hi {
+        if !in_bucket(x[g]) {
+            if in_bucket(x[g - 1]) {
+                return g;
+            }
+            end = g - 1;
+        } else {
+            lo = g + 1;
+            let mut step = 1;
+            while lo + step <= hi {
+                let p = lo + step - 1;
+                if !in_bucket(x[p]) {
+                    end = p;
+                    break;
+                }
+                lo = p + 1;
+                step *= 2;
+            }
+        }
+    }
+    lo + x[lo..end].partition_point(|&v| in_bucket(v))
+}
+
+/// Buckets shorter than this skip argminmax (uneven x leaves some in a dense
+/// window): its per-call cost outweighs the SIMD scan.
+const M4_SIMD_MIN_ROWS: usize = 32;
+
+/// Emit the {first, min, max, last} rows of bucket `[start, end)` over its
+/// finite y (NaN and ±∞ are non-plottable, §19). Ties keep the first row.
+/// `simd::try_argminmax` skips NaN but not ±∞, so a bucket with a ±∞ extreme
+/// takes the scalar scan too (as does `XY_SIMD=0`); both give the same rows.
+fn m4_bucket(y: &[f64], start: usize, end: usize, out: &mut Vec<u32>) {
+    if end - start >= M4_SIMD_MIN_ROWS {
+        let bucket = &y[start..end];
+        let Some(first) = bucket.iter().position(|v| v.is_finite()) else {
+            return;
+        };
+        let last = bucket
+            .iter()
+            .rposition(|v| v.is_finite())
+            .expect("first is finite");
+        let ys = &bucket[first..=last];
+        if let Some((min_i, max_i)) = crate::simd::try_argminmax(ys) {
+            if ys[min_i].is_finite() && ys[max_i].is_finite() {
+                let row = |i: usize| (start + first + i) as u32;
+                m4_flush(row(0), row(min_i), row(max_i), row(ys.len() - 1), out);
+                return;
+            }
+        }
+    }
+    let mut rows = (start..end).filter(|&i| y[i].is_finite());
+    let Some(first) = rows.next() else {
+        return;
+    };
+    let (mut min_i, mut max_i, mut last) = (first, first, first);
+    for i in rows {
+        if y[i] < y[min_i] {
+            min_i = i;
+        }
+        if y[i] > y[max_i] {
+            max_i = i;
+        }
+        last = i;
+    }
+    m4_flush(first as u32, min_i as u32, max_i as u32, last as u32, out);
+}
+
+fn m4_partition_bucket_end<F>(values: &[f64], bucket: usize, bucket_of: &F) -> usize
 where
     F: Fn(f64) -> usize,
 {
@@ -2621,7 +2757,7 @@ fn m4_indices_impl(
         let mut j = (start + t * chunk).min(end);
         if j < end {
             let bucket = bucket_of(x[j - 1]);
-            j += m4_bucket_end(&x[j..end], bucket, &bucket_of);
+            j += m4_partition_bucket_end(&x[j..end], bucket, &bucket_of);
         }
         let prev = *bounds.last().expect("bounds never empty");
         if j > prev && j < end {
@@ -6826,7 +6962,7 @@ mod tests {
             value as usize
         };
 
-        let offset = m4_bucket_end(&x[10_000..], 0, &bucket_of);
+        let offset = m4_partition_bucket_end(&x[10_000..], 0, &bucket_of);
 
         assert_eq!(offset, 184_000);
         assert!(
@@ -7334,6 +7470,40 @@ mod fuzz {
         }
     }
 
+    #[test]
+    fn fuzz_m4_bucket_walk_matches_per_row_loop() {
+        let mut rng = Rng(0x5EED_00A4);
+        for it in 0..2000 {
+            let n = 1 + (rng.next() % 3000) as usize;
+            // Ascending x: uniform, duplicate-heavy, or irregular.
+            let mut x: Vec<f64> = match it % 3 {
+                0 => (0..n).map(|i| i as f64).collect(),
+                1 => (0..n).map(|i| (i / (1 + it % 7)) as f64).collect(),
+                _ => (0..n).map(|_| rng.f01() * n as f64).collect(),
+            };
+            x.sort_by(f64::total_cmp);
+            let mut y = rng.hostile_vec(n, -1.0, 1.0);
+            for v in y.iter_mut() {
+                // Every other case drops ±∞ so long buckets reach argminmax;
+                // repeated values exercise first-row tie-breaking.
+                if it % 2 == 0 && v.is_infinite() {
+                    *v = 0.25;
+                }
+                if rng.next().is_multiple_of(8) {
+                    *v = 0.5;
+                }
+            }
+            let buckets = 1 + (rng.next() % 300) as usize;
+            let (x0, x1) = (-1.0, n as f64 + 1.0);
+            let inv = buckets as f64 / (x1 - x0);
+            assert_eq!(
+                m4_range_buckets(&x, &y, x0, inv, buckets, 0, n),
+                m4_range_rows(&x, &y, x0, inv, buckets, 0, n),
+                "it={it} n={n} buckets={buckets}"
+            );
+        }
+    }
+
     #[test]
     fn fuzz_m4_indices_window_and_bounds() {
         let mut rng = Rng(0x5EED_0004);