Revert back.

Cargo.toml

@@ -4,3 +4,5 @@ version = "0.1.0"
 authors = ["Anders Olsson <anders.e.olsson@gmail.com>"]
 
 [dependencies]
+statrs = "0.10"
+noisy_float = "0.1"

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_arena.get_mut(variable).unwrap().attach_factor(id);
-            variable.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) {
+        println!("PriorFactor: variable.id={}, msg={:?}", self.variable.index, self.gaussian);
-            variable.update_value(self.id, 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_arena.get_mut(mean).unwrap().attach_factor(id);
-            variable.attach_factor(id);
+        variable_arena.get_mut(value).unwrap().attach_factor(id);
-        }
-
-        if let Some(variable) = variable_arena.get_mut(value) {
-            variable.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) {
+        println!("LikelihoodFactor: mean.id={}, msg={:?}", self.mean.index, gaussian);
-            variable.update_message(self.id, 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<VariableId>,
         coeffs: Vec<f64>,
     ) -> SumFactor {
-        if let Some(variable) = variable_arena.get_mut(sum) {
+        variable_arena.get_mut(sum).unwrap().attach_factor(id);
-            variable.attach_factor(id);
-        }
 
         for term in &terms {
-            if let Some(variable) = variable_arena.get_mut(*term) {
+            variable_arena.get_mut(*term).unwrap().attach_factor(id);
-                variable.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_pi += da.powi(2) / (gy.pi() - gfy.pi());
-            sum_tau += da * (gy.tau - gfy.tau) / (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) {
+        if variable == self.sum {
-            variable.update_value(self.id, gaussian);
+            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_arena.get_mut(variable).unwrap().attach_factor(id);
-            variable.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 c = x.pi() - fx.pi();
-        let d = x.tau - fx.tau;
+        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_arena.get_mut(self.variable).unwrap().update_value(self.id, gaussian);
-            variable.update_value(self.id, gaussian);
-        }
     }
 }

src/gaussian.rs

@@ -1,29 +1,33 @@
 use std::ops;
 
-#[derive(Clone, Copy)]
+#[derive(Clone, Copy, Debug)]
 pub struct Gaussian {
-    pub pi: f64,
+    pi: f64,
-    pub tau: f64,
+    tau: f64,
 }
 
 impl Gaussian {
-    pub fn new() -> Gaussian {
+    pub fn from_pi_tau(pi: f64, tau: f64) -> Gaussian {
-        Gaussian::with_pi_tau(0.0, 0.0)
-    }
-
-    pub fn with_pi_tau(pi: f64, tau: f64) -> Gaussian {
         Gaussian { pi, tau }
     }
 
-    pub fn with_mu_sigma(mu: f64, sigma: f64) -> Gaussian {
+    pub fn from_mu_sigma(mu: f64, sigma: f64) -> Gaussian {
-        let pi = 1.0 / sigma.powi(2);
+        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<Gaussian> 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<Gaussian> 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)
     }
 }

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 fn rate(rating_groups: &[&[Gaussian]]) {
-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]]) {
     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,
+            ss[i],
-            value,
+            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 mu = value.mean();
-        let sigma = value.sigma();
+        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 mu = value.mean();
-        let sigma = value.sigma();
+        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 {
+        let alice = Gaussian::from_mean_std_dev(MU, SIGMA);
-            mu: 25.0,
+        let bob = Gaussian::from_mean_std_dev(MU, SIGMA);
-            sigma: SIGMA,
-        };
 
-        let bob = Rating {
+        assert_eq!(quality(&[&[alice], &[bob]]), 0.4472135954999579);
-            mu: 30.0,
-            sigma: SIGMA,
-        };
-
-        assert_eq!(quality(&[&[alice], &[bob]]), 0.41614606763952605);
     }
 
     #[test]
     fn test_rate_1vs1() {
-        let alice = Rating {
+        let alice = Gaussian::from_mean_std_dev(MU, SIGMA);
-            mu: 25.0,
+        let bob = Gaussian::from_mean_std_dev(MU, SIGMA);
-            sigma: SIGMA,
+        let chris = Gaussian::from_mean_std_dev(MU, SIGMA);
-        };
+        let darren = Gaussian::from_mean_std_dev(MU, SIGMA);
 
-        let bob = Rating {
+        // println!("alice: {:?}", alice);
-            mu: 30.0,
+        // println!("bob: {:?}", alice);
-            sigma: SIGMA,
-        };
 
-        // alice: mu=27.175, sigma=7.112
+        rate(&[&[alice], &[bob], &[chris], &[darren]]);
-        // bob:   mu=23.141, sigma=6.859
-
-        rate(&[&[alice], &[bob]]);
 
         assert_eq!(true, false);
     }

src/math.rs

@@ -1,30 +1,20 @@
-use std::f64;
+use statrs::distribution::{Normal, Univariate, Continuous};
-
-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
-    }
-}
 
 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);
+    }
 }