feat: added a Drift trait and a "default" ConstantDrift implementation

.gitignore

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

CLAUDE.md

@@ -0,0 +1,45 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Commands
+
+```bash
+cargo build                          # Build the library
+cargo test --lib                     # Run all library tests
+cargo test --lib <test_name>         # Run a single test by name
+cargo test --lib -- --nocapture      # Run tests with stdout output
+cargo clippy                         # Lint
+cargo bench                          # Run benchmarks (criterion)
+```
+
+The `approx` feature enables `approx::AbsDiffEq` for `Gaussian`:
+```bash
+cargo test --features approx
+```
+
+## Architecture
+
+This is a Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py) — a Bayesian skill rating system that tracks skill evolution over time using Gaussian message passing.
+
+### Data flow
+
+```
+History  →  Batch[]  →  Game[]  →  teams/players
+```
+
+- **`History`** (`history.rs`) — top-level container. Organizes games by time into `Batch`es, runs forward/backward message passing across batches, and exposes `learning_curves()` and `log_evidence()`.
+- **`Batch`** (`batch.rs`) — all games at a single time step. Runs `iteration()` to update skill estimates via `Game::posteriors()`, collecting `Skill` distributions per player.
+- **`Game`** (`game.rs`) — a single match. Given teams (slices of `Gaussian`), computes posterior skill distributions using Gaussian factor graphs and `message.rs` helpers.
+- **`Agent`** (`agent.rs`) — wraps a `Player` with temporal state (`last_time`, `message`). `receive()` applies time-decay (`gamma`) when the player reappears after a gap.
+- **`Player`** (`player.rs`) — static configuration: prior `Gaussian`, `beta` (performance noise), `gamma` (skill drift per time unit).
+- **`Gaussian`** (`gaussian.rs`) — core probability type. Stored as natural parameters (`pi = 1/sigma²`, `tau = mu/sigma²`). Arithmetic ops implement message multiplication/division in the factor graph.
+- **`message.rs`** — `TeamMessage` and `DiffMessage`: intermediate factor graph messages used inside `Game`.
+- **`lib.rs`** — exports the public API (`Game`, `Gaussian`, `History`, `Player`) and standalone functions (`quality()`, `pdf()`, `cdf()`, `erfc()`). Also defines global defaults: `MU=0.0`, `SIGMA=6.0`, `BETA=1.0`, `GAMMA=0.03`, `P_DRAW=0.0`, `EPSILON=1e-6`, `ITERATIONS=30`.
+
+### Key design points
+
+- `History` uses `IndexMap<K>` (defined in `lib.rs`) to map arbitrary player keys to `Agent` state.
+- Convergence is measured by the maximum `delta()` across all skill distributions; iteration stops when below `EPSILON` or after `ITERATIONS` rounds.
+- The `approx` feature gates `AbsDiffEq` on `Gaussian` for use in tests — the feature is optional and only needed for approximate equality assertions.
+- `time` in `History`/`Batch` is currently an `f64`; the README notes it needs to become an enum to support richer temporal states.

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "trueskill-tt"
 version = "0.1.0"
-edition = "2021"
+edition = "2024"
 
 [lib]
 bench = false

NOTEPAD.md

@@ -1,15 +1,16 @@
 # History
 
-```shell
-teams: [[player]]
-weights: [[f64]]
-results: [f64]
+```rust
+let mut history = History::new();
 
-player: (gaussian, f64, f64)
-
-players: [player]
-weights: [f64]
-
-
-teams: [([(player, weight)], result)]
+let agent_a = history.new_agent();
+let agent_b = history.new_agent_with_prior(Prior::new(Gaussian::default(), BETA, GAMMA));
+```
+
+```rust
+trait Team {
+    fn players(&self) -> impl Iterator<Item = P>;
+    fn weights(&self) -> impl Iterator<Item = f64>;
+    fn score(&self) -> u16;
+}
 ```

README.md

@@ -11,10 +11,71 @@ Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillTh
 - [TrueSkill Through Time: Revisiting the History of Chess](https://www.microsoft.com/en-us/research/wp-content/uploads/2008/01/NIPS2007_0931.pdf)
 - [TrueSkill Through Time. The full scientific documentation](https://glandfried.github.io/publication/landfried2021-learning/)
 
+## Drift
+
+Skill drift models how a player's true skill can change between appearances. Each time a player reappears after a gap, their skill uncertainty is widened by the drift model before the new evidence is incorporated.
+
+Drift is represented by the `Drift` trait:
+
+```rust
+pub trait Drift: Copy + Debug {
+    fn variance_delta(&self, elapsed: i64) -> f64;
+}
+```
+
+`variance_delta` returns the amount to add to `σ²` given the elapsed time since the player last played. Internally, `Gaussian::forget` uses this to compute the new sigma: `σ_new = sqrt(σ² + variance_delta)`.
+
+### ConstantDrift
+
+The built-in `ConstantDrift` implements a linear random walk — skill uncertainty grows proportionally to time:
+
+```
+variance_delta = elapsed * γ²
+```
+
+This is the standard TrueSkill Through Time model. Use it by passing a `ConstantDrift(gamma)` when constructing a `Player`:
+
+```rust
+use trueskill_tt::{Player, Gaussian, drift::ConstantDrift};
+
+// gamma = 0.1 means skill can shift ~0.1 per time unit
+let player = Player::new(Gaussian::from_ms(0.0, 6.0), 1.0, ConstantDrift(0.1));
+```
+
+### Custom drift
+
+Implement `Drift` to express any other model. For example, a drift that saturates after a long absence (uncertainty grows with the square root of elapsed time instead of linearly):
+
+```rust
+use trueskill_tt::drift::Drift;
+
+#[derive(Clone, Copy, Debug)]
+struct SqrtDrift {
+    gamma: f64,
+}
+
+impl Drift for SqrtDrift {
+    fn variance_delta(&self, elapsed: i64) -> f64 {
+        (elapsed as f64).sqrt() * self.gamma * self.gamma
+    }
+}
+
+let player = Player::new(Gaussian::from_ms(0.0, 6.0), 1.0, SqrtDrift { gamma: 0.5 });
+```
+
+To use a custom drift type with `History`, use the `.drift()` builder method instead of `.gamma()`:
+
+```rust
+let h = History::builder()
+    .drift(SqrtDrift { gamma: 0.5 })
+    .build();
+```
+
 ## Todo
 
 - [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

@@ -1,9 +1,9 @@
 use std::collections::HashMap;
 
-use criterion::{criterion_group, criterion_main, Criterion};
+use criterion::{Criterion, criterion_group, criterion_main};
 use trueskill_tt::{
-    agent::Agent, batch::Batch, gaussian::Gaussian, player::Player, IndexMap, BETA, GAMMA, MU,
-    P_DRAW, SIGMA,
+    BETA, GAMMA, IndexMap, MU, P_DRAW, SIGMA, agent::Agent, batch::Batch, drift::ConstantDrift,
+    gaussian::Gaussian, player::Player,
 };
 
 fn criterion_benchmark(criterion: &mut Criterion) {
@@ -19,21 +19,21 @@ fn criterion_benchmark(criterion: &mut Criterion) {
         map.insert(
             a,
             Agent {
-                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, GAMMA),
+                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, ConstantDrift(GAMMA)),
                 ..Default::default()
             },
         );
         map.insert(
             b,
             Agent {
-                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, GAMMA),
+                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, ConstantDrift(GAMMA)),
                 ..Default::default()
             },
         );
         map.insert(
             c,
             Agent {
-                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, GAMMA),
+                player: Player::new(Gaussian::from_ms(MU, SIGMA), BETA, ConstantDrift(GAMMA)),
                 ..Default::default()
             },
         );
@@ -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))

graph.d2

@@ -0,0 +1,64 @@
+vars: {
+  d2-config: {
+    layout-engine: elk
+    # Terminal theme code
+    theme-id: 300
+  }
+}
+
+History: {
+  shape: class
+
+  agents: "HashMap<Index, Agent>"
+  batches: "Vec<Batch>"
+}
+
+Batch: {
+  shape: class
+
+  skills: "HashMap<Index, Skill>"
+  events: "Vec<Event>"
+  time: "i64"
+  p_draw: "f64"
+}
+
+Event: {
+  shape: class
+
+  teams: "Vec<Team>"
+  weights: "Vec<Vec<f64>>"
+  evidence: "f64"
+}
+
+Team: {
+  shape: class
+
+  items: "Vec<Item>"
+  output: "f64"
+}
+
+Item: {
+  shape: class
+
+  agent: "Index"
+  likelihood: "Gaussian"
+}
+
+Skill: {
+  shape: class
+
+  forward: "Gaussian"
+  backward: "Gaussian"
+  likelihood: "Gaussian"
+  elapsed: "i64"
+  online: "Gaussian"
+}
+
+History -> Batch
+
+Batch -> Skill
+Batch -> Event
+
+Event -> Team
+
+Team -> Item

src/agent.rs

@@ -1,23 +1,29 @@
-use crate::{gaussian::Gaussian, player::Player, N_INF};
+use crate::{
+    N_INF,
+    drift::{ConstantDrift, Drift},
+    gaussian::Gaussian,
+    player::Player,
+};
 
 #[derive(Debug)]
-pub struct Agent {
-    pub player: Player,
+pub struct Agent<D: Drift = ConstantDrift> {
+    pub player: Player<D>,
     pub message: Gaussian,
     pub last_time: i64,
 }
 
-impl Agent {
+impl<D: Drift> Agent<D> {
     pub(crate) fn receive(&self, elapsed: i64) -> Gaussian {
         if self.message != N_INF {
-            self.message.forget(self.player.gamma, elapsed)
+            self.message
+                .forget(self.player.drift.variance_delta(elapsed))
         } else {
             self.player.prior
         }
     }
 }
 
-impl Default for Agent {
+impl Default for Agent<ConstantDrift> {
     fn default() -> Self {
         Self {
             player: Player::default(),
@@ -27,7 +33,10 @@ impl Default for Agent {
     }
 }
 
-pub(crate) fn clean<'a, A: Iterator<Item = &'a mut Agent>>(agents: A, last_time: bool) {
+pub(crate) fn clean<'a, D: Drift + 'a, A: Iterator<Item = &'a mut Agent<D>>>(
+    agents: A,
+    last_time: bool,
+) {
     for a in agents {
         a.message = N_INF;
 

src/batch.rs

@@ -1,7 +1,8 @@
-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,
+    Index, N_INF, agent::Agent, drift::Drift, game::Game, gaussian::Gaussian, player::Player,
+    tuple_gt, tuple_max,
 };
 
 #[derive(Debug)]
@@ -38,22 +39,22 @@ struct Item {
 }
 
 impl Item {
-    fn within_prior(
+    fn within_prior<D: Drift>(
         &self,
         online: bool,
         forward: bool,
         skills: &HashMap<Index, Skill>,
-        agents: &HashMap<Index, Agent>,
-    ) -> Player {
+        agents: &HashMap<Index, Agent<D>>,
+    ) -> Player<D> {
         let r = &agents[&self.agent].player;
         let skill = &skills[&self.agent];
 
         if online {
-            Player::new(skill.online, r.beta, r.gamma)
+            Player::new(skill.online, r.beta, r.drift)
         } else if forward {
-            Player::new(skill.forward, r.beta, r.gamma)
+            Player::new(skill.forward, r.beta, r.drift)
         } else {
-            Player::new(skill.posterior() / self.likelihood, r.beta, r.gamma)
+            Player::new(skill.posterior() / self.likelihood, r.beta, r.drift)
         }
     }
 }
@@ -79,13 +80,13 @@ impl Event {
             .collect::<Vec<_>>()
     }
 
-    pub(crate) fn within_priors(
+    pub(crate) fn within_priors<D: Drift>(
         &self,
         online: bool,
         forward: bool,
         skills: &HashMap<Index, Skill>,
-        agents: &HashMap<Index, Agent>,
-    ) -> Vec<Vec<Player>> {
+        agents: &HashMap<Index, Agent<D>>,
+    ) -> Vec<Vec<Player<D>>> {
         self.teams
             .iter()
             .map(|team| {
@@ -107,119 +108,43 @@ pub struct Batch {
 }
 
 impl Batch {
-    pub fn new(
-        composition: Vec<Vec<Vec<Index>>>,
-        results: Vec<Vec<f64>>,
-        weights: Vec<Vec<Vec<f64>>>,
-        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 {
+    pub fn new(time: i64, p_draw: f64) -> Self {
+        Self {
+            events: Vec::new(),
+            skills: HashMap::new(),
             time,
-            events,
-            skills,
             p_draw,
-        };
-
-        this.iteration(0, agents);
-
-        this
+        }
     }
 
-    pub(crate) fn add_events(
+    pub fn add_events<D: Drift>(
         &mut self,
         composition: Vec<Vec<Vec<Index>>>,
         results: Vec<Vec<f64>>,
         weights: Vec<Vec<Vec<f64>>>,
-        agents: &HashMap<Index, Agent>,
+        agents: &HashMap<Index, Agent<D>>,
     ) {
-        let this_agent = composition
-            .iter()
-            .flatten()
-            .flatten()
-            .cloned()
-            .collect::<HashSet<_>>();
+        let mut unique = Vec::with_capacity(10);
+
+        let this_agent = composition.iter().flatten().flatten().filter(|idx| {
+            if !unique.contains(idx) {
+                unique.push(*idx);
+
+                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 +178,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() {
-                    Vec::new()
-                } else {
-                    weights[e].clone()
-                },
+                weights,
             }
         });
 
@@ -278,12 +208,12 @@ impl Batch {
             .collect::<HashMap<_, _>>()
     }
 
-    pub fn iteration(&mut self, from: usize, agents: &HashMap<Index, Agent>) {
+    pub fn iteration<D: Drift>(&mut self, from: usize, agents: &HashMap<Index, Agent<D>>) {
         for event in self.events.iter_mut().skip(from) {
             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() {
@@ -300,7 +230,7 @@ impl Batch {
     }
 
     #[allow(dead_code)]
-    pub(crate) fn convergence(&mut self, agents: &HashMap<Index, Agent>) -> usize {
+    pub(crate) fn convergence<D: Drift>(&mut self, agents: &HashMap<Index, Agent<D>>) -> usize {
         let epsilon = 1e-6;
         let iterations = 20;
 
@@ -330,18 +260,18 @@ impl Batch {
         skill.forward * skill.likelihood
     }
 
-    pub(crate) fn backward_prior_out(
+    pub(crate) fn backward_prior_out<D: Drift>(
         &self,
         agent: &Index,
-        agents: &HashMap<Index, Agent>,
+        agents: &HashMap<Index, Agent<D>>,
     ) -> Gaussian {
         let skill = &self.skills[agent];
         let n = skill.likelihood * skill.backward;
 
-        n.forget(agents[agent].player.gamma, skill.elapsed)
+        n.forget(agents[agent].player.drift.variance_delta(skill.elapsed))
     }
 
-    pub(crate) fn new_backward_info(&mut self, agents: &HashMap<Index, Agent>) {
+    pub(crate) fn new_backward_info<D: Drift>(&mut self, agents: &HashMap<Index, Agent<D>>) {
         for (agent, skill) in self.skills.iter_mut() {
             skill.backward = agents[agent].message;
         }
@@ -349,7 +279,7 @@ impl Batch {
         self.iteration(0, agents);
     }
 
-    pub(crate) fn new_forward_info(&mut self, agents: &HashMap<Index, Agent>) {
+    pub(crate) fn new_forward_info<D: Drift>(&mut self, agents: &HashMap<Index, Agent<D>>) {
         for (agent, skill) in self.skills.iter_mut() {
             skill.forward = agents[agent].receive(skill.elapsed);
         }
@@ -357,12 +287,12 @@ impl Batch {
         self.iteration(0, agents);
     }
 
-    pub(crate) fn log_evidence(
+    pub(crate) fn log_evidence<D: Drift>(
         &self,
         online: bool,
         targets: &[Index],
         forward: bool,
-        agents: &HashMap<Index, Agent>,
+        agents: &HashMap<Index, Agent<D>>,
     ) -> f64 {
         if targets.is_empty() {
             if online || forward {
@@ -372,8 +302,8 @@ impl Batch {
                     .map(|(_, event)| {
                         Game::new(
                             event.within_priors(online, forward, &self.skills, agents),
-                            event.outputs(),
-                            event.weights.clone(),
+                            &event.outputs(),
+                            &event.weights,
                             self.p_draw,
                         )
                         .evidence
@@ -397,8 +327,8 @@ impl Batch {
                 .map(|(_, event)| {
                     Game::new(
                         event.within_priors(online, forward, &self.skills, agents),
-                        event.outputs(),
-                        event.weights.clone(),
+                        &event.outputs(),
+                        &event.weights,
                         self.p_draw,
                     )
                     .evidence
@@ -461,7 +391,7 @@ pub(crate) fn compute_elapsed(last_time: i64, actual_time: i64) -> i64 {
 mod tests {
     use approx::assert_ulps_eq;
 
-    use crate::{agent::Agent, player::Player, IndexMap};
+    use crate::{IndexMap, agent::Agent, drift::ConstantDrift, player::Player};
 
     use super::*;
 
@@ -485,14 +415,16 @@ mod tests {
                     player: Player::new(
                         Gaussian::from_ms(25.0, 25.0 / 3.0),
                         25.0 / 6.0,
-                        25.0 / 300.0,
+                        ConstantDrift(25.0 / 300.0),
                     ),
                     ..Default::default()
                 },
             );
         }
 
-        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 +432,6 @@ mod tests {
             ],
             vec![vec![1.0, 0.0], vec![0.0, 1.0], vec![1.0, 0.0]],
             vec![],
-            0,
-            0.0,
             &agents,
         );
 
@@ -561,14 +491,16 @@ mod tests {
                     player: Player::new(
                         Gaussian::from_ms(25.0, 25.0 / 3.0),
                         25.0 / 6.0,
-                        25.0 / 300.0,
+                        ConstantDrift(25.0 / 300.0),
                     ),
                     ..Default::default()
                 },
             );
         }
 
-        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 +508,6 @@ mod tests {
             ],
             vec![vec![1.0, 0.0], vec![0.0, 1.0], vec![1.0, 0.0]],
             vec![],
-            0,
-            0.0,
             &agents,
         );
 
@@ -640,14 +570,16 @@ mod tests {
                     player: Player::new(
                         Gaussian::from_ms(25.0, 25.0 / 3.0),
                         25.0 / 6.0,
-                        25.0 / 300.0,
+                        ConstantDrift(25.0 / 300.0),
                     ),
                     ..Default::default()
                 },
             );
         }
 
-        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 +587,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/drift.rs

@@ -0,0 +1,14 @@
+use std::fmt::Debug;
+
+pub trait Drift: Copy + Debug {
+    fn variance_delta(&self, elapsed: i64) -> f64;
+}
+
+#[derive(Clone, Copy, Debug)]
+pub struct ConstantDrift(pub f64);
+
+impl Drift for ConstantDrift {
+    fn variance_delta(&self, elapsed: i64) -> f64 {
+        elapsed as f64 * self.0 * self.0
+    }
+}

src/game.rs

@@ -1,72 +1,57 @@
 use crate::{
-    approx, compute_margin, evidence,
+    N_INF, N00, approx, compute_margin,
+    drift::Drift,
+    evidence,
     gaussian::Gaussian,
     message::{DiffMessage, TeamMessage},
     player::Player,
-    sort_perm, tuple_gt, tuple_max, N00, N_INF,
+    sort_perm, tuple_gt, tuple_max,
 };
 
 #[derive(Debug)]
-pub struct Game {
-    teams: Vec<Vec<Player>>,
-    result: Vec<f64>,
-    weights: Vec<Vec<f64>>,
+pub struct Game<'a, D: Drift> {
+    teams: Vec<Vec<Player<D>>>,
+    result: &'a [f64],
+    weights: &'a [Vec<f64>],
     p_draw: f64,
     pub(crate) likelihoods: Vec<Vec<Gaussian>>,
     pub(crate) evidence: f64,
 }
 
-impl Game {
+impl<'a, D: Drift> Game<'a, D> {
     pub fn new(
-        teams: Vec<Vec<Player>>,
-        mut result: Vec<f64>,
-        mut weights: Vec<Vec<f64>>,
+        teams: Vec<Vec<Player<D>>>,
+        result: &'a [f64],
+        weights: &'a [Vec<f64>],
         p_draw: f64,
     ) -> Self {
-        assert!(
-            (result.is_empty() || result.len() == teams.len()),
-            "result must be empty or the same length as teams"
+        debug_assert!(
+            (result.len() == teams.len()),
+            "result must have the same length as teams"
         );
 
-        assert!(
-            (weights.is_empty() || weights.len() == teams.len()),
-            "weights must be empty or the same length as teams"
+        debug_assert!(
+            weights
+                .iter()
+                .zip(teams.iter())
+                .all(|(w, t)| w.len() == t.len()),
+            "weights must have the same dimensions 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"
-        );
-
-        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 +67,7 @@ impl Game {
     }
 
     fn likelihoods(&mut self) {
-        let m_t_ft = self.likelihood_teams();
-
-        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 o = sort_perm(self.result, true);
 
         let mut team = o
             .iter()
@@ -143,14 +104,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();
+                    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())
-                    }
+                    compute_margin(self.p_draw, (a + b).sqrt())
                 })
                 .collect::<Vec<_>>()
         };
@@ -205,7 +162,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.powi(2))
+                    })
+                    .collect::<Vec<_>>()
+            })
+            .collect::<Vec<_>>();
     }
 
     pub fn posteriors(&self) -> Vec<Vec<Gaussian>> {
@@ -226,7 +203,7 @@ impl Game {
 mod tests {
     use ::approx::assert_ulps_eq;
 
-    use crate::{Gaussian, Player, GAMMA, N_INF};
+    use crate::{ConstantDrift, GAMMA, Gaussian, N_INF, Player};
 
     use super::*;
 
@@ -235,15 +212,16 @@ mod tests {
         let t_a = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(25.0 / 300.0),
         );
         let t_b = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(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];
@@ -252,10 +230,19 @@ mod tests {
         assert_ulps_eq!(a, Gaussian::from_ms(20.794779, 7.194481), epsilon = 1e-6);
         assert_ulps_eq!(b, Gaussian::from_ms(29.205220, 7.194481), epsilon = 1e-6);
 
-        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 t_a = Player::new(
+            Gaussian::from_ms(29.0, 1.0),
+            25.0 / 6.0,
+            ConstantDrift(GAMMA),
+        );
+        let t_b = Player::new(
+            Gaussian::from_ms(25.0, 25.0 / 3.0),
+            25.0 / 6.0,
+            ConstantDrift(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];
@@ -264,10 +251,11 @@ mod tests {
         assert_ulps_eq!(a, Gaussian::from_ms(28.896475, 0.996604), epsilon = 1e-6);
         assert_ulps_eq!(b, Gaussian::from_ms(32.189211, 6.062063), epsilon = 1e-6);
 
-        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 t_a = Player::new(Gaussian::from_ms(1.139, 0.531), 1.0, ConstantDrift(0.2125));
+        let t_b = Player::new(Gaussian::from_ms(15.568, 0.51), 1.0, ConstantDrift(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);
@@ -279,21 +267,22 @@ mod tests {
             vec![Player::new(
                 Gaussian::from_ms(25.0, 25.0 / 3.0),
                 25.0 / 6.0,
-                25.0 / 300.0,
+                ConstantDrift(25.0 / 300.0),
             )],
             vec![Player::new(
                 Gaussian::from_ms(25.0, 25.0 / 3.0),
                 25.0 / 6.0,
-                25.0 / 300.0,
+                ConstantDrift(25.0 / 300.0),
             )],
             vec![Player::new(
                 Gaussian::from_ms(25.0, 25.0 / 3.0),
                 25.0 / 6.0,
-                25.0 / 300.0,
+                ConstantDrift(25.0 / 300.0),
             )],
         ];
 
-        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 +291,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 +301,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];
@@ -328,15 +319,16 @@ mod tests {
         let t_a = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(25.0 / 300.0),
         );
         let t_b = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(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];
@@ -345,10 +337,19 @@ mod tests {
         assert_ulps_eq!(a, Gaussian::from_ms(24.999999, 6.469480), epsilon = 1e-6);
         assert_ulps_eq!(b, Gaussian::from_ms(24.999999, 6.469480), epsilon = 1e-6);
 
-        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 t_a = Player::new(
+            Gaussian::from_ms(25.0, 3.0),
+            25.0 / 6.0,
+            ConstantDrift(25.0 / 300.0),
+        );
+        let t_b = Player::new(
+            Gaussian::from_ms(29.0, 2.0),
+            25.0 / 6.0,
+            ConstantDrift(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];
@@ -363,23 +364,24 @@ mod tests {
         let t_a = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(25.0 / 300.0),
         );
         let t_b = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(25.0 / 300.0),
         );
         let t_c = Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(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],
-            vec![],
+            &[0.0, 0.0, 0.0],
+            &w,
             0.25,
         );
         let p = g.posteriors();
@@ -392,14 +394,27 @@ mod tests {
         assert_ulps_eq!(b, Gaussian::from_ms(25.000000, 5.707423), epsilon = 1e-6);
         assert_ulps_eq!(c, Gaussian::from_ms(24.999999, 5.729068), epsilon = 1e-6);
 
-        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(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 t_a = Player::new(
+            Gaussian::from_ms(25.0, 3.0),
+            25.0 / 6.0,
+            ConstantDrift(25.0 / 300.0),
+        );
+        let t_b = Player::new(
+            Gaussian::from_ms(25.0, 3.0),
+            25.0 / 6.0,
+            ConstantDrift(25.0 / 300.0),
+        );
+        let t_c = Player::new(
+            Gaussian::from_ms(29.0, 2.0),
+            25.0 / 6.0,
+            ConstantDrift(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],
-            vec![],
+            &[0.0, 0.0, 0.0],
+            &w,
             0.25,
         );
         let p = g.posteriors();
@@ -416,20 +431,37 @@ mod tests {
     #[test]
     fn test_2vs1vs2_mixed() {
         let t_a = vec![
-            Player::new(Gaussian::from_ms(12.0, 3.0), 25.0 / 6.0, 25.0 / 300.0),
-            Player::new(Gaussian::from_ms(18.0, 3.0), 25.0 / 6.0, 25.0 / 300.0),
+            Player::new(
+                Gaussian::from_ms(12.0, 3.0),
+                25.0 / 6.0,
+                ConstantDrift(25.0 / 300.0),
+            ),
+            Player::new(
+                Gaussian::from_ms(18.0, 3.0),
+                25.0 / 6.0,
+                ConstantDrift(25.0 / 300.0),
+            ),
         ];
         let t_b = vec![Player::new(
             Gaussian::from_ms(30.0, 3.0),
             25.0 / 6.0,
-            25.0 / 300.0,
+            ConstantDrift(25.0 / 300.0),
         )];
         let t_c = vec![
-            Player::new(Gaussian::from_ms(14.0, 3.0), 25.0 / 6.0, 25.0 / 300.0),
-            Player::new(Gaussian::from_ms(16., 3.0), 25.0 / 6.0, 25.0 / 300.0),
+            Player::new(
+                Gaussian::from_ms(14.0, 3.0),
+                25.0 / 6.0,
+                ConstantDrift(25.0 / 300.0),
+            ),
+            Player::new(
+                Gaussian::from_ms(16., 3.0),
+                25.0 / 6.0,
+                ConstantDrift(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);
@@ -447,15 +479,16 @@ mod tests {
         let t_a = vec![Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            0.0,
+            ConstantDrift(0.0),
         )];
         let t_b = vec![Player::new(
             Gaussian::from_ms(25.0, 25.0 / 3.0),
             25.0 / 6.0,
-            0.0,
+            ConstantDrift(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 +505,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 +523,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!(
@@ -506,10 +541,19 @@ mod tests {
         let w_a = vec![1.0];
         let w_b = vec![0.0];
 
-        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 t_a = vec![Player::new(
+            Gaussian::from_ms(2.0, 6.0),
+            1.0,
+            ConstantDrift(0.0),
+        )];
+        let t_b = vec![Player::new(
+            Gaussian::from_ms(2.0, 6.0),
+            1.0,
+            ConstantDrift(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!(
@@ -526,10 +570,19 @@ mod tests {
         let w_a = vec![1.0];
         let w_b = vec![-1.0];
 
-        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 t_a = vec![Player::new(
+            Gaussian::from_ms(2.0, 6.0),
+            1.0,
+            ConstantDrift(0.0),
+        )];
+        let t_b = vec![Player::new(
+            Gaussian::from_ms(2.0, 6.0),
+            1.0,
+            ConstantDrift(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);
@@ -538,18 +591,35 @@ mod tests {
     #[test]
     fn test_2vs2_weighted() {
         let t_a = vec![
-            Player::new(Gaussian::from_ms(25.0, 25.0 / 3.0), 25.0 / 6.0, 0.0),
-            Player::new(Gaussian::from_ms(25.0, 25.0 / 3.0), 25.0 / 6.0, 0.0),
+            Player::new(
+                Gaussian::from_ms(25.0, 25.0 / 3.0),
+                25.0 / 6.0,
+                ConstantDrift(0.0),
+            ),
+            Player::new(
+                Gaussian::from_ms(25.0, 25.0 / 3.0),
+                25.0 / 6.0,
+                ConstantDrift(0.0),
+            ),
         ];
         let w_a = vec![0.4, 0.8];
 
         let t_b = vec![
-            Player::new(Gaussian::from_ms(25.0, 25.0 / 3.0), 25.0 / 6.0, 0.0),
-            Player::new(Gaussian::from_ms(25.0, 25.0 / 3.0), 25.0 / 6.0, 0.0),
+            Player::new(
+                Gaussian::from_ms(25.0, 25.0 / 3.0),
+                25.0 / 6.0,
+                ConstantDrift(0.0),
+            ),
+            Player::new(
+                Gaussian::from_ms(25.0, 25.0 / 3.0),
+                25.0 / 6.0,
+                ConstantDrift(0.0),
+            ),
         ];
         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 +646,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 +674,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,17 +699,18 @@ mod tests {
             epsilon = 1e-6
         );
 
+        let w = [vec![1.0, 1.0], vec![1.0]];
         let g = Game::new(
             vec![
                 t_a.clone(),
                 vec![Player::new(
                     Gaussian::from_ms(25.0, 25.0 / 3.0),
                     25.0 / 6.0,
-                    0.0,
+                    ConstantDrift(0.0),
                 )],
             ],
-            vec![],
-            vec![],
+            &[1.0, 0.0],
+            &w,
             0.0,
         );
         let post_2vs1 = g.posteriors();
@@ -645,7 +718,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/gaussian.rs

@@ -40,10 +40,10 @@ impl Gaussian {
         }
     }
 
-    pub(crate) fn forget(&self, gamma: f64, t: i64) -> Self {
+    pub(crate) fn forget(&self, variance_delta: f64) -> Self {
         Self {
             mu: self.mu,
-            sigma: (self.sigma.powi(2) + t as f64 * gamma.powi(2)).sqrt(),
+            sigma: (self.sigma.powi(2) + variance_delta).sqrt(),
         }
     }
 }

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,25 +1,27 @@
-use std::collections::{HashMap, HashSet};
+use std::collections::HashMap;
 
 use crate::{
+    BETA, GAMMA, Index, MU, N_INF, P_DRAW, SIGMA,
     agent::{self, Agent},
     batch::{self, Batch},
+    drift::{ConstantDrift, Drift},
     gaussian::Gaussian,
     player::Player,
-    sort_time, tuple_gt, tuple_max, Index, BETA, GAMMA, MU, P_DRAW, SIGMA,
+    sort_time, tuple_gt, tuple_max,
 };
 
 #[derive(Clone)]
-pub struct HistoryBuilder {
+pub struct HistoryBuilder<D: Drift = ConstantDrift> {
     time: bool,
     mu: f64,
     sigma: f64,
     beta: f64,
-    gamma: f64,
+    drift: D,
     p_draw: f64,
     online: bool,
 }
 
-impl HistoryBuilder {
+impl<D: Drift> HistoryBuilder<D> {
     pub fn time(mut self, time: bool) -> Self {
         self.time = time;
         self
@@ -40,9 +42,16 @@ impl HistoryBuilder {
         self
     }
 
-    pub fn gamma(mut self, gamma: f64) -> Self {
-        self.gamma = gamma;
-        self
+    pub fn drift<D2: Drift>(self, drift: D2) -> HistoryBuilder<D2> {
+        HistoryBuilder {
+            drift,
+            time: self.time,
+            mu: self.mu,
+            sigma: self.sigma,
+            beta: self.beta,
+            p_draw: self.p_draw,
+            online: self.online,
+        }
     }
 
     pub fn p_draw(mut self, p_draw: f64) -> Self {
@@ -55,7 +64,7 @@ impl HistoryBuilder {
         self
     }
 
-    pub fn build(self) -> History {
+    pub fn build(self) -> History<D> {
         History {
             size: 0,
             batches: Vec::new(),
@@ -64,41 +73,48 @@ impl HistoryBuilder {
             mu: self.mu,
             sigma: self.sigma,
             beta: self.beta,
-            gamma: self.gamma,
+            drift: self.drift,
             p_draw: self.p_draw,
             online: self.online,
         }
     }
 }
 
-impl Default for HistoryBuilder {
+impl HistoryBuilder<ConstantDrift> {
+    pub fn gamma(mut self, gamma: f64) -> Self {
+        self.drift = ConstantDrift(gamma);
+        self
+    }
+}
+
+impl Default for HistoryBuilder<ConstantDrift> {
     fn default() -> Self {
         Self {
             time: true,
             mu: MU,
             sigma: SIGMA,
             beta: BETA,
-            gamma: GAMMA,
+            drift: ConstantDrift(GAMMA),
             p_draw: P_DRAW,
             online: false,
         }
     }
 }
 
-pub struct History {
+pub struct History<D: Drift = ConstantDrift> {
     size: usize,
     pub(crate) batches: Vec<Batch>,
-    agents: HashMap<Index, Agent>,
+    agents: HashMap<Index, Agent<D>>,
     time: bool,
     mu: f64,
     sigma: f64,
     beta: f64,
-    gamma: f64,
+    drift: D,
     p_draw: f64,
     online: bool,
 }
 
-impl Default for History {
+impl Default for History<ConstantDrift> {
     fn default() -> Self {
         Self {
             size: 0,
@@ -108,18 +124,20 @@ impl Default for History {
             mu: MU,
             sigma: SIGMA,
             beta: BETA,
-            gamma: GAMMA,
+            drift: ConstantDrift(GAMMA),
             p_draw: P_DRAW,
             online: false,
         }
     }
 }
 
-impl History {
-    pub fn builder() -> HistoryBuilder {
+impl History<ConstantDrift> {
+    pub fn builder() -> HistoryBuilder<ConstantDrift> {
         HistoryBuilder::default()
     }
+}
 
+impl<D: Drift> History<D> {
     fn iteration(&mut self) -> (f64, f64) {
         let mut step = (0.0, 0.0);
 
@@ -247,7 +265,7 @@ impl History {
         results: Vec<Vec<f64>>,
         times: Vec<i64>,
         weights: Vec<Vec<Vec<f64>>>,
-        priors: HashMap<Index, Player>,
+        mut priors: HashMap<Index, Player<D>>,
     ) {
         assert!(times.is_empty() || self.time, "length(times)>0 but !h.time");
         assert!(
@@ -267,33 +285,35 @@ impl History {
             "(length(weights) > 0) & (length(composition) != length(weights))"
         );
 
-        let this_agent = composition
-            .iter()
-            .flatten()
-            .flatten()
-            .cloned()
-            .collect::<HashSet<_>>();
+        agent::clean(self.agents.values_mut(), true);
+
+        let mut this_agent = Vec::with_capacity(1024);
+
+        for agent in composition.iter().flatten().flatten() {
+            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,
-                                self.gamma,
+                                self.drift,
                             )
                         }),
-                        ..Default::default()
+                        message: N_INF,
+                        last_time: i64::MIN,
                     },
                 );
             }
         }
 
-        agent::clean(self.agents.values_mut(), true);
-
         let n = composition.len();
         let o = if self.time {
             sort_time(&times, false)
@@ -315,26 +335,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
-                    .iter()
-                    .filter(|&agent| b.skills.contains_key(agent))
-                    .cloned()
-                    .collect::<Vec<_>>();
+                // TODO: Is it faster to iterate over agents in batch instead?
+                for agent_idx in &this_agent {
+                    if let Some(skill) = batch.skills.get_mut(agent_idx) {
+                        skill.elapsed =
+                            batch::compute_elapsed(self.agents[agent_idx].last_time, batch.time);
 
-                for agent in &intersect {
-                    b.skills.get_mut(agent).unwrap().elapsed =
-                        batch::compute_elapsed(self.agents[agent].last_time, b.time);
+                        let agent = self.agents.get_mut(agent_idx).unwrap();
 
-                    let a = self.agents.get_mut(agent).unwrap();
-
-                    a.last_time = if self.time { b.time } else { i64::MAX };
-                    a.message = b.forward_prior_out(agent);
+                        agent.last_time = if self.time { batch.time } else { i64::MAX };
+                        agent.message = batch.forward_prior_out(agent_idx);
+                    }
                 }
 
                 k += 1;
@@ -343,6 +360,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 +374,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 a in b.skills.keys() {
-                    let agent = self.agents.get_mut(a).unwrap();
+                for agent_idx in batch.skills.keys() {
+                    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);
                 }
             } 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() {
-                    let agent = self.agents.get_mut(a).unwrap();
+                for agent_idx in batch.skills.keys() {
+                    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 +405,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
-                .iter()
-                .filter(|&agent| b.skills.contains_key(agent))
-                .cloned()
-                .collect::<Vec<_>>();
+            // TODO: Is it faster to iterate over agents in batch instead?
+            for agent_idx in &this_agent {
+                if let Some(skill) = batch.skills.get_mut(agent_idx) {
+                    skill.elapsed =
+                        batch::compute_elapsed(self.agents[agent_idx].last_time, batch.time);
 
-            for agent in &intersect {
-                b.skills.get_mut(agent).unwrap().elapsed =
-                    batch::compute_elapsed(self.agents[agent].last_time, b.time);
+                    let agent = self.agents.get_mut(agent_idx).unwrap();
 
-                let a = self.agents.get_mut(agent).unwrap();
-
-                a.last_time = if self.time { b.time } else { i64::MAX };
-                a.message = b.forward_prior_out(agent);
+                    agent.last_time = if self.time { batch.time } else { i64::MAX };
+                    agent.message = batch.forward_prior_out(agent_idx);
+                }
             }
 
             k += 1;
@@ -418,7 +433,7 @@ impl History {
 mod tests {
     use approx::assert_ulps_eq;
 
-    use crate::{Game, Gaussian, IndexMap, Player, EPSILON, ITERATIONS, P_DRAW};
+    use crate::{ConstantDrift, EPSILON, Game, Gaussian, ITERATIONS, IndexMap, P_DRAW, Player};
 
     use super::*;
 
@@ -445,7 +460,7 @@ mod tests {
                 Player::new(
                     Gaussian::from_ms(25.0, 25.0 / 3.0),
                     25.0 / 6.0,
-                    0.15 * 25.0 / 3.0,
+                    ConstantDrift(0.15 * 25.0 / 3.0),
                 ),
             );
         }
@@ -470,10 +485,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],
-            vec![],
+            &[0.0, 1.0],
+            &w,
             P_DRAW,
         )
         .posteriors();
@@ -506,7 +522,7 @@ mod tests {
                 Player::new(
                     Gaussian::from_ms(25.0, 25.0 / 3.0),
                     25.0 / 6.0,
-                    0.15 * 25.0 / 3.0,
+                    ConstantDrift(0.15 * 25.0 / 3.0),
                 ),
             );
         }
@@ -555,7 +571,7 @@ mod tests {
                 Player::new(
                     Gaussian::from_ms(25.0, 25.0 / 3.0),
                     25.0 / 6.0,
-                    25.0 / 300.0,
+                    ConstantDrift(25.0 / 300.0),
                 ),
             );
         }
@@ -613,7 +629,7 @@ mod tests {
                 Player::new(
                     Gaussian::from_ms(25.0, 25.0 / 3.0),
                     25.0 / 6.0,
-                    25.0 / 300.0,
+                    ConstantDrift(25.0 / 300.0),
                 ),
             );
         }

src/lib.rs

@@ -8,16 +8,19 @@ pub mod agent;
 #[cfg(feature = "approx")]
 mod approx;
 pub mod batch;
+pub mod drift;
 mod game;
 pub mod gaussian;
-// mod gaussian2;
 mod history;
+mod matrix;
 mod message;
 pub mod player;
 
+pub use drift::{ConstantDrift, Drift};
 pub use game::Game;
 pub use gaussian::Gaussian;
 pub use history::History;
+use matrix::Matrix;
 use message::DiffMessage;
 pub use player::Player;
 
@@ -81,7 +84,7 @@ where
     pub fn key(&self, idx: Index) -> Option<&K> {
         self.0
             .iter()
-            .find(|(_, &value)| value == idx)
+            .find(|&(_, value)| *value == idx)
             .map(|(key, _)| key)
     }
 
@@ -114,11 +117,7 @@ fn erfc(x: f64) -> f64 {
 
     let r = t * (-z * z - 1.26551223 + t * h).exp();
 
-    if x >= 0.0 {
-        r
-    } else {
-        2.0 - r
-    }
+    if x >= 0.0 { r } else { 2.0 - r }
 }
 
 fn erfc_inv(mut y: f64) -> f64 {
@@ -146,11 +145,7 @@ fn erfc_inv(mut y: f64) -> f64 {
         x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err)
     }
 
-    if y < 1.0 {
-        x
-    } else {
-        -x
-    }
+    if y < 1.0 { x } else { -x }
 }
 
 fn ppf(p: f64, mu: f64, sigma: f64) -> f64 {
@@ -238,9 +233,9 @@ pub(crate) fn sort_time(xs: &[i64], reverse: bool) -> Vec<usize> {
     let mut x = xs.iter().enumerate().collect::<Vec<_>>();
 
     if reverse {
-        x.sort_by_key(|(_, &x)| Reverse(x));
+        x.sort_by_key(|&(_, x)| Reverse(x));
     } else {
-        x.sort_by_key(|(_, &x)| x);
+        x.sort_by_key(|&(_, x)| x);
     }
 
     x.into_iter().map(|(i, _)| i).collect()
@@ -255,8 +250,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 +327,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

@@ -1,5 +1,4 @@
-use crate::gaussian::Gaussian;
-use crate::N_INF;
+use crate::{N_INF, gaussian::Gaussian};
 
 pub(crate) struct TeamMessage {
     pub(crate) prior: Gaussian,
@@ -15,14 +14,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
     }

src/player.rs

@@ -1,35 +1,32 @@
-use crate::{gaussian::Gaussian, BETA, GAMMA};
+use crate::{
+    BETA, GAMMA,
+    drift::{ConstantDrift, Drift},
+    gaussian::Gaussian,
+};
 
 #[derive(Clone, Copy, Debug)]
-pub struct Player {
+pub struct Player<D: Drift = ConstantDrift> {
     pub(crate) prior: Gaussian,
     pub(crate) beta: f64,
-    pub(crate) gamma: f64,
-    // pub(crate) draw: Gaussian,
+    pub(crate) drift: D,
 }
 
-impl Player {
-    pub fn new(prior: Gaussian, beta: f64, gamma: f64) -> Self {
-        Self {
-            prior,
-            beta,
-            gamma,
-            // draw: N_INF,
-        }
+impl<D: Drift> Player<D> {
+    pub fn new(prior: Gaussian, beta: f64, drift: D) -> Self {
+        Self { prior, beta, drift }
     }
 
     pub(crate) fn performance(&self) -> Gaussian {
-        self.prior.forget(self.beta, 1)
+        self.prior.forget(self.beta.powi(2))
     }
 }
 
-impl Default for Player {
+impl Default for Player<ConstantDrift> {
     fn default() -> Self {
         Self {
             prior: Gaussian::default(),
             beta: BETA,
-            gamma: GAMMA,
-            // draw: N_INF,
+            drift: ConstantDrift(GAMMA),
         }
     }
 }