docs: refresh README and CLAUDE.md; add ingest benchmark

The CLAUDE.md architecture section still described the pre-redesign engine:
its data flow named `Batch`, `Agent`, `Player` and `message.rs`, none of
which have existed since T2, and the public API it listed did not match
`lib.rs`. It is the first thing a fresh session reads, so it was actively
misleading. Rewritten against the current module layout, with the invariants
that are easy to violate — ties needing a positive `p_draw`, NaN never being
convergence, log-space evidence, color contiguity, `forbid(unsafe_code)`,
and ingestion-order equivalence — written down.

The README Todo list had five entries that were already done, including
"Time needs to be an enum": `Time` has been a trait since T2, and the
`batch::compute_elapsed()` it pointed at no longer exists. The genuinely
open item — cross-checking `quality()` against sublee/trueskill — stays.

`benches/ingest.rs` measures one-event-per-call against a single batched
call. The rest of the suite only measured batched construction, which is why
the quadratic fixed earlier on this branch went unnoticed for so long.

`TimeSlice::log_evidence` also hashes its target set once instead of
scanning the slice per player per event, so `log_evidence_for` with many
keys is no longer quadratic.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DnsaJg74eNSva3PJjK2eej

.gitea/workflows/ci.yml

@@ -0,0 +1,87 @@
+name: CI
+
+on:
+  push:
+    branches: [main]
+  pull_request:
+
+env:
+  CARGO_TERM_COLOR: always
+  RUSTFLAGS: -D warnings
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          # The build most consumers get.
+          - name: default
+            features: ""
+            profile: ""
+          # Most numerical goldens need `approx` for assert_ulps_eq.
+          - name: approx
+            features: "--features approx"
+            profile: ""
+          # The parallel path, including tests/determinism.rs.
+          - name: rayon
+            features: "--features approx,rayon"
+            profile: ""
+          # Critical: debug_assert! is compiled out here, which is where the
+          # tie/p_draw and score_sigma validation actually has to hold.
+          - name: release
+            features: "--features approx"
+            profile: "--release"
+    name: test (${{ matrix.name }})
+    steps:
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@stable
+      - uses: Swatinem/rust-cache@v2
+      - run: cargo test ${{ matrix.profile }} ${{ matrix.features }}
+      - run: cargo test ${{ matrix.profile }} ${{ matrix.features }} --doc
+
+  determinism:
+    name: determinism across thread counts
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@stable
+      - uses: Swatinem/rust-cache@v2
+      # Posteriors must be bit-identical regardless of how many rayon workers
+      # run the color-group sweep.
+      - run: |
+          for threads in 1 2 4 8; do
+            echo "== RAYON_NUM_THREADS=$threads =="
+            RAYON_NUM_THREADS=$threads cargo test --release \
+              --features approx,rayon --test determinism
+          done
+
+  lint:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@stable
+        with:
+          components: clippy
+      - uses: Swatinem/rust-cache@v2
+      - run: cargo clippy --all-targets --all-features -- -D warnings
+
+  format:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      # rustfmt.toml uses nightly-only options (imports_granularity).
+      - uses: dtolnay/rust-toolchain@nightly
+        with:
+          components: rustfmt
+      - run: cargo +nightly fmt --check
+
+  msrv:
+    name: minimum supported Rust version
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: dtolnay/rust-toolchain@1.85.0
+      - uses: Swatinem/rust-cache@v2
+      - run: cargo check --all-targets --features approx,rayon

CLAUDE.md

@@ -5,42 +5,96 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 ## 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)
+just test                            # Full suite across every feature combination CI checks
+just check                           # Fast inner loop: cargo test --features approx
+just lint                            # clippy, warnings denied
+just fmt                             # ALWAYS nightly — rustfmt.toml uses nightly-only options
+just determinism                     # Bit-identical posteriors at RAYON_NUM_THREADS 1/2/4/8
+just ci                              # Everything CI runs
+cargo test --lib <test_name>         # A single test by name
+cargo bench                          # Criterion benchmarks
 ```
 
-The `approx` feature enables `approx::AbsDiffEq` for `Gaussian`:
-```bash
-cargo test --features approx
-```
+**Run tests in release too.** `debug_assert!` is compiled out there, and that
+is where several defects have hidden — a debug-only run is not evidence.
+`just test` includes a release job.
+
+### Feature flags
+
+- `approx` — `approx::AbsDiffEq` etc. for `Gaussian`. Most numerical goldens need it.
+- `rayon` — opt-in parallel within-slice sweep and per-slice query passes.
 
 ## 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.
+A Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py):
+Bayesian skill rating that infers skill at every point in time, propagating
+evidence both forward and backward across a history.
 
 ### Data flow
 
 ```
-History  →  Batch[]  →  Game[]  →  teams/players
+History  →  TimeSlice[]  →  Event[]  →  Team[] → Item[]
+                                  ↓
+                     Game (factor graph) → Schedule → BuiltinFactor[]
 ```
 
-- **`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`.
-- **`MarginFactor`** (`factor/margin.rs`) — Gaussian observation factor on a diff variable; engaged by `Outcome::Scored`.
-- **`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`.
+- **`History`** (`history.rs`) — top level. Interns keys, groups events into
+  `TimeSlice`s by time, runs the forward/backward sweep in `converge()`, and
+  answers `learning_curves()`, `current_skill()`, `log_evidence()`,
+  `predict_quality()`, `predict_outcome()`. Built via `HistoryBuilder`.
+- **`TimeSlice`** (`time_slice.rs`) — all events at one time. Owns a
+  `SkillStore` and a `ScratchArena`; `iteration()` sweeps its events, using
+  `ColorGroups` to partition independent ones.
+- **`Event`** (`time_slice.rs`) — one match. `compute()` runs inference reading
+  skills immutably; `apply()` folds the result back. The split is what lets a
+  color group run in parallel with no `unsafe`.
+- **`Game`** (`game.rs`) — a single match's factor graph. `run_chain` builds the
+  diff chain between rank-adjacent teams and drives it to convergence.
+- **`Gaussian`** (`gaussian.rs`) — natural parameters (`pi = 1/sigma²`,
+  `tau = mu/sigma²`). `Mul`/`Div` are the EP product/cavity: pure adds and
+  subtracts. Variance-space ops (`Add`, `Sub`, `exclude`, `forget`) go through
+  `from_mv`/`variance()` and take no square root.
+- **`factor/`** — `TeamSumFactor`, `RankDiffFactor`, `TruncFactor` (ranked),
+  `MarginFactor` (scored), over a flat `VarStore`. `BuiltinFactor` dispatches
+  by enum rather than `dyn`.
+- **`Schedule`** (`schedule.rs`) — drives factor propagation. `EpsilonOrMax` is
+  the only implementation.
+- **`Competitor`** (`competitor.rs`) — per-history temporal state (`message`,
+  `last_time`). **`Rating`** (`rating.rs`) — static config (prior, `beta`, drift).
+- **`storage/`** — `SkillStore` (per slice) and `CompetitorStore` (per history),
+  both dense `Vec`s indexed by `Index`.
+- **`KeyTable`** (`key_table.rs`) — user key ↔ `Index`, both directions O(1).
+- **`Drift`** (`drift.rs`) / **`Time`** (`time.rs`) — traits. `Time` is a *trait*
+  (`i64`, `Untimed`), not an enum.
+- **`lib.rs`** — public exports, global defaults (`MU`, `SIGMA`, `BETA`,
+  `GAMMA`, `P_DRAW`, `EPSILON`, `ITERATIONS`), and the standalone `quality()`,
+  `cdf()`, `erfc()`.
 
-### Key design points
+### Invariants worth knowing
 
-- `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.
+- **A tie needs `p_draw > 0`.** With `p_draw == 0.0` the truncation margin is
+  zero and the two-sided tie update evaluates `0/0`. Ingestion rejects such
+  events with `InferenceError::TieWithoutDrawProbability`. This includes
+  `Outcome::winner(w, n)` for `n >= 3`, which ties every loser.
+- **NaN is never convergence.** Comparisons against NaN are all false, so
+  `tuple_gt` reads NaN as "below epsilon". Use `step_converged` /
+  `step_is_finite`, never `!tuple_gt(..)` alone.
+- **Evidence accumulates in log space.** A linear product over a long diff
+  chain underflows to zero, and `ln(0)` is `-inf`.
+- **Colors are contiguous.** `recompute_color_groups` reorders events so each
+  color occupies one range; `ColorGroups::groups_are_contiguous` asserts it.
+- **The crate is `#![forbid(unsafe_code)]`.** Keep it that way.
+- **Ingestion order must not change the answer.** Events added one at a time
+  must converge to the same fixed point as the same events batched — see
+  `tests/ingestion_equivalence.rs`.
+
+### Testing notes
+
+- Numerical goldens are cross-validated against the Python/Julia reference.
+  Some are *convergence residuals*, not exact values; treat a small movement
+  as suspicious but check whether the new value is closer to the analytic
+  truth (symmetric fixtures converge to their prior mean exactly) before
+  assuming a regression.
+- `tests/degenerate_inputs.rs` covers empty/boundary/error paths,
+  `tests/ingestion_equivalence.rs` covers batching order, `tests/quality.rs`
+  covers N-group quality, `tests/determinism.rs` covers thread counts.

Cargo.toml

@@ -2,6 +2,17 @@
 name = "trueskill-tt"
 version = "0.1.2"
 edition = "2024"
+rust-version = "1.85"
+description = "TrueSkill Through Time: Bayesian skill rating that tracks how skill evolves over time, via Gaussian message passing"
+repository = "https://git.aceofba.se/logaritmisk/trueskill-tt"
+readme = "README.md"
+keywords = ["trueskill", "rating", "bayesian", "elo", "skill"]
+categories = ["algorithms", "science", "game-development"]
+# TODO: pick a licence before publishing. `cargo publish` rejects a crate
+# without `license` (or `license-file`), and without one the source carries no
+# stated terms. The Rust convention is `license = "MIT OR Apache-2.0"` plus the
+# matching LICENSE-MIT / LICENSE-APACHE files.
+exclude = ["/docs", "/benches/*.txt", "/temp", "/.gitea"]
 
 [lib]
 bench = false
@@ -22,6 +33,10 @@ harness = false
 name = "scored"
 harness = false
 
+[[bench]]
+name = "ingest"
+harness = false
+
 [dependencies]
 approx = { version = "0.5.1", optional = true }
 rayon = { version = "1", optional = true }

Justfile

@@ -1,4 +1,39 @@
 alias b := bench
+alias t := test
+
+# Run the full test suite across the feature combinations CI checks.
+test:
+    cargo test
+    cargo test --features approx
+    cargo test --features approx,rayon
+    cargo test --release --features approx
+
+# Fast inner-loop tests.
+check:
+    cargo test --features approx
+
+# Posteriors must be bit-identical across rayon worker counts.
+determinism:
+    #!/usr/bin/env bash
+    set -euo pipefail
+    for threads in 1 2 4 8; do
+        echo "== RAYON_NUM_THREADS=$threads =="
+        RAYON_NUM_THREADS=$threads cargo test --release \
+            --features approx,rayon --test determinism
+    done
+
+lint:
+    cargo clippy --all-targets --all-features -- -D warnings
+
+# Always nightly: rustfmt.toml uses nightly-only options.
+fmt:
+    cargo +nightly fmt
+
+fmt-check:
+    cargo +nightly fmt --check
+
+# Everything CI runs.
+ci: fmt-check lint test determinism
 
 store:
     cargo bench -- --save-baseline base

README.md

@@ -96,8 +96,8 @@ h.converge().unwrap();
 
 - [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))
-- [ ] Add Observer (see [argmin](https://docs.rs/argmin/latest/argmin/core/trait.Observe.html) for inspiration)
+- [x] Generalise a time axis — `Time` is now a trait (`Untimed`, `i64`), not an enum
+- [x] Add examples (`examples/atp.rs`, `examples/scored.rs`)
+- [x] Add Observer (`Observer` / `NullObserver`)
+- [x] Benchmark the inference loop (`benches/batch.rs`, `benches/history_converge.rs`, `benches/ingest.rs`)
+- [ ] Cross-check `quality()` against [sublee/trueskill](https://github.com/sublee/trueskill/tree/master) — N-group support works and is covered by invariants, but no reference values are asserted

benches/ingest.rs

@@ -0,0 +1,62 @@
+//! Ingestion cost: one event per call versus one batched call.
+//!
+//! The rest of the suite only measures batched construction, which is why a
+//! quadratic in the incremental path went unnoticed — `record_winner` and
+//! `event(..).commit()` each ingest a single event, so a caller looping over a
+//! match feed takes that path.
+
+use std::hint::black_box;
+
+use criterion::{BenchmarkId, Criterion, criterion_group, criterion_main};
+use smallvec::smallvec;
+use trueskill_tt::{Event, History, Member, Outcome, Team};
+
+fn events(n: usize, time: i64) -> Vec<Event<i64, String>> {
+    (0..n)
+        .map(|i| Event {
+            time,
+            teams: smallvec![
+                Team::with_members([Member::new(format!("p{}", 2 * i))]),
+                Team::with_members([Member::new(format!("p{}", 2 * i + 1))]),
+            ],
+            outcome: Outcome::winner(0, 2),
+        })
+        .collect()
+}
+
+fn bench_ingest(c: &mut Criterion) {
+    let mut group = c.benchmark_group("ingest");
+
+    for n in [250usize, 500, 1000] {
+        group.bench_with_input(BenchmarkId::new("one-at-a-time", n), &n, |b, &n| {
+            b.iter_batched(
+                || events(n, 0),
+                |evs| {
+                    let mut h: History<i64, _, _, String> = History::builder_with_key().build();
+                    for ev in evs {
+                        h.add_events(std::iter::once(ev)).unwrap();
+                    }
+                    black_box(h.time_slices_len())
+                },
+                criterion::BatchSize::SmallInput,
+            );
+        });
+
+        group.bench_with_input(BenchmarkId::new("single-batch", n), &n, |b, &n| {
+            b.iter_batched(
+                || events(n, 0),
+                |evs| {
+                    let mut h: History<i64, _, _, String> = History::builder_with_key().build();
+                    h.add_events(evs).unwrap();
+                    black_box(h.time_slices_len())
+                },
+                criterion::BatchSize::SmallInput,
+            );
+        });
+    }
+
+    group.finish();
+}
+
+criterion_group!(benches, bench_ingest);
+criterion_main!(benches);

src/color_group.rs

@@ -26,39 +26,75 @@ pub(crate) struct ColorGroups {
 }
 
 impl ColorGroups {
-    #[allow(dead_code)]
     pub(crate) fn new() -> Self {
         Self::default()
     }
 
-    #[allow(dead_code)]
-    pub(crate) fn n_colors(&self) -> usize {
-        self.groups.len()
-    }
-
-    #[allow(dead_code)]
     pub(crate) fn is_empty(&self) -> bool {
         self.groups.is_empty()
     }
 
-    /// Total event count across all colors.
-    #[allow(dead_code)]
+    /// Number of distinct colors in the partition. Test-only.
+    #[cfg(test)]
+    pub(crate) fn n_colors(&self) -> usize {
+        self.groups.len()
+    }
+
+    /// Total event count across all colors. Test-only.
+    #[cfg(test)]
     pub(crate) fn total_events(&self) -> usize {
         self.groups.iter().map(|g| g.len()).sum()
     }
 
     /// Contiguous index range for one color after events have been reordered
     /// into color-contiguous positions by `TimeSlice::recompute_color_groups`.
-    #[allow(dead_code)]
     pub(crate) fn color_range(&self, color_idx: usize) -> std::ops::Range<usize> {
         let group = &self.groups[color_idx];
         if group.is_empty() {
             return 0..0;
         }
+
         let start = *group.first().unwrap();
         let end = *group.last().unwrap() + 1;
+
+        debug_assert_eq!(
+            end - start,
+            group.len(),
+            "color {color_idx} is not contiguous; its range would overlap other colors"
+        );
+
         start..end
     }
+
+    /// Whether every color occupies a contiguous, ascending range of event
+    /// indices, and no two colors overlap.
+    ///
+    /// The parallel sweep derives one `&mut` sub-slice per color from these
+    /// ranges and relies on them being disjoint. That disjointness is what
+    /// makes concurrent writes to distinct skills sound, so it is checked
+    /// rather than assumed.
+    pub(crate) fn groups_are_contiguous(&self) -> bool {
+        let mut expected_start = 0;
+
+        for group in &self.groups {
+            if group.is_empty() {
+                continue;
+            }
+
+            let ascending_run = group
+                .iter()
+                .enumerate()
+                .all(|(offset, &idx)| idx == group[0] + offset);
+
+            if !ascending_run || group[0] != expected_start {
+                return false;
+            }
+
+            expected_start += group.len();
+        }
+
+        true
+    }
 }
 
 /// Compute color groups greedily.

src/error.rs

@@ -1,6 +1,7 @@
 use std::fmt;
 
 #[derive(Debug, Clone, PartialEq)]
+#[non_exhaustive]
 pub enum InferenceError {
     /// Expected and actual lengths of some array-shaped input differ.
     MismatchedShape {
@@ -8,15 +9,35 @@ pub enum InferenceError {
         expected: usize,
         got: usize,
     },
+    /// An `Outcome` of the wrong variant was supplied for the requested inference.
+    WrongOutcomeKind {
+        context: &'static str,
+        expected: &'static str,
+        got: &'static str,
+    },
     /// A probability value is outside `[0, 1]`.
     InvalidProbability { value: f64 },
     /// A scalar parameter is outside its valid range.
     InvalidParameter { name: &'static str, value: f64 },
+    /// An event contains tied teams, but the draw probability is zero.
+    ///
+    /// A zero draw probability asserts that draws cannot occur, so a tied
+    /// result has no representable likelihood. Configure a positive `p_draw`
+    /// (via `HistoryBuilder::p_draw` or `GameOptions::p_draw`) to admit ties.
+    TieWithoutDrawProbability { teams: (usize, usize) },
     /// Convergence exceeded `max_iter` without falling below `epsilon`.
     ConvergenceFailed {
         last_step: (f64, f64),
         iterations: usize,
     },
+    /// Inference produced a non-finite value (NaN or infinity).
+    ///
+    /// Indicates numerical breakdown; the resulting skills are meaningless
+    /// and must not be treated as a converged estimate.
+    NonFiniteResult {
+        context: &'static str,
+        step: (f64, f64),
+    },
     /// Negative precision: a Gaussian with `pi < 0` slipped into an API call.
     NegativePrecision { pi: f64 },
 }
@@ -31,9 +52,29 @@ impl fmt::Display for InferenceError {
             } => {
                 write!(f, "{kind}: expected length {expected}, got {got}")
             }
+            Self::WrongOutcomeKind {
+                context,
+                expected,
+                got,
+            } => {
+                write!(f, "{context}: expected {expected}, got {got}")
+            }
             Self::InvalidProbability { value } => {
                 write!(f, "probability must be in [0, 1]; got {value}")
             }
+            Self::TieWithoutDrawProbability { teams } => {
+                write!(
+                    f,
+                    "teams {} and {} are tied, but p_draw is 0.0; set a positive draw probability to admit ties",
+                    teams.0, teams.1
+                )
+            }
+            Self::NonFiniteResult { context, step } => {
+                write!(
+                    f,
+                    "{context}: inference produced a non-finite result (step = {step:?})"
+                )
+            }
             Self::InvalidParameter { name, value } => {
                 write!(f, "{name} is invalid: {value}")
             }

src/factor/margin.rs

@@ -64,9 +64,13 @@ impl Factor for MarginFactor {
     }
 }
 
+/// Density of the observed margin under the cavity, clamped to a positive
+/// floor so a far-out observation cannot underflow to `0.0` and make
+/// `log_evidence` `-inf`.
 fn cavity_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 {
     let combined_sigma = (cavity.sigma().powi(2) + sigma.powi(2)).sqrt();
-    pdf(m_obs, cavity.mu(), combined_sigma)
+
+    pdf(m_obs, cavity.mu(), combined_sigma).max(f64::MIN_POSITIVE)
 }
 
 #[cfg(test)]

src/factor/trunc.rs

@@ -72,12 +72,20 @@ impl Factor for TruncFactor {
 }
 
 /// P(diff > margin) for non-tie, P(|diff| < margin) for tie.
+///
+/// Clamped to a positive floor: for a near-certain outcome the tail rounds to
+/// exactly 0.0, and the `erfc` approximation used by `cdf` carries ~1e-7 error
+/// so it can even return slightly more than 1.0, making the difference
+/// negative. Either would send `log_evidence` to `-inf` or NaN and poison the
+/// sum across the whole history.
 fn cavity_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
-    if tie {
+    let raw = if tie {
         cdf(margin, diff.mu(), diff.sigma()) - cdf(-margin, diff.mu(), diff.sigma())
     } else {
         1.0 - cdf(margin, diff.mu(), diff.sigma())
-    }
+    };
+
+    raw.clamp(f64::MIN_POSITIVE, 1.0)
 }
 
 #[cfg(test)]

src/game.rs

@@ -37,10 +37,17 @@ impl DiffFactor {
         }
     }
 
-    pub(crate) fn evidence(&self) -> f64 {
+    /// Log of this link's cached evidence.
+    ///
+    /// Accumulating in log space keeps a long diff chain from underflowing:
+    /// each link contributes a probability in `(0, 1]`, so the linear product
+    /// over an n-team game decays geometrically and flushes to zero — and
+    /// `ln(0.0)` is `-inf` — well within the team counts a large free-for-all
+    /// reaches.
+    pub(crate) fn log_evidence(&self) -> f64 {
         match self {
-            Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0),
-            Self::Margin(f) => f.evidence_cached.unwrap_or(1.0),
+            Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0).ln(),
+            Self::Margin(f) => f.evidence_cached.unwrap_or(1.0).ln(),
         }
     }
 
@@ -81,18 +88,14 @@ impl Default for GameOptions {
 /// Owned variant of `Game` returned by public constructors.
 ///
 /// Unlike `Game<'a, T, D>` (which borrows its result/weights slices from
-/// History's internal state), `OwnedGame<T, D>` owns its inputs so it can
-/// be returned freely from public constructors.
+/// History's internal state), `OwnedGame<T, D>` owns the team ratings, so it
+/// can be returned freely from public constructors. The inference inputs
+/// themselves are not retained — nothing reads them back.
 #[derive(Debug)]
-#[allow(dead_code)]
 pub struct OwnedGame<T: Time, D: Drift<T>> {
     teams: Vec<Vec<Rating<T, D>>>,
-    result: Vec<f64>,
-    weights: Vec<Vec<f64>>,
-    p_draw: f64,
-    pub(crate) convergence: crate::ConvergenceOptions,
     pub(crate) likelihoods: Vec<Vec<Gaussian>>,
-    pub(crate) evidence: f64,
+    pub(crate) log_evidence: f64,
 }
 
 impl<T: Time, D: Drift<T>> OwnedGame<T, D> {
@@ -112,16 +115,10 @@ impl<T: Time, D: Drift<T>> OwnedGame<T, D> {
             convergence,
             &mut arena,
         );
-        let likelihoods = g.likelihoods;
-        let evidence = g.evidence;
         Self {
             teams,
-            result,
-            weights,
-            p_draw,
-            convergence,
-            likelihoods,
-            evidence,
+            likelihoods: g.likelihoods,
+            log_evidence: g.log_evidence,
         }
     }
 
@@ -141,16 +138,10 @@ impl<T: Time, D: Drift<T>> OwnedGame<T, D> {
             convergence,
             &mut arena,
         );
-        let likelihoods = g.likelihoods;
-        let evidence = g.evidence;
         Self {
             teams,
-            result: scores,
-            weights,
-            p_draw: 0.0,
-            convergence,
-            likelihoods,
-            evidence,
+            likelihoods: g.likelihoods,
+            log_evidence: g.log_evidence,
         }
     }
 
@@ -163,7 +154,7 @@ impl<T: Time, D: Drift<T>> OwnedGame<T, D> {
     }
 
     pub fn log_evidence(&self) -> f64 {
-        self.evidence.ln()
+        self.log_evidence
     }
 }
 
@@ -175,7 +166,7 @@ pub struct Game<'a, T: Time = i64, D: Drift<T> = crate::drift::ConstantDrift> {
     p_draw: f64,
     pub(crate) convergence: crate::ConvergenceOptions,
     pub(crate) likelihoods: Vec<Vec<Gaussian>>,
-    pub(crate) evidence: f64,
+    pub(crate) log_evidence: f64,
 }
 
 impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
@@ -222,7 +213,7 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             p_draw,
             convergence,
             likelihoods: Vec::new(),
-            evidence: 0.0,
+            log_evidence: 0.0,
         };
 
         this.likelihoods(arena);
@@ -261,7 +252,7 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             p_draw: 0.0,
             convergence,
             likelihoods: Vec::new(),
-            evidence: 0.0,
+            log_evidence: 0.0,
         };
 
         this.likelihoods_scored(arena, score_sigma);
@@ -355,7 +346,7 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             arena.lhood_lose[n_teams - 1] = pw_last - links[n_diffs - 1].msg();
         }
 
-        let evidence: f64 = links.iter().map(|l| l.evidence()).product();
+        let log_evidence: f64 = links.iter().map(DiffFactor::log_evidence).sum();
 
         // Inverse permutation: inv_buf[orig_i] = sorted_i.
         arena.inv_buf.resize(n_teams, 0);
@@ -371,10 +362,9 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             .map(|(orig_i, (players, weights))| {
                 let si = arena.inv_buf[orig_i];
                 let m = arena.lhood_win[si] * arena.lhood_lose[si];
-                let performance = players
-                    .iter()
-                    .zip(weights.iter())
-                    .fold(N00, |p, (player, &w)| p + (player.performance() * w));
+                // Already folded into `team_prior` at the top of the chain,
+                // indexed by sorted position.
+                let performance = arena.team_prior[si];
                 players
                     .iter()
                     .zip(weights.iter())
@@ -386,11 +376,11 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             })
             .collect::<Vec<_>>();
 
-        (evidence, likelihoods)
+        (log_evidence, likelihoods)
     }
 
     fn likelihoods(&mut self, arena: &mut ScratchArena) {
-        let (evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
+        let (log_evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
             let tie = self.result[sort_buf[i]] == self.result[sort_buf[i + 1]];
             let margin = if self.p_draw == 0.0 {
                 0.0
@@ -405,17 +395,17 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
             let vid = vars.alloc(N_INF);
             DiffFactor::Trunc(TruncFactor::new(vid, margin, tie))
         });
-        self.evidence = evidence;
+        self.log_evidence = log_evidence;
         self.likelihoods = likelihoods;
     }
 
     fn likelihoods_scored(&mut self, arena: &mut ScratchArena, score_sigma: f64) {
-        let (evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
+        let (log_evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
             let m_obs = self.result[sort_buf[i]] - self.result[sort_buf[i + 1]];
             let vid = vars.alloc(N_INF);
             DiffFactor::Margin(MarginFactor::new(vid, m_obs, score_sigma))
         });
-        self.evidence = evidence;
+        self.log_evidence = log_evidence;
         self.likelihoods = likelihoods;
     }
 
@@ -433,7 +423,7 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
     }
 
     pub fn log_evidence(&self) -> f64 {
-        self.evidence.ln()
+        self.log_evidence
     }
 }
 
@@ -458,11 +448,22 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
 
         let ranks = outcome
             .as_ranks()
-            .ok_or(crate::InferenceError::MismatchedShape {
-                kind: "Game::ranked requires Outcome::Ranked",
-                expected: 0,
-                got: 0,
+            .ok_or(crate::InferenceError::WrongOutcomeKind {
+                context: "Game::ranked",
+                expected: "Outcome::Ranked",
+                got: "Outcome::Scored",
             })?;
+
+        let tied = if options.p_draw == 0.0 {
+            crate::first_tied_pair(ranks)
+        } else {
+            None
+        };
+
+        if let Some(teams) = tied {
+            return Err(crate::InferenceError::TieWithoutDrawProbability { teams });
+        }
+
         let max_rank = ranks.iter().copied().max().unwrap_or(0) as f64;
         let result: Vec<f64> = ranks.iter().map(|&r| max_rank - r as f64).collect();
         let teams_owned: Vec<Vec<Rating<T, D>>> = teams.iter().map(|t| t.to_vec()).collect();
@@ -497,10 +498,10 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
         }
         let scores = outcome
             .as_scores()
-            .ok_or(crate::InferenceError::MismatchedShape {
-                kind: "Game::scored requires Outcome::Scored",
-                expected: 0,
-                got: 0,
+            .ok_or(crate::InferenceError::WrongOutcomeKind {
+                context: "Game::scored",
+                expected: "Outcome::Scored",
+                got: "Outcome::Ranked",
             })?
             .to_vec();
         let teams_owned: Vec<Vec<Rating<T, D>>> = teams.iter().map(|t| t.to_vec()).collect();
@@ -730,8 +731,12 @@ mod tests {
         let a = p[0][0];
         let b = p[1][0];
 
-        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);
+        // Two identical competitors drawing must land on their shared prior
+        // mean exactly, by symmetry. The reference transcription of 24.999999
+        // is that value rounded to six decimals; asserting it at epsilon 1e-6
+        // left no headroom. The root-free variance path now hits 25.0 exactly.
+        assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6);
+        assert_ulps_eq!(b, Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6);
 
         let t_a = R::new(
             Gaussian::from_ms(25.0, 3.0),
@@ -1124,7 +1129,10 @@ mod tests {
             &GameOptions::default(),
         )
         .unwrap_err();
-        assert!(matches!(err, crate::InferenceError::MismatchedShape { .. }));
+        assert!(matches!(
+            err,
+            crate::InferenceError::WrongOutcomeKind { .. }
+        ));
     }
 
     #[test]

src/gaussian.rs

@@ -35,6 +35,28 @@ impl Gaussian {
         }
     }
 
+    /// Construct from mean and *variance*, skipping the square-root round trip.
+    ///
+    /// `from_ms(mu, var.sqrt())` immediately squares the root away again to
+    /// recover `pi = 1/var`. Variance-combining operations (`Add`, `Sub`,
+    /// `exclude`, `forget`) work in variance space throughout, so they go
+    /// through here instead and never take a root.
+    #[inline]
+    pub(crate) fn from_mv(mu: f64, var: f64) -> Self {
+        if var == f64::INFINITY {
+            Self { pi: 0.0, tau: 0.0 }
+        } else if var == 0.0 {
+            // Point mass at mu; see `from_ms` for the tau convention.
+            Self {
+                pi: f64::INFINITY,
+                tau: if mu == 0.0 { 0.0 } else { f64::INFINITY },
+            }
+        } else {
+            let pi = 1.0 / var;
+            Self { pi, tau: mu * pi }
+        }
+    }
+
     /// Construct directly from natural parameters.
     #[inline]
     pub(crate) const fn from_natural(pi: f64, tau: f64) -> Self {
@@ -64,6 +86,21 @@ impl Gaussian {
         }
     }
 
+    /// Variance, `1 / pi`, without the root-and-square of `sigma().powi(2)`.
+    ///
+    /// Mirrors `sigma()`'s treatment of the improper (`pi <= 0`) and point-mass
+    /// (`pi == inf`) cases.
+    #[inline]
+    pub(crate) fn variance(&self) -> f64 {
+        if self.pi <= 0.0 {
+            f64::INFINITY
+        } else if self.pi.is_infinite() {
+            0.0
+        } else {
+            1.0 / self.pi
+        }
+    }
+
     #[inline]
     pub fn sigma(&self) -> f64 {
         // A non-positive precision is improper → infinite standard deviation. Guarding
@@ -86,22 +123,21 @@ impl Gaussian {
     }
 
     pub(crate) fn exclude(&self, other: Gaussian) -> Self {
-        let var = self.sigma().powi(2) - other.sigma().powi(2);
+        let var = self.variance() - other.variance();
         if var <= 0.0 {
             // When sigma_self ≈ sigma_other (including ULP-level rounding differences
             // from the pi→sigma accessor round-trip), the excluded contribution is N00.
             // Computing from_ms(tiny_mu, 0.0) would give {pi:inf, tau:inf}, whose
             // mu() = inf/inf = NaN.  Returning N00 is correct: when both Gaussians
             // carry the same variance, the residual is a point mass at 0.
-            return Gaussian::from_ms(0.0, 0.0);
+            return Gaussian::from_mv(0.0, 0.0);
         }
-        let mu = self.mu() - other.mu();
-        Self::from_ms(mu, var.sqrt())
+
+        Self::from_mv(self.mu() - other.mu(), var)
     }
 
     pub(crate) fn forget(&self, variance_delta: f64) -> Self {
-        let var = self.sigma().powi(2) + variance_delta;
-        Self::from_ms(self.mu(), var.sqrt())
+        Self::from_mv(self.mu(), self.variance() + variance_delta)
     }
 
     /// EP damping in natural-parameter space: `α·new + (1−α)·self`.
@@ -128,9 +164,7 @@ impl ops::Add<Gaussian> for Gaussian {
     /// Variance addition: (mu1 + mu2, sqrt(σ1² + σ2²)).
     /// Used for combining performance and noise; rare relative to mul/div.
     fn add(self, rhs: Gaussian) -> Self::Output {
-        let mu = self.mu() + rhs.mu();
-        let var = self.sigma().powi(2) + rhs.sigma().powi(2);
-        Self::from_ms(mu, var.sqrt())
+        Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance())
     }
 }
 
@@ -138,9 +172,7 @@ impl ops::Sub<Gaussian> for Gaussian {
     type Output = Gaussian;
     /// (mu1 - mu2, sqrt(σ1² + σ2²)). Same sigma combination as Add.
     fn sub(self, rhs: Gaussian) -> Self::Output {
-        let mu = self.mu() - rhs.mu();
-        let var = self.sigma().powi(2) + rhs.sigma().powi(2);
-        Self::from_ms(mu, var.sqrt())
+        Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance())
     }
 }
 
@@ -161,7 +193,7 @@ impl ops::Mul<f64> for Gaussian {
         if scalar == 0.0 {
             // Scaling by 0 collapses to a point mass at 0 (sigma' = 0, mu' = 0).
             // This is N00, the additive identity, NOT N_INF.
-            return Gaussian::from_ms(0.0, 0.0);
+            return Gaussian::from_mv(0.0, 0.0);
         }
         // sigma' = sigma * |scalar|  =>  pi' = pi / scalar²
         // mu'    = mu * scalar       =>  tau' = tau / scalar

src/history.rs

@@ -69,7 +69,20 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> HistoryBuilder<
         }
     }
 
+    /// Probability that two evenly-matched sides draw.
+    ///
+    /// Must be in `[0.0, 1.0)`. A zero draw probability asserts that draws
+    /// cannot occur, so ingesting a tied outcome then fails with
+    /// `InferenceError::TieWithoutDrawProbability`.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `p_draw` is outside `[0.0, 1.0)` or is NaN.
     pub fn p_draw(mut self, p_draw: f64) -> Self {
+        assert!(
+            (0.0..1.0).contains(&p_draw),
+            "p_draw must be in [0.0, 1.0) (got {p_draw})"
+        );
         self.p_draw = p_draw;
         self
     }
@@ -79,6 +92,11 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> HistoryBuilder<
         self
     }
 
+    /// Default observation noise for scored outcomes.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `score_sigma` is not strictly positive.
     pub fn score_sigma(mut self, score_sigma: f64) -> Self {
         assert!(
             score_sigma > 0.0,
@@ -88,7 +106,24 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> HistoryBuilder<
         self
     }
 
+    /// Convergence tolerance, iteration cap, and EP damping.
+    ///
+    /// # Panics
+    ///
+    /// Panics if `alpha` is outside `(0.0, 1.0]`, or if `epsilon` is negative
+    /// or NaN. An `alpha` of zero would leave every EP update unapplied, so
+    /// inference would silently return the priors.
     pub fn convergence(mut self, opts: ConvergenceOptions) -> Self {
+        assert!(
+            opts.alpha > 0.0 && opts.alpha <= 1.0,
+            "convergence alpha must be in (0.0, 1.0] (got {})",
+            opts.alpha
+        );
+        assert!(
+            opts.epsilon >= 0.0,
+            "convergence epsilon must be non-negative (got {})",
+            opts.epsilon
+        );
         self.convergence = opts;
         self
     }
@@ -220,6 +255,10 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     fn iteration(&mut self) -> (f64, f64) {
         let mut step = (0.0, 0.0);
 
+        if self.time_slices.is_empty() {
+            return step;
+        }
+
         competitor::clean(self.agents.values_mut(), false);
 
         for j in (0..self.time_slices.len() - 1).rev() {
@@ -273,6 +312,12 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         step
     }
 
+    /// Number of distinct time slices in the history.
+    #[must_use]
+    pub fn time_slices_len(&self) -> usize {
+        self.time_slices.len()
+    }
+
     /// Learning curves for all competitors, keyed by their user-facing key.
     ///
     /// Note: `key(idx)` is O(n) per lookup; this method is therefore O(n²)
@@ -385,7 +430,15 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
 
     /// Draw-probability quality metric for the given teams (key slices).
     ///
-    /// Values range roughly [0, 1]; 1 == perfectly matched.
+    /// Values range roughly [0, 1]; 1 == perfectly matched. Supports any
+    /// number of teams.
+    ///
+    /// # Panics
+    ///
+    /// Panics if fewer than two teams are supplied, or if a team resolves to
+    /// no known competitors — keys absent from the history, or competitors
+    /// with no recorded skill, are dropped, so a team of entirely-unknown
+    /// keys becomes empty. Use `lookup` to check keys first.
     pub fn predict_quality(&self, teams: &[&[&K]]) -> f64 {
         let groups: Vec<Vec<Gaussian>> = teams
             .iter()
@@ -435,6 +488,18 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         use smallvec::SmallVec;
 
         let opts = self.convergence;
+
+        if self.time_slices.is_empty() {
+            return Ok(ConvergenceReport {
+                iterations: 0,
+                final_step: (0.0, 0.0),
+                log_evidence: 0.0,
+                converged: true,
+                per_iteration_time: SmallVec::new(),
+                slices_skipped: 0,
+            });
+        }
+
         let mut step = (f64::INFINITY, f64::INFINITY);
         let mut i = 0;
         let mut per_iter: SmallVec<[std::time::Duration; 32]> = SmallVec::new();
@@ -444,8 +509,24 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             per_iter.push(t0.elapsed());
             i += 1;
             self.observer.on_iteration_end(i, step);
+
+            // A non-finite step means EP has broken down; further iterations
+            // cannot recover, and `tuple_gt` would read NaN as converged.
+            if !crate::step_is_finite(step) {
+                break;
+            }
         }
-        let converged = !tuple_gt(step, opts.epsilon);
+
+        if !crate::step_is_finite(step) {
+            self.observer.on_converged(i, step, false);
+
+            return Err(InferenceError::NonFiniteResult {
+                context: "History::converge",
+                step,
+            });
+        }
+
+        let converged = crate::step_converged(step, opts.epsilon);
         let log_evidence = self.log_evidence_internal(false, &[]);
         self.observer.on_converged(i, step, converged);
         Ok(ConvergenceReport {
@@ -498,6 +579,21 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
             });
         }
 
+        // Chokepoint for tie validation: every ingestion route lands here,
+        // including `record_draw`, which builds its results directly rather
+        // than going through `Outcome`.
+        if self.p_draw == 0.0 {
+            for (event_results, kind) in results.iter().zip(kinds.iter()) {
+                if !matches!(kind, EventKind::Ranked) {
+                    continue;
+                }
+
+                if let Some(teams) = crate::first_tied_output(event_results) {
+                    return Err(InferenceError::TieWithoutDrawProbability { teams });
+                }
+            }
+        }
+
         competitor::clean(self.agents.values_mut(), true);
 
         let mut this_agent = Vec::with_capacity(1024);
@@ -593,6 +689,8 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
                     agent.last_time = Some(t);
                     agent.message = time_slice.forward_prior_out(&agent_idx);
                 }
+
+                k += 1;
             } else {
                 let mut time_slice = TimeSlice::new(t, self.p_draw, self.convergence);
                 time_slice.add_events(composition, results, weights, kinds_chunk, &self.agents);
@@ -734,10 +832,13 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
                 }
                 crate::Outcome::Scored { scores, sigma } => {
                     let resolved = sigma.unwrap_or(self.score_sigma);
-                    debug_assert!(
-                        resolved > 0.0,
-                        "resolved score_sigma must be > 0.0 (got {resolved})"
-                    );
+                    if resolved <= 0.0 || resolved.is_nan() {
+                        return Err(InferenceError::InvalidParameter {
+                            name: "score_sigma",
+                            value: resolved,
+                        });
+                    }
+
                     kinds.push(EventKind::Scored {
                         score_sigma: resolved,
                     });

src/key_table.rs

@@ -12,59 +12,68 @@ use crate::Index;
 /// crate. Power users can promote `&K` to `Index` via `get_or_create` and
 /// skip the lookup on subsequent hot-path calls.
 #[derive(Debug)]
-pub struct KeyTable<K>(HashMap<K, Index>);
+pub struct KeyTable<K> {
+    forward: HashMap<K, Index>,
+    /// Reverse mapping, indexed by `Index.0`.
+    ///
+    /// Indices are handed out densely and sequentially, so position *is* the
+    /// index and `key()` is a lookup rather than a scan over every entry.
+    reverse: Vec<K>,
+}
 
 impl<K> KeyTable<K>
 where
-    K: Eq + Hash,
+    K: Eq + Hash + Clone,
 {
     pub fn new() -> Self {
-        Self(HashMap::new())
+        Self {
+            forward: HashMap::new(),
+            reverse: Vec::new(),
+        }
     }
 
     pub fn get<Q: ?Sized + Hash + Eq>(&self, k: &Q) -> Option<Index>
     where
         K: Borrow<Q>,
     {
-        self.0.get(k).cloned()
+        self.forward.get(k).cloned()
     }
 
     pub fn get_or_create<Q: ?Sized + Hash + Eq + ToOwned<Owned = K>>(&mut self, k: &Q) -> Index
     where
         K: Borrow<Q>,
     {
-        if let Some(idx) = self.0.get(k) {
+        if let Some(idx) = self.forward.get(k) {
             *idx
         } else {
-            let idx = Index::from(self.0.len());
-            self.0.insert(k.to_owned(), idx);
+            let idx = Index::from(self.reverse.len());
+            let owned = k.to_owned();
+            self.reverse.push(owned.clone());
+            self.forward.insert(owned, idx);
             idx
         }
     }
 
     pub fn key(&self, idx: Index) -> Option<&K> {
-        self.0
-            .iter()
-            .find(|&(_, value)| *value == idx)
-            .map(|(key, _)| key)
+        self.reverse.get(idx.0)
     }
 
     pub fn keys(&self) -> impl Iterator<Item = &K> {
-        self.0.keys()
+        self.forward.keys()
     }
 
     pub fn len(&self) -> usize {
-        self.0.len()
+        self.reverse.len()
     }
 
     pub fn is_empty(&self) -> bool {
-        self.0.is_empty()
+        self.reverse.is_empty()
     }
 }
 
 impl<K> Default for KeyTable<K>
 where
-    K: Eq + Hash,
+    K: Eq + Hash + Clone,
 {
     fn default() -> Self {
         KeyTable::new()

src/lib.rs

@@ -1,3 +1,91 @@
+//! TrueSkill Through Time — Bayesian skill rating over a time axis.
+//!
+//! Where plain TrueSkill gives each competitor one running estimate, TrueSkill
+//! Through Time treats a whole history as a single model and infers skill *at
+//! every point in time*. Evidence flows both directions: a result today
+//! sharpens the estimate of who someone was last year, so early estimates stop
+//! being frozen guesses and comparisons across eras become meaningful.
+//!
+//! This is a Rust port of
+//! [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py).
+//!
+//! # Getting started
+//!
+//! Record results, converge, then read off skills:
+//!
+//! ```
+//! use trueskill_tt::History;
+//!
+//! let mut history = History::default();
+//!
+//! history.record_winner(&"alice", &"bob", 1)?;
+//! history.record_winner(&"bob", &"carol", 2)?;
+//! history.record_winner(&"alice", &"carol", 3)?;
+//!
+//! let report = history.converge()?;
+//! assert!(report.converged);
+//!
+//! let alice = history.current_skill("alice").unwrap();
+//! assert!(alice.mu() > 0.0, "alice won every game she played");
+//! # Ok::<(), trueskill_tt::InferenceError>(())
+//! ```
+//!
+//! Teams, weights, explicit rankings and continuous scores go through the
+//! fluent event builder:
+//!
+//! ```
+//! use trueskill_tt::History;
+//!
+//! let mut history = History::builder().p_draw(0.1).build();
+//!
+//! history
+//!     .event(1)
+//!     .team(["alice", "bob"])
+//!     .team(["carol", "dave"])
+//!     .ranking([0, 1])
+//!     .commit()?;
+//!
+//! history.converge()?;
+//! # Ok::<(), trueskill_tt::InferenceError>(())
+//! ```
+//!
+//! # Draws need a draw probability
+//!
+//! A `p_draw` of zero asserts that draws cannot happen, so a tied result has
+//! no representable likelihood and is rejected:
+//!
+//! ```
+//! use trueskill_tt::{History, InferenceError};
+//!
+//! let mut history = History::default(); // p_draw defaults to 0.0
+//! let err = history.record_draw(&"alice", &"bob", 1).unwrap_err();
+//! assert!(matches!(err, InferenceError::TieWithoutDrawProbability { .. }));
+//! ```
+//!
+//! This also applies to [`Outcome::winner`] for three or more teams, which
+//! ties every loser. Configure a positive `p_draw` for those.
+//!
+//! # Core types
+//!
+//! - [`History`] — the top-level container: ingests events, runs
+//!   forward/backward message passing, and answers queries.
+//! - [`Gaussian`] — the probability type, stored in natural parameters
+//!   (`pi = 1/sigma²`, `tau = mu/sigma²`) so message passing is add/subtract.
+//! - [`Game`] — one match in isolation, for scoring a hypothetical without a
+//!   history.
+//! - [`Outcome`] — how a match ended: ranks, or continuous scores.
+//! - [`Rating`] — a competitor's static configuration (prior, `beta`, drift).
+//!
+//! # Feature flags
+//!
+//! - `approx` — implements [`approx`](https://docs.rs/approx) equality traits
+//!   for [`Gaussian`]. Useful in tests.
+//! - `rayon` — parallelises the within-slice sweep and the per-slice passes of
+//!   `learning_curves`/`log_evidence`. Opt-in; results stay bit-identical
+//!   regardless of worker count.
+
+#![forbid(unsafe_code)]
+
 use std::{
     cmp::Reverse,
     f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2},
@@ -37,7 +125,7 @@ pub use event::{Event, Member, Team};
 pub use event_builder::EventBuilder;
 pub use game::{Game, GameOptions, OwnedGame};
 pub use gaussian::Gaussian;
-pub use history::History;
+pub use history::{History, HistoryBuilder};
 pub use key_table::KeyTable;
 use matrix::Matrix;
 pub use observer::{NullObserver, Observer};
@@ -63,12 +151,29 @@ pub const N_INF: Gaussian = Gaussian::from_ms(0.0, f64::INFINITY);
 #[derive(Copy, Clone, Default, PartialEq, PartialOrd, Eq, Ord, Hash, Debug)]
 pub struct Index(usize);
 
+impl Index {
+    /// The underlying slot number.
+    ///
+    /// Indices are dense and assigned in interning order, so this is usable as
+    /// a key into a caller-side side table.
+    #[must_use]
+    pub fn get(self) -> usize {
+        self.0
+    }
+}
+
 impl From<usize> for Index {
     fn from(ix: usize) -> Self {
         Self(ix)
     }
 }
 
+impl From<Index> for usize {
+    fn from(idx: Index) -> Self {
+        idx.0
+    }
+}
+
 fn erfc(x: f64) -> f64 {
     let z = x.abs();
     let t = 1.0 / (1.0 + z / 2.0);
@@ -184,6 +289,56 @@ pub(crate) fn tuple_gt(t: (f64, f64), e: f64) -> bool {
     t.0 > e || t.1 > e
 }
 
+/// Whether a convergence step is finite in both components.
+///
+/// A NaN step means EP broke down numerically. Because every comparison
+/// against NaN is false, `tuple_gt` reads NaN as "below epsilon" — so
+/// convergence checks must test finiteness explicitly rather than inferring
+/// success from `!tuple_gt(..)`.
+pub(crate) fn step_is_finite(t: (f64, f64)) -> bool {
+    t.0.is_finite() && t.1.is_finite()
+}
+
+/// Whether a step counts as converged: finite *and* within `epsilon`.
+pub(crate) fn step_converged(t: (f64, f64), epsilon: f64) -> bool {
+    step_is_finite(t) && !tuple_gt(t, epsilon)
+}
+
+/// Indices of the first pair of teams sharing a rank, if any.
+///
+/// A tie is only representable when the draw probability is positive: with
+/// `p_draw == 0.0` the truncation margin collapses to zero and the two-sided
+/// tie update evaluates `0/0`. Callers use this to reject such events before
+/// they reach inference.
+pub(crate) fn first_tied_pair(ranks: &[u32]) -> Option<(usize, usize)> {
+    for (i, a) in ranks.iter().enumerate() {
+        for (j, b) in ranks.iter().enumerate().skip(i + 1) {
+            if a == b {
+                return Some((i, j));
+            }
+        }
+    }
+
+    None
+}
+
+/// As `first_tied_pair`, but over the engine's internal `f64` outputs.
+///
+/// Ranks reach the engine already converted to descending `f64` outputs, and
+/// `Game` decides a tie by exact equality of those values — so this mirrors
+/// the comparison inference itself performs.
+pub(crate) fn first_tied_output(outputs: &[f64]) -> Option<(usize, usize)> {
+    for (i, a) in outputs.iter().enumerate() {
+        for (j, b) in outputs.iter().enumerate().skip(i + 1) {
+            if a == b {
+                return Some((i, j));
+            }
+        }
+    }
+
+    None
+}
+
 pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
     let mut x: Vec<(usize, T)> = xs.iter().enumerate().map(|(i, &t)| (i, t)).collect();
 
@@ -197,7 +352,26 @@ pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
 }
 
 /// Calculates the match quality of the given rating groups. A result is the draw probability in the association
+///
+/// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a
+/// perfectly balanced match.
+///
+/// # Panics
+///
+/// Panics if fewer than two rating groups are supplied, or if any group is
+/// empty — match quality is a property of a contest between at least two
+/// non-empty sides.
 pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
+    assert!(
+        rating_groups.len() >= 2,
+        "quality() requires at least 2 rating groups, got {}",
+        rating_groups.len()
+    );
+    assert!(
+        rating_groups.iter().all(|group| !group.is_empty()),
+        "quality() requires every rating group to be non-empty"
+    );
+
     let flatten_ratings = rating_groups
         .iter()
         .flat_map(|group| group.iter())
@@ -221,8 +395,10 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
 
     let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length);
 
+    // Row `row` contrasts group `row` (+weight) against group `row + 1`
+    // (-weight). `t` is the column where the current group's players start;
+    // the negative block begins immediately after it.
     let mut t = 0;
-    let mut x = 0;
 
     for (row, group) in rating_groups.windows(2).enumerate() {
         let current = group[0];
@@ -230,17 +406,13 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
 
         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() {
+        for n in t..t + next.len() {
             rotated_a_matrix[(row, n)] = -flatten_weights[n];
         }
-
-        x += next.len();
     }
 
     let a_matrix = rotated_a_matrix.transpose();

src/matrix.rs

@@ -1,29 +1,13 @@
+//! Minimal dense matrix used by `quality()`.
+//!
+//! `determinant` and `inverse` go through one LU decomposition with partial
+//! pivoting — O(n³) and numerically stable. The previous implementation
+//! expanded cofactors recursively (O(n!), allocating a `Vec` per minor) and
+//! only implemented `inverse` for the 1×1 case, which limited `quality()` to
+//! exactly two rating groups.
+
 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]>,
@@ -31,6 +15,107 @@ pub struct Matrix {
     width: usize,
 }
 
+/// LU decomposition with partial pivoting: `PA = LU`, stored compactly.
+///
+/// `lu` holds `L` below the diagonal (unit diagonal implied) and `U` on and
+/// above it. `sign` is the determinant sign contributed by row swaps, or 0.0
+/// when the matrix is singular.
+struct Lu {
+    lu: Vec<f64>,
+    perm: Vec<usize>,
+    n: usize,
+    sign: f64,
+}
+
+impl Lu {
+    fn decompose(m: &Matrix) -> Self {
+        debug_assert_eq!(m.width, m.height, "LU requires a square matrix");
+
+        let n = m.width;
+        let mut lu = m.data.to_vec();
+        let mut perm: Vec<usize> = (0..n).collect();
+        let mut sign = 1.0;
+
+        for col in 0..n {
+            // Partial pivot: take the largest-magnitude candidate to limit
+            // growth of round-off in the elimination below.
+            let mut pivot_row = col;
+            let mut pivot_max = lu[col * n + col].abs();
+
+            for row in (col + 1)..n {
+                let candidate = lu[row * n + col].abs();
+                if candidate > pivot_max {
+                    pivot_max = candidate;
+                    pivot_row = row;
+                }
+            }
+
+            if pivot_max == 0.0 {
+                sign = 0.0;
+                continue;
+            }
+
+            if pivot_row != col {
+                for k in 0..n {
+                    lu.swap(col * n + k, pivot_row * n + k);
+                }
+                perm.swap(col, pivot_row);
+                sign = -sign;
+            }
+
+            let pivot = lu[col * n + col];
+
+            for row in (col + 1)..n {
+                let factor = lu[row * n + col] / pivot;
+                lu[row * n + col] = factor;
+
+                for k in (col + 1)..n {
+                    lu[row * n + k] -= factor * lu[col * n + k];
+                }
+            }
+        }
+
+        Self { lu, perm, n, sign }
+    }
+
+    fn determinant(&self) -> f64 {
+        if self.sign == 0.0 {
+            return 0.0;
+        }
+
+        let mut det = self.sign;
+        for i in 0..self.n {
+            det *= self.lu[i * self.n + i];
+        }
+
+        det
+    }
+
+    /// Solve `Ax = b` for a single column of the identity, giving one column
+    /// of the inverse.
+    fn solve_column(&self, col: usize, out: &mut [f64]) {
+        let n = self.n;
+
+        // Forward substitution through L, applying the row permutation.
+        for i in 0..n {
+            let mut sum = if self.perm[i] == col { 1.0 } else { 0.0 };
+            for (k, &solved) in out.iter().enumerate().take(i) {
+                sum -= self.lu[i * n + k] * solved;
+            }
+            out[i] = sum;
+        }
+
+        // Back substitution through U.
+        for i in (0..n).rev() {
+            let mut sum = out[i];
+            for (k, &solved) in out.iter().enumerate().skip(i + 1) {
+                sum -= self.lu[i * n + k] * solved;
+            }
+            out[i] = sum / self.lu[i * n + i];
+        }
+    }
+}
+
 impl Matrix {
     pub fn new(height: usize, width: usize) -> Matrix {
         Matrix {
@@ -52,73 +137,59 @@ impl Matrix {
         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
-    }
-
+    /// Determinant of a square matrix. The 0×0 determinant is 1 by convention
+    /// (the empty product).
+    ///
+    /// # Panics
+    ///
+    /// Panics if the matrix is not square.
     pub fn determinant(&self) -> f64 {
-        debug_assert!(self.width == self.height);
+        assert_eq!(
+            self.width, self.height,
+            "determinant requires a square matrix, got {}x{}",
+            self.height, self.width
+        );
 
-        det(&self.data, self.width)
+        if self.width == 0 {
+            return 1.0;
+        }
+
+        Lu::decompose(self).determinant()
     }
 
-    pub fn adjugate(&self) -> Matrix {
-        debug_assert!(self.width == self.height);
+    /// Matrix inverse via LU decomposition.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the matrix is not square or is singular.
+    pub fn inverse(&self) -> Matrix {
+        assert_eq!(
+            self.width, self.height,
+            "inverse requires a square matrix, got {}x{}",
+            self.height, self.width
+        );
 
-        let mut matrix = Matrix::new(self.height, self.width);
+        let n = self.width;
+        let mut inverse = Matrix::new(n, n);
 
-        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 };
+        if n == 0 {
+            return inverse;
+        }
 
-                    matrix[(r, c)] = self.minor(r, c).determinant() * sign;
-                }
+        let lu = Lu::decompose(self);
+        assert!(lu.sign != 0.0, "cannot invert a singular matrix");
+
+        let mut column = vec![0.0; n];
+
+        for c in 0..n {
+            lu.solve_column(c, &mut column);
+
+            for (r, &value) in column.iter().enumerate() {
+                inverse[(r, c)] = value;
             }
         }
 
-        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
+        inverse
     }
 }
 
@@ -126,20 +197,62 @@ impl ops::Index<(usize, usize)> for Matrix {
     type Output = f64;
 
     fn index(&self, pos: (usize, usize)) -> &Self::Output {
+        debug_assert!(
+            pos.0 < self.height && pos.1 < self.width,
+            "index ({}, {}) out of bounds for {}x{} matrix",
+            pos.0,
+            pos.1,
+            self.height,
+            self.width
+        );
+
         &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 {
+        debug_assert!(
+            pos.0 < self.height && pos.1 < self.width,
+            "index ({}, {}) out of bounds for {}x{} matrix",
+            pos.0,
+            pos.1,
+            self.height,
+            self.width
+        );
+
         &mut self.data[(self.width * pos.0) + pos.1]
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for f64 {
+fn multiply(lhs: &Matrix, rhs: &Matrix) -> Matrix {
+    assert_eq!(
+        lhs.width, rhs.height,
+        "cannot multiply {}x{} by {}x{}",
+        lhs.height, lhs.width, rhs.height, rhs.width
+    );
+
+    let mut matrix = Matrix::new(lhs.height, rhs.width);
+
+    for r in 0..matrix.height {
+        for c in 0..matrix.width {
+            let mut value = 0.0;
+
+            for x in 0..lhs.width {
+                value += lhs[(r, x)] * rhs[(x, c)];
+            }
+
+            matrix[(r, c)] = value;
+        }
+    }
+
+    matrix
+}
+
+impl ops::Mul<&Matrix> for f64 {
     type Output = Matrix;
 
-    fn mul(self, rhs: &'a Matrix) -> Matrix {
+    fn mul(self, rhs: &Matrix) -> Matrix {
         let mut matrix = Matrix::new(rhs.height, rhs.width);
 
         for r in 0..rhs.height {
@@ -152,54 +265,35 @@ impl<'a> ops::Mul<&'a Matrix> for f64 {
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for Matrix {
+impl ops::Mul<&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
+    fn mul(self, rhs: &Matrix) -> Matrix {
+        multiply(&self, rhs)
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for &'a Matrix {
+impl ops::Mul<&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
+    fn mul(self, rhs: &Matrix) -> Matrix {
+        multiply(self, rhs)
     }
 }
 
-impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
+impl ops::Add<&Matrix> for &Matrix {
     type Output = Matrix;
 
-    fn add(self, rhs: &'a Matrix) -> Matrix {
+    fn add(self, rhs: &Matrix) -> Matrix {
+        assert!(
+            self.height == rhs.height && self.width == rhs.width,
+            "cannot add {}x{} to {}x{}",
+            self.height,
+            self.width,
+            rhs.height,
+            rhs.width
+        );
+
         let mut matrix = Matrix::new(self.height, self.width);
 
         for r in 0..matrix.height {
@@ -211,3 +305,106 @@ impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
         matrix
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn from_rows(rows: &[&[f64]]) -> Matrix {
+        let mut m = Matrix::new(rows.len(), rows[0].len());
+        for (r, row) in rows.iter().enumerate() {
+            for (c, &v) in row.iter().enumerate() {
+                m[(r, c)] = v;
+            }
+        }
+        m
+    }
+
+    #[test]
+    fn determinant_1x1() {
+        assert!((from_rows(&[&[3.0]]).determinant() - 3.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_2x2() {
+        let m = from_rows(&[&[1.0, 2.0], &[3.0, 4.0]]);
+        assert!((m.determinant() - (-2.0)).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_3x3() {
+        let m = from_rows(&[&[6.0, 1.0, 1.0], &[4.0, -2.0, 5.0], &[2.0, 8.0, 7.0]]);
+        assert!((m.determinant() - (-306.0)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn determinant_requires_no_pivot_at_origin() {
+        // A zero in the top-left forces a row swap; the sign must follow.
+        let m = from_rows(&[&[0.0, 1.0], &[1.0, 0.0]]);
+        assert!((m.determinant() - (-1.0)).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_of_singular_is_zero() {
+        let m = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]);
+        assert!(m.determinant().abs() < 1e-12);
+    }
+
+    #[test]
+    fn inverse_1x1() {
+        let inv = from_rows(&[&[4.0]]).inverse();
+        assert!((inv[(0, 0)] - 0.25).abs() < 1e-12);
+    }
+
+    #[test]
+    fn inverse_times_original_is_identity() {
+        for rows in [
+            vec![vec![1.0, 2.0], vec![3.0, 4.0]],
+            vec![
+                vec![6.0, 1.0, 1.0],
+                vec![4.0, -2.0, 5.0],
+                vec![2.0, 8.0, 7.0],
+            ],
+            vec![
+                vec![2.0, 0.0, 1.0, 3.0],
+                vec![1.0, 5.0, 2.0, 0.0],
+                vec![0.0, 1.0, 4.0, 1.0],
+                vec![3.0, 2.0, 0.0, 6.0],
+            ],
+        ] {
+            let refs: Vec<&[f64]> = rows.iter().map(|r| r.as_slice()).collect();
+            let m = from_rows(&refs);
+            let product = &m * &m.inverse();
+
+            for r in 0..product.height {
+                for c in 0..product.width {
+                    let expected = if r == c { 1.0 } else { 0.0 };
+                    assert!(
+                        (product[(r, c)] - expected).abs() < 1e-9,
+                        "({r},{c}) = {} expected {expected}",
+                        product[(r, c)]
+                    );
+                }
+            }
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "singular")]
+    fn inverse_of_singular_panics() {
+        let _ = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]).inverse();
+    }
+
+    #[test]
+    fn empty_determinant_is_one() {
+        assert!((Matrix::new(0, 0).determinant() - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn transpose_round_trips() {
+        let m = from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
+        let t = m.transpose();
+        assert_eq!((t.height, t.width), (3, 2));
+        assert_eq!(t.transpose()[(1, 2)], m[(1, 2)]);
+    }
+}

src/outcome.rs

@@ -57,9 +57,11 @@ impl Outcome {
 
     /// Explicit per-team continuous scores with a per-event noise override.
     ///
-    /// `sigma` must be `> 0.0`; debug-asserts otherwise.
+    /// `sigma` must be `> 0.0`. Constructing an `Outcome` with a non-positive
+    /// or NaN sigma is allowed; the value is rejected with
+    /// `InferenceError::InvalidParameter` when the event is ingested, so
+    /// callers get an error rather than a panic.
     pub fn scores_with_sigma<I: IntoIterator<Item = f64>>(scores: I, sigma: f64) -> Self {
-        debug_assert!(sigma > 0.0, "score_sigma must be > 0.0 (got {sigma})");
         Self::Scored {
             scores: scores.into_iter().collect(),
             sigma: Some(sigma),
@@ -169,9 +171,15 @@ mod tests {
         }
     }
 
+    /// Construction accepts any sigma; the value is validated at ingestion so
+    /// callers receive an `InferenceError` rather than a panic. See
+    /// `tests/degenerate_inputs.rs::scored_event_rejects_non_positive_sigma`.
     #[test]
-    #[should_panic(expected = "score_sigma must be > 0.0")]
-    fn scores_with_sigma_rejects_zero() {
-        let _ = Outcome::scores_with_sigma([3.0, 1.0], 0.0);
+    fn scores_with_sigma_defers_validation_to_ingestion() {
+        let o = Outcome::scores_with_sigma([3.0, 1.0], 0.0);
+        match o {
+            Outcome::Scored { sigma, .. } => assert_eq!(sigma, Some(0.0)),
+            Outcome::Ranked(_) => panic!("expected Scored variant"),
+        }
     }
 }

src/rating.rs

@@ -29,6 +29,24 @@ impl<T: Time, D: Drift<T>> Rating<T, D> {
         }
     }
 
+    /// The configured prior skill estimate.
+    #[must_use]
+    pub fn prior(&self) -> Gaussian {
+        self.prior
+    }
+
+    /// Performance noise: how much a single showing varies around the skill.
+    #[must_use]
+    pub fn beta(&self) -> f64 {
+        self.beta
+    }
+
+    /// The drift model governing how skill may move between events.
+    #[must_use]
+    pub fn drift(&self) -> D {
+        self.drift
+    }
+
     pub(crate) fn performance(&self) -> Gaussian {
         self.prior.forget(self.beta.powi(2))
     }

src/schedule.rs

@@ -32,8 +32,17 @@ pub struct EpsilonOrMax {
 
 impl Default for EpsilonOrMax {
     fn default() -> Self {
-        // Matches today's hard-coded tolerance and iteration cap.
-        Self { eps: 1e-6, max: 10 }
+        // Derived from `ConvergenceOptions` so there is one source of truth for
+        // the tolerance and iteration cap. These previously disagreed: this
+        // default capped at 10 iterations while `ConvergenceOptions` allowed 30,
+        // and which applied depended on whether inference went through
+        // `run_chain` or a `Schedule`.
+        let defaults = crate::ConvergenceOptions::default();
+
+        Self {
+            eps: defaults.epsilon,
+            max: defaults.max_iter,
+        }
     }
 }
 
@@ -50,10 +59,16 @@ impl Schedule for EpsilonOrMax {
         }
 
         let mut iterations = 0;
-        let mut final_step = (f64::INFINITY, f64::INFINITY);
-        let mut converged = false;
+        // With no iterating factors the graph is already at its fixed point:
+        // the setup pass above is all there is to do. Reporting `converged:
+        // false` with an infinite step for that case gave callers a false
+        // negative.
+        let mut final_step = (0.0, 0.0);
+        let mut converged = true;
 
         if n_setup < factors.len() {
+            final_step = (f64::INFINITY, f64::INFINITY);
+            converged = false;
             for _ in 0..self.max {
                 let mut step = (0.0_f64, 0.0_f64);
 
@@ -113,7 +128,8 @@ mod tests {
 
     #[test]
     fn report_marks_converged_when_no_iterating_factors() {
-        // No iterating factors → 0 iterations, converged stays false (loop never ran).
+        // A graph of only setup factors has nothing to iterate, so it is at its
+        // fixed point after the setup pass: 0 iterations, and converged.
         let mut vars = VarStore::new();
         let out = vars.alloc(N_INF);
         let mut factors = vec![BuiltinFactor::TeamSum(TeamSumFactor {
@@ -122,5 +138,15 @@ mod tests {
         })];
         let report = EpsilonOrMax::default().run(&mut factors, &mut vars);
         assert_eq!(report.iterations, 0);
+        assert!(report.converged);
+        assert_eq!(report.final_step, (0.0, 0.0));
+    }
+
+    #[test]
+    fn default_matches_convergence_options() {
+        let schedule = EpsilonOrMax::default();
+        let options = crate::ConvergenceOptions::default();
+        assert_eq!(schedule.max, options.max_iter);
+        assert_eq!(schedule.eps, options.epsilon);
     }
 }

src/storage/skill_store.rs

@@ -41,6 +41,18 @@ impl SkillStore {
         }
     }
 
+    /// Whether a slot is occupied. Test-only.
+    #[cfg(test)]
+    pub fn contains(&self, idx: Index) -> bool {
+        idx.0 < self.present.len() && self.present[idx.0]
+    }
+
+    /// Number of occupied slots. Test-only.
+    #[cfg(test)]
+    pub fn len(&self) -> usize {
+        self.n_present
+    }
+
     pub fn get_mut(&mut self, idx: Index) -> Option<&mut Skill> {
         if idx.0 < self.present.len() && self.present[idx.0] {
             Some(&mut self.skills[idx.0])
@@ -49,21 +61,6 @@ impl SkillStore {
         }
     }
 
-    #[allow(dead_code)]
-    pub fn contains(&self, idx: Index) -> bool {
-        idx.0 < self.present.len() && self.present[idx.0]
-    }
-
-    #[allow(dead_code)]
-    pub fn len(&self) -> usize {
-        self.n_present
-    }
-
-    #[allow(dead_code)]
-    pub fn is_empty(&self) -> bool {
-        self.n_present == 0
-    }
-
     pub fn iter(&self) -> impl Iterator<Item = (Index, &Skill)> {
         self.present.iter().enumerate().filter_map(|(i, &p)| {
             if p {

src/time_slice.rs

@@ -14,7 +14,6 @@ use crate::{
     rating::Rating,
     storage::{CompetitorStore, SkillStore},
     time::Time,
-    tuple_gt, tuple_max,
 };
 
 #[derive(Debug)]
@@ -87,7 +86,7 @@ struct Team {
 #[derive(Debug)]
 pub(crate) struct Event {
     teams: Vec<Team>,
-    evidence: f64,
+    log_evidence: f64,
     weights: Vec<Vec<f64>>,
     kind: EventKind,
 }
@@ -124,18 +123,19 @@ impl Event {
             .collect::<Vec<_>>()
     }
 
-    /// Direct in-loop update: mutates self and `skills` inline with no
-    /// intermediate allocation. Used by both the sequential sweep path and,
-    /// via unsafe, by the parallel rayon path for events in the same color
-    /// group (which have disjoint agent sets — see `sweep_color_groups`).
-    fn iteration_direct<T: Time, D: Drift<T>>(
-        &mut self,
-        skills: &mut SkillStore,
+    /// Run inference for this event and return its per-item likelihoods.
+    ///
+    /// Reads `skills` immutably and does not touch `self`, so every event in
+    /// a color group can run concurrently without any aliasing question —
+    /// the mutation is deferred to `apply`.
+    fn compute<T: Time, D: Drift<T>>(
+        &self,
+        skills: &SkillStore,
         agents: &CompetitorStore<T, D>,
         p_draw: f64,
         convergence: crate::ConvergenceOptions,
         arena: &mut ScratchArena,
-    ) {
+    ) -> EventUpdate {
         let teams = self.within_priors(false, false, skills, agents);
         let result = self.outputs();
         let g = match self.kind {
@@ -152,17 +152,47 @@ impl Event {
             ),
         };
 
+        EventUpdate {
+            log_evidence: g.log_evidence,
+            likelihoods: g.likelihoods,
+        }
+    }
+
+    /// Fold a computed update into the skill store and cache it on the items.
+    fn apply(&mut self, skills: &mut SkillStore, update: EventUpdate) {
         for (t, team) in self.teams.iter_mut().enumerate() {
             for (i, item) in team.items.iter_mut().enumerate() {
+                let fresh = update.likelihoods[t][i];
                 let old_likelihood = skills.get(item.agent).unwrap().likelihood;
-                let new_likelihood = (old_likelihood / item.likelihood) * g.likelihoods[t][i];
+                let new_likelihood = (old_likelihood / item.likelihood) * fresh;
                 skills.get_mut(item.agent).unwrap().likelihood = new_likelihood;
-                item.likelihood = g.likelihoods[t][i];
+                item.likelihood = fresh;
             }
         }
 
-        self.evidence = g.evidence;
+        self.log_evidence = update.log_evidence;
     }
+
+    /// Compute and apply in one step — the sequential sweep.
+    fn iteration_direct<T: Time, D: Drift<T>>(
+        &mut self,
+        skills: &mut SkillStore,
+        agents: &CompetitorStore<T, D>,
+        p_draw: f64,
+        convergence: crate::ConvergenceOptions,
+        arena: &mut ScratchArena,
+    ) {
+        let update = self.compute(skills, agents, p_draw, convergence, arena);
+        self.apply(skills, update);
+    }
+}
+
+/// The result of running inference for one event, before it is folded back
+/// into the shared skill store.
+#[derive(Debug)]
+struct EventUpdate {
+    log_evidence: f64,
+    likelihoods: Vec<Vec<Gaussian>>,
 }
 
 #[derive(Debug)]
@@ -174,6 +204,14 @@ pub struct TimeSlice<T: Time = i64> {
     pub(crate) convergence: crate::ConvergenceOptions,
     arena: ScratchArena,
     pub(crate) color_groups: ColorGroups,
+    /// Whether `color_groups` still reflects `events`.
+    ///
+    /// Coloring is rebuilt lazily, on the first full sweep after an append,
+    /// rather than eagerly per append: the partition is thrown away and
+    /// recomputed wholesale either way, so doing it per append made ingesting
+    /// n events O(n^2) with no benefit — nothing reads the partition between
+    /// an append and the next full sweep.
+    color_groups_dirty: bool,
 }
 
 impl<T: Time> TimeSlice<T> {
@@ -186,6 +224,7 @@ impl<T: Time> TimeSlice<T> {
             convergence,
             arena: ScratchArena::new(),
             color_groups: ColorGroups::new(),
+            color_groups_dirty: false,
         }
     }
 
@@ -198,6 +237,7 @@ impl<T: Time> TimeSlice<T> {
         let n = self.events.len();
         if n == 0 {
             self.color_groups = ColorGroups::new();
+            self.color_groups_dirty = false;
             return;
         }
 
@@ -221,6 +261,12 @@ impl<T: Time> TimeSlice<T> {
 
         self.events = reordered;
         self.color_groups = ColorGroups { groups: new_groups };
+        self.color_groups_dirty = false;
+
+        debug_assert!(
+            self.color_groups.groups_are_contiguous(),
+            "color groups must occupy contiguous event ranges"
+        );
     }
 
     pub fn add_events<D: Drift<T>>(
@@ -296,7 +342,7 @@ impl<T: Time> TimeSlice<T> {
 
             Event {
                 teams,
-                evidence: 0.0,
+                log_evidence: 0.0,
                 weights,
                 kind: kinds[e],
             }
@@ -306,8 +352,9 @@ impl<T: Time> TimeSlice<T> {
 
         self.events.extend(events);
 
+        self.color_groups_dirty = true;
+
         self.iteration(from, agents);
-        self.recompute_color_groups();
     }
 
     pub(crate) fn posteriors(&self) -> HashMap<Index, Gaussian> {
@@ -318,6 +365,10 @@ impl<T: Time> TimeSlice<T> {
     }
 
     pub fn iteration<D: Drift<T>>(&mut self, from: usize, agents: &CompetitorStore<T, D>) {
+        if from == 0 && self.color_groups_dirty {
+            self.recompute_color_groups();
+        }
+
         if from > 0 || self.color_groups.is_empty() {
             // Initial pass (add_events) or no color groups yet: simple sequential sweep.
             for event in self.events.iter_mut().skip(from) {
@@ -353,7 +404,7 @@ impl<T: Time> TimeSlice<T> {
                     }
                 }
 
-                event.evidence = g.evidence;
+                event.log_evidence = g.log_evidence;
             }
         } else {
             self.sweep_color_groups(agents);
@@ -363,14 +414,13 @@ impl<T: Time> TimeSlice<T> {
     /// Full event sweep using the color-group partition. Colors are processed
     /// sequentially; within each color the inner loop is parallel under rayon.
     ///
-    /// Events within each color group touch disjoint agent sets (guaranteed by
-    /// the greedy coloring). This lets each rayon thread write directly to its
-    /// events' skill likelihoods without a deferred-apply step, matching the
-    /// sequential path's allocation profile. The unsafe block is sound because:
-    ///   1. `self.events[range]` and `self.skills` are separate fields → disjoint.
-    ///   2. Events in the same color group access disjoint `Index` values in
-    ///      `self.skills`, so concurrent writes land on different memory locations.
-    ///   3. Each event only writes to its own items' likelihoods (no sharing).
+    /// Events in one color group touch disjoint agent sets, so none of them
+    /// can observe another's writes. That makes the sweep separable: inference
+    /// runs concurrently over shared `&self.skills`, and the resulting updates
+    /// are folded in afterwards in index order. Splitting it this way needs no
+    /// `unsafe` and no aliasing argument, and it keeps results bit-identical
+    /// across thread counts because the apply order does not depend on which
+    /// worker finished first.
     #[cfg(feature = "rayon")]
     fn sweep_color_groups<D: Drift<T>>(&mut self, agents: &CompetitorStore<T, D>) {
         use rayon::prelude::*;
@@ -390,29 +440,28 @@ impl<T: Time> TimeSlice<T> {
             if group_len == 0 {
                 continue;
             }
+
             let range = self.color_groups.color_range(color_idx);
             let p_draw = self.p_draw;
             let convergence = self.convergence;
 
             if group_len >= RAYON_THRESHOLD {
-                // Obtain a raw pointer from the unique `&mut self.skills` reference.
-                // Casting back to `&mut` inside the closure is sound because:
-                //   1. The pointer originates from a `&mut` — no aliasing with shared refs.
-                //   2. Events in the same color group touch disjoint `Index` slots in the
-                //      underlying Vec, so concurrent writes from different threads land on
-                //      different memory locations — no data race.
-                //   3. `self.events[range]` and `self.skills` are separate struct fields,
-                //      so the borrow splits cleanly.
-                let skills_addr: usize = (&mut self.skills as *mut SkillStore) as usize;
-                self.events[range].par_iter_mut().for_each(move |ev| {
-                    // SAFETY: see above.
-                    let skills: &mut SkillStore = unsafe { &mut *(skills_addr as *mut SkillStore) };
-                    ARENA.with(|cell| {
-                        let mut arena = cell.borrow_mut();
-                        arena.reset();
-                        ev.iteration_direct(skills, agents, p_draw, convergence, &mut arena);
-                    });
-                });
+                let skills = &self.skills;
+                let updates: Vec<EventUpdate> = self.events[range.clone()]
+                    .par_iter()
+                    .map(|ev| {
+                        ARENA.with(|cell| {
+                            let mut arena = cell.borrow_mut();
+                            arena.reset();
+
+                            ev.compute(skills, agents, p_draw, convergence, &mut arena)
+                        })
+                    })
+                    .collect();
+
+                for (ev, update) in self.events[range].iter_mut().zip(updates) {
+                    ev.apply(&mut self.skills, update);
+                }
             } else {
                 for ev in &mut self.events[range] {
                     ev.iteration_direct(
@@ -454,18 +503,29 @@ impl<T: Time> TimeSlice<T> {
         }
     }
 
-    #[allow(dead_code)]
+    /// Iterate this slice alone until its posteriors stop moving, returning
+    /// the number of iterations taken.
+    ///
+    /// Only used by tests: production convergence is driven across slices by
+    /// `History::converge`.
+    ///
+    /// Honours `self.convergence`; it previously hard-coded an epsilon and a
+    /// 20-iteration cap that matched neither `ConvergenceOptions` nor the
+    /// schedule default.
+    #[cfg(test)]
     pub(crate) fn iterate_to_convergence<D: Drift<T>>(
         &mut self,
         agents: &CompetitorStore<T, D>,
     ) -> usize {
-        let epsilon = 1e-6;
-        let iterations = 20;
+        use crate::{tuple_gt, tuple_max};
+
+        let epsilon = self.convergence.epsilon;
+        let max_iter = self.convergence.max_iter;
 
         let mut step = (f64::INFINITY, f64::INFINITY);
         let mut i = 0;
 
-        while tuple_gt(step, epsilon) && i < iterations {
+        while tuple_gt(step, epsilon) && i < max_iter {
             let old = self.posteriors();
 
             self.iteration(0, agents);
@@ -477,6 +537,10 @@ impl<T: Time> TimeSlice<T> {
             });
 
             i += 1;
+
+            if !crate::step_is_finite(step) {
+                break;
+            }
         }
 
         i
@@ -523,6 +587,9 @@ impl<T: Time> TimeSlice<T> {
         forward: bool,
         agents: &CompetitorStore<T, D>,
     ) -> f64 {
+        // Hashed once rather than scanned per player per event, so a
+        // `log_evidence_for` with many keys is not quadratic.
+        let target_set: std::collections::HashSet<Index> = targets.iter().copied().collect();
         // log_evidence is infrequent; a local arena avoids needing &mut self.
         let mut arena = ScratchArena::new();
 
@@ -530,26 +597,28 @@ impl<T: Time> TimeSlice<T> {
             let teams = event.within_priors(online, forward, &self.skills, agents);
             let result = event.outputs();
             match event.kind {
-                EventKind::Ranked => Game::ranked_with_arena(
-                    teams,
-                    &result,
-                    &event.weights,
-                    self.p_draw,
-                    self.convergence,
-                    arena,
-                )
-                .evidence
-                .ln(),
-                EventKind::Scored { score_sigma } => Game::scored_with_arena(
-                    teams,
-                    &result,
-                    &event.weights,
-                    score_sigma,
-                    self.convergence,
-                    arena,
-                )
-                .evidence
-                .ln(),
+                EventKind::Ranked => {
+                    Game::ranked_with_arena(
+                        teams,
+                        &result,
+                        &event.weights,
+                        self.p_draw,
+                        self.convergence,
+                        arena,
+                    )
+                    .log_evidence
+                }
+                EventKind::Scored { score_sigma } => {
+                    Game::scored_with_arena(
+                        teams,
+                        &result,
+                        &event.weights,
+                        score_sigma,
+                        self.convergence,
+                        arena,
+                    )
+                    .log_evidence
+                }
             }
         };
 
@@ -560,7 +629,7 @@ impl<T: Time> TimeSlice<T> {
                     .map(|event| run_event(event, &mut arena))
                     .sum()
             } else {
-                self.events.iter().map(|event| event.evidence.ln()).sum()
+                self.events.iter().map(|event| event.log_evidence).sum()
             }
         } else if online || forward {
             self.events
@@ -570,7 +639,7 @@ impl<T: Time> TimeSlice<T> {
                         .teams
                         .iter()
                         .flat_map(|team| &team.items)
-                        .any(|item| targets.contains(&item.agent))
+                        .any(|item| target_set.contains(&item.agent))
                 })
                 .map(|event| run_event(event, &mut arena))
                 .sum()
@@ -582,9 +651,9 @@ impl<T: Time> TimeSlice<T> {
                         .teams
                         .iter()
                         .flat_map(|team| &team.items)
-                        .any(|item| targets.contains(&item.agent))
+                        .any(|item| target_set.contains(&item.agent))
                 })
-                .map(|event| event.evidence.ln())
+                .map(|event| event.log_evidence)
                 .sum()
         }
     }
@@ -866,19 +935,24 @@ mod tests {
 
         let post = time_slice.posteriors();
 
+        // These are convergence residuals, not exact values: by symmetry the
+        // true mean is 25.0 and the iteration approaches it from above. The
+        // previous expectation of 25.000003 was the residual after the
+        // hard-coded 20-iteration cap; honouring `ConvergenceOptions` runs to
+        // 30 and lands nearer the truth.
         assert_ulps_eq!(
             post[&a],
-            Gaussian::from_ms(25.000003, 3.880150),
+            Gaussian::from_ms(25.000001, 3.880150),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
             post[&b],
-            Gaussian::from_ms(25.000003, 3.880150),
+            Gaussian::from_ms(25.000001, 3.880150),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
             post[&c],
-            Gaussian::from_ms(25.000003, 3.880150),
+            Gaussian::from_ms(25.000001, 3.880150),
             epsilon = 1e-6
         );
     }

tests/degenerate_inputs.rs

@@ -0,0 +1,247 @@
+//! Degenerate, boundary, and error-path coverage.
+//!
+//! These run in both debug and release: the defects they pin were all
+//! guarded only by `debug_assert!`, so a debug-only suite never saw them.
+
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Game, GameOptions, Gaussian, History, InferenceError,
+    Outcome, Rating,
+};
+
+type R = Rating<i64, ConstantDrift>;
+
+fn rating() -> R {
+    R::new(
+        Gaussian::from_ms(25.0, 25.0 / 3.0),
+        25.0 / 6.0,
+        ConstantDrift(25.0 / 300.0),
+    )
+}
+
+fn assert_finite(g: Gaussian, what: &str) {
+    assert!(
+        g.mu().is_finite() && g.sigma().is_finite(),
+        "{what} must be finite, got mu={} sigma={}",
+        g.mu(),
+        g.sigma()
+    );
+}
+
+#[test]
+fn record_draw_without_draw_probability_is_rejected() {
+    let mut h = History::default();
+    let err = h.record_draw(&"a", &"b", 1).unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::TieWithoutDrawProbability { .. }
+    ));
+}
+
+#[test]
+fn builder_draw_without_draw_probability_is_rejected() {
+    let mut h = History::default();
+    let err = h
+        .event(1)
+        .team(["a"])
+        .team(["b"])
+        .draw()
+        .commit()
+        .unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::TieWithoutDrawProbability { .. }
+    ));
+}
+
+#[test]
+fn draw_with_positive_draw_probability_is_finite() {
+    let mut h = History::builder().p_draw(0.25).build();
+    h.record_draw(&"a", &"b", 1).unwrap();
+    let report = h.converge().unwrap();
+
+    assert_finite(h.current_skill("a").unwrap(), "drawn competitor skill");
+    assert_finite(h.current_skill("b").unwrap(), "drawn competitor skill");
+    assert!(report.log_evidence.is_finite());
+    assert!(report.converged);
+}
+
+#[test]
+fn game_ranked_rejects_tie_without_draw_probability() {
+    let a = [rating()];
+    let b = [rating()];
+    let teams: Vec<&[R]> = vec![&a, &b];
+    let err = Game::ranked(&teams, Outcome::draw(2), &GameOptions::default()).unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::TieWithoutDrawProbability { .. }
+    ));
+}
+
+/// `Outcome::winner(w, n)` ties every loser, so any n >= 3 free-for-all hits
+/// the tie path even though the caller never asked for a draw.
+#[test]
+fn winner_of_three_or_more_requires_draw_probability() {
+    let a = [rating()];
+    let b = [rating()];
+    let c = [rating()];
+    let teams: Vec<&[R]> = vec![&a, &b, &c];
+
+    let err = Game::ranked(&teams, Outcome::winner(0, 3), &GameOptions::default()).unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::TieWithoutDrawProbability { .. }
+    ));
+
+    let opts = GameOptions {
+        p_draw: 0.1,
+        ..GameOptions::default()
+    };
+    let game = Game::ranked(&teams, Outcome::winner(0, 3), &opts).unwrap();
+    for team in game.posteriors() {
+        for skill in team {
+            assert_finite(skill, "3-team winner posterior");
+        }
+    }
+}
+
+#[test]
+fn full_ranking_without_ties_needs_no_draw_probability() {
+    let a = [rating()];
+    let b = [rating()];
+    let c = [rating()];
+    let teams: Vec<&[R]> = vec![&a, &b, &c];
+    let game = Game::ranked(&teams, Outcome::ranking([0, 1, 2]), &GameOptions::default()).unwrap();
+
+    for team in game.posteriors() {
+        for skill in team {
+            assert_finite(skill, "strict ranking posterior");
+        }
+    }
+}
+
+#[test]
+fn empty_history_converges_trivially() {
+    let mut h = History::default();
+    let report = h.converge().unwrap();
+    assert_eq!(report.iterations, 0);
+    assert!(report.converged);
+}
+
+#[test]
+fn empty_event_stream_then_converge() {
+    let mut h = History::default();
+    h.add_events(std::iter::empty()).unwrap();
+    let report = h.converge().unwrap();
+    assert_eq!(report.iterations, 0);
+}
+
+#[test]
+fn empty_history_queries_do_not_panic() {
+    let h = History::default();
+    assert!(h.learning_curves().is_empty());
+    assert!(h.learning_curve("nobody").is_empty());
+    assert!(h.current_skill("nobody").is_none());
+}
+
+#[test]
+fn single_event_history_converges() {
+    let mut h = History::default();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let report = h.converge().unwrap();
+    assert!(report.converged);
+    assert_finite(h.current_skill("a").unwrap(), "single-event skill");
+}
+
+#[test]
+fn scored_event_rejects_non_positive_sigma() {
+    let mut h = History::builder().score_sigma(2.0).build();
+    let err = h
+        .event(1)
+        .team(["a"])
+        .team(["b"])
+        .scores_with_sigma([3.0, 1.0], f64::NAN)
+        .commit()
+        .unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::InvalidParameter {
+            name: "score_sigma",
+            ..
+        }
+    ));
+}
+
+#[test]
+fn convergence_reports_are_finite_across_many_teams() {
+    let opts = GameOptions {
+        p_draw: 0.1,
+        convergence: ConvergenceOptions::default(),
+        ..GameOptions::default()
+    };
+    let holders: Vec<[R; 1]> = (0..12).map(|_| [rating()]).collect();
+    let teams: Vec<&[R]> = holders.iter().map(|t| t.as_slice()).collect();
+    let game = Game::ranked(&teams, Outcome::ranking(0..12), &opts).unwrap();
+
+    assert!(
+        game.log_evidence().is_finite(),
+        "12-team log-evidence must be finite, got {}",
+        game.log_evidence()
+    );
+    for team in game.posteriors() {
+        for skill in team {
+            assert_finite(skill, "12-team posterior");
+        }
+    }
+}
+
+/// A long diff chain underflows a linear evidence product: each link
+/// contributes a probability in (0, 1], so ~1000 links flush the product to
+/// exactly 0.0 and `ln(0.0)` is `-inf`. Accumulating in log space keeps it
+/// finite.
+#[test]
+fn log_evidence_survives_a_long_diff_chain() {
+    let holders: Vec<[R; 1]> = (0..1200).map(|_| [rating()]).collect();
+    let teams: Vec<&[R]> = holders.iter().map(|t| t.as_slice()).collect();
+    let game = Game::ranked(
+        &teams,
+        Outcome::ranking(0..holders.len() as u32),
+        &GameOptions::default(),
+    )
+    .unwrap();
+
+    let log_evidence = game.log_evidence();
+    assert!(
+        log_evidence.is_finite(),
+        "1200-team log-evidence must be finite, got {log_evidence}"
+    );
+    assert!(
+        log_evidence < 0.0,
+        "log-evidence of a probability must be negative, got {log_evidence}"
+    );
+}
+
+/// A near-certain outcome rounds the losing tail to exactly zero in the
+/// `erfc` approximation; the evidence floor keeps `ln` finite.
+#[test]
+fn log_evidence_finite_for_near_certain_outcome() {
+    let overwhelming = R::new(Gaussian::from_ms(5_000.0, 0.5), 1.0, ConstantDrift(0.0));
+    let hopeless = R::new(Gaussian::from_ms(-5_000.0, 0.5), 1.0, ConstantDrift(0.0));
+    let a = [overwhelming];
+    let b = [hopeless];
+    let teams: Vec<&[R]> = vec![&a, &b];
+
+    let game = Game::ranked(&teams, Outcome::winner(0, 2), &GameOptions::default()).unwrap();
+    assert!(
+        game.log_evidence().is_finite(),
+        "got {}",
+        game.log_evidence()
+    );
+
+    // And the reverse — a colossal upset — must also stay finite.
+    let upset = Game::ranked(&teams, Outcome::winner(1, 2), &GameOptions::default()).unwrap();
+    assert!(
+        upset.log_evidence().is_finite(),
+        "upset log-evidence must be finite, got {}",
+        upset.log_evidence()
+    );
+}

tests/equivalence.rs

@@ -48,15 +48,9 @@ fn game_1v1_draw_golden() {
     )
     .unwrap();
     let p = g.posteriors();
-    // Historical golden from pre-T2 test_1vs1_draw:
-    assert_ulps_eq!(
-        p[0][0],
-        Gaussian::from_ms(24.999999, 6.469480),
-        epsilon = 1e-6
-    );
-    assert_ulps_eq!(
-        p[1][0],
-        Gaussian::from_ms(24.999999, 6.469480),
-        epsilon = 1e-6
-    );
+    // Historical golden from pre-T2 test_1vs1_draw. The mean is 25.0 exactly
+    // by symmetry — two identical competitors drawing cannot move apart — and
+    // the reference's 24.999999 is that value transcribed to six decimals.
+    assert_ulps_eq!(p[0][0], Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6);
+    assert_ulps_eq!(p[1][0], Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6);
 }

tests/ingestion_equivalence.rs

@@ -0,0 +1,147 @@
+//! Ingesting the same events must give the same answer however they were
+//! batched.
+//!
+//! The numerical goldens all ingest in a single call with one slice per
+//! timestamp, so they never exercise the "append to an existing slice" path.
+//! These do.
+
+use smallvec::smallvec;
+use trueskill_tt::{ConvergenceOptions, Event, Gaussian, History, Member, Outcome, Team};
+
+/// Converge tightly: the default cap of 30 iterations leaves a residual around
+/// 1e-6, which would swamp the comparison. Both paths must reach the same
+/// fixed point, so drive both well past it.
+fn tight() -> ConvergenceOptions {
+    ConvergenceOptions {
+        max_iter: 2_000,
+        epsilon: 1e-12,
+        ..ConvergenceOptions::default()
+    }
+}
+
+fn event(a: &str, b: &str, time: i64) -> Event<i64, String> {
+    Event {
+        time,
+        teams: smallvec![
+            Team::with_members([Member::new(a.to_string())]),
+            Team::with_members([Member::new(b.to_string())]),
+        ],
+        outcome: Outcome::winner(0, 2),
+    }
+}
+
+fn converged_skills(events: Vec<Event<i64, String>>, batched: bool) -> Vec<(String, Gaussian)> {
+    let mut h: History<i64, _, _, String> =
+        History::builder_with_key().convergence(tight()).build();
+
+    if batched {
+        h.add_events(events).unwrap();
+    } else {
+        for ev in events {
+            h.add_events(std::iter::once(ev)).unwrap();
+        }
+    }
+
+    let report = h.converge().unwrap();
+    assert!(
+        report.converged,
+        "fixture must converge before results can be compared; final step {:?}",
+        report.final_step
+    );
+
+    let mut skills: Vec<(String, Gaussian)> = h
+        .learning_curves()
+        .into_iter()
+        .map(|(key, curve)| (key, curve.last().unwrap().1))
+        .collect();
+    skills.sort_by(|a, b| a.0.cmp(&b.0));
+    skills
+}
+
+fn assert_same(batched: &[(String, Gaussian)], incremental: &[(String, Gaussian)], what: &str) {
+    assert_eq!(
+        batched.len(),
+        incremental.len(),
+        "{what}: competitor count differs"
+    );
+
+    for ((kb, gb), (ki, gi)) in batched.iter().zip(incremental.iter()) {
+        assert_eq!(kb, ki, "{what}: key order differs");
+        assert!(
+            (gb.mu() - gi.mu()).abs() < 1e-8 && (gb.sigma() - gi.sigma()).abs() < 1e-8,
+            "{what}: {kb} differs — batched mu={} sigma={}, incremental mu={} sigma={}",
+            gb.mu(),
+            gb.sigma(),
+            gi.mu(),
+            gi.sigma()
+        );
+    }
+}
+
+/// All events share one timestamp, so incremental ingestion repeatedly appends
+/// to an existing slice.
+#[test]
+fn same_slice_incremental_matches_batched() {
+    let events = vec![
+        event("a", "b", 1),
+        event("c", "d", 1),
+        event("e", "f", 1),
+        event("a", "c", 1),
+        event("b", "e", 1),
+    ];
+
+    let batched = converged_skills(events.clone(), true);
+    let incremental = converged_skills(events, false);
+    assert_same(&batched, &incremental, "single shared slice");
+}
+
+/// Distinct timestamps, so each append lands in a fresh slice appended after
+/// the existing ones.
+#[test]
+fn distinct_slices_incremental_matches_batched() {
+    let events = vec![
+        event("a", "b", 1),
+        event("b", "c", 2),
+        event("c", "a", 3),
+        event("a", "c", 4),
+    ];
+
+    let batched = converged_skills(events.clone(), true);
+    let incremental = converged_skills(events, false);
+    assert_same(&batched, &incremental, "distinct slices");
+}
+
+/// Several events per timestamp across several timestamps — appends to
+/// existing slices interleaved with new ones.
+#[test]
+fn mixed_slices_incremental_matches_batched() {
+    let events = vec![
+        event("a", "b", 1),
+        event("c", "d", 1),
+        event("a", "c", 2),
+        event("b", "d", 2),
+        event("a", "d", 3),
+        event("b", "c", 3),
+    ];
+
+    let batched = converged_skills(events.clone(), true);
+    let incremental = converged_skills(events, false);
+    assert_same(&batched, &incremental, "mixed slices");
+}
+
+/// Appending an event to a slice that is *not* the most recent one exercises
+/// the forward refresh of every later slice.
+#[test]
+fn back_dated_event_matches_batched() {
+    let events = vec![
+        event("a", "b", 1),
+        event("b", "c", 5),
+        event("c", "a", 9),
+        // arrives last, but belongs to the middle slice
+        event("a", "c", 5),
+    ];
+
+    let batched = converged_skills(events.clone(), true);
+    let incremental = converged_skills(events, false);
+    assert_same(&batched, &incremental, "back-dated event");
+}

tests/quality.rs

@@ -0,0 +1,119 @@
+//! `quality()` beyond two rating groups.
+//!
+//! The historical golden (two equal singletons) is asserted in
+//! `src/lib.rs::tests::test_quality`. These cover the N-group generalisation,
+//! which previously panicked with an out-of-bounds index at 3+ groups.
+
+use trueskill_tt::{Gaussian, quality};
+
+const BETA: f64 = 25.0 / 3.0 / 2.0;
+
+fn rating(mu: f64, sigma: f64) -> Gaussian {
+    Gaussian::from_ms(mu, sigma)
+}
+
+#[test]
+fn three_equal_groups_is_finite_and_in_range() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r], &[r], &[r]], BETA);
+
+    assert!(q.is_finite(), "quality must be finite, got {q}");
+    assert!((0.0..=1.0).contains(&q), "quality out of range: {q}");
+}
+
+#[test]
+fn quality_supports_many_groups() {
+    let r = rating(25.0, 3.0);
+    for n in 2..=8 {
+        let holders: Vec<[Gaussian; 1]> = (0..n).map(|_| [r]).collect();
+        let groups: Vec<&[Gaussian]> = holders.iter().map(|g| g.as_slice()).collect();
+        let q = quality(&groups, BETA);
+        assert!(q.is_finite(), "n={n}: quality must be finite, got {q}");
+        assert!((0.0..=1.0).contains(&q), "n={n}: out of range: {q}");
+    }
+}
+
+/// Equal-strength groups are the best-matched case: introducing a skill gap
+/// must lower quality.
+#[test]
+fn imbalance_lowers_quality() {
+    let strong = rating(40.0, 3.0);
+    let average = rating(25.0, 3.0);
+
+    let balanced = quality(&[&[average], &[average], &[average]], BETA);
+    let lopsided = quality(&[&[strong], &[average], &[average]], BETA);
+
+    assert!(
+        lopsided < balanced,
+        "expected imbalanced quality {lopsided} < balanced {balanced}"
+    );
+}
+
+/// Quality is a property of the multiset of groups, not their order.
+#[test]
+fn quality_is_permutation_invariant() {
+    let a = rating(30.0, 2.0);
+    let b = rating(25.0, 3.0);
+    let c = rating(20.0, 4.0);
+
+    let forward = quality(&[&[a], &[b], &[c]], BETA);
+    let reversed = quality(&[&[c], &[b], &[a]], BETA);
+
+    assert!(
+        (forward - reversed).abs() < 1e-9,
+        "permutation changed quality: {forward} vs {reversed}"
+    );
+}
+
+#[test]
+fn multi_player_groups_work() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r, r], &[r, r], &[r, r]], BETA);
+    assert!(q.is_finite());
+    assert!((0.0..=1.0).contains(&q));
+}
+
+#[test]
+fn uneven_group_sizes_work() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r, r], &[r], &[r, r, r]], BETA);
+    assert!(q.is_finite(), "got {q}");
+    assert!((0.0..=1.0).contains(&q), "got {q}");
+}
+
+#[test]
+#[should_panic(expected = "at least 2 rating groups")]
+fn single_group_panics_with_clear_message() {
+    let r = rating(25.0, 3.0);
+    let _ = quality(&[&[r]], BETA);
+}
+
+#[test]
+#[should_panic(expected = "at least 2 rating groups")]
+fn zero_groups_panics_with_clear_message() {
+    let _ = quality(&[], BETA);
+}
+
+#[test]
+#[should_panic(expected = "non-empty")]
+fn empty_group_panics_with_clear_message() {
+    let r = rating(25.0, 3.0);
+    let _ = quality(&[&[r], &[]], BETA);
+}
+
+#[test]
+fn history_predict_quality_supports_three_teams() {
+    use trueskill_tt::History;
+
+    let mut h = History::default();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    h.record_winner(&"b", &"c", 2).unwrap();
+    h.converge().unwrap();
+
+    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]);
+    assert!(
+        q.is_finite(),
+        "3-team predict_quality must be finite, got {q}"
+    );
+    assert!((0.0..=1.0).contains(&q), "out of range: {q}");
+}