From be74c1eac72fea14316df6756f1439e686151ff1 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 24 Oct 2018 11:12:34 +0200 Subject: [PATCH] It works! --- src/factor_graph.rs | 48 +++++------------ src/lib.rs | 66 +++++++++++------------ src/main.rs | 15 ++++++ src/math.rs | 129 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 190 insertions(+), 68 deletions(-) create mode 100644 src/main.rs diff --git a/src/factor_graph.rs b/src/factor_graph.rs index 73d8b85..1d4cc1b 100644 --- a/src/factor_graph.rs +++ b/src/factor_graph.rs @@ -23,7 +23,7 @@ impl VariableArena { let index = self.variables.len(); self.variables.push(Variable { - value: Gaussian::from_precision(0.0, 0.0), + value: Gaussian::from_pi_tau(0.0, 0.0), factors: HashMap::new(), }); @@ -46,26 +46,19 @@ pub struct Variable { impl Variable { pub fn attach_factor(&mut self, factor: usize) { - self.factors.insert(factor, Gaussian::from_precision(0.0, 0.0)); + self.factors.insert(factor, Gaussian::from_pi_tau(0.0, 0.0)); } pub fn update_value(&mut self, factor: usize, value: Gaussian) { let old = self.factors[&factor]; let intermediate = value * old; - let new = intermediate / self.value; + self.factors.insert(factor, intermediate / self.value); - /* - println!("update_value: old={}, value={:?}, new={:?}, self.value={:?}", - render_g(old), - value, - new, - self.value); - */ - self.value = new; + self.value = value; } - pub fn get_value(&self) -> Gaussian { + pub fn get_value(&self) -> Gaussian { self.value } @@ -77,8 +70,6 @@ impl Variable { self.value = value; - println!("update_message: old={:?} msg={:?}, new={:?}", old, message, value); - self.factors.insert(factor, message); } @@ -110,8 +101,6 @@ impl PriorFactor { } pub fn start(&self, variable_arena: &mut VariableArena) { - println!("PriorFactor: variable.id={}, msg={:?}", self.variable.index, self.gaussian); - variable_arena.get_mut(self.variable).unwrap().update_value(self.id, self.gaussian); } } @@ -153,11 +142,9 @@ impl LikelihoodFactor { .map(|variable| variable.get_message(self.id)) .unwrap(); - let a = 1.0 / (1.0 + self.variance * (x.precision_mean() - fx.precision_mean())); + let a = 1.0 / (1.0 + self.variance * (x.pi() - fx.pi())); - let gaussian = Gaussian::from_precision(a * (x.precision_mean() - fx.precision_mean()), a * (x.precision() - fx.precision())); - - println!("LikelihoodFactor: mean.id={}, msg={:?}", self.mean.index, gaussian); + let gaussian = Gaussian::from_pi_tau(a * (x.pi() - fx.pi()), a * (x.tau() - fx.tau())); variable_arena.get_mut(self.mean).unwrap().update_message(self.id, gaussian); } @@ -173,12 +160,9 @@ impl LikelihoodFactor { .map(|variable| variable.get_message(self.id)) .unwrap(); - let a = 1.0 / (1.0 + self.variance * (y.precision_mean() - fy.precision_mean())); + let a = 1.0 / (1.0 + self.variance * (y.pi() - fy.pi())); - let gaussian = Gaussian::from_precision(a * (y.precision_mean() - fy.precision_mean()), a * (y.precision() - fy.precision())); - - - println!("LikelihoodFactor: value.id={}, msg={:?}", self.value.index, gaussian); + let gaussian = Gaussian::from_pi_tau(a * (y.pi() - fy.pi()), a * (y.tau() - fy.tau())); variable_arena.get_mut(self.value).unwrap().update_message(self.id, gaussian); } @@ -240,13 +224,7 @@ impl SumFactor { let gaussian = Gaussian::from_pi_tau(new_pi, new_tau); - if variable == self.sum { - println!("SumFactor: sum.id={}, msg={:?}", variable.index, gaussian); - } else { - println!("SumFactor: term.id={}, msg={:?}", variable.index, gaussian); - } - - variable_arena.get_mut(variable).unwrap().update_value(self.id, gaussian); + variable_arena.get_mut(variable).unwrap().update_message(self.id, gaussian); } pub fn update_sum(&self, variable_arena: &mut VariableArena) { @@ -298,16 +276,16 @@ impl SumFactor { v[index] = self.sum; - for term in v { + for term in &v { let value = variable_arena - .get(term) + .get(*term) .map(|variable| variable.get_value()) .unwrap(); y.push(value); let value = variable_arena - .get(term) + .get(*term) .map(|variable| variable.get_message(self.id)) .unwrap(); diff --git a/src/lib.rs b/src/lib.rs index 3d26f38..ee02d74 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -7,14 +7,14 @@ mod math; mod matrix; use factor_graph::*; -use gaussian::Gaussian; +pub use gaussian::Gaussian; use matrix::Matrix; /// Default initial mean of ratings. -const MU: f64 = 25.0; +pub const MU: f64 = 25.0; /// Default initial standard deviation of ratings. -const SIGMA: f64 = MU / 3.0; +pub const SIGMA: f64 = MU / 3.0; /// Default distance that guarantees about 76% chance of winning. const BETA: f64 = SIGMA / 2.0; @@ -28,6 +28,22 @@ const DRAW_PROBABILITY: f64 = 0.10; /// A basis to check reliability of the result. const DELTA: f64 = 0.0001; + +fn draw_margin(p: f64, beta: f64, total_players: f64) -> f64 { + math::icdf((p + 1.0) / 2.0) * total_players.sqrt() * beta +} + +/* +Constants::Constants() { + double INITIAL_MU = 25.0; + double INITIAL_SIGMA = INITIAL_MU / 3.0; + double TOTAL_PLAYERS = 2.0; + + this->BETA = INITIAL_SIGMA / 2.0; + this->EPSILON = draw_margin(0.1, this->BETA, TOTAL_PLAYERS); + this->GAMMA = INITIAL_SIGMA / 100.0; +} +*/ pub fn rate(rating_groups: &[&[Gaussian]]) { let flatten_ratings = rating_groups .iter() @@ -62,7 +78,7 @@ pub fn rate(rating_groups: &[&[Gaussian]]) { for (i, rating) in flatten_ratings.iter().enumerate() { let variable = ss[i]; - let gaussian = Gaussian::from_mean_variance(rating.mean(), (rating.variance() + tau_sqr).sqrt()); + let gaussian = Gaussian::from_mu_sigma(rating.mu(), (rating.sigma().powi(2) + tau_sqr).sqrt()); skill.push(PriorFactor::new( &mut variable_arena, @@ -118,12 +134,14 @@ pub fn rate(rating_groups: &[&[Gaussian]]) { let mut trunc = Vec::new(); + let epsilon = draw_margin(0.1, BETA, 2.0); + for i in 0..size - 1 { trunc.push(TruncateFactor::new( &mut variable_arena, factor_id, ds[i], - DELTA, + epsilon, false, )); @@ -134,22 +152,6 @@ pub fn rate(rating_groups: &[&[Gaussian]]) { factor.start(&mut variable_arena); } - /* - println!("before:"); - - for i in 0..size { - let value = variable_arena - .get(ss[i]) - .map(|variable| variable.get_value()) - .unwrap(); - - let mu = value.mean(); - let sigma = value.std_dev(); - - println!("* player={}, mu={}, sigma={}", i, mu, sigma); - } - */ - for factor in &skill_to_perf { factor.update_value(&mut variable_arena); } @@ -181,7 +183,6 @@ pub fn rate(rating_groups: &[&[Gaussian]]) { factor.update_mean(&mut variable_arena); } - /* println!("after:"); for i in 0..size { @@ -190,12 +191,11 @@ pub fn rate(rating_groups: &[&[Gaussian]]) { .map(|variable| variable.get_value()) .unwrap(); - let mu = value.mean(); - let sigma = value.std_dev(); + let mu = value.mu(); + let sigma = value.sigma(); println!("* player={}, mu={}, sigma={}", i, mu, sigma); } - */ } pub fn quality(rating_groups: &[&[Gaussian]]) -> f64 { @@ -211,13 +211,13 @@ pub fn quality(rating_groups: &[&[Gaussian]]) -> f64 { let mut mean_matrix = Matrix::new(length, 1); for (i, rating) in flatten_ratings.iter().enumerate() { - mean_matrix[(i, 0)] = rating.mean(); + mean_matrix[(i, 0)] = rating.mu(); } let mut variance_matrix = Matrix::new(length, length); for (i, rating) in flatten_ratings.iter().enumerate() { - variance_matrix[(i, i)] = rating.std_dev().powi(2); + variance_matrix[(i, i)] = rating.sigma().powi(2); } let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length); @@ -261,18 +261,18 @@ mod tests { #[test] fn test_quality_1vs1() { - let alice = Gaussian::from_mean_std_dev(MU, SIGMA); - let bob = Gaussian::from_mean_std_dev(MU, SIGMA); + let alice = Gaussian::from_mu_sigma(MU, SIGMA); + let bob = Gaussian::from_mu_sigma(MU, SIGMA); assert_eq!(quality(&[&[alice], &[bob]]), 0.4472135954999579); } #[test] fn test_rate_1vs1() { - let alice = Gaussian::from_mean_std_dev(MU, SIGMA); - let bob = Gaussian::from_mean_std_dev(MU, SIGMA); - let chris = Gaussian::from_mean_std_dev(MU, SIGMA); - let darren = Gaussian::from_mean_std_dev(MU, SIGMA); + let alice = Gaussian::from_mu_sigma(MU, SIGMA); + let bob = Gaussian::from_mu_sigma(MU, SIGMA); + let chris = Gaussian::from_mu_sigma(MU, SIGMA); + let darren = Gaussian::from_mu_sigma(MU, SIGMA); // println!("alice: {:?}", alice); // println!("bob: {:?}", alice); diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..9991c5a --- /dev/null +++ b/src/main.rs @@ -0,0 +1,15 @@ +extern crate trueskill; + +use trueskill::{rate, Gaussian, MU, SIGMA}; + +fn main() { + let alice = Gaussian::from_mu_sigma(MU, SIGMA); + let bob = Gaussian::from_mu_sigma(MU, SIGMA); + let chris = Gaussian::from_mu_sigma(MU, SIGMA); + let darren = Gaussian::from_mu_sigma(MU, SIGMA); + + // println!("alice: {:?}", alice); + // println!("bob: {:?}", alice); + + rate(&[&[alice], &[bob], &[chris], &[darren]]); +} diff --git a/src/math.rs b/src/math.rs index 1701474..598297e 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,5 +1,130 @@ use statrs::distribution::{Normal, Univariate, Continuous}; +const S2PI: f64 = 2.50662827463100050242E0; + +const P0: [f64; 5] = [ + -5.99633501014107895267E1, + 9.80010754185999661536E1, + -5.66762857469070293439E1, + 1.39312609387279679503E1, + -1.23916583867381258016E0, +]; + +const Q0: [f64; 8] = [ + 1.95448858338141759834E0, + 4.67627912898881538453E0, + 8.63602421390890590575E1, + -2.25462687854119370527E2, + 2.00260212380060660359E2, + -8.20372256168333339912E1, + 1.59056225126211695515E1, + -1.18331621121330003142E0, +]; + +const P1: [f64; 9] = [ + 4.05544892305962419923E0, + 3.15251094599893866154E1, + 5.71628192246421288162E1, + 4.40805073893200834700E1, + 1.46849561928858024014E1, + 2.18663306850790267539E0, + -1.40256079171354495875E-1, + -3.50424626827848203418E-2, + -8.57456785154685413611E-4, +]; + +const Q1: [f64; 8] = [ + 1.57799883256466749731E1, + 4.53907635128879210584E1, + 4.13172038254672030440E1, + 1.50425385692907503408E1, + 2.50464946208309415979E0, + -1.42182922854787788574E-1, + -3.80806407691578277194E-2, + -9.33259480895457427372E-4, +]; + +const P2: [f64; 9] = [ + 3.23774891776946035970E0, + 6.91522889068984211695E0, + 3.93881025292474443415E0, + 1.33303460815807542389E0, + 2.01485389549179081538E-1, + 1.23716634817820021358E-2, + 3.01581553508235416007E-4, + 2.65806974686737550832E-6, + 6.23974539184983293730E-9, +]; + +const Q2: [f64; 8] = [ + 6.02427039364742014255E0, + 3.67983563856160859403E0, + 1.37702099489081330271E0, + 2.16236993594496635890E-1, + 1.34204006088543189037E-2, + 3.28014464682127739104E-4, + 2.89247864745380683936E-6, + 6.79019408009981274425E-9, +]; + +fn polevl(x: f64, coef: &[f64], n: usize) -> f64 { + let mut ans = coef[0]; + + for i in 1..n + 1 { + ans = ans * x + coef[i]; + } + + ans +} + +fn p1evl(x: f64, coef: &[f64], n: usize) -> f64 { + let mut ans = x + coef[0]; + + for i in 1..n { + ans = ans * x + coef[i]; + } + + ans +} + +fn ndtri(y0: f64) -> f64 { + let mut code = 1; + let mut y = y0; + + if y > (1.0 - 0.13533528323661269189) { + y = 1.0 - y; + code = 0; + } + + if y > 0.13533528323661269189 { + y = y - 0.5; + let y2 = y * y; + let x = y + y * (y2 * polevl(y2, &P0, 4) / p1evl(y2, &Q0, 8)); + let x = x * S2PI; + + return x; + } + + let x = (-2.0 * y.ln()).sqrt(); + let x0 = x - x.ln() / x; + + let z = 1.0 / x; + + let x1 = if x < 8.0 { + z * polevl(z, &P1, 8) / p1evl(z, &Q1, 8) + } else { + z * polevl(z, &P2, 8) / p1evl(z, &Q2, 8) + }; + + let mut x = x0 - x1; + + if code != 0 { + x = -x; + } + + x +} + pub fn cdf(x: f64) -> f64 { Normal::new(0.0, 1.0).unwrap().cdf(x) } @@ -8,6 +133,10 @@ pub fn pdf(x: f64) -> f64 { Normal::new(0.0, 1.0).unwrap().pdf(x) } +pub fn icdf(x: f64) -> f64 { + ndtri(x) +} + #[cfg(test)] mod tests {