fix(quality): support any number of rating groups
`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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DnsaJg74eNSva3PJjK2eej
This commit is contained in:
Claude Opus 5
+9
-1
@@ -389,7 +389,15 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
|
||||
|
||||
/// Draw-probability quality metric for the given teams (key slices).
|
||||
///
|
||||
/// Values range roughly [0, 1]; 1 == perfectly matched.
|
||||
/// Values range roughly [0, 1]; 1 == perfectly matched. Supports any
|
||||
/// number of teams.
|
||||
///
|
||||
/// # Panics
|
||||
///
|
||||
/// Panics if fewer than two teams are supplied, or if a team resolves to
|
||||
/// no known competitors — keys absent from the history, or competitors
|
||||
/// with no recorded skill, are dropped, so a team of entirely-unknown
|
||||
/// keys becomes empty. Use `lookup` to check keys first.
|
||||
pub fn predict_quality(&self, teams: &[&[&K]]) -> f64 {
|
||||
let groups: Vec<Vec<Gaussian>> = teams
|
||||
.iter()
|
||||
|
||||
+23
-6
@@ -247,7 +247,26 @@ pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
|
||||
}
|
||||
|
||||
/// 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();
|
||||
|
||||
+316
-119
@@ -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::<Vec<_>>();
|
||||
|
||||
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<f64>,
|
||||
perm: Vec<usize>,
|
||||
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<usize> = (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)]);
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user