From 1b840e737d7f31138442bc5282724f7899cfaa88 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Tue, 23 Oct 2018 14:58:30 +0200 Subject: [PATCH] It works?! --- src/factor_graph.rs | 400 ++++++++++++++++++++++++++++++++++++++++---- src/gaussian.rs | 43 ++++- src/lib.rs | 196 ++++++++++++++++++---- src/math.rs | 30 ++++ 4 files changed, 607 insertions(+), 62 deletions(-) create mode 100644 src/math.rs diff --git a/src/factor_graph.rs b/src/factor_graph.rs index e6f6ce7..124c366 100644 --- a/src/factor_graph.rs +++ b/src/factor_graph.rs @@ -1,54 +1,396 @@ -use std::cmp; +use std::collections::HashMap; use gaussian::Gaussian; +use math; + +#[derive(Clone, Copy)] +pub struct VariableId { + index: usize, +} + +pub struct VariableArena { + variables: Vec, +} + +impl VariableArena { + pub fn new() -> VariableArena { + VariableArena { + variables: Vec::new(), + } + } + + pub fn create(&mut self) -> VariableId { + let index = self.variables.len(); + + self.variables.push(Variable { + value: Gaussian::with_pi_tau(0.0, 0.0), + factors: HashMap::new(), + }); + + VariableId { index } + } + + pub fn get(&mut self, id: VariableId) -> Option<&Variable> { + self.variables.get(id.index) + } + + pub fn get_mut(&mut self, id: VariableId) -> Option<&mut Variable> { + self.variables.get_mut(id.index) + } +} pub struct Variable { - gaussian: Gaussian, + value: Gaussian, + factors: HashMap, } impl Variable { - pub fn new() -> Variable { - Variable { - gaussian: Gaussian::new(0.0, 0.0), - } + pub fn attach_factor(&mut self, factor: usize) { + self.factors.insert(factor, Gaussian::new()); } - fn delta(&self, other: &Variable) -> f32 { - let pi_delta = self.gaussian.pi - other.gaussian.pi; + pub fn update_value(&mut self, factor: usize, value: Gaussian) { + let old = self.factors[&factor]; - if pi_delta.is_infinite() { - 0.0 - } else { - let tau_delta = (self.gaussian.tau - other.gaussian.tau).abs(); + let intermediate = value * old; + let value = intermediate / self.value; - if pi_delta > tau_delta { - pi_delta - } else { - tau_delta - } - } + self.value = value; } -} -pub trait Factor { - fn down(&self) -> f32 { - 0.0 + pub fn get_value(&self) -> Gaussian { + self.value + } + + pub fn update_message(&mut self, factor: usize, message: Gaussian) { + let old = self.factors[&factor]; + + let intermediate = self.value / old; + let value = intermediate * message; + + self.value = value; + + self.factors.insert(factor, message); + } + + pub fn get_message(&self, factor: usize) -> Gaussian { + self.factors[&factor] } } pub struct PriorFactor { - variable: Variable, - dynamic: f32 + id: usize, + variable: VariableId, + gaussian: Gaussian, } impl PriorFactor { - pub fn new(variable: Variable, dynamic: f32) -> PriorFactor { - PriorFactor { variable, dynamic } + pub fn new( + variable_arena: &mut VariableArena, + id: usize, + variable: VariableId, + gaussian: Gaussian, + ) -> PriorFactor { + if let Some(variable) = variable_arena.get_mut(variable) { + variable.attach_factor(id); + } + + PriorFactor { + id, + variable, + gaussian, + } + } + + pub fn start(&self, variable_arena: &mut VariableArena) { + if let Some(variable) = variable_arena.get_mut(self.variable) { + variable.update_value(self.id, self.gaussian); + } } } -impl Factor for PriorFactor { - fn down(&self) -> f32 { - 0.0 +pub struct LikelihoodFactor { + id: usize, + mean: VariableId, + value: VariableId, + variance: f32, +} + +impl LikelihoodFactor { + pub fn new( + variable_arena: &mut VariableArena, + id: usize, + mean: VariableId, + value: VariableId, + variance: f32, + ) -> LikelihoodFactor { + if let Some(variable) = variable_arena.get_mut(mean) { + variable.attach_factor(id); + } + + if let Some(variable) = variable_arena.get_mut(value) { + variable.attach_factor(id); + } + + LikelihoodFactor { + id, + mean, + value, + variance, + } + } + + pub fn update_mean(&self, variable_arena: &mut VariableArena) { + let x = variable_arena + .get(self.value) + .map(|variable| variable.get_value()) + .unwrap(); + let fx = variable_arena + .get_mut(self.value) + .map(|variable| variable.get_message(self.id)) + .unwrap(); + + let a = 1.0 / (1.0 + self.variance * (x.pi - fx.pi)); + + let gaussian = Gaussian::with_pi_tau(a * (x.pi - fx.pi), a * (x.tau - fx.tau)); + + if let Some(variable) = variable_arena.get_mut(self.mean) { + variable.update_message(self.id, gaussian); + } + } + + pub fn update_value(&self, variable_arena: &mut VariableArena) { + let y = variable_arena + .get(self.mean) + .map(|variable| variable.get_value()) + .unwrap(); + let fy = variable_arena + .get(self.mean) + .map(|variable| variable.get_message(self.id)) + .unwrap(); + + let a = 1.0 / (1.0 + self.variance * (y.pi - fy.pi)); + + let gaussian = Gaussian::with_pi_tau(a * (y.pi - fy.pi), a * (y.tau - fy.tau)); + + if let Some(variable) = variable_arena.get_mut(self.value) { + variable.update_message(self.id, gaussian); + } + } +} + +pub struct SumFactor { + id: usize, + sum: VariableId, + terms: Vec, + coeffs: Vec, +} + +impl SumFactor { + pub fn new( + variable_arena: &mut VariableArena, + id: usize, + sum: VariableId, + terms: Vec, + coeffs: Vec, + ) -> SumFactor { + if let Some(variable) = variable_arena.get_mut(sum) { + variable.attach_factor(id); + } + + for term in &terms { + if let Some(variable) = variable_arena.get_mut(*term) { + variable.attach_factor(id); + } + } + + SumFactor { + id, + sum, + terms, + coeffs, + } + } + + fn internal_update( + &self, + variable_arena: &mut VariableArena, + variable: VariableId, + y: Vec, + fy: Vec, + a: &Vec, + ) { + let size = a.len(); + + let mut sum_pi = 0.0; + let mut sum_tau = 0.0; + + for i in 0..size { + let da = a[i]; + let gy = y[i]; + let gfy = fy[i]; + + sum_pi += da.powi(2) / (gy.pi - gfy.pi); + sum_tau += da * (gy.tau - gfy.tau) / (gy.pi - gfy.pi); + } + + let new_pi = 1.0 / sum_pi; + let new_tau = new_pi * sum_tau; + + let gaussian = Gaussian::with_pi_tau(new_pi, new_tau); + + if let Some(variable) = variable_arena.get_mut(variable) { + variable.update_value(self.id, gaussian); + } + } + + pub fn update_sum(&self, variable_arena: &mut VariableArena) { + let mut y = Vec::new(); + + for term in &self.terms { + let value = variable_arena + .get(*term) + .map(|variable| variable.get_value()) + .unwrap(); + + y.push(value); + } + + let mut fy = Vec::new(); + + for term in &self.terms { + let value = variable_arena + .get(*term) + .map(|variable| variable.get_message(self.id)) + .unwrap(); + + fy.push(value); + } + + self.internal_update(variable_arena, self.sum, y, fy, &self.coeffs); + } + + pub fn update_term(&self, variable_arena: &mut VariableArena, index: usize) { + let size = self.coeffs.len(); + let idx_coeff = self.coeffs[index]; + + let mut a = vec![0.0; size]; + + for i in 0..size { + if i != index { + a[i] = -self.coeffs[i] / idx_coeff; + } + } + + a[index] = 1.0 / idx_coeff; + + let idx_term = self.terms[index]; + + let mut y = Vec::new(); + let mut fy = Vec::new(); + + let mut v = self.terms.clone(); + + v[index] = self.sum; + + for term in v { + let value = variable_arena + .get(term) + .map(|variable| variable.get_value()) + .unwrap(); + + y.push(value); + + let value = variable_arena + .get(term) + .map(|variable| variable.get_message(self.id)) + .unwrap(); + + fy.push(value); + } + + self.internal_update(variable_arena, idx_term, y, fy, &a); + } +} + +fn v_win(t: f32, e: f32) -> f32 { + math::pdf(t - e) / math::cdf(t - e) +} + +fn w_win(t: f32, e: f32) -> f32 { + let vwin = v_win(t, e); + + vwin * (vwin + t - e) +} + +fn v_draw(t: f32, e: f32) -> f32 { + (math::pdf(-e - t) - math::pdf(e - t)) / (math::cdf(e - t) - math::cdf(-e - t)) +} + +fn w_draw(t: f32, e: f32) -> f32 { + let vdraw = v_draw(t, e); + let n = (vdraw * vdraw) + ((e - t) * math::pdf(e - t) + (e + t) * math::pdf(e + t)); + let d = math::cdf(e - t) - math::cdf(-e - t); + + n / d +} + +pub struct TruncateFactor { + id: usize, + variable: VariableId, + epsilon: f32, + draw: bool, +} + +impl TruncateFactor { + pub fn new( + variable_arena: &mut VariableArena, + id: usize, + variable: VariableId, + epsilon: f32, + draw: bool, + ) -> TruncateFactor { + if let Some(variable) = variable_arena.get_mut(variable) { + variable.attach_factor(id); + } + + TruncateFactor { + id, + variable, + epsilon, + draw, + } + } + + pub fn update(&self, variable_arena: &mut VariableArena) { + let x = variable_arena + .get(self.variable) + .map(|variable| variable.get_value()) + .unwrap(); + let fx = variable_arena + .get_mut(self.variable) + .map(|variable| variable.get_message(self.id)) + .unwrap(); + + let c = x.pi - fx.pi; + let d = x.tau - fx.tau; + let sqrt_c = c.sqrt(); + + let t = d / sqrt_c; + let e = self.epsilon * sqrt_c; + + let (v, w) = if self.draw { + (v_draw(t, e), w_draw(t, e)) + } else { + (v_win(t, e), w_win(t, e)) + }; + + let m_w = 1.0 - w; + + let gaussian = Gaussian::with_pi_tau(c / m_w, (d + sqrt_c * v) / m_w); + + if let Some(variable) = variable_arena.get_mut(self.variable) { + variable.update_value(self.id, gaussian); + } } } diff --git a/src/gaussian.rs b/src/gaussian.rs index 702b677..e5f3542 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -1,10 +1,51 @@ +use std::ops; + +#[derive(Clone, Copy)] pub struct Gaussian { pub pi: f32, pub tau: f32, } impl Gaussian { - pub fn new(pi: f32, tau: f32) -> Gaussian { + pub fn new() -> Gaussian { + Gaussian::with_pi_tau(0.0, 0.0) + } + + pub fn with_pi_tau(pi: f32, tau: f32) -> Gaussian { Gaussian { pi, tau } } + + pub fn with_mu_sigma(mu: f32, sigma: f32) -> Gaussian { + let pi = 1.0 / sigma.powi(2); + + Gaussian::with_pi_tau(pi, pi * mu) + } + + pub fn mu(&self) -> f32 { + if self.pi == 0.0 { + 0.0 + } else { + self.tau / self.pi + } + } + + pub fn sigma(&self) -> f32 { + (1.0 / self.pi).sqrt() + } +} + +impl ops::Mul for Gaussian { + type Output = Gaussian; + + fn mul(self, rhs: Gaussian) -> Gaussian { + Gaussian::with_pi_tau(self.pi + rhs.pi, self.tau + rhs.tau) + } +} + +impl ops::Div for Gaussian { + type Output = Gaussian; + + fn div(self, rhs: Gaussian) -> Gaussian { + Gaussian::with_pi_tau(self.pi - rhs.pi, self.tau - rhs.tau) + } } diff --git a/src/lib.rs b/src/lib.rs index 812f6c6..5519acf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,9 +1,11 @@ -mod matrix; mod factor_graph; mod gaussian; +mod math; +mod matrix; -use matrix::Matrix; use factor_graph::*; +use gaussian::Gaussian; +use matrix::Matrix; /// Default initial mean of ratings. const MU: f32 = 25.0; @@ -24,9 +26,9 @@ const DRAW_PROBABILITY: f32 = 0.10; const DELTA: f32 = 0.0001; #[derive(Debug, PartialEq)] -struct Rating { - mu: f32, - sigma: f32, +pub struct Rating { + pub mu: f32, + pub sigma: f32, } impl Default for Rating { @@ -38,42 +40,155 @@ impl Default for Rating { } } - -fn _team_sizes(rating_groups: &[&[Rating]]) -> Vec { - let mut team_sizes = Vec::new(); - - for group in rating_groups { - let last = team_sizes.last().map(|v| *v).unwrap_or(0); - - team_sizes.push(group.len() + last); - } - - team_sizes -} - -fn factor_graph_builders(rating_groups: &[&[Rating]]) { +fn rate(rating_groups: &[&[Rating]]) { let flatten_ratings = rating_groups .iter() .flat_map(|group| group.iter()) .collect::>(); - let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice(); - let size = flatten_ratings.len(); - let group_size = rating_groups.len(); - let rating_vars = (0..size).map(|_| Variable::new()).collect::>(); - let perf_vars = (0..size).map(|_| Variable::new()).collect::>(); - let team_perf_vars = (0..group_size).map(|_| Variable::new()).collect::>(); - let team_diff_vars = (0..group_size - 1).map(|_| Variable::new()).collect::>(); + let tau_sqr = TAU.powi(2); + let beta_sqr = BETA.powi(2); - let team_sizes = _team_sizes(rating_groups); -} + let mut variable_arena = VariableArena::new(); -fn rate(rating_groups: &[&[Rating]]) { - let ranks = (0..rating_groups.len()).collect::>(); - let weights = vec![1.0; rating_groups.len()]; - let min_delta = DELTA; + let mut ss = Vec::new(); + let mut ps = Vec::new(); + let mut ts = Vec::new(); + let mut ds = Vec::new(); + + for _ in 0..size { + ss.push(variable_arena.create()); + ps.push(variable_arena.create()); + ts.push(variable_arena.create()); + } + + for _ in 0..size - 1 { + ds.push(variable_arena.create()); + } + + let mut factor_id = 0; + + let mut skill = Vec::new(); + + for (i, rating) in flatten_ratings.iter().enumerate() { + let variable = ss[i]; + let gaussian = Gaussian::with_mu_sigma(rating.mu, (rating.sigma.powi(2) + tau_sqr).sqrt()); + + skill.push(PriorFactor::new( + &mut variable_arena, + factor_id, + variable, + gaussian, + )); + + factor_id += 1; + } + + let mut skill_to_perf = Vec::new();; + + for i in 0..size { + let mean = ss[i]; + let value = ps[i]; + + skill_to_perf.push(LikelihoodFactor::new( + &mut variable_arena, + factor_id, + mean, + value, + beta_sqr, + )); + + factor_id += 1; + } + + let mut perf_to_team = Vec::new();; + + for i in 0..size { + perf_to_team.push(SumFactor::new( + &mut variable_arena, + factor_id, + ts[i], + vec![ps[i]], + vec![1.0], + )); + + factor_id += 1; + } + + let mut team_diff = Vec::new();; + + for i in 0..size - 1 { + team_diff.push(SumFactor::new( + &mut variable_arena, + factor_id, + ds[i], + vec![ts[i], ts[i + 1]], + vec![1.0, -1.0], + )); + + factor_id += 1; + } + + let mut trunc = Vec::new(); + + for i in 0..size - 1 { + trunc.push(TruncateFactor::new( + &mut variable_arena, + factor_id, + ds[i], + DELTA, + false, + )); + + factor_id += 1; + } + + for factor in &skill { + factor.start(&mut variable_arena); + } + + for factor in &skill_to_perf { + factor.update_value(&mut variable_arena); + } + + for factor in &perf_to_team { + factor.update_sum(&mut variable_arena); + } + + for factor in &team_diff { + factor.update_sum(&mut variable_arena); + } + + for factor in &trunc { + factor.update(&mut variable_arena); + } + + for factor in &team_diff { + factor.update_term(&mut variable_arena, 0); + factor.update_term(&mut variable_arena, 1); + } + + for factor in &perf_to_team { + factor.update_term(&mut variable_arena, 0); + } + + for factor in &skill_to_perf { + factor.update_mean(&mut variable_arena); + } + + for i in 0..size { + let value = variable_arena + .get(ss[i]) + .map(|variable| variable.get_value()) + .unwrap(); + + let mu = value.mu(); + let sigma = value.sigma(); + + println!("player={}, mu={}, sigma={}", i, mu, sigma); + } } fn quality(rating_groups: &[&[Rating]]) -> f32 { @@ -151,4 +266,21 @@ mod tests { assert_eq!(quality(&[&[alice], &[bob]]), 0.41614607); } + + #[test] + fn test_rate_1vs1() { + let alice = Rating { + mu: 25.0, + sigma: SIGMA, + }; + + let bob = Rating { + mu: 30.0, + sigma: SIGMA, + }; + + rate(&[&[alice], &[bob]]); + + assert_eq!(true, false); + } } diff --git a/src/math.rs b/src/math.rs new file mode 100644 index 0000000..d96a404 --- /dev/null +++ b/src/math.rs @@ -0,0 +1,30 @@ +use std::f32; + +fn erfc(x: f32) -> f32 { + let z = x.abs(); + let t = 1.0 / (1.0 + z / 2.0); + let r = t + * (-z * z - 1.26551223 + + t * (1.00002368 + + t * (0.37409196 + + t * (0.09678418 + + t * (-0.18628806 + + t * (0.27886807 + + t * (-1.13520398 + + t * (1.48851587 + t * (-0.82215223 + t * 0.17087277))))))))) + .exp(); + + if x < 0.0 { + 2.0 - r + } else { + r + } +} + +pub fn cdf(x: f32) -> f32 { + 0.5 * erfc(-x / 2.0f32.sqrt()) +} + +pub fn pdf(x: f32) -> f32 { + 1.0 / (2.0 * f32::consts::PI).sqrt() * (-((x / 1.0).powi(2) / 2.0)).exp() +}