Ignore temp folder

.gitignore

@@ -1,5 +1,6 @@
 /target
 /Cargo.lock
 
+/temp
 .justfile
 *.svg

NOTEPAD.md

@@ -1,15 +1,8 @@
 # History
 
-```shell
+```rust
-teams: [[player]]
+let mut history = History::new();
-weights: [[f64]]
-results: [f64]
 
-player: (gaussian, f64, f64)
+let agent_a = history.new_agent();
-
+let agent_b = history.new_agent_with_prior(Prior::new(Gaussian::default(), BETA, GAMMA));
-players: [player]
-weights: [f64]
-
-
-teams: [([(player, weight)], result)]
 ```

README.md

@@ -15,6 +15,7 @@ Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillTh
 
 - [x] Implement approx for Gaussian
 - [x] Add more tests from `TrueSkillThroughTime.jl`
+- [ ] Add tests for `quality()` (Use [sublee/trueskill](https://github.com/sublee/trueskill/tree/master) as reference)
 - [ ] Benchmark Batch::iteration()
 - [ ] Time needs to be an enum so we can have multiple states (see `batch::compute_elapsed()`)
 - [ ] Add examples (use same TrueSkillThroughTime.(py|jl))

benches/batch.rs

@@ -51,7 +51,8 @@ fn criterion_benchmark(criterion: &mut Criterion) {
         weights.push(vec![vec![1.0], vec![1.0]]);
     }
 
-    let mut batch = Batch::new(composition, results, weights, 1, P_DRAW, &agents);
+    let mut batch = Batch::new(1, P_DRAW);
+    batch.add_events(composition, results, weights, &agents);
 
     criterion.bench_function("Batch::iteration", |b| {
         b.iter(|| batch.iteration(0, &agents))

src/batch.rs

@@ -1,4 +1,4 @@
-use std::collections::{HashMap, HashSet};
+use std::collections::HashMap;
 
 use crate::{
     agent::Agent, game::Game, gaussian::Gaussian, player::Player, tuple_gt, tuple_max, Index, N_INF,
@@ -107,119 +107,43 @@ pub struct Batch {
 }
 
 impl Batch {
-    pub fn new(
+    pub fn new(time: i64, p_draw: f64) -> Self {
-        composition: Vec<Vec<Vec<Index>>>,
+        Self {
-        results: Vec<Vec<f64>>,
+            events: Vec::new(),
-        weights: Vec<Vec<Vec<f64>>>,
+            skills: HashMap::new(),
-        time: i64,
-        p_draw: f64,
-        agents: &HashMap<Index, Agent>,
-    ) -> Self {
-        assert!(
-            results.is_empty() || results.len() == composition.len(),
-            "TODO: Add a comment here"
-        );
-        assert!(
-            weights.is_empty() || weights.len() == composition.len(),
-            "TODO: Add a comment here"
-        );
-
-        let this_agent = composition
-            .iter()
-            .flatten()
-            .flatten()
-            .collect::<HashSet<_>>();
-
-        let skills = this_agent
-            .iter()
-            .map(|&idx| {
-                let elapsed = compute_elapsed(agents[idx].last_time, time);
-
-                (
-                    *idx,
-                    Skill {
-                        forward: agents[idx].receive(elapsed),
-                        elapsed,
-                        ..Default::default()
-                    },
-                )
-            })
-            .collect::<HashMap<_, _>>();
-
-        let events = composition
-            .iter()
-            .enumerate()
-            .map(|(e, event)| {
-                let teams = event
-                    .iter()
-                    .enumerate()
-                    .map(|(t, team)| {
-                        let items = team
-                            .iter()
-                            .map(|&agent| Item {
-                                agent,
-                                likelihood: N_INF,
-                            })
-                            .collect::<Vec<_>>();
-
-                        Team {
-                            items,
-                            output: if results.is_empty() {
-                                (event.len() - (t + 1)) as f64
-                            } else {
-                                results[e][t]
-                            },
-                        }
-                    })
-                    .collect::<Vec<_>>();
-
-                Event {
-                    teams,
-                    evidence: 0.0,
-                    weights: if weights.is_empty() {
-                        Vec::new()
-                    } else {
-                        weights[e].clone()
-                    },
-                }
-            })
-            .collect::<Vec<_>>();
-
-        let mut this = Self {
             time,
-            events,
-            skills,
             p_draw,
-        };
+        }
-
-        this.iteration(0, agents);
-
-        this
     }
 
-    pub(crate) fn add_events(
+    pub fn add_events(
         &mut self,
         composition: Vec<Vec<Vec<Index>>>,
         results: Vec<Vec<f64>>,
         weights: Vec<Vec<Vec<f64>>>,
         agents: &HashMap<Index, Agent>,
     ) {
-        let this_agent = composition
+        let mut unique = Vec::with_capacity(10);
-            .iter()
+
-            .flatten()
+        let this_agent = composition.iter().flatten().flatten().filter(|idx| {
-            .flatten()
+            if !unique.contains(idx) {
-            .cloned()
+                unique.push(*idx);
-            .collect::<HashSet<_>>();
+
+                return true;
+            }
+
+            false
+        });
 
         for idx in this_agent {
             let elapsed = compute_elapsed(agents[&idx].last_time, self.time);
 
-            if let Some(skill) = self.skills.get_mut(&idx) {
+            if let Some(skill) = self.skills.get_mut(idx) {
                 skill.elapsed = elapsed;
                 skill.forward = agents[&idx].receive(elapsed);
             } else {
                 self.skills.insert(
-                    idx,
+                    *idx,
                     Skill {
                         forward: agents[&idx].receive(elapsed),
                         elapsed,
@@ -253,14 +177,19 @@ impl Batch {
                 })
                 .collect::<Vec<_>>();
 
+            let weights = if weights.is_empty() {
+                teams
+                    .iter()
+                    .map(|team| vec![1.0; team.items.len()])
+                    .collect::<Vec<_>>()
+            } else {
+                weights[e].clone()
+            };
+
             Event {
                 teams,
                 evidence: 0.0,
-                weights: if weights.is_empty() {
+                weights,
-                    Vec::new()
-                } else {
-                    weights[e].clone()
-                },
             }
         });
 
@@ -283,7 +212,7 @@ impl Batch {
             let teams = event.within_priors(false, false, &self.skills, agents);
             let result = event.outputs();
 
-            let g = Game::new(teams, result, event.weights.clone(), self.p_draw);
+            let g = Game::new(teams, &result, &event.weights, self.p_draw);
 
             for (t, team) in event.teams.iter_mut().enumerate() {
                 for (i, item) in team.items.iter_mut().enumerate() {
@@ -372,8 +301,8 @@ impl Batch {
                     .map(|(_, event)| {
                         Game::new(
                             event.within_priors(online, forward, &self.skills, agents),
-                            event.outputs(),
+                            &event.outputs(),
-                            event.weights.clone(),
+                            &event.weights,
                             self.p_draw,
                         )
                         .evidence
@@ -397,8 +326,8 @@ impl Batch {
                 .map(|(_, event)| {
                     Game::new(
                         event.within_priors(online, forward, &self.skills, agents),
-                        event.outputs(),
+                        &event.outputs(),
-                        event.weights.clone(),
+                        &event.weights,
                         self.p_draw,
                     )
                     .evidence
@@ -492,7 +421,9 @@ mod tests {
             );
         }
 
-        let mut batch = Batch::new(
+        let mut batch = Batch::new(0, 0.0);
+
+        batch.add_events(
             vec![
                 vec![vec![a], vec![b]],
                 vec![vec![c], vec![d]],
@@ -500,8 +431,6 @@ mod tests {
             ],
             vec![vec![1.0, 0.0], vec![0.0, 1.0], vec![1.0, 0.0]],
             vec![],
-            0,
-            0.0,
             &agents,
         );
 
@@ -568,7 +497,9 @@ mod tests {
             );
         }
 
-        let mut batch = Batch::new(
+        let mut batch = Batch::new(0, 0.0);
+
+        batch.add_events(
             vec![
                 vec![vec![a], vec![b]],
                 vec![vec![a], vec![c]],
@@ -576,8 +507,6 @@ mod tests {
             ],
             vec![vec![1.0, 0.0], vec![0.0, 1.0], vec![1.0, 0.0]],
             vec![],
-            0,
-            0.0,
             &agents,
         );
 
@@ -647,7 +576,9 @@ mod tests {
             );
         }
 
-        let mut batch = Batch::new(
+        let mut batch = Batch::new(0, 0.0);
+
+        batch.add_events(
             vec![
                 vec![vec![a], vec![b]],
                 vec![vec![a], vec![c]],
@@ -655,8 +586,6 @@ mod tests {
             ],
             vec![vec![1.0, 0.0], vec![0.0, 1.0], vec![1.0, 0.0]],
             vec![],
-            0,
-            0.0,
             &agents,
         );
 

src/game.rs

@@ -7,66 +7,49 @@ use crate::{
 };
 
 #[derive(Debug)]
-pub struct Game {
+pub struct Game<'a> {
     teams: Vec<Vec<Player>>,
-    result: Vec<f64>,
+    result: &'a [f64],
-    weights: Vec<Vec<f64>>,
+    weights: &'a [Vec<f64>],
     p_draw: f64,
     pub(crate) likelihoods: Vec<Vec<Gaussian>>,
     pub(crate) evidence: f64,
 }
 
-impl Game {
+impl<'a> Game<'a> {
     pub fn new(
         teams: Vec<Vec<Player>>,
-        mut result: Vec<f64>,
+        result: &'a [f64],
-        mut weights: Vec<Vec<f64>>,
+        weights: &'a [Vec<f64>],
         p_draw: f64,
     ) -> Self {
-        assert!(
+        debug_assert!(
-            (result.is_empty() || result.len() == teams.len()),
+            (result.len() == teams.len()),
-            "result must be empty or the same length as teams"
+            "result must have the same length as teams"
         );
 
-        assert!(
+        debug_assert!(
-            (weights.is_empty() || weights.len() == teams.len()),
+            weights
-            "weights must be empty or the same length as teams"
-        );
-
-        assert!(
-            weights.is_empty()
-                || weights
                 .iter()
                 .zip(teams.iter())
                 .all(|(w, t)| w.len() == t.len()),
-            "weights must be empty or has the same dimensions as teams"
+            "weights must have the same dimensions as teams"
         );
 
-        assert!(
+        debug_assert!(
             (0.0..1.0).contains(&p_draw),
             "draw probability.must be >= 0.0 and < 1.0"
         );
 
-        assert!(
+        debug_assert!(
             p_draw > 0.0 || {
-                let mut r = result.clone();
+                let mut r = result.to_vec();
                 r.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
                 r.windows(2).all(|w| w[0] != w[1])
             },
             "draw must be > 0.0 if there is teams with draw"
         );
 
-        if result.is_empty() {
-            result = (0..teams.len()).rev().map(|i| i as f64).collect::<Vec<_>>();
-        }
-
-        if weights.is_empty() {
-            weights = teams
-                .iter()
-                .map(|team| vec![1.0; team.len()])
-                .collect::<Vec<_>>();
-        }
-
         let mut this = Self {
             teams,
             result,
@@ -82,31 +65,7 @@ impl Game {
     }
 
     fn likelihoods(&mut self) {
-        let m_t_ft = self.likelihood_teams();
+        let o = sort_perm(self.result, true);
-
-        self.likelihoods = self
-            .teams
-            .iter()
-            .zip(self.weights.iter())
-            .zip(m_t_ft)
-            .map(|((p, w), m)| {
-                let performance = p.iter().zip(w.iter()).fold(N00, |p, (player, &weight)| {
-                    p + (player.performance() * weight)
-                });
-
-                p.iter()
-                    .zip(w.iter())
-                    .map(|(p, &w)| {
-                        ((m - performance.exclude(p.performance() * w)) * (1.0 / w))
-                            .forget(p.beta, 1)
-                    })
-                    .collect::<Vec<_>>()
-            })
-            .collect::<Vec<_>>();
-    }
-
-    fn likelihood_teams(&mut self) -> Vec<Gaussian> {
-        let o = sort_perm(&self.result, true);
 
         let mut team = o
             .iter()
@@ -143,14 +102,10 @@ impl Game {
         } else {
             o.windows(2)
                 .map(|w| {
-                    if self.p_draw == 0.0 {
-                        0.0
-                    } else {
                     let a: f64 = self.teams[w[0]].iter().map(|a| a.beta.powi(2)).sum();
                     let b: f64 = self.teams[w[1]].iter().map(|a| a.beta.powi(2)).sum();
 
                     compute_margin(self.p_draw, (a + b).sqrt())
-                    }
                 })
                 .collect::<Vec<_>>()
         };
@@ -205,7 +160,27 @@ impl Game {
         team[0].likelihood_win = team[1].posterior_lose() + diff[0].likelihood;
         team[t_end].likelihood_lose = team[t_end - 1].posterior_win() - diff[d_end].likelihood;
 
-        o.iter().map(|&e| team[e].likelihood()).collect::<Vec<_>>()
+        let m_t_ft = o.into_iter().map(|e| team[e].likelihood());
+
+        self.likelihoods = self
+            .teams
+            .iter()
+            .zip(self.weights.iter())
+            .zip(m_t_ft)
+            .map(|((p, w), m)| {
+                let performance = p.iter().zip(w.iter()).fold(N00, |p, (player, &weight)| {
+                    p + (player.performance() * weight)
+                });
+
+                p.iter()
+                    .zip(w.iter())
+                    .map(|(p, &w)| {
+                        ((m - performance.exclude(p.performance() * w)) * (1.0 / w))
+                            .forget(p.beta, 1)
+                    })
+                    .collect::<Vec<_>>()
+            })
+            .collect::<Vec<_>>();
     }
 
     pub fn posteriors(&self) -> Vec<Vec<Gaussian>> {
@@ -243,7 +218,8 @@ mod tests {
             25.0 / 300.0,
         );
 
-        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0.0, 1.0], vec![], 0.0);
+        let w = [vec![1.0], vec![1.0]];
+        let g = Game::new(vec![vec![t_a], vec![t_b]], &[0.0, 1.0], &w, 0.0);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -255,7 +231,8 @@ mod tests {
         let t_a = Player::new(Gaussian::from_ms(29.0, 1.0), 25.0 / 6.0, GAMMA);
         let t_b = Player::new(Gaussian::from_ms(25.0, 25.0 / 3.0), 25.0 / 6.0, GAMMA);
 
-        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0.0, 1.0], vec![], 0.0);
+        let w = [vec![1.0], vec![1.0]];
+        let g = Game::new(vec![vec![t_a], vec![t_b]], &[0.0, 1.0], &w, 0.0);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -267,7 +244,8 @@ mod tests {
         let t_a = Player::new(Gaussian::from_ms(1.139, 0.531), 1.0, 0.2125);
         let t_b = Player::new(Gaussian::from_ms(15.568, 0.51), 1.0, 0.2125);
 
-        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0.0, 1.0], vec![], 0.0);
+        let w = [vec![1.0], vec![1.0]];
+        let g = Game::new(vec![vec![t_a], vec![t_b]], &[0.0, 1.0], &w, 0.0);
 
         assert_eq!(g.likelihoods[0][0], N_INF);
         assert_eq!(g.likelihoods[1][0], N_INF);
@@ -293,7 +271,8 @@ mod tests {
             )],
         ];
 
-        let g = Game::new(teams.clone(), vec![1.0, 2.0, 0.0], vec![], 0.0);
+        let w = [vec![1.0], vec![1.0], vec![1.0]];
+        let g = Game::new(teams.clone(), &[1.0, 2.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -302,7 +281,8 @@ mod tests {
         assert_ulps_eq!(a, Gaussian::from_ms(25.000000, 6.238469), epsilon = 1e-6);
         assert_ulps_eq!(b, Gaussian::from_ms(31.311358, 6.698818), epsilon = 1e-6);
 
-        let g = Game::new(teams.clone(), vec![], vec![], 0.0);
+        let w = [vec![1.0], vec![1.0], vec![1.0]];
+        let g = Game::new(teams.clone(), &[2.0, 1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -311,7 +291,8 @@ mod tests {
         assert_ulps_eq!(a, Gaussian::from_ms(31.311358, 6.698818), epsilon = 1e-6);
         assert_ulps_eq!(b, Gaussian::from_ms(25.000000, 6.238469), epsilon = 1e-6);
 
-        let g = Game::new(teams, vec![1.0, 2.0, 0.0], vec![], 0.5);
+        let w = [vec![1.0], vec![1.0], vec![1.0]];
+        let g = Game::new(teams, &[1.0, 2.0, 0.0], &w, 0.5);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -336,7 +317,8 @@ mod tests {
             25.0 / 300.0,
         );
 
-        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0.0, 0.0], vec![], 0.25);
+        let w = [vec![1.0], vec![1.0]];
+        let g = Game::new(vec![vec![t_a], vec![t_b]], &[0.0, 0.0], &w, 0.25);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -348,7 +330,8 @@ mod tests {
         let t_a = Player::new(Gaussian::from_ms(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0);
         let t_b = Player::new(Gaussian::from_ms(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0);
 
-        let g = Game::new(vec![vec![t_a], vec![t_b]], vec![0.0, 0.0], vec![], 0.25);
+        let w = [vec![1.0], vec![1.0]];
+        let g = Game::new(vec![vec![t_a], vec![t_b]], &[0.0, 0.0], &w, 0.25);
         let p = g.posteriors();
 
         let a = p[0][0];
@@ -376,10 +359,11 @@ mod tests {
             25.0 / 300.0,
         );
 
+        let w = [vec![1.0], vec![1.0], vec![1.0]];
         let g = Game::new(
             vec![vec![t_a], vec![t_b], vec![t_c]],
-            vec![0.0, 0.0, 0.0],
+            &[0.0, 0.0, 0.0],
-            vec![],
+            &w,
             0.25,
         );
         let p = g.posteriors();
@@ -396,10 +380,11 @@ mod tests {
         let t_b = Player::new(Gaussian::from_ms(25.0, 3.0), 25.0 / 6.0, 25.0 / 300.0);
         let t_c = Player::new(Gaussian::from_ms(29.0, 2.0), 25.0 / 6.0, 25.0 / 300.0);
 
+        let w = [vec![1.0], vec![1.0], vec![1.0]];
         let g = Game::new(
             vec![vec![t_a], vec![t_b], vec![t_c]],
-            vec![0.0, 0.0, 0.0],
+            &[0.0, 0.0, 0.0],
-            vec![],
+            &w,
             0.25,
         );
         let p = g.posteriors();
@@ -429,7 +414,8 @@ mod tests {
             Player::new(Gaussian::from_ms(16., 3.0), 25.0 / 6.0, 25.0 / 300.0),
         ];
 
-        let g = Game::new(vec![t_a, t_b, t_c], vec![1.0, 0.0, 0.0], vec![], 0.25);
+        let w = [vec![1.0, 1.0], vec![1.0], vec![1.0, 1.0]];
+        let g = Game::new(vec![t_a, t_b, t_c], &[1.0, 0.0, 0.0], &w, 0.25);
         let p = g.posteriors();
 
         assert_ulps_eq!(p[0][0], Gaussian::from_ms(13.051, 2.864), epsilon = 1e-3);
@@ -455,7 +441,8 @@ mod tests {
             0.0,
         )];
 
-        let g = Game::new(vec![t_a.clone(), t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a.clone(), t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -472,7 +459,8 @@ mod tests {
         let w_a = vec![1.0];
         let w_b = vec![0.7];
 
-        let g = Game::new(vec![t_a.clone(), t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a.clone(), t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -489,7 +477,8 @@ mod tests {
         let w_a = vec![1.6];
         let w_b = vec![0.7];
 
-        let g = Game::new(vec![t_a, t_b], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a, t_b], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -509,7 +498,8 @@ mod tests {
         let t_a = vec![Player::new(Gaussian::from_ms(2.0, 6.0), 1.0, 0.0)];
         let t_b = vec![Player::new(Gaussian::from_ms(2.0, 6.0), 1.0, 0.0)];
 
-        let g = Game::new(vec![t_a, t_b], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a, t_b], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -529,7 +519,8 @@ mod tests {
         let t_a = vec![Player::new(Gaussian::from_ms(2.0, 6.0), 1.0, 0.0)];
         let t_b = vec![Player::new(Gaussian::from_ms(2.0, 6.0), 1.0, 0.0)];
 
-        let g = Game::new(vec![t_a, t_b], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a, t_b], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(p[0][0], p[1][0], epsilon = 1e-6);
@@ -549,7 +540,8 @@ mod tests {
         ];
         let w_b = vec![0.9, 0.6];
 
-        let g = Game::new(vec![t_a.clone(), t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a.clone(), t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -576,7 +568,8 @@ mod tests {
         let w_a = vec![1.3, 1.5];
         let w_b = vec![0.7, 0.4];
 
-        let g = Game::new(vec![t_a.clone(), t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a.clone(), t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -603,7 +596,8 @@ mod tests {
         let w_a = vec![1.6, 0.2];
         let w_b = vec![0.7, 2.4];
 
-        let g = Game::new(vec![t_a.clone(), t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a.clone(), t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(
@@ -627,6 +621,7 @@ mod tests {
             epsilon = 1e-6
         );
 
+        let w = [vec![1.0, 1.0], vec![1.0]];
         let g = Game::new(
             vec![
                 t_a.clone(),
@@ -636,8 +631,8 @@ mod tests {
                     0.0,
                 )],
             ],
-            vec![],
+            &[1.0, 0.0],
-            vec![],
+            &w,
             0.0,
         );
         let post_2vs1 = g.posteriors();
@@ -645,7 +640,8 @@ mod tests {
         let w_a = vec![1.0, 1.0];
         let w_b = vec![1.0, 0.0];
 
-        let g = Game::new(vec![t_a, t_b.clone()], vec![], vec![w_a, w_b], 0.0);
+        let w = [w_a, w_b];
+        let g = Game::new(vec![t_a, t_b.clone()], &[1.0, 0.0], &w, 0.0);
         let p = g.posteriors();
 
         assert_ulps_eq!(p[0][0], post_2vs1[0][0], epsilon = 1e-6);

src/gaussian2.rs

@@ -1,159 +0,0 @@
-use std::ops;
-
-#[derive(Clone, Copy, PartialEq, Debug)]
-pub struct Gaussian {
-    mu: f64,
-    sigma: f64,
-}
-
-impl Gaussian {
-    #[inline(always)]
-    pub const fn from_ms(mu: f64, sigma: f64) -> Self {
-        Self { mu, sigma }
-    }
-
-    #[inline(always)]
-    pub fn from_pt(pi: f64, tau: f64) -> Self {
-        Self::from_ms(tau / pi, (1.0 / pi).sqrt())
-    }
-
-    #[inline(always)]
-    pub fn mu(&self) -> f64 {
-        self.mu
-    }
-
-    #[inline(always)]
-    pub fn sigma(&self) -> f64 {
-        self.sigma
-    }
-
-    #[inline(always)]
-    pub fn pi(&self) -> f64 {
-        if self.sigma > 0.0 {
-            self.sigma.powi(-2)
-        } else {
-            f64::INFINITY
-        }
-    }
-
-    #[inline(always)]
-    pub fn tau(&self) -> f64 {
-        if self.sigma > 0.0 {
-            self.mu * self.pi()
-        } else {
-            f64::INFINITY
-        }
-    }
-}
-
-impl ops::Add<Gaussian> for Gaussian {
-    type Output = Gaussian;
-
-    fn add(self, rhs: Gaussian) -> Self::Output {
-        Self {
-            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 {
-        Self {
-            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 {
-        /*
-        if self.sigma == 0.0 || rhs.sigma == 0.0 {
-            let mu = self.mu / (self.sigma.powi(2) / rhs.sigma.powi(2) + 1.0)
-                + rhs.mu / (rhs.sigma.powi(2) / self.sigma.powi(2) + 1.0);
-
-            let sigma = (1.0 / ((1.0 / self.sigma.powi(2)) + (1.0 / rhs.sigma.powi(2)))).sqrt();
-
-            Self::from_ms(mu, sigma)
-        } else {
-            Self::from_pt(self.pi() + rhs.pi(), self.tau() + rhs.tau())
-        }
-        */
-
-        Self::from_pt(self.pi() + rhs.pi(), self.tau() + rhs.tau())
-    }
-}
-
-impl ops::Div<Gaussian> for Gaussian {
-    type Output = Gaussian;
-
-    fn div(self, rhs: Gaussian) -> Self::Output {
-        /*
-        let (mu, sigma) = if self.sigma == 0.0 || rhs.sigma == 0.0 {
-            let mu = self.mu / (1.0 - self.sigma.powi(2) / rhs.sigma.powi(2))
-                + rhs.mu / (rhs.sigma.powi(2) / self.sigma.powi(2) - 1.0);
-
-            let sigma = (1.0 / ((1.0 / self.sigma.powi(2)) - (1.0 / rhs.sigma.powi(2)))).sqrt();
-
-            Self::from_ms(mu, sigma)
-        } else {
-            Self::from_pt(self.pi() - rhs.pi(), self.tau() - rhs.tau())
-        }
-        */
-
-        Self::from_pt(self.pi() - rhs.pi(), self.tau() - rhs.tau())
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn test_add() {
-        let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
-        let m = Gaussian::from_ms(0.0, 1.0);
-
-        assert_eq!(n + m, Gaussian::from_ms(25.0, 8.393118874676116));
-    }
-
-    #[test]
-    fn test_sub() {
-        let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
-        let m = Gaussian::from_ms(1.0, 1.0);
-
-        assert_eq!(n - m, Gaussian::from_ms(24.0, 8.393118874676116));
-    }
-
-    #[test]
-    fn test_mul() {
-        let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
-        let m = Gaussian::from_ms(0.0, 1.0);
-
-        assert_eq!(
-            n * m,
-            Gaussian::from_ms(0.35488958990536273, 0.992876838486922)
-        );
-    }
-
-    #[test]
-    fn test_div() {
-        let n = Gaussian::from_ms(25.0, 25.0 / 3.0);
-        let m = Gaussian::from_ms(0.0, 1.0);
-
-        assert_eq!(
-            m / n,
-            Gaussian::from_ms(-0.3652597402597402, 1.0072787050317253)
-        );
-
-        assert_eq!(
-            n / m,
-            Gaussian::from_ms(-0.3652597402597402, 1.0072787050317253)
-        );
-    }
-}

src/history.rs

@@ -1,4 +1,4 @@
-use std::collections::{HashMap, HashSet};
+use std::collections::HashMap;
 
 use crate::{
     agent::{self, Agent},
@@ -247,7 +247,7 @@ impl History {
         results: Vec<Vec<f64>>,
         times: Vec<i64>,
         weights: Vec<Vec<Vec<f64>>>,
-        priors: HashMap<Index, Player>,
+        mut priors: HashMap<Index, Player>,
     ) {
         assert!(times.is_empty() || self.time, "length(times)>0 but !h.time");
         assert!(
@@ -267,19 +267,22 @@ impl History {
             "(length(weights) > 0) & (length(composition) != length(weights))"
         );
 
-        let this_agent = composition
+        agent::clean(self.agents.values_mut(), true);
-            .iter()
+
-            .flatten()
+        let mut this_agent = Vec::with_capacity(1024);
-            .flatten()
+
-            .cloned()
+        for agent in composition.iter().flatten().flatten() {
-            .collect::<HashSet<_>>();
+            if this_agent.contains(agent) {
+                continue;
+            }
+
+            this_agent.push(*agent);
 
-        for agent in &this_agent {
             if !self.agents.contains_key(agent) {
                 self.agents.insert(
                     *agent,
                     Agent {
-                        player: priors.get(agent).cloned().unwrap_or_else(|| {
+                        player: priors.remove(agent).unwrap_or_else(|| {
                             Player::new(
                                 Gaussian::from_ms(self.mu, self.sigma),
                                 self.beta,
@@ -292,8 +295,6 @@ impl History {
             }
         }
 
-        agent::clean(self.agents.values_mut(), true);
-
         let n = composition.len();
         let o = if self.time {
             sort_time(&times, false)
@@ -315,26 +316,23 @@ impl History {
             while (!self.time && (self.size > k))
                 || (self.time && self.batches.len() > k && self.batches[k].time < t)
             {
-                let b = &mut self.batches[k];
+                let batch = &mut self.batches[k];
 
                 if k > 0 {
-                    b.new_forward_info(&self.agents);
+                    batch.new_forward_info(&self.agents);
                 }
 
-                let intersect = this_agent
+                // TODO: Is it faster to iterate over agents in batch instead?
-                    .iter()
+                for agent_idx in &this_agent {
-                    .filter(|&agent| b.skills.contains_key(agent))
+                    if let Some(skill) = batch.skills.get_mut(agent_idx) {
-                    .cloned()
+                        skill.elapsed =
-                    .collect::<Vec<_>>();
+                            batch::compute_elapsed(self.agents[agent_idx].last_time, batch.time);
 
-                for agent in &intersect {
+                        let agent = self.agents.get_mut(agent_idx).unwrap();
-                    b.skills.get_mut(agent).unwrap().elapsed =
-                        batch::compute_elapsed(self.agents[agent].last_time, b.time);
 
-                    let a = self.agents.get_mut(agent).unwrap();
+                        agent.last_time = if self.time { batch.time } else { i64::MAX };
-
+                        agent.message = batch.forward_prior_out(agent_idx);
-                    a.last_time = if self.time { b.time } else { i64::MAX };
+                    }
-                    a.message = b.forward_prior_out(agent);
                 }
 
                 k += 1;
@@ -343,6 +341,7 @@ impl History {
             let composition = (i..j)
                 .map(|e| composition[o[e]].clone())
                 .collect::<Vec<_>>();
+
             let results = if results.is_empty() {
                 Vec::new()
             } else {
@@ -356,28 +355,28 @@ impl History {
             };
 
             if self.time && self.batches.len() > k && self.batches[k].time == t {
-                let b = &mut self.batches[k];
+                let batch = &mut self.batches[k];
+                batch.add_events(composition, results, weights, &self.agents);
 
-                b.add_events(composition, results, weights, &self.agents);
+                for agent_idx in batch.skills.keys() {
-
+                    let agent = self.agents.get_mut(agent_idx).unwrap();
-                for a in b.skills.keys() {
-                    let agent = self.agents.get_mut(a).unwrap();
 
                     agent.last_time = if self.time { t } else { i64::MAX };
-                    agent.message = b.forward_prior_out(a);
+                    agent.message = batch.forward_prior_out(agent_idx);
                 }
             } else {
-                let b = Batch::new(composition, results, weights, t, self.p_draw, &self.agents);
+                let mut batch: Batch = Batch::new(t, self.p_draw);
+                batch.add_events(composition, results, weights, &self.agents);
 
-                self.batches.insert(k, b);
+                self.batches.insert(k, batch);
 
-                let b = &self.batches[k];
+                let batch = &self.batches[k];
 
-                for a in b.skills.keys() {
+                for agent_idx in batch.skills.keys() {
-                    let agent = self.agents.get_mut(a).unwrap();
+                    let agent = self.agents.get_mut(agent_idx).unwrap();
 
                     agent.last_time = if self.time { t } else { i64::MAX };
-                    agent.message = b.forward_prior_out(a);
+                    agent.message = batch.forward_prior_out(agent_idx);
                 }
 
                 k += 1;
@@ -387,24 +386,21 @@ impl History {
         }
 
         while self.time && self.batches.len() > k {
-            let b = &mut self.batches[k];
+            let batch = &mut self.batches[k];
 
-            b.new_forward_info(&self.agents);
+            batch.new_forward_info(&self.agents);
 
-            let intersect = this_agent
+            // TODO: Is it faster to iterate over agents in batch instead?
-                .iter()
+            for agent_idx in &this_agent {
-                .filter(|&agent| b.skills.contains_key(agent))
+                if let Some(skill) = batch.skills.get_mut(agent_idx) {
-                .cloned()
+                    skill.elapsed =
-                .collect::<Vec<_>>();
+                        batch::compute_elapsed(self.agents[agent_idx].last_time, batch.time);
 
-            for agent in &intersect {
+                    let agent = self.agents.get_mut(agent_idx).unwrap();
-                b.skills.get_mut(agent).unwrap().elapsed =
-                    batch::compute_elapsed(self.agents[agent].last_time, b.time);
 
-                let a = self.agents.get_mut(agent).unwrap();
+                    agent.last_time = if self.time { batch.time } else { i64::MAX };
-
+                    agent.message = batch.forward_prior_out(agent_idx);
-                a.last_time = if self.time { b.time } else { i64::MAX };
+                }
-                a.message = b.forward_prior_out(agent);
             }
 
             k += 1;
@@ -470,10 +466,11 @@ mod tests {
 
         let observed = h.batches[1].skills[&a].posterior();
 
+        let w = [vec![1.0], vec![1.0]];
         let p = Game::new(
             h.batches[1].events[0].within_priors(false, false, &h.batches[1].skills, &h.agents),
-            vec![0.0, 1.0],
+            &[0.0, 1.0],
-            vec![],
+            &w,
             P_DRAW,
         )
         .posteriors();

src/lib.rs

@@ -10,14 +10,15 @@ mod approx;
 pub mod batch;
 mod game;
 pub mod gaussian;
-// mod gaussian2;
 mod history;
+mod matrix;
 mod message;
 pub mod player;
 
 pub use game::Game;
 pub use gaussian::Gaussian;
 pub use history::History;
+use matrix::Matrix;
 use message::DiffMessage;
 pub use player::Player;
 
@@ -255,8 +256,72 @@ pub(crate) fn evidence(d: &[DiffMessage], margin: &[f64], tie: &[bool], e: usize
     }
 }
 
+/// Calculates the match quality of the given rating groups. A result is the draw probability in the association
+pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
+    let flatten_ratings = rating_groups
+        .iter()
+        .flat_map(|group| group.iter())
+        .collect::<Vec<_>>();
+
+    let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice();
+
+    let length = flatten_ratings.len();
+
+    let mut mean_matrix = Matrix::new(length, 1);
+
+    for (i, rating) in flatten_ratings.iter().enumerate() {
+        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.sigma.powi(2);
+    }
+
+    let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length);
+
+    let mut t = 0;
+    let mut x = 0;
+
+    for (row, group) in rating_groups.windows(2).enumerate() {
+        let current = group[0];
+        let next = group[1];
+
+        for n in t..t + current.len() {
+            rotated_a_matrix[(row, n)] = flatten_weights[n];
+
+            x += 1;
+        }
+
+        t += current.len();
+
+        for n in x..x + next.len() {
+            rotated_a_matrix[(row, n)] = -flatten_weights[n];
+        }
+
+        x += next.len();
+    }
+
+    let a_matrix = rotated_a_matrix.transpose();
+
+    let ata = beta.powi(2) * &rotated_a_matrix * &a_matrix;
+    let atsa = &rotated_a_matrix * &variance_matrix * &a_matrix;
+
+    let start = mean_matrix.transpose() * &a_matrix;
+    let middle = &ata + &atsa;
+    let end = &rotated_a_matrix * &mean_matrix;
+
+    let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
+    let s_arg = ata.determinant() / middle.determinant();
+
+    e_arg.exp() * s_arg.sqrt()
+}
+
 #[cfg(test)]
 mod tests {
+    use ::approx::assert_ulps_eq;
+
     use super::*;
 
     #[test]
@@ -268,4 +333,14 @@ mod tests {
     fn test_sort_time() {
         assert_eq!(sort_time(&[0, 1, 2, 0], true), vec![2, 1, 0, 3]);
     }
+
+    #[test]
+    fn test_quality() {
+        let a = Gaussian::from_ms(25.0, 3.0);
+        let b = Gaussian::from_ms(25.0, 3.0);
+
+        let q = quality(&[&[a], &[b]], 25.0 / 3.0 / 2.0);
+
+        assert_ulps_eq!(q, 0.8115343414514944, epsilon = 1e-6)
+    }
 }

src/matrix.rs

@@ -0,0 +1,213 @@
+use std::ops;
+
+fn det(m: &[f64], x: usize) -> f64 {
+    if x == 1 {
+        m[0]
+    } else if x == 2 {
+        m[0] * m[3] - m[1] * m[2]
+    } else {
+        let mut d = 0.0;
+
+        for n in 0..x {
+            let ms = m
+                .iter()
+                .enumerate()
+                .skip(x)
+                .filter(|(i, _)| (i % x) != n)
+                .map(|(_, v)| *v)
+                .collect::<Vec<_>>();
+
+            d += (-1.0f64).powi(n as i32) * m[n] * det(&ms, x - 1);
+        }
+
+        d
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct Matrix {
+    data: Box<[f64]>,
+    height: usize,
+    width: usize,
+}
+
+impl Matrix {
+    pub fn new(height: usize, width: usize) -> Matrix {
+        Matrix {
+            data: vec![0.0; height * width].into_boxed_slice(),
+            height,
+            width,
+        }
+    }
+
+    pub fn transpose(&self) -> Matrix {
+        let mut matrix = Matrix::new(self.width, self.height);
+
+        for c in 0..self.width {
+            for r in 0..self.height {
+                matrix[(c, r)] = self[(r, c)];
+            }
+        }
+
+        matrix
+    }
+
+    pub fn minor(&self, row_n: usize, col_n: usize) -> Matrix {
+        let mut matrix = Matrix::new(self.height - 1, self.width - 1);
+
+        let mut nr = 0;
+
+        for r in 0..self.height {
+            if r == row_n {
+                continue;
+            }
+
+            let mut nc = 0;
+
+            for c in 0..self.width {
+                if c == col_n {
+                    continue;
+                }
+
+                matrix[(nr, nc)] = self[(r, c)];
+
+                nc += 1;
+            }
+
+            nr += 1;
+        }
+
+        matrix
+    }
+
+    pub fn determinant(&self) -> f64 {
+        debug_assert!(self.width == self.height);
+
+        det(&self.data, self.width)
+    }
+
+    pub fn adjugate(&self) -> Matrix {
+        debug_assert!(self.width == self.height);
+
+        let mut matrix = Matrix::new(self.height, self.width);
+
+        if matrix.height == 2 {
+            matrix[(0, 0)] = self[(1, 1)];
+            matrix[(0, 1)] = -self[(0, 1)];
+            matrix[(1, 0)] = -self[(1, 0)];
+            matrix[(1, 1)] = self[(0, 0)];
+        } else {
+            for r in 0..matrix.height {
+                for c in 0..matrix.width {
+                    let sign = if (r + c) % 2 == 0 { 1.0 } else { -1.0 };
+
+                    matrix[(r, c)] = self.minor(r, c).determinant() * sign;
+                }
+            }
+        }
+
+        matrix
+    }
+
+    pub fn inverse(&self) -> Matrix {
+        let mut matrix = Matrix::new(self.width, self.height);
+
+        if self.height == self.width && self.height == 1 {
+            matrix[(0, 0)] = 1.0 / self[(0, 0)];
+        } else {
+            panic!("eh, okey")
+        }
+
+        matrix
+    }
+}
+
+impl ops::Index<(usize, usize)> for Matrix {
+    type Output = f64;
+
+    fn index(&self, pos: (usize, usize)) -> &Self::Output {
+        &self.data[(self.width * pos.0) + pos.1]
+    }
+}
+
+impl ops::IndexMut<(usize, usize)> for Matrix {
+    fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output {
+        &mut self.data[(self.width * pos.0) + pos.1]
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for f64 {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(rhs.height, rhs.width);
+
+        for r in 0..rhs.height {
+            for c in 0..rhs.width {
+                matrix[(r, c)] = self * rhs[(r, c)];
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for Matrix {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, rhs.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                let mut value = 0.0;
+
+                for x in 0..self.width {
+                    value += self[(r, x)] * rhs[(x, c)];
+                }
+
+                matrix[(r, c)] = value;
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for &'a Matrix {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, rhs.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                let mut value = 0.0;
+
+                for x in 0..self.width {
+                    value += self[(r, x)] * rhs[(x, c)];
+                }
+
+                matrix[(r, c)] = value;
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
+    type Output = Matrix;
+
+    fn add(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, self.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                matrix[(r, c)] = self[(r, c)] + rhs[(r, c)];
+            }
+        }
+
+        matrix
+    }
+}

src/message.rs

@@ -15,14 +15,17 @@ impl TeamMessage {
     }
     */
 
+    #[inline]
     pub(crate) fn posterior_win(&self) -> Gaussian {
         self.prior * self.likelihood_lose * self.likelihood_draw
     }
 
+    #[inline]
     pub(crate) fn posterior_lose(&self) -> Gaussian {
         self.prior * self.likelihood_win * self.likelihood_draw
     }
 
+    #[inline]
     pub(crate) fn likelihood(&self) -> Gaussian {
         self.likelihood_win * self.likelihood_lose * self.likelihood_draw
     }