From 3f84990c7fcb5951b1561d8c10df8f19d52eeea9 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 8 Oct 2026 18:21:55 +0000 Subject: [PATCH] D-MHB-1: Mexican-hat response without a raster (probe + board) Measures whether popcount stacking, a quantised DoG and an exact early exit can answer a centre-surround query without per-point geometry or a materialised raster. Probe only; no production primitive. - Bit-sliced weight planes (F) are exact (asserted equal to direct geometry on every centre) at a constant ~600 ns per query. - Direct window geometry (D) wins below a density of about 0.1. - Equal-q ring buckets (E) are fast but wrong at every K tried. - The shipped mask-risc path (B) materialises a weight lane per query and is 26-150x slower. - The exact-bound early exit (H) saves about 20 % at the median threshold and equals the full decision on every query. - Rolling-Floor control and statistical early exit are parked: the deciding quantity is an exact popcount. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01EFw2WdKr1oxvaKCJC2ua2R --- .claude/board/STATUS_BOARD.md | 8 + .../2026-10-08-mexhat-bucket-cascade-probe.md | 221 ++++ .claude/board/entries/README.md | 3 +- .../examples/mexhat_bucket_probe.rs | 1128 +++++++++++++++++ 4 files changed, 1359 insertions(+), 1 deletion(-) create mode 100644 .claude/board/entries/2026-10-08-mexhat-bucket-cascade-probe.md create mode 100644 crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs diff --git a/.claude/board/STATUS_BOARD.md b/.claude/board/STATUS_BOARD.md index a103a47a7..290bff029 100644 --- a/.claude/board/STATUS_BOARD.md +++ b/.claude/board/STATUS_BOARD.md @@ -1,3 +1,11 @@ +## D-MHB — Centre-surround response without a raster (2026-10-08) + +Board entry: `entries/2026-10-08-mexhat-bucket-cascade-probe.md`. Popcount stacking, quantised DoG and exact early exit against direct window geometry; no wavefield, no raster. + +| D-id | scope | status | gate / falsifier | +|---|---|---|---| +| **D-MHB-1** | Mexican-hat response: direct geometry (D), ring buckets (E), bit-sliced weight planes (F), per-row choice (G), mask-risc lane (B), exact-bound early exit (H) | Shipped (probe; R-MHB-1 proposed) | D = F = G = B = full scan asserted on every centre; H equals the full decision on every query; one zero crossing, one minimum, symmetry; dropping the unvisited-row bound decides 52 / 512 wrongly | + ## D-RPF — Mask × fold over projections of resident bytes (2026-10-08) Plan: `.claude/plans/2026-10-08-resident-projection-fold-mask-v1.md`. Mask = admissibility, fold = one terminal, CE64 register = one instruction; no materialised hop between them. diff --git a/.claude/board/entries/2026-10-08-mexhat-bucket-cascade-probe.md b/.claude/board/entries/2026-10-08-mexhat-bucket-cascade-probe.md new file mode 100644 index 000000000..731981fd1 --- /dev/null +++ b/.claude/board/entries/2026-10-08-mexhat-bucket-cascade-probe.md @@ -0,0 +1,221 @@ +# 2026-10-08 — Mexican hat × popcount stacking × early exit, without a raster (D-MHB-1) + +**Status:** MEASURED. Probe only, no production primitive. +- **Probe:** `crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs`. +- **Revisions:** lance-graph `origin/main` d127f7dc; ndarray `origin/master` 30ce119. + +```text +CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example mexhat_bucket_probe +``` + +**The question.** Can popcount stacking, a quantised Mexican hat and an exact +early exit answer a centre-surround query without per-point geometry and +without a materialised raster? If so, is a Rolling-Floor controller worth +adding on top? + +**The query:** + +```text +S(c) = Σ_{p ∈ P, |p − c|² ≤ R²} w(|p − c|²) +``` + +- `P` is a presence bitmap over a row-major `W × H` grid. +- `w` is a normalised DoG, quantised to integers with `w(0) = 2^12`. +- The kernel is σc 3, σs 6, R 18: a 37 × 37 window holding 1,009 disk cells. + +## Inventory + +### SHIPPED (VERIFIED-IN-CODE) + +**ndarray:** +- **`RollingFloor`** (`hpc/rolling_floor.rs`): + - exact `(n, Σx, Σx²)` moments, deterministic reservoir; + - a quarter-σ lattice, thresholds derived on demand; + - drift read at the 1σ/3σ cuts. + - Lower is stronger (`shade`, l. 755). + - Not re-exported from `ndarray::simd`. +- **`hamming_distance_within`** (`bitwise.rs`): an exact monotone early exit per 256-byte block. Re-exported from `ndarray::simd`. +- **`Cascade`** (`hpc/cascade.rs`): + - Stroke 2 has an exact budget early exit (l. 346). + - **Stroke 1 is a statistical prune with no exact fallback.** A prefix whose estimate exceeds `threshold + 3σ` is dropped (l. 329–337). An unresolved candidate disappears instead of falling through to an exact check. + +**lance-graph:** +- **mask-risc** fused `Count` / `Any` over resident planes. +- **`MaskedSumI32`** with `execute_extent`. + +### PROPOSED, or present in name only + +- **`styles/lsi.rs` `mexican_hat`** is a five-step percentile function in f32 (1 / 0.5 / 0 / −0.5 / −1). The percentiles come from a normal fit. It is not a normalised DoG. +- **`pillar/mexican_hat.rs` (Pillar-15) is DEFERRED.** `prove_pillar_15` runs nothing and returns `passed: true`. No DoG kernel exists in ndarray. + +### MISSING + +- A DoG kernel in ndarray. +- A 2-D window or stencil operand in mask-risc: planes are 1-D row spaces. + +### Observation 6, reproduced + +`examples/rolling_floor_probe` (ndarray) still prints the same result: + +| phase | FIXED reject | SHIPPED reject | EWMA reject | +|---|---|---|---| +| 2 | 97.55 % | 97.55 % | 0.12 % | +| 4 | 64.10 % | 64.10 % | 0.07 % | + +FIXED and SHIPPED are identical because `Cascade`'s `ShiftAlert` never fires on a per-sample feed. That probe exercises **`Cascade`**, not `RollingFloor`. `RollingFloor`'s checkpoint drift rule is untested by it. + +## Semantics, kept apart + +- **Continuous** (arm A, f64): the reference for the quantisation error. +- **Quantised** (integer `w(q)` per lattice offset): arms Ascan, D, F, G and B must agree, and the probe **asserts exact equality** on every checked centre: + - 1,024 centres × 4 densities × 2 grid sizes; + - B on every timed query. + + Its distance from the continuous reference, as a mean of |S|: + + | ρ | quantised vs continuous | + |---|---| + | 0.01 | 0.07 % | + | 0.1 | 0.18 % | + | 0.5 | 0.5 % | + | 0.9 | 1.3 % | + +- **Bucketed** (arm E: K equal-q rings, one mean weight each): a further approximation, measured against the quantised answer, never asserted equal. + +## Kernel falsifiers + +All six `(σc, κ)` pairs pass, for κ ∈ {1.5, 2, 3}: +- a positive centre; +- a negative surround; +- **exactly one sign change**, at the analytic radius² `2σc²σs² ln(σs²/σc²)/(σs² − σc²)` within one q-step (e.g. 59.15 lands between lattice q 59 and 60); +- **one annular minimum**: no false extremum after quantisation; +- the 8 square symmetries of the template; +- rings that partition the disk; +- no overflow of the worst-case bounds. + +## Results + +**Setup:** median of 10 runs. +- Host: Intel Xeon @ 2.8 GHz, 4 vCPU (AVX-512 on the host). +- Build: `x86-64-v3` (AVX2, POPCNT), release, `debug = 0`, rustc 1.98.1. +- Single runs scatter up to 2×, which is why medians are reported. + +**ns per query:** + +| ρ (1M grid) | A f64 | Ascan | **D** window + LUT | E K=8 | E K=32 | **F** bit-sliced | G per-row | B mask-risc | +|---|---|---|---|---|---|---|---|---| +| 0.01 | 529 | 125 k | **246** | 489 | 1378 | 624 | 267 | 61.7 k | +| 0.1 | 2812 | 499 k | 619 | 464 | 1376 | **592** | 620 | 62.4 k | +| 0.5 | 11.8 k | 2.0 M | 1166 | 489 | 1413 | **624** | 894 | 80.4 k | +| 0.9 | 21.0 k | 3.5 M | 1516 | 473 | 1374 | **604** | 679 | 88.9 k | +| wavefront ring + 1 % noise | 560 | 135 k | **252** | 477 | 1334 | 591 | 270 | 64.0 k | + +The 64K grid gives the same orderings (D 230 / 620 / 1070 / 1516; F ≈ 600). + +| arm | work per query | memory | +|---|---|---| +| D | ~1.4 ns per present disk cell + ~210 ns base | | +| F | constant: 418 popcounts | | +| E, F, D | read about 57 `u64` words (the row-major window) | 0 allocations per query | +| B | | writes a 37.9 KB (64K) / 151.6 KB (1M) weight lane per query, over a 4 MB resident lane | + +## Exact-bound early exit (arm H) + +**The decision:** `S ≥ T`, using F's rows, centre rows first. + +**The bounds:** each row's suffix bound is `[Σ negative, Σ positive]` of the rows not yet visited. The early decision is: + +| condition | decision | +|---|---| +| `S + U < T` | reject | +| `S + L ≥ T` | accept | +| otherwise | visit the next row | + +**Correctness:** the decision is asserted equal to the full answer for every query: 4,096 queries × 4 thresholds × 5 fixtures. + +**Rows visited** (of 37; mean, with p95 in brackets): + +| T quantile | ρ 0.01 | ρ 0.5 | ρ 0.9 | +|---|---|---|---| +| 10 % | 30.3 | 21.6 | 15.8 | +| 50 % | 24.3 | 24.5 | 24.0 | +| 90 % | 14.2 | 21.7 | 26.1 | +| 99 % | 9.9 (p95 14) | 17.8 | 24.4 | + +**Timing at the median T:** H 500–534 ns against F 628–651 ns, about −20 %. On sparse input, D (~250 ns) beats H. + +## Counterexamples + +1. **The running sum peaks above T, then ends below it.** A positive core plus a late negative surround: peak 143,154 above T = 100,777, final S = 58,400. H stays undecided while the sum is above T and rejects after 18 of 37 rows. +2. **Cancellation.** A ring placed exactly at the zero crossing gives S = 116. Deciding `S ≥ 1` needs 35 of 37 rows. +3. **Dense input with T at the median:** 24.2 of 37 rows. The bound cannot help here. +4. **Disable run, in place.** An upper bound that drops the unvisited rows decides 52 of 512 queries wrongly, so the bound is load-bearing. + A second disable run, made after the commit: shifting F's negative planes one bit too far fails `F differs from D` at the first checked centre. The file was restored afterwards. +5. **Geometry buckets cannot see phase.** Two coherent sources, λ = 8, detector cells bucketed by `(r1, r2)`: + - buckets of width λ/2: 738 of 1,003 occupied buckets hold both a dark (I < 0.1) and a bright (I > 0.9) cell; + - buckets of width λ/8: 0 of 8,032 do. + + A geometric prune coarser than about λ/8 drops dark fringes with the noise. +6. **Locality.** A 37 × 37 window touches 56.7 row-major words and 30.2 Morton words. Morton touches fewer words, but its bits are not row-shaped. Arms D/E/F extract row windows, so Morton needs different (8 × 8 tile) templates. Word counts only, not timed. + +## Findings + +1. **The exact answer is already cheap; the bucketed one is wrong.** Equal-q rings (E): + - K = 8 is fast (~470 ns, constant) but wrong by 950–5,600 on |S| in the thousands; + - K = 32 is still wrong by 255–1,830 and slower than exact F. + + A popcount ring fold buys nothing an exact path does not. +2. **Bit-slicing makes popcount stacking exact.** For a static integer weight template: + + ```text + Σ_p w(p) = Σ_b 2^b · (popcount(P ∧ Pos_b) − popcount(P ∧ Neg_b)) + ``` + + - It costs 2 × 13 planes per row, independent of density, at ~600 ns. + - It is identical to direct geometry: asserted on every centre. + - It needs no value lane, no square root and no `exp`. +3. **Crossover at ρ ≈ 0.1** (≈ 100 present cells in a 1,009-cell disk): + - below it, direct per-point geometry D wins: 246 vs 624 ns at ρ 0.01, and on the wavefront fixture; + - above it, F wins: 624 vs 1,166 ns at ρ 0.5, 604 vs 1,516 at 0.9. +4. **Per-row adaptivity (G) does not reach min(D, F).** Choosing per row with one popcount is exact, but it loses up to 43 % at mid density (894 vs 624 ns at ρ 0.5). Choose per **query** instead, from the window's own popcount (37 popcounts ≈ 50 ns). That is an exact count, not a statistic. +5. **The shipped mask-risc path (B) is 26–40× slower than F on the 64K grid and 100–150× on the 1M grid** (the band is full width) and materialises a weight lane per query. The cause is that mask-risc has no 2-D window operand. +6. **Exact early exit is a modest, safe win.** About −20 % on a threshold decision over F, up to ~4× fewer rows at extreme quantiles on sparse input. It is never wrong, by construction and as asserted. + +## Recommendations + +| item | verdict | why | +|---|---|---| +| Bit-sliced weight fold (F) | **ADOPT as the dense path** (a separate PR, routed through `ndarray::simd` popcount) | exact, constant cost, no lane, no sqrt/exp | +| Direct window geometry (D) | **ADOPT as the sparse path** | fastest below ρ ≈ 0.1 | +| Per-query D/F choice from the window popcount | **PROBE** | exact count, ~50 ns; should reach min(D, F), not measured yet | +| Exact-bound early exit (H) | **PROBE** | −20 % at median T; needs a D-path variant for sparse input | +| Ring-bucket Mexican hat (E) | **REJECT** | wrong at every K tried and not faster than exact F | +| Rolling-Floor-controlled cascade | **PARK** | the deciding quantity (window density) is an exact popcount; a statistical controller has nothing to decide here | +| Statistical early exit + exact fallback | **PARK** | the exact bound already decides without error | +| mask-risc `MaskedSumI32` over a per-query weight lane (B) | **REJECT for this query** | 26–150× slower, materialises a lane | +| `lsi.rs` `mexican_hat` as a DoG | **REJECT the equivalence** | a percentile step function, not a DoG | +| Pillar-15 activation | **PARK** (separate from any planner change) | needs a real DoG kernel in ndarray first; the deferred stub reports `passed: true` | + +## Proposed rewrite (not implemented) + +**R-MHB-1, bit-sliced weighted fold.** +- **When:** a masked sum whose weights come from a static, narrow integer template. +- **Rewrite:** `MaskedSumI32` over that template → `Σ_b 2^b · (Count(P ∧ Pos_b) − Count(P ∧ Neg_b))` over resident planes. +- This turns a value-lane fold into fused `Count`s (D-RPF-9 finding 1: no slot, no bitmap). +- **Exact** for the quantised semantics. +- **Applies only** when the template planes are resident and aligned with `P`. For a moving window that needs the missing 2-D window operand. + +## Boundaries + +- No wavefield or raster is built. The only per-query state is a running sum and the window words. +- An interference intensity `|E1 + E2|²` stays a phase computation. A centre-surround score is not a fringe classifier (counterexample 5). +- Every popcount in the probe is `u64::count_ones`. A production path goes through `ndarray::simd`. + +## OPEN + +- Per-query D/F selection: not yet measured. +- An early exit on the D path for sparse input: not yet measured. +- Morton-shaped (8 × 8) templates: not timed. +- `RollingFloor`'s own drift rule under drifting density streams: untested. +- `Cascade` Stroke 1 silently drops unresolved candidates. It is an ndarray behaviour, reported here, not changed. +- One host class only. The ρ ≈ 0.1 crossover is a pin for this host. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index 3ffc42e98..b54dd5b7f 100644 --- a/.claude/board/entries/README.md +++ b/.claude/board/entries/README.md @@ -25,7 +25,7 @@ index row, (3) no duplicate entry id. Checks 1 and 2 are deliberately opposite directions; the stranding this convention prevents shows up in exactly one of them, never both. -273 entries, 2026-08-06 .. 2026-10-08. +274 entries, 2026-08-06 .. 2026-10-08. | date | entry id | finding | file | |---|---|---|---| @@ -35,6 +35,7 @@ exactly one of them, never both. | 2026-10-08 | `rbac-hotplug-socket` | | [2026-10-08-rbac-hotplug-socket.md](2026-10-08-rbac-hotplug-socket.md) | | 2026-10-08 | `moore-nars16-isa-visible-representation` | | [2026-10-08-moore-nars16-isa-visible-representation.md](2026-10-08-moore-nars16-isa-visible-representation.md) | | 2026-10-08 | `D-MOORE-NARS-0` | | [2026-10-08-moore-nars-0-recipe-learning-gomoku.md](2026-10-08-moore-nars-0-recipe-learning-gomoku.md) | +| 2026-10-08 | `D-MHB-1` | | [2026-10-08-mexhat-bucket-cascade-probe.md](2026-10-08-mexhat-bucket-cascade-probe.md) | | 2026-10-08 | `hhtl-nars-moore-value-tenants` | | [2026-10-08-hhtl-nars-moore-value-tenants.md](2026-10-08-hhtl-nars-moore-value-tenants.md) | | 2026-10-08 | `D-RPF-9` | | [2026-10-08-fold-join-deforestation-probe.md](2026-10-08-fold-join-deforestation-probe.md) | | 2026-10-08 | `D-RPF-0` | | [2026-10-08-d-rpf-0-resident-ce64-predicate.md](2026-10-08-d-rpf-0-resident-ce64-predicate.md) | diff --git a/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs b/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs new file mode 100644 index 000000000..7391b9788 --- /dev/null +++ b/crates/lance-graph-mask-risc/examples/mexhat_bucket_probe.rs @@ -0,0 +1,1128 @@ +//! D-MHB-1 — Mexican-hat response without a raster: does popcount stacking +//! (rings or bit-sliced weight planes) beat direct per-point geometry, and how +//! much does an exact-bound early exit save on a threshold decision? +//! +//! Nothing here is a production primitive. The question is answered for ONE +//! query shape: the centre-surround response at a query cell `c` of a +//! presence bitmap `P` over a row-major `W × H` grid, +//! +//! ```text +//! S(c) = Σ_{p ∈ P, |p − c|² ≤ R²} w(|p − c|²) +//! ``` +//! +//! where `w` is a normalised Difference-of-Gaussians quantised to integers +//! (`w(q) = round(DoG(√q) / DoG(0) · 2^SHIFT)`, evaluated on `q = r²`). +//! +//! Three semantics, never mixed: +//! +//! - **continuous**: `f64` DoG, the reference the quantisation is measured +//! against (arm A); +//! - **quantised**: integer `w(q)` per lattice offset. Arms Ascan, D, F, G and B +//! must agree on it EXACTLY, and the probe asserts that they do; +//! - **bucketed**: `K` radius² rings, one representative weight each (arm E). +//! A further approximation, measured against the quantised answer. +//! +//! | arm | execution | +//! |---|---| +//! | A | `f64` DoG over present points in the window (continuous reference) | +//! | Ascan | every set bit of the WHOLE grid, `q`, LUT (full candidate scan) | +//! | D | set bits in the window only, `q = dx² + dy²`, LUT | +//! | E | `K` ring templates: `Σ_k w_k · popcount(win ∧ ring_k)` | +//! | F | bit-sliced weight planes: `Σ_b 2^b · (popcnt(win ∧ Pos_b) − popcnt(win ∧ Neg_b))` | +//! | G | per row, the cheaper of D and F (both exact): one popcount decides | +//! | B | shipped mask-risc: per-query weight lane over the row band, `MaskedSumI32` over `P` | +//! | H | early exit for `S ≥ T` over F's rows, with exact suffix bounds `[L, U]` | +//! +//! Timings are printed, never asserted. Every equality is asserted. +//! +//! ```text +//! CARGO_PROFILE_RELEASE_DEBUG=0 cargo run --release -p lance-graph-mask-risc --example mexhat_bucket_probe +//! ``` + +// The template builders index window rows and columns as 2-D coordinates +// (`dy = yi − R`, `dx = xi − R`) into several tables at once; an iterator per +// table would hide the geometry the loops exist to express. +#![allow(clippy::needless_range_loop)] + +use std::alloc::{GlobalAlloc, Layout, System}; +use std::sync::atomic::{AtomicUsize, Ordering}; +use std::time::Instant; + +use lance_graph_mask_risc::exec::{execute_extent, Scratch}; +use lance_graph_mask_risc::{Foreign, LaneRef, Operand, Out, Planes, Program, Terminal, Value}; + +// ───────────────────────── allocation counter ───────────────────────── + +struct Counting; +static ALLOCS: AtomicUsize = AtomicUsize::new(0); + +// SAFETY: a pure pass-through to `System`; the counter is the only addition. +unsafe impl GlobalAlloc for Counting { + unsafe fn alloc(&self, layout: Layout) -> *mut u8 { + ALLOCS.fetch_add(1, Ordering::Relaxed); + // SAFETY: same layout, same contract as the caller's. + unsafe { System.alloc(layout) } + } + unsafe fn dealloc(&self, ptr: *mut u8, layout: Layout) { + // SAFETY: `ptr` came from `alloc` above with this `layout`. + unsafe { System.dealloc(ptr, layout) } + } +} + +#[global_allocator] +static GLOBAL: Counting = Counting; + +fn allocs() -> usize { + ALLOCS.load(Ordering::Relaxed) +} + +// ───────────────────────── rng ───────────────────────── + +struct Rng(u64); +impl Rng { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = self.0; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + fn unit(&mut self) -> f64 { + (self.next() >> 11) as f64 / (1u64 << 53) as f64 + } + fn below(&mut self, n: usize) -> usize { + (self.next() % n as u64) as usize + } +} + +// ───────────────────────── the kernel ───────────────────────── + +/// Weight scale: `w(0) = 2^SHIFT`. +const SHIFT: u32 = 12; +/// Bit planes needed for `|w| ≤ 2^SHIFT`. +const PLANES: usize = SHIFT as usize + 1; + +/// Normalised DoG at squared radius `q` (no square root needed). +fn dog(q: f64, sc: f64, ss: f64) -> f64 { + let g = |s: f64| (-q / (2.0 * s * s)).exp() / (2.0 * std::f64::consts::PI * s * s); + g(sc) - g(ss) +} + +/// Analytic zero crossing of the normalised DoG, as a squared radius: +/// `r0² = 2 σc² σs² ln(σs²/σc²) / (σs² − σc²)`. +fn zero_crossing_q(sc: f64, ss: f64) -> f64 { + let (a, b) = (sc * sc, ss * ss); + 2.0 * a * b * (b / a).ln() / (b - a) +} + +struct Kernel { + sc: f64, + ss: f64, + /// Cutoff radius; the window is `2R + 1` cells wide and must fit a `u64`. + r: usize, + r2: usize, + /// Quantised weight per squared radius `0..=R²`. + lut: Vec, + /// Per window row `dy + R`: the disk cells (`q ≤ R²`) as a `u64` over `dx + R`. + disk: Vec, + /// Per window row: bit-sliced positive and negative magnitude planes. + pos: Vec<[u64; PLANES]>, + neg: Vec<[u64; PLANES]>, + /// Per window row: the largest positive / most negative contribution the + /// row could still make (all its positive / negative cells present). + row_up: Vec, + row_lo: Vec, + /// Per window row: how many of its `2 · PLANES` planes are non-zero, + /// i.e. what arm F pays for the row in popcounts. + planes_used: Vec, +} + +impl Kernel { + fn new(sc: f64, kappa: f64) -> Self { + let ss = sc * kappa; + let r = (3.0 * ss).ceil() as usize; + assert!(2 * r < 64, "window must fit one u64 (R = {r})"); + let r2 = r * r; + let peak = dog(0.0, sc, ss); + let lut: Vec = (0..=r2) + .map(|q| (dog(q as f64, sc, ss) / peak * f64::from(1u32 << SHIFT)).round() as i32) + .collect(); + let wd = 2 * r + 1; + let (mut disk, mut pos, mut neg, mut row_up, mut row_lo) = ( + vec![0u64; wd], + vec![[0u64; PLANES]; wd], + vec![[0u64; PLANES]; wd], + vec![0i64; wd], + vec![0i64; wd], + ); + for yi in 0..wd { + let dy = yi as i64 - r as i64; + for xi in 0..wd { + let dx = xi as i64 - r as i64; + let q = (dx * dx + dy * dy) as usize; + if q > r2 { + continue; + } + disk[yi] |= 1 << xi; + let w = lut[q]; + let (planes, mag) = if w >= 0 { + (&mut pos[yi], w) + } else { + (&mut neg[yi], -w) + }; + for (b, p) in planes.iter_mut().enumerate() { + if (mag >> b) & 1 == 1 { + *p |= 1 << xi; + } + } + if w > 0 { + row_up[yi] += i64::from(w); + } else { + row_lo[yi] += i64::from(w); + } + } + } + let planes_used = (0..wd) + .map(|yi| { + (0..PLANES) + .map(|b| u32::from(pos[yi][b] != 0) + u32::from(neg[yi][b] != 0)) + .sum() + }) + .collect(); + Kernel { + sc, + ss, + r, + r2, + lut, + disk, + pos, + neg, + row_up, + row_lo, + planes_used, + } + } + + fn wd(&self) -> usize { + 2 * self.r + 1 + } +} + +/// `K` radius² rings over `0..=R²` with one representative weight each: the +/// mean of the quantised weights of the lattice cells the ring contains. +struct Rings { + /// `ring[k][row]`: ring `k`'s cells in window row `row`. + ring: Vec>, + w: Vec, +} + +fn rings(k: &Kernel, n: usize) -> Rings { + // Equal-width rings in q. Edge rule: cell with q goes to ring + // `min(q * n / (R² + 1), n − 1)` — one ring per q, deterministically. + let ring_of = |q: usize| (q * n / (k.r2 + 1)).min(n - 1); + let wd = k.wd(); + let mut ring = vec![vec![0u64; wd]; n]; + let (mut sum, mut cnt) = (vec![0i64; n], vec![0i64; n]); + for yi in 0..wd { + let dy = yi as i64 - k.r as i64; + for xi in 0..wd { + let dx = xi as i64 - k.r as i64; + let q = (dx * dx + dy * dy) as usize; + if q > k.r2 { + continue; + } + let b = ring_of(q); + ring[b][yi] |= 1 << xi; + sum[b] += i64::from(k.lut[q]); + cnt[b] += 1; + } + } + let w = sum + .iter() + .zip(&cnt) + .map(|(s, c)| { + if *c == 0 { + 0 + } else { + (*s as f64 / *c as f64).round() as i64 + } + }) + .collect(); + Rings { ring, w } +} + +// ───────────────────────── the grid ───────────────────────── + +struct Grid { + w: usize, + h: usize, + /// Row-major presence, `w / 64` words per row. + p: Vec, +} + +impl Grid { + fn words_per_row(&self) -> usize { + self.w / 64 + } + fn set(&mut self, x: usize, y: usize) { + let i = y * self.w + x; + self.p[i / 64] |= 1 << (i % 64); + } + /// `len ≤ 63` bits of row `y` starting at column `x0`, bit `j` = column `x0 + j`. + fn window(&self, x0: usize, y: usize, len: usize) -> u64 { + let base = y * self.words_per_row(); + let (wi, sh) = (x0 / 64, x0 % 64); + let lo = self.p[base + wi] >> sh; + let v = if sh != 0 && sh + len > 64 { + lo | (self.p[base + wi + 1] << (64 - sh)) + } else { + lo + }; + v & ((1u64 << len) - 1) + } +} + +fn random_grid(w: usize, h: usize, rho: f64, seed: u64) -> Grid { + let mut g = Grid { + w, + h, + p: vec![0; w * h / 64], + }; + let mut r = Rng(seed); + for y in 0..h { + for x in 0..w { + if r.unit() < rho { + g.set(x, y); + } + } + } + g +} + +/// A wavefront: one ring of radius `r0` and width ~1 cell around the grid +/// centre, plus `noise` density elsewhere. +fn wavefront_grid(w: usize, h: usize, r0: f64, noise: f64, seed: u64) -> Grid { + let mut g = random_grid(w, h, noise, seed); + let (cx, cy) = (w as f64 / 2.0, h as f64 / 2.0); + for y in 0..h { + for x in 0..w { + let d = ((x as f64 - cx).powi(2) + (y as f64 - cy).powi(2)).sqrt(); + if (d - r0).abs() < 0.5 { + g.set(x, y); + } + } + } + g +} + +// ───────────────────────── arms ───────────────────────── + +#[derive(Default, Clone, Copy)] +struct Work { + /// `q` computations (one per visited present point). + geometry: u64, + /// `exp` evaluations. + exps: u64, + popcounts: u64, + rows: u64, +} + +/// Arm A, continuous: `f64` DoG over present points of the window. +fn arm_a(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> f64 { + let peak = dog(0.0, k.sc, k.ss); + let mut s = 0.0; + for yi in 0..k.wd() { + let y = cy + yi - k.r; + let dy = yi as i64 - k.r as i64; + let mut m = g.window(cx - k.r, y, k.wd()); + while m != 0 { + let xi = m.trailing_zeros() as i64; + m &= m - 1; + let dx = xi - k.r as i64; + let q = (dx * dx + dy * dy) as usize; + wk.geometry += 1; + if q <= k.r2 { + wk.exps += 2; + s += dog(q as f64, k.sc, k.ss) / peak * f64::from(1u32 << SHIFT); + } + } + } + s +} + +/// Arm Ascan: every set bit of the whole grid. +fn arm_scan(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> i64 { + let mut s = 0i64; + for (wi, word) in g.p.iter().enumerate() { + let mut m = *word; + while m != 0 { + let i = wi * 64 + m.trailing_zeros() as usize; + m &= m - 1; + let (x, y) = ((i % g.w) as i64, (i / g.w) as i64); + let (dx, dy) = (x - cx as i64, y - cy as i64); + wk.geometry += 1; + let q = dx * dx + dy * dy; + if q <= k.r2 as i64 { + s += i64::from(k.lut[q as usize]); + } + } + } + s +} + +/// Arm D: set bits of the window, `q`, LUT. +fn arm_d(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> i64 { + let mut s = 0i64; + for yi in 0..k.wd() { + let dy = yi as i64 - k.r as i64; + let mut m = g.window(cx - k.r, cy + yi - k.r, k.wd()) & k.disk[yi]; + while m != 0 { + let dx = m.trailing_zeros() as i64 - k.r as i64; + m &= m - 1; + wk.geometry += 1; + s += i64::from(k.lut[(dx * dx + dy * dy) as usize]); + } + } + s +} + +/// Arm E: ring popcounts times the ring's representative weight. +fn arm_e(g: &Grid, k: &Kernel, rg: &Rings, cx: usize, cy: usize, wk: &mut Work) -> i64 { + let mut s = 0i64; + for yi in 0..k.wd() { + let win = g.window(cx - k.r, cy + yi - k.r, k.wd()); + for (ring, w) in rg.ring.iter().zip(&rg.w) { + let r = ring[yi]; + if r != 0 { + wk.popcounts += 1; + s += w * i64::from((win & r).count_ones()); + } + } + } + s +} + +/// One window row of arm F: the row's exact quantised contribution. +#[inline] +fn f_row(k: &Kernel, win: u64, yi: usize, wk: &mut Work) -> i64 { + let mut s = 0i64; + for b in 0..PLANES { + let (p, n) = (k.pos[yi][b], k.neg[yi][b]); + if p != 0 { + wk.popcounts += 1; + s += i64::from((win & p).count_ones()) << b; + } + if n != 0 { + wk.popcounts += 1; + s -= i64::from((win & n).count_ones()) << b; + } + } + s +} + +/// Arm F: bit-sliced weight planes, exact for the quantised semantics. +fn arm_f(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> i64 { + (0..k.wd()) + .map(|yi| { + wk.rows += 1; + f_row(k, g.window(cx - k.r, cy + yi - k.r, k.wd()), yi, wk) + }) + .sum() +} + +/// Arm G: per row, the cheaper of the two exact executions. A row with fewer +/// present disk cells than non-zero planes takes D's per-point path, +/// otherwise F's popcounts. Both are exact, so the choice cannot change the +/// answer; the probe asserts that it does not. +fn arm_g(g: &Grid, k: &Kernel, cx: usize, cy: usize, wk: &mut Work) -> i64 { + let mut s = 0i64; + for yi in 0..k.wd() { + let win = g.window(cx - k.r, cy + yi - k.r, k.wd()); + let mut m = win & k.disk[yi]; + wk.popcounts += 1; + if m.count_ones() <= k.planes_used[yi] { + let dy = yi as i64 - k.r as i64; + while m != 0 { + let dx = m.trailing_zeros() as i64 - k.r as i64; + m &= m - 1; + wk.geometry += 1; + s += i64::from(k.lut[(dx * dx + dy * dy) as usize]); + } + } else { + wk.rows += 1; + s += f_row(k, win, yi, wk); + } + } + s +} + +/// Rows in the order arm H visits them: largest `|contribution|` bound first. +fn h_order(k: &Kernel) -> Vec { + let mut o: Vec = (0..k.wd()).collect(); + o.sort_by_key(|&yi| std::cmp::Reverse(k.row_up[yi] - k.row_lo[yi])); + o +} + +/// Suffix bounds over `order`: `up[i]` / `lo[i]` bound what rows `order[i..]` +/// can still add. +fn h_suffix(k: &Kernel, order: &[usize]) -> (Vec, Vec) { + let n = order.len(); + let (mut up, mut lo) = (vec![0i64; n + 1], vec![0i64; n + 1]); + for i in (0..n).rev() { + up[i] = up[i + 1] + k.row_up[order[i]]; + lo[i] = lo[i + 1] + k.row_lo[order[i]]; + } + (up, lo) +} + +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +enum Decision { + Accept, + Reject, +} + +/// Arm H: decide `S ≥ t` over F's rows, stopping as soon as the exact suffix +/// bounds settle it. Returns the decision and the rows visited. +#[allow(clippy::too_many_arguments)] +fn arm_h( + g: &Grid, + k: &Kernel, + order: &[usize], + up: &[i64], + lo: &[i64], + cx: usize, + cy: usize, + t: i64, + wk: &mut Work, +) -> (Decision, usize) { + let mut s = 0i64; + for i in 0..=order.len() { + // Rows order[i..] can add at most up[i] and at least lo[i]. + if s + up[i] < t { + return (Decision::Reject, i); + } + if s + lo[i] >= t { + return (Decision::Accept, i); + } + if i == order.len() { + break; + } + let yi = order[i]; + wk.rows += 1; + s += f_row(k, g.window(cx - k.r, cy + yi - k.r, k.wd()), yi, wk); + } + unreachable!("with no rows left the bounds are [0, 0], so one of the two tests fired"); +} + +// ───────────────────────── harness ───────────────────────── + +fn centres(g: &Grid, k: &Kernel, n: usize, seed: u64) -> Vec<(usize, usize)> { + let mut r = Rng(seed); + (0..n) + .map(|_| (k.r + r.below(g.w - 2 * k.r), k.r + r.below(g.h - 2 * k.r))) + .collect() +} + +fn time(reps: usize, mut f: impl FnMut() -> T) -> (T, f64) { + let mut last = f(); + let t = Instant::now(); + for _ in 0..reps { + last = f(); + } + (last, t.elapsed().as_nanos() as f64 / reps as f64) +} + +fn pct(v: &mut [usize], p: f64) -> usize { + v.sort_unstable(); + v[((v.len() - 1) as f64 * p) as usize] +} + +// ───────────────────────── sections ───────────────────────── + +/// Kernel falsifiers: sign, single zero crossing at the analytic radius, +/// single annular minimum, symmetry, ring partition, no overflow. +fn kernel_checks() { + println!("== kernel falsifiers =="); + for &(sc, kappa) in &[ + (2.0, 1.5), + (3.0, 1.5), + (2.0, 2.0), + (4.0, 2.0), + (2.0, 3.0), + (3.0, 3.0), + ] { + let k = Kernel::new(sc, kappa); + // Positive centre, negative surround. + assert!(k.lut[0] > 0, "centre must be positive"); + assert!(k.lut.iter().any(|w| *w < 0), "surround must go negative"); + // Exactly one sign change from + to − over q (zeros are a plateau, not a crossing). + let signs: Vec = k + .lut + .iter() + .map(|w| w.signum()) + .filter(|s| *s != 0) + .collect(); + let changes = signs.windows(2).filter(|p| p[0] != p[1]).count(); + assert_eq!(changes, 1, "σc {sc} κ {kappa}: {changes} sign changes"); + // The crossing sits within one q-step of the analytic radius². + let q0 = zero_crossing_q(k.sc, k.ss); + let last_pos = k.lut.iter().rposition(|w| *w > 0).unwrap() as f64; + let first_neg = k.lut.iter().position(|w| *w < 0).unwrap() as f64; + assert!( + last_pos <= q0 + 1.0 && q0 - 1.0 <= first_neg, + "crossing {q0:.2} vs [{last_pos}, {first_neg}]" + ); + // One annular minimum: non-increasing to the minimum, non-decreasing after. + let qmin = (0..=k.r2).min_by_key(|q| (k.lut[*q], *q)).unwrap(); + let down = k.lut[..=qmin].windows(2).filter(|p| p[1] > p[0]).count(); + let up = k.lut[qmin..].windows(2).filter(|p| p[1] < p[0]).count(); + assert_eq!((down, up), (0, 0), "σc {sc} κ {kappa}: false extremum"); + // The 8 square symmetries: every template row is a palindrome, and + // row yi equals column yi. + let wd = k.wd(); + for yi in 0..wd { + assert_eq!(k.disk[yi].reverse_bits() >> (64 - wd), k.disk[yi]); + for xi in 0..wd { + assert_eq!((k.disk[yi] >> xi) & 1, (k.disk[xi] >> yi) & 1); + } + } + // Ring partition: disjoint, and their union is the disk. + for n in [4, 8, 16, 32] { + let rg = rings(&k, n); + for yi in 0..wd { + let mut union = 0u64; + for ring in &rg.ring { + assert_eq!(union & ring[yi], 0, "rings overlap"); + union |= ring[yi]; + } + assert_eq!(union, k.disk[yi], "rings do not cover the disk"); + } + } + // No overflow: the worst case is every positive or every negative cell. + let up: i64 = k.row_up.iter().sum(); + let lo: i64 = k.row_lo.iter().sum(); + assert!(up < i64::from(i32::MAX) && lo > i64::from(i32::MIN)); + println!( + " σc {sc} κ {kappa}: R {:>2} zero crossing q {q0:6.2} (lattice between {last_pos} and {first_neg}) minimum at q {qmin} bounds [{lo}, {up}]", + k.r + ); + } + println!(" sign, one crossing at the analytic radius, one minimum, symmetry, ring partition, no overflow: all hold"); +} + +/// Equality and approximation of every arm on one grid. +fn agreement(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { + let mut wk = Work::default(); + let rgs: Vec<(usize, Rings)> = [8, 16, 32].iter().map(|n| (*n, rings(k, *n))).collect(); + let (mut err_q, mut norm) = (0.0f64, 0.0f64); + let mut err_e = vec![0.0f64; rgs.len()]; + for &(cx, cy) in qs { + let d = arm_d(g, k, cx, cy, &mut wk); + assert_eq!( + arm_f(g, k, cx, cy, &mut wk), + d, + "F differs from D at ({cx},{cy})" + ); + assert_eq!( + arm_g(g, k, cx, cy, &mut wk), + d, + "G differs from D at ({cx},{cy})" + ); + let a = arm_a(g, k, cx, cy, &mut wk); + err_q += (a - d as f64).abs(); + norm += a.abs().max(1.0); + for (i, (_, rg)) in rgs.iter().enumerate() { + err_e[i] += (arm_e(g, k, rg, cx, cy, &mut wk) - d).abs() as f64; + } + } + // The full scan is O(grid) per query; check it on a few centres. + for &(cx, cy) in qs.iter().take(8) { + assert_eq!( + arm_scan(g, k, cx, cy, &mut wk), + arm_d(g, k, cx, cy, &mut wk) + ); + } + let n = qs.len() as f64; + println!( + " {name}: Ascan = D = F = G on every checked centre; |quantised − continuous| mean {:.2} ({:.3} % of |S|) |E − D| mean: {}", + err_q / n, + 100.0 * err_q / norm, + rgs.iter() + .zip(&err_e) + .map(|((n_r, _), e)| format!("K={n_r} {:.1}", e / n)) + .collect::>() + .join(" ") + ); +} + +struct Line { + name: String, + ns: f64, + wk: Work, + queries: usize, +} + +fn line(l: &Line) { + let q = l.queries as f64; + println!( + " {:<26} {:>10.0} ns/query geometry {:>8.1} exp {:>7.1} popcount {:>7.1} rows {:>5.1}", + l.name, + l.ns, + l.wk.geometry as f64 / q, + l.wk.exps as f64 / q, + l.wk.popcounts as f64 / q, + l.wk.rows as f64 / q, + ); +} + +fn bench(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)], scan_queries: usize) { + println!(" -- {name} --"); + let reps = 5; + let run = |label: &str, f: &mut dyn FnMut(&mut Work, usize, usize) -> i64, nq: usize| { + let mut wk = Work::default(); + let a0 = allocs(); + let (_, ns) = time(reps, || { + let mut acc = 0i64; + let mut w = Work::default(); + for &(cx, cy) in &qs[..nq] { + acc = acc.wrapping_add(f(&mut w, cx, cy)); + } + wk = w; + std::hint::black_box(acc) + }); + let a1 = allocs(); + assert_eq!(a1, a0, "{label} allocated"); + line(&Line { + name: label.into(), + ns: ns / nq as f64, + wk, + queries: nq, + }); + }; + let nq = qs.len(); + run( + "A f64 DoG (continuous)", + &mut |w, x, y| arm_a(g, k, x, y, w) as i64, + nq, + ); + run( + "Ascan full candidate scan", + &mut |w, x, y| arm_scan(g, k, x, y, w), + scan_queries, + ); + run( + "D window bits + LUT", + &mut |w, x, y| arm_d(g, k, x, y, w), + nq, + ); + for n in [8, 32] { + let rg = rings(k, n); + run( + &format!("E rings K={n}"), + &mut |w, x, y| arm_e(g, k, &rg, x, y, w), + nq, + ); + } + run( + "F bit-sliced planes", + &mut |w, x, y| arm_f(g, k, x, y, w), + nq, + ); + run( + "G per-row cheaper of D/F", + &mut |w, x, y| arm_g(g, k, x, y, w), + nq, + ); + + // B: the shipped executor. The weight lane is filled per query over the + // row band `cy − R ..= cy + R` (full width; zero outside the disk), then + // `MaskedSumI32` over the resident presence plane on that extent. + let n = g.w * g.h; + let mut lane = vec![0i32; n]; + let program = Program::new( + vec![], + Terminal::MaskedSumI32 { + mask: Operand::Plane(0), + lane: 0, + }, + ); + let mut scratch = Scratch::for_program(&program, n).expect("scratch"); + let mut written = 0u64; + let a0 = allocs(); + let t = Instant::now(); + for _ in 0..reps { + written = 0; + for &(cx, cy) in qs { + let (y0, y1) = (cy - k.r, cy + k.r + 1); + for y in y0..y1 { + let dy = y as i64 - cy as i64; + for x in 0..g.w { + let dx = x as i64 - cx as i64; + let q = dx * dx + dy * dy; + lane[y * g.w + x] = if q <= k.r2 as i64 { + k.lut[q as usize] + } else { + 0 + }; + } + } + written += ((y1 - y0) * g.w * 4) as u64; + let masks: [&[u64]; 1] = [&g.p]; + let lanes = [LaneRef::I32(&lane)]; + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let v = execute_extent( + &program, + &planes, + &Foreign::NONE, + &mut scratch, + Out::None, + y0 * g.w..y1 * g.w, + ); + let Ok(Value::SumI64(s)) = v else { + panic!("B: {v:?}") + }; + assert_eq!( + s, + arm_d(g, k, cx, cy, &mut Work::default()), + "B differs from D" + ); + } + } + let ns = t.elapsed().as_nanos() as f64 / (reps * qs.len()) as f64; + let a1 = allocs(); + println!( + " {:<26} {:>10.0} ns/query weight lane written {:>8} B/query resident lane {} B allocs/query {:.1}", + "B mask-risc MaskedSumI32", + ns, + written / qs.len() as u64, + n * 4, + (a1 - a0) as f64 / (reps * qs.len()) as f64 + ); +} + +/// Arm H against the full answer at several selectivities. +fn early_exit(name: &str, g: &Grid, k: &Kernel, qs: &[(usize, usize)]) { + let order = h_order(k); + let (up, lo) = h_suffix(k, &order); + let full: Vec = qs + .iter() + .map(|&(x, y)| arm_f(g, k, x, y, &mut Work::default())) + .collect(); + let mut sorted = full.clone(); + sorted.sort_unstable(); + println!( + " -- {name}: early exit for S ≥ T (rows in |bound| order, {} rows max) --", + k.wd() + ); + for p in [0.10, 0.50, 0.90, 0.99] { + let t = sorted[((sorted.len() - 1) as f64 * p) as usize]; + let mut rows = Vec::with_capacity(qs.len()); + let mut accepted = 0; + for (&(cx, cy), &s) in qs.iter().zip(&full) { + let (d, r) = arm_h(g, k, &order, &up, &lo, cx, cy, t, &mut Work::default()); + let want = if s >= t { + Decision::Accept + } else { + Decision::Reject + }; + assert_eq!(d, want, "H decided {d:?} at ({cx},{cy}), S {s}, T {t}"); + accepted += usize::from(d == Decision::Accept); + rows.push(r); + } + let mean = rows.iter().sum::() as f64 / rows.len() as f64; + println!( + " T at the {:>2.0} % quantile ({t:>7}): accepted {:>5.1} % rows visited mean {mean:5.1} p95 {:>2} of {}", + p * 100.0, + 100.0 * accepted as f64 / qs.len() as f64, + pct(&mut rows, 0.95), + k.wd() + ); + } + // Timing at the median threshold, against the full F. + let t = sorted[sorted.len() / 2]; + let (_, ns_h) = time(5, || { + qs.iter() + .map(|&(x, y)| arm_h(g, k, &order, &up, &lo, x, y, t, &mut Work::default()).1) + .sum::() + }); + let (_, ns_f) = time(5, || { + qs.iter() + .map(|&(x, y)| arm_f(g, k, x, y, &mut Work::default())) + .sum::() + }); + println!( + " median T: H {:.0} ns/query against F {:.0} ns/query", + ns_h / qs.len() as f64, + ns_f / qs.len() as f64 + ); +} + +/// Counterexamples the early exit must survive. +fn counterexamples(k: &Kernel) { + println!("== counterexamples =="); + let order = h_order(k); + let (up, lo) = h_suffix(k, &order); + let (w, h) = (256, 256); + let (cx, cy) = (128, 128); + + // 1. A partial sum that rises above T and ends below it. The positive + // core is present, and negative cells only in rows the visit order + // reaches late (|dy| > the zero-crossing radius), so the running sum + // climbs first and falls after. The bounds must keep H undecided + // while the sum sits above T. + let q0 = zero_crossing_q(k.sc, k.ss); + let mut g = Grid { + w, + h, + p: vec![0; w * h / 64], + }; + for yi in 0..k.wd() { + for xi in 0..k.wd() { + let (dx, dy) = (xi as i64 - k.r as i64, yi as i64 - k.r as i64); + let q = (dx * dx + dy * dy) as usize; + let core = q <= k.r2 && k.lut[q] > 0; + let late_surround = q <= k.r2 && k.lut[q] < 0 && (dy * dy) as f64 > q0; + if core || late_surround { + g.set(cx + xi - k.r, cy + yi - k.r); + } + } + } + let s = arm_f(&g, k, cx, cy, &mut Work::default()); + let (mut run, mut peak) = (0i64, i64::MIN); + for &yi in &order { + run += f_row( + k, + g.window(cx - k.r, cy + yi - k.r, k.wd()), + yi, + &mut Work::default(), + ); + peak = peak.max(run); + } + let t = (peak + s) / 2; + assert!( + s < t && t < peak, + "the fixture must rise above T and end below it (S {s}, T {t}, peak {peak})" + ); + let (d, rows) = arm_h(&g, k, &order, &up, &lo, cx, cy, t, &mut Work::default()); + assert_eq!(d, Decision::Reject); + println!( + " sign change during evaluation: running sum peaks at {peak}, above T = {t}, and ends at S = {s}; H rejects after {rows} of {} rows", + k.wd() + ); + + // 2. Cancellation: one ring exactly at the zero crossing. S ≈ 0. + let mut g = Grid { + w, + h, + p: vec![0; w * h / 64], + }; + for yi in 0..k.wd() { + for xi in 0..k.wd() { + let (dx, dy) = (xi as i64 - k.r as i64, yi as i64 - k.r as i64); + if ((dx * dx + dy * dy) as f64 - q0).abs() < q0.sqrt() { + g.set(cx + xi - k.r, cy + yi - k.r); + } + } + } + let s = arm_f(&g, k, cx, cy, &mut Work::default()); + let (_, rows) = arm_h(&g, k, &order, &up, &lo, cx, cy, 1, &mut Work::default()); + println!( + " cancellation ring at the zero crossing: S = {s}, decision S ≥ 1 needs {rows} of {} rows", + k.wd() + ); + + // 3. Dense, threshold at the dense mean: no early exit to be had. + let g = random_grid(w, h, 0.9, 0xD3); + let qs = centres(&g, k, 512, 0xC3); + let mut full: Vec = qs + .iter() + .map(|&(x, y)| arm_f(&g, k, x, y, &mut Work::default())) + .collect(); + full.sort_unstable(); + let t = full[full.len() / 2]; + let rows: usize = qs + .iter() + .map(|&(x, y)| arm_h(&g, k, &order, &up, &lo, x, y, t, &mut Work::default()).1) + .sum(); + println!( + " dense ρ = 0.9, T at the median: rows visited mean {:.1} of {} (the bound cannot help here)", + rows as f64 / qs.len() as f64, + k.wd() + ); + + // 4. Disable run, in place: an upper bound that ignores the unvisited + // rows must decide wrongly somewhere. Proves the bound is load-bearing. + let g = random_grid(w, h, 0.2, 0xD4); + let qs = centres(&g, k, 512, 0xC4); + let zero_up = vec![0i64; up.len()]; + let mut wrong = 0; + let s_all: Vec = qs + .iter() + .map(|&(x, y)| arm_f(&g, k, x, y, &mut Work::default())) + .collect(); + let mut sorted = s_all.clone(); + sorted.sort_unstable(); + let t = sorted[sorted.len() * 9 / 10]; + assert!(t > 0, "the disable run needs a positive threshold (T {t})"); + for (&(x, y), &s) in qs.iter().zip(&s_all) { + let (d, _) = arm_h(&g, k, &order, &zero_up, &lo, x, y, t, &mut Work::default()); + wrong += usize::from((d == Decision::Accept) != (s >= t)); + } + assert!( + wrong > 0, + "a bound that drops the unvisited rows must be caught" + ); + println!(" disable run: upper bound without the unvisited rows decides {wrong} of {} queries wrongly", qs.len()); +} + +/// A centre-surround response over geometry cannot see phase: two coherent +/// sources, detector cells bucketed by (r1, r2). With buckets as wide as λ/2, +/// bright and dark fringes share buckets, so pruning by bucket would drop +/// dark fringes along with "uninteresting" cells. +fn interference() { + println!("== interference: geometry buckets versus phase =="); + let lambda = 8.0; + let (s1, s2) = ((96.0, 128.0), (160.0, 128.0)); + for width in [lambda / 2.0, lambda / 8.0] { + let mut cls: std::collections::HashMap<(i64, i64), (bool, bool)> = + std::collections::HashMap::new(); + let (mut dark, mut bright) = (0usize, 0usize); + for y in 0..256 { + for x in 0..256 { + let (xf, yf) = (f64::from(x), f64::from(y)); + let r1 = ((xf - s1.0).powi(2) + (yf - s1.1).powi(2)).sqrt(); + let r2 = ((xf - s2.0).powi(2) + (yf - s2.1).powi(2)).sqrt(); + // Equal-amplitude coherent sources: I / 4I0 = cos²(πΔr/λ). + let i = (std::f64::consts::PI * (r1 - r2) / lambda).cos().powi(2); + let key = ((r1 / width) as i64, (r2 / width) as i64); + let e = cls.entry(key).or_default(); + if i < 0.1 { + e.0 = true; + dark += 1; + } + if i > 0.9 { + e.1 = true; + bright += 1; + } + } + } + let mixed = cls.values().filter(|(d, b)| *d && *b).count(); + let either = cls.values().filter(|(d, b)| *d || *b).count(); + println!( + " bucket width {:>4.1} (λ = {lambda}): {mixed} of {either} occupied (r1, r2) buckets hold both a dark and a bright cell ({dark} dark, {bright} bright cells)", + width + ); + if width >= lambda / 2.0 { + assert!(mixed > 0, "λ/2 buckets must mix fringes"); + } + } +} + +/// Words a (2R+1)² window touches: row-major against Morton order. +fn locality(k: &Kernel) { + fn morton(x: u32, y: u32) -> u64 { + let spread = |v: u32| { + let mut v = u64::from(v); + v = (v | (v << 16)) & 0x0000_FFFF_0000_FFFF; + v = (v | (v << 8)) & 0x00FF_00FF_00FF_00FF; + v = (v | (v << 4)) & 0x0F0F_0F0F_0F0F_0F0F; + v = (v | (v << 2)) & 0x3333_3333_3333_3333; + (v | (v << 1)) & 0x5555_5555_5555_5555 + }; + spread(x) | (spread(y) << 1) + } + let (w, h) = (1024usize, 1024usize); + let g = Grid { w, h, p: vec![] }; + let qs = centres(&g, k, 2048, 0x10C); + let (mut rm, mut mo) = (0usize, 0usize); + for &(cx, cy) in &qs { + let mut a = std::collections::HashSet::new(); + let mut b = std::collections::HashSet::new(); + for y in cy - k.r..=cy + k.r { + for x in cx - k.r..=cx + k.r { + a.insert((y * w + x) / 64); + b.insert(morton(x as u32, y as u32) / 64); + } + } + rm += a.len(); + mo += b.len(); + } + let n = qs.len() as f64; + println!("== locality: u64 words a {0}×{0} window touches ==", k.wd()); + println!( + " row-major {:.1} Morton {:.1} (the window cells are {})", + rm as f64 / n, + mo as f64 / n, + k.wd() * k.wd() + ); +} + +fn main() { + println!( + "D-MHB-1 Mexican-hat bucket probe avx512f={} avx2={}", + cfg!(target_feature = "avx512f"), + cfg!(target_feature = "avx2") + ); + kernel_checks(); + + let k = Kernel::new(3.0, 2.0); + println!( + "\nkernel σc {} σs {} R {} window {}×{} disk cells {}", + k.sc, + k.ss, + k.r, + k.wd(), + k.wd(), + k.disk.iter().map(|r| r.count_ones()).sum::() + ); + + println!("\n== agreement (quantised arms are asserted equal) =="); + for (w, h, label) in [(256usize, 256usize, "64K"), (1024, 1024, "1M")] { + for rho in [0.01, 0.1, 0.5, 0.9] { + let g = random_grid(w, h, rho, 0xA0 + (rho * 100.0) as u64); + let qs = centres(&g, &k, 1024, 0xB0); + agreement(&format!("{label} ρ {rho}"), &g, &k, &qs); + } + } + + println!("\n== speed (ns per query; work counted per query) =="); + for (w, h, label, scan) in [(256usize, 256usize, "64K", 64usize), (1024, 1024, "1M", 8)] { + for rho in [0.01, 0.1, 0.5, 0.9] { + let g = random_grid(w, h, rho, 0xA0 + (rho * 100.0) as u64); + let qs = centres(&g, &k, 2048, 0xB1); + bench(&format!("{label} ρ {rho}"), &g, &k, &qs, scan); + } + let g = wavefront_grid(w, h, w as f64 / 4.0, 0.01, 0xE0); + let qs = centres(&g, &k, 2048, 0xB2); + bench( + &format!("{label} wavefront ring + 1 % noise"), + &g, + &k, + &qs, + scan, + ); + } + + println!("\n== early exit =="); + for rho in [0.01, 0.1, 0.5, 0.9] { + let g = random_grid(1024, 1024, rho, 0xA0 + (rho * 100.0) as u64); + let qs = centres(&g, &k, 4096, 0xB3); + early_exit(&format!("1M ρ {rho}"), &g, &k, &qs); + } + let g = wavefront_grid(1024, 1024, 256.0, 0.01, 0xE0); + let qs = centres(&g, &k, 4096, 0xB4); + early_exit("1M wavefront ring + 1 % noise", &g, &k, &qs); + + println!(); + counterexamples(&k); + println!(); + interference(); + println!(); + locality(&k); + println!("\nall equalities and falsifiers hold"); +}