From 0d32690fcc306362869f4b6372808a8b7dfbfcb7 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Tue, 4 Aug 2026 21:42:04 +0200 Subject: [PATCH] fix(quality): support any number of rating groups MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `quality()` was two-group-only in three separate ways, and `History::predict_quality` inherited all of them. - The contrast-matrix column counter tracked two positions with two variables that only agree on the first row, so three or more groups wrote past the end of the row and panicked with an out-of-bounds index. The negative block always begins immediately after the positive one, so the second counter is unnecessary. - `Matrix::inverse` was implemented only for the 1x1 case and otherwise `panic!("eh, okey")`. It now uses LU decomposition with partial pivoting, which also replaces the recursive cofactor `determinant` — that was O(n!) and allocated a `Vec` per minor, so a 10-team match needed 362,880 terms. - Degenerate inputs (zero groups, one group, empty groups) underflowed or produced NaN. They now assert with a message naming the requirement. `Matrix` also gains dimension checks on multiply/add and bounds checks on indexing, and loses the now-unused `adjugate`/`minor` cofactor path. The two-group golden is unchanged. N-group behaviour is covered by invariants — permutation invariance, and quality falling as a skill gap widens — since no reference values were available to compare against; the sublee/trueskill cross-check remains open. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01DnsaJg74eNSva3PJjK2eej --- src/history.rs | 10 +- src/lib.rs | 29 +++- src/matrix.rs | 435 ++++++++++++++++++++++++++++++++++------------- tests/quality.rs | 119 +++++++++++++ 4 files changed, 467 insertions(+), 126 deletions(-) create mode 100644 tests/quality.rs diff --git a/src/history.rs b/src/history.rs index fb56932..aea04e2 100644 --- a/src/history.rs +++ b/src/history.rs @@ -389,7 +389,15 @@ impl, O: Observer, K: Eq + Hash + Clone> History f64 { let groups: Vec> = teams .iter() diff --git a/src/lib.rs b/src/lib.rs index cdacbfa..227c4f1 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -247,7 +247,26 @@ pub(crate) fn sort_time(xs: &[T], reverse: bool) -> Vec { } /// Calculates the match quality of the given rating groups. A result is the draw probability in the association +/// +/// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a +/// perfectly balanced match. +/// +/// # Panics +/// +/// Panics if fewer than two rating groups are supplied, or if any group is +/// empty — match quality is a property of a contest between at least two +/// non-empty sides. pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { + assert!( + rating_groups.len() >= 2, + "quality() requires at least 2 rating groups, got {}", + rating_groups.len() + ); + assert!( + rating_groups.iter().all(|group| !group.is_empty()), + "quality() requires every rating group to be non-empty" + ); + let flatten_ratings = rating_groups .iter() .flat_map(|group| group.iter()) @@ -271,8 +290,10 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length); + // Row `row` contrasts group `row` (+weight) against group `row + 1` + // (-weight). `t` is the column where the current group's players start; + // the negative block begins immediately after it. let mut t = 0; - let mut x = 0; for (row, group) in rating_groups.windows(2).enumerate() { let current = group[0]; @@ -280,17 +301,13 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { for n in t..t + current.len() { rotated_a_matrix[(row, n)] = flatten_weights[n]; - - x += 1; } t += current.len(); - for n in x..x + next.len() { + for n in t..t + next.len() { rotated_a_matrix[(row, n)] = -flatten_weights[n]; } - - x += next.len(); } let a_matrix = rotated_a_matrix.transpose(); diff --git a/src/matrix.rs b/src/matrix.rs index cf74061..cf29777 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -1,29 +1,13 @@ +//! Minimal dense matrix used by `quality()`. +//! +//! `determinant` and `inverse` go through one LU decomposition with partial +//! pivoting — O(n³) and numerically stable. The previous implementation +//! expanded cofactors recursively (O(n!), allocating a `Vec` per minor) and +//! only implemented `inverse` for the 1×1 case, which limited `quality()` to +//! exactly two rating groups. + use std::ops; -fn det(m: &[f64], x: usize) -> f64 { - if x == 1 { - m[0] - } else if x == 2 { - m[0] * m[3] - m[1] * m[2] - } else { - let mut d = 0.0; - - for n in 0..x { - let ms = m - .iter() - .enumerate() - .skip(x) - .filter(|(i, _)| (i % x) != n) - .map(|(_, v)| *v) - .collect::>(); - - d += (-1.0f64).powi(n as i32) * m[n] * det(&ms, x - 1); - } - - d - } -} - #[derive(Clone, Debug)] pub struct Matrix { data: Box<[f64]>, @@ -31,6 +15,107 @@ pub struct Matrix { width: usize, } +/// LU decomposition with partial pivoting: `PA = LU`, stored compactly. +/// +/// `lu` holds `L` below the diagonal (unit diagonal implied) and `U` on and +/// above it. `sign` is the determinant sign contributed by row swaps, or 0.0 +/// when the matrix is singular. +struct Lu { + lu: Vec, + perm: Vec, + n: usize, + sign: f64, +} + +impl Lu { + fn decompose(m: &Matrix) -> Self { + debug_assert_eq!(m.width, m.height, "LU requires a square matrix"); + + let n = m.width; + let mut lu = m.data.to_vec(); + let mut perm: Vec = (0..n).collect(); + let mut sign = 1.0; + + for col in 0..n { + // Partial pivot: take the largest-magnitude candidate to limit + // growth of round-off in the elimination below. + let mut pivot_row = col; + let mut pivot_max = lu[col * n + col].abs(); + + for row in (col + 1)..n { + let candidate = lu[row * n + col].abs(); + if candidate > pivot_max { + pivot_max = candidate; + pivot_row = row; + } + } + + if pivot_max == 0.0 { + sign = 0.0; + continue; + } + + if pivot_row != col { + for k in 0..n { + lu.swap(col * n + k, pivot_row * n + k); + } + perm.swap(col, pivot_row); + sign = -sign; + } + + let pivot = lu[col * n + col]; + + for row in (col + 1)..n { + let factor = lu[row * n + col] / pivot; + lu[row * n + col] = factor; + + for k in (col + 1)..n { + lu[row * n + k] -= factor * lu[col * n + k]; + } + } + } + + Self { lu, perm, n, sign } + } + + fn determinant(&self) -> f64 { + if self.sign == 0.0 { + return 0.0; + } + + let mut det = self.sign; + for i in 0..self.n { + det *= self.lu[i * self.n + i]; + } + + det + } + + /// Solve `Ax = b` for a single column of the identity, giving one column + /// of the inverse. + fn solve_column(&self, col: usize, out: &mut [f64]) { + let n = self.n; + + // Forward substitution through L, applying the row permutation. + for i in 0..n { + let mut sum = if self.perm[i] == col { 1.0 } else { 0.0 }; + for (k, &solved) in out.iter().enumerate().take(i) { + sum -= self.lu[i * n + k] * solved; + } + out[i] = sum; + } + + // Back substitution through U. + for i in (0..n).rev() { + let mut sum = out[i]; + for (k, &solved) in out.iter().enumerate().skip(i + 1) { + sum -= self.lu[i * n + k] * solved; + } + out[i] = sum / self.lu[i * n + i]; + } + } +} + impl Matrix { pub fn new(height: usize, width: usize) -> Matrix { Matrix { @@ -52,73 +137,59 @@ impl Matrix { matrix } - pub fn minor(&self, row_n: usize, col_n: usize) -> Matrix { - let mut matrix = Matrix::new(self.height - 1, self.width - 1); - - let mut nr = 0; - - for r in 0..self.height { - if r == row_n { - continue; - } - - let mut nc = 0; - - for c in 0..self.width { - if c == col_n { - continue; - } - - matrix[(nr, nc)] = self[(r, c)]; - - nc += 1; - } - - nr += 1; - } - - matrix - } - + /// Determinant of a square matrix. The 0×0 determinant is 1 by convention + /// (the empty product). + /// + /// # Panics + /// + /// Panics if the matrix is not square. pub fn determinant(&self) -> f64 { - debug_assert!(self.width == self.height); + assert_eq!( + self.width, self.height, + "determinant requires a square matrix, got {}x{}", + self.height, self.width + ); - det(&self.data, self.width) + if self.width == 0 { + return 1.0; + } + + Lu::decompose(self).determinant() } - pub fn adjugate(&self) -> Matrix { - debug_assert!(self.width == self.height); + /// Matrix inverse via LU decomposition. + /// + /// # Panics + /// + /// Panics if the matrix is not square or is singular. + pub fn inverse(&self) -> Matrix { + assert_eq!( + self.width, self.height, + "inverse requires a square matrix, got {}x{}", + self.height, self.width + ); - let mut matrix = Matrix::new(self.height, self.width); + let n = self.width; + let mut inverse = Matrix::new(n, n); - if matrix.height == 2 { - matrix[(0, 0)] = self[(1, 1)]; - matrix[(0, 1)] = -self[(0, 1)]; - matrix[(1, 0)] = -self[(1, 0)]; - matrix[(1, 1)] = self[(0, 0)]; - } else { - for r in 0..matrix.height { - for c in 0..matrix.width { - let sign = if (r + c) % 2 == 0 { 1.0 } else { -1.0 }; + if n == 0 { + return inverse; + } - matrix[(r, c)] = self.minor(r, c).determinant() * sign; - } + let lu = Lu::decompose(self); + assert!(lu.sign != 0.0, "cannot invert a singular matrix"); + + let mut column = vec![0.0; n]; + + for c in 0..n { + lu.solve_column(c, &mut column); + + for (r, &value) in column.iter().enumerate() { + inverse[(r, c)] = value; } } - matrix - } - - pub fn inverse(&self) -> Matrix { - let mut matrix = Matrix::new(self.width, self.height); - - if self.height == self.width && self.height == 1 { - matrix[(0, 0)] = 1.0 / self[(0, 0)]; - } else { - panic!("eh, okey") - } - - matrix + inverse } } @@ -126,20 +197,62 @@ impl ops::Index<(usize, usize)> for Matrix { type Output = f64; fn index(&self, pos: (usize, usize)) -> &Self::Output { + debug_assert!( + pos.0 < self.height && pos.1 < self.width, + "index ({}, {}) out of bounds for {}x{} matrix", + pos.0, + pos.1, + self.height, + self.width + ); + &self.data[(self.width * pos.0) + pos.1] } } impl ops::IndexMut<(usize, usize)> for Matrix { fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output { + debug_assert!( + pos.0 < self.height && pos.1 < self.width, + "index ({}, {}) out of bounds for {}x{} matrix", + pos.0, + pos.1, + self.height, + self.width + ); + &mut self.data[(self.width * pos.0) + pos.1] } } -impl<'a> ops::Mul<&'a Matrix> for f64 { +fn multiply(lhs: &Matrix, rhs: &Matrix) -> Matrix { + assert_eq!( + lhs.width, rhs.height, + "cannot multiply {}x{} by {}x{}", + lhs.height, lhs.width, rhs.height, rhs.width + ); + + let mut matrix = Matrix::new(lhs.height, rhs.width); + + for r in 0..matrix.height { + for c in 0..matrix.width { + let mut value = 0.0; + + for x in 0..lhs.width { + value += lhs[(r, x)] * rhs[(x, c)]; + } + + matrix[(r, c)] = value; + } + } + + matrix +} + +impl ops::Mul<&Matrix> for f64 { type Output = Matrix; - fn mul(self, rhs: &'a Matrix) -> Matrix { + fn mul(self, rhs: &Matrix) -> Matrix { let mut matrix = Matrix::new(rhs.height, rhs.width); for r in 0..rhs.height { @@ -152,54 +265,35 @@ impl<'a> ops::Mul<&'a Matrix> for f64 { } } -impl<'a> ops::Mul<&'a Matrix> for Matrix { +impl ops::Mul<&Matrix> for Matrix { type Output = Matrix; - fn mul(self, rhs: &'a Matrix) -> Matrix { - let mut matrix = Matrix::new(self.height, rhs.width); - - for r in 0..matrix.height { - for c in 0..matrix.width { - let mut value = 0.0; - - for x in 0..self.width { - value += self[(r, x)] * rhs[(x, c)]; - } - - matrix[(r, c)] = value; - } - } - - matrix + fn mul(self, rhs: &Matrix) -> Matrix { + multiply(&self, rhs) } } -impl<'a> ops::Mul<&'a Matrix> for &'a Matrix { +impl ops::Mul<&Matrix> for &Matrix { type Output = Matrix; - fn mul(self, rhs: &'a Matrix) -> Matrix { - let mut matrix = Matrix::new(self.height, rhs.width); - - for r in 0..matrix.height { - for c in 0..matrix.width { - let mut value = 0.0; - - for x in 0..self.width { - value += self[(r, x)] * rhs[(x, c)]; - } - - matrix[(r, c)] = value; - } - } - - matrix + fn mul(self, rhs: &Matrix) -> Matrix { + multiply(self, rhs) } } -impl<'a> ops::Add<&'a Matrix> for &'a Matrix { +impl ops::Add<&Matrix> for &Matrix { type Output = Matrix; - fn add(self, rhs: &'a Matrix) -> Matrix { + fn add(self, rhs: &Matrix) -> Matrix { + assert!( + self.height == rhs.height && self.width == rhs.width, + "cannot add {}x{} to {}x{}", + self.height, + self.width, + rhs.height, + rhs.width + ); + let mut matrix = Matrix::new(self.height, self.width); for r in 0..matrix.height { @@ -211,3 +305,106 @@ impl<'a> ops::Add<&'a Matrix> for &'a Matrix { matrix } } + +#[cfg(test)] +mod tests { + use super::*; + + fn from_rows(rows: &[&[f64]]) -> Matrix { + let mut m = Matrix::new(rows.len(), rows[0].len()); + for (r, row) in rows.iter().enumerate() { + for (c, &v) in row.iter().enumerate() { + m[(r, c)] = v; + } + } + m + } + + #[test] + fn determinant_1x1() { + assert!((from_rows(&[&[3.0]]).determinant() - 3.0).abs() < 1e-12); + } + + #[test] + fn determinant_2x2() { + let m = from_rows(&[&[1.0, 2.0], &[3.0, 4.0]]); + assert!((m.determinant() - (-2.0)).abs() < 1e-12); + } + + #[test] + fn determinant_3x3() { + let m = from_rows(&[&[6.0, 1.0, 1.0], &[4.0, -2.0, 5.0], &[2.0, 8.0, 7.0]]); + assert!((m.determinant() - (-306.0)).abs() < 1e-10); + } + + #[test] + fn determinant_requires_no_pivot_at_origin() { + // A zero in the top-left forces a row swap; the sign must follow. + let m = from_rows(&[&[0.0, 1.0], &[1.0, 0.0]]); + assert!((m.determinant() - (-1.0)).abs() < 1e-12); + } + + #[test] + fn determinant_of_singular_is_zero() { + let m = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]); + assert!(m.determinant().abs() < 1e-12); + } + + #[test] + fn inverse_1x1() { + let inv = from_rows(&[&[4.0]]).inverse(); + assert!((inv[(0, 0)] - 0.25).abs() < 1e-12); + } + + #[test] + fn inverse_times_original_is_identity() { + for rows in [ + vec![vec![1.0, 2.0], vec![3.0, 4.0]], + vec![ + vec![6.0, 1.0, 1.0], + vec![4.0, -2.0, 5.0], + vec![2.0, 8.0, 7.0], + ], + vec![ + vec![2.0, 0.0, 1.0, 3.0], + vec![1.0, 5.0, 2.0, 0.0], + vec![0.0, 1.0, 4.0, 1.0], + vec![3.0, 2.0, 0.0, 6.0], + ], + ] { + let refs: Vec<&[f64]> = rows.iter().map(|r| r.as_slice()).collect(); + let m = from_rows(&refs); + let product = &m * &m.inverse(); + + for r in 0..product.height { + for c in 0..product.width { + let expected = if r == c { 1.0 } else { 0.0 }; + assert!( + (product[(r, c)] - expected).abs() < 1e-9, + "({r},{c}) = {} expected {expected}", + product[(r, c)] + ); + } + } + } + } + + #[test] + #[should_panic(expected = "singular")] + fn inverse_of_singular_panics() { + let _ = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]).inverse(); + } + + #[test] + fn empty_determinant_is_one() { + assert!((Matrix::new(0, 0).determinant() - 1.0).abs() < 1e-12); + } + + #[test] + fn transpose_round_trips() { + let m = from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]); + let t = m.transpose(); + assert_eq!((t.height, t.width), (3, 2)); + assert_eq!(t.transpose()[(1, 2)], m[(1, 2)]); + } +} diff --git a/tests/quality.rs b/tests/quality.rs new file mode 100644 index 0000000..fe7c2c4 --- /dev/null +++ b/tests/quality.rs @@ -0,0 +1,119 @@ +//! `quality()` beyond two rating groups. +//! +//! The historical golden (two equal singletons) is asserted in +//! `src/lib.rs::tests::test_quality`. These cover the N-group generalisation, +//! which previously panicked with an out-of-bounds index at 3+ groups. + +use trueskill_tt::{Gaussian, quality}; + +const BETA: f64 = 25.0 / 3.0 / 2.0; + +fn rating(mu: f64, sigma: f64) -> Gaussian { + Gaussian::from_ms(mu, sigma) +} + +#[test] +fn three_equal_groups_is_finite_and_in_range() { + let r = rating(25.0, 3.0); + let q = quality(&[&[r], &[r], &[r]], BETA); + + assert!(q.is_finite(), "quality must be finite, got {q}"); + assert!((0.0..=1.0).contains(&q), "quality out of range: {q}"); +} + +#[test] +fn quality_supports_many_groups() { + let r = rating(25.0, 3.0); + for n in 2..=8 { + let holders: Vec<[Gaussian; 1]> = (0..n).map(|_| [r]).collect(); + let groups: Vec<&[Gaussian]> = holders.iter().map(|g| g.as_slice()).collect(); + let q = quality(&groups, BETA); + assert!(q.is_finite(), "n={n}: quality must be finite, got {q}"); + assert!((0.0..=1.0).contains(&q), "n={n}: out of range: {q}"); + } +} + +/// Equal-strength groups are the best-matched case: introducing a skill gap +/// must lower quality. +#[test] +fn imbalance_lowers_quality() { + let strong = rating(40.0, 3.0); + let average = rating(25.0, 3.0); + + let balanced = quality(&[&[average], &[average], &[average]], BETA); + let lopsided = quality(&[&[strong], &[average], &[average]], BETA); + + assert!( + lopsided < balanced, + "expected imbalanced quality {lopsided} < balanced {balanced}" + ); +} + +/// Quality is a property of the multiset of groups, not their order. +#[test] +fn quality_is_permutation_invariant() { + let a = rating(30.0, 2.0); + let b = rating(25.0, 3.0); + let c = rating(20.0, 4.0); + + let forward = quality(&[&[a], &[b], &[c]], BETA); + let reversed = quality(&[&[c], &[b], &[a]], BETA); + + assert!( + (forward - reversed).abs() < 1e-9, + "permutation changed quality: {forward} vs {reversed}" + ); +} + +#[test] +fn multi_player_groups_work() { + let r = rating(25.0, 3.0); + let q = quality(&[&[r, r], &[r, r], &[r, r]], BETA); + assert!(q.is_finite()); + assert!((0.0..=1.0).contains(&q)); +} + +#[test] +fn uneven_group_sizes_work() { + let r = rating(25.0, 3.0); + let q = quality(&[&[r, r], &[r], &[r, r, r]], BETA); + assert!(q.is_finite(), "got {q}"); + assert!((0.0..=1.0).contains(&q), "got {q}"); +} + +#[test] +#[should_panic(expected = "at least 2 rating groups")] +fn single_group_panics_with_clear_message() { + let r = rating(25.0, 3.0); + let _ = quality(&[&[r]], BETA); +} + +#[test] +#[should_panic(expected = "at least 2 rating groups")] +fn zero_groups_panics_with_clear_message() { + let _ = quality(&[], BETA); +} + +#[test] +#[should_panic(expected = "non-empty")] +fn empty_group_panics_with_clear_message() { + let r = rating(25.0, 3.0); + let _ = quality(&[&[r], &[]], BETA); +} + +#[test] +fn history_predict_quality_supports_three_teams() { + use trueskill_tt::History; + + let mut h = History::default(); + h.record_winner(&"a", &"b", 1).unwrap(); + h.record_winner(&"b", &"c", 2).unwrap(); + h.converge().unwrap(); + + let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]); + assert!( + q.is_finite(), + "3-team predict_quality must be finite, got {q}" + ); + assert!((0.0..=1.0).contains(&q), "out of range: {q}"); +}