From b39c446b37403704c5a9e15e00e7be46aedb6044 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 24 Oct 2018 07:11:38 +0200 Subject: [PATCH] Revert back. --- Cargo.toml | 2 + src/factor_graph.rs | 98 +++++++++++++++++++++++---------------------- src/gaussian.rs | 32 ++++++++------- src/lib.rs | 79 +++++++++++++----------------------- src/math.rs | 38 +++++++----------- 5 files changed, 112 insertions(+), 137 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index be29532..02b8bb4 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,3 +4,5 @@ version = "0.1.0" authors = ["Anders Olsson "] [dependencies] +statrs = "0.10" +noisy_float = "0.1" diff --git a/src/factor_graph.rs b/src/factor_graph.rs index c671ae8..73d8b85 100644 --- a/src/factor_graph.rs +++ b/src/factor_graph.rs @@ -3,7 +3,7 @@ use std::collections::HashMap; use gaussian::Gaussian; use math; -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug, PartialEq)] pub struct VariableId { index: usize, } @@ -23,7 +23,7 @@ impl VariableArena { let index = self.variables.len(); self.variables.push(Variable { - value: Gaussian::with_pi_tau(0.0, 0.0), + value: Gaussian::from_precision(0.0, 0.0), factors: HashMap::new(), }); @@ -46,16 +46,23 @@ pub struct Variable { impl Variable { pub fn attach_factor(&mut self, factor: usize) { - self.factors.insert(factor, Gaussian::new()); + self.factors.insert(factor, Gaussian::from_precision(0.0, 0.0)); } pub fn update_value(&mut self, factor: usize, value: Gaussian) { let old = self.factors[&factor]; let intermediate = value * old; - let value = intermediate / self.value; + let new = intermediate / self.value; - self.value = value; + /* + println!("update_value: old={}, value={:?}, new={:?}, self.value={:?}", + render_g(old), + value, + new, + self.value); + */ + self.value = new; } pub fn get_value(&self) -> Gaussian { @@ -70,6 +77,8 @@ impl Variable { self.value = value; + println!("update_message: old={:?} msg={:?}, new={:?}", old, message, value); + self.factors.insert(factor, message); } @@ -91,9 +100,7 @@ impl PriorFactor { variable: VariableId, gaussian: Gaussian, ) -> PriorFactor { - if let Some(variable) = variable_arena.get_mut(variable) { - variable.attach_factor(id); - } + variable_arena.get_mut(variable).unwrap().attach_factor(id); PriorFactor { id, @@ -103,9 +110,9 @@ impl PriorFactor { } 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); - } + println!("PriorFactor: variable.id={}, msg={:?}", self.variable.index, self.gaussian); + + variable_arena.get_mut(self.variable).unwrap().update_value(self.id, self.gaussian); } } @@ -124,13 +131,8 @@ impl LikelihoodFactor { value: VariableId, variance: f64, ) -> 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); - } + variable_arena.get_mut(mean).unwrap().attach_factor(id); + variable_arena.get_mut(value).unwrap().attach_factor(id); LikelihoodFactor { id, @@ -151,13 +153,13 @@ impl LikelihoodFactor { .map(|variable| variable.get_message(self.id)) .unwrap(); - let a = 1.0 / (1.0 + self.variance * (x.pi - fx.pi)); + let a = 1.0 / (1.0 + self.variance * (x.precision_mean() - fx.precision_mean())); - let gaussian = Gaussian::with_pi_tau(a * (x.pi - fx.pi), a * (x.tau - fx.tau)); + let gaussian = Gaussian::from_precision(a * (x.precision_mean() - fx.precision_mean()), a * (x.precision() - fx.precision())); - if let Some(variable) = variable_arena.get_mut(self.mean) { - variable.update_message(self.id, gaussian); - } + println!("LikelihoodFactor: mean.id={}, msg={:?}", self.mean.index, gaussian); + + variable_arena.get_mut(self.mean).unwrap().update_message(self.id, gaussian); } pub fn update_value(&self, variable_arena: &mut VariableArena) { @@ -165,18 +167,20 @@ impl LikelihoodFactor { .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 a = 1.0 / (1.0 + self.variance * (y.precision_mean() - fy.precision_mean())); - let gaussian = Gaussian::with_pi_tau(a * (y.pi - fy.pi), a * (y.tau - fy.tau)); + let gaussian = Gaussian::from_precision(a * (y.precision_mean() - fy.precision_mean()), a * (y.precision() - fy.precision())); - if let Some(variable) = variable_arena.get_mut(self.value) { - variable.update_message(self.id, gaussian); - } + + println!("LikelihoodFactor: value.id={}, msg={:?}", self.value.index, gaussian); + + variable_arena.get_mut(self.value).unwrap().update_message(self.id, gaussian); } } @@ -195,14 +199,10 @@ impl SumFactor { terms: Vec, coeffs: Vec, ) -> SumFactor { - if let Some(variable) = variable_arena.get_mut(sum) { - variable.attach_factor(id); - } + variable_arena.get_mut(sum).unwrap().attach_factor(id); for term in &terms { - if let Some(variable) = variable_arena.get_mut(*term) { - variable.attach_factor(id); - } + variable_arena.get_mut(*term).unwrap().attach_factor(id); } SumFactor { @@ -231,18 +231,22 @@ impl SumFactor { 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); + 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); + let gaussian = Gaussian::from_pi_tau(new_pi, new_tau); - if let Some(variable) = variable_arena.get_mut(variable) { - variable.update_value(self.id, gaussian); + 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); } pub fn update_sum(&self, variable_arena: &mut VariableArena) { @@ -351,9 +355,7 @@ impl TruncateFactor { epsilon: f64, draw: bool, ) -> TruncateFactor { - if let Some(variable) = variable_arena.get_mut(variable) { - variable.attach_factor(id); - } + variable_arena.get_mut(variable).unwrap().attach_factor(id); TruncateFactor { id, @@ -368,13 +370,15 @@ impl TruncateFactor { .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 c = x.pi() - fx.pi(); + let d = x.tau() - fx.tau(); + let sqrt_c = c.sqrt(); let t = d / sqrt_c; @@ -388,10 +392,8 @@ impl TruncateFactor { let m_w = 1.0 - w; - let gaussian = Gaussian::with_pi_tau(c / m_w, (d + sqrt_c * v) / m_w); + let gaussian = Gaussian::from_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); - } + variable_arena.get_mut(self.variable).unwrap().update_value(self.id, gaussian); } } diff --git a/src/gaussian.rs b/src/gaussian.rs index b046a7d..33166ac 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -1,29 +1,33 @@ use std::ops; -#[derive(Clone, Copy)] +#[derive(Clone, Copy, Debug)] pub struct Gaussian { - pub pi: f64, - pub tau: f64, + pi: f64, + tau: f64, } impl Gaussian { - pub fn new() -> Gaussian { - Gaussian::with_pi_tau(0.0, 0.0) - } - - pub fn with_pi_tau(pi: f64, tau: f64) -> Gaussian { + pub fn from_pi_tau(pi: f64, tau: f64) -> Gaussian { Gaussian { pi, tau } } - pub fn with_mu_sigma(mu: f64, sigma: f64) -> Gaussian { - let pi = 1.0 / sigma.powi(2); + pub fn from_mu_sigma(mu: f64, sigma: f64) -> Gaussian { + let pi = 1.0 / (sigma * sigma); - Gaussian::with_pi_tau(pi, pi * mu) + Self::from_pi_tau(pi, pi * mu) + } + + pub fn pi(&self) -> f64 { + self.pi + } + + pub fn tau(&self) -> f64 { + self.tau } pub fn mu(&self) -> f64 { if self.pi == 0.0 { - 0.0 + return 0.0; } else { self.tau / self.pi } @@ -38,7 +42,7 @@ 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) + Gaussian::from_pi_tau(self.pi + rhs.pi, self.tau + rhs.tau) } } @@ -46,6 +50,6 @@ 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) + Gaussian::from_pi_tau(self.pi - rhs.pi, self.tau - rhs.tau) } } diff --git a/src/lib.rs b/src/lib.rs index 47fa4b7..3d26f38 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,6 @@ +extern crate statrs; +extern crate noisy_float; + mod factor_graph; mod gaussian; mod math; @@ -25,22 +28,7 @@ const DRAW_PROBABILITY: f64 = 0.10; /// A basis to check reliability of the result. const DELTA: f64 = 0.0001; -#[derive(Debug, PartialEq)] -pub struct Rating { - pub mu: f64, - pub sigma: f64, -} - -impl Default for Rating { - fn default() -> Rating { - Rating { - mu: MU, - sigma: SIGMA, - } - } -} - -fn rate(rating_groups: &[&[Rating]]) { +pub fn rate(rating_groups: &[&[Gaussian]]) { let flatten_ratings = rating_groups .iter() .flat_map(|group| group.iter()) @@ -74,7 +62,7 @@ fn rate(rating_groups: &[&[Rating]]) { 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()); + let gaussian = Gaussian::from_mean_variance(rating.mean(), (rating.variance() + tau_sqr).sqrt()); skill.push(PriorFactor::new( &mut variable_arena, @@ -89,14 +77,11 @@ fn rate(rating_groups: &[&[Rating]]) { 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, + ss[i], + ps[i], beta_sqr, )); @@ -149,6 +134,7 @@ fn rate(rating_groups: &[&[Rating]]) { factor.start(&mut variable_arena); } + /* println!("before:"); for i in 0..size { @@ -157,11 +143,12 @@ fn rate(rating_groups: &[&[Rating]]) { .map(|variable| variable.get_value()) .unwrap(); - let mu = value.mu(); - let sigma = value.sigma(); + 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); @@ -194,6 +181,7 @@ fn rate(rating_groups: &[&[Rating]]) { factor.update_mean(&mut variable_arena); } + /* println!("after:"); for i in 0..size { @@ -202,14 +190,15 @@ fn rate(rating_groups: &[&[Rating]]) { .map(|variable| variable.get_value()) .unwrap(); - let mu = value.mu(); - let sigma = value.sigma(); + let mu = value.mean(); + let sigma = value.std_dev(); println!("* player={}, mu={}, sigma={}", i, mu, sigma); } + */ } -fn quality(rating_groups: &[&[Rating]]) -> f64 { +pub fn quality(rating_groups: &[&[Gaussian]]) -> f64 { let flatten_ratings = rating_groups .iter() .flat_map(|group| group.iter()) @@ -222,13 +211,13 @@ fn quality(rating_groups: &[&[Rating]]) -> f64 { let mut mean_matrix = Matrix::new(length, 1); for (i, rating) in flatten_ratings.iter().enumerate() { - mean_matrix[(i, 0)] = rating.mu; + mean_matrix[(i, 0)] = rating.mean(); } let mut variance_matrix = Matrix::new(length, length); for (i, rating) in flatten_ratings.iter().enumerate() { - variance_matrix[(i, i)] = rating.sigma.powi(2); + variance_matrix[(i, i)] = rating.std_dev().powi(2); } let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length); @@ -272,35 +261,23 @@ mod tests { #[test] fn test_quality_1vs1() { - let alice = Rating { - mu: 25.0, - sigma: SIGMA, - }; + let alice = Gaussian::from_mean_std_dev(MU, SIGMA); + let bob = Gaussian::from_mean_std_dev(MU, SIGMA); - let bob = Rating { - mu: 30.0, - sigma: SIGMA, - }; - - assert_eq!(quality(&[&[alice], &[bob]]), 0.41614606763952605); + assert_eq!(quality(&[&[alice], &[bob]]), 0.4472135954999579); } #[test] fn test_rate_1vs1() { - let alice = Rating { - mu: 25.0, - sigma: SIGMA, - }; + 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 bob = Rating { - mu: 30.0, - sigma: SIGMA, - }; + // println!("alice: {:?}", alice); + // println!("bob: {:?}", alice); - // alice: mu=27.175, sigma=7.112 - // bob: mu=23.141, sigma=6.859 - - rate(&[&[alice], &[bob]]); + rate(&[&[alice], &[bob], &[chris], &[darren]]); assert_eq!(true, false); } diff --git a/src/math.rs b/src/math.rs index 1986b8e..1701474 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,30 +1,20 @@ -use std::f64; - -fn erfc(x: f64) -> f64 { - 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 - } -} +use statrs::distribution::{Normal, Univariate, Continuous}; pub fn cdf(x: f64) -> f64 { - 0.5 * erfc(-x / 2.0f64.sqrt()) + Normal::new(0.0, 1.0).unwrap().cdf(x) } pub fn pdf(x: f64) -> f64 { - 1.0 / (2.0 * f64::consts::PI).sqrt() * (-((x / 1.0).powi(2) / 2.0)).exp() + Normal::new(0.0, 1.0).unwrap().pdf(x) +} + + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_cdf() { + assert_eq!(cdf(0.5), 0.6914624612740131); + } }