From 1234df3b5f3d226cb83a948ce5d2c6b425026973 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 19:43:38 +0000 Subject: [PATCH 01/14] mask-risc + jc: grouped moments terminal and one-way ANOVA from moments MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mask-risc: Terminal::GroupMomentsI32 { mask, key: GroupKey, val } folds (n, Σx, Σx²) per group into a new Out::Moments sink, one delegation per tile to ndarray::simd::masked_group_moments_i32{,_via,_pair}. It is its own terminal, not a GroupFold member: a GroupFold slot is one seeded i64. Validation refuses a wrong-width key/value lane, a missing or wrong-shaped sink, a plane past the 2^32-row exactness bound, and a partial extent. The independent oracle folds straight to i128; the key resolution it shares with GroupReduce is now one longhand helper (row_group), and the data type reaches it through value.rs so the oracle still never names the SIMD facade (law L4, guarded by the_oracle_has_no_facade_token). jc: anova_from_moments / eta_squared_from_moments project the same statistics as anova_one_way / eta_squared from grouped moments. The F / p / η² policy and every None case now live once, in anova_from_ss and eta_from_ss, shared by both paths. Sums of squares are formed from exact i128 quantities (W_g = n·Σx² − (Σx)², D_g = N·S_g − n_g·S), so no two large floats are subtracted and zero within-group variance is detected exactly. Tests: executor == oracle for GroupMomentsI32 on resident, VIA and pair keys across tiles; refusal cases; and in jc, mask -> fold -> moments ANOVA vs materialized -> anova_one_way over arbitrary masks, unequal groups, negative values, values near the i32 bound, chunk-merge identity, degenerate layouts, and an adversarial layout where the two-pass slice path loses the within-group variance. quack's duckdb_differential gains the unreachable arm for the new Value. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/jc/src/stats.rs | 403 +++++++++++++++++- crates/lance-graph-mask-risc/src/exec.rs | 38 +- crates/lance-graph-mask-risc/src/ir.rs | 24 ++ crates/lance-graph-mask-risc/src/lib.rs | 2 +- crates/lance-graph-mask-risc/src/reference.rs | 192 ++++++--- crates/lance-graph-mask-risc/src/value.rs | 11 + crates/lance-graph-mask-risc/tests/foreign.rs | 203 ++++++++- .../tests/duckdb_differential.rs | 1 + 8 files changed, 804 insertions(+), 70 deletions(-) diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index d72775926..9700becab 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -74,6 +74,7 @@ //! front by the same `all_finite` guard. use crate::reliability::{all_finite, mean, pearson}; +use ndarray::simd::GroupMoments; use std::collections::BTreeSet; // ─────────────────────────── local helpers ─────────────────────────── @@ -947,6 +948,12 @@ fn one_way_ss(groups: &[Vec]) -> Option<(f64, f64, usize, usize)> { /// ``` pub fn eta_squared(groups: &[Vec]) -> Option { let (ss_b, ss_t, _, _) = one_way_ss(groups)?; + eta_from_ss(ss_b, ss_t) +} + +/// η² from the two sums of squares — the ONE place its degeneracy policy +/// lives, shared by [`eta_squared`] and [`eta_squared_from_moments`]. +fn eta_from_ss(ss_b: f64, ss_t: f64) -> Option { if ss_t == 0.0 || !ss_t.is_finite() { return None; } @@ -1134,12 +1141,19 @@ pub fn t_test_student(a: &[f64], b: &[f64]) -> Option { /// ``` pub fn anova_one_way(groups: &[Vec]) -> Option { let (ss_b, ss_t, k, n_total) = one_way_ss(groups)?; + anova_from_ss(ss_b, ss_t - ss_b, ss_t, k, n_total) +} + +/// The F test from the three sums of squares — the ONE place the one-way +/// ANOVA's degeneracy policy lives, shared by [`anova_one_way`] (which +/// passes `ss_w = ss_t − ss_b`) and [`anova_from_moments`] (which forms +/// `ss_w` exactly and `ss_t = ss_b + ss_w`). +fn anova_from_ss(ss_b: f64, ss_w: f64, ss_t: f64, k: usize, n_total: usize) -> Option { if n_total <= k { return None; // no within-group df } let df_b = (k - 1) as f64; let df_w = (n_total - k) as f64; - let ss_w = ss_t - ss_b; if ss_w <= 0.0 || !ss_w.is_finite() { // ≤ 0 → no within-group variance (F undefined / degenerate). return None; @@ -1166,6 +1180,108 @@ pub fn anova_one_way(groups: &[Vec]) -> Option { }) } +// ──────────── one-way ANOVA from grouped sufficient statistics ──────────── +// +// The same statistics as `anova_one_way` / `eta_squared`, projected from +// per-group `(n, Σx, Σx²)` instead of from materialized group vectors — the +// shape `ndarray::simd::masked_group_moments_i32` folds out of a population +// mask in one pass. The statistical policy (every `None` case, the F / p / +// η² formulas) is shared with the slice path through `anova_from_ss` and +// `eta_from_ss`; only the sums of squares are formed differently. + +/// Between- and within-group sums of squares from grouped moments, as +/// `(ss_b, ss_w, ss_t, k, n_total)`. +/// +/// Formed in exact integer arithmetic as far as the division allows, so no +/// two large floats are ever subtracted: +/// +/// - per group, `W_g = n_g·Σx² − (Σx)² = n_g · Σ(x − x̄_g)²` is an exact `i128` +/// (never negative), and `ss_w = Σ W_g / n_g`; +/// - per group, `D_g = N·S_g − n_g·S` is an exact `i128`, and +/// `ss_b = Σ D_g² / (n_g·N²)` — the textbook `Σ n_g (x̄_g − x̄)²` with the +/// means cleared of their denominators; +/// - `ss_t = ss_b + ss_w`, a sum of two non-negative terms. +/// +/// So `ss_w == 0.0` exactly when every group is constant, rather than +/// whenever rounding happens to cancel. +/// +/// `None` on fewer than 2 groups, any empty group (`n == 0`, matching +/// [`anova_one_way`]'s empty-group rule), or moments too large for the exact +/// `i128` forms — i.e. outside the `2^32`-rows-per-group bound +/// `GroupMoments` is exact under. +fn one_way_ss_from_moments(groups: &[GroupMoments]) -> Option<(f64, f64, f64, usize, usize)> { + if groups.len() < 2 || groups.iter().any(|g| g.n == 0) { + return None; + } + let mut n_total: u64 = 0; + let mut s_total: i128 = 0; + for g in groups { + n_total = n_total.checked_add(g.n)?; + s_total = s_total.checked_add(i128::from(g.sum))?; + } + let big_n = i128::from(n_total); + let n_total_f = n_total as f64; + let mut ss_w = 0.0f64; + let mut ss_b = 0.0f64; + for g in groups { + let n_g = i128::from(g.n); + let s_g = i128::from(g.sum); + let w_g = n_g + .checked_mul(g.sum_sq)? + .checked_sub(s_g.checked_mul(s_g)?)?; + if w_g < 0 { + return None; // impossible for moments of real data: inconsistent input + } + let d_g = big_n + .checked_mul(s_g)? + .checked_sub(n_g.checked_mul(s_total)?)?; + let n_g_f = g.n as f64; + ss_w += w_g as f64 / n_g_f; + let d = d_g as f64; + ss_b += d * d / (n_g_f * n_total_f * n_total_f); + } + let n_total = usize::try_from(n_total).ok()?; + Some((ss_b, ss_w, ss_b + ss_w, groups.len(), n_total)) +} + +/// One-way ANOVA from grouped sufficient statistics — the same result as +/// [`anova_one_way`] on the materialized groups, without the groups. +/// +/// `groups[g]` holds group `g`'s `(n, Σx, Σx²)`, e.g. as folded out of a +/// population mask by `ndarray::simd::masked_group_moments_i32` (or +/// `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32`). Every degenerate +/// case is [`anova_one_way`]'s: fewer than 2 groups, an empty group, +/// `N ≤ k`, and zero within-group variance are all `None`. The sums of +/// squares are formed exactly (see `one_way_ss_from_moments`), so on data +/// where the slice path's two-pass floats cancel badly this is the more +/// accurate of the two, not merely an approximation of it. +/// +/// ``` +/// use jc::stats::{anova_from_moments, anova_one_way}; +/// use ndarray::simd::GroupMoments; +/// +/// let groups = vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]; +/// let moments = [ +/// GroupMoments { n: 3, sum: 6, sum_sq: 14 }, +/// GroupMoments { n: 3, sum: 15, sum_sq: 77 }, +/// ]; +/// let (a, b) = (anova_one_way(&groups).unwrap(), anova_from_moments(&moments).unwrap()); +/// assert!((a.f - b.f).abs() < 1e-12 && (a.p - b.p).abs() < 1e-12); +/// ``` +pub fn anova_from_moments(groups: &[GroupMoments]) -> Option { + let (ss_b, ss_w, ss_t, k, n_total) = one_way_ss_from_moments(groups)?; + anova_from_ss(ss_b, ss_w, ss_t, k, n_total) +} + +/// η² from grouped sufficient statistics — [`eta_squared`]'s value and +/// degeneracy policy (`None` only on fewer than 2 groups, an empty group, or +/// zero total variance; a perfectly separated layout with zero within-group +/// variance is a real `1.0`, unlike the F test). +pub fn eta_squared_from_moments(groups: &[GroupMoments]) -> Option { + let (ss_b, _, ss_t, _, _) = one_way_ss_from_moments(groups)?; + eta_from_ss(ss_b, ss_t) +} + // ── Fisher 2z (D-BLW-5 payload space) ── // // This is the "2z" of plan `cycle-loop-closure-driver-v1.md` §Stage A and the @@ -2246,4 +2362,289 @@ mod tests { ); } } + + // ─────────── one-way ANOVA: mask fold → moments vs materialized ─────────── + // + // The claim under test: a statistic computed from grouped sufficient + // statistics folded straight off a population mask equals the same + // statistic computed from the materialized groups. The fold is + // `ndarray::simd::masked_group_moments_i32` — the kernel + // `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32` delegates to. + // + // Tolerances. Both paths form the same sums of squares by different + // float orderings (two-pass deviations vs exact integers then one + // division). Each is a sum of at most a few thousand terms, so each is + // within ~n·ε ≈ 1e-12 relative of the true value on well-conditioned + // data; F is a ratio of two such sums (≤ 2× that), p goes through the + // regularised incomplete beta whose sensitivity near the observed F is + // bounded here by keeping F moderate. 1e-9 relative on F and η² and + // 1e-9 absolute on p leave three orders of margin and still fail on any + // real algebraic difference (a dropped row, a wrong df, a wrong group). + + mod moments_equivalence { + use super::super::*; + use ndarray::simd::masked_group_moments_i32; + + fn lcg(s: &mut u64) -> u64 { + *s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + *s >> 11 + } + + /// A population: `n` rows, a key lane over `0..k`, an `i32` value + /// lane, and a population mask of the given density (in 1/256ths). + struct Pop { + mask: Vec, + keys: Vec, + values: Vec, + } + + fn population( + n: usize, + k: u32, + density: u64, + seed: u64, + value: impl Fn(&mut u64, u32) -> i32, + ) -> Pop { + let mut s = seed; + let mut mask = vec![0u64; n.div_ceil(64)]; + let mut keys = Vec::with_capacity(n); + let mut values = Vec::with_capacity(n); + for i in 0..n { + // Skewed keys: group g is drawn with weight g+1, so group + // sizes are unequal by construction. + let total = u64::from(k * (k + 1) / 2); + let mut r = lcg(&mut s) % total; + let mut g = 0u32; + while r >= u64::from(g + 1) { + r -= u64::from(g + 1); + g += 1; + } + keys.push(g); + values.push(value(&mut s, g)); + if lcg(&mut s) % 256 < density { + mask[i / 64] |= 1 << (i % 64); + } + } + Pop { mask, keys, values } + } + + /// The conventional path: copy the selected observations out into + /// one vector per group. + fn materialize(p: &Pop, k: u32) -> Vec> { + let mut groups = vec![Vec::new(); k as usize]; + for i in 0..p.values.len() { + if p.mask[i / 64] >> (i % 64) & 1 == 1 { + groups[p.keys[i] as usize].push(f64::from(p.values[i])); + } + } + groups + } + + /// The substrate path: one masked fold, no observation copied out. + fn fold(p: &Pop, k: u32) -> Vec { + let mut out = vec![GroupMoments::EMPTY; k as usize]; + masked_group_moments_i32(&p.mask, &p.keys, &p.values, &mut out); + out + } + + fn rel(a: f64, b: f64) -> f64 { + (a - b).abs() / a.abs().max(b.abs()).max(f64::MIN_POSITIVE) + } + + fn assert_same(slice: Option, moments: Option, what: &str) { + match (slice, moments) { + (None, None) => {} + (Some(a), Some(b)) => { + assert_eq!(a.df_between, b.df_between, "{what}: df_between"); + assert_eq!(a.df_within, b.df_within, "{what}: df_within"); + assert!(rel(a.f, b.f) < 1e-9, "{what}: F {} vs {}", a.f, b.f); + assert!((a.p - b.p).abs() < 1e-9, "{what}: p {} vs {}", a.p, b.p); + assert!( + (a.eta_squared - b.eta_squared).abs() < 1e-9, + "{what}: eta² {} vs {}", + a.eta_squared, + b.eta_squared + ); + } + (a, b) => panic!("{what}: slice {a:?} vs moments {b:?}"), + } + } + + /// FAILS IF: the mask fold + moments projection disagrees with the + /// materialized path anywhere on ordinary data — arbitrary masks + /// (sparse to full), unequal group sizes, negative values, 2..6 + /// groups, several population sizes. + #[test] + fn anova_from_mask_fold_equals_anova_on_materialized_groups() { + let mut checked = 0; + for &n in &[40usize, 130, 1000, 5000] { + for &k in &[2u32, 3, 6] { + for &density in &[20u64, 128, 256] { + let seed = 0xA0A ^ (n as u64) << 8 ^ u64::from(k) << 32 ^ density; + // Group g is centred at 3·g with spread ±40, negatives included. + let p = population(n, k, density, seed, |s, g| { + (lcg(s) % 81) as i32 - 40 + 3 * g as i32 + }); + let groups = materialize(&p, k); + let moments = fold(&p, k); + let what = format!("n={n} k={k} density={density}"); + let a = anova_one_way(&groups); + assert_same(a, anova_from_moments(&moments), &what); + let (ea, eb) = (eta_squared(&groups), eta_squared_from_moments(&moments)); + match (ea, eb) { + (Some(x), Some(y)) => assert!((x - y).abs() < 1e-9, "{what}: η²"), + (x, y) => assert_eq!(x, y, "{what}: η² degeneracy"), + } + checked += usize::from(a.is_some()); + } + } + } + // Anti-vacuity: most configurations must be non-degenerate, or + // the comparison would be `None == None` everywhere. + assert!(checked >= 30, "only {checked} non-degenerate comparisons"); + } + + /// FAILS IF: the moments of the whole population differ from the + /// merged moments of any chunking of it — the fold must be a + /// monoid, so the ANOVA of chunked-then-merged moments is not just + /// close to the one-pass ANOVA but IDENTICAL. + #[test] + fn chunked_folds_merge_to_the_identical_anova() { + let n = 3000; + let k = 4; + let p = population(n, k, 180, 0xC4C4, |s, g| { + (lcg(s) % 1001) as i32 - 500 + 50 * g as i32 + }); + let whole = fold(&p, k); + for chunks in [2usize, 3, 7, 47] { + let mut merged = vec![GroupMoments::EMPTY; k as usize]; + let step = n.div_ceil(chunks); + for c in 0..chunks { + let (lo, hi) = (c * step, ((c + 1) * step).min(n)); + let mut m = vec![0u64; n.div_ceil(64)]; + for i in lo..hi { + m[i / 64] |= p.mask[i / 64] & (1 << (i % 64)); + } + let mut part = vec![GroupMoments::EMPTY; k as usize]; + masked_group_moments_i32(&m, &p.keys, &p.values, &mut part); + for (acc, x) in merged.iter_mut().zip(part) { + *acc = acc.checked_merge(x).expect("in bound"); + } + } + assert_eq!(merged, whole, "chunks={chunks}"); + assert_eq!(anova_from_moments(&merged), anova_from_moments(&whole)); + } + assert!(anova_from_moments(&whole).is_some()); + } + + /// FAILS IF: values near the i32 bound lose precision on the + /// moments path (an i64 square, an f64 accumulation of Σx²). + #[test] + fn values_near_the_i32_bound_agree() { + for (k, base) in [(3u32, i32::MAX - 5000), (3, i32::MIN + 5000)] { + let p = population(2000, k, 200, 0xB0B ^ u64::from(k), move |s, g| { + base + (lcg(s) % 2001) as i32 - 1000 + 400 * g as i32 * base.signum() + }); + let groups = materialize(&p, k); + let a = anova_one_way(&groups); + assert!(a.is_some(), "fixture must be non-degenerate"); + assert_same(a, anova_from_moments(&fold(&p, k)), &format!("base={base}")); + } + } + + /// FAILS IF: any degenerate layout is treated differently from + /// `anova_one_way` — one group, an empty/absent group, `N ≤ k`, or + /// zero within-group variance. + #[test] + fn degenerate_layouts_match_the_slice_contract() { + let m = |n: u64, xs: &[i64]| GroupMoments { + n, + sum: xs.iter().sum(), + sum_sq: xs.iter().map(|&x| i128::from(x * x)).sum(), + }; + let f = |xs: &[i64]| xs.iter().map(|&x| x as f64).collect::>(); + type Case = (&'static str, Vec>, Vec); + let cases: Vec = vec![ + ("one group", vec![f(&[1, 2, 3])], vec![m(3, &[1, 2, 3])]), + ( + "absent group", + vec![f(&[1, 2, 3]), vec![], f(&[4, 5])], + vec![m(3, &[1, 2, 3]), GroupMoments::EMPTY, m(2, &[4, 5])], + ), + ( + "N <= k", + vec![f(&[1]), f(&[2])], + vec![m(1, &[1]), m(1, &[2])], + ), + ( + "zero within-group variance", + vec![f(&[7, 7, 7]), f(&[9, 9])], + vec![m(3, &[7, 7, 7]), m(2, &[9, 9])], + ), + ( + "zero total variance", + vec![f(&[5, 5]), f(&[5, 5, 5])], + vec![m(2, &[5, 5]), m(3, &[5, 5, 5])], + ), + ]; + for (what, groups, moments) in cases { + assert_eq!(anova_one_way(&groups), None, "{what}: slice"); + assert_eq!(anova_from_moments(&moments), None, "{what}: moments"); + assert_eq!( + eta_squared(&groups), + eta_squared_from_moments(&moments), + "{what}: η²" + ); + } + // The one degenerate case where η² is NOT None: perfect separation. + let sep = [m(3, &[7, 7, 7]), m(2, &[9, 9])]; + assert_eq!(eta_squared_from_moments(&sep), Some(1.0)); + } + + /// FAILS IF: the moments path loses the within-group variance to + /// float cancellation. Group means ±2·10⁹ apart with a within-group + /// spread of 1: SS_between ≈ 4.8·10¹⁹, whose f64 ulp (8192) dwarfs + /// SS_within = 12. The expected F is computed independently in + /// exact integers (integer group means, integer grand mean). + #[test] + fn adversarial_large_nearly_equal_values_are_exact() { + let bases = [-2_000_000_000i64, 0, 2_000_000_000]; + let offsets = [-1i64, 0, 1, -1, 0, 1]; + let mut moments = Vec::new(); + let mut groups = Vec::new(); + for b in bases { + let xs: Vec = offsets.iter().map(|o| b + o).collect(); + moments.push(GroupMoments { + n: xs.len() as u64, + sum: xs.iter().sum(), + sum_sq: xs.iter().map(|&x| i128::from(x) * i128::from(x)).sum(), + }); + groups.push(xs.iter().map(|&x| x as f64).collect::>()); + } + // Independent exact oracle: SS_w = Σ offset² per group, SS_b = + // Σ n·(base − 0)², df = (2, 15). + let ss_w: i128 = 3 * offsets.iter().map(|&o| i128::from(o * o)).sum::(); + let ss_b: i128 = bases + .iter() + .map(|&b| 6 * i128::from(b) * i128::from(b)) + .sum(); + let f_exact = (ss_b as f64 / 2.0) / (ss_w as f64 / 15.0); + let got = anova_from_moments(&moments).expect("non-degenerate"); + assert_eq!((got.df_between, got.df_within), (2.0, 15.0)); + assert!( + rel(got.f, f_exact) < 1e-12, + "moments F {} vs exact {f_exact}", + got.f + ); + // Anti-vacuity: the fixture really is adversarial for two-pass + // floats — the slice path either refuses or misses by far more + // than the tolerance used everywhere else. + match anova_one_way(&groups) { + None => {} + Some(a) => assert!(rel(a.f, f_exact) > 1e-6, "fixture not adversarial: {}", a.f), + } + } + } } diff --git a/crates/lance-graph-mask-risc/src/exec.rs b/crates/lance-graph-mask-risc/src/exec.rs index 5ef160695..7a8860d56 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -31,13 +31,14 @@ use ndarray::simd::{ mask_scatter_or_u32, mask_set_range, mask_xor, mask_xor_assign, masked_group_count_u32, masked_group_count_u32_pair, masked_group_count_u32_via, masked_group_max_i32, masked_group_max_i32_pair, masked_group_max_i32_via, masked_group_min_i32, - masked_group_min_i32_pair, masked_group_min_i32_via, masked_group_sum_i32, + masked_group_min_i32_pair, masked_group_min_i32_via, masked_group_moments_i32, + masked_group_moments_i32_pair, masked_group_moments_i32_via, masked_group_sum_i32, masked_group_sum_i32_via, masked_group_sum_sym_i32, masked_group_sum_sym_i32_pair, masked_group_sum_sym_i32_via, masked_key_run_count_u32, masked_max_i32, masked_min_i32, masked_strided_group_sum, masked_sum_i32, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, popcount_batch_u64, ternary_match_strided_to_mask, ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, ternary_match_u64_to_mask, - ternary_match_u64_to_mask_under, KeyRunCarry, + ternary_match_u64_to_mask_under, GroupMoments, KeyRunCarry, }; use crate::ir::{ @@ -1419,6 +1420,7 @@ fn precheck( Terminal::GroupSumI32 { .. } => Some("GroupSumI32"), Terminal::GroupSumViaI32 { .. } => Some("GroupSumViaI32"), Terminal::GroupReduce { .. } => Some("GroupReduce"), + Terminal::GroupMomentsI32 { .. } => Some("GroupMomentsI32"), }; if let Some(what) = refused { return Err(ExecError::ExtentUnsupported { what }); @@ -1572,6 +1574,7 @@ pub fn execute_compiled( } (Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. }, Out::I64(o)) => o.fill(0), (Terminal::GroupReduce { fold, .. }, Out::I64(o)) => o.fill(fold.seed()), + (Terminal::GroupMomentsI32 { .. }, Out::Moments(o)) => o.fill(GroupMoments::EMPTY), _ => {} } let n_rows = planes.n_rows; @@ -1983,6 +1986,36 @@ pub fn execute_compiled( } } } + Terminal::GroupMomentsI32 { mask, key, val } => { + // `validate` already refused a missing/too-small `out`, every + // wrong-width lane and a plane past the 2^32-row exactness + // bound; one delegation per tile (law L3) into the sink + // seeded above with `GroupMoments::EMPTY`. + if let Out::Moments(o) = &mut out { + let m = read(planes, &slots, mask, t); + let v = lane_i32(planes, val, t); + match key { + GroupKey::Lane(k) => { + masked_group_moments_i32(m, lane_u32(planes, k, t), v, o) + } + GroupKey::Via { fk, key } => masked_group_moments_i32_via( + m, + lane_u32(planes, fk, t), + foreign_lane_u32(foreign, key), + v, + o, + ), + GroupKey::Pair { hi, lo, stride } => masked_group_moments_i32_pair( + m, + lane_u32(planes, hi, t), + lane_u32(planes, lo, t), + stride, + v, + o, + ), + } + } + } Terminal::Keep { mask } => { // The demanded mask, one tile at a time. With `Out::None` the // scratch is single-tile (checked above) and the slot IS the @@ -2017,6 +2050,7 @@ pub fn execute_compiled( Terminal::CountKeyRunsU32 { .. } => Value::Count(runs + run_carry.finish()), Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. } => Value::GroupSummed, Terminal::GroupReduce { .. } => Value::GroupReduced, + Terminal::GroupMomentsI32 { .. } => Value::GroupMoments, Terminal::Keep { mask } => Value::Mask(mask), }) } diff --git a/crates/lance-graph-mask-risc/src/ir.rs b/crates/lance-graph-mask-risc/src/ir.rs index 5cf112d6b..9cfdf146a 100644 --- a/crates/lance-graph-mask-risc/src/ir.rs +++ b/crates/lance-graph-mask-risc/src/ir.rs @@ -421,6 +421,29 @@ pub enum Terminal { key: GroupKey, fold: GroupFold, }, + /// Grouped sufficient statistics of an `I32` lane: for every row `i` + /// where `mask` holds, resolves the row's group through [`GroupKey`] and + /// folds `lanes[val][i]` into the caller's `Out::Moments` buffer as + /// `n += 1, Σx += x, Σx² += x²` ([`ndarray::simd::GroupMoments`]). One + /// pass over the selected rows; no selected value is copied out. The + /// buffer's length IS the group universe `K`, and the same zero-fallback + /// drops as [`Terminal::GroupReduce`] apply. + /// + /// A physical fold, not a statistic: consumers project means, variances, + /// F ratios or t statistics from the three sums. It is its own terminal + /// rather than a [`GroupFold`] member because a `GroupFold` slot is one + /// seeded `i64`, and three exact sums are not. + /// + /// Exact under the same row bound as [`Terminal::MaskedSumI32`] + /// ([`MASKED_SUM_I32_MAX_ROWS`], `2^32` rows): the `Σx` field is an `i64`, + /// and no group can exceed the plane. The executor refuses a wider plane + /// rather than wrap. The sink is seeded with `GroupMoments::EMPTY` before + /// the first tile, so an empty group reads `n == 0`. + GroupMomentsI32 { + mask: Operand, + key: GroupKey, + val: u16, + }, } /// Where a [`Terminal::GroupReduce`] reads each row's group. @@ -612,6 +635,7 @@ impl Program { | Terminal::GroupSumI32 { mask, .. } | Terminal::GroupSumViaI32 { mask, .. } | Terminal::GroupReduce { mask, .. } + | Terminal::GroupMomentsI32 { mask, .. } | Terminal::Keep { mask } => touch(mask), } Self { diff --git a/crates/lance-graph-mask-risc/src/lib.rs b/crates/lance-graph-mask-risc/src/lib.rs index 9c58ca92f..4ed2eef45 100644 --- a/crates/lance-graph-mask-risc/src/lib.rs +++ b/crates/lance-graph-mask-risc/src/lib.rs @@ -138,7 +138,7 @@ pub use reference::{ pub use ternlog_dispatch::{ ternlog_any_dispatch, ternlog_dispatch, ternlog_dispatch_assign, ternlog_popcount_dispatch, }; -pub use value::{ExecError, LaneKind, Out, Value}; +pub use value::{ExecError, GroupMoments, LaneKind, Out, Value}; /// Number of `u64` words a mask over `n_rows` occupies. #[inline] diff --git a/crates/lance-graph-mask-risc/src/reference.rs b/crates/lance-graph-mask-risc/src/reference.rs index e1fb7055b..975192612 100644 --- a/crates/lance-graph-mask-risc/src/reference.rs +++ b/crates/lance-graph-mask-risc/src/reference.rs @@ -20,6 +20,7 @@ use crate::ir::{ Foreign, GroupFold, GroupKey, LaneRef, MaskOp, Operand, Planes, Pred, Program, Terminal, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, MAX_SCRATCH_SLOTS, }; +use crate::value::GroupMoments; use crate::value::{ExecError, LaneKind, Out, Value}; use crate::words_for; @@ -34,6 +35,7 @@ pub(crate) enum OutShape { I32(usize), I64(usize), Mask(usize), + Moments(usize), } /// [`OutShape`] of a borrowed `out` — the caller keeps `out` itself to write @@ -44,6 +46,7 @@ pub(crate) fn out_shape(out: &Out<'_>) -> OutShape { Out::I32(v) => OutShape::I32(v.len()), Out::I64(v) => OutShape::I64(v.len()), Out::Mask(v) => OutShape::Mask(v.len()), + Out::Moments(v) => OutShape::Moments(v.len()), } } @@ -516,7 +519,9 @@ pub(crate) fn validate( found: len, }), OutShape::I32(_) => Ok(()), - OutShape::I64(_) | OutShape::Mask(_) => Err(ExecError::BlendNeedsOut), + OutShape::I64(_) | OutShape::Mask(_) | OutShape::Moments(_) => { + Err(ExecError::BlendNeedsOut) + } } } Terminal::ScatterOrU32 { @@ -544,7 +549,7 @@ pub(crate) fn validate( expected: want, found: len, }), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) => { + OutShape::None | OutShape::I32(_) | OutShape::I64(_) | OutShape::Moments(_) => { Err(ExecError::TerminalNeedsOut { what }) } } @@ -559,11 +564,13 @@ pub(crate) fn validate( } match out { OutShape::I64(len) if len >= 1 => Ok(()), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) => { - Err(ExecError::TerminalNeedsOut { - what: "GroupSumI32", - }) - } + OutShape::None + | OutShape::I32(_) + | OutShape::I64(_) + | OutShape::Mask(_) + | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + what: "GroupSumI32", + }), } } Terminal::GroupSumViaI32 { mask, fk, key, val } => { @@ -577,27 +584,19 @@ pub(crate) fn validate( } match out { OutShape::I64(len) if len >= 1 => Ok(()), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) => { - Err(ExecError::TerminalNeedsOut { - what: "GroupSumViaI32", - }) - } + OutShape::None + | OutShape::I32(_) + | OutShape::I64(_) + | OutShape::Mask(_) + | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + what: "GroupSumViaI32", + }), } } Terminal::GroupReduce { mask, key, fold } => { check_operand(p, planes, mask)?; written_slots.readable(mask)?; - match key { - GroupKey::Lane(k) => check_lane(planes, k, LaneKind::U32)?, - GroupKey::Via { fk, key } => { - check_lane(planes, fk, LaneKind::U32)?; - check_foreign_lane(foreign, key, LaneKind::U32)?; - } - GroupKey::Pair { hi, lo, .. } => { - check_lane(planes, hi, LaneKind::U32)?; - check_lane(planes, lo, LaneKind::U32)?; - } - } + check_group_key(planes, foreign, key)?; match fold { GroupFold::Count => {} GroupFold::MinI32(v) | GroupFold::MaxI32(v) => { @@ -612,16 +611,56 @@ pub(crate) fn validate( } match out { OutShape::I64(len) if len >= 1 => Ok(()), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) => { - Err(ExecError::TerminalNeedsOut { - what: "GroupReduce", - }) - } + OutShape::None + | OutShape::I32(_) + | OutShape::I64(_) + | OutShape::Mask(_) + | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + what: "GroupReduce", + }), + } + } + Terminal::GroupMomentsI32 { mask, key, val } => { + check_operand(p, planes, mask)?; + written_slots.readable(mask)?; + check_group_key(planes, foreign, key)?; + check_lane(planes, val, LaneKind::I32)?; + // `Σx` is an `i64`: exact for any group of at most 2^32 rows, and + // no group can hold more rows than the plane. + if n > MASKED_SUM_I32_MAX_ROWS { + return Err(ExecError::SumRowBound { n_rows: n }); + } + match out { + OutShape::Moments(len) if len >= 1 => Ok(()), + _ => Err(ExecError::TerminalNeedsOut { + what: "GroupMomentsI32", + }), } } } } +/// The lane checks every [`GroupKey`] needs: resident and pair keys are +/// `U32` lanes of this table; a VIA key is a `U32` fk here and a `U32` key +/// lane on the foreign table. +fn check_group_key( + planes: &Planes<'_>, + foreign: &Foreign<'_>, + key: GroupKey, +) -> Result<(), ExecError> { + match key { + GroupKey::Lane(k) => check_lane(planes, k, LaneKind::U32), + GroupKey::Via { fk, key } => { + check_lane(planes, fk, LaneKind::U32)?; + check_foreign_lane(foreign, key, LaneKind::U32) + } + GroupKey::Pair { hi, lo, .. } => { + check_lane(planes, hi, LaneKind::U32)?; + check_lane(planes, lo, LaneKind::U32) + } + } +} + fn plane_bit(planes: &Planes<'_>, i: u16, row: usize) -> bool { (planes.masks[usize::from(i)][row / 64] >> (row % 64)) & 1 == 1 } @@ -640,6 +679,41 @@ fn u32_at(planes: &Planes<'_>, lane: u16, row: usize) -> u32 { } } +/// The group of row `r` under `key` in a universe of `groups`, or `None` for +/// every zero-fallback drop: a VIA fk naming no foreign row, a pair minor key +/// at or past its stride, or a resolved key past the universe. Written out +/// longhand for the oracle — never through `ndarray::simd`'s walker. +fn row_group( + planes: &Planes<'_>, + foreign: &Foreign<'_>, + key: GroupKey, + r: usize, + groups: usize, +) -> Option { + let k = match key { + GroupKey::Lane(lane) => u32_at(planes, lane, r) as usize, + GroupKey::Via { fk, key } => { + let remap = match foreign.lanes.get(usize::from(key)) { + Some(LaneRef::U32(v)) => &v[..], + _ => &[][..], + }; + *remap.get(u32_at(planes, fk, r) as usize)? as usize + } + GroupKey::Pair { hi, lo, stride } => { + let lo_v = u32_at(planes, lo, r); + if lo_v >= stride { + return None; + } + // hi and lo are both u32, stride is u32: widen to u64 first so the + // multiply-add cannot overflow, matching ndarray::simd's + // GroupKeyAddr::Pair. + let hi_v = u32_at(planes, hi, r); + usize::try_from(u64::from(hi_v) * u64::from(stride) + u64::from(lo_v)).ok()? + } + }; + (k < groups).then_some(k) +} + /// One predicate, one row — the scalar SPEC of what each `Pred` means. /// /// This is the definition the executor's `ndarray::simd` delegation is diffed @@ -1035,42 +1109,10 @@ pub fn reference_execute_into( for x in o.iter_mut() { *x = seed; } - let remap = match key { - GroupKey::Via { key, .. } => match foreign.lanes.get(usize::from(key)) { - Some(LaneRef::U32(v)) => &v[..], - _ => &[][..], - }, - GroupKey::Lane(_) | GroupKey::Pair { .. } => &[][..], - }; for r in survivors(mask) { - let k = match key { - GroupKey::Lane(lane) => u32_at(planes, lane, r) as usize, - GroupKey::Via { fk, .. } => { - let idx = u32_at(planes, fk, r) as usize; - if idx >= remap.len() { - continue; - } - remap[idx] as usize - } - GroupKey::Pair { hi, lo, stride } => { - let lo_v = u32_at(planes, lo, r); - if lo_v >= stride { - continue; - } - // hi and lo are both u32, stride is u32: widen to - // u64 first so the multiply-add cannot overflow, - // matching ndarray::simd's GroupKeyAddr::Pair. - let hi_v = u32_at(planes, hi, r); - let composite = u64::from(hi_v) * u64::from(stride) + u64::from(lo_v); - match usize::try_from(composite) { - Ok(k) => k, - Err(_) => continue, - } - } - }; - if k >= o.len() { + let Some(k) = row_group(planes, foreign, key, r, o.len()) else { continue; - } + }; o[k] = match fold { GroupFold::Count => o[k] + 1, GroupFold::MinI32(v) => o[k].min(i64::from(i32_at(planes, v, r))), @@ -1089,6 +1131,32 @@ pub fn reference_execute_into( } Value::GroupReduced } + Terminal::GroupMomentsI32 { mask, key, val } => { + if let Out::Moments(o) = out { + // Independent formulation: seed every slot, then walk the + // survivors one row at a time, widening straight to i128 — + // no ndarray kernel and no `GroupMoments::observe` involved. + for x in o.iter_mut() { + *x = GroupMoments::EMPTY; + } + let mut sum = vec![0i128; o.len()]; + for r in survivors(mask) { + let Some(k) = row_group(planes, foreign, key, r, o.len()) else { + continue; + }; + let x = i128::from(i32_at(planes, val, r)); + o[k].n += 1; + sum[k] += x; + o[k].sum_sq += x * x; + } + for (slot, s) in o.iter_mut().zip(sum) { + // In range by the validated row bound; a failure here is + // an oracle bug, never a data condition. + slot.sum = i64::try_from(s).expect("validated 2^32-row bound keeps Σx in i64"); + } + } + Value::GroupMoments + } Terminal::Keep { mask } => { if let Out::Mask(o) = out { for w in o.iter_mut() { diff --git a/crates/lance-graph-mask-risc/src/value.rs b/crates/lance-graph-mask-risc/src/value.rs index 9eeea4bc8..2a9747866 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -5,6 +5,10 @@ //! rejects, or reject with a different reason. use crate::ir::Operand; +/// The per-group `(n, Σx, Σx²)` accumulator [`Out::Moments`] carries. Re-exported +/// here as shared result vocabulary so the independent oracle names the data +/// type through this crate, never through the SIMD facade it falsifies. +pub use ndarray::simd::GroupMoments; /// What a [`crate::Program`] produced. #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -32,6 +36,10 @@ pub enum Value { /// written, one slot per group; see [`crate::GroupFold::seed`] for what /// an empty group holds. GroupReduced, + /// [`crate::Terminal::GroupMomentsI32`]: the caller's `Out::Moments` + /// buffer was written, one [`GroupMoments`] per group; an empty group + /// holds [`GroupMoments::EMPTY`] (`n == 0`). + GroupMoments, /// [`crate::Terminal::MaskedStridedGroupSum`]: the widened sum, or `None` /// when it does not fit an `i64` (never a wrapped value). StridedSum(Option), @@ -52,6 +60,9 @@ pub enum Out<'a> { /// [`crate::Terminal::ScatterOrU32`]'s destination, `words_for(out_rows)` /// long. Mask(&'a mut [u64]), + /// [`crate::Terminal::GroupMomentsI32`]'s destination — one + /// [`GroupMoments`] per group, its length IS the group universe `K`. + Moments(&'a mut [GroupMoments]), } /// The lane width a predicate or terminal expects, for [`ExecError::LaneKind`]. diff --git a/crates/lance-graph-mask-risc/tests/foreign.rs b/crates/lance-graph-mask-risc/tests/foreign.rs index 5774bd965..f70384b91 100644 --- a/crates/lance-graph-mask-risc/tests/foreign.rs +++ b/crates/lance-graph-mask-risc/tests/foreign.rs @@ -2,12 +2,12 @@ //! against the row-at-a-time oracle — the same differential shape //! `tests/differential.rs` uses, extended over a SECOND, foreign row space. -use lance_graph_mask_risc::exec::{execute_into, Scratch}; +use lance_graph_mask_risc::exec::{execute_extent, execute_into, Scratch}; use lance_graph_mask_risc::reference::{reference_execute_into, reference_scratch_with_foreign}; use lance_graph_mask_risc::{ - scratch_words_for, words_for, ExecError, Foreign, ForeignPlane, GroupFold, GroupKey, LaneKind, - LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, Terminal, Value, GROUP_SUM_SYM_MAX_ROWS, - MASKED_SUM_I32_MAX_ROWS, + scratch_words_for, words_for, ExecError, Foreign, ForeignPlane, GroupFold, GroupKey, + GroupMoments, LaneKind, LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, Terminal, Value, + GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, }; fn lcg(seed: &mut u64) -> u64 { @@ -1568,3 +1568,198 @@ fn pair_key_drops_a_minor_key_at_stride() { "rows with lo >= stride must never be counted anywhere" ); } + +/// FAILS IF: `Terminal::GroupMomentsI32` disagrees with the row-at-a-time +/// oracle for any key address (resident, VIA, pair) — including across TILE +/// boundaries, where a sink re-seeded per tile, or a tile that re-counted a +/// row, would change `n`, `Σx` or `Σx²`. +/// +/// Anti-vacuity: the fixture must run more than one tile, drop rows at every +/// key hop, leave some group empty (`n == 0`) and fill some group with more +/// than one row. +#[test] +fn group_moments_match_the_oracle_for_every_key_address() { + let mut multi_tile = false; + let mut empty_group = false; + let mut multi_row_group = false; + for &n in &[1usize, 63, 64, 130, 1000, 5000] { + let fx = Fixture::new(n, 10, 4, 0x3033 ^ n as u64); + let (lanes, masks) = fx.planes(); + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let mut s = 0xA11Cu64 ^ n as u64; + let remap: Vec = (0..fx.foreign_rows) + .map(|_| match lcg(&mut s) % 5 { + 3 => 9, // past every universe below: a second-hop drop + k => k as u32 % 3, + }) + .collect(); + let foreign_lanes = [LaneRef::U32(&remap)]; + let foreign = Foreign { + planes: &[], + lanes: &foreign_lanes, + }; + let words = test_tile_words(n); + multi_tile |= words < words_for(n); + for (key, groups) in [ + (GroupKey::Lane(3), 4usize), + (GroupKey::Via { fk: 0, key: 0 }, 4), + // status (0..3) × key (0..9): minor keys >= stride 2 drop. + ( + GroupKey::Pair { + hi: 3, + lo: 1, + stride: 2, + }, + 12, + ), + ] { + let p = Program::new( + vec![MaskOp::Pred { + pred: Pred::NeU32 { lane: 1, v: 0 }, + under: None, + dst: 0, + }], + Terminal::GroupMomentsI32 { + mask: S0, + key, + val: 2, + }, + ); + let slots = p.scratch_slots as usize; + let mut buf = vec![0u64; scratch_words_for(words, slots).expect("sized")]; + let mut scratch = Scratch::over(&mut buf, words, slots).expect("carves"); + // Dirty sinks on purpose: the terminal must seed them itself. + let dirty = GroupMoments { + n: 7, + sum: -7, + sum_sq: 7, + }; + let mut got_out = vec![dirty; groups]; + let got = execute_into( + &p, + &planes, + &foreign, + &mut scratch, + Out::Moments(&mut got_out), + ) + .expect("runs"); + let mut want_out = vec![dirty; groups]; + let want = reference_execute_into(&p, &planes, &foreign, Out::Moments(&mut want_out)) + .expect("oracle runs"); + assert_eq!(got, Value::GroupMoments, "n={n} {key:?}"); + assert_eq!(got, want, "n={n} {key:?}: value"); + assert_eq!(got_out, want_out, "n={n} {key:?}: sink"); + empty_group |= got_out.iter().any(|g| g.n == 0); + multi_row_group |= got_out.iter().any(|g| g.n > 1); + } + } + assert!(multi_tile, "must exercise more than one tile"); + assert!(empty_group, "must leave some group empty"); + assert!(multi_row_group, "must fold several rows into one group"); +} + +/// FAILS IF: `GroupMomentsI32` accepts a wrong-width value or key lane, a +/// missing or wrong-shaped sink, or runs under an extent (a partial +/// population would silently report partial moments as whole). +#[test] +fn group_moments_refuses_malformed_programs() { + let n = 130; + let fx = Fixture::new(n, 10, 4, 0x5EF); + let (lanes, masks) = fx.planes(); + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let none = Foreign { + planes: &[], + lanes: &[], + }; + let run = |t: Terminal, out: Out<'_>| { + let p = Program::new( + vec![MaskOp::Pred { + pred: Pred::NeU32 { lane: 1, v: 0 }, + under: None, + dst: 0, + }], + t, + ); + reference_execute_into(&p, &planes, &none, out) + }; + let mut sink = vec![GroupMoments::EMPTY; 4]; + // `val` must be an I32 lane (lane 1 is U32). + assert!(matches!( + run( + Terminal::GroupMomentsI32 { + mask: S0, + key: GroupKey::Lane(3), + val: 1 + }, + Out::Moments(&mut sink) + ), + Err(ExecError::LaneKind { .. }) + )); + // The key must be a U32 lane (lane 2 is I32). + assert!(matches!( + run( + Terminal::GroupMomentsI32 { + mask: S0, + key: GroupKey::Lane(2), + val: 2 + }, + Out::Moments(&mut sink) + ), + Err(ExecError::LaneKind { .. }) + )); + // A moments sink is required; an i64 sink is the wrong shape. + let mut i64_sink = vec![0i64; 4]; + for out in [Out::None, Out::I64(&mut i64_sink), Out::Moments(&mut [])] { + assert_eq!( + run( + Terminal::GroupMomentsI32 { + mask: S0, + key: GroupKey::Lane(3), + val: 2 + }, + out + ), + Err(ExecError::TerminalNeedsOut { + what: "GroupMomentsI32" + }) + ); + } + // A partial extent is refused before anything is written: moments over + // part of the population must never be reported as moments of the whole. + let p = Program::new( + vec![], + Terminal::GroupMomentsI32 { + mask: Operand::Plane(0), + key: GroupKey::Lane(3), + val: 2, + }, + ); + let pl = vec![u64::MAX; words_for(n)]; + let masks: [&[u64]; 1] = [&pl]; + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let mut s = Scratch::for_program(&p, n).expect("scratch"); + let mut sink = vec![GroupMoments::EMPTY; 4]; + assert_eq!( + execute_extent(&p, &planes, &none, &mut s, Out::Moments(&mut sink), 10..20), + Err(ExecError::ExtentUnsupported { + what: "GroupMomentsI32" + }) + ); + assert_eq!( + sink, + vec![GroupMoments::EMPTY; 4], + "a refusal writes nothing" + ); +} diff --git a/crates/lance-graph-quack/tests/duckdb_differential.rs b/crates/lance-graph-quack/tests/duckdb_differential.rs index 96da4eadb..95d3b57cc 100644 --- a/crates/lance-graph-quack/tests/duckdb_differential.rs +++ b/crates/lance-graph-quack/tests/duckdb_differential.rs @@ -264,6 +264,7 @@ fn run_query(id: &str, planes: &Planes<'_>, filter: Filter, agg: Agg) -> (String Value::StridedSum(_) => { panic!("case {id}: no MaskedStridedGroupSum case in this suite") } + Value::GroupMoments => panic!("case {id}: no GroupMomentsI32 case in this suite"), }; let out_bytes = encoded.len(); ( From a9a4d1186aae29fadea7aec612c03529522155f2 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 19:45:33 +0000 Subject: [PATCH 02/14] board: grouped-moments fold + ANOVA-over-masks entry Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- ...9-29-grouped-moments-fold-anova-over-masks.md | 16 ++++++++++++++++ .claude/board/entries/README.md | 1 + 2 files changed, 17 insertions(+) create mode 100644 .claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md diff --git a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md new file mode 100644 index 000000000..76e07efaa --- /dev/null +++ b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md @@ -0,0 +1,16 @@ +# 2026-09-29 — One-way ANOVA over a population mask from grouped moments + +**Status:** MEASURED · DONE — ndarray `simd_masking_ops.rs`, `crates/lance-graph-mask-risc`, `crates/jc/src/stats.rs` + +## What landed +- ndarray: `GroupMoments { n: u64, sum: i64, sum_sq: i128 }` and `masked_group_moments_i32{,_via,_pair}` — `(n, Σx, Σx²)` per group in one pass over the rows a mask selects, over the existing `group_walk` (now generic over its slot type). `checked_merge` is exact integer addition; exact up to `GROUP_MOMENTS_MAX_ROWS = 2^32` rows per group (the `i64` sum is the binding field). +- mask-risc: `Terminal::GroupMomentsI32 { mask, key: GroupKey, val }` → `Out::Moments`. Own terminal, not a `GroupFold` member (a `GroupFold` slot is one seeded `i64`). Refuses planes past `MASKED_SUM_I32_MAX_ROWS` and partial extents. +- jc: `anova_from_moments` / `eta_squared_from_moments`. F/p/η² policy is shared with `anova_one_way` / `eta_squared` via `anova_from_ss` / `eta_from_ss`. Sums of squares come from exact `i128` quantities; no two large floats are subtracted. + +## Measured +Mask → fold → `anova_from_moments` vs materialized → `anova_one_way`, 42 non-degenerate layouts (n ≤ 100k, k ≤ 6, densities 20/128/256 of 256): max relative ΔF 4.2e-14, max |Δp| 4.4e-16, max |Δη²| 6.5e-16. Adversarial layout (group means ±2·10⁹, within spread 1): exact F = 3·10¹⁹; the moments path returns 3·10¹⁹, the slice path returns `None` (its `ss_t − ss_b` cancels the within-group SS to ≤ 0). A naive float `Q − S²/N` inside the moments path fails even the near-`i32::MAX` agreement test (ΔF ≈ 1e-3), so the exact forms are load-bearing. + +## Open +- Multi-membership is NOT one pass: `GroupKey` resolves one group per row, and `lance-graph-report`'s `CoordSpec::MaskSet` executes one `Filter::Plane` program per member. Seam: a key address that walks the selected rows once and folds each row into every member plane that holds it. +- Bivariate `(n, Σx, Σy, Σx², Σy², Σxy)`: same walk, a second value lane in the closure, a wider slot type; not built. +- No early exit, and no `CausalEdge64` commit. ANOVA is non-monotone under future rows. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index b0eaf97ad..cbbd0fb1f 100644 --- a/.claude/board/entries/README.md +++ b/.claude/board/entries/README.md @@ -61,6 +61,7 @@ exactly one of them, never both. | 2026-09-30 | `three-reference-sets-are-not-ordinal-aligned` | | [2026-09-30-three-reference-sets-are-not-ordinal-aligned.md](2026-09-30-three-reference-sets-are-not-ordinal-aligned.md) | | 2026-09-30 | `deepnsm-v2-coverage-bands` | | [2026-09-30-deepnsm-v2-coverage-bands.md](2026-09-30-deepnsm-v2-coverage-bands.md) | | 2026-09-30 | `cypher-mask-v2-is-a-replacement-not-a-phase` | | [2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md](2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md) | +| 2026-09-29 | `grouped-moments-fold-anova-over-masks` | | [2026-09-29-grouped-moments-fold-anova-over-masks.md](2026-09-29-grouped-moments-fold-anova-over-masks.md) | | 2026-09-29 | `deepnsm-v2-counted-pick-tag-deltas` | | [2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md](2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md) | | 2026-09-26 | `deepnsm-v2-lexical-evidence-survives-routing` | | [2026-09-26-deepnsm-v2-lexical-evidence-survives-routing.md](2026-09-26-deepnsm-v2-lexical-evidence-survives-routing.md) | | 2026-09-25 | `window-scheduling-and-two-level-ternlog` | | [2026-09-25-window-scheduling-and-two-level-ternlog.md](2026-09-25-window-scheduling-and-two-level-ternlog.md) | From 19e66e3e06717cdf0d562835189fb6fc12d6f971 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 20:20:10 +0000 Subject: [PATCH 03/14] =?UTF-8?q?mask-risc=20+=20jc:=20grouped=20cross=20m?= =?UTF-8?q?oments=20terminal;=20Pearson=20/=20cov=20/=20OLS=20/=20R=C2=B2?= =?UTF-8?q?=20from=20moments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit mask-risc: Terminal::GroupCrossMomentsI32 { mask, key: GroupKey, x, y } folds (n, Σx, Σy, Σx², Σy², Σxy) per group into Out::CrossMoments, one delegation per tile to ndarray::simd::masked_group_cross_moments_i32{,_via, _pair}. Same key addresses, drops, 2^32-row bound and extent refusal as GroupMomentsI32; the independent oracle accumulates all six fields in i128 row by row and names the type only through value.rs (law L4). jc: pearson_from_cross_moments, sample_covariance_from_cross_moments, simple_regression_from_cross_moments (SimpleRegression { slope, intercept }) and r_squared_from_cross_moments. The centred sums n·Σxy − Σx·Σy etc. are formed exactly in i128 (checked; None past the bound). Each projection shares its acceptance policy with the slice implementation it mirrors, extracted as pearson_from_centered (reliability), sample_cov_tail and r_squared_tail (stats) — no second definition of Pearson, the covariance convention (n−1) or the R² slack rule. R² follows multiple_r_squared's one-predictor contract (n ≥ 3, constant x or y → None). Recorded finding: unlike one-way ANOVA, the slice pearson / multiple_r_squared are two-pass and do NOT fail under a huge offset with a tiny spread; a naive f64 projection of the same moments does (NaN), which is what the exact i128 centring prevents. The regression test asserts both halves against an independently computed exact reference. quack's duckdb_differential gains the unreachable arm for the new Value. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/jc/src/reliability.rs | 8 + crates/jc/src/stats.rs | 538 +++++++++++++++++- crates/lance-graph-mask-risc/src/exec.rs | 40 +- crates/lance-graph-mask-risc/src/ir.rs | 17 + crates/lance-graph-mask-risc/src/lib.rs | 2 +- crates/lance-graph-mask-risc/src/reference.rs | 81 ++- crates/lance-graph-mask-risc/src/value.rs | 10 + crates/lance-graph-mask-risc/tests/foreign.rs | 192 ++++++- .../tests/duckdb_differential.rs | 3 + 9 files changed, 859 insertions(+), 32 deletions(-) diff --git a/crates/jc/src/reliability.rs b/crates/jc/src/reliability.rs index c73449ee2..34982499f 100644 --- a/crates/jc/src/reliability.rs +++ b/crates/jc/src/reliability.rs @@ -110,6 +110,14 @@ pub fn pearson(x: &[f64], y: &[f64]) -> Option { sxx += dx * dx; syy += dy * dy; } + pearson_from_centered(sxy, sxx, syy) +} + +/// Pearson's `r` from the centred co-moment and the two centred +/// second moments — the ONE place its degeneracy policy lives, shared by +/// [`pearson`] and `stats::pearson_from_cross_moments`. Any common positive +/// scale on all three (e.g. `n·S` instead of `S`) cancels. +pub(crate) fn pearson_from_centered(sxy: f64, sxx: f64, syy: f64) -> Option { let denom = (sxx * syy).sqrt(); if denom == 0.0 || !denom.is_finite() { // `denom == 0` → at least one series is constant. `denom == ∞` → the diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index 9700becab..bcf1a344c 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -73,8 +73,8 @@ //! [`crate::reliability`] contract exactly. Non-finite input is rejected up //! front by the same `all_finite` guard. -use crate::reliability::{all_finite, mean, pearson}; -use ndarray::simd::GroupMoments; +use crate::reliability::{all_finite, mean, pearson, pearson_from_centered}; +use ndarray::simd::{GroupCrossMoments, GroupMoments}; use std::collections::BTreeSet; // ─────────────────────────── local helpers ─────────────────────────── @@ -104,14 +104,20 @@ fn sample_cov(x: &[f64], y: &[f64]) -> Option { } let mx = mean(x)?; let my = mean(y)?; - let n = x.len() as f64; - Some( - x.iter() - .zip(y.iter()) - .map(|(&a, &b)| (a - mx) * (b - my)) - .sum::() - / (n - 1.0), - ) + let sxy = x + .iter() + .zip(y.iter()) + .map(|(&a, &b)| (a - mx) * (b - my)) + .sum::(); + sample_cov_tail(sxy, x.len()) +} + +/// The unbiased (divisor `n−1`) covariance from the centred co-moment +/// `Σ(x−x̄)(y−ȳ)` — the ONE place the convention lives, shared by +/// `sample_cov` and [`sample_covariance_from_cross_moments`]. +#[inline] +fn sample_cov_tail(sxy: f64, n: usize) -> Option { + (n >= 2).then(|| sxy / (n as f64 - 1.0)) } // ───────────────────── incomplete beta (p-values) ───────────────────── @@ -878,12 +884,17 @@ pub fn multiple_r_squared(y: &[f64], predictors: &[Vec]) -> Option { if !ss_res.is_finite() || ss_res < 0.0 { return None; } - let r2 = 1.0 - ss_res / ss_tot; - // `ss_res` is a sum of squares, so it cannot be negative; the only - // excursions possible are `ss_res` marginally exceeding `ss_tot` (R² - // slightly below 0) or cancellation pushing it a few ulps past 1. Clamp - // ONLY that rounding-scale band — a materially out-of-range value means - // the solve failed and must surface as `None`, not as a plausible 0 or 1. + r_squared_tail(1.0 - ss_res / ss_tot) +} + +/// The acceptance rule for a computed `R²` — the ONE place it lives, shared +/// by [`multiple_r_squared`] and [`r_squared_from_cross_moments`]. +/// +/// The only legitimate excursions are rounding-scale: `ss_res` marginally +/// exceeding `ss_tot` (R² slightly below 0) or cancellation pushing it a few +/// ulps past 1. Clamp ONLY that band — a materially out-of-range value means +/// the solve failed and must surface as `None`, not as a plausible 0 or 1. +fn r_squared_tail(r2: f64) -> Option { const R2_SLACK: f64 = 1e-9; if !(-R2_SLACK..=1.0 + R2_SLACK).contains(&r2) || !r2.is_finite() { return None; @@ -1282,6 +1293,113 @@ pub fn eta_squared_from_moments(groups: &[GroupMoments]) -> Option { eta_from_ss(ss_b, ss_t) } +// ─────────── correlation / covariance / simple OLS from cross moments ─────────── +// +// Pearson, the sample covariance, the simple least-squares line and its R², +// projected from one group's `(n, Σx, Σy, Σx², Σy², Σxy)` — the shape +// `ndarray::simd::masked_group_cross_moments_i32` folds off a population +// mask. Each shares its acceptance policy with the slice implementation it +// mirrors (`pearson_from_centered`, `sample_cov_tail`, `r_squared_tail`); only +// the centred sums are formed differently. + +/// The three centred sums scaled by `n`, formed exactly in `i128`: +/// `(C_xy, C_xx, C_yy) = (n·Σxy − Σx·Σy, n·Σx² − (Σx)², n·Σy² − (Σy)²)`, +/// i.e. `n·Σ(x−x̄)(y−ȳ)` and friends with the means cleared of their +/// denominators. No two large floats are ever subtracted, so a huge common +/// offset with a tiny spread loses nothing. +/// +/// Every product fits: within `GroupCrossMoments`' `2^32`-row bound each of +/// `n·Σx²`, `(Σx)²`, `n·Σxy`, `Σx·Σy` is at most `2^126` in magnitude, and by +/// Cauchy–Schwarz so is each difference. Past the bound the checked +/// arithmetic returns `None` rather than wrap. +fn centered_cross(m: &GroupCrossMoments) -> Option<(i128, i128, i128)> { + let n = i128::from(m.n); + let (sx, sy) = (i128::from(m.sum_x), i128::from(m.sum_y)); + let cxy = n.checked_mul(m.sum_xy)?.checked_sub(sx.checked_mul(sy)?)?; + let cxx = n.checked_mul(m.sum_x2)?.checked_sub(sx.checked_mul(sx)?)?; + let cyy = n.checked_mul(m.sum_y2)?.checked_sub(sy.checked_mul(sy)?)?; + // A negative centred square is impossible for moments of real data. + (cxx >= 0 && cyy >= 0).then_some((cxy, cxx, cyy)) +} + +/// Pearson's `r` from grouped cross moments — [`pearson`]'s value and +/// degeneracy policy (`None` for `n < 2` or a constant `x` or `y`). +/// +/// ``` +/// use jc::stats::pearson_from_cross_moments; +/// use ndarray::simd::GroupCrossMoments; +/// // x = [1,2,3], y = [2,4,6]: perfectly correlated. +/// let m = GroupCrossMoments { n: 3, sum_x: 6, sum_y: 12, sum_x2: 14, sum_y2: 56, sum_xy: 28 }; +/// assert!((pearson_from_cross_moments(&m).unwrap() - 1.0).abs() < 1e-12); +/// ``` +pub fn pearson_from_cross_moments(m: &GroupCrossMoments) -> Option { + if m.n < 2 { + return None; + } + let (cxy, cxx, cyy) = centered_cross(m)?; + pearson_from_centered(cxy as f64, cxx as f64, cyy as f64) +} + +/// The unbiased (divisor `n−1`) sample covariance from grouped cross moments +/// — the convention this module already uses (`None` for `n < 2`). +pub fn sample_covariance_from_cross_moments(m: &GroupCrossMoments) -> Option { + let (cxy, _, _) = centered_cross(m)?; + // `C_xy = n·Σ(x−x̄)(y−ȳ)`: divide the `n` back out before the tail. + let sxy = cxy as f64 / m.n.max(1) as f64; + sample_cov_tail(sxy, usize::try_from(m.n).ok()?) +} + +/// The ordinary least-squares line `y = intercept + slope·x`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct SimpleRegression { + /// `β = C_xy / C_xx`. + pub slope: f64, + /// `α = (Σy − β·Σx) / n`. + pub intercept: f64, +} + +/// The simple least-squares line from grouped cross moments. +/// +/// `None` for `n < 2` or a constant `x` (no line is identified). A constant +/// `y` is a real answer: slope `0`, intercept `ȳ`. +/// +/// **Conditioning of the intercept.** The slope is formed from exact +/// integers (one rounding each for `C_xy` and `C_xx`, then a division), so its +/// relative error is a few ulps whatever the offset. The intercept is an +/// extrapolation to `x = 0`: `α = ȳ − β·x̄` is computed in `f64`, so when +/// `|β·x̄|` is large and `α` is small its ABSOLUTE error is on the order of +/// `ε·(|β|·|x̄| + |ȳ|)`. That is the conditioning of the quantity itself, not +/// of the fold; centre `x` before regressing when `α` itself matters. +pub fn simple_regression_from_cross_moments(m: &GroupCrossMoments) -> Option { + if m.n < 2 { + return None; + } + let (cxy, cxx, _) = centered_cross(m)?; + if cxx == 0 { + return None; // constant x: no slope + } + let slope = cxy as f64 / cxx as f64; + let n = m.n as f64; + let intercept = (m.sum_y as f64 - slope * m.sum_x as f64) / n; + (slope.is_finite() && intercept.is_finite()).then_some(SimpleRegression { slope, intercept }) +} + +/// `R²` of the simple least-squares line from grouped cross moments — +/// [`multiple_r_squared`]'s one-predictor contract: `None` for `n < 3` (one +/// residual degree of freedom beyond the two coefficients), a constant `x` or +/// a constant `y`; otherwise `r²`, through the same acceptance rule. +pub fn r_squared_from_cross_moments(m: &GroupCrossMoments) -> Option { + if m.n < 3 { + return None; + } + let (cxy, cxx, cyy) = centered_cross(m)?; + if cxx == 0 || cyy == 0 { + return None; + } + let r = cxy as f64 / ((cxx as f64).sqrt() * (cyy as f64).sqrt()); + r_squared_tail(r * r) +} + // ── Fisher 2z (D-BLW-5 payload space) ── // // This is the "2z" of plan `cycle-loop-closure-driver-v1.md` §Stage A and the @@ -2647,4 +2765,392 @@ mod tests { } } } + + // ─────── Pearson / covariance / simple OLS: mask fold vs materialized ─────── + // + // Tolerances: both sides form the same centred sums, by two-pass floats + // (slice) or exact integers then one rounding (fold). On the data below + // each is within a few ulps × n of the truth, so 1e-9 relative on r, + // covariance, R² and slope leaves wide margin while still catching any + // real algebraic error. The intercept is compared with the ABSOLUTE + // bound its own conditioning allows (see + // `simple_regression_from_cross_moments`). + + mod cross_moments_equivalence { + use super::super::*; + use crate::reliability::pearson; + use ndarray::simd::masked_group_cross_moments_i32; + + fn lcg(s: &mut u64) -> u64 { + *s = s + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + *s >> 11 + } + + struct Pop { + mask: Vec, + keys: Vec, + xs: Vec, + ys: Vec, + } + + /// Skewed keys (group g drawn with weight g+1 → unequal sizes), an + /// arbitrary mask density, and `(x, y)` from `gen(state, group)`. + fn population( + n: usize, + k: u32, + density: u64, + seed: u64, + gen: impl Fn(&mut u64, u32) -> (i32, i32), + ) -> Pop { + let mut s = seed; + let mut p = Pop { + mask: vec![0u64; n.div_ceil(64)], + keys: Vec::new(), + xs: Vec::new(), + ys: Vec::new(), + }; + let total = u64::from(k * (k + 1) / 2); + for i in 0..n { + let mut r = lcg(&mut s) % total; + let mut g = 0u32; + while r >= u64::from(g + 1) { + r -= u64::from(g + 1); + g += 1; + } + let (x, y) = gen(&mut s, g); + p.keys.push(g); + p.xs.push(x); + p.ys.push(y); + if lcg(&mut s) % 256 < density { + p.mask[i / 64] |= 1 << (i % 64); + } + } + p + } + + fn fold(p: &Pop, k: u32) -> Vec { + let mut out = vec![GroupCrossMoments::EMPTY; k as usize]; + masked_group_cross_moments_i32(&p.mask, &p.keys, &p.xs, &p.ys, &mut out); + out + } + + /// The conventional path: the selected `(x, y)` pairs copied out, per group. + fn materialize(p: &Pop, k: u32) -> Vec<(Vec, Vec)> { + let mut g = vec![(Vec::new(), Vec::new()); k as usize]; + for i in 0..p.xs.len() { + if p.mask[i / 64] >> (i % 64) & 1 == 1 { + let slot = &mut g[p.keys[i] as usize]; + slot.0.push(f64::from(p.xs[i])); + slot.1.push(f64::from(p.ys[i])); + } + } + g + } + + /// The textbook two-pass OLS line on materialized data — an + /// independent formulation (JC has no slice slope/intercept API). + fn ols(x: &[f64], y: &[f64]) -> Option<(f64, f64)> { + if x.len() < 2 { + return None; + } + let (mx, my) = (mean(x)?, mean(y)?); + let sxx: f64 = x.iter().map(|v| (v - mx) * (v - mx)).sum(); + if sxx == 0.0 { + return None; + } + let sxy: f64 = x.iter().zip(y).map(|(a, b)| (a - mx) * (b - my)).sum(); + let slope = sxy / sxx; + Some((slope, my - slope * mx)) + } + + fn rel(a: f64, b: f64) -> f64 { + (a - b).abs() / a.abs().max(b.abs()).max(f64::MIN_POSITIVE) + } + + /// Compare every projection of one group against its slice + /// counterpart; both `None` or both within tolerance. + fn assert_group_agrees(x: &[f64], y: &[f64], m: &GroupCrossMoments, what: &str) { + match (pearson(x, y), pearson_from_cross_moments(m)) { + (Some(a), Some(b)) => assert!(rel(a, b) < 1e-9, "{what}: r {a} vs {b}"), + (a, b) => assert_eq!(a, b, "{what}: r degeneracy"), + } + match (sample_cov(x, y), sample_covariance_from_cross_moments(m)) { + (Some(a), Some(b)) => { + assert!( + (a - b).abs() <= 1e-9 * a.abs().max(1.0), + "{what}: cov {a} vs {b}" + ) + } + (a, b) => assert_eq!(a, b, "{what}: cov degeneracy"), + } + match ( + multiple_r_squared(y, &[x.to_vec()]), + r_squared_from_cross_moments(m), + ) { + (Some(a), Some(b)) => assert!((a - b).abs() < 1e-9, "{what}: R² {a} vs {b}"), + (a, b) => assert_eq!(a, b, "{what}: R² degeneracy"), + } + match (ols(x, y), simple_regression_from_cross_moments(m)) { + (Some((sl, ic)), Some(r)) => { + assert!(rel(sl, r.slope) < 1e-9, "{what}: slope {sl} vs {}", r.slope); + let mx = mean(x).unwrap().abs(); + let my = mean(y).unwrap().abs(); + let bound = 1e-12 * (sl.abs() * mx + my) + 1e-9; + assert!( + (ic - r.intercept).abs() <= bound, + "{what}: intercept {ic} vs {}", + r.intercept + ); + } + (a, b) => assert_eq!(a.is_some(), b.is_some(), "{what}: OLS degeneracy"), + } + } + + /// FAILS IF: any projection from the mask fold disagrees with its + /// materialized counterpart — arbitrary masks, unequal groups, + /// negative values, correlations of every sign and strength. + #[test] + fn cross_projections_from_mask_fold_equal_the_materialized_path() { + let mut compared = 0; + for &n in &[40usize, 130, 1000, 5000] { + for &k in &[1u32, 3, 6] { + for &density in &[20u64, 128, 256] { + for slope in [-3i32, 0, 1, 2] { + let seed = 0xC0 ^ (n as u64) << 8 ^ u64::from(k) << 24 ^ density << 32; + let p = population(n, k, density, seed ^ slope as u64, move |s, g| { + let x = (lcg(s) % 2001) as i32 - 1000 + 7 * g as i32; + let noise = (lcg(s) % 401) as i32 - 200; + (x, slope * x + noise - 50 * g as i32) + }); + let groups = materialize(&p, k); + let moments = fold(&p, k); + for (g, ((x, y), m)) in groups.iter().zip(&moments).enumerate() { + assert_group_agrees( + x, + y, + m, + &format!("n={n} k={k} d={density} slope={slope} g={g}"), + ); + compared += usize::from(x.len() >= 3); + } + } + } + } + } + assert!( + compared > 200, + "only {compared} non-trivial groups compared" + ); + } + + /// FAILS IF: chunked folds, merged in either order, do not reproduce + /// the one-pass moments — and therefore bit-identical projections. + #[test] + fn chunked_cross_folds_merge_to_identical_projections() { + let n = 3000; + let k = 3; + let p = population(n, k, 200, 0xDEC, |s, g| { + let x = (lcg(s) % 1001) as i32 - 500; + (x, 2 * x + (lcg(s) % 101) as i32 - 50 + g as i32) + }); + let whole = fold(&p, k); + for chunks in [2usize, 5, 47] { + let step = n.div_ceil(chunks); + let parts: Vec> = (0..chunks) + .map(|c| { + let mut m = vec![0u64; n.div_ceil(64)]; + for i in c * step..((c + 1) * step).min(n) { + m[i / 64] |= p.mask[i / 64] & (1 << (i % 64)); + } + let mut part = vec![GroupCrossMoments::EMPTY; k as usize]; + masked_group_cross_moments_i32(&m, &p.keys, &p.xs, &p.ys, &mut part); + part + }) + .collect(); + for order in [false, true] { + let mut merged = vec![GroupCrossMoments::EMPTY; k as usize]; + let seq: Vec<&Vec> = if order { + parts.iter().rev().collect() + } else { + parts.iter().collect() + }; + for part in seq { + for (acc, x) in merged.iter_mut().zip(part) { + *acc = acc.checked_merge(*x).expect("in bound"); + } + } + assert_eq!(merged, whole, "chunks={chunks} reversed={order}"); + for (a, b) in merged.iter().zip(&whole) { + assert_eq!(pearson_from_cross_moments(a), pearson_from_cross_moments(b)); + assert_eq!( + simple_regression_from_cross_moments(a), + simple_regression_from_cross_moments(b) + ); + } + } + } + } + + fn moments_of(xs: &[i64], ys: &[i64]) -> GroupCrossMoments { + let mut m = GroupCrossMoments::EMPTY; + for (&x, &y) in xs.iter().zip(ys) { + m.observe(x as i32, y as i32); + } + m + } + + /// FAILS IF: a degenerate or boundary layout is treated differently + /// from the slice contract: N = 0, N = 1, N = 2 (R² needs 3), + /// constant X, constant Y, perfect ±1 correlation. + #[test] + fn degenerate_layouts_match_the_slice_contract() { + let cases: [(&str, &[i64], &[i64]); 7] = [ + ("N=0", &[], &[]), + ("N=1", &[4], &[9]), + ("N=2", &[1, 3], &[2, 7]), + ("constant x", &[5, 5, 5, 5], &[1, 2, 3, 9]), + ("constant y", &[1, 2, 3, 9], &[5, 5, 5, 5]), + ("perfect +1", &[-3, 1, 4, 10], &[-5, 3, 9, 21]), + ("perfect -1", &[-3, 1, 4, 10], &[7, -1, -7, -19]), + ]; + for (what, xs, ys) in cases { + let m = moments_of(xs, ys); + let xf: Vec = xs.iter().map(|&v| v as f64).collect(); + let yf: Vec = ys.iter().map(|&v| v as f64).collect(); + assert_group_agrees(&xf, &yf, &m, what); + } + // The contract, stated positively rather than only as agreement. + let m = |x: &[i64], y: &[i64]| moments_of(x, y); + assert_eq!(pearson_from_cross_moments(&m(&[], &[])), None); + assert_eq!(sample_covariance_from_cross_moments(&m(&[4], &[9])), None); + assert_eq!(r_squared_from_cross_moments(&m(&[1, 3], &[2, 7])), None); + assert!(pearson_from_cross_moments(&m(&[1, 3], &[2, 7])).is_some()); + assert_eq!(pearson_from_cross_moments(&m(&[5, 5, 5], &[1, 2, 3])), None); + assert_eq!( + simple_regression_from_cross_moments(&m(&[5, 5, 5], &[1, 2, 3])), + None + ); + let flat = simple_regression_from_cross_moments(&m(&[1, 2, 3], &[5, 5, 5])).unwrap(); + assert_eq!((flat.slope, flat.intercept), (0.0, 5.0)); + let r = pearson_from_cross_moments(&m(&[-3, 1, 4, 10], &[7, -1, -7, -19])).unwrap(); + assert!((r + 1.0).abs() < 1e-15, "perfect -1: {r}"); + let line = + simple_regression_from_cross_moments(&m(&[-3, 1, 4, 10], &[-5, 3, 9, 21])).unwrap(); + assert_eq!((line.slope, line.intercept), (2.0, 1.0)); + } + + /// FAILS IF: values at the i32 extremes lose precision on the fold + /// path (an i32 product, an f64 Σxy). + #[test] + fn values_at_the_i32_extremes_agree() { + let p = population(2000, 2, 200, 0xE7, |s, g| { + let x = match lcg(s) % 5 { + 0 => i32::MIN, + 1 => i32::MAX, + _ => (lcg(s) % 2_000_001) as i32 - 1_000_000, + }; + let y = if g == 0 { + x / 2 + (lcg(s) % 1000) as i32 + } else { + i32::MAX - (lcg(s) % 7) as i32 + }; + (x, y) + }); + for (g, ((x, y), m)) in materialize(&p, 2).iter().zip(&fold(&p, 2)).enumerate() { + assert!(x.len() > 100, "group {g} too small"); + assert_group_agrees(x, y, m, &format!("extremes g={g}")); + } + } + + /// FAILS IF: the fold path loses a tiny spread under a huge offset. + /// + /// `x = B + d`, `y = C + d + e` with `Σd = Σe = 0` and `|d|, |e| ≤ 4`, + /// for `B ≈ ±2·10⁹`: the exact centred sums are the small integers + /// `Σde`, `Σd²`, … — an INDEPENDENT reference, never computed from the + /// moments. Recorded finding (not forced): the slice `pearson` / + /// `multiple_r_squared` are two-pass (the mean is subtracted before + /// any product), so unlike one-way ANOVA they do NOT fail here — both + /// paths land within an ulp of the exact value. What does fail is a + /// naive f64 projection of the same moments (`n·Σxy − Σx·Σy` in + /// floats): it cancels to garbage. That is the case the exact `i128` + /// centring exists for, and it is asserted as the anti-vacuity half. + #[test] + fn large_offsets_with_tiny_spread_are_exact() { + let d = [-3i64, -1, 0, 1, 3, -2, 2, -4, 4, 0]; + let e = [0i64, 1, -1, 0, 1, -1, 0, 1, -1, 0]; + for (b, c) in [ + (2_000_000_000i64, 1_000_000_000i64), + (-2_000_000_000, -1_900_000_000), + ] { + let xs: Vec = d.iter().map(|v| b + v).collect(); + let ys: Vec = d.iter().zip(e).map(|(v, w)| c + v + w).collect(); + let m = moments_of(&xs, &ys); + let sxy: i64 = d.iter().zip(e).map(|(v, w)| v * (v + w)).sum(); + let sxx: i64 = d.iter().map(|v| v * v).sum(); + let syy: i64 = d.iter().zip(e).map(|(v, w)| (v + w) * (v + w)).sum(); + let r_exact = sxy as f64 / ((sxx * syy) as f64).sqrt(); + let slope_exact = sxy as f64 / sxx as f64; + // α = C − β·B exactly, as the integer fraction (C·Sxx − B·Sxy) / Sxx. + let alpha_exact = (i128::from(c) * i128::from(sxx) + - i128::from(b) * i128::from(sxy)) as f64 + / sxx as f64; + let cov_exact = sxy as f64 / 9.0; + + let r = pearson_from_cross_moments(&m).unwrap(); + assert!(rel(r, r_exact) < 1e-15, "B={b}: r {r} vs exact {r_exact}"); + let r2 = r_squared_from_cross_moments(&m).unwrap(); + assert!(rel(r2, r_exact * r_exact) < 1e-15, "B={b}: R²"); + let cov = sample_covariance_from_cross_moments(&m).unwrap(); + assert!( + rel(cov, cov_exact) < 1e-15, + "B={b}: cov {cov} vs {cov_exact}" + ); + let line = simple_regression_from_cross_moments(&m).unwrap(); + assert!(rel(line.slope, slope_exact) < 1e-15, "B={b}: slope"); + // The documented intercept bound: ε·(|β|·|x̄| + |ȳ|). + let bound = + f64::EPSILON * 4.0 * (slope_exact.abs() * b.abs() as f64 + c.abs() as f64); + assert!( + (line.intercept - alpha_exact).abs() <= bound, + "B={b}: α {} vs {alpha_exact}", + line.intercept + ); + + // The slice path is two-pass and survives (recorded, not forced). + let xf: Vec = xs.iter().map(|&v| v as f64).collect(); + let yf: Vec = ys.iter().map(|&v| v as f64).collect(); + assert!(rel(pearson(&xf, &yf).unwrap(), r_exact) < 1e-12); + + // Anti-vacuity: the naive float projection of the SAME moments fails. + let n = m.n as f64; + let naive = (n * m.sum_xy as f64 - m.sum_x as f64 * m.sum_y as f64) + / ((n * m.sum_x2 as f64 - (m.sum_x as f64).powi(2)).sqrt() + * (n * m.sum_y2 as f64 - (m.sum_y as f64).powi(2)).sqrt()); + assert!( + !naive.is_finite() || rel(naive, r_exact) > 1e-3, + "fixture not adversarial: {naive}" + ); + } + } + + /// FAILS IF: moments past the exactness bound are projected instead + /// of refused — `centered_cross`'s checked products must say `None`. + #[test] + fn moments_past_the_bound_are_refused() { + let huge = GroupCrossMoments { + n: u64::MAX, + sum_x: 1, + sum_y: 1, + sum_x2: i128::MAX / 2, + sum_y2: i128::MAX / 2, + sum_xy: i128::MAX / 2, + }; + assert_eq!(pearson_from_cross_moments(&huge), None); + assert_eq!(sample_covariance_from_cross_moments(&huge), None); + assert_eq!(simple_regression_from_cross_moments(&huge), None); + assert_eq!(r_squared_from_cross_moments(&huge), None); + } + } } diff --git a/crates/lance-graph-mask-risc/src/exec.rs b/crates/lance-graph-mask-risc/src/exec.rs index 7a8860d56..72cc41f88 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -29,7 +29,8 @@ use ndarray::simd::{ lt_i32_to_mask_under, mask_all, mask_and, mask_and_assign, mask_andnot, mask_andnot_assign, mask_any, mask_gather_u32, mask_not, mask_not_assign, mask_or, mask_or_assign, mask_scatter_or_u32, mask_set_range, mask_xor, mask_xor_assign, masked_group_count_u32, - masked_group_count_u32_pair, masked_group_count_u32_via, masked_group_max_i32, + masked_group_count_u32_pair, masked_group_count_u32_via, masked_group_cross_moments_i32, + masked_group_cross_moments_i32_pair, masked_group_cross_moments_i32_via, masked_group_max_i32, masked_group_max_i32_pair, masked_group_max_i32_via, masked_group_min_i32, masked_group_min_i32_pair, masked_group_min_i32_via, masked_group_moments_i32, masked_group_moments_i32_pair, masked_group_moments_i32_via, masked_group_sum_i32, @@ -38,7 +39,7 @@ use ndarray::simd::{ masked_strided_group_sum, masked_sum_i32, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, popcount_batch_u64, ternary_match_strided_to_mask, ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, ternary_match_u64_to_mask, - ternary_match_u64_to_mask_under, GroupMoments, KeyRunCarry, + ternary_match_u64_to_mask_under, GroupCrossMoments, GroupMoments, KeyRunCarry, }; use crate::ir::{ @@ -1421,6 +1422,7 @@ fn precheck( Terminal::GroupSumViaI32 { .. } => Some("GroupSumViaI32"), Terminal::GroupReduce { .. } => Some("GroupReduce"), Terminal::GroupMomentsI32 { .. } => Some("GroupMomentsI32"), + Terminal::GroupCrossMomentsI32 { .. } => Some("GroupCrossMomentsI32"), }; if let Some(what) = refused { return Err(ExecError::ExtentUnsupported { what }); @@ -1575,6 +1577,9 @@ pub fn execute_compiled( (Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. }, Out::I64(o)) => o.fill(0), (Terminal::GroupReduce { fold, .. }, Out::I64(o)) => o.fill(fold.seed()), (Terminal::GroupMomentsI32 { .. }, Out::Moments(o)) => o.fill(GroupMoments::EMPTY), + (Terminal::GroupCrossMomentsI32 { .. }, Out::CrossMoments(o)) => { + o.fill(GroupCrossMoments::EMPTY) + } _ => {} } let n_rows = planes.n_rows; @@ -2016,6 +2021,36 @@ pub fn execute_compiled( } } } + Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { + // Same contract as GroupMomentsI32, both lanes read in place: + // one delegation per tile into the sink seeded above. + if let Out::CrossMoments(o) = &mut out { + let m = read(planes, &slots, mask, t); + let (xs, ys) = (lane_i32(planes, x, t), lane_i32(planes, y, t)); + match key { + GroupKey::Lane(k) => { + masked_group_cross_moments_i32(m, lane_u32(planes, k, t), xs, ys, o) + } + GroupKey::Via { fk, key } => masked_group_cross_moments_i32_via( + m, + lane_u32(planes, fk, t), + foreign_lane_u32(foreign, key), + xs, + ys, + o, + ), + GroupKey::Pair { hi, lo, stride } => masked_group_cross_moments_i32_pair( + m, + lane_u32(planes, hi, t), + lane_u32(planes, lo, t), + stride, + xs, + ys, + o, + ), + } + } + } Terminal::Keep { mask } => { // The demanded mask, one tile at a time. With `Out::None` the // scratch is single-tile (checked above) and the slot IS the @@ -2051,6 +2086,7 @@ pub fn execute_compiled( Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. } => Value::GroupSummed, Terminal::GroupReduce { .. } => Value::GroupReduced, Terminal::GroupMomentsI32 { .. } => Value::GroupMoments, + Terminal::GroupCrossMomentsI32 { .. } => Value::GroupCrossMoments, Terminal::Keep { mask } => Value::Mask(mask), }) } diff --git a/crates/lance-graph-mask-risc/src/ir.rs b/crates/lance-graph-mask-risc/src/ir.rs index 9cfdf146a..d22ce5e2d 100644 --- a/crates/lance-graph-mask-risc/src/ir.rs +++ b/crates/lance-graph-mask-risc/src/ir.rs @@ -444,6 +444,22 @@ pub enum Terminal { key: GroupKey, val: u16, }, + /// Grouped cross moments of TWO `I32` lanes: for every row `i` where + /// `mask` holds, folds `(lanes[x][i], lanes[y][i])` into the caller's + /// `Out::CrossMoments` buffer as `n, Σx, Σy, Σx², Σy², Σxy` + /// ([`ndarray::simd::GroupCrossMoments`]). The bivariate member of the + /// [`Terminal::GroupMomentsI32`] family: same key addresses, same drops, + /// same one pass with both lanes read in place, same `2^32`-row + /// exactness bound ([`MASKED_SUM_I32_MAX_ROWS`] — the `Σx`/`Σy` fields + /// are `i64`). A physical fold: covariance, correlation and simple + /// regression are consumer projections. `x == y` is legal (it folds the + /// univariate moments twice over). + GroupCrossMomentsI32 { + mask: Operand, + key: GroupKey, + x: u16, + y: u16, + }, } /// Where a [`Terminal::GroupReduce`] reads each row's group. @@ -636,6 +652,7 @@ impl Program { | Terminal::GroupSumViaI32 { mask, .. } | Terminal::GroupReduce { mask, .. } | Terminal::GroupMomentsI32 { mask, .. } + | Terminal::GroupCrossMomentsI32 { mask, .. } | Terminal::Keep { mask } => touch(mask), } Self { diff --git a/crates/lance-graph-mask-risc/src/lib.rs b/crates/lance-graph-mask-risc/src/lib.rs index 4ed2eef45..ba9454285 100644 --- a/crates/lance-graph-mask-risc/src/lib.rs +++ b/crates/lance-graph-mask-risc/src/lib.rs @@ -138,7 +138,7 @@ pub use reference::{ pub use ternlog_dispatch::{ ternlog_any_dispatch, ternlog_dispatch, ternlog_dispatch_assign, ternlog_popcount_dispatch, }; -pub use value::{ExecError, GroupMoments, LaneKind, Out, Value}; +pub use value::{ExecError, GroupCrossMoments, GroupMoments, LaneKind, Out, Value}; /// Number of `u64` words a mask over `n_rows` occupies. #[inline] diff --git a/crates/lance-graph-mask-risc/src/reference.rs b/crates/lance-graph-mask-risc/src/reference.rs index 975192612..b8652760c 100644 --- a/crates/lance-graph-mask-risc/src/reference.rs +++ b/crates/lance-graph-mask-risc/src/reference.rs @@ -20,8 +20,8 @@ use crate::ir::{ Foreign, GroupFold, GroupKey, LaneRef, MaskOp, Operand, Planes, Pred, Program, Terminal, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, MAX_SCRATCH_SLOTS, }; -use crate::value::GroupMoments; use crate::value::{ExecError, LaneKind, Out, Value}; +use crate::value::{GroupCrossMoments, GroupMoments}; use crate::words_for; /// The caller's terminal-out, described by SHAPE rather than borrowed — what @@ -36,6 +36,7 @@ pub(crate) enum OutShape { I64(usize), Mask(usize), Moments(usize), + CrossMoments(usize), } /// [`OutShape`] of a borrowed `out` — the caller keeps `out` itself to write @@ -47,6 +48,7 @@ pub(crate) fn out_shape(out: &Out<'_>) -> OutShape { Out::I64(v) => OutShape::I64(v.len()), Out::Mask(v) => OutShape::Mask(v.len()), Out::Moments(v) => OutShape::Moments(v.len()), + Out::CrossMoments(v) => OutShape::CrossMoments(v.len()), } } @@ -519,9 +521,10 @@ pub(crate) fn validate( found: len, }), OutShape::I32(_) => Ok(()), - OutShape::I64(_) | OutShape::Mask(_) | OutShape::Moments(_) => { - Err(ExecError::BlendNeedsOut) - } + OutShape::I64(_) + | OutShape::Mask(_) + | OutShape::Moments(_) + | OutShape::CrossMoments(_) => Err(ExecError::BlendNeedsOut), } } Terminal::ScatterOrU32 { @@ -549,9 +552,11 @@ pub(crate) fn validate( expected: want, found: len, }), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) | OutShape::Moments(_) => { - Err(ExecError::TerminalNeedsOut { what }) - } + OutShape::None + | OutShape::I32(_) + | OutShape::I64(_) + | OutShape::Moments(_) + | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { what }), } } Terminal::GroupSumI32 { mask, key, val } => { @@ -568,7 +573,8 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::Moments(_) + | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { what: "GroupSumI32", }), } @@ -588,7 +594,8 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::Moments(_) + | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { what: "GroupSumViaI32", }), } @@ -615,7 +622,8 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::Moments(_) + | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { what: "GroupReduce", }), } @@ -637,6 +645,24 @@ pub(crate) fn validate( }), } } + Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { + check_operand(p, planes, mask)?; + written_slots.readable(mask)?; + check_group_key(planes, foreign, key)?; + check_lane(planes, x, LaneKind::I32)?; + check_lane(planes, y, LaneKind::I32)?; + // `Σx` and `Σy` are `i64`: exact for any group of at most 2^32 + // rows. The i128 fields cannot overflow at any row count. + if n > MASKED_SUM_I32_MAX_ROWS { + return Err(ExecError::SumRowBound { n_rows: n }); + } + match out { + OutShape::CrossMoments(len) if len >= 1 => Ok(()), + _ => Err(ExecError::TerminalNeedsOut { + what: "GroupCrossMomentsI32", + }), + } + } } } @@ -1157,6 +1183,41 @@ pub fn reference_execute_into( } Value::GroupMoments } + Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { + if let Out::CrossMoments(o) = out { + // Independent formulation, as for GroupMomentsI32: every + // field accumulated in i128 row by row, no kernel involved. + let mut acc = vec![[0i128; 6]; o.len()]; + for r in survivors(mask) { + let Some(k) = row_group(planes, foreign, key, r, o.len()) else { + continue; + }; + let xv = i128::from(i32_at(planes, x, r)); + let yv = i128::from(i32_at(planes, y, r)); + let a = &mut acc[k]; + a[0] += 1; + a[1] += xv; + a[2] += yv; + a[3] += xv * xv; + a[4] += yv * yv; + a[5] += xv * yv; + } + for (slot, a) in o.iter_mut().zip(acc) { + let narrow = |v: i128| { + i64::try_from(v).expect("validated 2^32-row bound keeps a sum in i64") + }; + *slot = GroupCrossMoments { + n: u64::try_from(a[0]).expect("a count is non-negative"), + sum_x: narrow(a[1]), + sum_y: narrow(a[2]), + sum_x2: a[3], + sum_y2: a[4], + sum_xy: a[5], + }; + } + } + Value::GroupCrossMoments + } Terminal::Keep { mask } => { if let Out::Mask(o) = out { for w in o.iter_mut() { diff --git a/crates/lance-graph-mask-risc/src/value.rs b/crates/lance-graph-mask-risc/src/value.rs index 2a9747866..84b603012 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -5,6 +5,9 @@ //! rejects, or reject with a different reason. use crate::ir::Operand; +/// The per-group `(n, Σx, Σy, Σx², Σy², Σxy)` accumulator +/// [`Out::CrossMoments`] carries — re-exported for the same reason. +pub use ndarray::simd::GroupCrossMoments; /// The per-group `(n, Σx, Σx²)` accumulator [`Out::Moments`] carries. Re-exported /// here as shared result vocabulary so the independent oracle names the data /// type through this crate, never through the SIMD facade it falsifies. @@ -40,6 +43,10 @@ pub enum Value { /// buffer was written, one [`GroupMoments`] per group; an empty group /// holds [`GroupMoments::EMPTY`] (`n == 0`). GroupMoments, + /// [`crate::Terminal::GroupCrossMomentsI32`]: the caller's + /// `Out::CrossMoments` buffer was written, one [`GroupCrossMoments`] per + /// group; an empty group holds [`GroupCrossMoments::EMPTY`]. + GroupCrossMoments, /// [`crate::Terminal::MaskedStridedGroupSum`]: the widened sum, or `None` /// when it does not fit an `i64` (never a wrapped value). StridedSum(Option), @@ -63,6 +70,9 @@ pub enum Out<'a> { /// [`crate::Terminal::GroupMomentsI32`]'s destination — one /// [`GroupMoments`] per group, its length IS the group universe `K`. Moments(&'a mut [GroupMoments]), + /// [`crate::Terminal::GroupCrossMomentsI32`]'s destination — one + /// [`GroupCrossMoments`] per group, its length IS the group universe `K`. + CrossMoments(&'a mut [GroupCrossMoments]), } /// The lane width a predicate or terminal expects, for [`ExecError::LaneKind`]. diff --git a/crates/lance-graph-mask-risc/tests/foreign.rs b/crates/lance-graph-mask-risc/tests/foreign.rs index f70384b91..2ebbb07fa 100644 --- a/crates/lance-graph-mask-risc/tests/foreign.rs +++ b/crates/lance-graph-mask-risc/tests/foreign.rs @@ -5,9 +5,9 @@ use lance_graph_mask_risc::exec::{execute_extent, execute_into, Scratch}; use lance_graph_mask_risc::reference::{reference_execute_into, reference_scratch_with_foreign}; use lance_graph_mask_risc::{ - scratch_words_for, words_for, ExecError, Foreign, ForeignPlane, GroupFold, GroupKey, - GroupMoments, LaneKind, LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, Terminal, Value, - GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, + scratch_words_for, words_for, ExecError, Foreign, ForeignPlane, GroupCrossMoments, GroupFold, + GroupKey, GroupMoments, LaneKind, LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, + Terminal, Value, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, }; fn lcg(seed: &mut u64) -> u64 { @@ -1763,3 +1763,189 @@ fn group_moments_refuses_malformed_programs() { "a refusal writes nothing" ); } + +/// A second `I32` lane for the cross terminal (lane 4), independent of +/// `amount`, with `i32::MIN`/`i32::MAX` planted so extreme products occur. +fn y_lane(n: usize, seed: u64) -> Vec { + let mut s = seed; + (0..n) + .map(|i| match i % 11 { + 3 => i32::MIN, + 4 => i32::MAX, + _ => (lcg(&mut s) % 6000) as i32 - 3000, + }) + .collect() +} + +/// FAILS IF: `Terminal::GroupCrossMomentsI32` disagrees with the +/// row-at-a-time oracle for any key address, across tile boundaries, or +/// pairs `x[i]` with a `y` from another row (the oracle reads both at `r`). +/// Anti-vacuity: multi-tile, an empty group, a multi-row group. +#[test] +fn group_cross_moments_match_the_oracle_for_every_key_address() { + let mut multi_tile = false; + let mut empty_group = false; + let mut multi_row_group = false; + for &n in &[1usize, 63, 64, 130, 1000, 5000] { + let fx = Fixture::new(n, 10, 4, 0xC055 ^ n as u64); + let ys = y_lane(n, 0x7777 ^ n as u64); + let (mut lanes, masks) = fx.planes(); + lanes.push(LaneRef::I32(&ys)); + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let mut s = 0xA11Du64 ^ n as u64; + let remap: Vec = (0..fx.foreign_rows) + .map(|_| match lcg(&mut s) % 5 { + 3 => 9, + k => k as u32 % 3, + }) + .collect(); + let foreign_lanes = [LaneRef::U32(&remap)]; + let foreign = Foreign { + planes: &[], + lanes: &foreign_lanes, + }; + let words = test_tile_words(n); + multi_tile |= words < words_for(n); + for (key, groups) in [ + (GroupKey::Lane(3), 4usize), + (GroupKey::Via { fk: 0, key: 0 }, 4), + ( + GroupKey::Pair { + hi: 3, + lo: 1, + stride: 2, + }, + 12, + ), + ] { + let p = Program::new( + vec![MaskOp::Pred { + pred: Pred::NeU32 { lane: 1, v: 0 }, + under: None, + dst: 0, + }], + Terminal::GroupCrossMomentsI32 { + mask: S0, + key, + x: 2, + y: 4, + }, + ); + let slots = p.scratch_slots as usize; + let mut buf = vec![0u64; scratch_words_for(words, slots).expect("sized")]; + let mut scratch = Scratch::over(&mut buf, words, slots).expect("carves"); + let dirty = GroupCrossMoments { + n: 3, + sum_xy: -3, + ..GroupCrossMoments::EMPTY + }; + let mut got_out = vec![dirty; groups]; + let got = execute_into( + &p, + &planes, + &foreign, + &mut scratch, + Out::CrossMoments(&mut got_out), + ) + .expect("runs"); + let mut want_out = vec![dirty; groups]; + let want = + reference_execute_into(&p, &planes, &foreign, Out::CrossMoments(&mut want_out)) + .expect("oracle runs"); + assert_eq!(got, Value::GroupCrossMoments, "n={n} {key:?}"); + assert_eq!(got, want, "n={n} {key:?}: value"); + assert_eq!(got_out, want_out, "n={n} {key:?}: sink"); + empty_group |= got_out.iter().any(|g| g.n == 0); + multi_row_group |= got_out.iter().any(|g| g.n > 1); + } + } + assert!(multi_tile && empty_group && multi_row_group); +} + +/// FAILS IF: the cross terminal accepts a wrong-width `x` or `y` lane, a +/// missing or wrong-shaped sink (an `Out::Moments` is the WRONG shape), or +/// runs under a partial extent. +#[test] +fn group_cross_moments_refuses_malformed_programs() { + let n = 130; + let fx = Fixture::new(n, 10, 4, 0x5F0); + let ys = y_lane(n, 1); + let (mut lanes, _) = fx.planes(); + lanes.push(LaneRef::I32(&ys)); + // Every row selected, tail bits past `n` clear (a dirty tail is refused + // as `PlaneTail` before any terminal check runs). + let mut pl = vec![u64::MAX; words_for(n)]; + pl[words_for(n) - 1] = (1u64 << (n % 64)) - 1; + let masks: [&[u64]; 1] = [&pl]; + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let none = Foreign { + planes: &[], + lanes: &[], + }; + let term = |x: u16, y: u16| Terminal::GroupCrossMomentsI32 { + mask: Operand::Plane(0), + key: GroupKey::Lane(3), + x, + y, + }; + let run = |t: Terminal, out: Out<'_>| { + reference_execute_into(&Program::new(vec![], t), &planes, &none, out) + }; + let mut sink = vec![GroupCrossMoments::EMPTY; 4]; + for (x, y) in [(1u16, 4u16), (2, 1)] { + assert!( + matches!( + run(term(x, y), Out::CrossMoments(&mut sink)), + Err(ExecError::LaneKind { .. }) + ), + "x={x} y={y}" + ); + } + let mut univariate = vec![GroupMoments::EMPTY; 4]; + for out in [ + Out::None, + Out::Moments(&mut univariate), + Out::CrossMoments(&mut []), + ] { + assert_eq!( + run(term(2, 4), out), + Err(ExecError::TerminalNeedsOut { + what: "GroupCrossMomentsI32" + }) + ); + } + let p = Program::new(vec![], term(2, 4)); + let mut s = Scratch::for_program(&p, n).expect("scratch"); + assert_eq!( + execute_extent( + &p, + &planes, + &none, + &mut s, + Out::CrossMoments(&mut sink), + 10..20 + ), + Err(ExecError::ExtentUnsupported { + what: "GroupCrossMomentsI32" + }) + ); + assert_eq!( + sink, + vec![GroupCrossMoments::EMPTY; 4], + "a refusal writes nothing" + ); + // x == y is legal and folds the univariate moments. + let mut same = vec![GroupCrossMoments::EMPTY; 4]; + run(term(2, 2), Out::CrossMoments(&mut same)).expect("x == y is legal"); + assert!(same + .iter() + .all(|g| g.sum_x == g.sum_y && g.sum_xy == g.sum_x2)); +} diff --git a/crates/lance-graph-quack/tests/duckdb_differential.rs b/crates/lance-graph-quack/tests/duckdb_differential.rs index 95d3b57cc..8113fb011 100644 --- a/crates/lance-graph-quack/tests/duckdb_differential.rs +++ b/crates/lance-graph-quack/tests/duckdb_differential.rs @@ -265,6 +265,9 @@ fn run_query(id: &str, planes: &Planes<'_>, filter: Filter, agg: Agg) -> (String panic!("case {id}: no MaskedStridedGroupSum case in this suite") } Value::GroupMoments => panic!("case {id}: no GroupMomentsI32 case in this suite"), + Value::GroupCrossMoments => { + panic!("case {id}: no GroupCrossMomentsI32 case in this suite") + } }; let out_bytes = encoded.len(); ( From 8101f91e4bb085dc7d82a3a8ab123f65f88bdc60 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 20:21:38 +0000 Subject: [PATCH 04/14] jc: make the past-the-bound refusal test load-bearing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The first fixture (i128::MAX/2 sums) wrapped to negative centred squares, so it was refused by the >= 0 check and passed with every checked operation replaced by wrapping arithmetic. The new fixture wraps to a plausible positive value (n·Σx² = 2^128 + 2^33 -> 2^33), so a wrapping implementation would report a confident r = 1; verified by a disable run. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/jc/src/stats.rs | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index bcf1a344c..e96daf5d0 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -3137,15 +3137,25 @@ mod tests { /// FAILS IF: moments past the exactness bound are projected instead /// of refused — `centered_cross`'s checked products must say `None`. + /// + /// The fixture is chosen so that WRAPPING arithmetic would produce a + /// plausible answer rather than an obviously broken one: + /// `n = 2^33`, `Σx² = Σy² = Σxy = 2^95 + 1`, so `n·Σx² = 2^128 + 2^33` + /// wraps to the small positive `2^33` and a wrapping implementation + /// reports a confident `r = 1`. (A first version used `i128::MAX / 2`, + /// which wraps to NEGATIVE centred squares — refused by the `≥ 0` + /// check for the wrong reason, so the test passed with every checked + /// operation disabled. Measured by a disable run.) #[test] fn moments_past_the_bound_are_refused() { + let v = (1i128 << 95) + 1; let huge = GroupCrossMoments { - n: u64::MAX, - sum_x: 1, - sum_y: 1, - sum_x2: i128::MAX / 2, - sum_y2: i128::MAX / 2, - sum_xy: i128::MAX / 2, + n: 1 << 33, + sum_x: 0, + sum_y: 0, + sum_x2: v, + sum_y2: v, + sum_xy: v, }; assert_eq!(pearson_from_cross_moments(&huge), None); assert_eq!(sample_covariance_from_cross_moments(&huge), None); From 9b0030e4d08b8ecf7cccb9c50178097897572303 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 20:23:19 +0000 Subject: [PATCH 05/14] board: grouped cross moments entry Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- .../2026-09-29-grouped-cross-moments-fold.md | 17 +++++++++++++++++ .claude/board/entries/README.md | 1 + 2 files changed, 18 insertions(+) create mode 100644 .claude/board/entries/2026-09-29-grouped-cross-moments-fold.md diff --git a/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md new file mode 100644 index 000000000..c64b1ab6b --- /dev/null +++ b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md @@ -0,0 +1,17 @@ +# 2026-09-29 — Grouped cross moments: Pearson / covariance / OLS / R² over masks + +**Status:** MEASURED · DONE — ndarray `simd_masking_ops.rs`, `crates/lance-graph-mask-risc`, `crates/jc` + +## What landed +- ndarray: `GroupCrossMoments { n, sum_x, sum_y: i64, sum_x2, sum_y2, sum_xy: i128 }` + `masked_group_cross_moments_i32{,_via,_pair}` over the same generic `group_walk`. `x_moments()`/`y_moments()` return the exact univariate `GroupMoments`; a test pins them equal to the univariate fold. Bounds: squares ∈ [0, 2^62], products ∈ [−2^62+2^31, 2^62]; the `i64` sums bind at 2^32 rows/group; the `i128` sums cannot overflow at any `u64` row count. +- mask-risc: `Terminal::GroupCrossMomentsI32 { mask, key, x, y }` → `Out::CrossMoments`. +- jc: `pearson_from_cross_moments`, `sample_covariance_from_cross_moments` (n−1), `simple_regression_from_cross_moments` (`SimpleRegression { slope, intercept }`), `r_squared_from_cross_moments` (`multiple_r_squared`'s k=1 contract). Shared tails extracted: `pearson_from_centered`, `sample_cov_tail`, `r_squared_tail`. Centred sums formed exactly in `i128`, checked. + +## Measured +48 layouts (n ≤ 100k): |Δr| ≤ 1.3e-14, |ΔR²| ≤ 4.1e-14, relative Δslope ≤ 2.8e-14 vs the materialized path. Huge offset + tiny spread (x ≈ ±2·10⁹): fold path within 1e-15 of an independent exact reference. **Unlike ANOVA, the slice path does not fail here** — `pearson`/`multiple_r_squared` are two-pass. A naive f64 projection of the same moments returns NaN; that is what exact `i128` centring prevents (disable run: 2 tests fail). + +## Redundant guards (disable runs, kept as explicit contract) +`pearson_from_cross_moments`'s `n < 2`, the regression's `cxx == 0`, and R²'s `cyy == 0` are each subsumed by the shared tail's zero/non-finite rejection. The first past-bound refusal fixture was vacuous (wrapped to negative centred squares); replaced with one that wraps to a plausible `r = 1`. + +## Open +- Multi-membership in one pass: `lance-graph-report` lowers `CoordSpec::MaskSet` to one `Filter::Plane` pass per member tuple and per fold state. The single-group limit is in `GroupKeyAddr::group_of` (one `Option` per row), in the IR's `GroupKey`, and in report lowering. Seam: a word-level walker `for each 64-row word: sel & plane_m → fold hits into out[m]` over existing planes — same plane traffic as K passes, value-lane traffic of the union instead of the sum. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index cbbd0fb1f..bd5dae1fa 100644 --- a/.claude/board/entries/README.md +++ b/.claude/board/entries/README.md @@ -62,6 +62,7 @@ exactly one of them, never both. | 2026-09-30 | `deepnsm-v2-coverage-bands` | | [2026-09-30-deepnsm-v2-coverage-bands.md](2026-09-30-deepnsm-v2-coverage-bands.md) | | 2026-09-30 | `cypher-mask-v2-is-a-replacement-not-a-phase` | | [2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md](2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md) | | 2026-09-29 | `grouped-moments-fold-anova-over-masks` | | [2026-09-29-grouped-moments-fold-anova-over-masks.md](2026-09-29-grouped-moments-fold-anova-over-masks.md) | +| 2026-09-29 | `grouped-cross-moments-fold` | | [2026-09-29-grouped-cross-moments-fold.md](2026-09-29-grouped-cross-moments-fold.md) | | 2026-09-29 | `deepnsm-v2-counted-pick-tag-deltas` | | [2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md](2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md) | | 2026-09-26 | `deepnsm-v2-lexical-evidence-survives-routing` | | [2026-09-26-deepnsm-v2-lexical-evidence-survives-routing.md](2026-09-26-deepnsm-v2-lexical-evidence-survives-routing.md) | | 2026-09-25 | `window-scheduling-and-two-level-ternlog` | | [2026-09-25-window-scheduling-and-two-level-ternlog.md](2026-09-25-window-scheduling-and-two-level-ternlog.md) | From d3d6ee6ef02fdbad9e260b688495048a4e198e59 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 20:53:07 +0000 Subject: [PATCH 06/14] =?UTF-8?q?perturbation-sim:=20angle=5Fcov=20?= =?UTF-8?q?=E2=80=94=20push=20injection=20covariance=20through=20L?= =?UTF-8?q?=E2=81=BA=20via=20CovHighD::sandwich?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sigma_theta = L+ Sigma_p L+ is the stochastic twin of pseudo_apply and is exact (not first-order) for DC flow in a fixed topology. L+ is symmetric, so it is ndarray's Pillar-9 sandwich M·Σ·M with M = L+; no local kernel. New default-off feature 'pillar' (dep:ndarray + ndarray/pillar). Stated limits: f32 inside the sandwich (measured ~1e-7 rel vs an f64 triple product), const-generic N checked against the runtime bus count, and an asymmetric Sigma_p is refused rather than silently lower-triangled. Out of scope, documented: line trips (a rank-1 update of L+, not a sandwich) and line-flow covariance (non-symmetric rectangular Jacobian). Tests: rank-1 case against pseudo_apply's outer product (never forms L+), f64 dense triple product, 40k-sample Monte Carlo of the deterministic solver (rel err 0.0009), N-mismatch and asymmetry refusals. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/perturbation-sim/Cargo.toml | 4 + crates/perturbation-sim/src/angle_cov.rs | 241 +++++++++++++++++++++++ crates/perturbation-sim/src/lib.rs | 5 + 3 files changed, 250 insertions(+) create mode 100644 crates/perturbation-sim/src/angle_cov.rs diff --git a/crates/perturbation-sim/Cargo.toml b/crates/perturbation-sim/Cargo.toml index 37d522f83..6dd535a92 100644 --- a/crates/perturbation-sim/Cargo.toml +++ b/crates/perturbation-sim/Cargo.toml @@ -18,6 +18,10 @@ description = "Spectral + edge-propagation outage simulator: models the perturba # (or target-cpu=native locally). All SIMD comes from `ndarray::simd` per the # workspace rule — never raw intrinsics here. ndarray-simd = ["dep:ndarray"] +# Push a nodal injection covariance through the DC model: Σθ = L⁺ Σp L⁺, +# computed by ndarray's Pillar-9 `CovHighD::sandwich` (f32, dense O(N³), +# const-generic N). Default OFF. Pulls ndarray's `pillar` (= linalg + splat3d). +pillar = ["dep:ndarray", "ndarray/pillar"] [dependencies] # The AdaWorldAPI fork ("The Foundation"), optional + behind `ndarray-simd`, diff --git a/crates/perturbation-sim/src/angle_cov.rs b/crates/perturbation-sim/src/angle_cov.rs new file mode 100644 index 000000000..dc2dc25aa --- /dev/null +++ b/crates/perturbation-sim/src/angle_cov.rs @@ -0,0 +1,241 @@ +//! Injection-uncertainty push-forward: `Σθ = L⁺ · Σp · L⁺`. +//! +//! The stochastic twin of [`crate::eigen::Eigen::pseudo_apply`]. For a fixed +//! topology the DC model `θ = L⁺ p` is linear, so a nodal injection covariance +//! `Σp` maps to the angle covariance `Σθ = L⁺ Σp (L⁺)ᵀ` **exactly**, not as a +//! first-order approximation. `L⁺` is symmetric, so this is the sandwich +//! `M·Σ·M` with `M = L⁺`, and it is computed by ndarray's certified +//! [`CovHighD::sandwich`] (Pillar-9) rather than by a second kernel here. +//! +//! # What this does NOT cover +//! +//! - **A line trip.** That changes `L`, so it is a rank-1 update of `L⁺` +//! (Sherman–Morrison / LODF), not a sandwich of the old covariance. Push a +//! post-contingency covariance by re-decomposing the post-trip Laplacian. +//! - **Line-flow covariance.** `f = B·A·θ` has a rectangular, non-symmetric +//! Jacobian, and `sandwich` assumes a symmetric `M`. Per-line variance is a +//! 2-sparse quadratic form over `Σθ` and does not need a sandwich at all. +//! +//! # Precision and size +//! +//! `CovHighD` is `f32` and dense `O(N³)`; this crate is `f64`. Values are +//! narrowed to `f32` once on the way in and widened on the way out, so expect +//! roughly `1e-6` relative error against an `f64` triple product. `N` is a +//! compile-time constant in `CovHighD`, while a [`crate::graph::Grid`]'s bus +//! count is a runtime value: the caller picks `N` and the call panics when the +//! decomposition disagrees with it. A runtime-sized sandwich would be a change +//! to ndarray, not something to route around here. + +use ndarray::hpc::pillar::cov_high_d::CovHighD; + +use crate::eigen::Eigen; + +/// Relative asymmetry tolerated in `Σp` before it is rejected. +/// +/// `CovHighD::from_symmetric_fn` reads only the lower triangle, so an +/// asymmetric input would be silently replaced by its lower half. Refusing is +/// the honest alternative to that silent substitution. +pub const SYMMETRY_TOL: f64 = 1e-9; + +/// Angle covariance `Σθ = L⁺ Σp L⁺` for a nodal injection covariance `Σp`. +/// +/// `eig` is the decomposition of the Laplacian (see +/// [`crate::eigen::symmetric_eigen`]); `rel_tol` is the null-space cutoff +/// passed to [`Eigen::pseudo_inverse`], identical to the one `pseudo_apply` +/// uses. `sigma_p` and the result are row-major `N×N`. +/// +/// # Panics +/// If `eig.n != N`, if `sigma_p.len() != N*N`, or if `sigma_p` is not +/// symmetric within [`SYMMETRY_TOL`] (relative to its largest entry). +pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f64) -> Vec { + assert_eq!( + eig.n, N, + "decomposition has {} buses, CovHighD has N = {N}", + eig.n + ); + assert_eq!(sigma_p.len(), N * N, "sigma_p must be N*N"); + let scale = sigma_p + .iter() + .fold(0.0_f64, |m, v| m.max(v.abs())) + .max(f64::MIN_POSITIVE); + for i in 0..N { + for j in 0..i { + let d = (sigma_p[i * N + j] - sigma_p[j * N + i]).abs(); + assert!( + d <= SYMMETRY_TOL * scale, + "sigma_p is not symmetric at ({i},{j}): off by {d}" + ); + } + } + + let l_plus = eig.pseudo_inverse(rel_tol); + let m = CovHighD::::from_symmetric_fn(|i, j| l_plus[i * N + j] as f32); + let s = CovHighD::::from_symmetric_fn(|i, j| sigma_p[i * N + j] as f32); + let out = s.sandwich(&m); + + let mut dense = vec![0.0_f64; N * N]; + for i in 0..N { + for j in 0..N { + dense[i * N + j] = out.get(i, j) as f64; + } + } + dense +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::eigen::symmetric_eigen; + use crate::graph::{Edge, Grid}; + + const N: usize = 6; + const TOL: f64 = 1e-9; + + /// A 6-bus ring with two chords and unequal susceptances. + fn grid() -> Grid { + Grid::new( + N, + vec![ + Edge::new(0, 1, 4.0, 10.0), + Edge::new(1, 2, 2.5, 10.0), + Edge::new(2, 3, 3.0, 10.0), + Edge::new(3, 4, 1.5, 10.0), + Edge::new(4, 5, 5.0, 10.0), + Edge::new(5, 0, 2.0, 10.0), + Edge::new(0, 3, 1.0, 10.0), + Edge::new(1, 4, 0.7, 10.0), + ], + ) + } + + fn rel_err(a: &[f64], b: &[f64]) -> f64 { + let scale = b.iter().fold(0.0_f64, |m, v| m.max(v.abs())); + a.iter() + .zip(b) + .fold(0.0_f64, |m, (x, y)| m.max((x - y).abs())) + / scale + } + + /// `Σp = A Aᵀ` for a fixed non-trivial `A` — SPD and genuinely coupled. + fn sigma_p() -> Vec { + let a: Vec = (0..N * N) + .map(|k| ((k * 7 + 3) % 11) as f64 / 11.0 - 0.4) + .collect(); + let mut s = vec![0.0; N * N]; + for i in 0..N { + for j in 0..N { + s[i * N + j] = (0..N).map(|k| a[i * N + k] * a[j * N + k]).sum(); + } + } + s + } + + #[test] + fn rank_one_injection_matches_pseudo_apply_outer_product() { + // Σp = p pᵀ ⇒ Σθ = (L⁺p)(L⁺p)ᵀ. The right-hand side comes from the + // eigen-coefficient path `pseudo_apply`, which never forms L⁺ — so this + // cross-checks the dense L⁺ AND the sandwich against an independent route. + let eig = symmetric_eigen(&grid().laplacian(), N); + let p = [1.2, -0.4, 0.9, -1.5, 0.3, -0.5]; // balanced: sums to 0 + let sp: Vec = (0..N * N).map(|k| p[k / N] * p[k % N]).collect(); + let theta = eig.pseudo_apply(&p, TOL); + let want: Vec = (0..N * N).map(|k| theta[k / N] * theta[k % N]).collect(); + let got = angle_covariance::(&eig, &sp, TOL); + let e = rel_err(&got, &want); + eprintln!("rank-1 push-forward rel err {e:e}"); + assert!(e < 1e-5, "rank-1 push-forward rel err {e:e}"); + assert!( + want.iter().any(|v| v.abs() > 1e-3), + "fixture is vacuous: θ is ~0" + ); + } + + #[test] + fn matches_f64_dense_triple_product() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let sp = sigma_p(); + let l = eig.pseudo_inverse(TOL); + let mut want = vec![0.0; N * N]; + for i in 0..N { + for m in 0..N { + let mut acc = 0.0; + for j in 0..N { + for k in 0..N { + acc += l[i * N + j] * sp[j * N + k] * l[k * N + m]; + } + } + want[i * N + m] = acc; + } + } + let got = angle_covariance::(&eig, &sp, TOL); + let e = rel_err(&got, &want); + eprintln!("sandwich vs f64 triple product rel err {e:e}"); + assert!(e < 1e-5, "sandwich vs f64 triple product rel err {e:e}"); + // Non-trivial: off-diagonal coupling must be present, else a diagonal-only + // implementation would pass. + assert!( + want[1].abs() > 1e-3 * want[0].abs(), + "fixture has no off-diagonal mass" + ); + } + + #[test] + fn matches_monte_carlo_of_the_deterministic_solver() { + // Draw p ~ N(0, Σp) via p = A z, solve θ = L⁺p with pseudo_apply, and + // compare the empirical covariance to the sandwich. This checks the + // STATISTICAL claim (push-forward of a distribution), not just algebra. + let eig = symmetric_eigen(&grid().laplacian(), N); + let a: Vec = (0..N * N) + .map(|k| ((k * 7 + 3) % 11) as f64 / 11.0 - 0.4) + .collect(); + let mut state = 0x9E37_79B9_7F4A_7C15_u64; + let mut unif = move || { + state = state.wrapping_add(0x9E37_79B9_7F4A_7C15); + let mut z = state; + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + ((z ^ (z >> 31)) >> 11) as f64 / (1u64 << 53) as f64 + }; + let samples = 40_000; + let mut emp = vec![0.0; N * N]; + for _ in 0..samples { + let z: Vec = (0..N) + .map(|_| { + let (u1, u2) = (unif().max(1e-300), unif()); + (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos() + }) + .collect(); + let p: Vec = (0..N) + .map(|i| (0..N).map(|k| a[i * N + k] * z[k]).sum()) + .collect(); + let th = eig.pseudo_apply(&p, TOL); + for i in 0..N { + for j in 0..N { + emp[i * N + j] += th[i] * th[j]; + } + } + } + emp.iter_mut().for_each(|v| *v /= samples as f64); + let got = angle_covariance::(&eig, &sigma_p(), TOL); + let e = rel_err(&emp, &got); + // Sampling error of a second moment at n = 40k is ~1/√n·√2 ≈ 0.7%. + eprintln!("Monte-Carlo vs sandwich rel err {e:.4}"); + assert!(e < 0.03, "Monte-Carlo vs sandwich rel err {e:.4}"); + } + + #[test] + #[should_panic(expected = "CovHighD has N")] + fn refuses_a_const_n_that_disagrees_with_the_grid() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let _ = angle_covariance::<4>(&eig, &[0.0; 16], TOL); + } + + #[test] + #[should_panic(expected = "not symmetric")] + fn refuses_an_asymmetric_injection_covariance() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let mut sp = sigma_p(); + sp[N] += 0.5; // (1,0) no longer equals (0,1) + let _ = angle_covariance::(&eig, &sp, TOL); + } +} diff --git a/crates/perturbation-sim/src/lib.rs b/crates/perturbation-sim/src/lib.rs index 737d8c2a5..5b416d6b2 100644 --- a/crates/perturbation-sim/src/lib.rs +++ b/crates/perturbation-sim/src/lib.rs @@ -32,6 +32,9 @@ //! - `jc::ewa_sandwich` is a genuine covariance Σ-push-forward along multi-hop //! edge paths — the *uncertainty-propagation* sibling of the deterministic //! flow cascade here. +//! - `angle_cov` (feature `pillar`) is the nodal version of that sibling: +//! `Σθ = L⁺ Σp L⁺`, exact for a fixed topology, computed by ndarray's +//! Pillar-9 `CovHighD::sandwich` rather than a local kernel. //! //! ## Statistical hand-off //! @@ -50,6 +53,8 @@ //! buses) — exactly the regime of a regional transmission graph. pub mod acflow; +#[cfg(feature = "pillar")] +pub mod angle_cov; pub mod basin; pub mod buffer; pub mod cascade; From a54824a042611994892175b08b0d4b0e41313a66 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 29 Sep 2026 20:53:33 +0000 Subject: [PATCH 07/14] board: perturbation-sim angle covariance via CovHighD sandwich Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- ...ion-sim-angle-covariance-via-cov-high-d.md | 31 +++++++++++++++++++ .claude/board/entries/README.md | 3 +- 2 files changed, 33 insertions(+), 1 deletion(-) create mode 100644 .claude/board/entries/2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md diff --git a/.claude/board/entries/2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md b/.claude/board/entries/2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md new file mode 100644 index 000000000..831e13769 --- /dev/null +++ b/.claude/board/entries/2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md @@ -0,0 +1,31 @@ +# 2026-09-29 — perturbation-sim: Σθ = L⁺ Σp L⁺ through ndarray CovHighD::sandwich + +**STATUS:** MEASURED · **Scope:** `crates/perturbation-sim/src/angle_cov.rs` (feature `pillar`), ndarray `hpc::pillar::cov_high_d` + +## What landed +- ndarray `4d4ee17`: `CovHighD::from_symmetric_fn` + public `get`. Before this a consumer + could not build a CovHighD from data without re-spelling the packed index. +- lance-graph `bb6d044`: `angle_covariance::(eig, Σp, rel_tol)`, the stochastic twin of + `pseudo_apply`. It is exact for DC flow in a fixed topology. It calls the Pillar-9 sandwich; + there is no local kernel. + +## Measured (6-bus ring + chords) +- Against an f64 dense triple product: rel err 1.2e-7. This is the f32 floor. +- Rank-1 Σp = ppᵀ against `pseudo_apply` outer product: 2.3e-7. That path never forms L⁺. +- 40k-sample Monte Carlo of the deterministic solver: rel err 0.0009. +- Disable runs, each red under exactly its own test: + - M := I; + - symmetry guard removed; + - N-mismatch guard removed. + +## Finding in ndarray +Every existing `sandwich` test used M = I. Under a deliberately broken kernel (M·Σ·Σ), +`sandwich_identity_is_identity` stayed green. The new dense non-identity test is the first +that can see an index/transpose defect in the Pillar-9 kernel. + +## OPEN +- CovHighD is const-generic N; `Grid::n` is runtime. The caller picks N and it is checked. + A runtime-sized sandwich is an ndarray change, not taken. +- Line trips are a rank-1 update of L⁺ (LODF), not a sandwich. They are not wired. +- Line-flow covariance has a non-symmetric rectangular Jacobian, which sandwich cannot + express. Per-line variance is a 2-sparse quadratic form over Σθ, not built. diff --git a/.claude/board/entries/README.md b/.claude/board/entries/README.md index bd5dae1fa..39dbe7344 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. -208 entries, 2026-08-06 .. 2026-10-04. +211 entries, 2026-08-06 .. 2026-10-04. | date | entry id | finding | file | |---|---|---|---| @@ -61,6 +61,7 @@ exactly one of them, never both. | 2026-09-30 | `three-reference-sets-are-not-ordinal-aligned` | | [2026-09-30-three-reference-sets-are-not-ordinal-aligned.md](2026-09-30-three-reference-sets-are-not-ordinal-aligned.md) | | 2026-09-30 | `deepnsm-v2-coverage-bands` | | [2026-09-30-deepnsm-v2-coverage-bands.md](2026-09-30-deepnsm-v2-coverage-bands.md) | | 2026-09-30 | `cypher-mask-v2-is-a-replacement-not-a-phase` | | [2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md](2026-09-30-cypher-mask-v2-is-a-replacement-not-a-phase.md) | +| 2026-09-29 | `perturbation-sim-angle-covariance-via-cov-high-d` | | [2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md](2026-09-29-perturbation-sim-angle-covariance-via-cov-high-d.md) | | 2026-09-29 | `grouped-moments-fold-anova-over-masks` | | [2026-09-29-grouped-moments-fold-anova-over-masks.md](2026-09-29-grouped-moments-fold-anova-over-masks.md) | | 2026-09-29 | `grouped-cross-moments-fold` | | [2026-09-29-grouped-cross-moments-fold.md](2026-09-29-grouped-cross-moments-fold.md) | | 2026-09-29 | `deepnsm-v2-counted-pick-tag-deltas` | | [2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md](2026-09-29-deepnsm-v2-counted-pick-tag-deltas.md) | From 4d3e3901660010e4d5b4b3de9e56285574064420 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 30 Sep 2026 22:49:38 +0000 Subject: [PATCH 08/14] mask-risc + jc: consume ndarray's PowerSums / CrossPowerSums (drop the duplicate record) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ndarray master shipped PowerSums { n: u64, sum: i64, sum_sq: u128 } — the same record this branch had introduced as GroupMoments. The ndarray side of this branch was rebuilt on master to add only checked_merge and the joint fold (CrossPowerSums), so consumers move to the upstream names: - GroupMoments -> PowerSums, GroupCrossMoments -> CrossPowerSums - masked_group_(cross_)moments_i32* -> masked_group_(cross_)power_sums_i32* - sum_x2/sum_y2 -> sum_x_sq/sum_y_sq; x_moments()/y_moments() -> x()/y() - EMPTY -> default() Square sums are u128 upstream. jc's exact i128 centring converts them with i128::try_from; one past i128::MAX is already outside the documented 2^32-row bound and returns None, the existing policy. mask-risc's own IR names (Terminal::GroupMomentsI32, Value::GroupMoments) are unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- .../2026-09-29-grouped-cross-moments-fold.md | 2 + ...9-grouped-moments-fold-anova-over-masks.md | 2 + crates/jc/src/stats.rs | 132 ++++++++++-------- crates/lance-graph-mask-risc/src/exec.rs | 59 ++++---- crates/lance-graph-mask-risc/src/ir.rs | 6 +- crates/lance-graph-mask-risc/src/lib.rs | 2 +- crates/lance-graph-mask-risc/src/reference.rs | 14 +- crates/lance-graph-mask-risc/src/value.rs | 20 +-- crates/lance-graph-mask-risc/tests/foreign.rs | 28 ++-- 9 files changed, 145 insertions(+), 120 deletions(-) diff --git a/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md index c64b1ab6b..d22a68051 100644 --- a/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md +++ b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md @@ -15,3 +15,5 @@ ## Open - Multi-membership in one pass: `lance-graph-report` lowers `CoordSpec::MaskSet` to one `Filter::Plane` pass per member tuple and per fold state. The single-group limit is in `GroupKeyAddr::group_of` (one `Option` per row), in the IR's `GroupKey`, and in report lowering. Seam: a word-level walker `for each 64-row word: sel & plane_m → fold hits into out[m]` over existing planes — same plane traffic as K passes, value-lane traffic of the union instead of the sum. + +**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupMoments`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. diff --git a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md index 76e07efaa..5c2b0f749 100644 --- a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md +++ b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md @@ -14,3 +14,5 @@ Mask → fold → `anova_from_moments` vs materialized → `anova_one_way`, 42 n - Multi-membership is NOT one pass: `GroupKey` resolves one group per row, and `lance-graph-report`'s `CoordSpec::MaskSet` executes one `Filter::Plane` program per member. Seam: a key address that walks the selected rows once and folds each row into every member plane that holds it. - Bivariate `(n, Σx, Σy, Σx², Σy², Σxy)`: same walk, a second value lane in the closure, a wider slot type; not built. - No early exit, and no `CausalEdge64` commit. ANOVA is non-monotone under future rows. + +**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupMoments`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index e96daf5d0..738e8b454 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -74,7 +74,7 @@ //! front by the same `all_finite` guard. use crate::reliability::{all_finite, mean, pearson, pearson_from_centered}; -use ndarray::simd::{GroupCrossMoments, GroupMoments}; +use ndarray::simd::{CrossPowerSums, PowerSums}; use std::collections::BTreeSet; // ─────────────────────────── local helpers ─────────────────────────── @@ -1195,7 +1195,7 @@ fn anova_from_ss(ss_b: f64, ss_w: f64, ss_t: f64, k: usize, n_total: usize) -> O // // The same statistics as `anova_one_way` / `eta_squared`, projected from // per-group `(n, Σx, Σx²)` instead of from materialized group vectors — the -// shape `ndarray::simd::masked_group_moments_i32` folds out of a population +// shape `ndarray::simd::masked_group_power_sums_i32` folds out of a population // mask in one pass. The statistical policy (every `None` case, the F / p / // η² formulas) is shared with the slice path through `anova_from_ss` and // `eta_from_ss`; only the sums of squares are formed differently. @@ -1219,8 +1219,8 @@ fn anova_from_ss(ss_b: f64, ss_w: f64, ss_t: f64, k: usize, n_total: usize) -> O /// `None` on fewer than 2 groups, any empty group (`n == 0`, matching /// [`anova_one_way`]'s empty-group rule), or moments too large for the exact /// `i128` forms — i.e. outside the `2^32`-rows-per-group bound -/// `GroupMoments` is exact under. -fn one_way_ss_from_moments(groups: &[GroupMoments]) -> Option<(f64, f64, f64, usize, usize)> { +/// `PowerSums` is exact under. +fn one_way_ss_from_moments(groups: &[PowerSums]) -> Option<(f64, f64, f64, usize, usize)> { if groups.len() < 2 || groups.iter().any(|g| g.n == 0) { return None; } @@ -1238,7 +1238,7 @@ fn one_way_ss_from_moments(groups: &[GroupMoments]) -> Option<(f64, f64, f64, us let n_g = i128::from(g.n); let s_g = i128::from(g.sum); let w_g = n_g - .checked_mul(g.sum_sq)? + .checked_mul(i128::try_from(g.sum_sq).ok()?)? .checked_sub(s_g.checked_mul(s_g)?)?; if w_g < 0 { return None; // impossible for moments of real data: inconsistent input @@ -1259,7 +1259,7 @@ fn one_way_ss_from_moments(groups: &[GroupMoments]) -> Option<(f64, f64, f64, us /// [`anova_one_way`] on the materialized groups, without the groups. /// /// `groups[g]` holds group `g`'s `(n, Σx, Σx²)`, e.g. as folded out of a -/// population mask by `ndarray::simd::masked_group_moments_i32` (or +/// population mask by `ndarray::simd::masked_group_power_sums_i32` (or /// `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32`). Every degenerate /// case is [`anova_one_way`]'s: fewer than 2 groups, an empty group, /// `N ≤ k`, and zero within-group variance are all `None`. The sums of @@ -1269,17 +1269,17 @@ fn one_way_ss_from_moments(groups: &[GroupMoments]) -> Option<(f64, f64, f64, us /// /// ``` /// use jc::stats::{anova_from_moments, anova_one_way}; -/// use ndarray::simd::GroupMoments; +/// use ndarray::simd::PowerSums; /// /// let groups = vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]; /// let moments = [ -/// GroupMoments { n: 3, sum: 6, sum_sq: 14 }, -/// GroupMoments { n: 3, sum: 15, sum_sq: 77 }, +/// PowerSums { n: 3, sum: 6, sum_sq: 14 }, +/// PowerSums { n: 3, sum: 15, sum_sq: 77 }, /// ]; /// let (a, b) = (anova_one_way(&groups).unwrap(), anova_from_moments(&moments).unwrap()); /// assert!((a.f - b.f).abs() < 1e-12 && (a.p - b.p).abs() < 1e-12); /// ``` -pub fn anova_from_moments(groups: &[GroupMoments]) -> Option { +pub fn anova_from_moments(groups: &[PowerSums]) -> Option { let (ss_b, ss_w, ss_t, k, n_total) = one_way_ss_from_moments(groups)?; anova_from_ss(ss_b, ss_w, ss_t, k, n_total) } @@ -1288,7 +1288,7 @@ pub fn anova_from_moments(groups: &[GroupMoments]) -> Option { /// degeneracy policy (`None` only on fewer than 2 groups, an empty group, or /// zero total variance; a perfectly separated layout with zero within-group /// variance is a real `1.0`, unlike the F test). -pub fn eta_squared_from_moments(groups: &[GroupMoments]) -> Option { +pub fn eta_squared_from_moments(groups: &[PowerSums]) -> Option { let (ss_b, _, ss_t, _, _) = one_way_ss_from_moments(groups)?; eta_from_ss(ss_b, ss_t) } @@ -1297,7 +1297,7 @@ pub fn eta_squared_from_moments(groups: &[GroupMoments]) -> Option { // // Pearson, the sample covariance, the simple least-squares line and its R², // projected from one group's `(n, Σx, Σy, Σx², Σy², Σxy)` — the shape -// `ndarray::simd::masked_group_cross_moments_i32` folds off a population +// `ndarray::simd::masked_group_cross_power_sums_i32` folds off a population // mask. Each shares its acceptance policy with the slice implementation it // mirrors (`pearson_from_centered`, `sample_cov_tail`, `r_squared_tail`); only // the centred sums are formed differently. @@ -1308,16 +1308,23 @@ pub fn eta_squared_from_moments(groups: &[GroupMoments]) -> Option { /// denominators. No two large floats are ever subtracted, so a huge common /// offset with a tiny spread loses nothing. /// -/// Every product fits: within `GroupCrossMoments`' `2^32`-row bound each of +/// The square sums arrive as `u128` (ndarray folds them unsigned); one past +/// `i128::MAX` is already outside the bound below and yields `None`. +/// +/// Every product fits: within `CrossPowerSums`' `2^32`-row bound each of /// `n·Σx²`, `(Σx)²`, `n·Σxy`, `Σx·Σy` is at most `2^126` in magnitude, and by /// Cauchy–Schwarz so is each difference. Past the bound the checked /// arithmetic returns `None` rather than wrap. -fn centered_cross(m: &GroupCrossMoments) -> Option<(i128, i128, i128)> { +fn centered_cross(m: &CrossPowerSums) -> Option<(i128, i128, i128)> { let n = i128::from(m.n); let (sx, sy) = (i128::from(m.sum_x), i128::from(m.sum_y)); let cxy = n.checked_mul(m.sum_xy)?.checked_sub(sx.checked_mul(sy)?)?; - let cxx = n.checked_mul(m.sum_x2)?.checked_sub(sx.checked_mul(sx)?)?; - let cyy = n.checked_mul(m.sum_y2)?.checked_sub(sy.checked_mul(sy)?)?; + let cxx = n + .checked_mul(i128::try_from(m.sum_x_sq).ok()?)? + .checked_sub(sx.checked_mul(sx)?)?; + let cyy = n + .checked_mul(i128::try_from(m.sum_y_sq).ok()?)? + .checked_sub(sy.checked_mul(sy)?)?; // A negative centred square is impossible for moments of real data. (cxx >= 0 && cyy >= 0).then_some((cxy, cxx, cyy)) } @@ -1327,12 +1334,12 @@ fn centered_cross(m: &GroupCrossMoments) -> Option<(i128, i128, i128)> { /// /// ``` /// use jc::stats::pearson_from_cross_moments; -/// use ndarray::simd::GroupCrossMoments; +/// use ndarray::simd::CrossPowerSums; /// // x = [1,2,3], y = [2,4,6]: perfectly correlated. -/// let m = GroupCrossMoments { n: 3, sum_x: 6, sum_y: 12, sum_x2: 14, sum_y2: 56, sum_xy: 28 }; +/// let m = CrossPowerSums { n: 3, sum_x: 6, sum_y: 12, sum_x_sq: 14, sum_y_sq: 56, sum_xy: 28 }; /// assert!((pearson_from_cross_moments(&m).unwrap() - 1.0).abs() < 1e-12); /// ``` -pub fn pearson_from_cross_moments(m: &GroupCrossMoments) -> Option { +pub fn pearson_from_cross_moments(m: &CrossPowerSums) -> Option { if m.n < 2 { return None; } @@ -1342,7 +1349,7 @@ pub fn pearson_from_cross_moments(m: &GroupCrossMoments) -> Option { /// The unbiased (divisor `n−1`) sample covariance from grouped cross moments /// — the convention this module already uses (`None` for `n < 2`). -pub fn sample_covariance_from_cross_moments(m: &GroupCrossMoments) -> Option { +pub fn sample_covariance_from_cross_moments(m: &CrossPowerSums) -> Option { let (cxy, _, _) = centered_cross(m)?; // `C_xy = n·Σ(x−x̄)(y−ȳ)`: divide the `n` back out before the tail. let sxy = cxy as f64 / m.n.max(1) as f64; @@ -1370,7 +1377,7 @@ pub struct SimpleRegression { /// `|β·x̄|` is large and `α` is small its ABSOLUTE error is on the order of /// `ε·(|β|·|x̄| + |ȳ|)`. That is the conditioning of the quantity itself, not /// of the fold; centre `x` before regressing when `α` itself matters. -pub fn simple_regression_from_cross_moments(m: &GroupCrossMoments) -> Option { +pub fn simple_regression_from_cross_moments(m: &CrossPowerSums) -> Option { if m.n < 2 { return None; } @@ -1388,7 +1395,7 @@ pub fn simple_regression_from_cross_moments(m: &GroupCrossMoments) -> Option Option { +pub fn r_squared_from_cross_moments(m: &CrossPowerSums) -> Option { if m.n < 3 { return None; } @@ -2486,7 +2493,7 @@ mod tests { // The claim under test: a statistic computed from grouped sufficient // statistics folded straight off a population mask equals the same // statistic computed from the materialized groups. The fold is - // `ndarray::simd::masked_group_moments_i32` — the kernel + // `ndarray::simd::masked_group_power_sums_i32` — the kernel // `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32` delegates to. // // Tolerances. Both paths form the same sums of squares by different @@ -2501,7 +2508,7 @@ mod tests { mod moments_equivalence { use super::super::*; - use ndarray::simd::masked_group_moments_i32; + use ndarray::simd::masked_group_power_sums_i32; fn lcg(s: &mut u64) -> u64 { *s = s @@ -2561,9 +2568,9 @@ mod tests { } /// The substrate path: one masked fold, no observation copied out. - fn fold(p: &Pop, k: u32) -> Vec { - let mut out = vec![GroupMoments::EMPTY; k as usize]; - masked_group_moments_i32(&p.mask, &p.keys, &p.values, &mut out); + fn fold(p: &Pop, k: u32) -> Vec { + let mut out = vec![PowerSums::default(); k as usize]; + masked_group_power_sums_i32(&p.mask, &p.keys, &p.values, &mut out); out } @@ -2637,7 +2644,7 @@ mod tests { }); let whole = fold(&p, k); for chunks in [2usize, 3, 7, 47] { - let mut merged = vec![GroupMoments::EMPTY; k as usize]; + let mut merged = vec![PowerSums::default(); k as usize]; let step = n.div_ceil(chunks); for c in 0..chunks { let (lo, hi) = (c * step, ((c + 1) * step).min(n)); @@ -2645,8 +2652,8 @@ mod tests { for i in lo..hi { m[i / 64] |= p.mask[i / 64] & (1 << (i % 64)); } - let mut part = vec![GroupMoments::EMPTY; k as usize]; - masked_group_moments_i32(&m, &p.keys, &p.values, &mut part); + let mut part = vec![PowerSums::default(); k as usize]; + masked_group_power_sums_i32(&m, &p.keys, &p.values, &mut part); for (acc, x) in merged.iter_mut().zip(part) { *acc = acc.checked_merge(x).expect("in bound"); } @@ -2677,19 +2684,19 @@ mod tests { /// zero within-group variance. #[test] fn degenerate_layouts_match_the_slice_contract() { - let m = |n: u64, xs: &[i64]| GroupMoments { + let m = |n: u64, xs: &[i64]| PowerSums { n, sum: xs.iter().sum(), - sum_sq: xs.iter().map(|&x| i128::from(x * x)).sum(), + sum_sq: xs.iter().map(|&x| (x * x) as u128).sum(), }; let f = |xs: &[i64]| xs.iter().map(|&x| x as f64).collect::>(); - type Case = (&'static str, Vec>, Vec); + type Case = (&'static str, Vec>, Vec); let cases: Vec = vec![ ("one group", vec![f(&[1, 2, 3])], vec![m(3, &[1, 2, 3])]), ( "absent group", vec![f(&[1, 2, 3]), vec![], f(&[4, 5])], - vec![m(3, &[1, 2, 3]), GroupMoments::EMPTY, m(2, &[4, 5])], + vec![m(3, &[1, 2, 3]), PowerSums::default(), m(2, &[4, 5])], ), ( "N <= k", @@ -2734,10 +2741,13 @@ mod tests { let mut groups = Vec::new(); for b in bases { let xs: Vec = offsets.iter().map(|o| b + o).collect(); - moments.push(GroupMoments { + moments.push(PowerSums { n: xs.len() as u64, sum: xs.iter().sum(), - sum_sq: xs.iter().map(|&x| i128::from(x) * i128::from(x)).sum(), + sum_sq: xs + .iter() + .map(|&x| (i128::from(x) * i128::from(x)) as u128) + .sum(), }); groups.push(xs.iter().map(|&x| x as f64).collect::>()); } @@ -2779,7 +2789,7 @@ mod tests { mod cross_moments_equivalence { use super::super::*; use crate::reliability::pearson; - use ndarray::simd::masked_group_cross_moments_i32; + use ndarray::simd::masked_group_cross_power_sums_i32; fn lcg(s: &mut u64) -> u64 { *s = s @@ -2830,9 +2840,9 @@ mod tests { p } - fn fold(p: &Pop, k: u32) -> Vec { - let mut out = vec![GroupCrossMoments::EMPTY; k as usize]; - masked_group_cross_moments_i32(&p.mask, &p.keys, &p.xs, &p.ys, &mut out); + fn fold(p: &Pop, k: u32) -> Vec { + let mut out = vec![CrossPowerSums::default(); k as usize]; + masked_group_cross_power_sums_i32(&p.mask, &p.keys, &p.xs, &p.ys, &mut out); out } @@ -2871,7 +2881,7 @@ mod tests { /// Compare every projection of one group against its slice /// counterpart; both `None` or both within tolerance. - fn assert_group_agrees(x: &[f64], y: &[f64], m: &GroupCrossMoments, what: &str) { + fn assert_group_agrees(x: &[f64], y: &[f64], m: &CrossPowerSums, what: &str) { match (pearson(x, y), pearson_from_cross_moments(m)) { (Some(a), Some(b)) => assert!(rel(a, b) < 1e-9, "{what}: r {a} vs {b}"), (a, b) => assert_eq!(a, b, "{what}: r degeneracy"), @@ -2958,20 +2968,20 @@ mod tests { let whole = fold(&p, k); for chunks in [2usize, 5, 47] { let step = n.div_ceil(chunks); - let parts: Vec> = (0..chunks) + let parts: Vec> = (0..chunks) .map(|c| { let mut m = vec![0u64; n.div_ceil(64)]; for i in c * step..((c + 1) * step).min(n) { m[i / 64] |= p.mask[i / 64] & (1 << (i % 64)); } - let mut part = vec![GroupCrossMoments::EMPTY; k as usize]; - masked_group_cross_moments_i32(&m, &p.keys, &p.xs, &p.ys, &mut part); + let mut part = vec![CrossPowerSums::default(); k as usize]; + masked_group_cross_power_sums_i32(&m, &p.keys, &p.xs, &p.ys, &mut part); part }) .collect(); for order in [false, true] { - let mut merged = vec![GroupCrossMoments::EMPTY; k as usize]; - let seq: Vec<&Vec> = if order { + let mut merged = vec![CrossPowerSums::default(); k as usize]; + let seq: Vec<&Vec> = if order { parts.iter().rev().collect() } else { parts.iter().collect() @@ -2993,12 +3003,20 @@ mod tests { } } - fn moments_of(xs: &[i64], ys: &[i64]) -> GroupCrossMoments { - let mut m = GroupCrossMoments::EMPTY; - for (&x, &y) in xs.iter().zip(ys) { - m.observe(x as i32, y as i32); - } - m + fn moments_of(xs: &[i64], ys: &[i64]) -> CrossPowerSums { + // Through the real fold: one group, every row selected. + let xs: Vec = xs.iter().map(|&x| x as i32).collect(); + let ys: Vec = ys.iter().map(|&y| y as i32).collect(); + let mut m = [CrossPowerSums::default()]; + let mask = vec![u64::MAX; xs.len().div_ceil(64)]; + ndarray::simd::masked_group_cross_power_sums_i32( + &mask, + &vec![0; xs.len()], + &xs, + &ys, + &mut m, + ); + m[0] } /// FAILS IF: a degenerate or boundary layout is treated differently @@ -3126,8 +3144,8 @@ mod tests { // Anti-vacuity: the naive float projection of the SAME moments fails. let n = m.n as f64; let naive = (n * m.sum_xy as f64 - m.sum_x as f64 * m.sum_y as f64) - / ((n * m.sum_x2 as f64 - (m.sum_x as f64).powi(2)).sqrt() - * (n * m.sum_y2 as f64 - (m.sum_y as f64).powi(2)).sqrt()); + / ((n * m.sum_x_sq as f64 - (m.sum_x as f64).powi(2)).sqrt() + * (n * m.sum_y_sq as f64 - (m.sum_y as f64).powi(2)).sqrt()); assert!( !naive.is_finite() || rel(naive, r_exact) > 1e-3, "fixture not adversarial: {naive}" @@ -3149,12 +3167,12 @@ mod tests { #[test] fn moments_past_the_bound_are_refused() { let v = (1i128 << 95) + 1; - let huge = GroupCrossMoments { + let huge = CrossPowerSums { n: 1 << 33, sum_x: 0, sum_y: 0, - sum_x2: v, - sum_y2: v, + sum_x_sq: v as u128, + sum_y_sq: v as u128, sum_xy: v, }; assert_eq!(pearson_from_cross_moments(&huge), None); diff --git a/crates/lance-graph-mask-risc/src/exec.rs b/crates/lance-graph-mask-risc/src/exec.rs index 72cc41f88..e951ee2b7 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -29,17 +29,18 @@ use ndarray::simd::{ lt_i32_to_mask_under, mask_all, mask_and, mask_and_assign, mask_andnot, mask_andnot_assign, mask_any, mask_gather_u32, mask_not, mask_not_assign, mask_or, mask_or_assign, mask_scatter_or_u32, mask_set_range, mask_xor, mask_xor_assign, masked_group_count_u32, - masked_group_count_u32_pair, masked_group_count_u32_via, masked_group_cross_moments_i32, - masked_group_cross_moments_i32_pair, masked_group_cross_moments_i32_via, masked_group_max_i32, - masked_group_max_i32_pair, masked_group_max_i32_via, masked_group_min_i32, - masked_group_min_i32_pair, masked_group_min_i32_via, masked_group_moments_i32, - masked_group_moments_i32_pair, masked_group_moments_i32_via, masked_group_sum_i32, - masked_group_sum_i32_via, masked_group_sum_sym_i32, masked_group_sum_sym_i32_pair, - masked_group_sum_sym_i32_via, masked_key_run_count_u32, masked_max_i32, masked_min_i32, - masked_strided_group_sum, masked_sum_i32, ne_i32_to_mask, ne_i32_to_mask_under, ne_u32_to_mask, - ne_u32_to_mask_under, popcount_batch_u64, ternary_match_strided_to_mask, - ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, ternary_match_u64_to_mask, - ternary_match_u64_to_mask_under, GroupCrossMoments, GroupMoments, KeyRunCarry, + masked_group_count_u32_pair, masked_group_count_u32_via, masked_group_cross_power_sums_i32, + masked_group_cross_power_sums_i32_pair, masked_group_cross_power_sums_i32_via, + masked_group_max_i32, masked_group_max_i32_pair, masked_group_max_i32_via, + masked_group_min_i32, masked_group_min_i32_pair, masked_group_min_i32_via, + masked_group_power_sums_i32, masked_group_power_sums_i32_pair, masked_group_power_sums_i32_via, + masked_group_sum_i32, masked_group_sum_i32_via, masked_group_sum_sym_i32, + masked_group_sum_sym_i32_pair, masked_group_sum_sym_i32_via, masked_key_run_count_u32, + masked_max_i32, masked_min_i32, masked_strided_group_sum, masked_sum_i32, ne_i32_to_mask, + ne_i32_to_mask_under, ne_u32_to_mask, ne_u32_to_mask_under, popcount_batch_u64, + ternary_match_strided_to_mask, ternary_match_u32_to_mask, ternary_match_u32_to_mask_under, + ternary_match_u64_to_mask, ternary_match_u64_to_mask_under, CrossPowerSums, KeyRunCarry, + PowerSums, }; use crate::ir::{ @@ -1576,9 +1577,9 @@ pub fn execute_compiled( } (Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. }, Out::I64(o)) => o.fill(0), (Terminal::GroupReduce { fold, .. }, Out::I64(o)) => o.fill(fold.seed()), - (Terminal::GroupMomentsI32 { .. }, Out::Moments(o)) => o.fill(GroupMoments::EMPTY), + (Terminal::GroupMomentsI32 { .. }, Out::Moments(o)) => o.fill(PowerSums::default()), (Terminal::GroupCrossMomentsI32 { .. }, Out::CrossMoments(o)) => { - o.fill(GroupCrossMoments::EMPTY) + o.fill(CrossPowerSums::default()) } _ => {} } @@ -1995,22 +1996,22 @@ pub fn execute_compiled( // `validate` already refused a missing/too-small `out`, every // wrong-width lane and a plane past the 2^32-row exactness // bound; one delegation per tile (law L3) into the sink - // seeded above with `GroupMoments::EMPTY`. + // seeded above with `PowerSums::default()`. if let Out::Moments(o) = &mut out { let m = read(planes, &slots, mask, t); let v = lane_i32(planes, val, t); match key { GroupKey::Lane(k) => { - masked_group_moments_i32(m, lane_u32(planes, k, t), v, o) + masked_group_power_sums_i32(m, lane_u32(planes, k, t), v, o) } - GroupKey::Via { fk, key } => masked_group_moments_i32_via( + GroupKey::Via { fk, key } => masked_group_power_sums_i32_via( m, lane_u32(planes, fk, t), foreign_lane_u32(foreign, key), v, o, ), - GroupKey::Pair { hi, lo, stride } => masked_group_moments_i32_pair( + GroupKey::Pair { hi, lo, stride } => masked_group_power_sums_i32_pair( m, lane_u32(planes, hi, t), lane_u32(planes, lo, t), @@ -2029,9 +2030,9 @@ pub fn execute_compiled( let (xs, ys) = (lane_i32(planes, x, t), lane_i32(planes, y, t)); match key { GroupKey::Lane(k) => { - masked_group_cross_moments_i32(m, lane_u32(planes, k, t), xs, ys, o) + masked_group_cross_power_sums_i32(m, lane_u32(planes, k, t), xs, ys, o) } - GroupKey::Via { fk, key } => masked_group_cross_moments_i32_via( + GroupKey::Via { fk, key } => masked_group_cross_power_sums_i32_via( m, lane_u32(planes, fk, t), foreign_lane_u32(foreign, key), @@ -2039,15 +2040,17 @@ pub fn execute_compiled( ys, o, ), - GroupKey::Pair { hi, lo, stride } => masked_group_cross_moments_i32_pair( - m, - lane_u32(planes, hi, t), - lane_u32(planes, lo, t), - stride, - xs, - ys, - o, - ), + GroupKey::Pair { hi, lo, stride } => { + masked_group_cross_power_sums_i32_pair( + m, + lane_u32(planes, hi, t), + lane_u32(planes, lo, t), + stride, + xs, + ys, + o, + ) + } } } } diff --git a/crates/lance-graph-mask-risc/src/ir.rs b/crates/lance-graph-mask-risc/src/ir.rs index d22ce5e2d..f33fb1db8 100644 --- a/crates/lance-graph-mask-risc/src/ir.rs +++ b/crates/lance-graph-mask-risc/src/ir.rs @@ -424,7 +424,7 @@ pub enum Terminal { /// Grouped sufficient statistics of an `I32` lane: for every row `i` /// where `mask` holds, resolves the row's group through [`GroupKey`] and /// folds `lanes[val][i]` into the caller's `Out::Moments` buffer as - /// `n += 1, Σx += x, Σx² += x²` ([`ndarray::simd::GroupMoments`]). One + /// `n += 1, Σx += x, Σx² += x²` ([`ndarray::simd::PowerSums`]). One /// pass over the selected rows; no selected value is copied out. The /// buffer's length IS the group universe `K`, and the same zero-fallback /// drops as [`Terminal::GroupReduce`] apply. @@ -437,7 +437,7 @@ pub enum Terminal { /// Exact under the same row bound as [`Terminal::MaskedSumI32`] /// ([`MASKED_SUM_I32_MAX_ROWS`], `2^32` rows): the `Σx` field is an `i64`, /// and no group can exceed the plane. The executor refuses a wider plane - /// rather than wrap. The sink is seeded with `GroupMoments::EMPTY` before + /// rather than wrap. The sink is seeded with `PowerSums::default()` before /// the first tile, so an empty group reads `n == 0`. GroupMomentsI32 { mask: Operand, @@ -447,7 +447,7 @@ pub enum Terminal { /// Grouped cross moments of TWO `I32` lanes: for every row `i` where /// `mask` holds, folds `(lanes[x][i], lanes[y][i])` into the caller's /// `Out::CrossMoments` buffer as `n, Σx, Σy, Σx², Σy², Σxy` - /// ([`ndarray::simd::GroupCrossMoments`]). The bivariate member of the + /// ([`ndarray::simd::CrossPowerSums`]). The bivariate member of the /// [`Terminal::GroupMomentsI32`] family: same key addresses, same drops, /// same one pass with both lanes read in place, same `2^32`-row /// exactness bound ([`MASKED_SUM_I32_MAX_ROWS`] — the `Σx`/`Σy` fields diff --git a/crates/lance-graph-mask-risc/src/lib.rs b/crates/lance-graph-mask-risc/src/lib.rs index ba9454285..bf90bb0bf 100644 --- a/crates/lance-graph-mask-risc/src/lib.rs +++ b/crates/lance-graph-mask-risc/src/lib.rs @@ -138,7 +138,7 @@ pub use reference::{ pub use ternlog_dispatch::{ ternlog_any_dispatch, ternlog_dispatch, ternlog_dispatch_assign, ternlog_popcount_dispatch, }; -pub use value::{ExecError, GroupCrossMoments, GroupMoments, LaneKind, Out, Value}; +pub use value::{CrossPowerSums, ExecError, LaneKind, Out, PowerSums, Value}; /// Number of `u64` words a mask over `n_rows` occupies. #[inline] diff --git a/crates/lance-graph-mask-risc/src/reference.rs b/crates/lance-graph-mask-risc/src/reference.rs index b8652760c..df8b4a336 100644 --- a/crates/lance-graph-mask-risc/src/reference.rs +++ b/crates/lance-graph-mask-risc/src/reference.rs @@ -20,8 +20,8 @@ use crate::ir::{ Foreign, GroupFold, GroupKey, LaneRef, MaskOp, Operand, Planes, Pred, Program, Terminal, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, MAX_SCRATCH_SLOTS, }; +use crate::value::{CrossPowerSums, PowerSums}; use crate::value::{ExecError, LaneKind, Out, Value}; -use crate::value::{GroupCrossMoments, GroupMoments}; use crate::words_for; /// The caller's terminal-out, described by SHAPE rather than borrowed — what @@ -1161,9 +1161,9 @@ pub fn reference_execute_into( if let Out::Moments(o) = out { // Independent formulation: seed every slot, then walk the // survivors one row at a time, widening straight to i128 — - // no ndarray kernel and no `GroupMoments::observe` involved. + // no ndarray kernel involved. for x in o.iter_mut() { - *x = GroupMoments::EMPTY; + *x = PowerSums::default(); } let mut sum = vec![0i128; o.len()]; for r in survivors(mask) { @@ -1173,7 +1173,7 @@ pub fn reference_execute_into( let x = i128::from(i32_at(planes, val, r)); o[k].n += 1; sum[k] += x; - o[k].sum_sq += x * x; + o[k].sum_sq += (x * x) as u128; // a square is never negative } for (slot, s) in o.iter_mut().zip(sum) { // In range by the validated row bound; a failure here is @@ -1206,12 +1206,12 @@ pub fn reference_execute_into( let narrow = |v: i128| { i64::try_from(v).expect("validated 2^32-row bound keeps a sum in i64") }; - *slot = GroupCrossMoments { + *slot = CrossPowerSums { n: u64::try_from(a[0]).expect("a count is non-negative"), sum_x: narrow(a[1]), sum_y: narrow(a[2]), - sum_x2: a[3], - sum_y2: a[4], + sum_x_sq: u128::try_from(a[3]).expect("a square sum is non-negative"), + sum_y_sq: u128::try_from(a[4]).expect("a square sum is non-negative"), sum_xy: a[5], }; } diff --git a/crates/lance-graph-mask-risc/src/value.rs b/crates/lance-graph-mask-risc/src/value.rs index 84b603012..a5c01b266 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -7,11 +7,11 @@ use crate::ir::Operand; /// The per-group `(n, Σx, Σy, Σx², Σy², Σxy)` accumulator /// [`Out::CrossMoments`] carries — re-exported for the same reason. -pub use ndarray::simd::GroupCrossMoments; +pub use ndarray::simd::CrossPowerSums; /// The per-group `(n, Σx, Σx²)` accumulator [`Out::Moments`] carries. Re-exported /// here as shared result vocabulary so the independent oracle names the data /// type through this crate, never through the SIMD facade it falsifies. -pub use ndarray::simd::GroupMoments; +pub use ndarray::simd::PowerSums; /// What a [`crate::Program`] produced. #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -40,12 +40,12 @@ pub enum Value { /// an empty group holds. GroupReduced, /// [`crate::Terminal::GroupMomentsI32`]: the caller's `Out::Moments` - /// buffer was written, one [`GroupMoments`] per group; an empty group - /// holds [`GroupMoments::EMPTY`] (`n == 0`). + /// buffer was written, one [`PowerSums`] per group; an empty group + /// holds [`PowerSums::default()`] (`n == 0`). GroupMoments, /// [`crate::Terminal::GroupCrossMomentsI32`]: the caller's - /// `Out::CrossMoments` buffer was written, one [`GroupCrossMoments`] per - /// group; an empty group holds [`GroupCrossMoments::EMPTY`]. + /// `Out::CrossMoments` buffer was written, one [`CrossPowerSums`] per + /// group; an empty group holds [`CrossPowerSums::default()`]. GroupCrossMoments, /// [`crate::Terminal::MaskedStridedGroupSum`]: the widened sum, or `None` /// when it does not fit an `i64` (never a wrapped value). @@ -68,11 +68,11 @@ pub enum Out<'a> { /// long. Mask(&'a mut [u64]), /// [`crate::Terminal::GroupMomentsI32`]'s destination — one - /// [`GroupMoments`] per group, its length IS the group universe `K`. - Moments(&'a mut [GroupMoments]), + /// [`PowerSums`] per group, its length IS the group universe `K`. + Moments(&'a mut [PowerSums]), /// [`crate::Terminal::GroupCrossMomentsI32`]'s destination — one - /// [`GroupCrossMoments`] per group, its length IS the group universe `K`. - CrossMoments(&'a mut [GroupCrossMoments]), + /// [`CrossPowerSums`] per group, its length IS the group universe `K`. + CrossMoments(&'a mut [CrossPowerSums]), } /// The lane width a predicate or terminal expects, for [`ExecError::LaneKind`]. diff --git a/crates/lance-graph-mask-risc/tests/foreign.rs b/crates/lance-graph-mask-risc/tests/foreign.rs index 2ebbb07fa..fee842ee0 100644 --- a/crates/lance-graph-mask-risc/tests/foreign.rs +++ b/crates/lance-graph-mask-risc/tests/foreign.rs @@ -5,9 +5,9 @@ use lance_graph_mask_risc::exec::{execute_extent, execute_into, Scratch}; use lance_graph_mask_risc::reference::{reference_execute_into, reference_scratch_with_foreign}; use lance_graph_mask_risc::{ - scratch_words_for, words_for, ExecError, Foreign, ForeignPlane, GroupCrossMoments, GroupFold, - GroupKey, GroupMoments, LaneKind, LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, - Terminal, Value, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, + scratch_words_for, words_for, CrossPowerSums, ExecError, Foreign, ForeignPlane, GroupFold, + GroupKey, LaneKind, LaneRef, MaskOp, Operand, Out, Planes, PowerSums, Pred, Program, Terminal, + Value, GROUP_SUM_SYM_MAX_ROWS, MASKED_SUM_I32_MAX_ROWS, }; fn lcg(seed: &mut u64) -> u64 { @@ -1633,7 +1633,7 @@ fn group_moments_match_the_oracle_for_every_key_address() { let mut buf = vec![0u64; scratch_words_for(words, slots).expect("sized")]; let mut scratch = Scratch::over(&mut buf, words, slots).expect("carves"); // Dirty sinks on purpose: the terminal must seed them itself. - let dirty = GroupMoments { + let dirty = PowerSums { n: 7, sum: -7, sum_sq: 7, @@ -1690,7 +1690,7 @@ fn group_moments_refuses_malformed_programs() { ); reference_execute_into(&p, &planes, &none, out) }; - let mut sink = vec![GroupMoments::EMPTY; 4]; + let mut sink = vec![PowerSums::default(); 4]; // `val` must be an I32 lane (lane 1 is U32). assert!(matches!( run( @@ -1750,7 +1750,7 @@ fn group_moments_refuses_malformed_programs() { lanes: &lanes, }; let mut s = Scratch::for_program(&p, n).expect("scratch"); - let mut sink = vec![GroupMoments::EMPTY; 4]; + let mut sink = vec![PowerSums::default(); 4]; assert_eq!( execute_extent(&p, &planes, &none, &mut s, Out::Moments(&mut sink), 10..20), Err(ExecError::ExtentUnsupported { @@ -1759,7 +1759,7 @@ fn group_moments_refuses_malformed_programs() { ); assert_eq!( sink, - vec![GroupMoments::EMPTY; 4], + vec![PowerSums::default(); 4], "a refusal writes nothing" ); } @@ -1838,10 +1838,10 @@ fn group_cross_moments_match_the_oracle_for_every_key_address() { let slots = p.scratch_slots as usize; let mut buf = vec![0u64; scratch_words_for(words, slots).expect("sized")]; let mut scratch = Scratch::over(&mut buf, words, slots).expect("carves"); - let dirty = GroupCrossMoments { + let dirty = CrossPowerSums { n: 3, sum_xy: -3, - ..GroupCrossMoments::EMPTY + ..CrossPowerSums::default() }; let mut got_out = vec![dirty; groups]; let got = execute_into( @@ -1899,7 +1899,7 @@ fn group_cross_moments_refuses_malformed_programs() { let run = |t: Terminal, out: Out<'_>| { reference_execute_into(&Program::new(vec![], t), &planes, &none, out) }; - let mut sink = vec![GroupCrossMoments::EMPTY; 4]; + let mut sink = vec![CrossPowerSums::default(); 4]; for (x, y) in [(1u16, 4u16), (2, 1)] { assert!( matches!( @@ -1909,7 +1909,7 @@ fn group_cross_moments_refuses_malformed_programs() { "x={x} y={y}" ); } - let mut univariate = vec![GroupMoments::EMPTY; 4]; + let mut univariate = vec![PowerSums::default(); 4]; for out in [ Out::None, Out::Moments(&mut univariate), @@ -1939,13 +1939,13 @@ fn group_cross_moments_refuses_malformed_programs() { ); assert_eq!( sink, - vec![GroupCrossMoments::EMPTY; 4], + vec![CrossPowerSums::default(); 4], "a refusal writes nothing" ); // x == y is legal and folds the univariate moments. - let mut same = vec![GroupCrossMoments::EMPTY; 4]; + let mut same = vec![CrossPowerSums::default(); 4]; run(term(2, 2), Out::CrossMoments(&mut same)).expect("x == y is legal"); assert!(same .iter() - .all(|g| g.sum_x == g.sum_y && g.sum_xy == g.sum_x2)); + .all(|g| g.sum_x == g.sum_y && u128::try_from(g.sum_xy) == Ok(g.sum_x_sq))); } From b13be86298f0ffb4e9d0565bd59d068de74639dd Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 17:50:16 +0000 Subject: [PATCH 09/14] mask-risc + jc: name the grouped folds power sums, matching ndarray Terminal::GroupMomentsI32 / GroupCrossMomentsI32 become GroupPowerSumsI32 / GroupCrossPowerSumsI32, Out::{Moments, CrossMoments} become Out::{PowerSums, CrossPowerSums}, Value::{GroupMoments, GroupCrossMoments} become Value::{GroupPowerSums, GroupCrossPowerSums}, and the jc finishes become *_from_power_sums / *_from_cross_power_sums. The records were already ndarray's PowerSums / CrossPowerSums; only the names lagged. No behaviour change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- .../2026-09-29-grouped-cross-moments-fold.md | 10 +- ...9-grouped-moments-fold-anova-over-masks.md | 10 +- crates/jc/src/reliability.rs | 2 +- crates/jc/src/stats.rs | 140 ++++++++++-------- crates/lance-graph-mask-risc/src/exec.rs | 22 +-- crates/lance-graph-mask-risc/src/ir.rs | 14 +- crates/lance-graph-mask-risc/src/reference.rs | 54 +++---- crates/lance-graph-mask-risc/src/value.rs | 22 +-- crates/lance-graph-mask-risc/tests/foreign.rs | 75 +++++----- .../tests/duckdb_differential.rs | 6 +- 10 files changed, 190 insertions(+), 165 deletions(-) diff --git a/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md index d22a68051..8d450d1a6 100644 --- a/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md +++ b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md @@ -3,17 +3,17 @@ **Status:** MEASURED · DONE — ndarray `simd_masking_ops.rs`, `crates/lance-graph-mask-risc`, `crates/jc` ## What landed -- ndarray: `GroupCrossMoments { n, sum_x, sum_y: i64, sum_x2, sum_y2, sum_xy: i128 }` + `masked_group_cross_moments_i32{,_via,_pair}` over the same generic `group_walk`. `x_moments()`/`y_moments()` return the exact univariate `GroupMoments`; a test pins them equal to the univariate fold. Bounds: squares ∈ [0, 2^62], products ∈ [−2^62+2^31, 2^62]; the `i64` sums bind at 2^32 rows/group; the `i128` sums cannot overflow at any `u64` row count. -- mask-risc: `Terminal::GroupCrossMomentsI32 { mask, key, x, y }` → `Out::CrossMoments`. -- jc: `pearson_from_cross_moments`, `sample_covariance_from_cross_moments` (n−1), `simple_regression_from_cross_moments` (`SimpleRegression { slope, intercept }`), `r_squared_from_cross_moments` (`multiple_r_squared`'s k=1 contract). Shared tails extracted: `pearson_from_centered`, `sample_cov_tail`, `r_squared_tail`. Centred sums formed exactly in `i128`, checked. +- ndarray: `GroupCrossPowerSums { n, sum_x, sum_y: i64, sum_x2, sum_y2, sum_xy: i128 }` + `masked_group_cross_moments_i32{,_via,_pair}` over the same generic `group_walk`. `x_moments()`/`y_moments()` return the exact univariate `GroupPowerSums`; a test pins them equal to the univariate fold. Bounds: squares ∈ [0, 2^62], products ∈ [−2^62+2^31, 2^62]; the `i64` sums bind at 2^32 rows/group; the `i128` sums cannot overflow at any `u64` row count. +- mask-risc: `Terminal::GroupCrossPowerSumsI32 { mask, key, x, y }` → `Out::CrossPowerSums`. +- jc: `pearson_from_cross_power_sums`, `sample_covariance_from_cross_power_sums` (n−1), `simple_regression_from_cross_power_sums` (`SimpleRegression { slope, intercept }`), `r_squared_from_cross_power_sums` (`multiple_r_squared`'s k=1 contract). Shared tails extracted: `pearson_from_centered`, `sample_cov_tail`, `r_squared_tail`. Centred sums formed exactly in `i128`, checked. ## Measured 48 layouts (n ≤ 100k): |Δr| ≤ 1.3e-14, |ΔR²| ≤ 4.1e-14, relative Δslope ≤ 2.8e-14 vs the materialized path. Huge offset + tiny spread (x ≈ ±2·10⁹): fold path within 1e-15 of an independent exact reference. **Unlike ANOVA, the slice path does not fail here** — `pearson`/`multiple_r_squared` are two-pass. A naive f64 projection of the same moments returns NaN; that is what exact `i128` centring prevents (disable run: 2 tests fail). ## Redundant guards (disable runs, kept as explicit contract) -`pearson_from_cross_moments`'s `n < 2`, the regression's `cxx == 0`, and R²'s `cyy == 0` are each subsumed by the shared tail's zero/non-finite rejection. The first past-bound refusal fixture was vacuous (wrapped to negative centred squares); replaced with one that wraps to a plausible `r = 1`. +`pearson_from_cross_power_sums`'s `n < 2`, the regression's `cxx == 0`, and R²'s `cyy == 0` are each subsumed by the shared tail's zero/non-finite rejection. The first past-bound refusal fixture was vacuous (wrapped to negative centred squares); replaced with one that wraps to a plausible `r = 1`. ## Open - Multi-membership in one pass: `lance-graph-report` lowers `CoordSpec::MaskSet` to one `Filter::Plane` pass per member tuple and per fold state. The single-group limit is in `GroupKeyAddr::group_of` (one `Option` per row), in the IR's `GroupKey`, and in report lowering. Seam: a word-level walker `for each 64-row word: sel & plane_m → fold hits into out[m]` over existing planes — same plane traffic as K passes, value-lane traffic of the union instead of the sum. -**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupMoments`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. +**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupPowerSums`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. diff --git a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md index 5c2b0f749..c32426630 100644 --- a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md +++ b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md @@ -3,16 +3,16 @@ **Status:** MEASURED · DONE — ndarray `simd_masking_ops.rs`, `crates/lance-graph-mask-risc`, `crates/jc/src/stats.rs` ## What landed -- ndarray: `GroupMoments { n: u64, sum: i64, sum_sq: i128 }` and `masked_group_moments_i32{,_via,_pair}` — `(n, Σx, Σx²)` per group in one pass over the rows a mask selects, over the existing `group_walk` (now generic over its slot type). `checked_merge` is exact integer addition; exact up to `GROUP_MOMENTS_MAX_ROWS = 2^32` rows per group (the `i64` sum is the binding field). -- mask-risc: `Terminal::GroupMomentsI32 { mask, key: GroupKey, val }` → `Out::Moments`. Own terminal, not a `GroupFold` member (a `GroupFold` slot is one seeded `i64`). Refuses planes past `MASKED_SUM_I32_MAX_ROWS` and partial extents. -- jc: `anova_from_moments` / `eta_squared_from_moments`. F/p/η² policy is shared with `anova_one_way` / `eta_squared` via `anova_from_ss` / `eta_from_ss`. Sums of squares come from exact `i128` quantities; no two large floats are subtracted. +- ndarray: `GroupPowerSums { n: u64, sum: i64, sum_sq: i128 }` and `masked_group_moments_i32{,_via,_pair}` — `(n, Σx, Σx²)` per group in one pass over the rows a mask selects, over the existing `group_walk` (now generic over its slot type). `checked_merge` is exact integer addition; exact up to `GROUP_MOMENTS_MAX_ROWS = 2^32` rows per group (the `i64` sum is the binding field). +- mask-risc: `Terminal::GroupPowerSumsI32 { mask, key: GroupKey, val }` → `Out::PowerSums`. Own terminal, not a `GroupFold` member (a `GroupFold` slot is one seeded `i64`). Refuses planes past `MASKED_SUM_I32_MAX_ROWS` and partial extents. +- jc: `anova_from_power_sums` / `eta_squared_from_power_sums`. F/p/η² policy is shared with `anova_one_way` / `eta_squared` via `anova_from_ss` / `eta_from_ss`. Sums of squares come from exact `i128` quantities; no two large floats are subtracted. ## Measured -Mask → fold → `anova_from_moments` vs materialized → `anova_one_way`, 42 non-degenerate layouts (n ≤ 100k, k ≤ 6, densities 20/128/256 of 256): max relative ΔF 4.2e-14, max |Δp| 4.4e-16, max |Δη²| 6.5e-16. Adversarial layout (group means ±2·10⁹, within spread 1): exact F = 3·10¹⁹; the moments path returns 3·10¹⁹, the slice path returns `None` (its `ss_t − ss_b` cancels the within-group SS to ≤ 0). A naive float `Q − S²/N` inside the moments path fails even the near-`i32::MAX` agreement test (ΔF ≈ 1e-3), so the exact forms are load-bearing. +Mask → fold → `anova_from_power_sums` vs materialized → `anova_one_way`, 42 non-degenerate layouts (n ≤ 100k, k ≤ 6, densities 20/128/256 of 256): max relative ΔF 4.2e-14, max |Δp| 4.4e-16, max |Δη²| 6.5e-16. Adversarial layout (group means ±2·10⁹, within spread 1): exact F = 3·10¹⁹; the moments path returns 3·10¹⁹, the slice path returns `None` (its `ss_t − ss_b` cancels the within-group SS to ≤ 0). A naive float `Q − S²/N` inside the moments path fails even the near-`i32::MAX` agreement test (ΔF ≈ 1e-3), so the exact forms are load-bearing. ## Open - Multi-membership is NOT one pass: `GroupKey` resolves one group per row, and `lance-graph-report`'s `CoordSpec::MaskSet` executes one `Filter::Plane` program per member. Seam: a key address that walks the selected rows once and folds each row into every member plane that holds it. - Bivariate `(n, Σx, Σy, Σx², Σy², Σxy)`: same walk, a second value lane in the closure, a wider slot type; not built. - No early exit, and no `CausalEdge64` commit. ANOVA is non-monotone under future rows. -**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupMoments`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. +**Renamed 2026-09-30 (rebase onto ndarray master):** upstream shipped `PowerSums { n, sum, sum_sq: u128 }` — the same record as `GroupPowerSums`. The duplicate was dropped: consumers use `PowerSums` / `CrossPowerSums` and `masked_group_(cross_)power_sums_i32*`; `checked_merge` moved onto the upstream types. Square sums are now `u128`; jc converts with `i128::try_from`, and a value past `i128::MAX` is refused as past the bound. diff --git a/crates/jc/src/reliability.rs b/crates/jc/src/reliability.rs index 34982499f..9f6b1de5c 100644 --- a/crates/jc/src/reliability.rs +++ b/crates/jc/src/reliability.rs @@ -115,7 +115,7 @@ pub fn pearson(x: &[f64], y: &[f64]) -> Option { /// Pearson's `r` from the centred co-moment and the two centred /// second moments — the ONE place its degeneracy policy lives, shared by -/// [`pearson`] and `stats::pearson_from_cross_moments`. Any common positive +/// [`pearson`] and `stats::pearson_from_cross_power_sums`. Any common positive /// scale on all three (e.g. `n·S` instead of `S`) cancels. pub(crate) fn pearson_from_centered(sxy: f64, sxx: f64, syy: f64) -> Option { let denom = (sxx * syy).sqrt(); diff --git a/crates/jc/src/stats.rs b/crates/jc/src/stats.rs index 738e8b454..39eec4a72 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -114,7 +114,7 @@ fn sample_cov(x: &[f64], y: &[f64]) -> Option { /// The unbiased (divisor `n−1`) covariance from the centred co-moment /// `Σ(x−x̄)(y−ȳ)` — the ONE place the convention lives, shared by -/// `sample_cov` and [`sample_covariance_from_cross_moments`]. +/// `sample_cov` and [`sample_covariance_from_cross_power_sums`]. #[inline] fn sample_cov_tail(sxy: f64, n: usize) -> Option { (n >= 2).then(|| sxy / (n as f64 - 1.0)) @@ -888,7 +888,7 @@ pub fn multiple_r_squared(y: &[f64], predictors: &[Vec]) -> Option { } /// The acceptance rule for a computed `R²` — the ONE place it lives, shared -/// by [`multiple_r_squared`] and [`r_squared_from_cross_moments`]. +/// by [`multiple_r_squared`] and [`r_squared_from_cross_power_sums`]. /// /// The only legitimate excursions are rounding-scale: `ss_res` marginally /// exceeding `ss_tot` (R² slightly below 0) or cancellation pushing it a few @@ -963,7 +963,7 @@ pub fn eta_squared(groups: &[Vec]) -> Option { } /// η² from the two sums of squares — the ONE place its degeneracy policy -/// lives, shared by [`eta_squared`] and [`eta_squared_from_moments`]. +/// lives, shared by [`eta_squared`] and [`eta_squared_from_power_sums`]. fn eta_from_ss(ss_b: f64, ss_t: f64) -> Option { if ss_t == 0.0 || !ss_t.is_finite() { return None; @@ -1157,7 +1157,7 @@ pub fn anova_one_way(groups: &[Vec]) -> Option { /// The F test from the three sums of squares — the ONE place the one-way /// ANOVA's degeneracy policy lives, shared by [`anova_one_way`] (which -/// passes `ss_w = ss_t − ss_b`) and [`anova_from_moments`] (which forms +/// passes `ss_w = ss_t − ss_b`) and [`anova_from_power_sums`] (which forms /// `ss_w` exactly and `ss_t = ss_b + ss_w`). fn anova_from_ss(ss_b: f64, ss_w: f64, ss_t: f64, k: usize, n_total: usize) -> Option { if n_total <= k { @@ -1220,7 +1220,7 @@ fn anova_from_ss(ss_b: f64, ss_w: f64, ss_t: f64, k: usize, n_total: usize) -> O /// [`anova_one_way`]'s empty-group rule), or moments too large for the exact /// `i128` forms — i.e. outside the `2^32`-rows-per-group bound /// `PowerSums` is exact under. -fn one_way_ss_from_moments(groups: &[PowerSums]) -> Option<(f64, f64, f64, usize, usize)> { +fn one_way_ss_from_power_sums(groups: &[PowerSums]) -> Option<(f64, f64, f64, usize, usize)> { if groups.len() < 2 || groups.iter().any(|g| g.n == 0) { return None; } @@ -1260,15 +1260,15 @@ fn one_way_ss_from_moments(groups: &[PowerSums]) -> Option<(f64, f64, f64, usize /// /// `groups[g]` holds group `g`'s `(n, Σx, Σx²)`, e.g. as folded out of a /// population mask by `ndarray::simd::masked_group_power_sums_i32` (or -/// `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32`). Every degenerate +/// `lance-graph-mask-risc`'s `Terminal::GroupPowerSumsI32`). Every degenerate /// case is [`anova_one_way`]'s: fewer than 2 groups, an empty group, /// `N ≤ k`, and zero within-group variance are all `None`. The sums of -/// squares are formed exactly (see `one_way_ss_from_moments`), so on data +/// squares are formed exactly (see `one_way_ss_from_power_sums`), so on data /// where the slice path's two-pass floats cancel badly this is the more /// accurate of the two, not merely an approximation of it. /// /// ``` -/// use jc::stats::{anova_from_moments, anova_one_way}; +/// use jc::stats::{anova_from_power_sums, anova_one_way}; /// use ndarray::simd::PowerSums; /// /// let groups = vec![vec![1.0, 2.0, 3.0], vec![4.0, 5.0, 6.0]]; @@ -1276,11 +1276,11 @@ fn one_way_ss_from_moments(groups: &[PowerSums]) -> Option<(f64, f64, f64, usize /// PowerSums { n: 3, sum: 6, sum_sq: 14 }, /// PowerSums { n: 3, sum: 15, sum_sq: 77 }, /// ]; -/// let (a, b) = (anova_one_way(&groups).unwrap(), anova_from_moments(&moments).unwrap()); +/// let (a, b) = (anova_one_way(&groups).unwrap(), anova_from_power_sums(&moments).unwrap()); /// assert!((a.f - b.f).abs() < 1e-12 && (a.p - b.p).abs() < 1e-12); /// ``` -pub fn anova_from_moments(groups: &[PowerSums]) -> Option { - let (ss_b, ss_w, ss_t, k, n_total) = one_way_ss_from_moments(groups)?; +pub fn anova_from_power_sums(groups: &[PowerSums]) -> Option { + let (ss_b, ss_w, ss_t, k, n_total) = one_way_ss_from_power_sums(groups)?; anova_from_ss(ss_b, ss_w, ss_t, k, n_total) } @@ -1288,8 +1288,8 @@ pub fn anova_from_moments(groups: &[PowerSums]) -> Option { /// degeneracy policy (`None` only on fewer than 2 groups, an empty group, or /// zero total variance; a perfectly separated layout with zero within-group /// variance is a real `1.0`, unlike the F test). -pub fn eta_squared_from_moments(groups: &[PowerSums]) -> Option { - let (ss_b, _, ss_t, _, _) = one_way_ss_from_moments(groups)?; +pub fn eta_squared_from_power_sums(groups: &[PowerSums]) -> Option { + let (ss_b, _, ss_t, _, _) = one_way_ss_from_power_sums(groups)?; eta_from_ss(ss_b, ss_t) } @@ -1333,13 +1333,13 @@ fn centered_cross(m: &CrossPowerSums) -> Option<(i128, i128, i128)> { /// degeneracy policy (`None` for `n < 2` or a constant `x` or `y`). /// /// ``` -/// use jc::stats::pearson_from_cross_moments; +/// use jc::stats::pearson_from_cross_power_sums; /// use ndarray::simd::CrossPowerSums; /// // x = [1,2,3], y = [2,4,6]: perfectly correlated. /// let m = CrossPowerSums { n: 3, sum_x: 6, sum_y: 12, sum_x_sq: 14, sum_y_sq: 56, sum_xy: 28 }; -/// assert!((pearson_from_cross_moments(&m).unwrap() - 1.0).abs() < 1e-12); +/// assert!((pearson_from_cross_power_sums(&m).unwrap() - 1.0).abs() < 1e-12); /// ``` -pub fn pearson_from_cross_moments(m: &CrossPowerSums) -> Option { +pub fn pearson_from_cross_power_sums(m: &CrossPowerSums) -> Option { if m.n < 2 { return None; } @@ -1349,7 +1349,7 @@ pub fn pearson_from_cross_moments(m: &CrossPowerSums) -> Option { /// The unbiased (divisor `n−1`) sample covariance from grouped cross moments /// — the convention this module already uses (`None` for `n < 2`). -pub fn sample_covariance_from_cross_moments(m: &CrossPowerSums) -> Option { +pub fn sample_covariance_from_cross_power_sums(m: &CrossPowerSums) -> Option { let (cxy, _, _) = centered_cross(m)?; // `C_xy = n·Σ(x−x̄)(y−ȳ)`: divide the `n` back out before the tail. let sxy = cxy as f64 / m.n.max(1) as f64; @@ -1377,7 +1377,7 @@ pub struct SimpleRegression { /// `|β·x̄|` is large and `α` is small its ABSOLUTE error is on the order of /// `ε·(|β|·|x̄| + |ȳ|)`. That is the conditioning of the quantity itself, not /// of the fold; centre `x` before regressing when `α` itself matters. -pub fn simple_regression_from_cross_moments(m: &CrossPowerSums) -> Option { +pub fn simple_regression_from_cross_power_sums(m: &CrossPowerSums) -> Option { if m.n < 2 { return None; } @@ -1395,7 +1395,7 @@ pub fn simple_regression_from_cross_moments(m: &CrossPowerSums) -> Option Option { +pub fn r_squared_from_cross_power_sums(m: &CrossPowerSums) -> Option { if m.n < 3 { return None; } @@ -2494,7 +2494,7 @@ mod tests { // statistics folded straight off a population mask equals the same // statistic computed from the materialized groups. The fold is // `ndarray::simd::masked_group_power_sums_i32` — the kernel - // `lance-graph-mask-risc`'s `Terminal::GroupMomentsI32` delegates to. + // `lance-graph-mask-risc`'s `Terminal::GroupPowerSumsI32` delegates to. // // Tolerances. Both paths form the same sums of squares by different // float orderings (two-pass deviations vs exact integers then one @@ -2506,7 +2506,7 @@ mod tests { // 1e-9 absolute on p leave three orders of margin and still fail on any // real algebraic difference (a dropped row, a wrong df, a wrong group). - mod moments_equivalence { + mod power_sums_equivalence { use super::super::*; use ndarray::simd::masked_group_power_sums_i32; @@ -2616,8 +2616,9 @@ mod tests { let moments = fold(&p, k); let what = format!("n={n} k={k} density={density}"); let a = anova_one_way(&groups); - assert_same(a, anova_from_moments(&moments), &what); - let (ea, eb) = (eta_squared(&groups), eta_squared_from_moments(&moments)); + assert_same(a, anova_from_power_sums(&moments), &what); + let (ea, eb) = + (eta_squared(&groups), eta_squared_from_power_sums(&moments)); match (ea, eb) { (Some(x), Some(y)) => assert!((x - y).abs() < 1e-9, "{what}: η²"), (x, y) => assert_eq!(x, y, "{what}: η² degeneracy"), @@ -2659,9 +2660,12 @@ mod tests { } } assert_eq!(merged, whole, "chunks={chunks}"); - assert_eq!(anova_from_moments(&merged), anova_from_moments(&whole)); + assert_eq!( + anova_from_power_sums(&merged), + anova_from_power_sums(&whole) + ); } - assert!(anova_from_moments(&whole).is_some()); + assert!(anova_from_power_sums(&whole).is_some()); } /// FAILS IF: values near the i32 bound lose precision on the @@ -2675,7 +2679,11 @@ mod tests { let groups = materialize(&p, k); let a = anova_one_way(&groups); assert!(a.is_some(), "fixture must be non-degenerate"); - assert_same(a, anova_from_moments(&fold(&p, k)), &format!("base={base}")); + assert_same( + a, + anova_from_power_sums(&fold(&p, k)), + &format!("base={base}"), + ); } } @@ -2716,16 +2724,16 @@ mod tests { ]; for (what, groups, moments) in cases { assert_eq!(anova_one_way(&groups), None, "{what}: slice"); - assert_eq!(anova_from_moments(&moments), None, "{what}: moments"); + assert_eq!(anova_from_power_sums(&moments), None, "{what}: moments"); assert_eq!( eta_squared(&groups), - eta_squared_from_moments(&moments), + eta_squared_from_power_sums(&moments), "{what}: η²" ); } // The one degenerate case where η² is NOT None: perfect separation. let sep = [m(3, &[7, 7, 7]), m(2, &[9, 9])]; - assert_eq!(eta_squared_from_moments(&sep), Some(1.0)); + assert_eq!(eta_squared_from_power_sums(&sep), Some(1.0)); } /// FAILS IF: the moments path loses the within-group variance to @@ -2759,7 +2767,7 @@ mod tests { .map(|&b| 6 * i128::from(b) * i128::from(b)) .sum(); let f_exact = (ss_b as f64 / 2.0) / (ss_w as f64 / 15.0); - let got = anova_from_moments(&moments).expect("non-degenerate"); + let got = anova_from_power_sums(&moments).expect("non-degenerate"); assert_eq!((got.df_between, got.df_within), (2.0, 15.0)); assert!( rel(got.f, f_exact) < 1e-12, @@ -2784,9 +2792,9 @@ mod tests { // covariance, R² and slope leaves wide margin while still catching any // real algebraic error. The intercept is compared with the ABSOLUTE // bound its own conditioning allows (see - // `simple_regression_from_cross_moments`). + // `simple_regression_from_cross_power_sums`). - mod cross_moments_equivalence { + mod cross_power_sums_equivalence { use super::super::*; use crate::reliability::pearson; use ndarray::simd::masked_group_cross_power_sums_i32; @@ -2882,11 +2890,11 @@ mod tests { /// Compare every projection of one group against its slice /// counterpart; both `None` or both within tolerance. fn assert_group_agrees(x: &[f64], y: &[f64], m: &CrossPowerSums, what: &str) { - match (pearson(x, y), pearson_from_cross_moments(m)) { + match (pearson(x, y), pearson_from_cross_power_sums(m)) { (Some(a), Some(b)) => assert!(rel(a, b) < 1e-9, "{what}: r {a} vs {b}"), (a, b) => assert_eq!(a, b, "{what}: r degeneracy"), } - match (sample_cov(x, y), sample_covariance_from_cross_moments(m)) { + match (sample_cov(x, y), sample_covariance_from_cross_power_sums(m)) { (Some(a), Some(b)) => { assert!( (a - b).abs() <= 1e-9 * a.abs().max(1.0), @@ -2897,12 +2905,12 @@ mod tests { } match ( multiple_r_squared(y, &[x.to_vec()]), - r_squared_from_cross_moments(m), + r_squared_from_cross_power_sums(m), ) { (Some(a), Some(b)) => assert!((a - b).abs() < 1e-9, "{what}: R² {a} vs {b}"), (a, b) => assert_eq!(a, b, "{what}: R² degeneracy"), } - match (ols(x, y), simple_regression_from_cross_moments(m)) { + match (ols(x, y), simple_regression_from_cross_power_sums(m)) { (Some((sl, ic)), Some(r)) => { assert!(rel(sl, r.slope) < 1e-9, "{what}: slope {sl} vs {}", r.slope); let mx = mean(x).unwrap().abs(); @@ -2993,17 +3001,20 @@ mod tests { } assert_eq!(merged, whole, "chunks={chunks} reversed={order}"); for (a, b) in merged.iter().zip(&whole) { - assert_eq!(pearson_from_cross_moments(a), pearson_from_cross_moments(b)); assert_eq!( - simple_regression_from_cross_moments(a), - simple_regression_from_cross_moments(b) + pearson_from_cross_power_sums(a), + pearson_from_cross_power_sums(b) + ); + assert_eq!( + simple_regression_from_cross_power_sums(a), + simple_regression_from_cross_power_sums(b) ); } } } } - fn moments_of(xs: &[i64], ys: &[i64]) -> CrossPowerSums { + fn power_sums_of(xs: &[i64], ys: &[i64]) -> CrossPowerSums { // Through the real fold: one group, every row selected. let xs: Vec = xs.iter().map(|&x| x as i32).collect(); let ys: Vec = ys.iter().map(|&y| y as i32).collect(); @@ -3034,28 +3045,35 @@ mod tests { ("perfect -1", &[-3, 1, 4, 10], &[7, -1, -7, -19]), ]; for (what, xs, ys) in cases { - let m = moments_of(xs, ys); + let m = power_sums_of(xs, ys); let xf: Vec = xs.iter().map(|&v| v as f64).collect(); let yf: Vec = ys.iter().map(|&v| v as f64).collect(); assert_group_agrees(&xf, &yf, &m, what); } // The contract, stated positively rather than only as agreement. - let m = |x: &[i64], y: &[i64]| moments_of(x, y); - assert_eq!(pearson_from_cross_moments(&m(&[], &[])), None); - assert_eq!(sample_covariance_from_cross_moments(&m(&[4], &[9])), None); - assert_eq!(r_squared_from_cross_moments(&m(&[1, 3], &[2, 7])), None); - assert!(pearson_from_cross_moments(&m(&[1, 3], &[2, 7])).is_some()); - assert_eq!(pearson_from_cross_moments(&m(&[5, 5, 5], &[1, 2, 3])), None); + let m = |x: &[i64], y: &[i64]| power_sums_of(x, y); + assert_eq!(pearson_from_cross_power_sums(&m(&[], &[])), None); + assert_eq!( + sample_covariance_from_cross_power_sums(&m(&[4], &[9])), + None + ); + assert_eq!(r_squared_from_cross_power_sums(&m(&[1, 3], &[2, 7])), None); + assert!(pearson_from_cross_power_sums(&m(&[1, 3], &[2, 7])).is_some()); + assert_eq!( + pearson_from_cross_power_sums(&m(&[5, 5, 5], &[1, 2, 3])), + None + ); assert_eq!( - simple_regression_from_cross_moments(&m(&[5, 5, 5], &[1, 2, 3])), + simple_regression_from_cross_power_sums(&m(&[5, 5, 5], &[1, 2, 3])), None ); - let flat = simple_regression_from_cross_moments(&m(&[1, 2, 3], &[5, 5, 5])).unwrap(); + let flat = simple_regression_from_cross_power_sums(&m(&[1, 2, 3], &[5, 5, 5])).unwrap(); assert_eq!((flat.slope, flat.intercept), (0.0, 5.0)); - let r = pearson_from_cross_moments(&m(&[-3, 1, 4, 10], &[7, -1, -7, -19])).unwrap(); + let r = pearson_from_cross_power_sums(&m(&[-3, 1, 4, 10], &[7, -1, -7, -19])).unwrap(); assert!((r + 1.0).abs() < 1e-15, "perfect -1: {r}"); let line = - simple_regression_from_cross_moments(&m(&[-3, 1, 4, 10], &[-5, 3, 9, 21])).unwrap(); + simple_regression_from_cross_power_sums(&m(&[-3, 1, 4, 10], &[-5, 3, 9, 21])) + .unwrap(); assert_eq!((line.slope, line.intercept), (2.0, 1.0)); } @@ -3104,7 +3122,7 @@ mod tests { ] { let xs: Vec = d.iter().map(|v| b + v).collect(); let ys: Vec = d.iter().zip(e).map(|(v, w)| c + v + w).collect(); - let m = moments_of(&xs, &ys); + let m = power_sums_of(&xs, &ys); let sxy: i64 = d.iter().zip(e).map(|(v, w)| v * (v + w)).sum(); let sxx: i64 = d.iter().map(|v| v * v).sum(); let syy: i64 = d.iter().zip(e).map(|(v, w)| (v + w) * (v + w)).sum(); @@ -3116,16 +3134,16 @@ mod tests { / sxx as f64; let cov_exact = sxy as f64 / 9.0; - let r = pearson_from_cross_moments(&m).unwrap(); + let r = pearson_from_cross_power_sums(&m).unwrap(); assert!(rel(r, r_exact) < 1e-15, "B={b}: r {r} vs exact {r_exact}"); - let r2 = r_squared_from_cross_moments(&m).unwrap(); + let r2 = r_squared_from_cross_power_sums(&m).unwrap(); assert!(rel(r2, r_exact * r_exact) < 1e-15, "B={b}: R²"); - let cov = sample_covariance_from_cross_moments(&m).unwrap(); + let cov = sample_covariance_from_cross_power_sums(&m).unwrap(); assert!( rel(cov, cov_exact) < 1e-15, "B={b}: cov {cov} vs {cov_exact}" ); - let line = simple_regression_from_cross_moments(&m).unwrap(); + let line = simple_regression_from_cross_power_sums(&m).unwrap(); assert!(rel(line.slope, slope_exact) < 1e-15, "B={b}: slope"); // The documented intercept bound: ε·(|β|·|x̄| + |ȳ|). let bound = @@ -3165,7 +3183,7 @@ mod tests { /// check for the wrong reason, so the test passed with every checked /// operation disabled. Measured by a disable run.) #[test] - fn moments_past_the_bound_are_refused() { + fn power_sums_past_the_bound_are_refused() { let v = (1i128 << 95) + 1; let huge = CrossPowerSums { n: 1 << 33, @@ -3175,10 +3193,10 @@ mod tests { sum_y_sq: v as u128, sum_xy: v, }; - assert_eq!(pearson_from_cross_moments(&huge), None); - assert_eq!(sample_covariance_from_cross_moments(&huge), None); - assert_eq!(simple_regression_from_cross_moments(&huge), None); - assert_eq!(r_squared_from_cross_moments(&huge), None); + assert_eq!(pearson_from_cross_power_sums(&huge), None); + assert_eq!(sample_covariance_from_cross_power_sums(&huge), None); + assert_eq!(simple_regression_from_cross_power_sums(&huge), None); + assert_eq!(r_squared_from_cross_power_sums(&huge), None); } } } diff --git a/crates/lance-graph-mask-risc/src/exec.rs b/crates/lance-graph-mask-risc/src/exec.rs index e951ee2b7..51d2311c6 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -1422,8 +1422,8 @@ fn precheck( Terminal::GroupSumI32 { .. } => Some("GroupSumI32"), Terminal::GroupSumViaI32 { .. } => Some("GroupSumViaI32"), Terminal::GroupReduce { .. } => Some("GroupReduce"), - Terminal::GroupMomentsI32 { .. } => Some("GroupMomentsI32"), - Terminal::GroupCrossMomentsI32 { .. } => Some("GroupCrossMomentsI32"), + Terminal::GroupPowerSumsI32 { .. } => Some("GroupPowerSumsI32"), + Terminal::GroupCrossPowerSumsI32 { .. } => Some("GroupCrossPowerSumsI32"), }; if let Some(what) = refused { return Err(ExecError::ExtentUnsupported { what }); @@ -1577,8 +1577,8 @@ pub fn execute_compiled( } (Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. }, Out::I64(o)) => o.fill(0), (Terminal::GroupReduce { fold, .. }, Out::I64(o)) => o.fill(fold.seed()), - (Terminal::GroupMomentsI32 { .. }, Out::Moments(o)) => o.fill(PowerSums::default()), - (Terminal::GroupCrossMomentsI32 { .. }, Out::CrossMoments(o)) => { + (Terminal::GroupPowerSumsI32 { .. }, Out::PowerSums(o)) => o.fill(PowerSums::default()), + (Terminal::GroupCrossPowerSumsI32 { .. }, Out::CrossPowerSums(o)) => { o.fill(CrossPowerSums::default()) } _ => {} @@ -1992,12 +1992,12 @@ pub fn execute_compiled( } } } - Terminal::GroupMomentsI32 { mask, key, val } => { + Terminal::GroupPowerSumsI32 { mask, key, val } => { // `validate` already refused a missing/too-small `out`, every // wrong-width lane and a plane past the 2^32-row exactness // bound; one delegation per tile (law L3) into the sink // seeded above with `PowerSums::default()`. - if let Out::Moments(o) = &mut out { + if let Out::PowerSums(o) = &mut out { let m = read(planes, &slots, mask, t); let v = lane_i32(planes, val, t); match key { @@ -2022,10 +2022,10 @@ pub fn execute_compiled( } } } - Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { - // Same contract as GroupMomentsI32, both lanes read in place: + Terminal::GroupCrossPowerSumsI32 { mask, key, x, y } => { + // Same contract as GroupPowerSumsI32, both lanes read in place: // one delegation per tile into the sink seeded above. - if let Out::CrossMoments(o) = &mut out { + if let Out::CrossPowerSums(o) = &mut out { let m = read(planes, &slots, mask, t); let (xs, ys) = (lane_i32(planes, x, t), lane_i32(planes, y, t)); match key { @@ -2088,8 +2088,8 @@ pub fn execute_compiled( Terminal::CountKeyRunsU32 { .. } => Value::Count(runs + run_carry.finish()), Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. } => Value::GroupSummed, Terminal::GroupReduce { .. } => Value::GroupReduced, - Terminal::GroupMomentsI32 { .. } => Value::GroupMoments, - Terminal::GroupCrossMomentsI32 { .. } => Value::GroupCrossMoments, + Terminal::GroupPowerSumsI32 { .. } => Value::GroupPowerSums, + Terminal::GroupCrossPowerSumsI32 { .. } => Value::GroupCrossPowerSums, Terminal::Keep { mask } => Value::Mask(mask), }) } diff --git a/crates/lance-graph-mask-risc/src/ir.rs b/crates/lance-graph-mask-risc/src/ir.rs index f33fb1db8..0745cf2f4 100644 --- a/crates/lance-graph-mask-risc/src/ir.rs +++ b/crates/lance-graph-mask-risc/src/ir.rs @@ -423,7 +423,7 @@ pub enum Terminal { }, /// Grouped sufficient statistics of an `I32` lane: for every row `i` /// where `mask` holds, resolves the row's group through [`GroupKey`] and - /// folds `lanes[val][i]` into the caller's `Out::Moments` buffer as + /// folds `lanes[val][i]` into the caller's `Out::PowerSums` buffer as /// `n += 1, Σx += x, Σx² += x²` ([`ndarray::simd::PowerSums`]). One /// pass over the selected rows; no selected value is copied out. The /// buffer's length IS the group universe `K`, and the same zero-fallback @@ -439,22 +439,22 @@ pub enum Terminal { /// and no group can exceed the plane. The executor refuses a wider plane /// rather than wrap. The sink is seeded with `PowerSums::default()` before /// the first tile, so an empty group reads `n == 0`. - GroupMomentsI32 { + GroupPowerSumsI32 { mask: Operand, key: GroupKey, val: u16, }, /// Grouped cross moments of TWO `I32` lanes: for every row `i` where /// `mask` holds, folds `(lanes[x][i], lanes[y][i])` into the caller's - /// `Out::CrossMoments` buffer as `n, Σx, Σy, Σx², Σy², Σxy` + /// `Out::CrossPowerSums` buffer as `n, Σx, Σy, Σx², Σy², Σxy` /// ([`ndarray::simd::CrossPowerSums`]). The bivariate member of the - /// [`Terminal::GroupMomentsI32`] family: same key addresses, same drops, + /// [`Terminal::GroupPowerSumsI32`] family: same key addresses, same drops, /// same one pass with both lanes read in place, same `2^32`-row /// exactness bound ([`MASKED_SUM_I32_MAX_ROWS`] — the `Σx`/`Σy` fields /// are `i64`). A physical fold: covariance, correlation and simple /// regression are consumer projections. `x == y` is legal (it folds the /// univariate moments twice over). - GroupCrossMomentsI32 { + GroupCrossPowerSumsI32 { mask: Operand, key: GroupKey, x: u16, @@ -651,8 +651,8 @@ impl Program { | Terminal::GroupSumI32 { mask, .. } | Terminal::GroupSumViaI32 { mask, .. } | Terminal::GroupReduce { mask, .. } - | Terminal::GroupMomentsI32 { mask, .. } - | Terminal::GroupCrossMomentsI32 { mask, .. } + | Terminal::GroupPowerSumsI32 { mask, .. } + | Terminal::GroupCrossPowerSumsI32 { mask, .. } | Terminal::Keep { mask } => touch(mask), } Self { diff --git a/crates/lance-graph-mask-risc/src/reference.rs b/crates/lance-graph-mask-risc/src/reference.rs index df8b4a336..394b04d29 100644 --- a/crates/lance-graph-mask-risc/src/reference.rs +++ b/crates/lance-graph-mask-risc/src/reference.rs @@ -35,8 +35,8 @@ pub(crate) enum OutShape { I32(usize), I64(usize), Mask(usize), - Moments(usize), - CrossMoments(usize), + PowerSums(usize), + CrossPowerSums(usize), } /// [`OutShape`] of a borrowed `out` — the caller keeps `out` itself to write @@ -47,8 +47,8 @@ pub(crate) fn out_shape(out: &Out<'_>) -> OutShape { Out::I32(v) => OutShape::I32(v.len()), Out::I64(v) => OutShape::I64(v.len()), Out::Mask(v) => OutShape::Mask(v.len()), - Out::Moments(v) => OutShape::Moments(v.len()), - Out::CrossMoments(v) => OutShape::CrossMoments(v.len()), + Out::PowerSums(v) => OutShape::PowerSums(v.len()), + Out::CrossPowerSums(v) => OutShape::CrossPowerSums(v.len()), } } @@ -523,8 +523,8 @@ pub(crate) fn validate( OutShape::I32(_) => Ok(()), OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) - | OutShape::CrossMoments(_) => Err(ExecError::BlendNeedsOut), + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::BlendNeedsOut), } } Terminal::ScatterOrU32 { @@ -555,8 +555,8 @@ pub(crate) fn validate( OutShape::None | OutShape::I32(_) | OutShape::I64(_) - | OutShape::Moments(_) - | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { what }), + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { what }), } } Terminal::GroupSumI32 { mask, key, val } => { @@ -573,8 +573,8 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) - | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { what: "GroupSumI32", }), } @@ -594,8 +594,8 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) - | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { what: "GroupSumViaI32", }), } @@ -622,13 +622,13 @@ pub(crate) fn validate( | OutShape::I32(_) | OutShape::I64(_) | OutShape::Mask(_) - | OutShape::Moments(_) - | OutShape::CrossMoments(_) => Err(ExecError::TerminalNeedsOut { + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { what: "GroupReduce", }), } } - Terminal::GroupMomentsI32 { mask, key, val } => { + Terminal::GroupPowerSumsI32 { mask, key, val } => { check_operand(p, planes, mask)?; written_slots.readable(mask)?; check_group_key(planes, foreign, key)?; @@ -639,13 +639,13 @@ pub(crate) fn validate( return Err(ExecError::SumRowBound { n_rows: n }); } match out { - OutShape::Moments(len) if len >= 1 => Ok(()), + OutShape::PowerSums(len) if len >= 1 => Ok(()), _ => Err(ExecError::TerminalNeedsOut { - what: "GroupMomentsI32", + what: "GroupPowerSumsI32", }), } } - Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { + Terminal::GroupCrossPowerSumsI32 { mask, key, x, y } => { check_operand(p, planes, mask)?; written_slots.readable(mask)?; check_group_key(planes, foreign, key)?; @@ -657,9 +657,9 @@ pub(crate) fn validate( return Err(ExecError::SumRowBound { n_rows: n }); } match out { - OutShape::CrossMoments(len) if len >= 1 => Ok(()), + OutShape::CrossPowerSums(len) if len >= 1 => Ok(()), _ => Err(ExecError::TerminalNeedsOut { - what: "GroupCrossMomentsI32", + what: "GroupCrossPowerSumsI32", }), } } @@ -1157,8 +1157,8 @@ pub fn reference_execute_into( } Value::GroupReduced } - Terminal::GroupMomentsI32 { mask, key, val } => { - if let Out::Moments(o) = out { + Terminal::GroupPowerSumsI32 { mask, key, val } => { + if let Out::PowerSums(o) = out { // Independent formulation: seed every slot, then walk the // survivors one row at a time, widening straight to i128 — // no ndarray kernel involved. @@ -1181,11 +1181,11 @@ pub fn reference_execute_into( slot.sum = i64::try_from(s).expect("validated 2^32-row bound keeps Σx in i64"); } } - Value::GroupMoments + Value::GroupPowerSums } - Terminal::GroupCrossMomentsI32 { mask, key, x, y } => { - if let Out::CrossMoments(o) = out { - // Independent formulation, as for GroupMomentsI32: every + Terminal::GroupCrossPowerSumsI32 { mask, key, x, y } => { + if let Out::CrossPowerSums(o) = out { + // Independent formulation, as for GroupPowerSumsI32: every // field accumulated in i128 row by row, no kernel involved. let mut acc = vec![[0i128; 6]; o.len()]; for r in survivors(mask) { @@ -1216,7 +1216,7 @@ pub fn reference_execute_into( }; } } - Value::GroupCrossMoments + Value::GroupCrossPowerSums } Terminal::Keep { mask } => { if let Out::Mask(o) = out { diff --git a/crates/lance-graph-mask-risc/src/value.rs b/crates/lance-graph-mask-risc/src/value.rs index a5c01b266..9480faf41 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -6,9 +6,9 @@ use crate::ir::Operand; /// The per-group `(n, Σx, Σy, Σx², Σy², Σxy)` accumulator -/// [`Out::CrossMoments`] carries — re-exported for the same reason. +/// [`Out::CrossPowerSums`] carries — re-exported for the same reason. pub use ndarray::simd::CrossPowerSums; -/// The per-group `(n, Σx, Σx²)` accumulator [`Out::Moments`] carries. Re-exported +/// The per-group `(n, Σx, Σx²)` accumulator [`Out::PowerSums`] carries. Re-exported /// here as shared result vocabulary so the independent oracle names the data /// type through this crate, never through the SIMD facade it falsifies. pub use ndarray::simd::PowerSums; @@ -39,14 +39,14 @@ pub enum Value { /// written, one slot per group; see [`crate::GroupFold::seed`] for what /// an empty group holds. GroupReduced, - /// [`crate::Terminal::GroupMomentsI32`]: the caller's `Out::Moments` + /// [`crate::Terminal::GroupPowerSumsI32`]: the caller's `Out::PowerSums` /// buffer was written, one [`PowerSums`] per group; an empty group /// holds [`PowerSums::default()`] (`n == 0`). - GroupMoments, - /// [`crate::Terminal::GroupCrossMomentsI32`]: the caller's - /// `Out::CrossMoments` buffer was written, one [`CrossPowerSums`] per + GroupPowerSums, + /// [`crate::Terminal::GroupCrossPowerSumsI32`]: the caller's + /// `Out::CrossPowerSums` buffer was written, one [`CrossPowerSums`] per /// group; an empty group holds [`CrossPowerSums::default()`]. - GroupCrossMoments, + GroupCrossPowerSums, /// [`crate::Terminal::MaskedStridedGroupSum`]: the widened sum, or `None` /// when it does not fit an `i64` (never a wrapped value). StridedSum(Option), @@ -67,12 +67,12 @@ pub enum Out<'a> { /// [`crate::Terminal::ScatterOrU32`]'s destination, `words_for(out_rows)` /// long. Mask(&'a mut [u64]), - /// [`crate::Terminal::GroupMomentsI32`]'s destination — one + /// [`crate::Terminal::GroupPowerSumsI32`]'s destination — one /// [`PowerSums`] per group, its length IS the group universe `K`. - Moments(&'a mut [PowerSums]), - /// [`crate::Terminal::GroupCrossMomentsI32`]'s destination — one + PowerSums(&'a mut [PowerSums]), + /// [`crate::Terminal::GroupCrossPowerSumsI32`]'s destination — one /// [`CrossPowerSums`] per group, its length IS the group universe `K`. - CrossMoments(&'a mut [CrossPowerSums]), + CrossPowerSums(&'a mut [CrossPowerSums]), } /// The lane width a predicate or terminal expects, for [`ExecError::LaneKind`]. diff --git a/crates/lance-graph-mask-risc/tests/foreign.rs b/crates/lance-graph-mask-risc/tests/foreign.rs index fee842ee0..370acd1dd 100644 --- a/crates/lance-graph-mask-risc/tests/foreign.rs +++ b/crates/lance-graph-mask-risc/tests/foreign.rs @@ -1569,7 +1569,7 @@ fn pair_key_drops_a_minor_key_at_stride() { ); } -/// FAILS IF: `Terminal::GroupMomentsI32` disagrees with the row-at-a-time +/// FAILS IF: `Terminal::GroupPowerSumsI32` disagrees with the row-at-a-time /// oracle for any key address (resident, VIA, pair) — including across TILE /// boundaries, where a sink re-seeded per tile, or a tile that re-counted a /// row, would change `n`, `Σx` or `Σx²`. @@ -1578,7 +1578,7 @@ fn pair_key_drops_a_minor_key_at_stride() { /// key hop, leave some group empty (`n == 0`) and fill some group with more /// than one row. #[test] -fn group_moments_match_the_oracle_for_every_key_address() { +fn group_power_sums_match_the_oracle_for_every_key_address() { let mut multi_tile = false; let mut empty_group = false; let mut multi_row_group = false; @@ -1623,7 +1623,7 @@ fn group_moments_match_the_oracle_for_every_key_address() { under: None, dst: 0, }], - Terminal::GroupMomentsI32 { + Terminal::GroupPowerSumsI32 { mask: S0, key, val: 2, @@ -1644,13 +1644,13 @@ fn group_moments_match_the_oracle_for_every_key_address() { &planes, &foreign, &mut scratch, - Out::Moments(&mut got_out), + Out::PowerSums(&mut got_out), ) .expect("runs"); let mut want_out = vec![dirty; groups]; - let want = reference_execute_into(&p, &planes, &foreign, Out::Moments(&mut want_out)) + let want = reference_execute_into(&p, &planes, &foreign, Out::PowerSums(&mut want_out)) .expect("oracle runs"); - assert_eq!(got, Value::GroupMoments, "n={n} {key:?}"); + assert_eq!(got, Value::GroupPowerSums, "n={n} {key:?}"); assert_eq!(got, want, "n={n} {key:?}: value"); assert_eq!(got_out, want_out, "n={n} {key:?}: sink"); empty_group |= got_out.iter().any(|g| g.n == 0); @@ -1662,11 +1662,11 @@ fn group_moments_match_the_oracle_for_every_key_address() { assert!(multi_row_group, "must fold several rows into one group"); } -/// FAILS IF: `GroupMomentsI32` accepts a wrong-width value or key lane, a +/// FAILS IF: `GroupPowerSumsI32` accepts a wrong-width value or key lane, a /// missing or wrong-shaped sink, or runs under an extent (a partial /// population would silently report partial moments as whole). #[test] -fn group_moments_refuses_malformed_programs() { +fn group_power_sums_refuses_malformed_programs() { let n = 130; let fx = Fixture::new(n, 10, 4, 0x5EF); let (lanes, masks) = fx.planes(); @@ -1694,33 +1694,33 @@ fn group_moments_refuses_malformed_programs() { // `val` must be an I32 lane (lane 1 is U32). assert!(matches!( run( - Terminal::GroupMomentsI32 { + Terminal::GroupPowerSumsI32 { mask: S0, key: GroupKey::Lane(3), val: 1 }, - Out::Moments(&mut sink) + Out::PowerSums(&mut sink) ), Err(ExecError::LaneKind { .. }) )); // The key must be a U32 lane (lane 2 is I32). assert!(matches!( run( - Terminal::GroupMomentsI32 { + Terminal::GroupPowerSumsI32 { mask: S0, key: GroupKey::Lane(2), val: 2 }, - Out::Moments(&mut sink) + Out::PowerSums(&mut sink) ), Err(ExecError::LaneKind { .. }) )); // A moments sink is required; an i64 sink is the wrong shape. let mut i64_sink = vec![0i64; 4]; - for out in [Out::None, Out::I64(&mut i64_sink), Out::Moments(&mut [])] { + for out in [Out::None, Out::I64(&mut i64_sink), Out::PowerSums(&mut [])] { assert_eq!( run( - Terminal::GroupMomentsI32 { + Terminal::GroupPowerSumsI32 { mask: S0, key: GroupKey::Lane(3), val: 2 @@ -1728,7 +1728,7 @@ fn group_moments_refuses_malformed_programs() { out ), Err(ExecError::TerminalNeedsOut { - what: "GroupMomentsI32" + what: "GroupPowerSumsI32" }) ); } @@ -1736,7 +1736,7 @@ fn group_moments_refuses_malformed_programs() { // part of the population must never be reported as moments of the whole. let p = Program::new( vec![], - Terminal::GroupMomentsI32 { + Terminal::GroupPowerSumsI32 { mask: Operand::Plane(0), key: GroupKey::Lane(3), val: 2, @@ -1752,9 +1752,16 @@ fn group_moments_refuses_malformed_programs() { let mut s = Scratch::for_program(&p, n).expect("scratch"); let mut sink = vec![PowerSums::default(); 4]; assert_eq!( - execute_extent(&p, &planes, &none, &mut s, Out::Moments(&mut sink), 10..20), + execute_extent( + &p, + &planes, + &none, + &mut s, + Out::PowerSums(&mut sink), + 10..20 + ), Err(ExecError::ExtentUnsupported { - what: "GroupMomentsI32" + what: "GroupPowerSumsI32" }) ); assert_eq!( @@ -1777,12 +1784,12 @@ fn y_lane(n: usize, seed: u64) -> Vec { .collect() } -/// FAILS IF: `Terminal::GroupCrossMomentsI32` disagrees with the +/// FAILS IF: `Terminal::GroupCrossPowerSumsI32` disagrees with the /// row-at-a-time oracle for any key address, across tile boundaries, or /// pairs `x[i]` with a `y` from another row (the oracle reads both at `r`). /// Anti-vacuity: multi-tile, an empty group, a multi-row group. #[test] -fn group_cross_moments_match_the_oracle_for_every_key_address() { +fn group_cross_power_sums_match_the_oracle_for_every_key_address() { let mut multi_tile = false; let mut empty_group = false; let mut multi_row_group = false; @@ -1828,7 +1835,7 @@ fn group_cross_moments_match_the_oracle_for_every_key_address() { under: None, dst: 0, }], - Terminal::GroupCrossMomentsI32 { + Terminal::GroupCrossPowerSumsI32 { mask: S0, key, x: 2, @@ -1849,14 +1856,14 @@ fn group_cross_moments_match_the_oracle_for_every_key_address() { &planes, &foreign, &mut scratch, - Out::CrossMoments(&mut got_out), + Out::CrossPowerSums(&mut got_out), ) .expect("runs"); let mut want_out = vec![dirty; groups]; let want = - reference_execute_into(&p, &planes, &foreign, Out::CrossMoments(&mut want_out)) + reference_execute_into(&p, &planes, &foreign, Out::CrossPowerSums(&mut want_out)) .expect("oracle runs"); - assert_eq!(got, Value::GroupCrossMoments, "n={n} {key:?}"); + assert_eq!(got, Value::GroupCrossPowerSums, "n={n} {key:?}"); assert_eq!(got, want, "n={n} {key:?}: value"); assert_eq!(got_out, want_out, "n={n} {key:?}: sink"); empty_group |= got_out.iter().any(|g| g.n == 0); @@ -1867,10 +1874,10 @@ fn group_cross_moments_match_the_oracle_for_every_key_address() { } /// FAILS IF: the cross terminal accepts a wrong-width `x` or `y` lane, a -/// missing or wrong-shaped sink (an `Out::Moments` is the WRONG shape), or +/// missing or wrong-shaped sink (an `Out::PowerSums` is the WRONG shape), or /// runs under a partial extent. #[test] -fn group_cross_moments_refuses_malformed_programs() { +fn group_cross_power_sums_refuses_malformed_programs() { let n = 130; let fx = Fixture::new(n, 10, 4, 0x5F0); let ys = y_lane(n, 1); @@ -1890,7 +1897,7 @@ fn group_cross_moments_refuses_malformed_programs() { planes: &[], lanes: &[], }; - let term = |x: u16, y: u16| Terminal::GroupCrossMomentsI32 { + let term = |x: u16, y: u16| Terminal::GroupCrossPowerSumsI32 { mask: Operand::Plane(0), key: GroupKey::Lane(3), x, @@ -1903,7 +1910,7 @@ fn group_cross_moments_refuses_malformed_programs() { for (x, y) in [(1u16, 4u16), (2, 1)] { assert!( matches!( - run(term(x, y), Out::CrossMoments(&mut sink)), + run(term(x, y), Out::CrossPowerSums(&mut sink)), Err(ExecError::LaneKind { .. }) ), "x={x} y={y}" @@ -1912,13 +1919,13 @@ fn group_cross_moments_refuses_malformed_programs() { let mut univariate = vec![PowerSums::default(); 4]; for out in [ Out::None, - Out::Moments(&mut univariate), - Out::CrossMoments(&mut []), + Out::PowerSums(&mut univariate), + Out::CrossPowerSums(&mut []), ] { assert_eq!( run(term(2, 4), out), Err(ExecError::TerminalNeedsOut { - what: "GroupCrossMomentsI32" + what: "GroupCrossPowerSumsI32" }) ); } @@ -1930,11 +1937,11 @@ fn group_cross_moments_refuses_malformed_programs() { &planes, &none, &mut s, - Out::CrossMoments(&mut sink), + Out::CrossPowerSums(&mut sink), 10..20 ), Err(ExecError::ExtentUnsupported { - what: "GroupCrossMomentsI32" + what: "GroupCrossPowerSumsI32" }) ); assert_eq!( @@ -1944,7 +1951,7 @@ fn group_cross_moments_refuses_malformed_programs() { ); // x == y is legal and folds the univariate moments. let mut same = vec![CrossPowerSums::default(); 4]; - run(term(2, 2), Out::CrossMoments(&mut same)).expect("x == y is legal"); + run(term(2, 2), Out::CrossPowerSums(&mut same)).expect("x == y is legal"); assert!(same .iter() .all(|g| g.sum_x == g.sum_y && u128::try_from(g.sum_xy) == Ok(g.sum_x_sq))); diff --git a/crates/lance-graph-quack/tests/duckdb_differential.rs b/crates/lance-graph-quack/tests/duckdb_differential.rs index 8113fb011..d15aaa3bc 100644 --- a/crates/lance-graph-quack/tests/duckdb_differential.rs +++ b/crates/lance-graph-quack/tests/duckdb_differential.rs @@ -264,9 +264,9 @@ fn run_query(id: &str, planes: &Planes<'_>, filter: Filter, agg: Agg) -> (String Value::StridedSum(_) => { panic!("case {id}: no MaskedStridedGroupSum case in this suite") } - Value::GroupMoments => panic!("case {id}: no GroupMomentsI32 case in this suite"), - Value::GroupCrossMoments => { - panic!("case {id}: no GroupCrossMomentsI32 case in this suite") + Value::GroupPowerSums => panic!("case {id}: no GroupPowerSumsI32 case in this suite"), + Value::GroupCrossPowerSums => { + panic!("case {id}: no GroupCrossPowerSumsI32 case in this suite") } }; let out_bytes = encoded.len(); From ba2e48fb3d0020ab10a53682eeadd3fb9a2d24fc Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 19:02:52 +0000 Subject: [PATCH 10/14] mask-risc: admit grouped power-sums terminals under a partial extent GroupPowerSumsI32 and GroupCrossPowerSumsI32 were refused by execute_extent although their per-extent results merge exactly by PowerSums::checked_merge / CrossPowerSums::checked_merge. Admit them: each call seeds its own sink, and the edge word's mask is clipped to the extent like the other mergeable folds. tests/extent.rs pins the law: whole-population execution equals the group-by-group checked_merge of any partition, in any order, for the resident, VIA and pair key forms, at two tile widths. The two refusal assertions in tests/foreign.rs are removed with this change. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- ...9-grouped-moments-fold-anova-over-masks.md | 2 +- crates/lance-graph-mask-risc/src/exec.rs | 13 +- crates/lance-graph-mask-risc/tests/extent.rs | 225 +++++++++++++++++- crates/lance-graph-mask-risc/tests/foreign.rs | 69 +----- 4 files changed, 237 insertions(+), 72 deletions(-) diff --git a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md index c32426630..79149a567 100644 --- a/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md +++ b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md @@ -4,7 +4,7 @@ ## What landed - ndarray: `GroupPowerSums { n: u64, sum: i64, sum_sq: i128 }` and `masked_group_moments_i32{,_via,_pair}` — `(n, Σx, Σx²)` per group in one pass over the rows a mask selects, over the existing `group_walk` (now generic over its slot type). `checked_merge` is exact integer addition; exact up to `GROUP_MOMENTS_MAX_ROWS = 2^32` rows per group (the `i64` sum is the binding field). -- mask-risc: `Terminal::GroupPowerSumsI32 { mask, key: GroupKey, val }` → `Out::PowerSums`. Own terminal, not a `GroupFold` member (a `GroupFold` slot is one seeded `i64`). Refuses planes past `MASKED_SUM_I32_MAX_ROWS` and partial extents. +- mask-risc: `Terminal::GroupPowerSumsI32 { mask, key: GroupKey, val }` → `Out::PowerSums`. Own terminal, not a `GroupFold` member (a `GroupFold` slot is one seeded `i64`). Refuses planes past `MASKED_SUM_I32_MAX_ROWS`. Partial extents: admitted since #1323's extent wiring — per-extent sinks merge by `checked_merge` (pinned in `tests/extent.rs`). - jc: `anova_from_power_sums` / `eta_squared_from_power_sums`. F/p/η² policy is shared with `anova_one_way` / `eta_squared` via `anova_from_ss` / `eta_from_ss`. Sums of squares come from exact `i128` quantities; no two large floats are subtracted. ## Measured diff --git a/crates/lance-graph-mask-risc/src/exec.rs b/crates/lance-graph-mask-risc/src/exec.rs index 51d2311c6..60b783a79 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -1342,7 +1342,10 @@ pub fn execute_into( /// `execute_into` is this call with `0..n_rows`, which accepts every /// terminal. A partial extent accepts the terminals whose per-extent results /// merge by a shipped law — `Count` (sum), `Any` (or), `All` (and), -/// `MaskedSumI32` (sum), `MaskedMinI32` / `MaskedMaxI32` (min / max) — plus +/// `MaskedSumI32` (sum), `MaskedMinI32` / `MaskedMaxI32` (min / max), +/// `GroupPowerSumsI32` / `GroupCrossPowerSumsI32` (each extent gets its own +/// sink, seeded fresh; partial sinks combine group-by-group with +/// `PowerSums::checked_merge` / `CrossPowerSums::checked_merge`) — plus /// `Keep`, which writes only the in-extent bits of its population-addressed /// [`Out::Mask`] and leaves every other bit as the caller holds it (so /// disjoint extents compose into one buffer in any SEQUENTIAL order). Anything else is @@ -1414,6 +1417,8 @@ fn precheck( | Terminal::MaskedMinI32 { .. } | Terminal::MaskedMaxI32 { .. } | Terminal::MaskedStridedGroupSum { .. } + | Terminal::GroupPowerSumsI32 { .. } + | Terminal::GroupCrossPowerSumsI32 { .. } | Terminal::Keep { .. } => None, Terminal::BlendI32 { .. } => Some("BlendI32"), Terminal::ScatterOrU32 { .. } => Some("ScatterOrU32"), @@ -1422,8 +1427,6 @@ fn precheck( Terminal::GroupSumI32 { .. } => Some("GroupSumI32"), Terminal::GroupSumViaI32 { .. } => Some("GroupSumViaI32"), Terminal::GroupReduce { .. } => Some("GroupReduce"), - Terminal::GroupPowerSumsI32 { .. } => Some("GroupPowerSumsI32"), - Terminal::GroupCrossPowerSumsI32 { .. } => Some("GroupCrossPowerSumsI32"), }; if let Some(what) = refused { return Err(ExecError::ExtentUnsupported { what }); @@ -1998,7 +2001,7 @@ pub fn execute_compiled( // bound; one delegation per tile (law L3) into the sink // seeded above with `PowerSums::default()`. if let Out::PowerSums(o) = &mut out { - let m = read(planes, &slots, mask, t); + let m = clip(read(planes, &slots, mask, t), edge, &mut eb, false); let v = lane_i32(planes, val, t); match key { GroupKey::Lane(k) => { @@ -2026,7 +2029,7 @@ pub fn execute_compiled( // Same contract as GroupPowerSumsI32, both lanes read in place: // one delegation per tile into the sink seeded above. if let Out::CrossPowerSums(o) = &mut out { - let m = read(planes, &slots, mask, t); + let m = clip(read(planes, &slots, mask, t), edge, &mut eb, false); let (xs, ys) = (lane_i32(planes, x, t), lane_i32(planes, y, t)); match key { GroupKey::Lane(k) => { diff --git a/crates/lance-graph-mask-risc/tests/extent.rs b/crates/lance-graph-mask-risc/tests/extent.rs index cd4cf049a..c112f6a36 100644 --- a/crates/lance-graph-mask-risc/tests/extent.rs +++ b/crates/lance-graph-mask-risc/tests/extent.rs @@ -9,7 +9,9 @@ //! The gates: //! - the edge matrix: every result equals a scalar oracle that knows nothing //! about tiles, words or edges; -//! - split composition: whole == the merge of any partition, in any order; +//! - split composition: whole == the merge of any partition, in any order — +//! including the grouped power-sums sinks, merged group-by-group with +//! `checked_merge`; //! - the non-rebasing falsifier: an extent that does not start at row 0 must //! read the absolute rows, and the rebased reading is shown to differ; //! - the structural gate (tiles visited scale with the extent, not with @@ -19,8 +21,8 @@ use lance_graph_mask_risc::exec::{execute_extent, execute_into, Scratch}; use lance_graph_mask_risc::{ - ExecError, Foreign, ForeignPlane, LaneRef, MaskOp, Operand, Out, Planes, Pred, Program, - Terminal, Value, + scratch_words_for, CrossPowerSums, ExecError, Foreign, ForeignPlane, GroupKey, LaneRef, MaskOp, + Operand, Out, Planes, PowerSums, Pred, Program, Terminal, Value, }; fn lcg(seed: &mut u64) -> u64 { @@ -584,3 +586,220 @@ fn the_extent_never_slices_a_foreign_plane() { let v = execute_extent(&p, &planes, &foreign, &mut s, Out::None, lo..hi).expect("gather"); assert_eq!(v, Value::Count(want)); } + +/// The grouped power-sums terminals over one extent, into a FRESH sink the +/// caller owns. The sink starts dirty on purpose: the terminal must seed it. +fn power_sums_part( + p: &Program, + planes: &Planes<'_>, + foreign: &Foreign<'_>, + tile_words: usize, + groups: usize, + ext: std::ops::Range, +) -> Vec { + let slots = p.scratch_slots as usize; + let mut buf = vec![0u64; scratch_words_for(tile_words, slots).expect("sized")]; + let mut s = Scratch::over(&mut buf, tile_words, slots).expect("carves"); + let dirty = PowerSums { + n: 7, + sum: -7, + sum_sq: 7, + }; + let mut sink = vec![dirty; groups]; + let v = execute_extent(p, planes, foreign, &mut s, Out::PowerSums(&mut sink), ext) + .expect("a partial extent is admitted"); + assert_eq!(v, Value::GroupPowerSums); + sink +} + +fn cross_part( + p: &Program, + planes: &Planes<'_>, + foreign: &Foreign<'_>, + tile_words: usize, + groups: usize, + ext: std::ops::Range, +) -> Vec { + let slots = p.scratch_slots as usize; + let mut buf = vec![0u64; scratch_words_for(tile_words, slots).expect("sized")]; + let mut s = Scratch::over(&mut buf, tile_words, slots).expect("carves"); + let dirty = CrossPowerSums { + n: 7, + sum_x: -7, + sum_y: 7, + sum_x_sq: 7, + sum_y_sq: 7, + sum_xy: -7, + }; + let mut sink = vec![dirty; groups]; + let v = execute_extent( + p, + planes, + foreign, + &mut s, + Out::CrossPowerSums(&mut sink), + ext, + ) + .expect("a partial extent is admitted"); + assert_eq!(v, Value::GroupCrossPowerSums); + sink +} + +/// FAILS IF: a grouped power-sums terminal over a partial extent counts a +/// row outside it (an unclipped edge word), re-counts a row two extents both +/// claim, or its sink is not re-seeded per call — i.e. whole-population +/// execution must equal the group-by-group `checked_merge` of ANY partition, +/// in any order, for every key address (resident, VIA, pair). +/// +/// Anti-vacuity: some partition must put rows of ONE group in two different +/// extents (else the merge never adds), and some cut must split a 64-row +/// word with selected rows on both sides of it (else edge clipping is never +/// exercised). +#[test] +fn grouped_power_sums_partials_merge_to_the_whole() { + let mut merged_nontrivially = false; + let mut split_a_live_word = false; + for n in [1317usize, 4096 + 37] { + let mut seed = 0x9_0E5 ^ n as u64; + let sel: Vec = (0..n).map(|_| !lcg(&mut seed).is_multiple_of(3)).collect(); + let pl = plane(n, |r| sel[r]); + let key: Vec = (0..n).map(|_| (lcg(&mut seed) % 6) as u32).collect(); // 5 = drop + let x: Vec = (0..n) + .map(|i| match i % 17 { + 3 => i32::MIN, + 9 => i32::MAX, + _ => (lcg(&mut seed) % 200_001) as i32 - 100_000, + }) + .collect(); + let y: Vec = (0..n) + .map(|i| match i % 13 { + 5 => i32::MAX, + _ => (lcg(&mut seed) % 2001) as i32 - 1000, + }) + .collect(); + let fk: Vec = (0..n).map(|_| (lcg(&mut seed) % 8) as u32).collect(); // 7 = past table + let hi: Vec = (0..n).map(|_| (lcg(&mut seed) % 3) as u32).collect(); + let lo: Vec = (0..n).map(|_| (lcg(&mut seed) % 5) as u32).collect(); // 4 = >= stride + let table: Vec = vec![0, 4, 2, 9, 1, 3, 4]; // 9 = second-hop drop + let masks: [&[u64]; 1] = [&pl]; + let lanes = [ + LaneRef::U32(&key), + LaneRef::I32(&x), + LaneRef::I32(&y), + LaneRef::U32(&fk), + LaneRef::U32(&hi), + LaneRef::U32(&lo), + ]; + let planes = Planes { + n_rows: n, + masks: &masks, + lanes: &lanes, + }; + let flanes = [LaneRef::U32(&table)]; + let foreign = Foreign { + planes: &[], + lanes: &flanes, + }; + let mut partitions: Vec> = [0, 1, 63, 64, 65, 127, 129, n / 2, n - 1, n] + .into_iter() + .map(|k| vec![0, k, n]) + .collect(); + for _ in 0..20 { + let a = (lcg(&mut seed) as usize) % (n + 1); + let b = (lcg(&mut seed) as usize) % (n + 1); + partitions.push(vec![0, a.min(b), a.max(b), n]); + } + for cuts in &partitions { + split_a_live_word |= cuts[1..cuts.len() - 1].iter().any(|&c| { + c % 64 != 0 + && (c - c % 64..c).any(|r| bit(&pl, r)) + && (c..(c - c % 64 + 64).min(n)).any(|r| bit(&pl, r)) + }); + } + let keys = [ + (GroupKey::Lane(0), 5usize), + (GroupKey::Via { fk: 3, key: 0 }, 5), + ( + GroupKey::Pair { + hi: 4, + lo: 5, + stride: 4, + }, + 12, + ), + ]; + let words = words(n); + for (gk, groups) in keys { + for tile_words in [2usize, words] { + let ps = Program::new( + vec![], + Terminal::GroupPowerSumsI32 { + mask: Operand::Plane(0), + key: gk, + val: 1, + }, + ); + let cs = Program::new( + vec![], + Terminal::GroupCrossPowerSumsI32 { + mask: Operand::Plane(0), + key: gk, + x: 1, + y: 2, + }, + ); + let whole_ps = power_sums_part(&ps, &planes, &foreign, tile_words, groups, 0..n); + let whole_cs = cross_part(&cs, &planes, &foreign, tile_words, groups, 0..n); + for cuts in &partitions { + let spans: Vec<_> = cuts.windows(2).map(|w| w[0]..w[1]).collect(); + let ps_parts: Vec> = spans + .iter() + .map(|e| { + power_sums_part(&ps, &planes, &foreign, tile_words, groups, e.clone()) + }) + .collect(); + let cs_parts: Vec> = spans + .iter() + .map(|e| cross_part(&cs, &planes, &foreign, tile_words, groups, e.clone())) + .collect(); + merged_nontrivially |= + (0..groups).any(|g| ps_parts.iter().filter(|p| p[g].n > 0).count() > 1); + let k = spans.len(); + for order in [ + (0..k).collect::>(), + (0..k).rev().collect(), + (0..k).map(|i| (i + 1) % k).collect(), + ] { + let ps_merged: Vec = (0..groups) + .map(|g| { + order + .iter() + .map(|&i| ps_parts[i][g]) + .try_fold(PowerSums::default(), PowerSums::checked_merge) + .expect("in-width merge") + }) + .collect(); + let cs_merged: Vec = (0..groups) + .map(|g| { + order + .iter() + .map(|&i| cs_parts[i][g]) + .try_fold( + CrossPowerSums::default(), + CrossPowerSums::checked_merge, + ) + .expect("in-width merge") + }) + .collect(); + let at = + format!("n={n} {gk:?} tile={tile_words} cuts={cuts:?} order={order:?}"); + assert_eq!(ps_merged, whole_ps, "power sums {at}"); + assert_eq!(cs_merged, whole_cs, "cross power sums {at}"); + } + } + } + } + } + assert!(merged_nontrivially, "some group must span two extents"); + assert!(split_a_live_word, "some cut must split a live 64-row word"); +} diff --git a/crates/lance-graph-mask-risc/tests/foreign.rs b/crates/lance-graph-mask-risc/tests/foreign.rs index 370acd1dd..e131c98ca 100644 --- a/crates/lance-graph-mask-risc/tests/foreign.rs +++ b/crates/lance-graph-mask-risc/tests/foreign.rs @@ -2,7 +2,7 @@ //! against the row-at-a-time oracle — the same differential shape //! `tests/differential.rs` uses, extended over a SECOND, foreign row space. -use lance_graph_mask_risc::exec::{execute_extent, execute_into, Scratch}; +use lance_graph_mask_risc::exec::{execute_into, Scratch}; use lance_graph_mask_risc::reference::{reference_execute_into, reference_scratch_with_foreign}; use lance_graph_mask_risc::{ scratch_words_for, words_for, CrossPowerSums, ExecError, Foreign, ForeignPlane, GroupFold, @@ -1662,9 +1662,9 @@ fn group_power_sums_match_the_oracle_for_every_key_address() { assert!(multi_row_group, "must fold several rows into one group"); } -/// FAILS IF: `GroupPowerSumsI32` accepts a wrong-width value or key lane, a -/// missing or wrong-shaped sink, or runs under an extent (a partial -/// population would silently report partial moments as whole). +/// FAILS IF: `GroupPowerSumsI32` accepts a wrong-width value or key lane, or +/// a missing or wrong-shaped sink. (Partial extents are admitted; their merge +/// law is pinned in `tests/extent.rs`.) #[test] fn group_power_sums_refuses_malformed_programs() { let n = 130; @@ -1732,43 +1732,6 @@ fn group_power_sums_refuses_malformed_programs() { }) ); } - // A partial extent is refused before anything is written: moments over - // part of the population must never be reported as moments of the whole. - let p = Program::new( - vec![], - Terminal::GroupPowerSumsI32 { - mask: Operand::Plane(0), - key: GroupKey::Lane(3), - val: 2, - }, - ); - let pl = vec![u64::MAX; words_for(n)]; - let masks: [&[u64]; 1] = [&pl]; - let planes = Planes { - n_rows: n, - masks: &masks, - lanes: &lanes, - }; - let mut s = Scratch::for_program(&p, n).expect("scratch"); - let mut sink = vec![PowerSums::default(); 4]; - assert_eq!( - execute_extent( - &p, - &planes, - &none, - &mut s, - Out::PowerSums(&mut sink), - 10..20 - ), - Err(ExecError::ExtentUnsupported { - what: "GroupPowerSumsI32" - }) - ); - assert_eq!( - sink, - vec![PowerSums::default(); 4], - "a refusal writes nothing" - ); } /// A second `I32` lane for the cross terminal (lane 4), independent of @@ -1874,8 +1837,8 @@ fn group_cross_power_sums_match_the_oracle_for_every_key_address() { } /// FAILS IF: the cross terminal accepts a wrong-width `x` or `y` lane, a -/// missing or wrong-shaped sink (an `Out::PowerSums` is the WRONG shape), or -/// runs under a partial extent. +/// missing or wrong-shaped sink (an `Out::PowerSums` is the WRONG shape). +/// (Partial extents are admitted; see `tests/extent.rs`.) #[test] fn group_cross_power_sums_refuses_malformed_programs() { let n = 130; @@ -1929,26 +1892,6 @@ fn group_cross_power_sums_refuses_malformed_programs() { }) ); } - let p = Program::new(vec![], term(2, 4)); - let mut s = Scratch::for_program(&p, n).expect("scratch"); - assert_eq!( - execute_extent( - &p, - &planes, - &none, - &mut s, - Out::CrossPowerSums(&mut sink), - 10..20 - ), - Err(ExecError::ExtentUnsupported { - what: "GroupCrossPowerSumsI32" - }) - ); - assert_eq!( - sink, - vec![CrossPowerSums::default(); 4], - "a refusal writes nothing" - ); // x == y is legal and folds the univariate moments. let mut same = vec![CrossPowerSums::default(); 4]; run(term(2, 2), Out::CrossPowerSums(&mut same)).expect("x == y is legal"); From 898179fc24a3ab9e526a66f858738ecd7d79f9a7 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 19:25:26 +0000 Subject: [PATCH 11/14] perturbation-sim: rescale before narrowing the angle-covariance sandwich to f32 An f64 entry of sigma_p or L+ past f32::MAX became infinity on the cast (and one below f32's smallest normal flushed to zero), so a finite f64 product could come back infinite or NaN. Each input is now divided by its largest absolute entry before narrowing and the result multiplied back in f64; the sandwich is bilinear, so this is exact up to rounding. Non-finite sigma_p entries are refused. Tests: magnitudes 1e45 and 1e-45 match the f64 triple product within 1e-5; a non-finite entry panics. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/perturbation-sim/src/angle_cov.rs | 79 ++++++++++++++++++++---- 1 file changed, 67 insertions(+), 12 deletions(-) diff --git a/crates/perturbation-sim/src/angle_cov.rs b/crates/perturbation-sim/src/angle_cov.rs index dc2dc25aa..b59cf0862 100644 --- a/crates/perturbation-sim/src/angle_cov.rs +++ b/crates/perturbation-sim/src/angle_cov.rs @@ -18,9 +18,13 @@ //! //! # Precision and size //! -//! `CovHighD` is `f32` and dense `O(N³)`; this crate is `f64`. Values are -//! narrowed to `f32` once on the way in and widened on the way out, so expect -//! roughly `1e-6` relative error against an `f64` triple product. `N` is a +//! `CovHighD` is `f32` and dense `O(N³)`; this crate is `f64`. Each input is +//! divided by its largest absolute entry before it is narrowed to `f32`, and +//! the result is multiplied back in `f64` (the sandwich is bilinear, so this +//! is exact up to rounding). Magnitude therefore never decides whether an +//! entry survives the narrowing: a `Σp` past `f32::MAX` or below its smallest +//! normal is as accurate as one near `1`. Expect roughly `1e-6` error relative +//! to the largest output entry. `N` is a //! compile-time constant in `CovHighD`, while a [`crate::graph::Grid`]'s bus //! count is a runtime value: the caller picks `N` and the call panics when the //! decomposition disagrees with it. A runtime-sized sandwich would be a change @@ -45,8 +49,9 @@ pub const SYMMETRY_TOL: f64 = 1e-9; /// uses. `sigma_p` and the result are row-major `N×N`. /// /// # Panics -/// If `eig.n != N`, if `sigma_p.len() != N*N`, or if `sigma_p` is not -/// symmetric within [`SYMMETRY_TOL`] (relative to its largest entry). +/// If `eig.n != N`, if `sigma_p.len() != N*N`, if any entry of `sigma_p` is +/// not finite, or if `sigma_p` is not symmetric within [`SYMMETRY_TOL`] +/// (relative to its largest entry). pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f64) -> Vec { assert_eq!( eig.n, N, @@ -54,6 +59,10 @@ pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f eig.n ); assert_eq!(sigma_p.len(), N * N, "sigma_p must be N*N"); + assert!( + sigma_p.iter().all(|v| v.is_finite()), + "sigma_p has a non-finite entry" + ); let scale = sigma_p .iter() .fold(0.0_f64, |m, v| m.max(v.abs())) @@ -69,14 +78,22 @@ pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f } let l_plus = eig.pseudo_inverse(rel_tol); - let m = CovHighD::::from_symmetric_fn(|i, j| l_plus[i * N + j] as f32); - let s = CovHighD::::from_symmetric_fn(|i, j| sigma_p[i * N + j] as f32); + let l_scale = l_plus + .iter() + .fold(0.0_f64, |m, v| m.max(v.abs())) + .max(f64::MIN_POSITIVE); + // Both inputs land in [-1, 1] before narrowing, so no entry overflows to + // infinity or flushes to zero merely because of its magnitude; the f32 + // sandwich then sums at most N² products bounded by 1. + let m = CovHighD::::from_symmetric_fn(|i, j| (l_plus[i * N + j] / l_scale) as f32); + let s = CovHighD::::from_symmetric_fn(|i, j| (sigma_p[i * N + j] / scale) as f32); let out = s.sandwich(&m); + let back = scale * l_scale * l_scale; let mut dense = vec![0.0_f64; N * N]; for i in 0..N { for j in 0..N { - dense[i * N + j] = out.get(i, j) as f64; + dense[i * N + j] = out.get(i, j) as f64 * back; } } dense @@ -150,10 +167,8 @@ mod tests { ); } - #[test] - fn matches_f64_dense_triple_product() { - let eig = symmetric_eigen(&grid().laplacian(), N); - let sp = sigma_p(); + /// `L⁺ Σp L⁺` in plain `f64`, the reference every test compares against. + fn triple_product(eig: &Eigen, sp: &[f64]) -> Vec { let l = eig.pseudo_inverse(TOL); let mut want = vec![0.0; N * N]; for i in 0..N { @@ -167,6 +182,14 @@ mod tests { want[i * N + m] = acc; } } + want + } + + #[test] + fn matches_f64_dense_triple_product() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let sp = sigma_p(); + let want = triple_product(&eig, &sp); let got = angle_covariance::(&eig, &sp, TOL); let e = rel_err(&got, &want); eprintln!("sandwich vs f64 triple product rel err {e:e}"); @@ -223,6 +246,38 @@ mod tests { assert!(e < 0.03, "Monte-Carlo vs sandwich rel err {e:.4}"); } + /// FAILS IF: an input is narrowed to `f32` at its raw magnitude, so a + /// `Σp` past `f32::MAX` becomes infinity (or one below the smallest + /// normal flushes to zero) although the `f64` product is finite. + #[test] + fn magnitude_outside_f32_range_is_exact_up_to_rounding() { + let eig = symmetric_eigen(&grid().laplacian(), N); + for k in [1e45_f64, 1e-45] { + let sp: Vec = sigma_p().iter().map(|v| v * k).collect(); + assert!( + sp.iter().any(|v| (*v as f32).is_infinite()) + || sp + .iter() + .all(|v| *v == 0.0 || (*v as f32).abs() < f32::MIN_POSITIVE), + "fixture must leave f32's normal range at k = {k:e}" + ); + let want = triple_product(&eig, &sp); + let got = angle_covariance::(&eig, &sp, TOL); + assert!(got.iter().all(|v| v.is_finite()), "k = {k:e}: non-finite"); + let e = rel_err(&got, &want); + assert!(e < 1e-5, "k = {k:e}: rel err {e:e}"); + } + } + + #[test] + #[should_panic(expected = "non-finite")] + fn refuses_a_non_finite_injection_covariance() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let mut sp = sigma_p(); + sp[0] = f64::INFINITY; + let _ = angle_covariance::(&eig, &sp, TOL); + } + #[test] #[should_panic(expected = "CovHighD has N")] fn refuses_a_const_n_that_disagrees_with_the_grid() { From 401265aff9ec0be85df7788316771be05a1f66fb Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 19:26:13 +0000 Subject: [PATCH 12/14] mask-risc: list the grouped power-sums terminals in the ExtentUnsupported doc Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/lance-graph-mask-risc/src/value.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/crates/lance-graph-mask-risc/src/value.rs b/crates/lance-graph-mask-risc/src/value.rs index 9480faf41..30dff5bf6 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -216,7 +216,8 @@ pub enum ExecError { /// write). The whole-population extent `[0, n_rows)` accepts every /// terminal; a partial one accepts `Count`, `Any`, `All`, /// `MaskedSumI32`, `MaskedMinI32`, `MaskedMaxI32`, `MaskedStridedGroupSum` - /// (a sum merges by addition) and `Keep`. `what` - /// names the refused terminal. + /// (a sum merges by addition), `GroupPowerSumsI32` / + /// `GroupCrossPowerSumsI32` (fresh per-extent sinks, merged group-by-group + /// with `checked_merge`) and `Keep`. `what` names the refused terminal. ExtentUnsupported { what: &'static str }, } From 74e621b3644e1891cb078b3d152e5e35698c1582 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 19:32:20 +0000 Subject: [PATCH 13/14] perturbation-sim: refuse a sigma_p wider than f32, rescale per entry Two review findings on the f32 sandwich: - A huge sigma_p entry (e.g. in L+'s null space) set the global scale, so the entries that determine the output flushed to zero in f32 and the result was silently wrong. A nonzero entry that would fall below f32::MIN_POSITIVE relative to the max is now refused. The module doc now states the error bound relative to the input magnitude, not the output. - The combined factor scale * l_scale^2 could overflow to infinity and turn an exact-zero entry into NaN. Each entry is now rescaled factor by factor. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/perturbation-sim/src/angle_cov.rs | 77 +++++++++++++++++++++--- 1 file changed, 69 insertions(+), 8 deletions(-) diff --git a/crates/perturbation-sim/src/angle_cov.rs b/crates/perturbation-sim/src/angle_cov.rs index b59cf0862..e03834c6b 100644 --- a/crates/perturbation-sim/src/angle_cov.rs +++ b/crates/perturbation-sim/src/angle_cov.rs @@ -21,10 +21,17 @@ //! `CovHighD` is `f32` and dense `O(N³)`; this crate is `f64`. Each input is //! divided by its largest absolute entry before it is narrowed to `f32`, and //! the result is multiplied back in `f64` (the sandwich is bilinear, so this -//! is exact up to rounding). Magnitude therefore never decides whether an -//! entry survives the narrowing: a `Σp` past `f32::MAX` or below its smallest -//! normal is as accurate as one near `1`. Expect roughly `1e-6` error relative -//! to the largest output entry. `N` is a +//! is exact up to rounding). A uniformly huge or tiny `Σp` is therefore as +//! accurate as one near `1`. +//! +//! The error bound is relative to the INPUT magnitude `max|Σp| · max|L⁺|²`, +//! not to the output: expect roughly `1e-6` of that. When most of `Σp`'s mass +//! lies in `L⁺`'s null space (e.g. a large variance on an isolated bus), the +//! output is small next to that bound and its relative error is large. A +//! `Σp` whose dynamic range exceeds `f32`'s — a nonzero entry that would +//! flush below `f32`'s smallest normal after the division — is refused rather +//! than silently zeroed. Both limits are `f32` limits of `CovHighD`; lifting +//! them is an `f64` sandwich in ndarray, not a second kernel here. `N` is a //! compile-time constant in `CovHighD`, while a [`crate::graph::Grid`]'s bus //! count is a runtime value: the caller picks `N` and the call panics when the //! decomposition disagrees with it. A runtime-sized sandwich would be a change @@ -50,8 +57,10 @@ pub const SYMMETRY_TOL: f64 = 1e-9; /// /// # Panics /// If `eig.n != N`, if `sigma_p.len() != N*N`, if any entry of `sigma_p` is -/// not finite, or if `sigma_p` is not symmetric within [`SYMMETRY_TOL`] -/// (relative to its largest entry). +/// not finite, if `sigma_p` is not symmetric within [`SYMMETRY_TOL`] +/// (relative to its largest entry), or if a nonzero entry of `sigma_p` is +/// below `f32::MIN_POSITIVE` relative to its largest entry (its dynamic range +/// does not fit `f32`, so that entry would silently become zero). pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f64) -> Vec { assert_eq!( eig.n, N, @@ -67,6 +76,15 @@ pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f .iter() .fold(0.0_f64, |m, v| m.max(v.abs())) .max(f64::MIN_POSITIVE); + if let Some(v) = sigma_p + .iter() + .find(|v| **v != 0.0 && (*v / scale).abs() < f64::from(f32::MIN_POSITIVE)) + { + panic!( + "sigma_p dynamic range exceeds f32: entry {v:e} next to max {scale:e} \ + would be zeroed by the f32 sandwich" + ); + } for i in 0..N { for j in 0..i { let d = (sigma_p[i * N + j] - sigma_p[j * N + i]).abs(); @@ -88,17 +106,25 @@ pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f let m = CovHighD::::from_symmetric_fn(|i, j| (l_plus[i * N + j] / l_scale) as f32); let s = CovHighD::::from_symmetric_fn(|i, j| (sigma_p[i * N + j] / scale) as f32); let out = s.sandwich(&m); - let back = scale * l_scale * l_scale; let mut dense = vec![0.0_f64; N * N]; for i in 0..N { for j in 0..N { - dense[i * N + j] = out.get(i, j) as f64 * back; + dense[i * N + j] = unscale(out.get(i, j), scale, l_scale); } } dense } +/// One sandwich entry back to `f64` units: `o · scale · l_scale²`. +/// +/// Applied factor by factor on the entry, never through a precomputed +/// `scale · l_scale²`: that product can overflow to infinity while the entry +/// itself does not, and an exact-zero entry would then become `0 · ∞ = NaN`. +fn unscale(o: f32, scale: f64, l_scale: f64) -> f64 { + f64::from(o) * scale * l_scale * l_scale +} + #[cfg(test)] mod tests { use super::*; @@ -269,6 +295,41 @@ mod tests { } } + /// FAILS IF: the result is rescaled through the combined factor + /// `scale · l_scale²`, which overflows here and turns an exact zero into + /// NaN. + #[test] + fn an_exact_zero_stays_zero_when_the_combined_factor_overflows() { + let (scale, l_scale) = (1e300_f64, 1e10_f64); + assert!( + (scale * l_scale * l_scale).is_infinite(), + "fixture must overflow" + ); + assert_eq!(unscale(0.0, scale, l_scale), 0.0); + // A nonzero entry whose true value is finite stays finite. + let v = unscale(1e-20, scale, l_scale); + let want = f64::from(1e-20_f32) * 1e300 * 1e20; + assert!( + v.is_finite() && (v - want).abs() <= 1e-12 * want.abs(), + "{v:e} vs {want:e}" + ); + } + + /// FAILS IF: an entry too small to survive the f32 narrowing is silently + /// zeroed instead of refused (a huge variance elsewhere, e.g. on a bus in + /// `L⁺`'s null space, sets the scale). + #[test] + #[should_panic(expected = "dynamic range exceeds f32")] + fn refuses_a_dynamic_range_wider_than_f32() { + let eig = symmetric_eigen(&grid().laplacian(), N); + let mut sp = vec![0.0; N * N]; + sp[0] = 1e100; + for k in 1..N { + sp[k * N + k] = 1.0; + } + let _ = angle_covariance::(&eig, &sp, TOL); + } + #[test] #[should_panic(expected = "non-finite")] fn refuses_a_non_finite_injection_covariance() { From 21376e45c76fc2d16c4c189b28d8297bd1a86709 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 4 Oct 2026 19:35:58 +0000 Subject: [PATCH 14/14] perturbation-sim: unscale in an order that cannot underflow or overflow early A fixed factor order fails one way or the other: tiny scale first underflows to zero before a large l_scale^2 would restore it, and the reverse overflows in the mirrored case. Each step now takes the factor that moves the running value toward 1, so only the final product can leave the representable range. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01X1YcYMRSFvfczXoP748wtB --- crates/perturbation-sim/src/angle_cov.rs | 54 +++++++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/crates/perturbation-sim/src/angle_cov.rs b/crates/perturbation-sim/src/angle_cov.rs index e03834c6b..838ad3031 100644 --- a/crates/perturbation-sim/src/angle_cov.rs +++ b/crates/perturbation-sim/src/angle_cov.rs @@ -121,8 +121,38 @@ pub fn angle_covariance(eig: &Eigen, sigma_p: &[f64], rel_tol: f /// Applied factor by factor on the entry, never through a precomputed /// `scale · l_scale²`: that product can overflow to infinity while the entry /// itself does not, and an exact-zero entry would then become `0 · ∞ = NaN`. +/// +/// No FIXED order is safe either: multiplying by a tiny `scale` first can +/// underflow to zero before a large `l_scale²` would have brought the value +/// back, and the reverse order overflows in the mirrored case. So each step +/// takes the factor that moves the running value toward 1 — the smallest +/// remaining factor while it is at least 1, the largest while it is below — +/// and only the final product can leave the representable range, which it +/// does only when the true result does. fn unscale(o: f32, scale: f64, l_scale: f64) -> f64 { - f64::from(o) * scale * l_scale * l_scale + let mut acc = f64::from(o); + let mut rest = [scale, l_scale, l_scale]; + let mut left = rest.len(); + while left > 0 { + let pick = (0..left) + .reduce(|a, b| { + let take_b = if acc.abs() >= 1.0 { + rest[b].abs() < rest[a].abs() + } else { + rest[b].abs() > rest[a].abs() + }; + if take_b { + b + } else { + a + } + }) + .unwrap_or(0); + acc *= rest[pick]; + rest.swap(pick, left - 1); + left -= 1; + } + acc } #[cfg(test)] @@ -315,6 +345,28 @@ mod tests { ); } + /// FAILS IF: `unscale` multiplies by a tiny `scale` before the large + /// `l_scale²`, so the intermediate underflows to zero although the final + /// value (~1e-30) is representable. The reverse case (huge `scale`, small + /// `l_scale`) must not overflow either. + #[test] + fn unscale_keeps_a_representable_result_when_one_factor_order_would_not() { + // L⁺ entries up to 1e150, sigma_p max 1e-310: entry (1e140/1e150)² = 1e-20. + let v = unscale(1e-20, 1e-310, 1e150); + assert!( + f64::from(1e-20_f32) * 1e-310 == 0.0, + "fixture must underflow scale-first" + ); + assert!(v != 0.0 && (v - 1e-30).abs() <= 1e-6 * 1e-30, "{v:e}"); + // Mirror: scale-first overflows (1e10 * 1e305), the true value is 1e-5. + assert!( + (1e10_f64 * 1e305).is_infinite(), + "fixture must overflow scale-first" + ); + let v = unscale(1e10, 1e305, 1e-160); + assert!(v.is_finite() && (v - 1e-5).abs() <= 1e-6 * 1e-5, "{v:e}"); + } + /// FAILS IF: an entry too small to survive the f32 narrowing is silently /// zeroed instead of refused (a huge variance elsewhere, e.g. on a bus in /// `L⁺`'s null space, sets the scale).