diff --git a/src/trueskill/factor_graph.rs b/src/trueskill/factor_graph.rs new file mode 100644 index 0000000..adc3ca3 --- /dev/null +++ b/src/trueskill/factor_graph.rs @@ -0,0 +1,276 @@ +use std::cell::RefCell; +use std::collections::HashMap; +use std::rc::Rc; + +use super::gaussian::Gaussian; + +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Variable { + pub gaussian: Gaussian, + messages: HashMap, +} + +impl Variable { + pub fn new() -> Self { + Self { + gaussian: Gaussian::default(), + messages: HashMap::new(), + } + } + + fn set(&mut self, val: Gaussian) -> f64 { + let delta = self.delta(val); + + self.gaussian.pi = val.pi; + self.gaussian.tau = val.tau; + + delta + } + + fn update_message(&mut self, factor_id: u32, message: Gaussian) -> f64 { + let old_message = self.messages[&factor_id]; + let v = self.messages.entry(factor_id).or_default(); + *v = message; + + self.set(self.gaussian / old_message * message) + } + + fn update_value(&mut self, factor_id: u32, val: Gaussian) -> f64 { + let old_message = self.messages[&factor_id]; + let v = self.messages.entry(factor_id).or_default(); + *v = val * old_message / self.gaussian; + + self.set(val) + } + + fn delta(&self, other: Gaussian) -> f64 { + let pi_delta = (self.gaussian.pi - other.pi).abs(); + if pi_delta.is_infinite() { + return 0.0; + } + + (self.gaussian.tau - other.tau).abs().max(pi_delta.sqrt()) + } +} + +pub struct PriorFactor { + id: u32, + pub variable: Rc>, + val: Gaussian, + dynamic: f64, +} + +impl PriorFactor { + pub fn new(id: u32, variable: Rc>, val: Gaussian, dynamic: f64) -> Self { + variable.borrow_mut().messages.entry(id).or_default(); + + Self { + id, + variable, + val, + dynamic, + } + } + + pub fn down(&mut self) -> f64 { + let sigma = self.val.sigma().hypot(self.dynamic); + let value = Gaussian::with_mu_sigma(self.val.mu(), sigma); + self.variable.borrow_mut().update_value(self.id, value) + } +} + +pub struct LikelihoodFactor { + id: u32, + mean: Rc>, + value: Rc>, + variance: f64, +} + +impl LikelihoodFactor { + pub fn new( + id: u32, + mean: Rc>, + value: Rc>, + variance: f64, + ) -> Self { + mean.borrow_mut().messages.entry(id).or_default(); + value.borrow_mut().messages.entry(id).or_default(); + + Self { + id, + mean, + value, + variance, + } + } + + pub fn down(&mut self) -> f64 { + let msg = { + let mean = self.mean.borrow(); + mean.gaussian / mean.messages[&self.id] + }; + let a = self.calc_a(msg); + self.value + .borrow_mut() + .update_message(self.id, Gaussian::with_pi_tau(a * msg.pi, a * msg.tau)) + } + + pub fn up(&mut self) -> f64 { + let msg = { + let value = self.value.borrow(); + value.gaussian / value.messages[&self.id] + }; + let a = self.calc_a(msg); + self.mean + .borrow_mut() + .update_message(self.id, Gaussian::with_pi_tau(a * msg.pi, a * msg.tau)) + } + + fn calc_a(&self, gaussian: Gaussian) -> f64 { + self.variance.mul_add(gaussian.pi, 1.0).recip() + } +} + +pub struct SumFactor { + id: u32, + sum: Rc>, + terms: Vec>>, + coeffs: Vec, +} + +impl SumFactor { + pub fn new( + id: u32, + sum: Rc>, + terms: Vec>>, + coeffs: Vec, + ) -> Self { + sum.borrow_mut().messages.entry(id).or_default(); + for term in &terms { + term.borrow_mut().messages.entry(id).or_default(); + } + + Self { + id, + sum, + terms, + coeffs, + } + } + + pub fn down(&mut self) -> f64 { + let msgs: Vec = self + .terms + .iter() + .map(|term| term.borrow().messages[&self.id]) + .collect(); + self.update(&self.sum, &self.terms, &msgs, &self.coeffs) + } + + pub fn up(&mut self, index: usize) -> f64 { + let coeff = self.coeffs[index]; + let mut coeffs = Vec::new(); + for (x, c) in self.coeffs.iter().enumerate() { + if coeff == 0.0 { + coeffs.push(0.0); + } else if x == index { + coeffs.push(coeff.recip()); + } else { + coeffs.push(-(*c) / coeff); + } + } + + let mut vals = self.terms.clone(); + vals[index] = self.sum.clone(); + let msgs: Vec = vals + .iter() + .map(|val| val.borrow().messages[&self.id]) + .collect(); + + self.update(&self.terms[index], &vals, &msgs, &coeffs) + } + + #[inline] + pub fn terms_len(&self) -> usize { + self.terms.len() + } + + fn update( + &self, + var: &Rc>, + vals: &[Rc>], + msgs: &[Gaussian], + coeffs: &[f64], + ) -> f64 { + let mut pi_inv = 0.0_f64; + let mut mu = 0.0; + + for ((val, msg), coeff) in vals.iter().zip(msgs).zip(coeffs) { + let div = val.borrow().gaussian / *msg; + mu += coeff * div.mu(); + if pi_inv.is_infinite() { + continue; + } + + if div.pi == 0.0 { + pi_inv = f64::INFINITY; + } else { + pi_inv += coeff.powi(2) / div.pi; + } + } + + let pi = pi_inv.recip(); + let tau = pi * mu; + + var.borrow_mut() + .update_message(self.id, Gaussian::with_pi_tau(pi, tau)) + } +} + +pub struct TruncateFactor { + id: u32, + variable: Rc>, + v_func: Box f64>, + w_func: Box f64>, + draw_margin: f64, +} + +impl TruncateFactor { + pub fn new( + id: u32, + variable: Rc>, + v_func: Box f64>, + w_func: Box f64>, + draw_margin: f64, + ) -> Self { + variable.borrow_mut().messages.entry(id).or_default(); + + Self { + id, + variable, + v_func, + w_func, + draw_margin, + } + } + + pub fn up(&mut self) -> f64 { + let div = { + let variable = self.variable.borrow(); + variable.gaussian / variable.messages[&self.id] + }; + let pi_sqrt = div.pi.sqrt(); + let arg_1 = div.tau / pi_sqrt; + let arg_2 = self.draw_margin * pi_sqrt; + let v = (self.v_func)(arg_1, arg_2); + let w = (self.w_func)(arg_1, arg_2); + let denom = 1.0 - w; + + let pi = div.pi / denom; + let tau = pi_sqrt.mul_add(v, div.tau) / denom; + + self.variable + .borrow_mut() + .update_value(self.id, Gaussian::with_pi_tau(pi, tau)) + } +} diff --git a/src/trueskill/gaussian.rs b/src/trueskill/gaussian.rs new file mode 100644 index 0000000..8bfb7e5 --- /dev/null +++ b/src/trueskill/gaussian.rs @@ -0,0 +1,65 @@ +use std::cmp::Ordering; +use std::ops::{Div, Mul}; + +#[derive(Copy, Clone, Debug, Default, PartialEq)] +pub struct Gaussian { + pub pi: f64, + pub tau: f64, +} + +impl Gaussian { + pub fn with_mu_sigma(mu: f64, sigma: f64) -> Self { + assert_ne!(sigma, 0.0, "sigma^2 needs to be greater than 0"); + + let pi = sigma.powi(-2); + Self { pi, tau: pi * mu } + } + + pub const fn with_pi_tau(pi: f64, tau: f64) -> Self { + Self { pi, tau } + } + + pub fn mu(&self) -> f64 { + if self.pi == 0.0 { + return 0.0; + } + + self.tau / self.pi + } + + pub fn sigma(&self) -> f64 { + if self.pi == 0.0 { + return f64::INFINITY; + } + + self.pi.recip().sqrt() + } +} + +impl Mul for Gaussian { + type Output = Self; + + fn mul(self, rhs: Self) -> Self::Output { + Self { + pi: self.pi + rhs.pi, + tau: self.tau + rhs.tau, + } + } +} + +impl Div for Gaussian { + type Output = Self; + + fn div(self, rhs: Self) -> Self::Output { + Self { + pi: self.pi - rhs.pi, + tau: self.tau - rhs.tau, + } + } +} + +impl PartialOrd for Gaussian { + fn partial_cmp(&self, other: &Self) -> Option { + self.mu().partial_cmp(&other.mu()) + } +} diff --git a/src/trueskill/matrix.rs b/src/trueskill/matrix.rs new file mode 100644 index 0000000..d6ef084 --- /dev/null +++ b/src/trueskill/matrix.rs @@ -0,0 +1,238 @@ +use super::TrueSkillRating; + +// This Matrix could have been imported, but we implement it ourselves, since we only have to use some basic things here. +#[derive(Clone, Debug)] +pub struct Matrix { + data: Vec, + rows: usize, + cols: usize, +} + +impl Matrix { + pub fn set(&mut self, row: usize, col: usize, val: f64) { + self.data[row * self.cols + col] = val; + } + + pub fn get(&self, row: usize, col: usize) -> f64 { + self.data[row * self.cols + col] + } + + pub fn new(rows: usize, cols: usize) -> Self { + Self { + data: vec![0.0; rows * cols], + rows, + cols, + } + } + + pub fn new_from_data(data: &[f64], rows: usize, cols: usize) -> Self { + Self { + data: data.to_vec(), + rows, + cols, + } + } + + pub fn new_diagonal(data: &[f64]) -> Self { + let mut matrix = Self::new(data.len(), data.len()); + + for (i, val) in data.iter().enumerate() { + matrix.set(i, i, *val); + } + + matrix + } + + pub fn create_rotated_a_matrix(teams: &[&[TrueSkillRating]]) -> Self { + let total_players = teams.iter().map(|team| team.len()).sum::(); + + let mut player_assignments: Vec = vec![]; + + let mut total_previous_players = 0; + + let team_assignments_list_count = teams.len(); + + for current_column in 0..team_assignments_list_count - 1 { + let current_team = teams[current_column]; + + player_assignments.append(&mut vec![0.0; total_previous_players]); + + for _current_player in current_team { + player_assignments.push(1.0); // TODO: Replace 1.0 by partial play weighting + total_previous_players += 1; + } + + let mut rows_remaining = total_players - total_previous_players; + let next_team = teams[current_column + 1]; + + for _next_player in next_team { + player_assignments.push(-1.0 * 1.0); // TODO: Replace 1.0 by partial play weighting + rows_remaining -= 1; + } + + player_assignments.append(&mut vec![0.0; rows_remaining]); + } + + Self::new_from_data( + &player_assignments, + team_assignments_list_count - 1, + total_players, + ) + } + + pub fn transpose(&self) -> Self { + let mut matrix = Self::new(self.cols, self.rows); + + for i in 0..self.rows { + for j in 0..self.cols { + matrix.set(j, i, self.get(i, j)); + } + } + + matrix + } + + #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)] + pub fn determinant(&self) -> f64 { + assert_eq!(self.rows, self.cols, "Matrix must be square"); + + if self.rows == 1 { + return self.get(0, 0); + } + + let mut sum = 0.0; + + for i in 0..self.rows { + sum += self.get(0, i) * self.minor(0, i).determinant() * (-1.0_f64).powi(i as i32); + } + + sum + } + + pub fn minor(&self, row: usize, col: usize) -> Self { + let mut matrix = Self::new(self.rows - 1, self.cols - 1); + + for i in 0..self.rows { + for j in 0..self.cols { + if i != row && j != col { + matrix.set( + if i > row { i - 1 } else { i }, + if j > col { j - 1 } else { j }, + self.get(i, j), + ); + } + } + } + + matrix + } + + #[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)] + pub fn adjugate(&self) -> Self { + let mut matrix = Self::new(self.rows, self.cols); + + for i in 0..self.rows { + for j in 0..self.cols { + matrix.set( + i, + j, + self.minor(j, i).determinant() * (-1.0_f64).powi((i + j) as i32), + ); + } + } + + matrix + } + + pub fn inverse(&self) -> Self { + let det = self.determinant(); + + // Avoiding 1/0 + assert!((det != 0.0), "Matrix is not invertible"); + + self.adjugate() * det.recip() + } +} + +impl std::ops::Mul for Matrix { + type Output = Self; + + fn mul(self, rhs: Self) -> Self::Output { + if self.cols == rhs.rows { + let mut matrix = Self::new(self.rows, rhs.cols); + + for i in 0..self.rows { + for j in 0..rhs.cols { + let mut sum = 0.0; + + for k in 0..self.cols { + sum += self.get(i, k) * rhs.get(k, j); + } + + matrix.set(i, j, sum); + } + } + + matrix + } else if self.rows == rhs.cols { + let mut matrix = Self::new(self.cols, rhs.rows); + + for i in 0..self.cols { + for j in 0..rhs.rows { + let mut sum = 0.0; + + for k in 0..self.rows { + sum += self.get(k, i) * rhs.get(j, k); + } + + matrix.set(i, j, sum); + } + } + + matrix + } else { + panic!("Cannot multiply matrices with incompatible dimensions"); + } + } +} + +impl std::ops::Mul for Matrix { + type Output = Self; + + fn mul(self, rhs: f64) -> Self::Output { + let mut matrix = Self::new(self.rows, self.cols); + + for i in 0..self.rows { + for j in 0..self.cols { + matrix.set(i, j, self.get(i, j) * rhs); + } + } + + matrix + } +} + +impl std::ops::Add for Matrix { + type Output = Self; + + fn add(self, rhs: Self) -> Self::Output { + assert_eq!( + self.rows, rhs.rows, + "Cannot add matrices with different row counts" + ); + assert_eq!( + self.cols, rhs.cols, + "Cannot add matrices with different column counts" + ); + + let mut matrix = Self::new(self.rows, self.cols); + + for i in 0..self.rows { + for j in 0..self.cols { + matrix.set(i, j, self.get(i, j) + rhs.get(i, j)); + } + } + + matrix + } +} diff --git a/src/trueskill.rs b/src/trueskill/mod.rs similarity index 65% rename from src/trueskill.rs rename to src/trueskill/mod.rs index edd1783..9a3554d 100644 --- a/src/trueskill.rs +++ b/src/trueskill/mod.rs @@ -69,13 +69,24 @@ //! - [Moserware: Computing Your Skill](http://www.moserware.com/2010/03/computing-your-skill.html) //! - [TrueSkill Calculator](https://trueskill-calculator.vercel.app/) +mod factor_graph; +mod gaussian; +mod matrix; + +use std::cell::RefCell; use std::f64::consts::{FRAC_1_SQRT_2, PI, SQRT_2}; +use std::rc::Rc; +use factor_graph::{LikelihoodFactor, PriorFactor, SumFactor, TruncateFactor, Variable}; +use gaussian::Gaussian; +use matrix::Matrix; #[cfg(feature = "serde")] use serde::{Deserialize, Serialize}; use crate::{weng_lin::WengLinRating, Outcomes}; -use crate::{Rating, RatingPeriodSystem, RatingSystem, TeamRatingSystem}; +use crate::{MultiTeamOutcome, Rating, RatingPeriodSystem, RatingSystem, TeamRatingSystem}; + +const MIN_DELTA: f64 = 0.0001; #[derive(Copy, Clone, Debug, PartialEq)] #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] @@ -615,13 +626,211 @@ pub fn trueskill_two_teams( (new_team_one, new_team_two) } +#[must_use] +/// Calculates the [`TrueSkillRating`] of multiple teams based on their ratings, uncertainties, and the outcome of the game. +/// +/// Takes in a slice, which contains tuples of teams, which are just slices of [`TrueSkillRating`]s, +/// as well the rank of the team as an [`MultiTeamOutcome`] and a [`TrueSkillConfig`]. +/// +/// Ties are represented by several teams having the same rank. +/// +/// Similar to [`trueskill_two_teams`]. +/// +/// **Caution regarding usage of TrueSkill**: +/// Microsoft permits only Xbox Live games or non-commercial projects to use TrueSkill(TM). +/// If your project is commercial, you should use another rating system included here. +/// +/// # Examples +/// ``` +/// use skillratings::{ +/// trueskill::{trueskill_multi_team, TrueSkillConfig, TrueSkillRating}, +/// MultiTeamOutcome, +/// }; +/// +/// let team_one = vec![ +/// TrueSkillRating::new(), +/// TrueSkillRating { +/// rating: 30.0, +/// uncertainty: 1.2, +/// }, +/// TrueSkillRating { +/// rating: 21.0, +/// uncertainty: 6.5, +/// }, +/// ]; +/// +/// let team_two = vec![ +/// TrueSkillRating::default(), +/// TrueSkillRating { +/// rating: 41.0, +/// uncertainty: 1.4, +/// }, +/// TrueSkillRating { +/// rating: 19.2, +/// uncertainty: 4.3, +/// }, +/// ]; +/// +/// let team_three = vec![ +/// TrueSkillRating::default(), +/// TrueSkillRating { +/// rating: 29.4, +/// uncertainty: 1.6, +/// }, +/// TrueSkillRating { +/// rating: 17.2, +/// uncertainty: 2.1, +/// }, +/// ]; +/// +/// let teams_and_ranks = vec![ +/// (&team_one[..], MultiTeamOutcome::new(2)), // Team 1 takes the second place. +/// (&team_two[..], MultiTeamOutcome::new(1)), // Team 2 takes the first place. +/// (&team_three[..], MultiTeamOutcome::new(3)), // Team 3 takes the third place. +/// ]; +/// +/// let new_teams = trueskill_multi_team(&teams_and_ranks, &TrueSkillConfig::new()); +/// +/// assert_eq!(new_teams.len(), 3); +/// +/// let new_one = &new_teams[0]; +/// let new_two = &new_teams[1]; +/// let new_three = &new_teams[2]; +/// +/// assert!((new_one[0].rating - 25.622_306_739_859_763).abs() < f64::EPSILON); +/// assert!((new_one[1].rating - 30.012_965_086_723_046).abs() < f64::EPSILON); +/// assert!((new_one[2].rating - 21.378635787625903).abs() < f64::EPSILON); +/// +/// assert!((new_two[0].rating - 28.246_057_397_676_047).abs() < f64::EPSILON); +/// assert!((new_two[1].rating - 41.091_932_136_518_125).abs() < f64::EPSILON); +/// assert!((new_two[2].rating - 20.064_520_412_174_183).abs() < f64::EPSILON); +/// +/// assert!((new_three[0].rating - 21.131_635_862_464_19).abs() < f64::EPSILON); +/// assert!((new_three[1].rating - 29.257_024_085_611_56).abs() < f64::EPSILON); +/// assert!((new_three[2].rating - 16.953_981_169_279_245).abs() < f64::EPSILON); +/// ``` +pub fn trueskill_multi_team( + teams_and_ranks: &[(&[TrueSkillRating], MultiTeamOutcome)], + config: &TrueSkillConfig, +) -> Vec> { + if teams_and_ranks.is_empty() { + return Vec::new(); + } + + // Just returning the original teams if a team is empty. + for (team, _) in teams_and_ranks { + if team.is_empty() { + return teams_and_ranks + .iter() + .map(|(team, _)| team.to_vec()) + .collect(); + } + } + + let mut sorted_teams_and_ranks_with_pos = Vec::new(); + for (pos, (team, outcome)) in teams_and_ranks.iter().enumerate() { + sorted_teams_and_ranks_with_pos.push((pos, (*team, *outcome))); + } + sorted_teams_and_ranks_with_pos.sort_by_key(|v| v.1 .1); + + let teams_and_ranks: Vec<(&[TrueSkillRating], MultiTeamOutcome)> = + sorted_teams_and_ranks_with_pos + .iter() + .map(|v| v.1) + .collect(); + + let mut flattened_ratings = Vec::new(); + for (team, _) in &teams_and_ranks { + for player in *team { + flattened_ratings.push(*player); + } + } + + let rating_vars = { + let mut v = Vec::with_capacity(flattened_ratings.len()); + for _ in 0..flattened_ratings.len() { + v.push(Rc::new(RefCell::new(Variable::new()))); + } + + v + }; + let perf_vars = { + let mut v = Vec::with_capacity(flattened_ratings.len()); + for _ in 0..flattened_ratings.len() { + v.push(Rc::new(RefCell::new(Variable::new()))); + } + + v + }; + let team_perf_vars = { + let mut v = Vec::with_capacity(teams_and_ranks.len()); + for _ in 0..teams_and_ranks.len() { + v.push(Rc::new(RefCell::new(Variable::new()))); + } + + v + }; + let team_diff_vars = { + let mut v = Vec::with_capacity(teams_and_ranks.len() - 1); + for _ in 0..(teams_and_ranks.len() - 1) { + v.push(Rc::new(RefCell::new(Variable::new()))); + } + + v + }; + let team_sizes = team_sizes(&teams_and_ranks); + + let rating_layer = run_schedule( + &rating_vars, + &perf_vars, + &team_perf_vars, + &team_diff_vars, + &team_sizes, + &teams_and_ranks, + &flattened_ratings, + config.default_dynamics, + config.beta, + config.draw_probability, + MIN_DELTA, + ); + + let mut transformed_groups = Vec::new(); + let mut iter_team_sizes = vec![0]; + iter_team_sizes.extend_from_slice(&team_sizes[..(team_sizes.len() - 1)]); + + for (start, end) in iter_team_sizes.into_iter().zip(&team_sizes) { + let mut group = Vec::new(); + for f in &rating_layer[start..*end] { + let gaussian = f.variable.borrow().gaussian; + let mu = gaussian.mu(); + let sigma = gaussian.sigma(); + + group.push(TrueSkillRating { + rating: mu, + uncertainty: sigma, + }); + } + + transformed_groups.push(group); + } + + let mut unsorted_with_pos = sorted_teams_and_ranks_with_pos + .iter() + .map(|v| v.0) + .zip(transformed_groups) + .collect::>(); + unsorted_with_pos.sort_by_key(|v| v.0); + + unsorted_with_pos.into_iter().map(|v| v.1).collect() +} + #[must_use] /// Gets the quality of the match, which is equal to the probability that the match will end in a draw. /// The higher the Value, the better the quality of the match. /// /// Takes in two players as [`TrueSkillRating`]s and returns the probability of a draw occurring as an [`f64`] between 1.0 and 0.0. /// -/// Similar to [`match_quality_two_teams`]. +/// Similar to [`match_quality_two_teams`] and [`match_quality_multi_team`]. /// /// # Examples /// ``` @@ -668,7 +877,7 @@ pub fn match_quality( /// /// Takes in two teams as a Slice of [`TrueSkillRating`]s and returns the probability of a draw occurring as an [`f64`] between 1.0 and 0.0. /// -/// Similar to [`match_quality`]. +/// Similar to [`match_quality`] and [`match_quality_multi_team`]. /// /// # Examples /// ``` @@ -727,6 +936,103 @@ pub fn match_quality_two_teams( a * b } +#[must_use] +/// Gets the quality of the match, which is equal to the probability that the match will end in a draw. +/// The higher the Value, the better the quality of the match. +/// +/// Takes in multiple teams as a Slices of [`TrueSkillRating`]s, a [`TrueSkillConfig`] +/// and returns the probability of a draw occurring as an [`f64`] between 1.0 and 0.0. +/// +/// Similar to [`match_quality`] and [`match_quality_two_teams`]. +/// +/// # Examples +/// ``` +/// use skillratings::trueskill::{match_quality_multi_team, TrueSkillConfig, TrueSkillRating}; +/// +/// let team_one = vec![ +/// TrueSkillRating { +/// rating: 20.0, +/// uncertainty: 2.0, +/// }, +/// TrueSkillRating { +/// rating: 25.0, +/// uncertainty: 2.0, +/// }, +/// ]; +/// let team_two = vec![ +/// TrueSkillRating { +/// rating: 35.0, +/// uncertainty: 2.0, +/// }, +/// TrueSkillRating { +/// rating: 20.0, +/// uncertainty: 3.0, +/// }, +/// ]; +/// let team_three = vec![ +/// TrueSkillRating { +/// rating: 20.0, +/// uncertainty: 2.0, +/// }, +/// TrueSkillRating { +/// rating: 22.0, +/// uncertainty: 1.0, +/// }, +/// ]; +/// +/// let quality = match_quality_multi_team( +/// &[&team_one, &team_two, &team_three], +/// &TrueSkillConfig::new(), +/// ); +/// +/// // There is a ~28.6% chance of a draw occurring. +/// assert_eq!(quality, 0.285_578_468_347_742_1); +/// ``` +pub fn match_quality_multi_team(teams: &[&[TrueSkillRating]], config: &TrueSkillConfig) -> f64 { + if teams.is_empty() { + return 0.0; + } + + for team in teams { + if team.is_empty() { + return 0.0; + } + } + + let total_players = teams.iter().map(|t| t.len()).sum::(); + + let team_uncertainties_sq_flatten = teams + .iter() + .flat_map(|team| { + team.iter() + .map(|p| p.uncertainty.powi(2)) + .collect::>() + }) + .collect::>(); + let team_ratings_flatten = teams + .iter() + .flat_map(|team| team.iter().map(|p| p.rating).collect::>()) + .collect::>(); + + let mean_matrix = Matrix::new_from_data(&team_ratings_flatten, total_players, 1); + let variance_matrix = Matrix::new_diagonal(&team_uncertainties_sq_flatten); + + let rotated_a_matrix = Matrix::create_rotated_a_matrix(teams); + let a_matrix = rotated_a_matrix.transpose(); + + let a_ta = rotated_a_matrix.clone() * a_matrix.clone() * config.beta.powi(2); + let atsa = rotated_a_matrix.clone() * variance_matrix * a_matrix.clone(); + let start = a_matrix * mean_matrix.transpose(); + let middle = a_ta.clone() + atsa; + + let end = rotated_a_matrix * mean_matrix; + + let e_arg = (start * middle.inverse() * end * -0.5).determinant(); + let s_arg = a_ta.determinant() / middle.determinant(); + + e_arg.exp() * s_arg.sqrt() +} + #[must_use] /// Calculates the expected outcome of two players based on TrueSkill. /// @@ -735,7 +1041,7 @@ pub fn match_quality_two_teams( /// 1.0 means a certain victory for the player, 0.0 means certain loss. /// Values near 0.5 mean a draw is likely to occur. /// -/// Similar to [`expected_score_two_teams`]. +/// Similar to [`expected_score_two_teams`] and [`expected_score_multi_team`]. /// /// To see the actual chances of a draw occurring, please use [`match_quality`]. /// @@ -791,7 +1097,7 @@ pub fn expected_score( /// 1.0 means a certain victory for the player, 0.0 means certain loss. /// Values near 0.5 mean a draw is likely to occur. /// -/// Similar to [`expected_score`]. +/// Similar to [`expected_score`] and [`expected_score_multi_team`]. /// /// To see the actual chances of a draw occurring, please use [`match_quality_two_teams`]. /// @@ -854,6 +1160,116 @@ pub fn expected_score_two_teams( (exp_one, exp_two) } +#[must_use] +/// Calculates the expected outcome of multiple teams based on TrueSkill. +/// +/// Takes in multiple teams as Slices of [`TrueSkillRating`]s, a [`TrueSkillConfig`] +/// and returns the probability of victory for each team as an [`f64`] between 1.0 and 0.0. +/// 1.0 means a certain victory for the player, 0.0 means certain loss. +/// Values near `1 / Number of Teams` mean a draw is likely to occur. +/// +/// Similar to [`expected_score`] and [`expected_score_multi_team`]. +/// +/// To see the actual chances of a draw occurring, please use [`match_quality_multi_team`]. +/// +/// # Examples +/// ``` +/// use skillratings::trueskill::{expected_score_multi_team, TrueSkillConfig, TrueSkillRating}; +/// +/// let team_one = vec![ +/// TrueSkillRating { +/// rating: 38.0, +/// uncertainty: 3.0, +/// }, +/// TrueSkillRating { +/// rating: 38.0, +/// uncertainty: 3.0, +/// }, +/// ]; +/// +/// let team_two = vec![ +/// TrueSkillRating { +/// rating: 44.0, +/// uncertainty: 3.0, +/// }, +/// TrueSkillRating { +/// rating: 44.0, +/// uncertainty: 3.0, +/// }, +/// ]; +/// +/// let team_three = vec![ +/// TrueSkillRating { +/// rating: 50.0, +/// uncertainty: 3.0, +/// }, +/// TrueSkillRating { +/// rating: 50.0, +/// uncertainty: 3.0, +/// }, +/// ]; +/// +/// let exp = expected_score_multi_team( +/// &[&team_one, &team_two, &team_three], +/// &TrueSkillConfig::new(), +/// ); +/// +/// assert!((exp.iter().sum::() - 1.0).abs() < f64::EPSILON); +/// +/// // Team one has a 6% chance of winning, Team two a 33% and Team three a 61% chance. +/// assert!((exp[0] * 100.0 - 6.0).round().abs() < f64::EPSILON); +/// assert!((exp[1] * 100.0 - 33.0).round().abs() < f64::EPSILON); +/// assert!((exp[2] * 100.0 - 61.0).round().abs() < f64::EPSILON); +/// ``` +pub fn expected_score_multi_team( + teams: &[&[TrueSkillRating]], + config: &TrueSkillConfig, +) -> Vec { + let player_count = teams.iter().map(|t| t.len()).sum::() as f64; + + let mut win_probabilities = Vec::with_capacity(teams.len()); + let mut total_probability = 0.0; + + for (i, team_one) in teams.iter().enumerate() { + // We are calculating the probability of team_one winning against all other teams. + // We do this for every team, sum up the probabilities + // and then divide by the total probability to get the probability of winning for each team. + let mut current_team_probabilities = Vec::with_capacity(teams.len() - 1); + let team_one_ratings = team_one.iter().map(|p| p.rating).sum::(); + let team_one_uncertainties = team_one.iter().map(|p| p.uncertainty.powi(2)).sum::(); + + for (j, team_two) in teams.iter().enumerate() { + if i == j { + continue; + } + + let team_two_ratings = team_two.iter().map(|p| p.rating).sum::(); + let team_two_uncertainties = + team_two.iter().map(|p| p.uncertainty.powi(2)).sum::(); + + let delta = team_one_ratings - team_two_ratings; + let denom = (team_two_uncertainties + + player_count.mul_add(config.beta.powi(2), team_one_uncertainties)) + .sqrt(); + + let result = cdf(delta / denom, 0.0, 1.0); + + current_team_probabilities.push(result); + total_probability += result; + } + + win_probabilities.push(current_team_probabilities); + } + + let mut expected_scores = Vec::new(); + + for probability in win_probabilities { + expected_scores.push(probability.iter().sum::() / total_probability); + } + + expected_scores +} + #[must_use] /// Calculates the expected outcome of a player in a rating period or tournament. /// @@ -942,7 +1358,7 @@ pub fn get_rank(player: &TrueSkillRating) -> f64 { } fn draw_margin(draw_probability: f64, beta: f64, total_players: f64) -> f64 { - inverse_cdf(0.5 * (draw_probability + 1.0), 0.0, 1.0) * total_players.sqrt() * beta + inverse_cdf((draw_probability + 1.0) / 2.0, 0.0, 1.0) * total_players.sqrt() * beta } fn v_non_draw(difference: f64, draw_margin: f64, c: f64) -> f64 { @@ -1134,8 +1550,281 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { ((2.0 * PI).sqrt() * sigma.abs()).recip() * (-(((x - mu) / sigma.abs()).powi(2) / 2.0)).exp() } +#[allow(clippy::too_many_arguments)] +fn run_schedule( + rating_vars: &[Rc>], + perf_vars: &[Rc>], + team_perf_vars: &[Rc>], + team_diff_vars: &[Rc>], + team_sizes: &[usize], + sorted_teams_and_ranks: &[(&[TrueSkillRating], MultiTeamOutcome)], + flattened_ratings: &[TrueSkillRating], + tau: f64, + beta: f64, + draw_probability: f64, + min_delta: f64, +) -> Vec { + assert!(!(min_delta <= 0.0), "min_delta must be greater than 0"); + + let mut id = 0; + + let mut rating_layer = build_rating_layer(rating_vars, flattened_ratings, tau, id); + id += >::try_into(rating_layer.len()).unwrap(); + let mut perf_layer = build_perf_layer(rating_vars, perf_vars, beta, id); + id += >::try_into(perf_layer.len()).unwrap(); + let mut team_perf_layer = build_team_perf_layer(team_perf_vars, perf_vars, team_sizes, id); + id += >::try_into(team_perf_layer.len()).unwrap(); + + for factor in &mut rating_layer { + factor.down(); + } + for factor in &mut perf_layer { + factor.down(); + } + for factor in &mut team_perf_layer { + factor.down(); + } + + let mut team_diff_layer = build_team_diff_layer(team_diff_vars, team_perf_vars, id); + let team_diff_len = team_diff_layer.len(); + id += >::try_into(team_diff_len).unwrap(); + let mut trunc_layer = build_trunc_layer( + team_diff_vars, + sorted_teams_and_ranks, + draw_probability, + beta, + id, + ); + + let mut delta: f64; + for _ in 0..10 { + if team_diff_len == 1 { + team_diff_layer[0].down(); + delta = trunc_layer[0].up(); + } else { + delta = 0.0; + for x in 0..(team_diff_len - 1) { + team_diff_layer[x].down(); + delta = delta.max(trunc_layer[x].up()); + team_diff_layer[x].up(1); + } + for x in (1..team_diff_len).rev() { + team_diff_layer[x].down(); + delta = delta.max(trunc_layer[x].up()); + team_diff_layer[x].up(0); + } + } + if delta <= min_delta { + break; + } + } + + team_diff_layer[0].up(0); + team_diff_layer[team_diff_len - 1].up(1); + for f in &mut team_perf_layer { + for x in 0..f.terms_len() { + f.up(x); + } + } + for f in &mut perf_layer { + f.up(); + } + + rating_layer +} + +fn build_rating_layer( + rating_vars: &[Rc>], + flattened_ratings: &[TrueSkillRating], + tau: f64, + starting_id: u32, +) -> Vec { + let mut v = Vec::with_capacity(rating_vars.len()); + let mut i = starting_id; + for (var, rating) in rating_vars.iter().zip(flattened_ratings) { + v.push(PriorFactor::new( + i, + Rc::clone(var), + Gaussian::with_mu_sigma(rating.rating, rating.uncertainty), + tau, + )); + i += 1; + } + + v +} + +fn build_perf_layer( + rating_vars: &[Rc>], + perf_vars: &[Rc>], + beta: f64, + starting_id: u32, +) -> Vec { + let beta_sq = beta.powi(2); + let mut v = Vec::with_capacity(rating_vars.len()); + let mut i = starting_id; + for (rating_var, perf_var) in rating_vars.iter().zip(perf_vars) { + v.push(LikelihoodFactor::new( + i, + Rc::clone(rating_var), + Rc::clone(perf_var), + beta_sq, + )); + i += 1; + } + + v +} + +fn build_team_perf_layer( + team_perf_vars: &[Rc>], + perf_vars: &[Rc>], + team_sizes: &[usize], + starting_id: u32, +) -> Vec { + let mut v = Vec::with_capacity(team_perf_vars.len()); + let mut i = starting_id; + for (team, team_perf_var) in team_perf_vars.iter().enumerate() { + let start = if team > 0 { team_sizes[team - 1] } else { 0 }; + + let end = team_sizes[team]; + let child_perf_vars = perf_vars[start..end].to_vec(); + let coeffs = vec![1.0; child_perf_vars.len()]; + + v.push(SumFactor::new( + i, + Rc::clone(team_perf_var), + child_perf_vars, + coeffs, + )); + i += 1; + } + + v +} + +fn build_team_diff_layer( + team_diff_vars: &[Rc>], + team_perf_vars: &[Rc>], + starting_id: u32, +) -> Vec { + let mut v = Vec::with_capacity(team_diff_vars.len()); + let mut i = starting_id; + for (team, team_diff_var) in team_diff_vars.iter().enumerate() { + v.push(SumFactor::new( + i, + Rc::clone(team_diff_var), + team_perf_vars[team..(team + 2)].to_vec(), + vec![1.0, -1.0], + )); + i += 1; + } + + v +} + +fn build_trunc_layer( + team_diff_vars: &[Rc>], + sorted_teams_and_ranks: &[(&[TrueSkillRating], MultiTeamOutcome)], + draw_probability: f64, + beta: f64, + starting_id: u32, +) -> Vec { + fn v_w(diff: f64, draw_margin: f64) -> f64 { + let x = diff - draw_margin; + let denom = cdf(x, 0.0, 1.0); + + if denom == 0.0 { + -x + } else { + pdf(x, 0.0, 1.0) / denom + } + } + + fn v_d(diff: f64, draw_margin: f64) -> f64 { + let abs_diff = diff.abs(); + let a = draw_margin - abs_diff; + let b = -draw_margin - abs_diff; + let denom = cdf(a, 0.0, 1.0) - cdf(b, 0.0, 1.0); + let numer = pdf(b, 0.0, 1.0) - pdf(a, 0.0, 1.0); + + let lhs = if denom == 0.0 { a } else { numer / denom }; + let rhs = if diff < 0.0 { -1.0 } else { 1.0 }; + + lhs * rhs + } + + fn w_w(diff: f64, draw_margin: f64) -> f64 { + let x = diff - draw_margin; + let v = v_w(diff, draw_margin); + let w = v * (v + x); + if 0.0 < w && w < 1.0 { + return w; + } + + panic!("floating point error"); + } + + fn w_d(diff: f64, draw_margin: f64) -> f64 { + let abs_diff = diff.abs(); + let a = draw_margin - abs_diff; + let b = -draw_margin - abs_diff; + let denom = cdf(a, 0.0, 1.0) - cdf(b, 0.0, 1.0); + + assert!(!(denom == 0.0), "floating point error"); + + let v = v_d(abs_diff, draw_margin); + v.mul_add(v, (a * pdf(a, 0.0, 1.0) - b * pdf(b, 0.0, 1.0)) / denom) + } + + let mut v = Vec::with_capacity(team_diff_vars.len()); + let mut i = starting_id; + for (x, team_diff_var) in team_diff_vars.iter().enumerate() { + let size = sorted_teams_and_ranks[x..(x + 2)] + .iter() + .map(|v| v.0.len() as f64) + .sum(); + let draw_margin = draw_margin(draw_probability, beta, size); + let v_func: Box f64>; + let w_func: Box f64>; + if sorted_teams_and_ranks[x].1 == sorted_teams_and_ranks[x + 1].1 { + v_func = Box::new(v_d); + w_func = Box::new(w_d); + } else { + v_func = Box::new(v_w); + w_func = Box::new(w_w); + }; + + v.push(TruncateFactor::new( + i, + Rc::clone(team_diff_var), + v_func, + w_func, + draw_margin, + )); + i += 1; + } + + v +} + +fn team_sizes(teams_and_ranks: &[(&[TrueSkillRating], MultiTeamOutcome)]) -> Vec { + let mut team_sizes = Vec::new(); + for (team, _) in teams_and_ranks { + if team_sizes.is_empty() { + team_sizes.push(team.len()); + } else { + team_sizes.push(team.len() + team_sizes[team_sizes.len() - 1]); + } + } + + team_sizes +} + #[cfg(test)] mod tests { + use crate::MultiTeamOutcome; + use super::*; use std::f64::{INFINITY, NEG_INFINITY}; @@ -1544,6 +2233,114 @@ mod tests { assert!(exp2.mul_add(100.0, -20.0).round().abs() < f64::EPSILON); assert!((exp1.mul_add(100.0, exp2 * 100.0).round() - 100.0).abs() < f64::EPSILON); + + // Testing if the other functions give the same result. + let team_one = [TrueSkillRating::from((44.0, 3.0))]; + let team_two = [TrueSkillRating::from((38.0, 3.0))]; + + let (e0, e1) = expected_score_two_teams(&team_one, &team_two, &TrueSkillConfig::new()); + let e = expected_score_multi_team(&[&team_one, &team_two], &TrueSkillConfig::new()); + + assert!((e0 - e[0]).abs() < f64::EPSILON); + assert!((e1 - e[1]).abs() < f64::EPSILON); + assert!((exp1 - e[0]).abs() < f64::EPSILON); + assert!((exp2 - e[1]).abs() < f64::EPSILON); + } + + #[test] + fn test_match_quality_multi_team() { + let team_one = vec![TrueSkillRating::new(); 2]; + let team_two = vec![TrueSkillRating::from((30.0, 3.0)); 2]; + let team_three = vec![TrueSkillRating::from((40.0, 2.0)); 2]; + + let exp = match_quality_multi_team( + &[&team_one, &team_two, &team_three], + &TrueSkillConfig::new(), + ); + + // Double checked this with the most popular python implementation. + assert!((exp - 0.017_538_349_223_941_27).abs() < f64::EPSILON); + + let exp = match_quality_multi_team(&[], &TrueSkillConfig::default()); + + assert!(exp < f64::EPSILON); + + let exp = match_quality_multi_team(&[&team_one, &[]], &TrueSkillConfig::default()); + + assert!(exp < f64::EPSILON); + } + + #[test] + fn test_multi_team_expected() { + let team_one = vec![ + TrueSkillRating { + rating: 38.0, + uncertainty: 3.0, + }, + TrueSkillRating { + rating: 38.0, + uncertainty: 3.0, + }, + ]; + + let team_two = vec![ + TrueSkillRating { + rating: 44.0, + uncertainty: 3.0, + }, + TrueSkillRating { + rating: 44.0, + uncertainty: 3.0, + }, + ]; + + let team_three = vec![ + TrueSkillRating { + rating: 50.0, + uncertainty: 3.0, + }, + TrueSkillRating { + rating: 50.0, + uncertainty: 3.0, + }, + ]; + + let exp = expected_score_multi_team( + &[&team_one, &team_two, &team_three], + &TrueSkillConfig::new(), + ); + + assert!((exp.iter().sum::() - 1.0).abs() < f64::EPSILON); + + assert_eq!( + exp, + vec![ + 0.058_904_655_169_257_615, + 0.333_333_333_333_333_3, + 0.607_762_011_497_409 + ] + ); + + let team_one = vec![TrueSkillRating::new(); 10]; + let team_two = vec![TrueSkillRating::new(); 10]; + let team_three = vec![TrueSkillRating::new(); 10]; + let team_four = vec![TrueSkillRating::new(); 10]; + + let exp = expected_score_multi_team( + &[&team_one, &team_two, &team_three, &team_four], + &TrueSkillConfig::new(), + ); + + assert!((exp.iter().sum::() - 1.0).abs() < f64::EPSILON); + assert_eq!( + exp, + vec![ + 0.249_999_999_999_999_97, + 0.249_999_999_999_999_97, + 0.249_999_999_999_999_97, + 0.249_999_999_999_999_97 + ] + ); } #[test] @@ -1640,6 +2437,26 @@ mod tests { assert!(v3 == NEG_INFINITY); } + #[test] + fn test_matrix_panics() { + use std::panic::catch_unwind; + + let result = catch_unwind(|| Matrix::new(2, 3).determinant()); + assert!(result.is_err()); + + let result = catch_unwind(|| Matrix::new(2, 2).inverse()); + assert!(result.is_err()); + + let result = catch_unwind(|| Matrix::new(2, 2) * Matrix::new(3, 3)); + assert!(result.is_err()); + + let result = catch_unwind(|| Matrix::new(3, 2) + Matrix::new(2, 2)); + assert!(result.is_err()); + + let result = catch_unwind(|| Matrix::new(2, 2) + Matrix::new(2, 3)); + assert!(result.is_err()); + } + #[test] #[allow(clippy::clone_on_copy)] fn test_misc_stuff() { @@ -1652,6 +2469,8 @@ mod tests { assert!(!format!("{player_one:?}").is_empty()); assert!(!format!("{config:?}").is_empty()); + assert!(!format!("{:?}", Matrix::new(2, 3)).is_empty()); + assert_eq!(player_one, TrueSkillRating::from((25.0, 25.0 / 3.0))); } @@ -1705,4 +2524,77 @@ mod tests { assert!((exp1 + exp2 - 1.0).abs() < f64::EPSILON); } + + #[test] + /// This example is taken from the Python TrueSkill package: + /// https://github.com/sublee/trueskill/blob/3ff78b26c8374b30cc0d6bae32cb45c9bee46a1d/trueskilltest.py#L311 + fn test_trueskill_multi_team() { + let t1p1 = TrueSkillRating { + rating: 40.0, + uncertainty: 4.0, + }; + let t1p2 = TrueSkillRating { + rating: 45.0, + uncertainty: 3.0, + }; + + let t2p1 = TrueSkillRating { + rating: 20.0, + uncertainty: 7.0, + }; + let t2p2 = TrueSkillRating { + rating: 19.0, + uncertainty: 6.0, + }; + let t2p3 = TrueSkillRating { + rating: 30.0, + uncertainty: 9.0, + }; + let t2p4 = TrueSkillRating { + rating: 10.0, + uncertainty: 4.0, + }; + + let t3p1 = TrueSkillRating { + rating: 50.0, + uncertainty: 5.0, + }; + let t3p2 = TrueSkillRating { + rating: 30.0, + uncertainty: 2.0, + }; + + let t1 = [t1p1, t1p2]; + let t2 = [t2p1, t2p2, t2p3, t2p4]; + let t3 = [t3p1, t3p2]; + + let teams_and_ranks = [ + (&t1[..], MultiTeamOutcome::new(0)), + (&t2[..], MultiTeamOutcome::new(1)), + (&t3[..], MultiTeamOutcome::new(1)), + ]; + let results = trueskill_multi_team(&teams_and_ranks, &TrueSkillConfig::new()); + + assert!((results[0][0].rating - 40.876_849_177_315_655).abs() < f64::EPSILON); + assert!((results[0][1].rating - 45.493_394_092_398_44).abs() < f64::EPSILON); + + assert!((results[1][0].rating - 19.608_650_920_845_236).abs() < f64::EPSILON); + assert!((results[1][1].rating - 18.712_463_514_890_54).abs() < f64::EPSILON); + assert!((results[1][2].rating - 29.353_112_227_810_637).abs() < f64::EPSILON); + assert!((results[1][3].rating - 9.872_175_198_037_164).abs() < f64::EPSILON); + + assert!((results[2][0].rating - 48.829_832_201_455_32).abs() < f64::EPSILON); + assert!((results[2][1].rating - 29.812_500_188_903_005).abs() < f64::EPSILON); + + assert!((results[0][0].uncertainty - 3.839_527_589_355_37).abs() < f64::EPSILON); + assert!((results[0][1].uncertainty - 2.933_671_613_522_051).abs() < f64::EPSILON); + + assert!((results[1][0].uncertainty - 6.396_044_310_523_897).abs() < f64::EPSILON); + assert!((results[1][1].uncertainty - 5.624_556_429_622_889).abs() < f64::EPSILON); + assert!((results[1][2].uncertainty - 7.673_456_361_986_594).abs() < f64::EPSILON); + assert!((results[1][3].uncertainty - 3.891_408_425_994_520_3).abs() < f64::EPSILON); + + assert!((results[2][0].uncertainty - 4.590_018_525_151_38).abs() < f64::EPSILON); + assert!((results[2][1].uncertainty - 1.976_314_792_712_798_2).abs() < f64::EPSILON); + } }