From 695bb822efee965b32a4f4db8ee1c7635e142e03 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Thu, 10 Sep 2026 07:08:18 +0200 Subject: [PATCH] perf!: sparse Cholesky with an AMD ordering for the joint MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 745 ms -> 1.11 ms on the fixture #52 was opened about. The joint precision matrix is 0.19% dense at scale and gets sparser as the history grows. We allocated all n^2 entries — 31 MB at n = 1976, 128 MB at ustat's ~4000 appearances — filled 99.8% of it with zeros, and ran an O(n^3) factorisation over the whole thing. Two measurements shaped the fix, and the first killed the plan #52 proposed. **Ordering alone does nothing to a dense factorisation.** Its inner loops run over every k whether the entry is zero or not, so a permutation changes which entries are zero and not how many multiplications happen. A 700x700 banded matrix at 0.43% density: 30.196 ms in band order, 29.544 ms under a scramble that destroyed the band. Identical, as the flop count says it must be. #52's step 1 — "reorder with AMD, keep our own Cholesky, and measure" — could not have worked, and measuring said so before any of it was written. **Sparsity and AMD together are worth four orders of magnitude.** Symbolic factorisation on the n = 1976 fixture, against 2.572e9 dense flops: sparse in natural order needs 5.597e7 (46x), sparse under AMD needs 8.656e4 — 29,710x. AMD is worth 646x on top of sparsity and nothing without it. Natural order fills in badly for exactly the reason #52 predicted about bandwidth: nnz(L) is 292,437 against A's 7,504, because a competitor idle from slice 0 to slice 75 links across the whole matrix. Measured end to end, factorising through `History::joint`: n = 480 215 us (bench: 9.11 ms -> 167 us, 54x) n = 1976 1.112 ms (was ~745 ms, 670x) n = 7800 4.616 ms (dense would be 1.58e11 flops) Scaling is near-linear now rather than cubic: 16x the variables costs 21x the time, where dense would cost 4096x. The factorisation is the up-looking sparse Cholesky of Davis's *Direct Methods for Sparse Linear Systems*, written here rather than taken from a crate. The scouting in #52 still holds and got one addition: `feral` itself pulls `pulp`, so it has the same runtime CPU-dispatch problem that ruled out `faer` — results could differ between an AVX-512 host and an AVX2 one, the drift the libm-over-std decision was made to avoid. `sprs-ldl` is still LGPL and `nalgebra-sparse` still disclaims fill-reduction in its own docs. Only the ordering is a dependency: `feral-amd`, two crates, both `#![forbid(unsafe_code)]`. The matrix is accumulated into a `BTreeMap`, not a hash map: the iteration order becomes the summation order, and a hash map's varies per process. `tests/cross_process_determinism.rs` exists because that has bitten before. `whiten` returns its result in the permuted order and leaves it there — a dot product does not care, as long as both operands were permuted the same way — so `bilinear` is unchanged. Correctness: the existing analytic goldens are 2x2 and 3x3, too small to permute or fill in, so they could not have caught a symbolic-pass bug. `agrees_with_a_dense_reference_on_random_sparse_systems` checks every bilinear form against a deliberately naive dense factorisation that shares no code with the thing it is checking, on chain-plus-long-range matrices up to n = 60. Closes #52. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- Cargo.toml | 3 + src/history.rs | 43 +++- src/joint.rs | 452 ++++++++++++++++++++++++++++++---- tests/sparsity_measurement.rs | 178 +++++++++++++ 4 files changed, 619 insertions(+), 57 deletions(-) create mode 100644 tests/sparsity_measurement.rs diff --git a/Cargo.toml b/Cargo.toml index 88fae9b..21e951d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -47,12 +47,15 @@ harness = false [dependencies] approx = { version = "0.5.1", optional = true } +feral-amd = "0.2" libm = "0.2.16" rayon = { version = "1", optional = true } smallvec = "1" [features] approx = ["dep:approx"] +# Exposes the joint sparsity pattern for the #52 measurement. Test-only. +measure-sparsity = [] rayon = ["dep:rayon"] [dev-dependencies] diff --git a/src/history.rs b/src/history.rs index 9b0c395..e243b10 100644 --- a/src/history.rs +++ b/src/history.rs @@ -386,8 +386,8 @@ pub(crate) struct CompetitorConfig { /// The joint precision over a history's appearances, with the maps needed to /// address a competitor either at their latest appearance or at a given slice. struct TimeExpanded { - /// Row-major precision matrix over appearances. - lambda: Vec, + /// Precision matrix over appearances, accumulated sparsely. + lambda: crate::joint::SymmetricBuilder, /// `(row, slice)` of each competitor's latest appearance. latest: HashMap, /// Row of each `(competitor, slice)` appearance. @@ -1415,6 +1415,23 @@ impl, O: Observer, K: Eq + Hash + Clone> History (usize, Vec>) { + let te = self.time_expanded_joint(); + let n = te.width; + let mut adj = vec![std::collections::HashSet::new(); n]; + for (i, j) in te.lambda.pattern() { + if i != j { + adj[i].insert(j); + } + } + (n, adj) + } + fn time_expanded_joint(&self) -> TimeExpanded { let mut latest: HashMap = HashMap::new(); let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new(); @@ -1455,17 +1472,23 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History, +} + +impl SymmetricBuilder { + pub(crate) fn new() -> Self { + Self::default() + } + + /// Add `value` to entry `(row, col)`. Both triangles must be supplied. + pub(crate) fn add(&mut self, row: usize, col: usize, value: f64) { + *self.entries.entry((row, col)).or_insert(0.0) += value; + } + + /// The `(row, col)` positions that hold a nonzero. For the #52 measurement. + #[cfg(feature = "measure-sparsity")] + pub(crate) fn pattern(&self) -> impl Iterator + '_ { + self.entries + .iter() + .filter(|(_, v)| **v != 0.0) + .map(|(&rc, _)| rc) + } +} /// A factorised symmetric positive-definite matrix, reusable across queries. pub(crate) struct Cholesky { - /// Lower triangle of `L`, row-major `n * n`. The upper triangle is - /// leftover scratch from the factorisation and is never read. - l: Vec, n: usize, + /// `inv[old] = new`: where each original row sits after the AMD reorder. + inv: Vec, + /// `L` in compressed-column form, permuted. Within a column the diagonal + /// is first and the rest ascend by row. + col_ptr: Vec, + row_idx: Vec, + val: Vec, } impl Cholesky { - /// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`. - /// - /// `a` is consumed as scratch. + /// Factorise the accumulated matrix into `L L^T`, under a fill-reducing + /// permutation. /// /// Returns `None` if the matrix is not positive-definite, which for a /// precision matrix means the model is improper — a competitor with - /// neither a proper prior nor any evidence. - pub(crate) fn factor(mut a: Vec, n: usize) -> Option { - debug_assert_eq!(a.len(), n * n); + /// neither a proper prior nor any evidence — or if the ordering fails. + pub(crate) fn factor(built: SymmetricBuilder, n: usize) -> Option { + if n == 0 { + return Some(Self { + n: 0, + inv: Vec::new(), + col_ptr: vec![0], + row_idx: Vec::new(), + val: Vec::new(), + }); + } - for j in 0..n { - let mut d = a[j * n + j]; - for k in 0..j { - d -= a[j * n + k] * a[j * n + k]; + let inv = Self::amd_permutation(n, &built)?; + + // Upper triangle of the permuted matrix, column-major: column `c` + // holds the rows `r <= c`. Exactly one of a symmetric pair survives + // the `r <= c` filter, so nothing is double-counted. + let mut cols: Vec> = vec![Vec::new(); n]; + for (&(old_r, old_c), &v) in &built.entries { + if v == 0.0 { + continue; } + let (r, c) = (inv[old_r], inv[old_c]); + if r <= c { + cols[c].push((r, v)); + } + } + let mut a_ptr = Vec::with_capacity(n + 1); + let mut a_row = Vec::new(); + let mut a_val = Vec::new(); + a_ptr.push(0usize); + for col in &mut cols { + col.sort_unstable_by_key(|&(r, _)| r); + for &(r, v) in col.iter() { + a_row.push(r); + a_val.push(v); + } + a_ptr.push(a_row.len()); + } + + let parent = Self::etree(n, &a_ptr, &a_row); + + // Symbolic pass: how many entries each column of L will hold. Running + // `ereach` per column costs O(nnz(L)) in total, which is the same order + // as the numeric pass it sizes. + let mut counts = vec![0usize; n]; + let mut stack = vec![0usize; n]; + let mut mark = vec![false; n]; + for k in 0..n { + let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark); + for &i in &stack[top..] { + counts[i] += 1; + } + counts[k] += 1; // the diagonal + } + + let mut col_ptr = Vec::with_capacity(n + 1); + col_ptr.push(0usize); + for &c in &counts { + col_ptr.push(col_ptr[col_ptr.len() - 1] + c); + } + let nnz = col_ptr[n]; + let mut row_idx = vec![0usize; nnz]; + let mut val = vec![0.0f64; nnz]; + + // `next[i]` is the slot column `i` will fill next. Column `i`'s + // diagonal lands first, at `col_ptr[i]`, because nothing is written to + // a column before its own iteration. + let mut next: Vec = col_ptr[..n].to_vec(); + let mut x = vec![0.0f64; n]; + + for k in 0..n { + let top = Self::ereach(k, &a_ptr, &a_row, &parent, &mut stack, &mut mark); + + for p in a_ptr[k]..a_ptr[k + 1] { + if a_row[p] <= k { + x[a_row[p]] = a_val[p]; + } + } + + let mut d = x[k]; + x[k] = 0.0; + + for &i in &stack[top..] { + let lki = x[i] / val[col_ptr[i]]; + x[i] = 0.0; + for p in col_ptr[i] + 1..next[i] { + x[row_idx[p]] -= val[p] * lki; + } + d -= lki * lki; + let p = next[i]; + next[i] += 1; + row_idx[p] = k; + val[p] = lki; + } + // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here // too, and a negated comparison would let it through as "not // positive". if d.is_nan() || d <= 0.0 { return None; } - let d = d.sqrt(); - a[j * n + j] = d; - - for i in j + 1..n { - let mut s = a[i * n + j]; - for k in 0..j { - s -= a[i * n + k] * a[j * n + k]; - } - a[i * n + j] = s / d; - } + let p = next[k]; + next[k] += 1; + row_idx[p] = k; + val[p] = d.sqrt(); } - Some(Self { l: a, n }) + Some(Self { + n, + inv, + col_ptr, + row_idx, + val, + }) } - /// Whiten a contrast: `y = L^-1 b`. + /// AMD fill-reducing order, as `inv[old] = new`. + fn amd_permutation(n: usize, built: &SymmetricBuilder) -> Option> { + let mut cols: Vec> = vec![Vec::new(); n]; + for (&(r, c), &v) in &built.entries { + if v != 0.0 { + cols[c].push(i32::try_from(r).ok()?); + } + } + let mut col_ptr = Vec::with_capacity(n + 1); + let mut row_idx = Vec::new(); + col_ptr.push(0i32); + for (j, col) in cols.iter_mut().enumerate() { + col.push(i32::try_from(j).ok()?); + col.sort_unstable(); + col.dedup(); + row_idx.extend_from_slice(col); + col_ptr.push(i32::try_from(row_idx.len()).ok()?); + } + + let pattern = feral_amd::CscPattern::new(n, &col_ptr, &row_idx)?; + // `perm[new] = old`; we want the inverse. + let perm = feral_amd::amd_order(&pattern).ok()?; + let mut inv = vec![0usize; n]; + for (new, &old) in perm.iter().enumerate() { + inv[usize::try_from(old).ok()?] = new; + } + Some(inv) + } + + /// Elimination tree of the upper-triangular pattern. `usize::MAX` is "no + /// parent", i.e. a root. + fn etree(n: usize, col_ptr: &[usize], row_idx: &[usize]) -> Vec { + let mut parent = vec![usize::MAX; n]; + let mut ancestor = vec![usize::MAX; n]; + for k in 0..n { + for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] { + let mut i = row; + while i != usize::MAX && i < k { + let next = ancestor[i]; + ancestor[i] = k; + if next == usize::MAX { + parent[i] = k; + } + i = next; + } + } + } + parent + } + + /// Nonzero pattern of row `k` of `L`, written into `stack[top..n]` in + /// topological order. Returns `top`. + /// + /// `stack` is used from both ends — a scratch region from `0` while walking + /// each path up the tree, and the result from `n` downwards. They cannot + /// collide because every node is pushed at most once across the whole call. + fn ereach( + k: usize, + col_ptr: &[usize], + row_idx: &[usize], + parent: &[usize], + stack: &mut [usize], + mark: &mut [bool], + ) -> usize { + let n = mark.len(); + let mut top = n; + mark[k] = true; + for &row in &row_idx[col_ptr[k]..col_ptr[k + 1]] { + let mut i = row; + if i > k { + continue; + } + let mut len = 0usize; + while i != usize::MAX && !mark[i] { + stack[len] = i; + len += 1; + mark[i] = true; + i = parent[i]; + } + // Reverse the path onto the output end, so the result stays in + // topological order overall. + while len > 0 { + len -= 1; + top -= 1; + stack[top] = stack[len]; + } + } + for &i in &stack[top..] { + mark[i] = false; + } + mark[k] = false; + top + } + + /// Whiten a contrast: `y = L^-1 P b`. /// /// The point of the result is the dot product, not the vector: for two /// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs. + /// + /// The result is in the permuted order, and stays there — a dot product + /// does not care, as long as both operands were permuted the same way. pub(crate) fn whiten(&self, b: &[f64]) -> Vec { debug_assert_eq!(b.len(), self.n); let n = self.n; - let mut y = b.to_vec(); - for i in 0..n { - // Folded from `y[i]` rather than summed and subtracted once, so the - // accumulation order matches a plain substitution loop exactly. - let row = &self.l[i * n..i * n + i]; - let s = row - .iter() - .zip(&y[..i]) - .fold(y[i], |acc, (l, v)| acc - l * v); - y[i] = s / self.l[i * n + i]; + let mut y = vec![0.0f64; n]; + for (old, &v) in b.iter().enumerate() { + y[self.inv[old]] = v; + } + for j in 0..n { + y[j] /= self.val[self.col_ptr[j]]; + let yj = y[j]; + for p in self.col_ptr[j] + 1..self.col_ptr[j + 1] { + y[self.row_idx[p]] -= self.val[p] * yj; + } } y } @@ -98,11 +350,24 @@ pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 { mod tests { use super::*; + /// Factorise a dense row-major matrix, for the goldens below. + fn dense(a: &[f64], n: usize) -> Option { + let mut b = SymmetricBuilder::new(); + for i in 0..n { + for j in 0..n { + if a[i * n + j] != 0.0 { + b.add(i, j, a[i * n + j]); + } + } + } + Cholesky::factor(b, n) + } + /// `[[4, 1], [1, 3]] z = [1, 2]` has `z = [1/11, 7/11]`, so the quadratic /// form `b^T A^-1 b` is `1 * 1/11 + 2 * 7/11 = 15/11`. #[test] fn reproduces_a_known_quadratic_form() { - let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap(); + let c = dense(&[4.0, 1.0, 1.0, 3.0], 2).unwrap(); let y = c.whiten(&[1.0, 2.0]); assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12); } @@ -113,8 +378,8 @@ mod tests { fn recovers_the_inverse_diagonal() { // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is // [0.75, 1.0, 0.75]. - let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; - let c = Cholesky::factor(a, 3).unwrap(); + let a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; + let c = dense(&a, 3).unwrap(); for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() { let mut e = vec![0.0; 3]; e[i] = 1.0; @@ -127,8 +392,8 @@ mod tests { #[test] fn recovers_an_off_diagonal_covariance() { // Same A; (A^-1)_{0,1} = 0.5. - let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; - let c = Cholesky::factor(a, 3).unwrap(); + let a = [2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; + let c = dense(&a, 3).unwrap(); let y0 = c.whiten(&[1.0, 0.0, 0.0]); let y1 = c.whiten(&[0.0, 1.0, 0.0]); assert!((bilinear(&y0, &y1) - 0.5).abs() < 1e-12); @@ -138,15 +403,108 @@ mod tests { /// A variance can never come out negative, because it is a sum of squares. #[test] fn a_quadratic_form_is_never_negative() { - let a = vec![1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12]; - let c = Cholesky::factor(a, 2).unwrap(); + let a = [1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12]; + let c = dense(&a, 2).unwrap(); let y = c.whiten(&[1.0, -1.0]); assert!(bilinear(&y, &y) >= 0.0); } + /// Against an independent dense reference, on random sparse SPD matrices. + /// + /// The goldens above are 2x2 and 3x3 — small enough that AMD does nothing + /// and no fill-in occurs, so they cannot catch a symbolic-pass bug. This + /// builds matrices big enough to permute and fill in, and checks every + /// bilinear form against a textbook dense factorisation of the *same* + /// matrix in its original order. + #[test] + fn agrees_with_a_dense_reference_on_random_sparse_systems() { + /// Dense Cholesky and quadratic form, deliberately naive: this is the + /// reference, so it must not share code with what it is checking. + fn dense_quadratic_form(a: &[f64], n: usize, b: &[f64], c: &[f64]) -> f64 { + let mut l = a.to_vec(); + for j in 0..n { + let mut d = l[j * n + j]; + for k in 0..j { + d -= l[j * n + k] * l[j * n + k]; + } + let d = d.sqrt(); + l[j * n + j] = d; + for i in j + 1..n { + let mut sum = l[i * n + j]; + for k in 0..j { + sum -= l[i * n + k] * l[j * n + k]; + } + l[i * n + j] = sum / d; + } + } + let solve = |rhs: &[f64]| -> Vec { + let mut y = rhs.to_vec(); + for i in 0..n { + for k in 0..i { + y[i] -= l[i * n + k] * y[k]; + } + y[i] /= l[i * n + i]; + } + y + }; + let (yb, yc) = (solve(b), solve(c)); + yb.iter().zip(&yc).map(|(x, y)| x * y).sum() + } + + // A cheap deterministic generator; no dependency, and reproducible. + let mut seed = 0x2545_F491_4F6C_DD1Du64; + let mut rand = move || { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + (seed >> 11) as f64 / (1u64 << 53) as f64 + }; + + for n in [7usize, 23, 60] { + let mut a = vec![0.0f64; n * n]; + // A chain plus scattered long-range couplings: the shape of a + // time-expanded joint, where a competitor's drift link can span + // the whole matrix. + for i in 0..n { + a[i * n + i] = 4.0 + rand(); + if i + 1 < n { + let v = -(0.5 + rand() * 0.5); + a[i * n + i + 1] = v; + a[(i + 1) * n + i] = v; + } + } + for step in 0..n / 3 { + let i = (step * 7) % n; + let j = (step * 29 + 3) % n; + if i != j { + let v = -(0.1 + rand() * 0.2); + a[i * n + j] = v; + a[j * n + i] = v; + // Keep it diagonally dominant, hence positive-definite. + a[i * n + i] += 0.6; + a[j * n + j] += 0.6; + } + } + + let sparse = dense(&a, n).expect("spd"); + + for trial in 0..8 { + let b: Vec = (0..n).map(|_| rand() * 2.0 - 1.0).collect(); + let c: Vec = (0..n).map(|_| rand() * 2.0 - 1.0).collect(); + let got = bilinear(&sparse.whiten(&b), &sparse.whiten(&c)); + let want = dense_quadratic_form(&a, n, &b, &c); + assert!( + (got - want).abs() <= 1e-10 * want.abs().max(1.0), + "n={n} trial={trial}: sparse {got} vs dense {want}" + ); + } + } + } + + /// A permutation must not change which matrices are rejected. #[test] fn rejects_a_non_positive_definite_matrix() { // Singular: the second row is a multiple of the first. - assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none()); + assert!(dense(&[1.0, 2.0, 2.0, 4.0], 2).is_none()); } } diff --git a/tests/sparsity_measurement.rs b/tests/sparsity_measurement.rs new file mode 100644 index 0000000..169f596 --- /dev/null +++ b/tests/sparsity_measurement.rs @@ -0,0 +1,178 @@ +//! What a sparse factorisation of the joint would actually buy (#52). +//! +//! Run explicitly: +//! +//! ```text +//! cargo test --release --features approx,measure-sparsity \ +//! --test sparsity_measurement -- --ignored --nocapture +//! ``` +//! +//! The whole file is gated: it reaches for the joint's sparsity pattern, which +//! is exposed only under `measure-sparsity`. +#![cfg(feature = "measure-sparsity")] + +use std::collections::HashSet; + +use trueskill_tt::{ConvergenceOptions, History}; + +/// A history shaped like the issue's fixture: many slices, scored duels, +/// competitors reappearing across slices so the drift links are long. +fn fitted(slices: i64, duels: usize, competitors: usize) -> History { + let mut h: History = History::builder() + .key_type::() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .score_sigma(2.0) + .gamma(0.05) + .convergence(ConvergenceOptions { + max_iter: trueskill_tt::ITERATIONS, + epsilon: 1e-8, + alpha: 1.0, + }) + .build(); + + let mut k = 0usize; + for t in 0..slices { + for _ in 0..duels { + k += 1; + h.event(t) + .team([format!("p{}", k % competitors)]) + .team([format!("p{}", (k + 37) % competitors)]) + .scores([ + (k as f64 * 0.3).sin().abs() * 20.0, + (k as f64 * 0.3).cos().abs() * 20.0, + ]) + .commit() + .expect("ingests"); + } + } + h.converge().expect("converges"); + h +} + +/// Symbolic Cholesky by row-merge: returns (nnz(L), flops). +/// +/// Fill-in is simulated directly — for each column, the set of rows below the +/// diagonal that are nonzero — which is exact and easily checked, at the cost +/// of being O(n * nnz(L)) rather than the linear elimination-tree method. +fn symbolic(n: usize, adj: &[HashSet], perm_of: &[usize]) -> (usize, f64) { + // `perm_of[old] = new`. Build the permuted lower-triangle pattern. + let mut cols: Vec> = vec![HashSet::new(); n]; + for (old, nbrs) in adj.iter().enumerate() { + let i = perm_of[old]; + for &old_j in nbrs { + let j = perm_of[old_j]; + if j < i { + cols[j].insert(i); + } + } + } + + let mut nnz = 0usize; + let mut flops = 0.0f64; + for j in 0..n { + // Column j's pattern is final once every earlier column has merged in. + let rows: Vec = cols[j].iter().copied().collect(); + let c = rows.len(); + nnz += c + 1; // below-diagonal entries plus the diagonal + // Cholesky work for this column: one outer product over its pattern. + flops += (c as f64 + 1.0) * (c as f64 + 1.0); + // Fill-in: every pair in column j becomes an edge in the remaining graph. + for (a_idx, &a) in rows.iter().enumerate() { + for &b in &rows[a_idx + 1..] { + let (lo, hi) = if a < b { (a, b) } else { (b, a) }; + cols[lo].insert(hi); + } + } + } + (nnz, flops) +} + +#[test] +#[ignore = "measurement, run explicitly"] +fn what_sparsity_would_buy() { + for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200)] { + let h = fitted(slices, duels, competitors); + let (n, pattern) = h.joint_pattern_for_measurement(); + + let nnz_a: usize = pattern.iter().map(HashSet::len).sum::() + n; + let dense_flops = (n as f64).powi(3) / 3.0; + + let natural: Vec = (0..n).collect(); + let (nnz_nat, flops_nat) = symbolic(n, &pattern, &natural); + + // AMD returns `perm[new] = old`; invert it. + let (col_ptr, row_idx) = csc(n, &pattern); + let p = feral_amd::amd_order( + &feral_amd::CscPattern::new(n, &col_ptr, &row_idx).expect("valid pattern"), + ) + .expect("amd"); + let mut perm_of = vec![0usize; n]; + for (new, &old) in p.iter().enumerate() { + perm_of[old as usize] = new; + } + let (nnz_amd, flops_amd) = symbolic(n, &pattern, &perm_of); + + println!( + "\n=== {slices} slices x {duels} duels, {competitors} competitors ===\n\ + n = {n}\n\ + nnz(A) = {nnz_a} ({:.4}% dense)\n\ + dense flops = {:.3e}\n\ + nnz(L) natural = {nnz_nat} flops = {:.3e} ({:.1}x vs dense)\n\ + nnz(L) AMD = {nnz_amd} flops = {:.3e} ({:.1}x vs dense)", + 100.0 * nnz_a as f64 / (n * n) as f64, + dense_flops, + flops_nat, + dense_flops / flops_nat, + flops_amd, + dense_flops / flops_amd, + ); + } +} + +/// Full symmetric pattern to CSC, as `feral-amd` wants it. +fn csc(n: usize, adj: &[HashSet]) -> (Vec, Vec) { + let mut col_ptr = Vec::with_capacity(n + 1); + let mut row_idx = Vec::new(); + col_ptr.push(0i32); + for (j, nbrs) in adj.iter().enumerate() { + let mut rows: Vec = nbrs.iter().map(|&i| i as i32).collect(); + rows.push(j as i32); + rows.sort_unstable(); + rows.dedup(); + row_idx.extend_from_slice(&rows); + col_ptr.push(row_idx.len() as i32); + } + (col_ptr, row_idx) +} + +/// End-to-end factorisation time at the scale #52 was opened about. +#[test] +#[ignore = "measurement, run explicitly"] +fn factorisation_time_at_scale() { + use std::time::Instant; + + for (slices, duels, competitors) in [(30, 8, 100), (76, 13, 200), (150, 26, 400)] { + let h = fitted(slices, duels, competitors); + let (n, _) = h.joint_pattern_for_measurement(); + + // Warm, then time. + let _ = h.joint().expect("scored history"); + let t = Instant::now(); + let joint = h.joint().expect("scored history"); + let factor = t.elapsed(); + + let a = "p0".to_string(); + let b = "p1".to_string(); + let t = Instant::now(); + let _ = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).expect("known"); + let query = t.elapsed(); + + println!( + "n = {n:5} factorise = {factor:>12?} query = {query:>10?} \ + (dense was O(n^3): {:.3e} flops)", + (n as f64).powi(3) / 3.0 + ); + } +}