diff --git a/Cargo.lock b/Cargo.lock index 8afdf0f5..a150db50 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -8,6 +8,21 @@ version = "2.0.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" +[[package]] +name = "argminmax" +version = "0.7.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "088a9c28ca36500c581f9f54b63460492d2b4f71e1e61d176e8fb79565109b01" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" + [[package]] name = "bitflags" version = "2.13.0" @@ -58,6 +73,15 @@ dependencies = [ "simd-adler32", ] +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + [[package]] name = "png" version = "0.18.1" @@ -81,5 +105,6 @@ checksum = "703d5c7ef118737c72f1af64ad2f6f8c5e1921f818cdcb97b8fe6fc69bf66214" name = "xy-core" version = "0.1.0" dependencies = [ + "argminmax", "png", ] diff --git a/Cargo.toml b/Cargo.toml index 1a2946b4..372f6ee6 100644 --- a/Cargo.toml +++ b/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 diff --git a/news/535.performance.md b/news/535.performance.md new file mode 100644 index 00000000..a9cd0a7a --- /dev/null +++ b/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. diff --git a/spec/design-dossier.md b/spec/design-dossier.md index ac0fcc82..adc01306 100644 --- a/spec/design-dossier.md +++ b/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-` wheel covers every supported Python version on that platform without PyO3 or `abi3`. Wheel matrix in CI, release-blocking: manylinux diff --git a/spec/design/rust-engine.md b/spec/design/rust-engine.md index 1bb29a42..a9cc974a 100644 --- a/spec/design/rust-engine.md +++ b/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 diff --git a/src/kernels.rs b/src/kernels.rs index e56cff3d..f0417969 100644 --- a/src/kernels.rs +++ b/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 { + 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 { let mut out: Vec = Vec::with_capacity((n_buckets * 4).min((hi - lo) * 4)); @@ -2583,7 +2614,112 @@ fn m4_range( out } -fn m4_bucket_end(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 { + let bucket_of = |v: f64| (((v - x0) * inv_bucket_w) as usize).min(n_buckets - 1); + let mut out: Vec = 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) { + 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(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 = 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); diff --git a/src/simd.rs b/src/simd.rs index 96c7bb0a..2d195229 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -9,9 +9,13 @@ //! Only the kernels where the restructure measurably wins live here (this //! box, 10M points): min_max +89%, bin_2d +~50%, local-density gather +~42%, //! range scan +20–50% at ≤1M (wash at 10M where it's bandwidth-bound). -//! M4 and histogram were tried and REVERTED — measured slower (the bucket -//! state machine / scatter increment dominates and the extra block pass costs -//! more than the vector math saves); see the notes on their scalar kernels. +//! Histogram was tried and REVERTED — measured slower (the scatter increment +//! dominates and the extra block pass costs more than the vector math saves); +//! see the note on its scalar kernel. The same two-phase restructure lost on +//! M4 too, so M4 windows that average 8+ rows per bucket scan whole buckets +//! instead: [`try_argminmax`] hands the y of each bucket of 32+ rows to the +//! `argminmax` crate, which picks its own SIMD width at run time +//! (SSE4.1/AVX2/AVX-512, NEON): M4 at 1M points runs 3–9× faster. //! //! Rules (rust-engine doc §3.4): //! - Only kernels whose outputs are **order-independent** are accelerated: @@ -27,7 +31,8 @@ //! before/after benchmark and a debugging escape hatch. //! //! Non-x86_64 (e.g. aarch64) needs nothing here: NEON is part of the aarch64 -//! baseline, so the scalar kernels already autovectorize at full width there. +//! baseline, so the scalar kernels already autovectorize at full width there +//! (argminmax uses NEON on its own). #![allow( dead_code, @@ -44,14 +49,16 @@ const BLOCK: usize = 1024; /// bucket indices are always well below this (grids are screen-sized). const SKIP: u32 = u32::MAX; -#[cfg(target_arch = "x86_64")] -pub(crate) fn use_avx2() -> bool { +/// The `XY_SIMD=0` kill switch, read once. +fn simd_enabled() -> bool { use std::sync::OnceLock; static ON: OnceLock = OnceLock::new(); - *ON.get_or_init(|| { - std::env::var_os("XY_SIMD").is_none_or(|v| v != "0") - && std::arch::is_x86_feature_detected!("avx2") - }) + *ON.get_or_init(|| std::env::var_os("XY_SIMD").is_none_or(|v| v != "0")) +} + +#[cfg(target_arch = "x86_64")] +pub(crate) fn use_avx2() -> bool { + simd_enabled() && std::arch::is_x86_feature_detected!("avx2") } #[cfg(not(target_arch = "x86_64"))] @@ -276,6 +283,13 @@ unsafe fn density_gather_avx2( // every `unsafe` call site inside this module. // --------------------------------------------------------------------------- +/// Rows of the first minimum and the first maximum of a non-empty `data`, +/// NaN skipped (±∞ count as values), from the `argminmax` crate. It detects +/// the CPU itself, so this needs only the kill switch. +pub(crate) fn try_argminmax(data: &[f64]) -> Option<(usize, usize)> { + simd_enabled().then(|| argminmax::ArgMinMax::argminmax(&data)) +} + pub(crate) fn try_min_max(data: &[f64]) -> Option<(f64, f64)> { #[cfg(target_arch = "x86_64")] if use_avx2() {