diff --git a/Cargo.toml b/Cargo.toml index 11bb8905..596da8ee 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -424,6 +424,9 @@ exclude = [ "crates/burn", "crates/wasm-simd-parity", "crates/neon-simd-parity", + # Cross-repo: its dev-dep is a PATH into a lance-graph sibling checkout. + # In-workspace, a missing sibling would fail resolution for EVERY member. + "crates/sigker-parity", "vendor/chacha20", ] default-members = [ diff --git a/crates/sigker-parity/.gitignore b/crates/sigker-parity/.gitignore new file mode 100644 index 00000000..2c96eb1b --- /dev/null +++ b/crates/sigker-parity/.gitignore @@ -0,0 +1,2 @@ +target/ +Cargo.lock diff --git a/crates/sigker-parity/Cargo.toml b/crates/sigker-parity/Cargo.toml new file mode 100644 index 00000000..79ca4448 --- /dev/null +++ b/crates/sigker-parity/Cargo.toml @@ -0,0 +1,32 @@ +# sigker-parity — the cross-repo gate for the Pillar-11 signature lanes. +# +# `ndarray::hpc::pillar::signature::signature_d2_deg3` (hardware: f32, fixed +# d=2/deg-3, Chen accumulation) and lance-graph `sigker::signature_truncated` +# (reference: f64, any d/depth) compute the SAME iterated integrals and had +# zero cross-checks — census finding F-4 of +# `pillar11-signature-certification-unification-v1`. The workspace's own +# architecture rule (ndarray = hardware, lance-graph = thinking) blesses the +# split but demands the parity test it never got. This crate is that test +# (W1), plus the depth-infinity PSD leg (W4) that needs the same sibling. +# +# EXCLUDED from the workspace (see root Cargo.toml `exclude`) because its +# `sigker` dep is a PATH into a sibling checkout. An unconditional path dep +# whose target is absent fails manifest resolution for the WHOLE workspace — +# ndarray CI does not check lance-graph out, so an in-workspace dep here +# would break every ndarray build on a fresh clone. Excluded, it costs +# nothing when the sibling is missing and runs on demand: +# +# cargo test --manifest-path crates/sigker-parity/Cargo.toml +# +# Same shape as `crates/wasm-simd-parity` and `crates/neon-simd-parity`. +[package] +name = "sigker-parity" +version = "0.0.0" +edition = "2021" +publish = false + +[dependencies] +ndarray = { path = "../..", default-features = false, features = ["std", "hpc-extras", "pillar"] } + +[dev-dependencies] +sigker = { path = "../../../lance-graph/crates/sigker" } diff --git a/crates/sigker-parity/examples/w1_diagnose.rs b/crates/sigker-parity/examples/w1_diagnose.rs new file mode 100644 index 00000000..477f2d63 --- /dev/null +++ b/crates/sigker-parity/examples/w1_diagnose.rs @@ -0,0 +1,37 @@ +//! Diagnostic only: is the W1 gap a FORMULA difference or f32 accumulation? +use ndarray::hpc::pillar::signature::signature_d2_deg3; +use sigker::signature_truncated; + +const NAMES: [&str; 15] = [ + "s0", "1x", "1y", "2xx", "2xy", "2yx", "2yy", "3xxx", "3xxy", "3xyx", "3xyy", "3yxx", "3yxy", "3yyx", "3yyy", +]; + +fn cmp(tag: &str, flat: &[f32], n: usize) { + let hw = signature_d2_deg3(flat, n); + let pts: Vec> = (0..n) + .map(|k| vec![flat[2 * k] as f64, flat[2 * k + 1] as f64]) + .collect(); + let refr: Vec = signature_truncated(&pts, 3) + .levels + .iter() + .flat_map(|l| l.iter().copied()) + .collect(); + println!("--- {tag} (n={n}) ---"); + for i in 0..15 { + let (h, r) = (hw[i] as f64, refr[i]); + let d = (h - r).abs(); + let rel = if r.abs() > 1e-12 { d / r.abs() } else { d }; + if rel > 1e-6 { + println!(" {:>5}: hw {:+.9} ref {:+.9} rel {:.3e} <== DIFFERS", NAMES[i], h, r, rel); + } + } +} + +fn main() { + // Single segment: closed form, zero accumulation — any gap here is FORMULA. + cmp("one segment", &[0.0, 0.0, 1.0, 0.5], 2); + // Two segments: Chen composition enters. + cmp("two segments", &[0.0, 0.0, 1.0, 0.5, 1.3, -0.2], 3); + // Three, exact small values (representable in f32) — still formula-only. + cmp("three segments", &[0.0, 0.0, 0.5, 0.25, 0.75, -0.5, 0.25, 0.125], 4); +} diff --git a/crates/sigker-parity/examples/w1_sweep.rs b/crates/sigker-parity/examples/w1_sweep.rs new file mode 100644 index 00000000..09634298 --- /dev/null +++ b/crates/sigker-parity/examples/w1_sweep.rs @@ -0,0 +1,60 @@ +//! Pre-registration sweep: how does the hardware-vs-reference error scale, +//! and under WHICH normalization is it a stable gate? +use ndarray::hpc::pillar::signature::signature_d2_deg3; +use sigker::signature_truncated; + +struct Rng(u64); +impl Rng { + fn f(&mut self) -> f32 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + ((x.wrapping_mul(0x2545_F491_4F6C_DD1D) >> 40) as f32 / (1u32 << 24) as f32) - 0.5 + } +} + +// level of each of the 15 coefficients +const LEVEL: [usize; 15] = [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3]; + +fn main() { + println!("{:>6} {:>12} {:>14} {:>16}", "N", "worst |abs|", "worst /coeff", "worst /levelmax"); + for &n in &[16usize, 32, 64, 128, 256] { + let mut rng = Rng(0x9E37_79B9_7F4A_7C15); + let (mut wa, mut wc, mut wl) = (0.0f64, 0.0f64, 0.0f64); + for _ in 0..1000 { + let (mut x, mut y) = (0.0f32, 0.0f32); + let mut flat = Vec::with_capacity(n * 2); + let mut pts = Vec::with_capacity(n); + for _ in 0..n { + flat.push(x); + flat.push(y); + pts.push(vec![x as f64, y as f64]); + x += rng.f(); + y += rng.f(); + } + let hw = signature_d2_deg3(&flat, n); + let refr: Vec = signature_truncated(&pts, 3) + .levels + .iter() + .flat_map(|l| l.iter().copied()) + .collect(); + // characteristic magnitude per level, from the REFERENCE + let mut lvmax = [0.0f64; 4]; + for i in 0..15 { + lvmax[LEVEL[i]] = lvmax[LEVEL[i]].max(refr[i].abs()); + } + for i in 0..15 { + let d = (hw[i] as f64 - refr[i]).abs(); + wa = wa.max(d); + if refr[i].abs() > 1e-12 { + wc = wc.max(d / refr[i].abs()); + } + let s = lvmax[LEVEL[i]].max(1e-12); + wl = wl.max(d / s); + } + } + println!("{n:>6} {wa:>12.3e} {wc:>14.3e} {wl:>16.3e}"); + } +} diff --git a/crates/sigker-parity/examples/w4_concentration_sweep.rs b/crates/sigker-parity/examples/w4_concentration_sweep.rs new file mode 100644 index 00000000..5d1685e1 --- /dev/null +++ b/crates/sigker-parity/examples/w4_concentration_sweep.rs @@ -0,0 +1,63 @@ +//! Does depth-inf self-kernel concentration shrink like 1/sqrt(N) (a sample- +//! size effect) or plateau (a genuine heavy tail)? Measure, do not assume. +use ndarray::hpc::pillar::signature::brownian_path_d2; +use ndarray::hpc::pillar::SplitMix64; +use sigker::{signature_kernel_pde, signature_truncated}; + +const SEED: u64 = 0x5EED_1111_5164_A7AB; +const N_STEPS: usize = 50; + +fn pool(n: usize) -> Vec>> { + let mut rng = SplitMix64::new(SEED); + (0..n) + .map(|_| { + let p = brownian_path_d2(&mut rng, N_STEPS); + (0..=N_STEPS) + .map(|k| vec![p[2 * k] as f64, p[2 * k + 1] as f64]) + .collect() + }) + .collect() +} + +fn stats(v: &[f64]) -> (f64, f64, f64) { + let n = v.len(); + let h = n / 2; + let m1 = v[..h].iter().sum::() / h as f64; + let m2 = v[h..].iter().sum::() / (n - h) as f64; + let mean = v.iter().sum::() / n as f64; + let var = v.iter().map(|x| (x - mean).powi(2)).sum::() / n as f64; + // half-mean gap, and the coefficient of variation that predicts it + ((m1 - m2).abs() / mean, var.sqrt() / mean, mean) +} + +fn main() { + println!("depth-INFINITY (Goursat PDE)"); + println!("{:>6} {:>12} {:>10} {:>14} {:>12}", "N", "concentr", "CV", "predicted", "mean K"); + for &n in &[64usize, 128, 256, 512, 1000] { + let p = pool(n); + let k: Vec = p.iter().map(|x| signature_kernel_pde(x, x)).collect(); + let (c, cv, mean) = stats(&k); + // For independent samples the expected half-mean gap ~ CV * sqrt(8/(pi*N)) + let pred = cv * (8.0 / (core::f64::consts::PI * n as f64)).sqrt(); + println!("{n:>6} {c:>12.4} {cv:>10.3} {pred:>14.4} {mean:>12.4e}"); + } + println!("\ndepth-3 TRUNCATED (the existing battery's kernel, f64 reference)"); + println!("{:>6} {:>12} {:>10} {:>14} {:>12}", "N", "concentr", "CV", "predicted", "mean K"); + for &n in &[64usize, 1000] { + let p = pool(n); + let k: Vec = p + .iter() + .map(|x| { + let s = signature_truncated(x, 3); + s.levels + .iter() + .flat_map(|l| l.iter()) + .map(|v| v * v) + .sum::() + }) + .collect(); + let (c, cv, mean) = stats(&k); + let pred = cv * (8.0 / (core::f64::consts::PI * n as f64)).sqrt(); + println!("{n:>6} {c:>12.4} {cv:>10.3} {pred:>14.4} {mean:>12.4e}"); + } +} diff --git a/crates/sigker-parity/src/lib.rs b/crates/sigker-parity/src/lib.rs new file mode 100644 index 00000000..3b3425ff --- /dev/null +++ b/crates/sigker-parity/src/lib.rs @@ -0,0 +1,3 @@ +pub fn sibling_is_wired() -> bool { + true +} diff --git a/crates/sigker-parity/tests/w1_signature_parity.rs b/crates/sigker-parity/tests/w1_signature_parity.rs new file mode 100644 index 00000000..f5b24e0b --- /dev/null +++ b/crates/sigker-parity/tests/w1_signature_parity.rs @@ -0,0 +1,231 @@ +//! W1 — the parity bridge (census F-4). +//! +//! `ndarray::…::signature_d2_deg3` (hardware: f32, fixed d=2/deg-3) and +//! `sigker::signature_truncated` (reference: f64, any d/depth) implement the +//! same iterated integrals in two repos with, until this file, zero +//! cross-checks. +//! +//! # Layout correspondence (the thing that makes the comparison meaningful) +//! +//! ndarray returns one flat `[f32; 15]`: +//! `[s0, s1x, s1y, s2xx, s2xy, s2yx, s2yy, s3xxx, s3xxy, s3xyx, s3xyy, +//! s3yxx, s3yxy, s3yyx, s3yyy]`. +//! +//! sigker returns per-level flat storage, `levels[k]` of length `d^k`, +//! row-major in the index tuple. For d = 2 that is `[1]`, `[x, y]`, +//! `[xx, xy, yx, yy]`, `[xxx, xxy, xyx, xyy, yxx, yxy, yyx, yyy]` — the same +//! order, concatenated. So the comparison is index-for-index, and +//! `concat_levels` asserts the total length is 15 rather than assuming it. +//! +//! # Why the gate is normalized by LEVEL, not by coefficient +//! +//! The obvious gate — per-coefficient relative error — is the cancellation +//! trap this plan's §6 already names as law, one level down from where it was +//! first met (the D-SK kernel-scalar finding). Measured by +//! `examples/w1_sweep.rs` over 1000 paths at each length: +//! +//! ```text +//! N worst |abs| worst /coeff worst /levelmax +//! 16 8.741e-7 2.010e2 6.452e-6 +//! 32 3.865e-6 1.827e3 2.439e-6 +//! 64 1.727e-5 7.911e2 4.154e-6 +//! 128 7.587e-5 7.779e2 4.070e-6 +//! 256 2.988e-4 7.926e3 9.032e-6 +//! ``` +//! +//! Per-coefficient relative error swings over 2e2..8e3 with no trend — a +//! level-3 coefficient of a random walk passes through zero by cancellation, +//! and relative error against a near-zero denominator is unbounded no matter +//! how correct the implementation is. Absolute error is no better: it grows +//! ~N² with the signature's own scale. Normalized by the characteristic +//! magnitude of the coefficient's OWN LEVEL, the error is flat at +//! 2.4e-6..9.0e-6 across a 16x range of path length — that is a property of +//! the implementation, so that is what the gate binds on. +//! +//! `REL_TOL = 1e-4` therefore sits ~11x above the worst measured value. It is +//! a pre-registered bound with margin, not a fitted line; the sweep is +//! committed alongside so a future tightening has its evidence. +//! +//! And the formula itself is EXACT: on a single segment (closed form, zero +//! accumulation) the two implementations agree to the last f32 bit — see +//! `examples/w1_diagnose.rs`. What this gate bounds is f32 accumulation +//! drift, nothing else. + +use ndarray::hpc::pillar::signature::{signature_d2_deg3, SIG_D2_DEG3_LEN}; +use sigker::signature_truncated; + +const N_PATHS: usize = 1000; +const N_POINTS: usize = 64; +/// Error bound, normalized by the characteristic magnitude of each +/// coefficient's own level. Pre-registered from the sweep above (worst +/// measured 9.03e-6 at N=256) with ~11x margin. +const REL_TOL: f64 = 1e-4; +/// Which signature level each of the 15 coefficients belongs to. +const LEVEL: [usize; SIG_D2_DEG3_LEN] = [0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3]; + +/// Deterministic path generator — a tiny xorshift so the fixture needs no +/// dev-dep and reproduces byte-for-byte across runs and machines. +struct Rng(u64); +impl Rng { + fn next_f32(&mut self) -> f32 { + // xorshift64* + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + let v = x.wrapping_mul(0x2545_F491_4F6C_DD1D) >> 40; // 24 bits + (v as f32 / (1u32 << 24) as f32) - 0.5 + } +} + +/// One random d=2 walk, in both repos' input shapes. +fn make_path(rng: &mut Rng, n_points: usize) -> (Vec, Vec>) { + let mut flat = Vec::with_capacity(n_points * 2); + let mut pts = Vec::with_capacity(n_points); + let (mut x, mut y) = (0.0f32, 0.0f32); + for _ in 0..n_points { + flat.push(x); + flat.push(y); + pts.push(vec![x as f64, y as f64]); + x += rng.next_f32(); + y += rng.next_f32(); + } + (flat, pts) +} + +/// sigker's per-level storage, flattened into ndarray's single-array order. +fn concat_levels(sig: &sigker::Signature) -> Vec { + let out: Vec = sig.levels.iter().flat_map(|l| l.iter().copied()).collect(); + assert_eq!( + out.len(), + SIG_D2_DEG3_LEN, + "layout correspondence broken: sigker d=2/depth=3 must flatten to \ + exactly the {SIG_D2_DEG3_LEN} coefficients ndarray returns" + ); + out +} + +/// Worst level-normalized error over `n_paths` random walks. +/// +/// The denominator is the largest reference coefficient WITHIN the same +/// level, so a coefficient that cancels toward zero is measured against the +/// scale its level actually carries rather than against its own vanishing +/// magnitude. +fn worst_level_error(n_paths: usize, seed: u64) -> (f64, usize) { + let mut rng = Rng(seed); + let mut worst = 0.0f64; + let mut worst_idx = 0usize; + for _ in 0..n_paths { + let (flat, pts) = make_path(&mut rng, N_POINTS); + let hw = signature_d2_deg3(&flat, N_POINTS); + let refr = concat_levels(&signature_truncated(&pts, 3)); + for (i, (&h, &r)) in hw.iter().zip(refr.iter()).enumerate() { + let err = (h as f64 - r).abs() / level_scale(&refr, LEVEL[i]); + if err > worst { + worst = err; + worst_idx = i; + } + } + } + (worst, worst_idx) +} + +/// Characteristic magnitude of one signature level, from the f64 reference. +fn level_scale(refr: &[f64], level: usize) -> f64 { + refr.iter() + .enumerate() + .filter(|(i, _)| LEVEL[*i] == level) + .map(|(_, v)| v.abs()) + .fold(0.0f64, f64::max) + .max(1e-12) +} + +#[test] +fn hardware_signature_matches_the_sigker_reference() { + let (worst, idx) = worst_level_error(N_PATHS, 0x9E37_79B9_7F4A_7C15); + assert!( + worst < REL_TOL, + "W1 parity FAILED: worst level-normalized error {worst:.3e} \ + (coefficient index {idx}) exceeds {REL_TOL:.0e} over {N_PATHS} paths" + ); +} + +/// The margin is REPORTED, not merely asserted — a gate whose measured value +/// is invisible cannot be re-tightened later, and a silent 100x margin is +/// indistinguishable from a test that compares nothing. +#[test] +fn parity_margin_is_reported() { + let (worst, idx) = worst_level_error(N_PATHS, 0xD1B5_4A32_D192_ED03); + println!( + "W1 parity margin: worst level-normalized err {worst:.3e} at coefficient {idx}, \ + bound {REL_TOL:.0e}, margin {:.1}x", + REL_TOL / worst.max(f64::MIN_POSITIVE) + ); + assert!(worst < REL_TOL, "second seed must hold the same bound"); +} + +// ════════════════════════════════════════════════════════════════════════════ +// Anti-vacuity half — the gate must be able to FAIL. +// +// A parity assertion is only evidence if some reachable implementation +// violates it. `sabotaged_signature` is the same Chen accumulation with ONE +// term removed: the `½ dx_i dx_j` self-term in the level-2 update (and its +// level-3 companions). That is the single most plausible way to get this +// wrong — it is exactly the term that distinguishes the true iterated +// integral from the naive "outer product of increments" — and dropping it +// still produces finite, same-magnitude, plausible-looking coefficients. +// If REL_TOL still passed on that, the test above would be measuring nothing. +// ════════════════════════════════════════════════════════════════════════════ + +fn sabotaged_signature(path: &[f32], n_points: usize) -> [f32; SIG_D2_DEG3_LEN] { + let (mut s1x, mut s1y) = (0.0f32, 0.0f32); + let (mut s2xx, mut s2xy, mut s2yx, mut s2yy) = (0.0f32, 0.0f32, 0.0f32, 0.0f32); + let mut s3 = [0.0f32; 8]; + for k in 0..n_points - 1 { + let dx = path[2 * k + 2] - path[2 * k]; + let dy = path[2 * k + 3] - path[2 * k + 1]; + // Level 3 — unchanged from the real update. + s3[0] += s2xx * dx + 0.5 * s1x * dx * dx + (1.0 / 6.0) * dx * dx * dx; + s3[1] += s2xx * dy + 0.5 * s1x * dx * dy + (1.0 / 6.0) * dx * dx * dy; + s3[2] += s2xy * dx + 0.5 * s1x * dy * dx + (1.0 / 6.0) * dx * dy * dx; + s3[3] += s2xy * dy + 0.5 * s1x * dy * dy + (1.0 / 6.0) * dx * dy * dy; + s3[4] += s2yx * dx + 0.5 * s1y * dx * dx + (1.0 / 6.0) * dy * dx * dx; + s3[5] += s2yx * dy + 0.5 * s1y * dx * dy + (1.0 / 6.0) * dy * dx * dy; + s3[6] += s2yy * dx + 0.5 * s1y * dy * dx + (1.0 / 6.0) * dy * dy * dx; + s3[7] += s2yy * dy + 0.5 * s1y * dy * dy + (1.0 / 6.0) * dy * dy * dy; + // THE SABOTAGE: the `+ 0.5 * d_i * d_j` self-term is dropped here. + s2xx += s1x * dx; + s2xy += s1x * dy; + s2yx += s1y * dx; + s2yy += s1y * dy; + s1x += dx; + s1y += dy; + } + [1.0, s1x, s1y, s2xx, s2xy, s2yx, s2yy, s3[0], s3[1], s3[2], s3[3], s3[4], s3[5], s3[6], s3[7]] +} + +#[test] +fn a_wrong_chen_accumulation_is_caught_by_the_same_bound() { + let mut rng = Rng(0x9E37_79B9_7F4A_7C15); + let mut worst = 0.0f64; + for _ in 0..N_PATHS { + let (flat, pts) = make_path(&mut rng, N_POINTS); + let bad = sabotaged_signature(&flat, N_POINTS); + let refr = concat_levels(&signature_truncated(&pts, 3)); + for (i, (&h, &r)) in bad.iter().zip(refr.iter()).enumerate() { + worst = worst.max((h as f64 - r).abs() / level_scale(&refr, LEVEL[i])); + } + } + assert!( + worst > REL_TOL, + "ANTI-VACUITY FAILED: the sabotaged accumulation passed the parity \ + bound (worst level-normalized err {worst:.3e} < {REL_TOL:.0e}) — the \ + gate above therefore proves nothing" + ); + println!( + "W1 anti-vacuity: sabotaged accumulation reaches level-normalized err \ + {worst:.3e}, {:.0}x the bound {REL_TOL:.0e}", + worst / REL_TOL + ); +} diff --git a/crates/sigker-parity/tests/w4_depth_infinity_psd.rs b/crates/sigker-parity/tests/w4_depth_infinity_psd.rs new file mode 100644 index 00000000..e7f1578a --- /dev/null +++ b/crates/sigker-parity/tests/w4_depth_infinity_psd.rs @@ -0,0 +1,246 @@ +//! W4 — PSD at depth-infinity (census M-3). +//! +//! ndarray's `prove_pillar_11` certifies the TRUNCATED (d=2, deg-3) kernel's +//! stability. The property kernel machines actually rely on — PSD-ness of the +//! Gram — is uncertified for the depth-infinity kernel everywhere. This leg +//! re-runs that machinery over `sigker::signature_kernel_pde` values on the +//! SAME Brownian pool, per ruling Q2 (cross-repo call, not a second f32 +//! Goursat port; a port waits for W5's trigger). +//! +//! # Strengthened relative to the truncated battery +//! +//! `prove_pillar_11`'s "PSD" criteria are diagonal positivity plus +//! Cauchy-Schwarz. Both are NECESSARY conditions, neither is sufficient: a +//! matrix can satisfy them and still have a negative eigenvalue. This leg +//! adds the sufficient test — a Cholesky factorization, which exists iff the +//! matrix is positive definite — and pairs it with the falsifier that the +//! weaker criteria cannot supply (see `a_non_psd_gram_is_rejected`). +//! +//! # Tolerance +//! +//! The Gram is built in f64 from a first-order solver, so the diagonal is not +//! exact; Cholesky is run with a relative jitter of `JITTER` times the mean +//! diagonal, the standard numerical allowance. The falsifier below confirms +//! that this jitter does not swallow a genuinely indefinite matrix. + +use ndarray::hpc::pillar::signature::brownian_path_d2; +use ndarray::hpc::pillar::SplitMix64; +use sigker::signature_kernel_pde; + +/// Same seed as the truncated battery, so both certify the same pool. +const PILLAR_11_SEED: u64 = 0x5EED_1111_5164_A7AB; +/// Gram pool — O(N^2) kernel solves, so kept small. +const N_PATHS: usize = 64; +/// Concentration pool — O(N) solves, matched to the truncated battery's 1000 +/// so the two numbers are directly comparable (see that leg's doc). +const N_CONC: usize = 1000; +const N_STEPS: usize = 50; +const JITTER: f64 = 1e-9; + +/// The f32 hardware path buffer, in sigker's point-list shape. +fn as_points(flat: &[f32], n_points: usize) -> Vec> { + (0..n_points) + .map(|k| vec![flat[2 * k] as f64, flat[2 * k + 1] as f64]) + .collect() +} + +fn brownian_pool(seed: u64, n_paths: usize) -> Vec>> { + let mut rng = SplitMix64::new(seed); + (0..n_paths) + .map(|_| { + let p = brownian_path_d2(&mut rng, N_STEPS); + as_points(&p, N_STEPS + 1) + }) + .collect() +} + +fn gram(pool: &[Vec>]) -> Vec> { + let n = pool.len(); + let mut k = vec![vec![0.0f64; n]; n]; + for i in 0..n { + for j in i..n { + let v = signature_kernel_pde(&pool[i], &pool[j]); + k[i][j] = v; + k[j][i] = v; // the kernel is symmetric by construction + } + } + k +} + +/// Cholesky with relative jitter. `Some(_)` iff the matrix is positive +/// definite to that tolerance — the SUFFICIENT PSD test the truncated +/// battery lacks. +fn cholesky(k: &[Vec], jitter: f64) -> Option { + let n = k.len(); + let mean_diag = (0..n).map(|i| k[i][i]).sum::() / n as f64; + let eps = jitter * mean_diag.abs().max(1.0); + let mut l = vec![vec![0.0f64; n]; n]; + for i in 0..n { + for j in 0..=i { + let mut s = k[i][j]; + if i == j { + s += eps; + } + // Dot of the two already-computed row prefixes. Both borrows are + // immutable and end before the write below. + let dot: f64 = l[i][..j].iter().zip(&l[j][..j]).map(|(a, b)| a * b).sum(); + s -= dot; + if i == j { + if s <= 0.0 { + return Some(i); // failed at this leading minor + } + l[i][i] = s.sqrt(); + } else { + l[i][j] = s / l[j][j]; + } + } + } + None +} + +#[test] +fn the_depth_infinity_gram_is_positive_definite() { + let pool = brownian_pool(PILLAR_11_SEED, N_PATHS); + let k = gram(&pool); + + // Necessary conditions first — the same two the truncated battery uses, + // so a failure here is directly comparable to that report. + let mut cs_violations = 0usize; + for i in 0..N_PATHS { + assert!(k[i][i] > 0.0, "depth-inf self-kernel K[{i},{i}] = {} is not positive", k[i][i]); + for j in i + 1..N_PATHS { + if k[i][j] * k[i][j] > k[i][i] * k[j][j] * 1.001 { + cs_violations += 1; + } + } + } + assert_eq!(cs_violations, 0, "Cauchy-Schwarz violated at depth-infinity"); + + // The sufficient one. + assert!( + cholesky(&k, JITTER).is_none(), + "depth-infinity Gram is NOT positive definite (Cholesky failed at \ + leading minor {:?}) — M-3 would be falsified", + cholesky(&k, JITTER) + ); + println!("W4: depth-inf Gram over {N_PATHS} Brownian paths — diag > 0, Cauchy-Schwarz clean, Cholesky OK"); +} + +/// Concentration is measured at `N_CONC`, NOT at `N_PATHS`. +/// +/// The 0.20 bound comes from the truncated battery, which runs 1000 paths. +/// Applying it to the 64-path Gram pool would silently be a DIFFERENT gate — +/// measured, at N = 64 the truncated kernel itself concentrates to 0.3461 and +/// would fail its own bound. Half-mean agreement is a sample-size statistic +/// before it is a kernel property, so the comparison is only meaningful at +/// matched N (`examples/w4_concentration_sweep.rs`): +/// +/// ```text +/// depth-INFINITY depth-3 TRUNCATED +/// N concentr N concentr +/// 64 0.2892 64 0.3461 +/// 128 0.2156 1000 0.0481 +/// 256 0.1884 +/// 512 0.0053 +/// 1000 0.0038 +/// ``` +/// +/// At matched N = 1000 the depth-infinity kernel concentrates BETTER than the +/// truncated one it extends (0.0038 vs 0.0481). The Gram/Cholesky leg above +/// stays at 64 paths because it costs O(N^2) solves; this leg is O(N). +#[test] +fn the_depth_infinity_self_kernel_concentrates() { + let pool = brownian_pool(PILLAR_11_SEED, N_CONC); + let self_k: Vec = pool.iter().map(|p| signature_kernel_pde(p, p)).collect(); + let concentration = half_mean_gap(&self_k); + println!( + "W4: depth-inf self-kernel concentration {concentration:.4} over {N_CONC} paths \ + (bound 0.20, the truncated battery's; truncated scores 0.0481 at the same N)" + ); + assert!( + concentration < 0.20, + "depth-inf self-kernel half-means disagree by {concentration:.3} at N = {N_CONC}" + ); +} + +/// Can-stay-silent's partner: a pool that genuinely does NOT concentrate must +/// be caught at the SAME N, so the green above is a property of the kernel +/// and not of the bound being loose. The fixture rescales the second half of +/// the pool, which is exactly the drift a half-mean statistic exists to see. +#[test] +fn a_non_concentrating_pool_is_caught() { + let pool = brownian_pool(PILLAR_11_SEED, N_CONC); + let self_k: Vec = pool + .iter() + .enumerate() + .map(|(i, p)| { + if i < N_CONC / 2 { + signature_kernel_pde(p, p) + } else { + let scaled: Vec> = p + .iter() + .map(|pt| pt.iter().map(|v| v * 1.35).collect()) + .collect(); + signature_kernel_pde(&scaled, &scaled) + } + }) + .collect(); + let concentration = half_mean_gap(&self_k); + assert!( + concentration >= 0.20, + "ANTI-VACUITY FAILED: a pool whose second half is rescaled 1.35x still \ + concentrated to {concentration:.4} — the bound cannot see drift" + ); + println!("W4 anti-vacuity: rescaled-half pool reaches concentration {concentration:.4} (bound 0.20)"); +} + +/// Relative gap between the two half-means — the truncated battery's statistic. +fn half_mean_gap(v: &[f64]) -> f64 { + let h = v.len() / 2; + let m1 = v[..h].iter().sum::() / h as f64; + let m2 = v[h..].iter().sum::() / (v.len() - h) as f64; + let combined = v.iter().sum::() / v.len() as f64; + if combined > 0.0 { + (m1 - m2).abs() / combined + } else { + f64::INFINITY + } +} + +// ════════════════════════════════════════════════════════════════════════════ +// Falsifier pair for the PSD gate. +// ════════════════════════════════════════════════════════════════════════════ + +/// Can-fire: an indefinite matrix of the same shape and scale must be +/// REJECTED. Without this, `cholesky` returning `None` would be evidence of +/// nothing — and it is precisely the test the truncated battery cannot run, +/// because diagonal-positivity + Cauchy-Schwarz both HOLD on this matrix. +#[test] +fn a_non_psd_gram_is_rejected() { + let pool = brownian_pool(PILLAR_11_SEED, 8); + let mut k = gram(&pool); + // Flip the sign of one off-diagonal pair. Diagonals are untouched (still + // positive) and the magnitudes are unchanged (so Cauchy-Schwarz still + // holds) — only definiteness is destroyed. + let scale = (k[0][0] * k[1][1]).sqrt(); + k[0][1] = -0.999 * scale; + k[1][0] = -0.999 * scale; + k[0][2] = 0.999 * (k[0][0] * k[2][2]).sqrt(); + k[2][0] = k[0][2]; + k[1][2] = 0.999 * (k[1][1] * k[2][2]).sqrt(); + k[2][1] = k[1][2]; + + // The weak criteria still pass — that is the point. + for i in 0..3 { + assert!(k[i][i] > 0.0, "diagonal must stay positive in the fixture"); + for j in i + 1..3 { + assert!(k[i][j] * k[i][j] <= k[i][i] * k[j][j] * 1.001, "Cauchy-Schwarz must still hold in the fixture"); + } + } + // The strong one must catch it. + assert!(cholesky(&k, JITTER).is_some(), "ANTI-VACUITY FAILED: an indefinite Gram passed the Cholesky gate"); + println!( + "W4 anti-vacuity: indefinite Gram passes diag>0 AND Cauchy-Schwarz, rejected by Cholesky at minor {:?}", + cholesky(&k, JITTER) + ); +}