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..8d450d1a6 --- /dev/null +++ b/.claude/board/entries/2026-09-29-grouped-cross-moments-fold.md @@ -0,0 +1,19 @@ +# 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: `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_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 `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 new file mode 100644 index 000000000..79149a567 --- /dev/null +++ b/.claude/board/entries/2026-09-29-grouped-moments-fold-anova-over-masks.md @@ -0,0 +1,18 @@ +# 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: `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`. 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 +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 `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-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 b0eaf97ad..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,9 @@ 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) | | 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) | diff --git a/crates/jc/src/reliability.rs b/crates/jc/src/reliability.rs index c73449ee2..9f6b1de5c 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_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(); 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 d72775926..39eec4a72 100644 --- a/crates/jc/src/stats.rs +++ b/crates/jc/src/stats.rs @@ -73,7 +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 crate::reliability::{all_finite, mean, pearson, pearson_from_centered}; +use ndarray::simd::{CrossPowerSums, PowerSums}; use std::collections::BTreeSet; // ─────────────────────────── local helpers ─────────────────────────── @@ -103,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_power_sums`]. +#[inline] +fn sample_cov_tail(sxy: f64, n: usize) -> Option { + (n >= 2).then(|| sxy / (n as f64 - 1.0)) } // ───────────────────── incomplete beta (p-values) ───────────────────── @@ -877,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_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 +/// 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; @@ -947,6 +959,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_power_sums`]. +fn eta_from_ss(ss_b: f64, ss_t: f64) -> Option { if ss_t == 0.0 || !ss_t.is_finite() { return None; } @@ -1134,12 +1152,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_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 { 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 +1191,222 @@ 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_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. + +/// 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 +/// `PowerSums` is exact under. +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; + } + 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(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 + } + 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_power_sums_i32` (or +/// `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_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_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]]; +/// let moments = [ +/// 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_power_sums(&moments).unwrap()); +/// assert!((a.f - b.f).abs() < 1e-12 && (a.p - b.p).abs() < 1e-12); +/// ``` +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) +} + +/// η² 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_power_sums(groups: &[PowerSums]) -> Option { + let (ss_b, _, ss_t, _, _) = one_way_ss_from_power_sums(groups)?; + 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_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. + +/// 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. +/// +/// 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: &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(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)) +} + +/// 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_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_power_sums(&m).unwrap() - 1.0).abs() < 1e-12); +/// ``` +pub fn pearson_from_cross_power_sums(m: &CrossPowerSums) -> 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_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; + 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_power_sums(m: &CrossPowerSums) -> 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_power_sums(m: &CrossPowerSums) -> 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 @@ -2246,4 +2487,716 @@ 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_power_sums_i32` — the kernel + // `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 + // 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 power_sums_equivalence { + use super::super::*; + use ndarray::simd::masked_group_power_sums_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![PowerSums::default(); k as usize]; + masked_group_power_sums_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_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"), + } + 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![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)); + 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![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"); + } + } + assert_eq!(merged, whole, "chunks={chunks}"); + assert_eq!( + anova_from_power_sums(&merged), + anova_from_power_sums(&whole) + ); + } + assert!(anova_from_power_sums(&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_power_sums(&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]| PowerSums { + n, + sum: xs.iter().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); + 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]), PowerSums::default(), 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_power_sums(&moments), None, "{what}: moments"); + assert_eq!( + eta_squared(&groups), + 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_power_sums(&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(PowerSums { + n: xs.len() as u64, + sum: xs.iter().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::>()); + } + // 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_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, + "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), + } + } + } + + // ─────── 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_power_sums`). + + mod cross_power_sums_equivalence { + use super::super::*; + use crate::reliability::pearson; + use ndarray::simd::masked_group_cross_power_sums_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![CrossPowerSums::default(); k as usize]; + masked_group_cross_power_sums_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: &CrossPowerSums, what: &str) { + 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_power_sums(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_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_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(); + 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![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![CrossPowerSums::default(); 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_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 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(); + 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 + /// 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 = 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]| 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_power_sums(&m(&[5, 5, 5], &[1, 2, 3])), + None + ); + 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_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_power_sums(&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 = 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(); + 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_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_power_sums(&m).unwrap(); + assert!(rel(r2, r_exact * r_exact) < 1e-15, "B={b}: R²"); + 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_power_sums(&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_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}" + ); + } + } + + /// 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 power_sums_past_the_bound_are_refused() { + let v = (1i128 << 95) + 1; + let huge = CrossPowerSums { + n: 1 << 33, + sum_x: 0, + sum_y: 0, + sum_x_sq: v as u128, + sum_y_sq: v as u128, + sum_xy: v, + }; + 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 5ef160695..60b783a79 100644 --- a/crates/lance-graph-mask-risc/src/exec.rs +++ b/crates/lance-graph-mask-risc/src/exec.rs @@ -29,15 +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_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_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, + 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::{ @@ -1339,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 @@ -1411,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"), @@ -1572,6 +1580,10 @@ 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::GroupPowerSumsI32 { .. }, Out::PowerSums(o)) => o.fill(PowerSums::default()), + (Terminal::GroupCrossPowerSumsI32 { .. }, Out::CrossPowerSums(o)) => { + o.fill(CrossPowerSums::default()) + } _ => {} } let n_rows = planes.n_rows; @@ -1983,6 +1995,68 @@ pub fn execute_compiled( } } } + 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::PowerSums(o) = &mut out { + let m = clip(read(planes, &slots, mask, t), edge, &mut eb, false); + let v = lane_i32(planes, val, t); + match key { + GroupKey::Lane(k) => { + masked_group_power_sums_i32(m, lane_u32(planes, k, t), v, o) + } + 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_power_sums_i32_pair( + m, + lane_u32(planes, hi, t), + lane_u32(planes, lo, t), + stride, + v, + o, + ), + } + } + } + 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::CrossPowerSums(o) = &mut out { + 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) => { + masked_group_cross_power_sums_i32(m, lane_u32(planes, k, t), xs, ys, o) + } + GroupKey::Via { fk, key } => masked_group_cross_power_sums_i32_via( + m, + lane_u32(planes, fk, t), + foreign_lane_u32(foreign, key), + 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, + ) + } + } + } + } 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 +2091,8 @@ pub fn execute_compiled( Terminal::CountKeyRunsU32 { .. } => Value::Count(runs + run_carry.finish()), Terminal::GroupSumI32 { .. } | Terminal::GroupSumViaI32 { .. } => Value::GroupSummed, Terminal::GroupReduce { .. } => Value::GroupReduced, + 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 5cf112d6b..0745cf2f4 100644 --- a/crates/lance-graph-mask-risc/src/ir.rs +++ b/crates/lance-graph-mask-risc/src/ir.rs @@ -421,6 +421,45 @@ 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::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 + /// 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 `PowerSums::default()` before + /// the first tile, so an empty group reads `n == 0`. + 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::CrossPowerSums` buffer as `n, Σx, Σy, Σx², Σy², Σxy` + /// ([`ndarray::simd::CrossPowerSums`]). The bivariate member of the + /// [`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). + GroupCrossPowerSumsI32 { + mask: Operand, + key: GroupKey, + x: u16, + y: u16, + }, } /// Where a [`Terminal::GroupReduce`] reads each row's group. @@ -612,6 +651,8 @@ impl Program { | Terminal::GroupSumI32 { mask, .. } | Terminal::GroupSumViaI32 { mask, .. } | Terminal::GroupReduce { mask, .. } + | Terminal::GroupPowerSumsI32 { mask, .. } + | Terminal::GroupCrossPowerSumsI32 { 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..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, 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 e1fb7055b..394b04d29 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::{CrossPowerSums, PowerSums}; use crate::value::{ExecError, LaneKind, Out, Value}; use crate::words_for; @@ -34,6 +35,8 @@ pub(crate) enum OutShape { I32(usize), I64(usize), Mask(usize), + PowerSums(usize), + CrossPowerSums(usize), } /// [`OutShape`] of a borrowed `out` — the caller keeps `out` itself to write @@ -44,6 +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::PowerSums(v) => OutShape::PowerSums(v.len()), + Out::CrossPowerSums(v) => OutShape::CrossPowerSums(v.len()), } } @@ -516,7 +521,10 @@ pub(crate) fn validate( found: len, }), OutShape::I32(_) => Ok(()), - OutShape::I64(_) | OutShape::Mask(_) => Err(ExecError::BlendNeedsOut), + OutShape::I64(_) + | OutShape::Mask(_) + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::BlendNeedsOut), } } Terminal::ScatterOrU32 { @@ -544,9 +552,11 @@ pub(crate) fn validate( expected: want, found: len, }), - OutShape::None | OutShape::I32(_) | OutShape::I64(_) => { - Err(ExecError::TerminalNeedsOut { what }) - } + OutShape::None + | OutShape::I32(_) + | OutShape::I64(_) + | OutShape::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { what }), } } Terminal::GroupSumI32 { mask, key, val } => { @@ -559,11 +569,14 @@ 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::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { + what: "GroupSumI32", + }), } } Terminal::GroupSumViaI32 { mask, fk, key, val } => { @@ -577,27 +590,20 @@ 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::PowerSums(_) + | OutShape::CrossPowerSums(_) => 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,12 +618,71 @@ 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::PowerSums(_) + | OutShape::CrossPowerSums(_) => Err(ExecError::TerminalNeedsOut { + what: "GroupReduce", + }), + } + } + Terminal::GroupPowerSumsI32 { 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::PowerSums(len) if len >= 1 => Ok(()), + _ => Err(ExecError::TerminalNeedsOut { + what: "GroupPowerSumsI32", + }), + } + } + Terminal::GroupCrossPowerSumsI32 { 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::CrossPowerSums(len) if len >= 1 => Ok(()), + _ => Err(ExecError::TerminalNeedsOut { + what: "GroupCrossPowerSumsI32", + }), + } + } + } +} + +/// 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) } } } @@ -640,6 +705,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 +1135,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 +1157,67 @@ pub fn reference_execute_into( } Value::GroupReduced } + 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. + for x in o.iter_mut() { + *x = PowerSums::default(); + } + 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) 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 + // an oracle bug, never a data condition. + slot.sum = i64::try_from(s).expect("validated 2^32-row bound keeps Σx in i64"); + } + } + Value::GroupPowerSums + } + 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) { + 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 = CrossPowerSums { + n: u64::try_from(a[0]).expect("a count is non-negative"), + sum_x: narrow(a[1]), + sum_y: narrow(a[2]), + 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], + }; + } + } + Value::GroupCrossPowerSums + } 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..30dff5bf6 100644 --- a/crates/lance-graph-mask-risc/src/value.rs +++ b/crates/lance-graph-mask-risc/src/value.rs @@ -5,6 +5,13 @@ //! rejects, or reject with a different reason. use crate::ir::Operand; +/// The per-group `(n, Σx, Σy, Σx², Σy², Σxy)` accumulator +/// [`Out::CrossPowerSums`] carries — re-exported for the same reason. +pub use ndarray::simd::CrossPowerSums; +/// 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; /// What a [`crate::Program`] produced. #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -32,6 +39,14 @@ pub enum Value { /// written, one slot per group; see [`crate::GroupFold::seed`] for what /// an empty group holds. GroupReduced, + /// [`crate::Terminal::GroupPowerSumsI32`]: the caller's `Out::PowerSums` + /// buffer was written, one [`PowerSums`] per group; an empty group + /// holds [`PowerSums::default()`] (`n == 0`). + GroupPowerSums, + /// [`crate::Terminal::GroupCrossPowerSumsI32`]: the caller's + /// `Out::CrossPowerSums` buffer was written, one [`CrossPowerSums`] per + /// group; an empty group holds [`CrossPowerSums::default()`]. + GroupCrossPowerSums, /// [`crate::Terminal::MaskedStridedGroupSum`]: the widened sum, or `None` /// when it does not fit an `i64` (never a wrapped value). StridedSum(Option), @@ -52,6 +67,12 @@ pub enum Out<'a> { /// [`crate::Terminal::ScatterOrU32`]'s destination, `words_for(out_rows)` /// long. Mask(&'a mut [u64]), + /// [`crate::Terminal::GroupPowerSumsI32`]'s destination — one + /// [`PowerSums`] per group, its length IS the group universe `K`. + PowerSums(&'a mut [PowerSums]), + /// [`crate::Terminal::GroupCrossPowerSumsI32`]'s destination — one + /// [`CrossPowerSums`] per group, its length IS the group universe `K`. + CrossPowerSums(&'a mut [CrossPowerSums]), } /// The lane width a predicate or terminal expects, for [`ExecError::LaneKind`]. @@ -195,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 }, } 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 5774bd965..e131c98ca 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_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, 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 { @@ -1568,3 +1568,334 @@ fn pair_key_drops_a_minor_key_at_stride() { "rows with lo >= stride must never be counted anywhere" ); } + +/// 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²`. +/// +/// 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_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; + 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::GroupPowerSumsI32 { + 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 = PowerSums { + n: 7, + sum: -7, + sum_sq: 7, + }; + let mut got_out = vec![dirty; groups]; + let got = execute_into( + &p, + &planes, + &foreign, + &mut scratch, + Out::PowerSums(&mut got_out), + ) + .expect("runs"); + let mut want_out = vec![dirty; groups]; + let want = reference_execute_into(&p, &planes, &foreign, Out::PowerSums(&mut want_out)) + .expect("oracle runs"); + 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); + 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: `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; + 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![PowerSums::default(); 4]; + // `val` must be an I32 lane (lane 1 is U32). + assert!(matches!( + run( + Terminal::GroupPowerSumsI32 { + mask: S0, + key: GroupKey::Lane(3), + val: 1 + }, + Out::PowerSums(&mut sink) + ), + Err(ExecError::LaneKind { .. }) + )); + // The key must be a U32 lane (lane 2 is I32). + assert!(matches!( + run( + Terminal::GroupPowerSumsI32 { + mask: S0, + key: GroupKey::Lane(2), + val: 2 + }, + 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::PowerSums(&mut [])] { + assert_eq!( + run( + Terminal::GroupPowerSumsI32 { + mask: S0, + key: GroupKey::Lane(3), + val: 2 + }, + out + ), + Err(ExecError::TerminalNeedsOut { + what: "GroupPowerSumsI32" + }) + ); + } +} + +/// 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::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_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; + 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::GroupCrossPowerSumsI32 { + 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 = CrossPowerSums { + n: 3, + sum_xy: -3, + ..CrossPowerSums::default() + }; + let mut got_out = vec![dirty; groups]; + let got = execute_into( + &p, + &planes, + &foreign, + &mut scratch, + Out::CrossPowerSums(&mut got_out), + ) + .expect("runs"); + let mut want_out = vec![dirty; groups]; + let want = + reference_execute_into(&p, &planes, &foreign, Out::CrossPowerSums(&mut want_out)) + .expect("oracle runs"); + 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); + 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::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; + 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::GroupCrossPowerSumsI32 { + 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![CrossPowerSums::default(); 4]; + for (x, y) in [(1u16, 4u16), (2, 1)] { + assert!( + matches!( + run(term(x, y), Out::CrossPowerSums(&mut sink)), + Err(ExecError::LaneKind { .. }) + ), + "x={x} y={y}" + ); + } + let mut univariate = vec![PowerSums::default(); 4]; + for out in [ + Out::None, + Out::PowerSums(&mut univariate), + Out::CrossPowerSums(&mut []), + ] { + assert_eq!( + run(term(2, 4), out), + Err(ExecError::TerminalNeedsOut { + what: "GroupCrossPowerSumsI32" + }) + ); + } + // 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"); + 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 96da4eadb..d15aaa3bc 100644 --- a/crates/lance-graph-quack/tests/duckdb_differential.rs +++ b/crates/lance-graph-quack/tests/duckdb_differential.rs @@ -264,6 +264,10 @@ fn run_query(id: &str, planes: &Planes<'_>, filter: Filter, agg: Agg) -> (String Value::StridedSum(_) => { panic!("case {id}: no MaskedStridedGroupSum 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(); ( 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..838ad3031 --- /dev/null +++ b/crates/perturbation-sim/src/angle_cov.rs @@ -0,0 +1,409 @@ +//! 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`. 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). 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 +//! 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`, if any entry of `sigma_p` is +/// 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, + "decomposition has {} buses, CovHighD has N = {N}", + 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())) + .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(); + 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 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 mut dense = vec![0.0_f64; N * N]; + for i in 0..N { + for j in 0..N { + 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`. +/// +/// 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 { + 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)] +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" + ); + } + + /// `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 { + 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; + } + } + 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}"); + 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}"); + } + + /// 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}"); + } + } + + /// 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: `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). + #[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() { + 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() { + 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;