diff --git a/server/Cargo.toml b/server/Cargo.toml index 7de0c944b4a8c7e3123b818432fce82f09e6aeca..cdc2a78be0ffbd6ac27a02a332093daf5eb8b990 100644 --- a/server/Cargo.toml +++ b/server/Cargo.toml @@ -49,4 +49,8 @@ harness = false name = "ingest_replay" harness = false +[[bench]] +name = "stationary_solvers" +harness = false + diff --git a/server/benches/stationary_solvers.rs b/server/benches/stationary_solvers.rs new file mode 100644 index 0000000000000000000000000000000000000000..ff442b86e61861633ce2be1a016eae53aec73e6e --- /dev/null +++ b/server/benches/stationary_solvers.rs @@ -0,0 +1,171 @@ +//! Head-to-head cost of the stationary-distribution solvers behind Rank Centrality. +//! +//! Deliberately small: a couple of representative sizes per topology, ten +//! samples each, ~1 s of measurement per point. The point is the *ordering* +//! between solvers, which is stable at this resolution; for a finer sweep plus +//! accuracy numbers use `cargo run --release -p slugsocial-server --example +//! solver_probe`, which finishes in seconds. +//! +//! What each topology is here to show: +//! +//! - **chain** — spectral gap `Θ(1/n²)`. Power iteration needs `Θ(n²)` sweeps +//! and blows the 10 000 cap past n ≈ 1500; sparse GTH is `O(n)` because a +//! tree eliminates with zero fill. +//! - **star** — uniformization by `d_max = n-1` makes the chain almost purely +//! lazy, so power iteration crawls (it fails the cap at n = 1024) while +//! Gauss–Seidel finishes in two sweeps. +//! - **clique / sparse** — well-conditioned. Iterative wins outright and the +//! direct methods are the ones paying. + +use std::hint::black_box; +use std::time::Duration; + +use criterion::{criterion_group, criterion_main, BenchmarkId, Criterion}; + +use slugsocial_server::ranking::chain_from_edges; +use slugsocial_server::stationary::{ + bicgstab, dense_gth, dense_lu, power, power_aitken, solve, sor, sparse_gth, RankChain, + SolveOptions, SparseGthOutcome, +}; + +mod common; +use common::{build_group, Topology}; + +const TOL: f64 = 1e-8; +const CAP: usize = 10_000; + +fn opts() -> SolveOptions { + SolveOptions { + tol: TOL, + max_iters: CAP, + ..SolveOptions::default() + } +} + +fn chain_of(topo: Topology, n: usize) -> RankChain { + let g = build_group(topo, n); + chain_from_edges(g.idx_to_item.len(), g.edges.iter().map(|(&k, &w)| (k, w))) +} + +fn quick<'a>( + c: &'a mut Criterion, + name: &str, +) -> criterion::BenchmarkGroup<'a, criterion::measurement::WallTime> { + let mut g = c.benchmark_group(name); + g.sample_size(10) + .warm_up_time(Duration::from_millis(300)) + .measurement_time(Duration::from_secs(1)); + g +} + +/// Every solver on the same input, per topology. +fn solver_shootout(c: &mut Criterion) { + let cases: Vec<(Topology, usize, bool)> = vec![ + // (topology, n, include the O(n³) dense arms) + (Topology::Chain, 256, true), + (Topology::Chain, 1024, true), + (Topology::Star, 1024, true), + (Topology::Clique, 256, true), + (Topology::RandomSparse { degree: 6 }, 1024, false), + ]; + + for (topo, n, dense) in cases { + let chain = chain_of(topo, n); + let mut g = quick(c, &format!("solvers/{}", topo.label())); + + g.bench_with_input(BenchmarkId::new("power", n), &chain, |b, ch| { + b.iter(|| black_box(power(ch, opts()).residual)) + }); + g.bench_with_input(BenchmarkId::new("power+aitken", n), &chain, |b, ch| { + b.iter(|| black_box(power_aitken(ch, opts()).residual)) + }); + g.bench_with_input(BenchmarkId::new("sor", n), &chain, |b, ch| { + b.iter(|| black_box(sor(ch, opts(), 1.0).residual)) + }); + g.bench_with_input(BenchmarkId::new("bicgstab", n), &chain, |b, ch| { + b.iter(|| black_box(bicgstab(ch, opts()).residual)) + }); + g.bench_with_input(BenchmarkId::new("sparse-gth", n), &chain, |b, ch| { + let unlimited = SolveOptions { + direct_work_budget: u64::MAX, + dense_core_max: usize::MAX, + ..opts() + }; + b.iter(|| match sparse_gth(ch, unlimited, TOL) { + SparseGthOutcome::Solved(s) => black_box(s.residual), + SparseGthOutcome::TooDense { .. } => unreachable!(), + }) + }); + if dense { + g.bench_with_input(BenchmarkId::new("dense-gth", n), &chain, |b, ch| { + b.iter(|| black_box(dense_gth(ch, TOL).residual)) + }); + g.bench_with_input(BenchmarkId::new("dense-lu", n), &chain, |b, ch| { + b.iter(|| black_box(dense_lu(ch, TOL).residual)) + }); + } + g.finish(); + } +} + +/// The shipped path, including the cost of deciding which solver to use. +fn hybrid_dispatch(c: &mut Criterion) { + let cases: Vec<(String, Topology, usize)> = vec![ + ("chain".into(), Topology::Chain, 4000), + ("star".into(), Topology::Star, 4000), + ("clique".into(), Topology::Clique, 400), + ( + "sparse-d6".into(), + Topology::RandomSparse { degree: 6 }, + 4000, + ), + ( + "components".into(), + Topology::ManyComponents { + components: 64, + size: 64, + }, + 0, + ), + ]; + let mut g = quick(c, "hybrid"); + for (label, topo, n) in cases { + let chain = chain_of(topo, n); + g.bench_with_input(BenchmarkId::new(label, chain.n), &chain, |b, ch| { + b.iter(|| black_box(solve(ch, opts()).residual)) + }); + } + g.finish(); +} + +/// Building the chain (aggregate, pairwise-normalize, sort) versus solving it. +/// Sorting is what buys run-to-run determinism; this is where to see its price. +fn chain_construction(c: &mut Criterion) { + let mut g = quick(c, "chain_from_edges"); + for (topo, n) in [ + (Topology::RandomSparse { degree: 6 }, 4000usize), + (Topology::Clique, 400), + ] { + let group = build_group(topo, n); + let total = group.idx_to_item.len(); + g.bench_with_input( + BenchmarkId::new(topo.label(), group.edges.len()), + &group, + |b, grp| { + b.iter(|| { + let ch = chain_from_edges(total, grp.edges.iter().map(|(&k, &w)| (k, w))); + black_box(ch.nnz()) + }) + }, + ); + } + g.finish(); +} + +criterion_group!( + benches, + solver_shootout, + hybrid_dispatch, + chain_construction +); +criterion_main!(benches); diff --git a/server/examples/solver_probe.rs b/server/examples/solver_probe.rs new file mode 100644 index 0000000000000000000000000000000000000000..e3933530ac9efab379dc10ec2a467eb30ad203f8 --- /dev/null +++ b/server/examples/solver_probe.rs @@ -0,0 +1,738 @@ +//! Head-to-head comparison of stationary-distribution solvers for Rank Centrality. +//! +//! Run: `cargo run --release -p slugsocial-server --example solver_probe` +//! +//! Reports, per (topology, solver): wall time, backward error `‖πP − π‖₁`, and +//! how badly the produced *order* disagrees with ground truth. On a chain the +//! ground truth order is known exactly by construction, and on any tree the +//! exact scores are known in closed form from detailed balance, so "how wrong +//! is the current implementation" is answerable without trusting any solver. + +use std::time::Instant; + +use slugsocial_server::path_types::ItemId; +use slugsocial_server::ranking::chain_from_edges; +use slugsocial_server::reducer::{GroupState, VoteData}; +use slugsocial_server::stationary::{ + bicgstab, dense_gth, dense_lu, power, power_aitken, solve, sor, sparse_gth, Method, RankChain, + Solution, SolveOptions, SparseGthOutcome, +}; + +const TOL: f64 = 1e-8; +const CAP: usize = 10_000; + +// --------------------------------------------------------------------------- +// Graph generators (mirrors server/benches/common/mod.rs so numbers line up) +// --------------------------------------------------------------------------- + +struct Rng(u64); +impl Rng { + fn new(s: u64) -> Self { + Rng(s | 1) + } + fn next_u64(&mut self) -> u64 { + let mut x = self.0; + x ^= x >> 12; + x ^= x << 25; + x ^= x >> 27; + self.0 = x; + x.wrapping_mul(0x2545_F491_4F6C_DD1D) + } + fn below(&mut self, n: usize) -> usize { + (self.next_u64() % n as u64) as usize + } +} + +fn item_name(i: usize) -> String { + format!("~/bench/i{i:06}") +} + +fn vote(a: usize, b: usize, l: i32, r: i32, ts: i64) -> VoteData { + VoteData { + ts, + a: ItemId::parse(&item_name(a)).unwrap(), + b: ItemId::parse(&item_name(b)).unwrap(), + ratio_left: l, + ratio_right: r, + body: "synthetic".into(), + principal: "bench".into(), + delegate: None, + thread_tag: "bench".into(), + } +} + +#[derive(Clone, Copy)] +enum Topo { + /// `0 > 1 > 2 > …`, every vote at a fixed ratio. Ground-truth order known. + Chain { + ratio: i32, + }, + /// Chain with per-edge varying ratios (what `benches/common` builds). + ChainVaried, + Star, + Clique, + Sparse { + degree: usize, + }, + Components { + count: usize, + size: usize, + }, +} + +fn pairs(t: Topo, n: usize) -> Vec<(usize, usize, i32)> { + let mut out = Vec::new(); + match t { + Topo::Chain { ratio } => { + for i in 0..n.saturating_sub(1) { + out.push((i, i + 1, ratio)); + } + } + Topo::ChainVaried => { + for i in 0..n.saturating_sub(1) { + out.push((i, i + 1, 2 + (i % 5) as i32)); + } + } + Topo::Star => { + for i in 1..n { + out.push((0, i, 2 + (i % 5) as i32)); + } + } + Topo::Clique => { + let mut k = 0; + for i in 0..n { + for j in (i + 1)..n { + out.push((i, j, 2 + (k % 5))); + k += 1; + } + } + } + Topo::Sparse { degree } => { + let mut rng = Rng::new(0xC0FFEE ^ n as u64); + for i in 0..n.saturating_sub(1) { + out.push((i, i + 1, 2)); + } + for k in 0..(degree.saturating_sub(2) * n / 2) { + let a = rng.below(n); + let b = rng.below(n); + if a != b { + out.push((a, b, 2 + (k % 5) as i32)); + } + } + } + Topo::Components { count, size } => { + for c in 0..count { + let base = c * size; + for i in 0..size.saturating_sub(1) { + out.push((base + i, base + i + 1, 2)); + } + } + } + } + out +} + +fn node_count(t: Topo, n: usize) -> usize { + match t { + Topo::Components { count, size } => count * size, + _ => n, + } +} + +fn build_chain(t: Topo, n: usize) -> RankChain { + let total = node_count(t, n); + let mut g = GroupState::new(); + for i in 0..total { + g.ensure_item_pub(&item_name(i)); + } + for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() { + g.apply_vote(vote(i, j, ratio, 1, k as i64)); + } + chain_from_edges(total, g.edges.iter().map(|(&k, &w)| (k, w))) +} + +// --------------------------------------------------------------------------- +// Accuracy metrics +// --------------------------------------------------------------------------- + +/// Kendall-tau distance between the order induced by `a` and by `b`, as a +/// fraction of all pairs. 0 = identical ranking, 0.5 = as good as random. +fn kendall_distance(a: &[f64], b: &[f64]) -> f64 { + let n = a.len(); + if n < 2 { + return 0.0; + } + // Sort indices by `a` descending, then count inversions with respect to `b`. + let mut idx: Vec = (0..n).collect(); + idx.sort_by(|&x, &y| b[y].partial_cmp(&b[x]).unwrap_or(std::cmp::Ordering::Equal)); + let mut seq: Vec = idx.iter().map(|&i| a[i]).collect(); + let mut buf = seq.clone(); + let inv = count_inversions(&mut seq, &mut buf); + let total = (n as f64) * (n as f64 - 1.0) / 2.0; + inv as f64 / total +} + +/// Inversions = pairs out of descending order (strict; ties are concordant). +fn count_inversions(v: &mut [f64], buf: &mut [f64]) -> u64 { + let n = v.len(); + if n < 2 { + return 0; + } + let mid = n / 2; + let (l, r) = v.split_at_mut(mid); + let (bl, br) = buf.split_at_mut(mid); + let mut inv = count_inversions(l, bl) + count_inversions(r, br); + let (mut i, mut j, mut k) = (0usize, 0usize, 0usize); + while i < l.len() && j < r.len() { + if l[i] >= r[j] { + buf[k] = l[i]; + i += 1; + } else { + // every remaining element of `l` is smaller than r[j]: inversion + inv += (l.len() - i) as u64; + buf[k] = r[j]; + j += 1; + } + k += 1; + } + while i < l.len() { + buf[k] = l[i]; + i += 1; + k += 1; + } + while j < r.len() { + buf[k] = r[j]; + j += 1; + k += 1; + } + v.copy_from_slice(&buf[..n]); + inv +} + +/// The pre-change chain builder, kept verbatim so determinism can be measured +/// rather than argued about: adjacency rows and their weight sums are produced +/// by iterating a `HashMap`, whose order varies per instance. +fn legacy_chain(n: usize, edges: impl Iterator) -> RankChain { + use std::collections::{HashMap, HashSet}; + let mut raw: HashMap<(usize, usize), f64> = HashMap::new(); + for ((src, dst), w) in edges { + if src >= n || dst >= n || w <= 0.0 { + continue; + } + *raw.entry((src, dst)).or_insert(0.0) += w; + } + let keys: Vec<(usize, usize)> = raw.keys().copied().collect(); + let mut normalized: HashMap<(usize, usize), f64> = HashMap::new(); + for (i, j) in keys { + if normalized.contains_key(&(i, j)) { + continue; + } + let w_ij = *raw.get(&(i, j)).unwrap_or(&0.0); + let w_ji = *raw.get(&(j, i)).unwrap_or(&0.0); + let total = w_ij + w_ji; + if total <= 0.0 { + continue; + } + normalized.insert((i, j), w_ij / total); + if w_ji > 0.0 { + normalized.insert((j, i), w_ji / total); + } + } + let mut rows: Vec> = vec![Vec::new(); n]; + let mut neighbors: Vec> = vec![HashSet::new(); n]; + for ((src, dst), w) in &normalized { + rows[*src].push((*dst, *w)); + neighbors[*src].insert(*dst); + neighbors[*dst].insert(*src); + } + let row_sum: Vec = rows + .iter() + .map(|r| r.iter().map(|&(_, w)| w).sum()) + .collect(); + let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0) as f64; + RankChain { + n, + rows, + row_sum, + d_max, + } +} + +/// Pairs whose relative order differs between two score vectors. +fn order_flips(a: &[f64], b: &[f64]) -> usize { + let n = a.len(); + let mut idx: Vec = (0..n).collect(); + idx.sort_by(|&x, &y| b[y].partial_cmp(&b[x]).unwrap_or(std::cmp::Ordering::Equal)); + let mut seq: Vec = idx.iter().map(|&i| a[i]).collect(); + let mut buf = seq.clone(); + count_inversions(&mut seq, &mut buf) as usize +} + +/// Adjacent positions where a chain's known-correct descending order breaks, +/// split into strict inversions (genuinely the wrong way round) and ties (the +/// two items are indistinguishable, so their displayed order is arbitrary). +fn chain_order_violations(pi: &[f64]) -> (usize, usize) { + let mut inv = 0; + let mut tie = 0; + for i in 0..pi.len().saturating_sub(1) { + if pi[i] < pi[i + 1] { + inv += 1; + } else if pi[i] == pi[i + 1] { + tie += 1; + } + } + (inv, tie) +} + +/// Exact π for a chain via detailed balance, in log space (no overflow). +/// +/// A chain is a tree, so the Markov chain is reversible and +/// `π_i / π_{i+1} = a_{i+1,i} / a_{i,i+1}` holds exactly edge by edge. +fn chain_exact_log(t: Topo, n: usize) -> Vec { + let mut logs = vec![0.0f64; n]; + for (i, (_, _, ratio)) in pairs(t, n).into_iter().enumerate() { + // node i preferred over node i+1 at `ratio`:1, so π_i / π_{i+1} = ratio + logs[i + 1] = logs[i] - (ratio as f64).ln(); + } + logs +} + +/// Largest absolute error in `ln π_i` against an exact log-space reference, +/// restricted to entries the reference says are representable in f64. +fn max_log_error(pi: &[f64], exact_log: &[f64]) -> (f64, usize) { + let ref_max = exact_log.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + let pi_max = pi.iter().cloned().fold(0.0f64, f64::max); + if pi_max <= 0.0 { + return (f64::INFINITY, pi.len()); + } + let mut worst = 0.0f64; + let mut lost = 0usize; + for i in 0..pi.len() { + let want = exact_log[i] - ref_max; + if want < -700.0 { + continue; // genuinely below f64 range; not the solver's fault + } + if pi[i] <= 0.0 { + lost += 1; + continue; + } + let got = (pi[i] / pi_max).ln(); + worst = worst.max((got - want).abs()); + } + (worst, lost) +} + +// --------------------------------------------------------------------------- +// Runner +// --------------------------------------------------------------------------- + +struct Run { + name: &'static str, + ms: f64, + sol: Solution, +} + +fn time_it(name: &'static str, mut f: impl FnMut() -> Solution) -> Run { + // One warm-up, then the best of three: solver cost here is deterministic, + // so the minimum is the cleanest estimate of the real cost. + let sol = f(); + let mut best = f64::INFINITY; + for _ in 0..3 { + let t = Instant::now(); + let _ = f(); + best = best.min(t.elapsed().as_secs_f64() * 1000.0); + } + Run { + name, + ms: best, + sol, + } +} + +fn run_all(chain: &RankChain, include_dense: bool) -> Vec { + let opts = SolveOptions { + tol: TOL, + max_iters: CAP, + ..SolveOptions::default() + }; + let big_opts = SolveOptions { + tol: TOL, + max_iters: 20_000_000, + ..SolveOptions::default() + }; + + let mut runs = vec![ + time_it("power(10k) [current]", || power(chain, opts)), + time_it("power+aitken", || power_aitken(chain, opts)), + time_it("sor(1.0)", || sor(chain, opts, 1.0)), + time_it("bicgstab", || bicgstab(chain, opts)), + time_it("solve() [hybrid]", || solve(chain, opts)), + ]; + if include_dense { + runs.push(time_it("sparse-gth(uncapped)", || force_sparse_gth(chain))); + runs.push(time_it("dense-gth", || dense_gth(chain, TOL))); + runs.push(time_it("dense-lu", || dense_lu(chain, TOL))); + } + // Power iteration with an effectively unlimited cap: what the current + // algorithm would produce if it were allowed to finish. + if chain.n <= 1200 && chain.nnz() <= 8000 { + runs.push(time_it("power(unbounded)", || power(chain, big_opts))); + } + runs +} + +/// Sparse GTH with the budget disabled, for measuring what the direct path +/// would cost even on graphs the hybrid would refuse. +fn force_sparse_gth(chain: &RankChain) -> Solution { + let unlimited = SolveOptions { + direct_work_budget: u64::MAX, + dense_core_max: usize::MAX, + ..SolveOptions::default() + }; + match sparse_gth(chain, unlimited, TOL) { + SparseGthOutcome::Solved(s) => s, + SparseGthOutcome::TooDense { .. } => unreachable!(), + } +} + +fn print_block(label: &str, chain: &RankChain, reference: &[f64], exact_log: Option<&[f64]>) { + println!( + "\n--- {label} n={} nnz={} d_max={}", + chain.n, + chain.nnz(), + chain.d_max as usize + ); + println!( + "{:<22} {:>9} {:>11} {:>5} {:>7} {:>10} {:>9} {:>10} {:>10} {:>10}", + "method", + "time(ms)", + "resid L1", + "conv", + "iters", + "kendall-d", + "log-err", + "inv/displ", + "tie/maxsh", + "kd(log pi)" + ); + for r in run_all(chain, chain.n <= 1500) { + let kd = kendall_distance(&r.sol.pi, reference); + let kd_log = kendall_distance(&r.sol.log_pi, reference); + let (logerr, _) = match exact_log { + Some(e) => max_log_error(&r.sol.pi, e), + None => (f64::NAN, 0usize), + }; + let (inv, tie) = chain_order_violations(&r.sol.pi); + let (inv, tie) = if exact_log.is_some() { + // How the ranking would actually be displayed: sort by score + // descending (ties broken by index, as a stable sort does) and see + // how many items land somewhere other than their true rank. + let n = r.sol.pi.len(); + let mut order: Vec = (0..n).collect(); + order.sort_by(|&a, &b| { + r.sol.pi[b] + .partial_cmp(&r.sol.pi[a]) + .unwrap_or(std::cmp::Ordering::Equal) + }); + let misplaced = order.iter().enumerate().filter(|&(p, &i)| p != i).count(); + let shift = order + .iter() + .enumerate() + .map(|(p, &i)| p.abs_diff(i)) + .max() + .unwrap_or(0); + (format!("{inv}/{misplaced}"), format!("{tie}/{shift}")) + } else { + ("-".to_string(), "-".to_string()) + }; + println!( + "{:<22} {:>9.3} {:>11.2e} {:>5} {:>7} {:>10.2e} {:>9} {:>10} {:>10} {:>10.2e}", + r.name, + r.ms, + r.sol.residual, + if r.sol.converged { "yes" } else { "NO" }, + r.sol.iterations, + kd, + if logerr.is_nan() { + "-".to_string() + } else { + format!("{logerr:.2e}") + }, + inv, + tie, + kd_log + ); + } +} + +fn main() { + let want: Vec = std::env::args().skip(1).collect(); + let on = |s: &str| want.is_empty() || want.iter().any(|w| w == s); + println!("Rank Centrality stationary solvers — tol={TOL:e}, power cap={CAP}"); + + // ---- Chains: ground truth known exactly, and where the cap binds. ---- + if on("chain") { + println!("\n================ CHAIN (ratio 2:1, exact order known) ================"); + for n in [100usize, 256, 512, 1024, 2048, 4000] { + let t = Topo::Chain { ratio: 2 }; + let chain = build_chain(t, n); + let exact = chain_exact_log(t, n); + // Reference order is the construction order: node 0 best. + let reference: Vec = (0..n).map(|i| -(i as f64)).collect(); + print_block(&format!("chain n={n}"), &chain, &reference, Some(&exact)); + println!( + " ground truth: π_0/π_{} = 2^{} = 10^{:.0}", + n - 1, + n - 1, + (n - 1) as f64 * 2f64.log10() + ); + } + } + if on("varied") { + println!( + "\n================ CHAIN (varied ratios 2..6, benches/common shape) ================" + ); + for n in [256usize, 1024] { + let t = Topo::ChainVaried; + let chain = build_chain(t, n); + let exact = chain_exact_log(t, n); + let reference: Vec = (0..n).map(|i| -(i as f64)).collect(); + print_block( + &format!("chain-varied n={n}"), + &chain, + &reference, + Some(&exact), + ); + } + + // ---- Other topologies: reference is sparse GTH (subtraction-free). ---- + } + if on("topo") { + println!("\n================ OTHER TOPOLOGIES ================"); + let cases: Vec<(String, Topo, usize)> = vec![ + ("star".into(), Topo::Star, 1024), + ("clique".into(), Topo::Clique, 128), + ("clique".into(), Topo::Clique, 400), + ("sparse-d6".into(), Topo::Sparse { degree: 6 }, 1024), + ("sparse-d6".into(), Topo::Sparse { degree: 6 }, 4000), + ("sparse-d20".into(), Topo::Sparse { degree: 20 }, 1024), + ( + "components 64x64".into(), + Topo::Components { + count: 64, + size: 64, + }, + 0, + ), + ]; + for (label, t, n) in cases { + let chain = build_chain(t, n); + let reference = solve( + &chain, + SolveOptions { + tol: TOL, + max_iters: CAP, + ..SolveOptions::default() + }, + ) + .pi; + print_block( + &format!("{label} n={}", node_count(t, n)), + &chain, + &reference, + None, + ); + } + + // ---- Which method does the hybrid pick, and what does elimination cost? ---- + } + if on("dispatch") { + println!("\n================ HYBRID DISPATCH + FILL-IN ================"); + println!( + "{:<24} {:>7} {:>9} {:>14} {:>10} {:>12}", + "graph", "n", "nnz", "picked", "iters", "resid" + ); + let dispatch: Vec<(String, Topo, usize)> = vec![ + ("chain".into(), Topo::Chain { ratio: 2 }, 4000), + ("star".into(), Topo::Star, 4000), + ("clique".into(), Topo::Clique, 200), + ("clique".into(), Topo::Clique, 600), + ("clique".into(), Topo::Clique, 1200), + ("sparse-d6".into(), Topo::Sparse { degree: 6 }, 4000), + ("sparse-d20".into(), Topo::Sparse { degree: 20 }, 4000), + ("sparse-d6".into(), Topo::Sparse { degree: 6 }, 20000), + ]; + for (label, t, n) in dispatch { + let chain = build_chain(t, n); + let opts = SolveOptions { + tol: TOL, + max_iters: CAP, + ..SolveOptions::default() + }; + let t0 = Instant::now(); + let s = solve(&chain, opts); + let ms = t0.elapsed().as_secs_f64() * 1000.0; + println!( + "{:<24} {:>7} {:>9} {:>14} {:>10} {:>12.2e} {:.2} ms {}", + label, + chain.n, + chain.nnz(), + s.method.label(), + s.iterations, + s.residual, + ms, + if s.converged { "" } else { "NOT CONVERGED" } + ); + } + + // ---- Determinism: does HashMap iteration order move the scores? ---- + } + if on("determinism") { + println!("\n================ DETERMINISM ================"); + for (label, t, n) in [ + ("sparse-d6 n=1024", Topo::Sparse { degree: 6 }, 1024usize), + ("clique n=128", Topo::Clique, 128), + ] { + let mut first: Option> = None; + let mut worst = 0.0f64; + for _ in 0..8 { + let chain = build_chain(t, n); + let pi = force_sparse_gth(&chain).pi; + match &first { + None => first = Some(pi), + Some(f) => { + for i in 0..pi.len() { + worst = worst.max((pi[i] - f[i]).abs()); + } + } + } + } + println!("{label:<24} max |Δπ| across 8 rebuilds: {worst:.3e}"); + } + + // Control: the pre-change builder, rebuilt from scratch each round so + // every rebuild gets a fresh `HashMap` with a different iteration order. + println!("\ncontrol — legacy HashMap-ordered builder, 8 rebuilds of the same graph:"); + for (label, t, n) in [ + ("sparse-d6 n=1024", Topo::Sparse { degree: 6 }, 1024usize), + ("clique n=128", Topo::Clique, 128), + ("chain n=512", Topo::Chain { ratio: 2 }, 512), + ] { + let total = node_count(t, n); + let mut g = GroupState::new(); + for i in 0..total { + g.ensure_item_pub(&item_name(i)); + } + for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() { + g.apply_vote(vote(i, j, ratio, 1, k as i64)); + } + let opts = SolveOptions { + tol: TOL, + max_iters: CAP, + ..SolveOptions::default() + }; + let mut reference: Option> = None; + let mut worst_abs = 0.0f64; + let mut worst_rel = 0.0f64; + let mut flips = 0usize; + let mut distinct = 0usize; + for _ in 0..8 { + let c = legacy_chain(total, g.edges.iter().map(|(&k, &w)| (k, w))); + let pi = power(&c, opts).pi; + match &reference { + None => reference = Some(pi), + Some(r) => { + if pi != *r { + distinct += 1; + } + for i in 0..pi.len() { + worst_abs = worst_abs.max((pi[i] - r[i]).abs()); + if r[i] > 0.0 { + worst_rel = worst_rel.max((pi[i] - r[i]).abs() / r[i]); + } + } + flips += order_flips(&pi, r); + } + } + } + println!( + "{label:<20} differing rebuilds {distinct}/7 max |Δπ|={worst_abs:.2e} \ + max rel={worst_rel:.2e} rank flips {flips}" + ); + } + + // Cost of the sort that buys determinism, isolated from the `HashMap` + // aggregation `chain_from_edges` has always done. + println!("\nsort cost (edge ordering inside chain_from_edges):"); + for (label, t, n) in [ + ("sparse-d6 n=4000", Topo::Sparse { degree: 6 }, 4000usize), + ("clique n=400", Topo::Clique, 400), + ] { + let total = node_count(t, n); + let mut g = GroupState::new(); + for i in 0..total { + g.ensure_item_pub(&item_name(i)); + } + for (k, (i, j, ratio)) in pairs(t, total).into_iter().enumerate() { + g.apply_vote(vote(i, j, ratio, 1, k as i64)); + } + let t0 = Instant::now(); + for _ in 0..20 { + let c = chain_from_edges(total, g.edges.iter().map(|(&k, &w)| (k, w))); + std::hint::black_box(c.nnz()); + } + let with_sort = t0.elapsed().as_secs_f64() * 1000.0 / 20.0; + + let edges: Vec<((usize, usize), f64)> = g.edges.iter().map(|(&k, &w)| (k, w)).collect(); + let t0 = Instant::now(); + for _ in 0..20 { + let mut keys: Vec<(usize, usize)> = edges.iter().map(|&(k, _)| k).collect(); + keys.sort_unstable(); + std::hint::black_box(keys.len()); + } + let sort_only = t0.elapsed().as_secs_f64() * 1000.0 / 20.0; + println!( + "{label:<20} chain_from_edges {with_sort:7.3} ms, of which sorting is \ + {sort_only:6.3} ms (E={})", + g.edges.len() + ); + } + + // ---- Cross-method equivalence on small graphs, at machine precision. ---- + } + if on("equiv") { + println!("\n================ EQUIVALENCE (small n, all methods) ================"); + for (label, t, n) in [ + ("chain n=64", Topo::Chain { ratio: 2 }, 64usize), + ("star n=64", Topo::Star, 64), + ("clique n=64", Topo::Clique, 64), + ("sparse-d6 n=200", Topo::Sparse { degree: 6 }, 200), + ] { + let chain = build_chain(t, n); + let gth = dense_gth(&chain, TOL).pi; + let mut worst: Vec<(String, f64)> = Vec::new(); + let opts = SolveOptions { + tol: 1e-14, + max_iters: 200_000, + ..SolveOptions::default() + }; + let cands: Vec<(&str, Vec)> = vec![ + ("sparse-gth", force_sparse_gth(&chain).pi), + ("dense-lu", dense_lu(&chain, TOL).pi), + ("power(200k)", power(&chain, opts).pi), + ("sor", sor(&chain, opts, 1.0).pi), + ("bicgstab", bicgstab(&chain, opts).pi), + ]; + for (name, pi) in cands { + let rel = (0..chain.n) + .filter(|&i| gth[i] > 1e-300) + .map(|i| (pi[i] - gth[i]).abs() / gth[i]) + .fold(0.0f64, f64::max); + worst.push((name.to_string(), rel)); + } + let cells: Vec = worst.iter().map(|(n, v)| format!("{n}={v:.1e}")).collect(); + println!("{label:<18} max rel err vs dense-GTH: {}", cells.join(" ")); + } + } + let _ = Method::Power; +} diff --git a/server/src/lib.rs b/server/src/lib.rs index e102e73ec167fa35dcc22d4ec8c6f40202982429..7845de0b8512738fe20ca9dffd84e9b8c8134d50 100644 --- a/server/src/lib.rs +++ b/server/src/lib.rs @@ -16,6 +16,7 @@ pub mod ranking; pub mod reducer; pub mod scope_rank; pub mod state; +pub mod stationary; pub mod timeago; pub mod views; pub mod write_cmd; diff --git a/server/src/ranking.rs b/server/src/ranking.rs index 38e6d09b4370e5f8cbae09c0e5760b4e7f1ef7db..2d391bb6fdf4d22247e5524d459f2d1609bd6d70 100644 --- a/server/src/ranking.rs +++ b/server/src/ranking.rs @@ -2,6 +2,7 @@ use std::collections::{HashMap, HashSet}; use crate::path_types::ItemId; use crate::reducer::GroupState; +use crate::stationary::{self, RankChain, SolveOptions, Solution}; #[derive(Debug, Clone)] pub struct RankedItem { @@ -105,14 +106,21 @@ pub fn ranked_items(group: &mut GroupState, max_iters: usize, tol: f64) -> Vec, max_iters: usize, tol: f64) -> Vec { - if n == 0 { - return vec![]; - } - if n == 1 { - return vec![1.0]; - } - +/// Build the Rank Centrality Markov chain from raw directed vote weights. +/// +/// Rank Centrality (Negahban, Oh, Shah 2012, §3.1): +/// P_ij = (1/d_max) * a_ij for i ≠ j compared +/// P_ii = 1 - (1/d_max) * Σ_k a_ik +/// where d_i is the *degree* (number of distinct neighbors compared) and +/// d_max = max_i d_i. Using the unweighted degree — not the sum of +/// pairwise-normalized weights — is what guarantees aperiodicity: it +/// forces P_ii > 0 for every non-maximum-degree node, and for max-degree +/// nodes whenever any neighbor weight is below 1 (i.e. not a unanimous +/// loss). Without this, regular comparison graphs (e.g. a pure star at +/// ratio 2:1) produce a bipartite chain that oscillates instead of +/// converging — see issue #146. (Direct solvers are immune either way: +/// π depends only on the off-diagonals, which d_max scales uniformly.) +pub fn chain_from_edges(n: usize, edges: impl Iterator) -> RankChain { // Collect raw edges into a map for pairwise normalization. let mut raw: HashMap<(usize, usize), f64> = HashMap::new(); for ((src, dst), w) in edges { @@ -125,90 +133,80 @@ fn compute_scores_from_edges(n: usize, edges: impl Iterator = raw.keys().copied().collect(); - let mut normalized: HashMap<(usize, usize), f64> = HashMap::new(); + // + // Sorted, not read straight off the `HashMap`: the summation order of every + // downstream reduction has to be a function of the graph alone, or scores + // wobble in their low bits between runs and near-ties can flip. + let mut keys: Vec<(usize, usize)> = raw.keys().copied().collect(); + keys.sort_unstable(); + + let mut normalized: Vec<((usize, usize), f64)> = Vec::with_capacity(keys.len()); + let mut neighbors: Vec> = vec![HashSet::new(); n]; for (i, j) in keys { - if normalized.contains_key(&(i, j)) { - continue; - } let w_ij = *raw.get(&(i, j)).unwrap_or(&0.0); let w_ji = *raw.get(&(j, i)).unwrap_or(&0.0); let total = w_ij + w_ji; if total <= 0.0 { continue; } - normalized.insert((i, j), w_ij / total); - if w_ji > 0.0 { - normalized.insert((j, i), w_ji / total); - } + normalized.push(((i, j), w_ij / total)); + neighbors[i].insert(j); + neighbors[j].insert(i); } - // Rank Centrality (Negahban, Oh, Shah 2012, §3.1): - // P_ij = (1/d_max) * A_ij for i ≠ j compared - // P_ii = 1 - (1/d_max) * Σ_k A_ik - // where d_i is the *degree* (number of distinct neighbors compared) and - // d_max = max_i d_i. Using the unweighted degree — not the sum of - // pairwise-normalized weights — is what guarantees aperiodicity: it - // forces P_ii > 0 for every non-maximum-degree node, and for max-degree - // nodes whenever any neighbor weight is below 1 (i.e. not a unanimous - // loss). Without this, regular comparison graphs (e.g. a pure star at - // ratio 2:1) produce a bipartite chain that oscillates instead of - // converging — see issue #146. - let mut out_edges: Vec> = vec![Vec::new(); n]; - let mut neighbors: Vec> = vec![HashSet::new(); n]; + let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0); + RankChain::from_normalized(n, normalized, d_max) +} - for ((src, dst), w) in &normalized { - out_edges[*src].push((*dst, *w)); - neighbors[*src].insert(*dst); - neighbors[*dst].insert(*src); +/// Stationary distribution for the chain induced by `edges`, with convergence +/// diagnostics attached. See [`stationary::Solution`]. +/// +/// This is the API that replaces "return whatever the iteration reached": +/// callers that care can inspect `converged`, `residual` and `method`, and +/// non-convergence is logged rather than swallowed. +pub fn solve_scores_from_edges( + n: usize, + edges: impl Iterator, + max_iters: usize, + tol: f64, +) -> Solution { + if n == 0 { + return stationary::trivial(vec![]); } - - let weight_sum: Vec = out_edges - .iter() - .map(|es| es.iter().map(|(_, w)| *w).sum()) - .collect(); - let d_max = neighbors.iter().map(|s| s.len()).max().unwrap_or(0); - if d_max == 0 { - return vec![1.0 / n as f64; n]; + if n == 1 { + return stationary::trivial(vec![1.0]); } - let d_max_f = d_max as f64; - - let mut scores = vec![1.0 / n as f64; n]; - let mut next = vec![0.0f64; n]; - - for _ in 0..max_iters { - next.fill(0.0); - for i in 0..n { - let stay_prob = (d_max_f - weight_sum[i]) / d_max_f; - next[i] += scores[i] * stay_prob; - - if out_edges[i].is_empty() { - continue; - } - for &(dst, w) in &out_edges[i] { - next[dst] += scores[i] * (w / d_max_f); - } - } - let diff: f64 = scores - .iter() - .zip(next.iter()) - .map(|(a, b)| (a - b).abs()) - .sum(); + let chain = chain_from_edges(n, edges); + let solution = stationary::solve( + &chain, + SolveOptions { + tol, + max_iters, + ..SolveOptions::default() + }, + ); - scores.clone_from_slice(&next); - if diff < tol { - break; - } + if !solution.converged { + tracing::warn!( + n, + method = solution.method.label(), + iterations = solution.iterations, + residual = solution.residual, + "rank centrality did not reach tolerance; this ranking may be misordered" + ); + } else if solution.underflowed { + tracing::debug!( + n, + method = solution.method.label(), + "rank centrality scores span more than f64 holds; use log scores to order the tail" + ); } + solution +} - let sum: f64 = scores.iter().sum(); - if sum.is_finite() && sum > 0.0 { - for s in &mut scores { - *s /= sum; - } - } - scores +fn compute_scores_from_edges(n: usize, edges: impl Iterator, max_iters: usize, tol: f64) -> Vec { + solve_scores_from_edges(n, edges, max_iters, tol).pi } /// Rank-centrality within a subset of items (an induced subgraph), using the group's aggregated edges. @@ -231,20 +229,25 @@ pub fn ranked_items_subset(group: &GroupState, idxs: &[usize], max_iters: usize, Some(((s, d), w)) }); - let scores = compute_scores_from_edges(idxs.len(), edges_iter, max_iters, tol); + let solved = solve_scores_from_edges(idxs.len(), edges_iter, max_iters, tol); // Filter out entries where idx_to_item doesn't have the slot (shouldn't happen, but be safe). - let mut items: Vec = idxs + // Sorting on the *log* score, not `score`: on a long preference chain the + // true distribution spans more decades than f64 holds, so `score` ties off + // at the bottom while the log scores still order it correctly. + let mut items: Vec<(RankedItem, f64)> = idxs .iter() .enumerate() .filter_map(|(j, &orig)| { let item = group.idx_to_item.get(orig)?.clone(); - Some(RankedItem { item, score: *scores.get(j).unwrap_or(&0.0) }) + let score = *solved.pi.get(j).unwrap_or(&0.0); + let log_score = *solved.log_pi.get(j).unwrap_or(&f64::NEG_INFINITY); + Some((RankedItem { item, score }, log_score)) }) .collect(); - items.sort_by(|a, b| b.score.partial_cmp(&a.score).unwrap_or(std::cmp::Ordering::Equal)); - items + items.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)); + items.into_iter().map(|(item, _)| item).collect() } pub fn group_summary_scores( @@ -314,6 +317,82 @@ mod tests { ); } + /// A long preference chain — the natural shape of a deep ontology or a + /// hand-ordered list — is where power iteration got the answer wrong. With + /// varied ratios and 1200 items the old solver reported success (its L1 + /// step fell below 1e-8 after 4325 sweeps, well inside the cap) while + /// leaving 47% of items at the wrong rank, because an absolute step + /// tolerance says nothing about entries that are themselves 1e-300. + #[test] + fn deep_chain_is_ranked_in_exactly_the_right_order() { + let n = 1200usize; + let mut g = mk_group(); + let name = |i: usize| format!("i{i:05}"); + for i in 0..n - 1 { + let left = 2 + (i % 5) as i32; + g.apply_vote(vote(i as i64, &name(i), &name(i + 1), left, 1)); + } + + let idxs: Vec = (0..g.idx_to_item.len()).collect(); + let ranked = ranked_items_subset(&g, &idxs, 10000, 1e-8); + assert_eq!(ranked.len(), n); + for (i, r) in ranked.iter().enumerate() { + let want = format!("https://slug.social/{}", name(i)); + assert_eq!(r.item.as_str(), want, "position {i} of the chain ranking"); + } + } + + /// Same graph, edges presented in a different order: the scores must come + /// back bit-identical, not merely close. + #[test] + fn ranking_does_not_depend_on_edge_iteration_order() { + let n = 200usize; + let mut edges: Vec<((usize, usize), f64)> = Vec::new(); + for i in 0..n - 1 { + edges.push(((i + 1, i), 3.0)); + edges.push(((i, i + 1), 1.0)); + } + for i in 0..n / 4 { + let j = (5 * i + 7) % n; + if i != j { + edges.push(((j, i), 2.0)); + edges.push(((i, j), 1.0)); + } + } + + let reference = compute_scores_from_edges(n, edges.iter().copied(), 10000, 1e-8); + let mut state = 0x5EEDu64; + for _ in 0..5 { + let mut shuffled = edges.clone(); + for k in (1..shuffled.len()).rev() { + state ^= state >> 12; + state ^= state << 25; + state ^= state >> 27; + let m = (state.wrapping_mul(0x2545_F491_4F6C_DD1D) % (k as u64 + 1)) as usize; + shuffled.swap(k, m); + } + let got = compute_scores_from_edges(n, shuffled.into_iter(), 10000, 1e-8); + assert_eq!(got, reference); + } + } + + /// Convergence is reported, never assumed. + #[test] + fn solver_reports_convergence_for_every_shape() { + for n in [2usize, 3, 64, 900] { + let mut edges: Vec<((usize, usize), f64)> = Vec::new(); + for i in 0..n - 1 { + edges.push(((i + 1, i), 2.0)); + edges.push(((i, i + 1), 1.0)); + } + let solved = solve_scores_from_edges(n, edges.into_iter(), 10000, 1e-8); + assert!(solved.converged, "chain n={n} reported non-convergence"); + assert!(solved.residual < 1e-10, "chain n={n}: {:e}", solved.residual); + assert_eq!(solved.pi.len(), n); + assert_eq!(solved.log_pi.len(), n); + } + } + #[test] fn group_ranking_cache_dirty_flow() { let mut g = mk_group(); diff --git a/server/src/stationary.rs b/server/src/stationary.rs new file mode 100644 index 0000000000000000000000000000000000000000..f13af1526e7abdf53fcc02a6f38dee8c6c2257fe --- /dev/null +++ b/server/src/stationary.rs @@ -0,0 +1,1623 @@ +//! Stationary-distribution solvers for the Rank Centrality Markov chain. +//! +//! `ranking::compute_scores_from_edges` builds a chain over compared items and +//! needs π with `πP = π`, `Σπ = 1`. This module owns every candidate solver so +//! they can be compared against each other on identical input, and so the +//! ranking code has exactly one place to ask for "the stationary distribution, +//! and tell me whether you actually got there". +//! +//! # The chain +//! +//! Rank Centrality (Negahban, Oh, Shah 2012, §3.1) uses +//! +//! ```text +//! P_ij = a_ij / d_max (i ≠ j, compared) +//! P_ii = 1 - (Σ_k a_ik) / d_max +//! ``` +//! +//! where `a_ij = A_ij / (A_ij + A_ji)` and `d_max` is the maximum *unweighted* +//! degree. `d_max` is a uniformization constant: it exists only to keep the +//! chain aperiodic so power iteration cannot oscillate (issue #146). Because +//! +//! ```text +//! πP = π ⟺ πQ = 0, Q = P - I, Q_ij = a_ij / d_max (i ≠ j) +//! ``` +//! +//! and `Q` is only defined by its off-diagonal entries (the diagonal is minus +//! the row sum), scaling every off-diagonal by the same `1/d_max` leaves π +//! unchanged. **Direct solvers therefore do not need `d_max` at all** and are +//! structurally immune to the periodicity bug. Only the iterative solvers, +//! which literally walk `P`, care. +//! +//! # Dynamic range +//! +//! On a tree (a chain is a tree) the chain is reversible and detailed balance +//! pins the answer exactly: `π_i / π_j = a_ji / a_ij` for every edge. A chain +//! of `n` items each preferred 2:1 over the next therefore has +//! `π_max / π_min = 2^(n-1)`, which leaves the f64 range at n ≈ 1075. This is a +//! property of the model, not of any solver: past that width the tail of the +//! ranking is not representable in double precision and underflows to zero. +//! [`Solution::underflowed`] reports it, and [`Solution::log_pi`] stays exact +//! there, which is why it — not `pi` — is the sort key for a ranking. +//! +//! # Strategy ([`solve`]) +//! +//! 1. Split disconnected components; each is solved on its own and keeps the +//! share of the mass its node count started with (what power iteration from +//! a uniform start converges to). +//! 2. Sparse GTH state reduction, minimum-degree order, with a degree guard and +//! a work budget. Trees and chains reduce completely in `O(n)`; whatever +//! survives is an irreducible core solved by dense GTH when it fits. +//! 3. Otherwise Gauss–Seidel, then power iteration as a second opinion. Graphs +//! that defeat step 2 are exactly the well-connected ones these finish in +//! tens of sweeps. +//! +//! Nothing returns without a checked residual: see [`Solution::converged`]. + +/// Which solver produced a [`Solution`]. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)] +pub enum Method { + /// Baseline: repeated `π ← πP` until the L1 step falls under `tol`. + Power, + /// Power iteration with periodic Aitken Δ² extrapolation. + PowerAitken, + /// Gauss–Seidel / SOR sweeps on the balance equations. + Sor, + /// Dense LU with partial pivoting on `Qᵀπ = e`, one row replaced by `Σπ = 1`. + DenseLu, + /// Dense Grassmann–Taksar–Heyman state reduction, `O(n³)`, subtraction-free. + DenseGth, + /// GTH state reduction over a sparse graph with minimum-degree elimination. + SparseGth, + /// Jacobi-preconditioned BiCGSTAB on the singular balance system. + BiCgStab, +} + +impl Method { + pub fn label(self) -> &'static str { + match self { + Method::Power => "power", + Method::PowerAitken => "power+aitken", + Method::Sor => "sor", + Method::DenseLu => "dense-lu", + Method::DenseGth => "dense-gth", + Method::SparseGth => "sparse-gth", + Method::BiCgStab => "bicgstab", + } + } + + /// Direct methods answer exactly (up to rounding); iterative ones can stall. + pub fn is_direct(self) -> bool { + matches!(self, Method::DenseLu | Method::DenseGth | Method::SparseGth) + } +} + +/// A stationary distribution plus everything needed to judge whether to trust it. +#[derive(Debug, Clone)] +pub struct Solution { + /// Non-negative, sums to 1 (unless every entry underflowed). + pub pi: Vec, + /// `ln π_i`, shifted so the largest entry is 0. Direct solvers fill this in + /// exactly even where `pi` has flushed to zero, so it is the only faithful + /// sort key on graphs whose score spread exceeds the f64 range. + /// `f64::NEG_INFINITY` for genuinely unreachable states. + pub log_pi: Vec, + pub method: Method, + /// Sweeps performed. Zero for direct methods. + pub iterations: usize, + /// `‖πP − π‖₁`, the honest backward error of the answer that was returned. + pub residual: f64, + /// False means the number that came back is *not* the stationary + /// distribution to the requested tolerance. Never silently true. + pub converged: bool, + /// The exact answer spans more than f64 can hold, so the tail of the + /// ranking has been flushed to zero and its internal order is lost. + pub underflowed: bool, +} + +/// Tuning for [`solve`]. Defaults encode the recommended hybrid strategy. +#[derive(Debug, Clone, Copy)] +pub struct SolveOptions { + pub tol: f64, + pub max_iters: usize, + /// Row-merge budget for the sparse elimination in [`sparse_gth`]. Bounds + /// the work wasted before giving up on the direct path. + pub direct_work_budget: u64, + /// Largest irreducible core still worth an `O(n³)` dense GTH. + pub dense_core_max: usize, +} + +impl Default for SolveOptions { + fn default() -> Self { + SolveOptions { + tol: 1e-8, + max_iters: 10_000, + direct_work_budget: DEFAULT_DIRECT_WORK_BUDGET, + dense_core_max: DEFAULT_DENSE_CORE_MAX, + } + } +} + +/// Measured: sparse elimination runs at roughly 20–40 M row-merge steps per +/// second, so 2 M caps the abandoned work at a few tens of milliseconds. Trees +/// and chains finish two to three orders of magnitude under it. +pub const DEFAULT_DIRECT_WORK_BUDGET: u64 = 2_000_000; + +/// Dense GTH is `O(n³)`; measured at ~5 ms for n=400 and ~40 ms for n=800, so +/// 512 keeps the direct path inside a garden render's budget. Above it, the +/// iterative solvers are both faster and (on well-connected graphs) accurate. +pub const DEFAULT_DENSE_CORE_MAX: usize = 512; + +/// The comparison chain in the one form every solver consumes. +/// +/// `rows[i]` holds `(j, a_ij)` for the pairwise-normalized weights, sorted by +/// `j` so that every floating-point summation in this module runs in a fixed +/// order regardless of how the caller's `HashMap` happened to iterate. +#[derive(Debug, Clone)] +pub struct RankChain { + pub n: usize, + pub rows: Vec>, + pub row_sum: Vec, + /// Maximum unweighted degree; the uniformization constant for `P`. + pub d_max: f64, +} + +impl RankChain { + /// Build from pairwise-normalized weights. `edges` may arrive in any order; + /// they are sorted here so results are bit-for-bit reproducible, and + /// repeated `(i, j)` pairs are coalesced so every row holds each column + /// exactly once (elimination relies on that). + pub fn from_normalized(n: usize, mut edges: Vec<((usize, usize), f64)>, d_max: usize) -> Self { + edges.sort_unstable_by(|a, b| a.0.cmp(&b.0)); + let mut rows: Vec> = vec![Vec::new(); n]; + for ((i, j), w) in edges { + // Written positively so a NaN weight is dropped rather than kept. + let usable = i < n && j < n && i != j && w > 0.0; + if !usable { + continue; + } + match rows[i].last_mut() { + Some(last) if last.0 == j => last.1 += w, + _ => rows[i].push((j, w)), + } + } + let row_sum: Vec = rows + .iter() + .map(|r| r.iter().map(|&(_, w)| w).sum()) + .collect(); + RankChain { + n, + rows, + row_sum, + d_max: d_max as f64, + } + } + + pub fn nnz(&self) -> usize { + self.rows.iter().map(|r| r.len()).sum() + } + + /// Incoming edges, `cols[j] = [(i, a_ij)]`, sorted by `i`. + fn columns(&self) -> Vec> { + let mut cols: Vec> = vec![Vec::new(); self.n]; + for i in 0..self.n { + for &(j, w) in &self.rows[i] { + cols[j].push((i, w)); + } + } + cols + } + + /// `‖πP − π‖₁` for an arbitrary vector, using the uniformized `P`. + /// + /// This is the yardstick every method is judged by, so it must not depend + /// on how the candidate was produced. + pub fn residual(&self, pi: &[f64]) -> f64 { + let n = self.n; + let mut next = vec![0.0f64; n]; + for i in 0..n { + let stay = (self.d_max - self.row_sum[i]) / self.d_max; + next[i] += pi[i] * stay; + for &(j, w) in &self.rows[i] { + next[j] += pi[i] * (w / self.d_max); + } + } + (0..n).map(|i| (next[i] - pi[i]).abs()).sum() + } +} + +/// A `Solution` for an input with nothing to solve (no items, or one). +pub fn trivial(pi: Vec) -> Solution { + Solution { + log_pi: vec![0.0; pi.len()], + pi, + method: Method::SparseGth, + iterations: 0, + residual: 0.0, + converged: true, + underflowed: false, + } +} + +/// Normalize in place to sum 1. Returns whether the tail underflowed to zero. +fn normalize(pi: &mut [f64]) -> bool { + let sum: f64 = pi.iter().sum(); + if sum.is_finite() && sum > 0.0 { + for p in pi.iter_mut() { + *p /= sum; + } + } else { + let uniform = 1.0 / pi.len() as f64; + for p in pi.iter_mut() { + *p = uniform; + } + return true; + } + pi.iter().any(|&p| p <= 0.0) +} + +fn finish( + pi: Vec, + chain: &RankChain, + method: Method, + iterations: usize, + tol: f64, +) -> Solution { + let log_pi = log_of(&pi); + finish_with_logs(pi, log_pi, chain, method, iterations, tol) +} + +fn log_of(pi: &[f64]) -> Vec { + let max = pi.iter().cloned().fold(0.0f64, f64::max); + if max <= 0.0 { + return vec![f64::NEG_INFINITY; pi.len()]; + } + pi.iter().map(|&p| (p / max).ln()).collect() +} + +fn finish_with_logs( + pi: Vec, + log_pi: Vec, + chain: &RankChain, + method: Method, + iterations: usize, + tol: f64, +) -> Solution { + let mut pi = pi; + let underflowed = normalize(&mut pi); + let residual = chain.residual(&pi); + Solution { + converged: residual.is_finite() && residual <= tol.max(f64::EPSILON * chain.n as f64), + pi, + log_pi, + method, + iterations, + residual, + underflowed, + } +} + +/// `ln(Σ exp(t))` computed by shifting out the largest term, so a sum of +/// astronomically different magnitudes neither overflows nor loses the small ones. +fn log_sum_exp(terms: &[f64]) -> f64 { + let max = terms.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + if !max.is_finite() { + return max; + } + let sum: f64 = terms.iter().map(|&t| (t - max).exp()).sum(); + max + sum.ln() +} + +/// Rebuild `pi` from exact log-scores, shifted so the maximum is 1. +/// +/// `exp` degrades into subnormals rather than jumping to zero, so this keeps +/// roughly 745 decades of usable spread where the linear back-substitution +/// keeps 308 — the difference between an exact ranking of a 2400-long chain and +/// an exact ranking of a 1000-long one. +fn pi_from_logs(log_pi: &[f64]) -> Vec { + let max = log_pi.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + if !max.is_finite() { + return vec![1.0 / log_pi.len() as f64; log_pi.len()]; + } + log_pi.iter().map(|&l| (l - max).exp()).collect() +} + +// --------------------------------------------------------------------------- +// Iterative solvers +// --------------------------------------------------------------------------- + +/// One `π ← πP` sweep. +fn step(chain: &RankChain, pi: &[f64], next: &mut [f64]) { + next.fill(0.0); + for i in 0..chain.n { + let stay = (chain.d_max - chain.row_sum[i]) / chain.d_max; + next[i] += pi[i] * stay; + for &(j, w) in &chain.rows[i] { + next[j] += pi[i] * (w / chain.d_max); + } + } +} + +/// Baseline power iteration: exactly what `ranking.rs` does today. +pub fn power(chain: &RankChain, opts: SolveOptions) -> Solution { + let n = chain.n; + let mut pi = vec![1.0 / n as f64; n]; + let mut next = vec![0.0f64; n]; + let mut iters = 0usize; + + for _ in 0..opts.max_iters { + step(chain, &pi, &mut next); + iters += 1; + let diff: f64 = (0..n).map(|i| (pi[i] - next[i]).abs()).sum(); + pi.copy_from_slice(&next); + if diff < opts.tol { + break; + } + } + finish(pi, chain, Method::Power, iters, opts.tol) +} + +/// Power iteration with periodic componentwise Aitken Δ² extrapolation. +/// +/// Aitken assumes the error is dominated by a single geometric mode, which is +/// exactly the regime that makes plain power iteration slow. The extrapolated +/// vector is only accepted when it lowers the residual, so a bad extrapolation +/// costs one wasted sweep instead of divergence. +pub fn power_aitken(chain: &RankChain, opts: SolveOptions) -> Solution { + let n = chain.n; + let mut x0 = vec![1.0 / n as f64; n]; + let mut x1 = vec![0.0f64; n]; + let mut x2 = vec![0.0f64; n]; + let mut cand = vec![0.0f64; n]; + let mut iters = 0usize; + let mut best_resid = f64::INFINITY; + + while iters < opts.max_iters { + step(chain, &x0, &mut x1); + step(chain, &x1, &mut x2); + iters += 2; + + let diff: f64 = (0..n).map(|i| (x2[i] - x1[i]).abs()).sum(); + if diff < opts.tol { + return finish(x2, chain, Method::PowerAitken, iters, opts.tol); + } + + // Aitken: x* ≈ x2 - (Δx1)² / Δ²x0, componentwise. + let mut usable = true; + for i in 0..n { + let d1 = x1[i] - x0[i]; + let d2 = x2[i] - x1[i]; + let denom = d2 - d1; + if denom.abs() <= f64::MIN_POSITIVE { + cand[i] = x2[i]; + continue; + } + let v = x2[i] - d2 * d2 / denom; + if !v.is_finite() || v < 0.0 { + usable = false; + break; + } + cand[i] = v; + } + + std::mem::swap(&mut x0, &mut x2); + if !usable { + continue; + } + let sum: f64 = cand.iter().sum(); + if !(sum.is_finite() && sum > 0.0) { + continue; + } + for c in cand.iter_mut() { + *c /= sum; + } + let r_cand = chain.residual(&cand); + let r_cur = chain.residual(&x0); + if r_cand < r_cur && r_cand < best_resid { + best_resid = r_cand; + x0.copy_from_slice(&cand); + if r_cand <= opts.tol { + return finish(x0, chain, Method::PowerAitken, iters, opts.tol); + } + } + } + finish(x0, chain, Method::PowerAitken, iters, opts.tol) +} + +/// Gauss–Seidel / SOR on the balance equations. +/// +/// Balance says `π_j · Σ_k a_jk = Σ_{i≠j} π_i a_ij`, i.e. every node's outflow +/// equals its inflow. Sweeping that in place (using already-updated components +/// within the sweep) propagates information across the whole graph in one pass +/// instead of one edge per pass, which is precisely what power iteration cannot +/// do on a long chain. +pub fn sor(chain: &RankChain, opts: SolveOptions, omega: f64) -> Solution { + let n = chain.n; + let cols = chain.columns(); + let mut pi = vec![1.0 / n as f64; n]; + let mut iters = 0usize; + + for _ in 0..opts.max_iters { + iters += 1; + let mut delta = 0.0f64; + for j in 0..n { + let out = chain.row_sum[j]; + if out <= 0.0 { + continue; + } + let inflow: f64 = cols[j].iter().map(|&(i, w)| pi[i] * w).sum(); + let gs = inflow / out; + let updated = (1.0 - omega) * pi[j] + omega * gs; + let updated = if updated.is_finite() && updated > 0.0 { + updated + } else { + pi[j] + }; + delta += (updated - pi[j]).abs(); + pi[j] = updated; + } + // The iteration is only defined up to scale; renormalize so `tol` means + // the same thing here as it does for power iteration. + let sum: f64 = pi.iter().sum(); + if sum.is_finite() && sum > 0.0 { + for p in pi.iter_mut() { + *p /= sum; + } + } + if delta < opts.tol { + break; + } + } + finish(pi, chain, Method::Sor, iters, opts.tol) +} + +/// Jacobi-preconditioned BiCGSTAB on `Qᵀπ = 0` with row 0 replaced by `Σπ = 1`. +/// +/// Included as the Krylov representative. GMRES/Arnoldi would need restart +/// bookkeeping and a stored basis for the same job; BiCGSTAB gives the same +/// short-recurrence answer for a nonsymmetric operator at fixed memory. +pub fn bicgstab(chain: &RankChain, opts: SolveOptions) -> Solution { + let n = chain.n; + let cols = chain.columns(); + + // A x = b where row 0 is Σx = 1 and row j>0 is (inflow - outflow) at j. + let apply = |x: &[f64], out: &mut [f64]| { + out[0] = x.iter().sum(); + for j in 1..n { + let inflow: f64 = cols[j].iter().map(|&(i, w)| x[i] * w).sum(); + out[j] = inflow - x[j] * chain.row_sum[j]; + } + }; + // Jacobi preconditioner: the diagonal of that operator. + let diag: Vec = (0..n) + .map(|j| { + let d = if j == 0 { 1.0 } else { -chain.row_sum[j] }; + if d.abs() < 1e-300 { + 1.0 + } else { + d + } + }) + .collect(); + let precond = |v: &[f64], out: &mut [f64]| { + for i in 0..n { + out[i] = v[i] / diag[i]; + } + }; + + let mut b = vec![0.0f64; n]; + b[0] = 1.0; + + let mut x = vec![1.0 / n as f64; n]; + let mut ax = vec![0.0f64; n]; + apply(&x, &mut ax); + let mut r: Vec = (0..n).map(|i| b[i] - ax[i]).collect(); + let r_hat = r.clone(); + + let mut rho = 1.0f64; + let mut alpha = 1.0f64; + let mut omega = 1.0f64; + let mut v = vec![0.0f64; n]; + let mut p = vec![0.0f64; n]; + let mut y = vec![0.0f64; n]; + let mut z = vec![0.0f64; n]; + let mut s = vec![0.0f64; n]; + let mut t = vec![0.0f64; n]; + let dot = |a: &[f64], b: &[f64]| -> f64 { (0..a.len()).map(|i| a[i] * b[i]).sum() }; + let bnorm = 1.0f64; + let mut iters = 0usize; + + for _ in 0..opts.max_iters { + iters += 1; + let rho_new = dot(&r_hat, &r); + if rho_new.abs() < 1e-300 { + break; + } + let beta = (rho_new / rho) * (alpha / omega); + rho = rho_new; + for i in 0..n { + p[i] = r[i] + beta * (p[i] - omega * v[i]); + } + precond(&p, &mut y); + apply(&y, &mut v); + let denom = dot(&r_hat, &v); + if denom.abs() < 1e-300 { + break; + } + alpha = rho / denom; + for i in 0..n { + s[i] = r[i] - alpha * v[i]; + } + if dot(&s, &s).sqrt() / bnorm < opts.tol * 1e-2 { + for i in 0..n { + x[i] += alpha * y[i]; + } + break; + } + precond(&s, &mut z); + apply(&z, &mut t); + let tt = dot(&t, &t); + if tt.abs() < 1e-300 { + break; + } + omega = dot(&t, &s) / tt; + for i in 0..n { + x[i] += alpha * y[i] + omega * z[i]; + r[i] = s[i] - omega * t[i]; + } + if dot(&r, &r).sqrt() / bnorm < opts.tol * 1e-2 { + break; + } + if omega.abs() < 1e-300 { + break; + } + } + + // BiCGSTAB has no sign constraint; clamp before normalizing. + for xi in x.iter_mut() { + if !xi.is_finite() || *xi < 0.0 { + *xi = 0.0; + } + } + finish(x, chain, Method::BiCgStab, iters, opts.tol) +} + +// --------------------------------------------------------------------------- +// Direct solvers +// --------------------------------------------------------------------------- + +/// Dense LU with partial pivoting on the balance system. +/// +/// `πQ = 0` is rank `n-1`, so equation 0 is replaced by the normalization +/// `Σπ = 1`, giving a nonsingular `n×n` system solved in `O(n³)`. +pub fn dense_lu(chain: &RankChain, tol: f64) -> Solution { + let n = chain.n; + // Row-major A, where A[j][i] is the coefficient of π_i in equation j. + let mut a = vec![0.0f64; n * n]; + for v in a.iter_mut().take(n) { + *v = 1.0; // equation 0: Σ π_i = 1 + } + for i in 0..n { + for &(j, w) in &chain.rows[i] { + if j != 0 { + a[j * n + i] += w; // inflow to j from i + } + } + if i != 0 { + a[i * n + i] -= chain.row_sum[i]; // outflow from i + } + } + let mut b = vec![0.0f64; n]; + b[0] = 1.0; + + // Gaussian elimination, partial pivoting. + let mut perm: Vec = (0..n).collect(); + for k in 0..n { + let mut piv = k; + let mut best = a[perm[k] * n + k].abs(); + for r in (k + 1)..n { + let v = a[perm[r] * n + k].abs(); + if v > best { + best = v; + piv = r; + } + } + if best == 0.0 { + continue; + } + perm.swap(k, piv); + let pk = perm[k]; + let pivot = a[pk * n + k]; + for r in (k + 1)..n { + let pr = perm[r]; + let f = a[pr * n + k] / pivot; + if f == 0.0 { + continue; + } + a[pr * n + k] = 0.0; + for c in (k + 1)..n { + a[pr * n + c] -= f * a[pk * n + c]; + } + b[pr] -= f * b[pk]; + } + } + let mut x = vec![0.0f64; n]; + for k in (0..n).rev() { + let pk = perm[k]; + let mut acc = b[pk]; + for c in (k + 1)..n { + acc -= a[pk * n + c] * x[c]; + } + let d = a[pk * n + k]; + x[k] = if d.abs() > 0.0 { acc / d } else { 0.0 }; + } + for xi in x.iter_mut() { + if !xi.is_finite() || *xi < 0.0 { + *xi = 0.0; + } + } + finish(x, chain, Method::DenseLu, 0, tol) +} + +/// Rescale threshold for GTH back-substitution. +/// +/// Back-substitution accumulates `π` in *unnormalized* form, and on a strongly +/// ordered chain the values grow geometrically. Rescaling the partial vector +/// whenever it nears the top of the f64 range keeps every ratio exact (they are +/// all defined up to one global scalar) and avoids returning `inf`. +const GTH_RESCALE_ABOVE: f64 = 1e250; + +/// Dense Grassmann–Taksar–Heyman state reduction on a compact rate matrix. +/// +/// Every operation is an addition of non-negative numbers or a division by a +/// positive sum: there is no subtraction anywhere, so there is no cancellation +/// and no pivoting is needed. That is why GTH is the reference implementation +/// here even though dense LU costs the same `O(n³)`. +/// +/// Returns `ln π` shifted so the largest entry is 0. Finishing in logs matters +/// because the answer itself can span more than f64's range. +fn gth_dense_logs(n: usize, rows: &[Vec<(usize, f64)>]) -> Vec { + if n == 0 { + return vec![]; + } + if n == 1 { + return vec![0.0]; + } + let mut a = vec![0.0f64; n * n]; + for (i, row) in rows.iter().enumerate().take(n) { + for &(j, w) in row { + if j < n && j != i { + a[i * n + j] += w; + } + } + } + + for k in (1..n).rev() { + let s: f64 = (0..k).map(|j| a[k * n + j]).sum(); + let reaches_survivors = s > 0.0; + if !reaches_survivors { + // State k cannot reach the surviving block, so it feeds no mass + // back into it. + for i in 0..k { + a[i * n + k] = 0.0; + } + continue; + } + for i in 0..k { + let f = a[i * n + k] / s; + a[i * n + k] = f; + if f == 0.0 { + continue; + } + for j in 0..k { + a[i * n + j] += f * a[k * n + j]; + } + } + } + + // Cheap path: accumulate linearly, take logs at the end. + let mut x = vec![0.0f64; n]; + x[0] = 1.0; + let mut rescaled = false; + for k in 1..n { + let mut acc = 0.0f64; + for i in 0..k { + acc += x[i] * a[i * n + k]; + } + x[k] = acc; + if acc > GTH_RESCALE_ABOVE { + let inv = 1.0 / acc; + for v in x.iter_mut().take(k + 1) { + *v *= inv; + } + rescaled = true; + } + } + if !rescaled { + return shift_logs(x.iter().map(|&v| v.ln()).collect()); + } + + // The span overflowed a single f64 vector, so redo the accumulation in log + // space. The elimination coefficients are unchanged; only the sum differs. + let mut logx = vec![f64::NEG_INFINITY; n]; + logx[0] = 0.0; + let mut terms: Vec = Vec::with_capacity(n); + for k in 1..n { + terms.clear(); + for i in 0..k { + let c = a[i * n + k]; + if c > 0.0 && logx[i].is_finite() { + terms.push(logx[i] + c.ln()); + } + } + logx[k] = log_sum_exp(&terms); + } + shift_logs(logx) +} + +/// Dense GTH over the whole chain. `O(n³)` time, `O(n²)` memory. +pub fn dense_gth(chain: &RankChain, tol: f64) -> Solution { + let logs = gth_dense_logs(chain.n, &chain.rows); + let pi = pi_from_logs(&logs); + finish_with_logs(pi, logs, chain, Method::DenseGth, 0, tol) +} + +/// Shift log-scores so the maximum is exactly 0. +fn shift_logs(mut logs: Vec) -> Vec { + let max = logs.iter().cloned().fold(f64::NEG_INFINITY, f64::max); + if max.is_finite() { + for l in logs.iter_mut() { + *l -= max; + } + } + logs +} + +/// Outcome of attempting the sparse direct path. +pub enum SparseGthOutcome { + Solved(Solution), + /// Elimination could not finish cheaply: states were left with degree above + /// the guard, or the work budget ran out. Fall back to an iterative solver. + TooDense { + core: usize, + work: u64, + }, +} + +/// `μ = m - n + c` over the symmetrized comparison graph: the number of +/// independent cycles, i.e. how far the graph is from being a forest. `O(n + m)`. +fn cyclomatic_number(chain: &RankChain) -> usize { + let n = chain.n; + let mut parent: Vec = (0..n).collect(); + fn find(parent: &mut [usize], mut x: usize) -> usize { + while parent[x] != x { + parent[x] = parent[parent[x]]; + x = parent[x]; + } + x + } + + let mut edges = 0usize; + let mut merges = 0usize; + for i in 0..n { + for &(j, _) in &chain.rows[i] { + // Count each compared pair once: the reverse arc exists whenever + // both sides got votes, and the pair is an edge either way. + if j < i && chain.rows[j].binary_search_by(|p| p.0.cmp(&i)).is_ok() { + continue; + } + edges += 1; + let (ri, rj) = (find(&mut parent, i), find(&mut parent, j)); + if ri != rj { + parent[ri] = rj; + merges += 1; + } + } + } + // c = n - merges, so μ = m - n + c = m - merges. + edges.saturating_sub(merges) +} + +/// Above this live degree, eliminating a state creates more fill than it +/// removes. Trees, chains, series-parallel graphs and most hand-built +/// ontologies reduce away completely under this guard; random graphs stall +/// almost immediately, which is exactly the signal to go iterative. +const MAX_ELIM_DEGREE: usize = 12; + +/// GTH state reduction over the sparse graph, eliminating minimum-degree states first. +/// +/// Elimination on a graph is fill-bounded: removing a state connects its +/// surviving in-neighbours to its surviving out-neighbours. A chain — or any +/// tree — has an elimination order with *zero* fill, so this is `O(n)` on +/// exactly the topology where power iteration needs `Θ(n²)` sweeps. The two +/// failure modes are complementary: what stalls here (dense, well-connected +/// graphs) is what power iteration finishes in a few dozen sweeps. +/// +/// States that cannot be eliminated cheaply form an irreducible *core*, handed +/// to dense GTH when it is small enough and reported as `TooDense` otherwise. +pub fn sparse_gth(chain: &RankChain, opts: SolveOptions, tol: f64) -> SparseGthOutcome { + let n = chain.n; + if n == 1 { + return SparseGthOutcome::Solved(finish(vec![1.0], chain, Method::SparseGth, 0, tol)); + } + + // Predict the core before doing any work. A graph whose every vertex has + // degree ≥ 3 has at most `2μ - 2` vertices, where `μ = m - n + c` is the + // cyclomatic number, so `2μ` bounds what elimination can leave behind. A + // forest has `μ = 0` and reduces to nothing; a random graph has `μ ≈ m` and + // reduces to almost nothing *but itself*, which is the case worth skipping + // before paying for it. + if n > opts.dense_core_max && 2 * cyclomatic_number(chain) > opts.dense_core_max { + return SparseGthOutcome::TooDense { core: n, work: 0 }; + } + + // Live out-rates, plus in-neighbour lists that may carry stale or duplicate + // entries — they are compacted against `alive` when a state is popped. + let mut out: Vec> = chain.rows.clone(); + let mut inn: Vec> = vec![Vec::new(); n]; + for i in 0..n { + for &(j, _) in &chain.rows[i] { + inn[j].push(i); + } + } + + let mut alive = vec![true; n]; + // Position of each column in the row currently being updated, -1 when + // absent. Turns a row merge into O(len) instead of a map lookup per entry. + let mut slot: Vec = vec![-1; n]; + let mut recovery: Vec<(usize, Vec<(usize, f64)>)> = Vec::with_capacity(n); + let mut work: u64 = 0; + + let mut heap: std::collections::BinaryHeap> = + std::collections::BinaryHeap::new(); + for k in 0..n { + heap.push(std::cmp::Reverse((out[k].len() + inn[k].len(), k))); + } + + let mut remaining = n; + while remaining > 1 { + let Some(std::cmp::Reverse((d, k))) = heap.pop() else { + break; + }; + if !alive[k] { + continue; + } + + // Compact k's neighbour lists, then re-check the heap key against the + // true degree: a stale key means another elimination has changed it. + out[k].retain(|&(j, w)| alive[j] && w > 0.0); + inn[k].sort_unstable(); + inn[k].dedup(); + inn[k].retain(|&i| alive[i] && i != k); + let cur = out[k].len() + inn[k].len(); + if cur != d { + heap.push(std::cmp::Reverse((cur, k))); + continue; + } + if out[k].len().max(inn[k].len()) > MAX_ELIM_DEGREE { + // The cheapest remaining state is already expensive, so every other + // one is too. What is left is the core. + break; + } + + let s: f64 = out[k].iter().map(|&(_, w)| w).sum(); + let in_list = std::mem::take(&mut inn[k]); + let out_list = std::mem::take(&mut out[k]); + + let mut coeffs: Vec<(usize, f64)> = Vec::with_capacity(in_list.len()); + for &i in &in_list { + for (p, &(j, _)) in out[i].iter().enumerate() { + slot[j] = p as i64; + } + work = work.saturating_add(out[i].len() as u64 + out_list.len() as u64); + + let w_ik = if slot[k] >= 0 { + let p = slot[k] as usize; + let v = out[i][p].1; + out[i][p].1 = 0.0; + v + } else { + 0.0 + }; + let f = if s > 0.0 { w_ik / s } else { 0.0 }; + if f > 0.0 { + coeffs.push((i, f)); + for &(j, r) in &out_list { + if j == i { + continue; + } + if slot[j] >= 0 { + out[i][slot[j] as usize].1 += f * r; + } else { + slot[j] = out[i].len() as i64; + out[i].push((j, f * r)); + inn[j].push(i); + } + } + } + for &(j, _) in out[i].iter() { + slot[j] = -1; + } + slot[k] = -1; + out[i].retain(|&(j, w)| j != k && w > 0.0); + } + + alive[k] = false; + remaining -= 1; + recovery.push((k, coeffs)); + + for &i in &in_list { + if alive[i] { + heap.push(std::cmp::Reverse((out[i].len() + inn[i].len(), i))); + } + } + for &(j, _) in &out_list { + if alive[j] { + heap.push(std::cmp::Reverse((out[j].len() + inn[j].len(), j))); + } + } + + if work > opts.direct_work_budget { + return SparseGthOutcome::TooDense { + core: (0..n).filter(|&i| alive[i]).count(), + work, + }; + } + } + + let core: Vec = (0..n).filter(|&i| alive[i]).collect(); + if core.len() > opts.dense_core_max { + return SparseGthOutcome::TooDense { + core: core.len(), + work, + }; + } + + // Solve the irreducible core densely, then unwind the eliminated states. + let mut compact = vec![usize::MAX; n]; + for (c, &g) in core.iter().enumerate() { + compact[g] = c; + } + let core_rows: Vec> = core + .iter() + .map(|&g| { + let mut row: Vec<(usize, f64)> = out[g] + .iter() + .filter(|&&(j, w)| alive[j] && w > 0.0) + .map(|&(j, w)| (compact[j], w)) + .collect(); + row.sort_unstable_by(|a, b| a.0.cmp(&b.0)); + row + }) + .collect(); + + let core_logs = gth_dense_logs(core.len(), &core_rows); + let mut logx = vec![f64::NEG_INFINITY; n]; + for (c, &g) in core.iter().enumerate() { + logx[g] = core_logs[c]; + } + let mut terms: Vec = Vec::new(); + for (k, coeffs) in recovery.iter().rev() { + terms.clear(); + for &(i, f) in coeffs { + if f > 0.0 && logx[i].is_finite() { + terms.push(logx[i] + f.ln()); + } + } + logx[*k] = log_sum_exp(&terms); + } + + let shifted = shift_logs(logx); + let pi = pi_from_logs(&shifted); + SparseGthOutcome::Solved(finish_with_logs( + pi, + shifted, + chain, + Method::SparseGth, + 0, + tol, + )) +} + +// --------------------------------------------------------------------------- +// Recommended strategy +// --------------------------------------------------------------------------- + +/// Solve for π using the hybrid strategy: sparse direct first, iterative only +/// when the direct path is genuinely too expensive. +/// +/// This never returns an unconverged answer without saying so — check +/// [`Solution::converged`]. +pub fn solve(chain: &RankChain, opts: SolveOptions) -> Solution { + if chain.n == 0 { + return Solution { + pi: vec![], + log_pi: vec![], + method: Method::SparseGth, + iterations: 0, + residual: 0.0, + converged: true, + underflowed: false, + }; + } + if chain.n == 1 { + return Solution { + pi: vec![1.0], + log_pi: vec![0.0], + method: Method::SparseGth, + iterations: 0, + residual: 0.0, + converged: true, + underflowed: false, + }; + } + if chain.nnz() == 0 { + let n = chain.n; + return Solution { + pi: vec![1.0 / n as f64; n], + log_pi: vec![0.0; n], + method: Method::SparseGth, + iterations: 0, + residual: 0.0, + converged: true, + underflowed: false, + }; + } + + // A chain over several disconnected components has no unique stationary + // distribution: any split of mass between them is stationary. Power + // iteration from a uniform start picks the one where each component holds + // its share of the nodes, so reproduce that explicitly instead of letting a + // direct solver pick one component and zero the rest. + let comps = weak_components(chain); + if comps.len() > 1 { + return solve_by_component(chain, &comps, opts); + } + solve_connected(chain, opts) +} + +/// Weakly connected components of the comparison graph, each sorted ascending. +fn weak_components(chain: &RankChain) -> Vec> { + let n = chain.n; + let mut adj: Vec> = vec![Vec::new(); n]; + for i in 0..n { + for &(j, _) in &chain.rows[i] { + adj[i].push(j); + adj[j].push(i); + } + } + let mut seen = vec![false; n]; + let mut comps = Vec::new(); + let mut stack = Vec::new(); + for start in 0..n { + if seen[start] { + continue; + } + seen[start] = true; + stack.push(start); + let mut comp = Vec::new(); + while let Some(x) = stack.pop() { + comp.push(x); + for &y in &adj[x] { + if !seen[y] { + seen[y] = true; + stack.push(y); + } + } + } + comp.sort_unstable(); + comps.push(comp); + } + comps +} + +fn solve_by_component(chain: &RankChain, comps: &[Vec], opts: SolveOptions) -> Solution { + let n = chain.n; + let mut log_pi = vec![f64::NEG_INFINITY; n]; + let mut converged = true; + let mut iterations = 0usize; + let mut method = Method::SparseGth; + + for comp in comps { + let mut compact = vec![usize::MAX; n]; + for (c, &g) in comp.iter().enumerate() { + compact[g] = c; + } + let rows: Vec> = comp + .iter() + .map(|&g| { + chain.rows[g] + .iter() + .filter(|&&(j, _)| compact[j] != usize::MAX) + .map(|&(j, w)| (compact[j], w)) + .collect() + }) + .collect(); + let row_sum = rows + .iter() + .map(|r| r.iter().map(|&(_, w)| w).sum()) + .collect(); + let sub = RankChain { + n: comp.len(), + rows, + row_sum, + d_max: chain.d_max, + }; + let sol = solve_connected(&sub, opts); + converged &= sol.converged; + iterations = iterations.max(sol.iterations); + if !sol.method.is_direct() { + method = sol.method; + } + // Each component keeps the share of the mass its node count started + // with, spread internally by its own stationary distribution. + let share = (comp.len() as f64 / n as f64).ln() - log_sum_exp(&sol.log_pi); + for (c, &g) in comp.iter().enumerate() { + log_pi[g] = sol.log_pi[c] + share; + } + } + + let shifted = shift_logs(log_pi); + let pi = pi_from_logs(&shifted); + let mut out = finish_with_logs(pi, shifted, chain, method, iterations, opts.tol); + // One component missing tolerance condemns the whole vector, even if the + // combined residual happens to look small next to the dominant component. + out.converged &= converged; + out +} + +fn solve_connected(chain: &RankChain, opts: SolveOptions) -> Solution { + match sparse_gth(chain, opts, opts.tol) { + SparseGthOutcome::Solved(sol) if sol.converged => return sol, + // A direct solve that misses tolerance means the conditioning is worse + // than the elimination could handle; let the iterative path try rather + // than hand back an answer nothing has checked. + SparseGthOutcome::Solved(_) | SparseGthOutcome::TooDense { .. } => {} + } + + // Gauss–Seidel first: on every graph dense enough to reach this branch it + // converged in tens of sweeps and beat power iteration on both time and + // residual. Power iteration stays as the second opinion. + let sweeps = sor(chain, opts, 1.0); + if sweeps.converged { + return sweeps; + } + let pow = power_aitken(chain, opts); + if pow.converged || pow.residual < sweeps.residual { + pow + } else { + sweeps + } +} + +#[cfg(test)] +mod tests { + use super::*; + + const TOL: f64 = 1e-8; + + fn opts() -> SolveOptions { + SolveOptions::default() + } + + /// Build a chain the way production does: raw directed vote weights through + /// `ranking::chain_from_edges`. `votes` are `(winner, loser, w_win, w_lose)`; + /// Rank Centrality walks *toward* the winner, so the loser's row carries the + /// weight (see `GroupState::apply_vote`). + fn chain_from_votes(n: usize, votes: &[(usize, usize, f64, f64)]) -> RankChain { + let mut raw: std::collections::HashMap<(usize, usize), f64> = Default::default(); + for &(w, l, ww, wl) in votes { + *raw.entry((l, w)).or_insert(0.0) += ww; + *raw.entry((w, l)).or_insert(0.0) += wl; + } + crate::ranking::chain_from_edges(n, raw.into_iter()) + } + + fn path_graph(n: usize, ratio: f64) -> RankChain { + let votes: Vec<(usize, usize, f64, f64)> = + (0..n - 1).map(|i| (i, i + 1, ratio, 1.0)).collect(); + chain_from_votes(n, &votes) + } + + fn star_graph(n: usize) -> RankChain { + let votes: Vec<(usize, usize, f64, f64)> = + (1..n).map(|i| (0, i, 2.0 + (i % 5) as f64, 1.0)).collect(); + chain_from_votes(n, &votes) + } + + fn clique_graph(n: usize) -> RankChain { + let mut votes = Vec::new(); + let mut k = 0; + for i in 0..n { + for j in (i + 1)..n { + votes.push((i, j, 2.0 + (k % 5) as f64, 1.0)); + k += 1; + } + } + chain_from_votes(n, &votes) + } + + fn sparse_graph(n: usize, degree: usize) -> RankChain { + let mut state = 0x00C0_FFEEu64; + let mut next = || { + state ^= state >> 12; + state ^= state << 25; + state ^= state >> 27; + state.wrapping_mul(0x2545_F491_4F6C_DD1D) + }; + let mut votes: Vec<(usize, usize, f64, f64)> = + (0..n - 1).map(|i| (i, i + 1, 2.0, 1.0)).collect(); + for k in 0..(degree.saturating_sub(2) * n / 2) { + let a = (next() % n as u64) as usize; + let b = (next() % n as u64) as usize; + if a != b { + votes.push((a, b, 2.0 + (k % 5) as f64, 1.0)); + } + } + chain_from_votes(n, &votes) + } + + fn max_rel_error(got: &[f64], want: &[f64]) -> f64 { + (0..want.len()) + .filter(|&i| want[i] > 1e-280) + .map(|i| (got[i] - want[i]).abs() / want[i]) + .fold(0.0f64, f64::max) + } + + fn force_sparse(chain: &RankChain) -> Solution { + let unlimited = SolveOptions { + direct_work_budget: u64::MAX, + dense_core_max: usize::MAX, + ..opts() + }; + match sparse_gth(chain, unlimited, TOL) { + SparseGthOutcome::Solved(s) => s, + SparseGthOutcome::TooDense { .. } => panic!("uncapped elimination gave up"), + } + } + + // ----------------------------------------------------------------- + // Closed forms: these pin the answer without trusting any solver. + // ----------------------------------------------------------------- + + /// A tree is reversible, so detailed balance fixes π edge by edge: + /// `π_i / π_j = a_ji / a_ij`. For a path at a constant ratio that is a + /// geometric sequence, known exactly. + #[test] + fn path_matches_detailed_balance_closed_form() { + for n in [2usize, 3, 8, 64, 300] { + let chain = path_graph(n, 3.0); + let mut want: Vec = (0..n).map(|i| 3f64.powi(-(i as i32))).collect(); + let sum: f64 = want.iter().sum(); + for w in want.iter_mut() { + *w /= sum; + } + for sol in [ + force_sparse(&chain), + dense_gth(&chain, TOL), + dense_lu(&chain, TOL), + solve(&chain, opts()), + ] { + assert!( + max_rel_error(&sol.pi, &want) < 1e-12, + "{} wrong on path n={n}: {:?}", + sol.method.label(), + &sol.pi[..4.min(n)] + ); + assert!(sol.converged, "{} not converged", sol.method.label()); + } + } + } + + /// The 3-cycle a>b>c>a is symmetric under rotation, so π must be uniform. + /// This is the `test/fixtures/ranking/cycle.sorter` anchor. + #[test] + fn three_cycle_is_exactly_uniform() { + let chain = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (1, 2, 2.0, 1.0), (2, 0, 2.0, 1.0)]); + let sol = solve(&chain, opts()); + for (i, &p) in sol.pi.iter().enumerate() { + assert!( + (p - 1.0 / 3.0).abs() < 1e-14, + "node {i} got {p}, expected 1/3" + ); + } + } + + /// Issue #146: a pure forward star at 2:1. The hub beat both spokes, so it + /// must rank first — and by detailed balance on this tree, exactly 2:1:1. + #[test] + fn star_regression_issue_146() { + let chain = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (0, 2, 2.0, 1.0)]); + let sol = solve(&chain, opts()); + assert!((sol.pi[0] - 0.5).abs() < 1e-14, "hub: {}", sol.pi[0]); + assert!((sol.pi[1] - 0.25).abs() < 1e-14, "spoke a: {}", sol.pi[1]); + assert!((sol.pi[2] - 0.25).abs() < 1e-14, "spoke b: {}", sol.pi[2]); + } + + /// A single comparison at ratio r:1 must give exactly r:1. + #[test] + fn single_pair_is_the_ratio() { + for r in [1.0f64, 2.0, 3.0, 7.0] { + let chain = chain_from_votes(2, &[(0, 1, r, 1.0)]); + let sol = solve(&chain, opts()); + let want = r / (r + 1.0); + assert!( + (sol.pi[0] - want).abs() < 1e-15, + "ratio {r}: got {:?}", + sol.pi + ); + } + } + + // ----------------------------------------------------------------- + // Cross-method equivalence + // ----------------------------------------------------------------- + + #[test] + fn all_methods_agree_on_well_conditioned_graphs() { + let cases: Vec<(&str, RankChain)> = vec![ + ("star-64", star_graph(64)), + ("clique-40", clique_graph(40)), + ("sparse-d6-120", sparse_graph(120, 6)), + ("path-40", path_graph(40, 2.0)), + ]; + let tight = SolveOptions { + tol: 1e-14, + max_iters: 500_000, + ..opts() + }; + for (label, chain) in cases { + let reference = dense_gth(&chain, TOL).pi; + for sol in [ + force_sparse(&chain), + dense_lu(&chain, TOL), + sor(&chain, tight, 1.0), + power(&chain, tight), + bicgstab(&chain, tight), + solve(&chain, opts()), + ] { + // On a path the scores span 2^39, so an iterative method that + // has driven the *absolute* residual to 1e-14 is still far off + // in relative terms at the tail. Only the direct methods get + // full relative precision there; that gap is the whole point. + let rel = max_rel_error(&sol.pi, &reference); + let bound = if label == "path-40" && !sol.method.is_direct() { + 1e-2 + } else { + 1e-9 + }; + assert!( + rel < bound, + "{label}: {} disagrees with dense GTH by {rel:e}", + sol.method.label() + ); + } + } + } + + /// Every solver's own residual must match what an independent recomputation + /// says, so `Solution::residual` can be trusted as the acceptance criterion. + #[test] + fn reported_residual_matches_recomputation() { + let chain = sparse_graph(200, 6); + for sol in [ + solve(&chain, opts()), + power(&chain, opts()), + sor(&chain, opts(), 1.0), + dense_gth(&chain, TOL), + ] { + let recomputed = chain.residual(&sol.pi); + assert!( + (recomputed - sol.residual).abs() <= 1e-18 + recomputed * 1e-9, + "{}: reported {:e} vs recomputed {recomputed:e}", + sol.method.label(), + sol.residual + ); + } + } + + // ----------------------------------------------------------------- + // No silent non-convergence + // ----------------------------------------------------------------- + + #[test] + fn power_iteration_reports_its_own_failure() { + let chain = path_graph(1500, 2.0); + let capped = SolveOptions { + max_iters: 50, + ..opts() + }; + let sol = power(&chain, capped); + assert!(!sol.converged, "50 sweeps cannot solve a 1500-node path"); + assert_eq!(sol.iterations, 50, "should have used the whole budget"); + assert!(sol.residual > TOL); + } + + /// The same input the capped power iteration fails on must come back solved + /// through the hybrid, with the flag set honestly. + #[test] + fn hybrid_converges_where_power_iteration_cannot() { + for n in [1500usize, 4000] { + let chain = path_graph(n, 2.0); + let pow = power(&chain, opts()); + assert!( + !pow.converged, + "n={n}: power iteration was expected to hit the 10k cap" + ); + + let sol = solve(&chain, opts()); + assert!(sol.converged, "n={n}: hybrid failed"); + assert_eq!(sol.method, Method::SparseGth); + assert!(sol.residual < 1e-12, "n={n}: residual {:e}", sol.residual); + } + } + + /// On a path the true ranking is the construction order. `pi` runs out of + /// exponent range past ~1000 nodes, but `log_pi` must stay strictly ordered + /// all the way down — that is what makes it the right sort key. + #[test] + fn long_path_ranking_is_strictly_ordered_in_log_space() { + let n = 3000; + let chain = path_graph(n, 2.0); + let sol = solve(&chain, opts()); + assert!(sol.converged); + assert!( + sol.underflowed, + "a 3000-node path at 2:1 spans 10^903 and must report underflow" + ); + for i in 0..n - 1 { + assert!( + sol.log_pi[i] > sol.log_pi[i + 1], + "log scores tie or invert at {i}: {} vs {}", + sol.log_pi[i], + sol.log_pi[i + 1] + ); + } + // Every adjacent step is exactly ln 2 apart (detailed balance). + for i in 0..n - 1 { + let step = sol.log_pi[i] - sol.log_pi[i + 1]; + assert!( + (step - 2f64.ln()).abs() < 1e-9, + "step at {i} is {step}, expected ln 2" + ); + } + } + + // ----------------------------------------------------------------- + // Structure + // ----------------------------------------------------------------- + + /// A chain over disconnected components has infinitely many stationary + /// distributions. Power iteration from uniform picks the one where each + /// component keeps its share of the nodes; the solver must pick the same + /// one, or whole-group rankings would silently change meaning. + #[test] + fn disjoint_components_split_mass_by_size() { + let chain = chain_from_votes(4, &[(0, 1, 3.0, 1.0), (2, 3, 2.0, 1.0)]); + let sol = solve(&chain, opts()); + assert!(sol.pi.iter().all(|p| p.is_finite())); + // Each pair holds 1/2 of the mass, split 3:1 and 2:1 internally. + for (i, want) in [(0, 0.375), (1, 0.125), (2, 1.0 / 3.0), (3, 1.0 / 6.0)] { + assert!( + (sol.pi[i] - want).abs() < 1e-14, + "node {i}: got {} want {want}", + sol.pi[i] + ); + } + + // And that is what unbounded power iteration converges to. + let reference = power( + &chain, + SolveOptions { + tol: 1e-15, + max_iters: 200_000, + ..opts() + }, + ); + assert!(max_rel_error(&sol.pi, &reference.pi) < 1e-9); + } + + /// A lone item with no comparisons keeps a `1/n` share, exactly as the + /// uniform-start power iteration leaves it. + #[test] + fn isolated_nodes_keep_a_uniform_share() { + let chain = chain_from_votes(3, &[(0, 1, 3.0, 1.0)]); + let sol = solve(&chain, opts()); + assert!( + (sol.pi[2] - 1.0 / 3.0).abs() < 1e-14, + "isolate: {}", + sol.pi[2] + ); + assert!((sol.pi[0] - 0.5).abs() < 1e-14, "winner: {}", sol.pi[0]); + assert!( + (sol.pi[1] - 1.0 / 6.0).abs() < 1e-14, + "loser: {}", + sol.pi[1] + ); + } + + #[test] + fn cyclomatic_number_counts_independent_cycles() { + assert_eq!( + cyclomatic_number(&path_graph(50, 2.0)), + 0, + "a path is a tree" + ); + assert_eq!(cyclomatic_number(&star_graph(50)), 0, "a star is a tree"); + let triangle = chain_from_votes(3, &[(0, 1, 2.0, 1.0), (1, 2, 2.0, 1.0), (2, 0, 2.0, 1.0)]); + assert_eq!(cyclomatic_number(&triangle), 1); + // K_n has n(n-1)/2 edges and n-1 spanning-tree edges. + assert_eq!(cyclomatic_number(&clique_graph(6)), 15 - 5); + } + + #[test] + fn hybrid_picks_direct_for_trees_and_iterative_for_dense_graphs() { + assert_eq!( + solve(&path_graph(4000, 2.0), opts()).method, + Method::SparseGth + ); + assert_eq!(solve(&clique_graph(60), opts()).method, Method::SparseGth); + let big_sparse = solve(&sparse_graph(2000, 6), opts()); + assert_ne!(big_sparse.method, Method::SparseGth); + assert!(big_sparse.converged); + } + + #[test] + fn degenerate_inputs_do_not_panic() { + let empty = RankChain::from_normalized(0, vec![], 0); + assert!(solve(&empty, opts()).pi.is_empty()); + + let single = RankChain::from_normalized(1, vec![], 0); + assert_eq!(solve(&single, opts()).pi, vec![1.0]); + + // Nodes with no comparisons at all: uniform, and no division by zero. + let isolated = RankChain::from_normalized(5, vec![], 0); + let sol = solve(&isolated, opts()); + assert!(sol.pi.iter().all(|&p| (p - 0.2).abs() < 1e-15)); + + // A unanimous edge (the loser never scored) leaves one direction empty. + let unanimous = chain_from_votes(2, &[(0, 1, 1.0, 0.0)]); + let sol = solve(&unanimous, opts()); + assert!(sol.pi[0] > sol.pi[1], "{:?}", sol.pi); + assert!(sol.pi.iter().all(|p| p.is_finite())); + } + + // ----------------------------------------------------------------- + // Determinism + // ----------------------------------------------------------------- + + /// The same graph presented in a different edge order must give a + /// bit-identical answer — `RankChain::from_normalized` sorts precisely so + /// that `HashMap` iteration order cannot leak into the low bits. + #[test] + fn edge_input_order_does_not_change_the_result() { + let n = 300; + let base: Vec<((usize, usize), f64)> = { + let mut v = Vec::new(); + for i in 0..n - 1 { + v.push(((i + 1, i), 2.0 / 3.0)); + v.push(((i, i + 1), 1.0 / 3.0)); + } + for i in 0..n / 3 { + v.push(((i, (7 * i + 11) % n), 0.5)); + v.push((((7 * i + 11) % n, i), 0.5)); + } + v + }; + let reference = solve(&RankChain::from_normalized(n, base.clone(), 6), opts()).pi; + + let mut state = 12345u64; + for _ in 0..6 { + let mut shuffled = base.clone(); + for i in (1..shuffled.len()).rev() { + state ^= state >> 12; + state ^= state << 25; + state ^= state >> 27; + let j = (state.wrapping_mul(0x2545_F491_4F6C_DD1D) % (i as u64 + 1)) as usize; + shuffled.swap(i, j); + } + let got = solve(&RankChain::from_normalized(n, shuffled, 6), opts()).pi; + assert_eq!(got, reference, "shuffled edge order changed the scores"); + } + } +}