Initial commit.

.gitignore

@@ -0,0 +1,2 @@
+/target
+/Cargo.lock

Cargo.toml

@@ -0,0 +1,8 @@
+[package]
+name = "trueskill-tt"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]

README.md

@@ -0,0 +1,3 @@
+# TrueSkill - Through Time
+
+Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py).

src/game.rs

@@ -0,0 +1,540 @@
+use std::cmp::Reverse;
+use std::collections::HashSet;
+
+use crate::{message::DiffMessages, utils, variable::TeamVariable, Gaussian, Player, N00};
+
+pub struct Game {
+    teams: Vec<Vec<Player>>,
+    result: Vec<u16>,
+    p_draw: f64,
+    pub likelihoods: Vec<Vec<Gaussian>>,
+    pub evidence: f64,
+}
+
+impl Game {
+    pub fn new(teams: Vec<Vec<Player>>, result: Vec<u16>, p_draw: f64) -> Self {
+        if !result.is_empty() {
+            assert!(
+                teams.len() == result.len(),
+                "len(result) and (len(teams) != len(result))"
+            );
+        }
+
+        assert!(p_draw >= 0.0 && p_draw < 1.0, "0.0 <= p_draw < 1.0");
+
+        if p_draw == 0.0 {
+            assert!(
+                result.iter().collect::<HashSet<_>>().len() == result.len(),
+                "(p_draw == 0.0) and (len(result) > 0) and (len(set(result)) != len(result))"
+            );
+        }
+
+        let mut this = Self {
+            teams,
+            result,
+            p_draw,
+            likelihoods: Vec::new(),
+            evidence: 0.0,
+        };
+
+        this.compute_likelihoods();
+
+        this
+    }
+
+    fn performance(&self, index: usize) -> Gaussian {
+        self.teams[index]
+            .iter()
+            .fold(N00, |sum, p| sum + p.performance())
+    }
+
+    fn partial_evidence(&mut self, d: &[DiffMessages], margin: &[f64], tie: &[bool], e: usize) {
+        let mu = d[e].prior.mu();
+        let sigma = d[e].prior.sigma();
+
+        if tie[e] {
+            self.evidence *= utils::cdf(margin[e], mu, sigma) - utils::cdf(-margin[e], mu, sigma)
+        } else {
+            self.evidence *= 1.0 - utils::cdf(margin[e], mu, sigma);
+        }
+    }
+
+    fn graphical_model(
+        &mut self,
+    ) -> (
+        Vec<usize>,
+        Vec<TeamVariable>,
+        Vec<DiffMessages>,
+        Vec<bool>,
+        Vec<f64>,
+    ) {
+        if self.result.is_empty() {
+            self.result = (0..self.teams.len() as u16).rev().collect::<Vec<_>>();
+        }
+
+        let r = &self.result;
+        let o = sortperm(r);
+
+        let t = (0..self.teams.len())
+            .map(|e| TeamVariable {
+                prior: self.teams[o[e]]
+                    .iter()
+                    .fold(N00, |sum, p| sum + p.performance()),
+                ..Default::default()
+            })
+            .collect::<Vec<_>>();
+
+        let d = t
+            .windows(2)
+            .map(|window| DiffMessages {
+                prior: window[0].prior - window[1].prior,
+                ..Default::default()
+            })
+            .collect::<Vec<_>>();
+
+        let tie = (0..d.len())
+            .map(|e| r[o[e]] == r[o[e + 1]])
+            .collect::<Vec<_>>();
+
+        let margin = (0..d.len())
+            .map(|e| {
+                if self.p_draw == 0.0 {
+                    0.0
+                } else {
+                    let a: f64 = self.teams[o[e]].iter().map(|a| a.beta.powi(2)).sum();
+                    let b: f64 = self.teams[o[e + 1]].iter().map(|a| a.beta.powi(2)).sum();
+
+                    utils::compute_margin(self.p_draw, (a + b).sqrt())
+                }
+            })
+            .collect::<Vec<_>>();
+
+        self.evidence = 1.0;
+
+        (o, t, d, tie, margin)
+    }
+
+    fn likelihood_analitico(&mut self) -> Vec<Vec<Gaussian>> {
+        let (o, t, d, tie, margin) = self.graphical_model();
+
+        self.partial_evidence(&d, &margin, &tie, 0);
+
+        let d = d[0].prior;
+        let (mu_trunc, sigma_trunc) = utils::trunc(d.mu(), d.sigma(), margin[0], tie[0]);
+
+        let (delta_div, theta_div_pow2) = if d.sigma() == sigma_trunc {
+            (
+                d.sigma().powi(2) * mu_trunc - sigma_trunc.powi(2) * d.mu(),
+                f64::INFINITY,
+            )
+        } else {
+            (
+                (d.sigma().powi(2) * mu_trunc - sigma_trunc.powi(2) * d.mu())
+                    / (d.sigma().powi(2) - sigma_trunc.powi(2)),
+                (sigma_trunc.powi(2) * d.sigma().powi(2))
+                    / (d.sigma().powi(2) - sigma_trunc.powi(2)),
+            )
+        };
+
+        let mut res = Vec::new();
+
+        for i in 0..t.len() {
+            let mut team = Vec::new();
+
+            for j in 0..self.teams[o[i]].len() {
+                //
+                let mu = if d.sigma() == sigma_trunc {
+                    0.0
+                } else {
+                    self.teams[o[i]][j].prior.mu()
+                        + (delta_div - d.mu()) * (-1.0f64).powi(if i == 1 { 1 } else { 0 })
+                };
+
+                let sigma_analitico = (theta_div_pow2 + d.sigma().powi(2)
+                    - self.teams[o[i]][j].prior.sigma().powi(2))
+                .sqrt();
+
+                team.push(Gaussian::new(mu, sigma_analitico));
+            }
+
+            res.push(team);
+        }
+
+        if o[0] >= o[1] {
+            res.swap(0, 1);
+        }
+
+        res
+    }
+
+    fn likelihood_teams(&mut self) -> Vec<Gaussian> {
+        let (o, mut t, mut d, tie, margin) = self.graphical_model();
+
+        let mut step = (f64::INFINITY, f64::INFINITY);
+        let mut i = 0;
+
+        while ((step.0 > 1e-6) || (step.1 > 1e-6)) && i < 10 {
+            step = (0.0, 0.0);
+
+            for e in 0..d.len() - 1 {
+                d[e].prior = t[e].posterior_win() - t[e + 1].posterior_lose();
+
+                if i == 0 {
+                    let mu = d[e].prior.mu();
+                    let sigma = d[e].prior.sigma();
+
+                    if tie[e] {
+                        self.evidence *=
+                            utils::cdf(margin[e], mu, sigma) - utils::cdf(-margin[e], mu, sigma)
+                    } else {
+                        self.evidence *= 1.0 - utils::cdf(margin[e], mu, sigma);
+                    }
+                }
+
+                d[e].likelihood = utils::approx(d[e].prior, margin[e], tie[e]) / d[e].prior;
+
+                let likelihood_lose = t[e].posterior_win() - d[e].likelihood;
+                let delta = t[e + 1].likelihood_lose.delta(likelihood_lose);
+
+                step = (
+                    if step.0 > delta.0 { step.0 } else { delta.0 },
+                    if step.1 > delta.1 { step.1 } else { delta.1 },
+                );
+
+                t[e + 1].likelihood_lose = likelihood_lose;
+            }
+
+            for e in (1..d.len()).rev() {
+                d[e].prior = t[e].posterior_win() - t[e + 1].posterior_lose();
+
+                if i == 0 && e == d.len() - 1 {
+                    self.partial_evidence(&d, &margin, &tie, e);
+                }
+
+                d[e].likelihood = utils::approx(d[e].prior, margin[e], tie[e]) / d[e].prior;
+
+                let likelihood_win = t[e + 1].posterior_lose() + d[e].likelihood;
+                let delta = t[e].likelihood_win.delta(likelihood_win);
+
+                step = (
+                    if step.0 > delta.0 { step.0 } else { delta.0 },
+                    if step.1 > delta.1 { step.1 } else { delta.1 },
+                );
+
+                t[e].likelihood_win = likelihood_win;
+            }
+
+            i += 1;
+        }
+
+        if d.len() == 1 {
+            self.partial_evidence(&d, &margin, &tie, 0);
+
+            d[0].prior = t[0].posterior_win() - t[1].posterior_lose();
+            d[0].likelihood = utils::approx(d[0].prior, margin[0], tie[0]) / d[0].prior;
+        }
+
+        let t_e = t.len();
+        let d_e = d.len();
+
+        t[0].likelihood_win = t[1].posterior_lose() + d[0].likelihood;
+        t[t_e - 1].likelihood_lose = t[t_e - 2].posterior_win() - d[d_e - 1].likelihood;
+
+        (0..t.len())
+            .map(|e| t[o[e]].likelihood())
+            .collect::<Vec<_>>()
+    }
+
+    fn compute_likelihoods(&mut self) {
+        if self.teams.len() > 2 {
+            let m_t_ft = self.likelihood_teams();
+
+            self.likelihoods = (0..self.teams.len())
+                .map(|e| {
+                    (0..self.teams[e].len())
+                        .map(|i| m_t_ft[e] - self.performance(e).exclude(self.teams[e][i].prior))
+                        .collect::<Vec<_>>()
+                })
+                .collect::<Vec<_>>();
+        } else {
+            self.likelihoods = self.likelihood_analitico();
+        }
+    }
+
+    pub fn posteriors(&self) -> Vec<Vec<Gaussian>> {
+        (0..self.teams.len())
+            .map(|e| {
+                (0..self.teams[e].len())
+                    .map(|i| self.likelihoods[e][i] * self.teams[e][i].prior)
+                    .collect::<Vec<_>>()
+            })
+            .collect::<Vec<_>>()
+    }
+}
+
+fn sortperm(xs: &[u16]) -> Vec<usize> {
+    let mut x = xs.iter().enumerate().collect::<Vec<_>>();
+    x.sort_unstable_by_key(|(_, x)| Reverse(*x));
+    x.into_iter().map(|(i, _)| i).collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{Gaussian, Player, GAMMA, N_INF};
+
+    use super::*;
+
+    #[test]
+    fn test_sortperm() {
+        assert_eq!(sortperm(&[0, 1, 2, 0]), vec![2, 1, 0, 3]);
+    }
+
+    #[test]
+    fn test_1vs1() {
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+
+        assert_eq!(a.mu(), 20.79477925612302);
+        assert_eq!(b.mu(), 29.205220743876975);
+        assert_eq!(a.sigma(), 7.194481422570443);
+
+        let t_a = Player::new(Gaussian::new(29.0, 1.0), 25.0 / 6.0, GAMMA, N_INF);
+        let t_b = Player::new(Gaussian::new(25.0, 25.0 / 3.0), 25.0 / 6.0, GAMMA, N_INF);
+
+        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+
+        assert_eq!(a.mu(), 28.896475351225412);
+        assert_eq!(a.sigma(), 0.9966043313004235);
+        assert_eq!(b.mu(), 32.18921172045737);
+        assert_eq!(b.sigma(), 6.062063735879715);
+
+        let t_a = Player::new(Gaussian::new(1.139, 0.531), 1.0, 0.2125, N_INF);
+        let t_b = Player::new(Gaussian::new(15.568, 0.51), 1.0, 0.2125, N_INF);
+
+        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0);
+
+        assert_eq!(g.likelihoods[0][0].sigma(), f64::INFINITY);
+        assert_eq!(g.likelihoods[1][0].sigma(), f64::INFINITY);
+        assert_eq!(g.likelihoods[0][0].mu(), 0.0);
+    }
+
+    #[test]
+    fn test_1vs1vs1() {
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_c = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![1, 2, 0], 0.0);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+        let c = p[2][0];
+
+        assert_eq!(a.mu(), 25.00000000000592);
+        assert_eq!(a.sigma(), 6.238469796269066);
+        assert_eq!(b.mu(), 31.31135822129149);
+        assert_eq!(b.sigma(), 6.69881865477675);
+        assert_eq!(c.mu(), 18.688641778702593);
+        assert_eq!(c.sigma(), 6.698818654778007);
+
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_c = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![2, 1, 0], 0.0);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+        let c = p[2][0];
+
+        assert_eq!(a.mu(), 31.31135822129149);
+        assert_eq!(a.sigma(), 6.69881865477675);
+        assert_eq!(b.mu(), 25.00000000000592);
+        assert_eq!(b.sigma(), 6.238469796269066);
+        assert_eq!(c.mu(), 18.688641778702593);
+        assert_eq!(c.sigma(), 6.698818654778007);
+
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_c = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![1, 2, 0], 0.5);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+        let c = p[2][0];
+
+        assert_eq!(a.mu(), 24.999999999511545);
+        assert_eq!(a.sigma(), 6.092561128305945);
+        assert_eq!(b.mu(), 33.37931495595287);
+        assert_eq!(b.sigma(), 6.483575782278924);
+        assert_eq!(c.mu(), 16.62068504453558);
+        assert_eq!(c.sigma(), 6.483575782198122);
+    }
+
+    #[test]
+    fn test_1vs1_draw() {
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 0], 0.25);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+
+        assert_eq!(a.mu(), 25.0);
+        assert_eq!(a.sigma(), 6.469480769842277);
+        assert_eq!(b.mu(), 25.0);
+        assert_eq!(b.sigma(), 6.469480769842277);
+
+        let t_a = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF);
+
+        let t_b = Player::new(Gaussian::new(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0, N_INF);
+
+        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 0], 0.25);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+
+        assert_eq!(a.mu(), 25.736001810566616);
+        assert_eq!(a.sigma(), 2.709956162204711);
+        assert_eq!(b.mu(), 28.67288808419261);
+        assert_eq!(b.sigma(), 1.9164711604544398);
+    }
+
+    #[test]
+    fn test_1vs1vs1_draw() {
+        let t_a = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_b = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+        let t_c = Player::new(
+            Gaussian::new(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            25.0 / 300.0,
+            N_INF,
+        );
+
+        let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![0, 0, 0], 0.25);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+        let c = p[2][0];
+
+        assert_eq!(a.mu(), 24.999999999999996);
+        assert_eq!(a.sigma(), 5.729068664890827);
+        assert_eq!(b.mu(), 25.000000000000004);
+        assert_eq!(b.sigma(), 5.707423522433266);
+        assert_eq!(c.mu(), 24.999999999999996);
+        assert_eq!(c.sigma(), 5.729068664890825);
+
+        let t_a = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF);
+        let t_b = Player::new(Gaussian::new(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0, N_INF);
+        let t_c = Player::new(Gaussian::new(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0, N_INF);
+
+        let g = Game::new(vec![vec![t_a], vec![t_b], vec![t_c]], vec![0, 0, 0], 0.25);
+        let p = g.posteriors();
+
+        let a = p[0][0];
+        let b = p[1][0];
+        let c = p[2][0];
+
+        assert_eq!(a.mu(), 25.48850755025261);
+        assert_eq!(a.sigma(), 2.638208444298423);
+        assert_eq!(b.mu(), 25.51067170990121);
+        assert_eq!(b.sigma(), 2.6287517663583633);
+        assert_eq!(c.mu(), 28.555920328820523);
+        assert_eq!(c.sigma(), 1.8856891308577184);
+    }
+}

src/gaussian.rs

@@ -0,0 +1,220 @@
+use std::fmt;
+use std::ops;
+
+use crate::utils;
+
+pub const BETA: f64 = 1.0;
+pub const MU: f64 = 0.0;
+pub const SIGMA: f64 = BETA * 6.0;
+pub const GAMMA: f64 = BETA * 0.03;
+
+pub const N01: Gaussian = Gaussian::new(0.0, 1.0);
+pub const N00: Gaussian = Gaussian::new(0.0, 0.0);
+pub const N_INF: Gaussian = Gaussian::new(0.0, f64::INFINITY);
+pub const N_MS: Gaussian = Gaussian::new(MU, SIGMA);
+
+#[derive(Clone, Copy, PartialEq, Debug)]
+pub struct Gaussian {
+    mu: f64,
+    sigma: f64,
+}
+
+impl Gaussian {
+    pub const fn new(mu: f64, sigma: f64) -> Self {
+        Gaussian { mu, sigma }
+    }
+
+    pub fn mu(&self) -> f64 {
+        self.mu
+    }
+
+    pub fn sigma(&self) -> f64 {
+        self.sigma
+    }
+
+    pub fn tau(&self) -> f64 {
+        self.mu * self.pi()
+    }
+
+    pub fn pi(&self) -> f64 {
+        self.sigma.powi(-2)
+    }
+
+    pub fn forget(&self, gamma: f64, t: u32) -> Self {
+        Self::new(
+            self.mu,
+            (self.sigma().powi(2) + t as f64 * gamma.powi(2)).sqrt(),
+        )
+    }
+
+    pub fn delta(&self, m: Gaussian) -> (f64, f64) {
+        ((self.mu() - m.mu()).abs(), (self.sigma() - m.sigma()).abs())
+    }
+
+    pub fn exclude(&self, m: Gaussian) -> Gaussian {
+        Gaussian::new(
+            self.mu() - m.mu(),
+            (self.sigma().powi(2) - m.sigma().powi(2)).sqrt(),
+        )
+    }
+
+    /*
+    def forget(self,gamma,t):
+        return Gaussian(self.mu, math.sqrt(self.sigma**2 + t*gamma**2))
+
+    def delta(self, M):
+        return abs(self.mu - M.mu) , abs(self.sigma - M.sigma)
+
+    def exclude(self, M):
+        return Gaussian(self.mu - M.mu, math.sqrt(self.sigma**2 - M.sigma**2) )
+
+    def isapprox(self, M, tol=1e-4):
+        return (abs(self.mu - M.mu) < tol) and (abs(self.sigma - M.sigma) < tol)
+    */
+}
+
+impl Default for Gaussian {
+    fn default() -> Self {
+        Gaussian {
+            mu: MU,
+            sigma: SIGMA,
+        }
+    }
+}
+
+impl fmt::Display for Gaussian {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(f, "N(mu={:.3}, sigma={:.3})", self.mu, self.sigma)
+    }
+}
+
+impl ops::Add<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn add(self, rhs: Gaussian) -> Self::Output {
+        Gaussian {
+            mu: self.mu + rhs.mu,
+            sigma: (self.sigma.powi(2) + rhs.sigma.powi(2)).sqrt(),
+        }
+    }
+}
+
+impl ops::Sub<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn sub(self, rhs: Gaussian) -> Self::Output {
+        Gaussian {
+            mu: self.mu - rhs.mu,
+            sigma: (self.sigma.powi(2) + rhs.sigma.powi(2)).sqrt(),
+        }
+    }
+}
+
+impl ops::Mul<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn mul(self, rhs: Gaussian) -> Self::Output {
+        let (mu, sigma) = utils::mu_sigma(self.tau() + rhs.tau(), self.pi() + rhs.pi());
+
+        Gaussian { mu, sigma }
+    }
+}
+
+impl ops::Div<Gaussian> for Gaussian {
+    type Output = Gaussian;
+
+    fn div(self, rhs: Gaussian) -> Self::Output {
+        let (mu, sigma) = utils::mu_sigma(self.tau() - rhs.tau(), self.pi() - rhs.pi());
+
+        Gaussian { mu, sigma }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_add() {
+        let n = Gaussian {
+            mu: 25.0,
+            sigma: 25.0 / 3.0,
+        };
+
+        let m = Gaussian {
+            mu: 0.0,
+            sigma: 1.0,
+        };
+
+        assert_eq!(
+            n + m,
+            Gaussian {
+                mu: 25.0,
+                sigma: 8.393118874676116
+            }
+        );
+    }
+
+    #[test]
+    fn test_sub() {
+        let n = Gaussian {
+            mu: 25.0,
+            sigma: 25.0 / 3.0,
+        };
+
+        let m = Gaussian {
+            mu: 1.0,
+            sigma: 1.0,
+        };
+
+        assert_eq!(
+            n - m,
+            Gaussian {
+                mu: 24.0,
+                sigma: 8.393118874676116
+            }
+        );
+    }
+
+    #[test]
+    fn test_mul() {
+        let n = Gaussian {
+            mu: 25.0,
+            sigma: 25.0 / 3.0,
+        };
+
+        let m = Gaussian {
+            mu: 0.0,
+            sigma: 1.0,
+        };
+
+        assert_eq!(
+            n * m,
+            Gaussian {
+                mu: 0.35488958990536273,
+                sigma: 0.992876838486922
+            }
+        );
+    }
+
+    #[test]
+    fn test_div() {
+        let n = Gaussian {
+            mu: 25.0,
+            sigma: 25.0 / 3.0,
+        };
+
+        let m = Gaussian {
+            mu: 0.0,
+            sigma: 1.0,
+        };
+
+        assert_eq!(
+            m / n,
+            Gaussian {
+                mu: -0.3652597402597402,
+                sigma: 1.0072787050317253
+            }
+        );
+    }
+}

src/history.rs

@@ -0,0 +1,5 @@
+pub struct History {
+    mu: f64,
+    sigma: f64,
+    gamma: f64,
+}

src/lib.rs

@@ -0,0 +1,12 @@
+mod game;
+mod gaussian;
+mod history;
+mod message;
+mod player;
+mod utils;
+mod variable;
+
+pub use game::*;
+pub use gaussian::*;
+pub use history::*;
+pub use player::*;

src/main.rs

@@ -0,0 +1,9 @@
+use trueskill_tt::*;
+
+fn main() {
+    let t_a = Player::new(Gaussian::new(29.0, 1.0), 25.0 / 6.0, GAMMA, N_INF);
+    let t_b = Player::new(Gaussian::new(25.0, 25.0 / 3.0), 25.0 / 6.0, GAMMA, N_INF);
+
+    let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0, 1], 0.0);
+    let p = g.posteriors();
+}

src/message.rs

@@ -0,0 +1,22 @@
+use crate::{Gaussian, N_INF};
+
+#[derive(Debug)]
+pub struct DiffMessages {
+    pub prior: Gaussian,
+    pub likelihood: Gaussian,
+}
+
+impl DiffMessages {
+    pub fn p(&self) -> Gaussian {
+        self.prior * self.likelihood
+    }
+}
+
+impl Default for DiffMessages {
+    fn default() -> Self {
+        Self {
+            prior: N_INF,
+            likelihood: N_INF,
+        }
+    }
+}

src/player.rs

@@ -0,0 +1,47 @@
+use std::fmt;
+
+use crate::{Gaussian, BETA, GAMMA, N_INF};
+
+#[derive(Debug)]
+pub struct Player {
+    pub prior: Gaussian,
+    pub beta: f64,
+    gamma: f64,
+    prior_draw: Gaussian,
+}
+
+impl Player {
+    pub fn new(prior: Gaussian, beta: f64, gamma: f64, prior_draw: Gaussian) -> Self {
+        Player {
+            prior,
+            beta,
+            gamma,
+            prior_draw,
+        }
+    }
+}
+
+impl Player {
+    pub fn performance(&self) -> Gaussian {
+        Gaussian::new(
+            self.prior.mu(),
+            (self.prior.sigma().powi(2) + self.beta.powi(2)).sqrt(),
+        )
+    }
+}
+
+impl Default for Player {
+    fn default() -> Self {
+        Player::new(Gaussian::default(), BETA, GAMMA, N_INF)
+    }
+}
+
+impl fmt::Display for Player {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        write!(
+            f,
+            "Player({}, beta={:.3}, gamma={:.3})",
+            self.prior, self.beta, self.gamma
+        )
+    }
+}

src/utils.rs

@@ -0,0 +1,209 @@
+use std::f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2};
+
+use crate::Gaussian;
+
+const SQRT_TAU: f64 = 2.5066282746310002;
+
+fn erfc(x: f64) -> f64 {
+    let z = x.abs();
+    let t = 1.0 / (1.0 + z / 2.0);
+
+    let a = -0.82215223 + t * 0.17087277;
+    let b = 1.48851587 + t * a;
+    let c = -1.13520398 + t * b;
+    let d = 0.27886807 + t * c;
+    let e = -0.18628806 + t * d;
+    let f = 0.09678418 + t * e;
+    let g = 0.37409196 + t * f;
+    let h = 1.00002368 + t * g;
+
+    let r = t * (-z * z - 1.26551223 + t * h).exp();
+
+    if x >= 0.0 {
+        r
+    } else {
+        2.0 - r
+    }
+}
+
+fn erfc_inv(mut y: f64) -> f64 {
+    if y >= 2.0 {
+        return f64::NEG_INFINITY;
+    }
+
+    debug_assert!(y >= 0.0, "argument must be nonnegative");
+
+    if y == 0.0 {
+        return f64::INFINITY;
+    }
+
+    if y >= 1.0 {
+        y = 2.0 - y;
+    }
+
+    let t = (-2.0 * (y / 2.0).ln()).sqrt();
+
+    let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
+
+    for _ in 0..3 {
+        let err = erfc(x) - y;
+
+        x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err)
+    }
+
+    if y < 1.0 {
+        x
+    } else {
+        -x
+    }
+}
+
+fn ppf(p: f64, mu: f64, sigma: f64) -> f64 {
+    mu - sigma * SQRT_2 * erfc_inv(2.0 * p)
+}
+
+pub(crate) fn mu_sigma(tau: f64, pi: f64) -> (f64, f64) {
+    if pi > 0.0 {
+        return (tau / pi, (1.0 / pi).sqrt());
+    }
+
+    if pi + 1e-5 < 0.0 {
+        panic!("sigma should be greater than 0");
+    }
+
+    (0.0, f64::INFINITY)
+}
+
+pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let z = -(x - mu) / (sigma * SQRT_2);
+
+    0.5 * erfc(z)
+}
+
+fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let normalizer = (SQRT_TAU * sigma).powi(-1);
+    let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp();
+
+    normalizer * functional
+}
+
+/*
+def ppf(p, mu, sigma):
+    return mu - sigma * sqrt2  * erfcinv(2 * p)
+*/
+
+fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
+    if !tie {
+        let alpha = (margin - mu) / sigma;
+
+        let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
+        let w = v * (v + (-alpha));
+
+        (v, w)
+    } else {
+        let alpha = (-margin - mu) / sigma;
+        let beta = (margin - mu) / sigma;
+
+        let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0))
+            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
+        let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0))
+            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
+        let w = -(u - v.powi(2));
+
+        (v, w)
+    }
+}
+
+pub(crate) fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
+    let (v, w) = v_w(mu, sigma, margin, tie);
+
+    let mu_trunc = mu + sigma * v;
+    let sigma_trunc = sigma * (1.0 - w).sqrt();
+
+    (mu_trunc, sigma_trunc)
+}
+
+pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian {
+    let (mu, sigma) = trunc(n.mu(), n.sigma(), margin, tie);
+
+    Gaussian::new(mu, sigma)
+}
+
+pub(crate) fn compute_margin(p_draw: f64, sd: f64) -> f64 {
+    ppf(0.5 - p_draw / 2.0, 0.0, sd).abs()
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::{Gaussian, N01};
+
+    use super::*;
+
+    #[test]
+    fn test_ppf() {
+        assert_eq!(ppf(0.3, N01.mu(), N01.sigma()), -0.5244004458961101);
+
+        let n23 = Gaussian::new(2.0, 3.0);
+        assert_eq!(ppf(0.3, n23.mu(), n23.sigma()), 0.4267986623116695);
+    }
+
+    #[test]
+    fn test_cdf() {
+        assert_eq!(cdf(0.3, N01.mu(), N01.sigma()), 0.6179114097962345);
+
+        let n23 = Gaussian::new(2.0, 3.0);
+        assert_eq!(cdf(0.3, n23.mu(), n23.sigma()), 0.28547031198297773);
+    }
+
+    #[test]
+    fn test_pdf() {
+        assert_eq!(pdf(0.3, N01.mu(), N01.sigma()), 0.38138781546052414);
+
+        let n23 = Gaussian::new(2.0, 3.0);
+        assert_eq!(pdf(0.3, n23.mu(), n23.sigma()), 0.11325579143491937);
+    }
+
+    #[test]
+    fn test_compute_margin() {
+        assert_eq!(
+            compute_margin(0.25, 2.0f64.sqrt() * (25.0 / 6.0)),
+            1.8776005988640154
+        );
+        assert_eq!(
+            compute_margin(0.25, 3.0f64.sqrt() * (25.0 / 6.0)),
+            2.2995817039804787
+        );
+        assert_eq!(
+            compute_margin(0.0, 3.0f64.sqrt() * (25.0 / 6.0)),
+            2.71348758713328e-7
+        );
+        assert_eq!(
+            compute_margin(1.0, 3.0f64.sqrt() * (25.0 / 6.0)),
+            f64::INFINITY
+        );
+    }
+
+    #[test]
+    fn test_trunc() {
+        let g = Gaussian::new(0.0, 1.0);
+
+        assert_eq!(
+            trunc(g.mu(), g.sigma(), 0.0, false),
+            (0.7978845368663289, 0.6028103066716792)
+        );
+
+        let g = Gaussian::new(0.0, SQRT_2 * (25.0 / 6.0));
+
+        assert_eq!(
+            trunc(g.mu(), g.sigma(), 1.8776005988, true),
+            (0.0, 1.0767055018086311)
+        );
+
+        let g = Gaussian::new(12.0, SQRT_2 * (25.0 / 6.0));
+
+        assert_eq!(
+            trunc(g.mu(), g.sigma(), 1.8776005988, true),
+            (0.39009949143595435, 1.034397855300721)
+        );
+    }
+}

src/variable.rs

@@ -0,0 +1,38 @@
+use crate::{Gaussian, N_INF};
+
+#[derive(Debug)]
+pub struct TeamVariable {
+    pub prior: Gaussian,
+    pub likelihood_lose: Gaussian,
+    pub likelihood_win: Gaussian,
+    pub likelihood_draw: Gaussian,
+}
+
+impl TeamVariable {
+    pub fn p(&self) -> Gaussian {
+        self.prior * self.likelihood_lose * self.likelihood_win * self.likelihood_draw
+    }
+
+    pub fn posterior_win(&self) -> Gaussian {
+        self.prior * self.likelihood_lose * self.likelihood_draw
+    }
+
+    pub fn posterior_lose(&self) -> Gaussian {
+        self.prior * self.likelihood_win * self.likelihood_draw
+    }
+
+    pub fn likelihood(&self) -> Gaussian {
+        self.likelihood_win * self.likelihood_lose * self.likelihood_draw
+    }
+}
+
+impl Default for TeamVariable {
+    fn default() -> Self {
+        Self {
+            prior: N_INF,
+            likelihood_lose: N_INF,
+            likelihood_win: N_INF,
+            likelihood_draw: N_INF,
+        }
+    }
+}