chore: Release trueskill-tt version 0.5.0

Two decisions taken before cutting 0.5.0.

#42 — the `Schedule` trait was public API the engine never called. Its
only call site was `Game::custom`, itself `#[doc(hidden)]`, and
`EpsilonOrMax` was never constructed anywhere. Removing that surface
showed the problem was larger than the issue described: with `custom`
gone, the compiler found `Factor`, `BuiltinFactor`, `RankDiffFactor` and
`TeamSumFactor` all dead too.

`Game::run_chain` drives a local `DiffFactor` enum and bypasses the
whole T1 abstraction — it has done since it was written. So this is not
just an unused extension point but the machinery it was built on, and
`CLAUDE.md` was documenting it as live architecture.

Removed: `graph` module, `Schedule`, `EpsilonOrMax`, `ScheduleReport`,
`Game::custom`, `Factor`, `BuiltinFactor`, `RankDiffFactor`,
`TeamSumFactor`. `TruncFactor`, `MarginFactor`, `VarStore` and `VarId`
stay — inference uses those. The measurement behind choosing removal
over wiring is in #42: the within-game loop converges in 1 to 8
iterations against a cap of 30, so a `Residual` schedule has no headroom
to reclaim, and `Damped` already shipped as `ConvergenceOptions::alpha`.

#20 — `Outcome::winner` panicking on an out-of-range index. Kept, and
the reasoning is now on the method. It is the only constructor here that
validates, which looks inconsistent until you try deferring like its
siblings: `winner(5, 2)` produces ranks `[1, 1]`, an all-tied draw that
ingestion accepts without complaint when `p_draw > 0`. Asking "team 5
won" and silently getting "everyone drew" is the exact failure this
crate keeps removing, so the check belongs where the mistake is.

Adds `Outcome::try_winner` for indices that are computed or parsed
rather than written literally, following the `new`/`try_new` convention.
That is additive; the panicking form stays because every call site in
this repo, its tests and its README passes literals, where a `?` would
be noise.

BREAKING CHANGE: the `graph` module and everything it exported are
removed, as are `Game::custom` and `ScheduleReport`.

Closes #42. Closes #20.

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

CHANGELOG.md

@@ -2,6 +2,95 @@
 
 All notable changes to this project will be documented in this file.
 
+## 0.5.0 - 2026-09-08
+
+### Breaking Changes
+
+- feat!: name the unknown key, expose tail probabilities, flag short fits
+- refactor!: remove the factor-graph surface nothing used, add try_winner
+
+### Bug Fixes
+
+- fix(test): the ingestion-order property was comparing two truncated fits
+
+### Documentation
+
+- docs: record that the event log is the source of truth, and why
+
+### Features
+
+- feat: add UnknownKeys::Prior, and explain why there is no Skip
+- feat: add History::posterior_of for a linear combination of competitors
+- feat: add History::predict_margin for scored matchups
+- feat: add expected_variance_reduction for scored active learning
+
+### Styling
+
+- style: factor the event-pair type out of the reconvergence fixture
+- style: use arrays rather than vec! in the calibration fixture
+
+### Testing
+
+- test: pin that re-convergence is path-independent
+- test: calibrate the marginals against the exact posterior
+- test: pin what an additive model does to combined uncertainty
+
+## 0.4.2 - 2026-09-07
+
+### Bug Fixes
+
+- fix: replace the erfc approximation with libm, for free
+- fix: route every transcendental through libm, and combine sigmas with hypot
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.4.2
+
+### Testing
+
+- test: localise the erfc_inv tail residual to the caller's argument
+
+## 0.4.1 - 2026-09-07
+
+### Bug Fixes
+
+- fix: correct erfc_inv's sign error and keep evidence in log space
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.4.1
+
+### Testing
+
+- test: pin quality()'s N-group closed form, closing the README cross-check
+
+## 0.4.0 - 2026-09-07
+
+### Breaking Changes
+
+- feat!: N-team outcome prediction with draw mass, replacing the 2-team panic
+- refactor!: close the remaining API gaps from #21
+- fix!: apply competitor configuration whenever it is supplied
+
+### Bug Fixes
+
+- fix(release): skip the changelog hook during a dry run
+- fix: stop destroying tail precision in evidence and truncation
+- fix: reject convergence options that silently disable inference
+
+### Documentation
+
+- docs: correct drifted documentation and compile the README in CI
+
+### Features
+
+- feat: add expected information gain for active matchup selection
+- feat: let observers be shared, boxed, or borrowed
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.4.0
+
 ## 0.3.0 - 2026-09-01
 
 ### Breaking Changes
@@ -25,6 +114,7 @@ All notable changes to this project will be documented in this file.
 ### Miscellaneous Tasks
 
 - chore: ignore proptest regression seed files
+- chore: Release trueskill-tt version 0.3.0
 
 ### Performance
 

CLAUDE.md

@@ -24,6 +24,20 @@ is where several defects have hidden — a debug-only run is not evidence.
 - `approx` — `approx::AbsDiffEq` etc. for `Gaussian`. Most numerical goldens need it.
 - `rayon` — opt-in parallel within-slice sweep and per-slice query passes.
 
+## Working rules
+
+- **Investigate before implementing.** Measure the actual behaviour first —
+  against an analytic reference where one exists. Several "obvious" fixes in
+  this repo turned out to be wrong in sign or unnecessary, and the measurement
+  is what caught them.
+- **Fix the root issue, not the symptom.** A clamp that hides an underflow, or
+  a tolerance loosened to make a test pass, is a defect deferred.
+- **Scout crates.io before hand-rolling numerics.** Check accuracy against an
+  independent reference rather than trusting downloads: `puruspe` has 1.4M
+  downloads and is 346 ULP off in the tail, where `libm` is 1. Fewer
+  dependencies is preferable, not mandatory — take the dependency when it is
+  measurably better.
+
 ## Architecture
 
 A Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py):
@@ -32,10 +46,19 @@ evidence both forward and backward across a history.
 
 ### Data flow
 
+Ingestion (public types, `event.rs`):
+
 ```
-History  →  TimeSlice[]  →  Event[]  →  Team[] → Item[]
-                                  ↓
-                     Game (factor graph) → Schedule → BuiltinFactor[]
+Event<T, K>  →  Team<K>[]  →  Member<K>[]
+```
+
+`History::add_events` flattens that into indices; teams survive only as
+grouping, not as a value. Inference then runs on the internal shapes:
+
+```
+History  →  TimeSlice[]  →  Event[]  →  Item[]
+                                ↓
+                   Game (factor graph) → Schedule → BuiltinFactor[]
 ```
 
 - **`History`** (`history.rs`) — top level. Interns keys, groups events into
@@ -45,30 +68,35 @@ History  →  TimeSlice[]  →  Event[]  →  Team[] → Item[]
 - **`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`.
+- **`Event`** — two distinct types, do not confuse them. The *public* ingestion
+  `Event<T, K>` is in `event.rs` (with `Team`/`Member`); the *internal*
+  `pub(crate) Event` in `time_slice.rs` is one match during inference, where
+  `compute()` runs inference reading skills immutably and `apply()` folds the
+  result back. That 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.
+- **`factor/`** — `TruncFactor` (ranked) and `MarginFactor` (scored) over a
+  flat `VarStore`. `Game::run_chain` drives them directly through a local
+  `DiffFactor` enum; there is no `Schedule` indirection and no generic `Factor`
+  trait. Both were removed once measurement showed nothing had ever used them
+  — see #42.
 - **`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`.
+- **`storage/`** — `SkillStore` (per slice, `pub(crate)`) and `CompetitorStore`
+  (per history, public), both indexed by `Index`. The module is `pub`, but only
+  `CompetitorStore` is reachable from outside the crate.
 - **`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()`.
+  `GAMMA`, `P_DRAW`, `EPSILON`, `ITERATIONS`), and the standalone `quality()`.
+  The `cdf()` / `erfc()` helpers live here too but are `pub(crate)` and private
+  respectively — not public API.
 
 ### Invariants worth knowing
 
@@ -83,6 +111,12 @@ History  →  TimeSlice[]  →  Event[]  →  Team[] → Item[]
   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.
+- **Transcendentals go through `libm`, not `std`.** IEEE 754 pins the basic
+  operations and `sqrt` but says nothing about `exp`/`log`/`erf`, and `std`
+  delegates to the *system* math library — measured, `f64::exp` and `libm::exp`
+  disagree on 9.7% of inputs by one ULP. Since inference is an iterative fixed
+  point, one ULP can change an iteration count. Use `libm::exp` / `libm::log` in
+  inference code; `f64::sqrt` is fine (IEEE specifies it). Tests may use either.
 - **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

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "trueskill-tt"
-version = "0.3.0"
+version = "0.5.0"
 edition = "2024"
 rust-version = "1.85"
 description = "TrueSkill Through Time: Bayesian skill rating that tracks how skill evolves over time, via Gaussian message passing"
@@ -51,6 +51,7 @@ harness = false
 
 [dependencies]
 approx = { version = "0.5.1", optional = true }
+libm = "0.2.16"
 rayon = { version = "1", optional = true }
 smallvec = "1"
 

README.md

@@ -13,63 +13,92 @@ Rust port of [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillTh
 
 ## Drift
 
-Skill drift models how a player's true skill can change between appearances. Each time a player reappears after a gap, their skill uncertainty is widened by the drift model before the new evidence is incorporated.
+Skill drift models how a competitor's true skill can change between appearances.
+Each time they reappear after a gap, their skill uncertainty is widened by the
+drift model before the new evidence is incorporated.
 
-Drift is represented by the `Drift` trait:
+Drift is represented by the `Drift` trait (`src/drift.rs`), generic over the
+history's time type:
 
-```rust
-pub trait Drift: Copy + Debug {
-    fn variance_delta(&self, elapsed: i64) -> f64;
+```text
+pub trait Drift<T: Time>: Copy + Debug + Send + Sync {
+    fn variance_delta(&self, from: &T, to: &T) -> f64;
+    fn variance_for_elapsed(&self, elapsed: i64) -> f64;
 }
 ```
 
-`variance_delta` returns the amount to add to `σ²` given the elapsed time since the player last played. Internally, `Gaussian::forget` uses this to compute the new sigma: `σ_new = sqrt(σ² + variance_delta)`.
+Both methods return the amount to add to `σ²`, not to `σ`. `variance_delta`
+works from two timestamps; `variance_for_elapsed` takes an already-computed
+elapsed count, and is used on the paths that cache it. `Gaussian::forget`
+applies the result entirely in variance space — `from_mv(mu, variance() +
+variance_delta)` — taking no square root.
+
+That block is a quotation rather than a doctest. The custom-drift example below
+is compiled by CI, so it is what actually pins the signature.
 
 ### ConstantDrift
 
-The built-in `ConstantDrift` implements a linear random walk — skill uncertainty grows proportionally to time:
+The built-in `ConstantDrift` implements a linear random walk — skill uncertainty
+grows proportionally to time:
 
-```
+```text
 variance_delta = elapsed * γ²
 ```
 
-This is the standard TrueSkill Through Time model. Use it by passing a `ConstantDrift(gamma)` when constructing a `Player`:
+This is the standard TrueSkill Through Time model. Pass a `ConstantDrift(gamma)`
+when constructing a `Rating`:
 
 ```rust
-use trueskill_tt::{Player, Gaussian, drift::ConstantDrift};
+use trueskill_tt::{ConstantDrift, Gaussian, Rating};
 
-// gamma = 0.1 means skill can shift ~0.1 per time unit
-let player = Player::new(Gaussian::from_ms(0.0, 6.0), 1.0, ConstantDrift(0.1));
+// gamma = 0.1 means skill can shift ~0.1 per time unit.
+let rating: Rating<i64, ConstantDrift> =
+    Rating::new(Gaussian::from_ms(0.0, 6.0), 1.0, ConstantDrift(0.1));
+
+assert_eq!(rating.drift().0, 0.1);
 ```
 
+The type annotation is load-bearing: `ConstantDrift` implements `Drift<T>` for
+every `T: Time`, so without it `T` is ambiguous.
+
 ### Custom drift
 
-Implement `Drift` to express any other model. For example, a drift that saturates after a long absence (uncertainty grows with the square root of elapsed time instead of linearly):
+Implement `Drift<T>` to express any other model. For example, a drift that
+saturates after a long absence, with uncertainty growing as the square root of
+elapsed time instead of linearly:
 
 ```rust
-use trueskill_tt::drift::Drift;
+use trueskill_tt::{Drift, Gaussian, History, Rating, Time};
 
 #[derive(Clone, Copy, Debug)]
 struct SqrtDrift {
     gamma: f64,
 }
 
-impl Drift for SqrtDrift {
-    fn variance_delta(&self, elapsed: i64) -> f64 {
-        (elapsed as f64).sqrt() * self.gamma * self.gamma
+impl<T: Time> Drift<T> for SqrtDrift {
+    fn variance_delta(&self, from: &T, to: &T) -> f64 {
+        let elapsed = from.elapsed_to(to).max(0) as f64;
+        elapsed.sqrt() * self.gamma * self.gamma
+    }
+
+    fn variance_for_elapsed(&self, elapsed: i64) -> f64 {
+        (elapsed.max(0) as f64).sqrt() * self.gamma * self.gamma
     }
 }
 
-let player = Player::new(Gaussian::from_ms(0.0, 6.0), 1.0, SqrtDrift { gamma: 0.5 });
+// On a single Rating:
+let rating: Rating<i64, SqrtDrift> =
+    Rating::new(Gaussian::from_ms(0.0, 6.0), 1.0, SqrtDrift { gamma: 0.5 });
+
+// Or for a whole History, via the builder:
+let history = History::builder().drift(SqrtDrift { gamma: 0.5 }).build();
+
+assert_eq!(rating.beta(), 1.0);
+assert_eq!(history.log_evidence(), 0.0);
 ```
 
-To use a custom drift type with `History`, use the `.drift()` builder method instead of `.gamma()`:
-
-```rust
-let h = History::builder()
-    .drift(SqrtDrift { gamma: 0.5 })
-    .build();
-```
+`HistoryBuilder::drift` is the only way to set a history's drift model; there is
+no `gamma()` shorthand. The default is `ConstantDrift(GAMMA)`.
 
 ### Per-competitor drift
 
@@ -84,16 +113,25 @@ expressible in the same graph as moving competitors — a bot at a known
 strength, a rating floor, a course difficulty:
 
 ```rust
-let events = vec![Event {
+use trueskill_tt::{ConstantDrift, Event, History, Member, Outcome, Team};
+
+let mut h = History::builder().drift(ConstantDrift(0.1)).build();
+
+h.add_events(vec![Event {
     time: 0,
-    teams: smallvec![
+    teams: [
         Team::with_members([Member::new("player")]),
         // A course does not improve. Pin it, and the round's evidence
         // lands on the player instead of being split between the two.
         Team::with_members([Member::new("layout_7").with_drift_scale(0.0)]),
-    ],
+    ]
+    .into_iter()
+    .collect(),
     outcome: Outcome::winner(0, 2),
-}];
+}])
+.unwrap();
+
+h.converge().unwrap();
 ```
 
 Like `with_prior`, the scale is **competitor configuration captured at first
@@ -101,6 +139,10 @@ appearance** — setting it on a key the history already knows has no effect. It
 must be finite and non-negative; ingestion otherwise fails with
 `InferenceError::InvalidParameter`.
 
+Note that the fluent `EventBuilder` (`h.event(t).team([...])`) sets weights but
+not `drift_scale` or `prior`; those need the typed `Event` / `Team` / `Member`
+shape shown above.
+
 ## Scored outcomes
 
 Use `Outcome::scores([...])` when you have continuous per-team scores rather
@@ -110,7 +152,7 @@ soft Gaussian evidence about the latent performance diff. Configure
 (smaller σ = more trust).
 
 ```rust
-use trueskill_tt::{History, Outcome};
+use trueskill_tt::History;
 
 let mut h = History::builder().score_sigma(2.0).build();
 h.event(1)
@@ -122,6 +164,114 @@ h.event(1)
 h.converge().unwrap();
 ```
 
+## Prediction
+
+`predict_outcome` gives the full distribution over finishing orders. Each entry
+is a rank vector in the same shape `Outcome::ranking` takes — equal ranks mean a
+tie — so an outcome feeds straight back into inference.
+
+```rust
+use trueskill_tt::History;
+
+let mut h = History::builder().p_draw(0.1).build();
+h.record_winner(&"alice", &"bob", 1).unwrap();
+h.converge().unwrap();
+
+let p = h.predict_outcome(&[&[&"alice"], &[&"bob"]]).unwrap();
+
+// Probabilities are exhaustive and disjoint, so they sum to one.
+assert!((p.total() - 1.0).abs() < 1e-6);
+
+let (best, likelihood) = p.most_likely().unwrap();
+println!("most likely: {best:?} at {likelihood:.3}");
+println!("draw:        {:.3}", p.probability_of(&[0, 0]));
+```
+
+Supports any number of teams. Because the outcome space grows factorially, the
+full distribution is capped at `MAX_PREDICTED_TEAMS`; two cheaper entry points
+stay available at any size:
+
+- `predict_win_probabilities(teams)` — `P(team i finishes strictly first)`,
+  quadratic in team count.
+- `predict_ranking(teams, ranks)` — one specific finishing order.
+
+Unknown keys are an error by default, not a silent omission: a team the history
+has never seen cannot produce a confident-looking probability. The error names
+the key, and every key must already be known — pre-filter with `lookup` or
+`current_skill` if your caller cannot guarantee that.
+
+If predicting for competitors you have never seen is the point rather than a
+mistake, say so once:
+
+```rust
+use trueskill_tt::{History, UnknownKeys};
+
+let h = History::builder().unknown_keys(UnknownKeys::Prior).build();
+```
+
+An unknown competitor is then answered from the configured prior, which is the
+honest reading — you have no evidence about them — and correctly *widens* a team
+that contains one. There is deliberately no "skip the member" mode: a team's
+performance is the sum of its members, so dropping one would make the model more
+certain because it knows less.
+
+### Asking about one competitor
+
+`Gaussian` answers tail questions directly, which is what a stopping rule needs:
+
+```rust
+use trueskill_tt::History;
+
+let mut h = History::builder().build();
+h.record_winner(&"alice", &"bob", 1).unwrap();
+let _ = h.converge().unwrap();
+
+let skill = h.current_skill(&"alice").unwrap();
+
+// "How sure am I that this is below the cutoff?" — a probability, not a
+// `mu + z * sigma` band whose confidence drifts as sigma changes.
+let _ = skill.probability_below(20.0);
+
+// Use this rather than `1.0 - probability_below(x)`: the complement cancels
+// away every digit in the upper tail, which is where a stopping rule lives.
+let _ = skill.probability_above(30.0);
+```
+
+## Which match to play next
+
+`quality()` measures whether a matchup is *fair*. That is not the same as
+whether it is *informative*, and the two only coincide for two evenly matched
+competitors. When each observation costs something, ask
+`expected_information_gain` instead — the outcome-weighted divergence between
+what you believe now and what you would believe afterwards.
+
+```rust
+use trueskill_tt::History;
+
+let mut h = History::builder().build();
+for t in 1..=10 {
+    h.record_winner(&"veteran", &"regular", t).unwrap();
+    h.record_winner(&"regular", &"veteran", t + 100).unwrap();
+}
+h.record_winner(&"veteran", &"newcomer", 500).unwrap();
+h.converge().unwrap();
+
+let settled = h.expected_information_gain(&[&[&"veteran"], &[&"regular"]]).unwrap();
+let unknown = h.expected_information_gain(&[&[&"veteran"], &[&"newcomer"]]).unwrap();
+
+// Playing the newcomer teaches you more than replaying a settled rivalry.
+assert!(unknown > settled);
+```
+
+The result is in nats, and is bounded by the entropy of the outcome: at most
+`ln 2 ≈ 0.693` for a two-way result, `ln 3` once draws are possible, `ln k` for
+`k` outcomes. A value near zero means you already know how it ends.
+
+This costs one full inference pass **per possible outcome**, so it is far more
+expensive than `quality()`. Scoring every pairing among `n` competitors is
+`O(n² × outcomes)` passes — shortlist with `quality()` or
+`predict_win_probabilities` first, then score only the shortlist.
+
 ## Todo
 
 - [x] Implement approx for Gaussian
@@ -130,7 +280,8 @@ h.converge().unwrap();
 - [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
+- [x] N-team `predict_outcome` with draw mass, and `expected_information_gain`
+- [x] Cross-check `quality()` against [sublee/trueskill](https://github.com/sublee/trueskill/tree/master) — N identical teams follow the closed form `(1/5)^((n-1)/2)` for the conventional parameters, asserted for n = 2..10, and the n=3/n=5 values (0.200, 0.040) match the reference package
 
 ## License
 

benches/history_converge.rs

@@ -82,7 +82,7 @@ fn bench_converge(c: &mut Criterion) {
         b.iter_batched(
             || build_history_1v1(500, 100, 10, 42),
             |mut h| {
-                h.converge().unwrap();
+                let _ = h.converge().unwrap();
             },
             BatchSize::SmallInput,
         );
@@ -92,7 +92,7 @@ fn bench_converge(c: &mut Criterion) {
         b.iter_batched(
             || build_history_1v1(2000, 200, 20, 42),
             |mut h| {
-                h.converge().unwrap();
+                let _ = h.converge().unwrap();
             },
             BatchSize::SmallInput,
         );
@@ -106,7 +106,7 @@ fn bench_converge(c: &mut Criterion) {
         b.iter_batched(
             || build_history_1v1(5000, 50000, 5000, 42),
             |mut h| {
-                h.converge().unwrap();
+                let _ = h.converge().unwrap();
             },
             BatchSize::SmallInput,
         );

benches/scored.rs

@@ -29,7 +29,7 @@ fn bench_scored_history(c: &mut Criterion) {
                 });
             }
             h.add_events(events).unwrap();
-            h.converge().unwrap();
+            let _ = h.converge().unwrap();
         });
     });
 }

examples/atp.rs

@@ -46,14 +46,33 @@ fn main() {
         .sigma(1.6)
         .drift(ConstantDrift(0.036))
         .convergence(trueskill_tt::ConvergenceOptions {
-            max_iter: 10,
+            // This history needs 30 sweeps to reach the epsilon below. It was
+            // capped at 10 until the `#[must_use]` on `ConvergenceReport`
+            // surfaced that the example had been shipping a short fit.
+            max_iter: 100,
             epsilon: 0.01,
             alpha: 1.0,
         })
         .build();
 
     hist.add_events(events).unwrap();
-    hist.converge().unwrap();
+
+    // Read the report rather than discarding it. A fit that hits `max_iter`
+    // without reaching `epsilon` is not an error and does not look wrong — every
+    // rating comes back finite and sensibly ordered — so this flag is the only
+    // thing that says the numbers were still moving when the sweep stopped.
+    let report = hist.converge().unwrap();
+    eprintln!(
+        "converged={} after {} sweeps, final step {:?}",
+        report.converged, report.iterations, report.final_step
+    );
+    if !report.converged {
+        eprintln!(
+            "warning: stopped after {} sweeps with a final step of {:?}, \
+             short of epsilon — raise ConvergenceOptions::max_iter",
+            report.iterations, report.final_step
+        );
+    }
 
     let players = [
         ("aggasi", "a092", 38800i64),

release.toml

@@ -3,4 +3,12 @@ publish = true
 # Hold off pushing until tags and publish have both succeeded; `just release`
 # pushes last.
 push = false
-pre-release-hook = ["sh", "-c", "git cliff -o CHANGELOG.md --tag {{version}} && git add CHANGELOG.md"]
+# Regenerate the changelog and stage it so it lands in the release commit.
+#
+# Guarded on DRY_RUN because cargo-release runs pre-release hooks during a dry
+# run too (verified against cargo-release 1.1.5, which exports DRY_RUN=true,
+# CRATE_NAME, PREV_VERSION and NEW_VERSION to the hook). Without the guard,
+# `just release-plan` — documented as a preview that writes nothing — writes and
+# `git add`s CHANGELOG.md, and the clean-tree check in `just release` then
+# refuses to run. That check is load-bearing: publishing is irreversible.
+pre-release-hook = ["sh", "-c", '[ "$DRY_RUN" = "true" ] || (git cliff -o CHANGELOG.md --tag {{version}} && git add CHANGELOG.md)']

src/acquisition.rs

@@ -0,0 +1,352 @@
+//! Active learning: which comparison teaches you the most.
+//!
+//! [`quality`](crate::quality) answers "is this matchup *fair*". That is a
+//! different question from "is this matchup *informative*", and the two
+//! coincide only for two evenly matched competitors. When each observation
+//! costs something — a human click, a scheduled fixture — the question worth
+//! asking is the second one.
+//!
+//! The quantity here is expected information gain: the outcome-weighted
+//! divergence between what you believe now and what you would believe after
+//! seeing the result.
+//!
+//! ```text
+//! EIG(matchup) = SUM  P(outcome) * KL( posterior_after(outcome) || prior )
+//!              outcome
+//! ```
+//!
+//! It is the mutual information between the observed outcome and the skills,
+//! which is worth remembering because it pins the scale: information gain
+//! cannot exceed the entropy of the thing you are about to observe. A contest
+//! with `k` distinguishable outcomes can teach you at most `ln k` nats,
+//! whatever the ratings. That ceiling is the sharpest available test of an
+//! implementation — see [`expected_information_gain`].
+
+use crate::{
+    GameOptions, Gaussian, InferenceError, Outcome, Rating, drift::Drift, predict, time::Time,
+};
+
+/// Outcomes below this probability contribute nothing measurable and are not
+/// worth an inference pass.
+///
+/// The contribution of an outcome is `P * KL`, and `KL` is bounded in practice
+/// by tens of nats, so a probability this small moves the total by less than
+/// the quadrature error already present in `P` itself.
+const NEGLIGIBLE: f64 = 1e-12;
+
+/// `KL(q || p)` for two univariate Gaussians, in nats.
+///
+/// Both arguments are proper posteriors from inference, so the degenerate
+/// cases guarded here (zero or infinite variance) indicate that inference has
+/// broken down rather than anything a caller did.
+fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 {
+    let (var_q, var_p) = (q.sigma().powi(2), p.sigma().powi(2));
+
+    if !(var_q.is_finite() && var_p.is_finite()) || var_q <= 0.0 || var_p <= 0.0 {
+        return 0.0;
+    }
+
+    let mean_gap = q.mu() - p.mu();
+    0.5 * (libm::log(var_p / var_q) + (var_q + mean_gap * mean_gap) / var_p - 1.0)
+}
+
+/// Expected information gain of a hypothetical matchup, in nats.
+///
+/// Enumerates the outcomes this matchup could have, runs inference for each to
+/// get the belief it would produce, and weights the resulting divergence by
+/// that outcome's probability. A higher value means the result would teach you
+/// more.
+///
+/// # Interpreting the value
+///
+/// Nats. The upper bound is the entropy of the outcome variable: at most
+/// `ln 2 ≈ 0.693` for a two-way result, `ln 3 ≈ 1.099` once draws are
+/// possible, `ln k` for `k` outcomes. A value near the ceiling means the
+/// result is close to a coin flip *and* would move the posteriors a long way;
+/// a value near zero means you already know what will happen, or that the
+/// result would barely change your beliefs if you saw it.
+///
+/// This is not a monotone transform of [`quality`](crate::quality). A lopsided
+/// matchup between two uncertain competitors scores well on quality-times-
+/// variance heuristics and poorly here, because the near-certain outcome
+/// carries almost no information.
+///
+/// # Cost
+///
+/// One full inference pass per possible outcome, so this is far more expensive
+/// than `quality()` — which is one closed-form evaluation. The outcome count
+/// grows quickly with team count (3 outcomes for two teams that can draw, 13
+/// for three, 75 for four), and scoring every candidate pairing among `n`
+/// competitors is `O(n² × outcomes)` inference passes.
+///
+/// For a selector over many candidates, shortlist with the cheap
+/// [`quality`](crate::quality) or
+/// [`predict_win_probabilities`](crate::History::predict_win_probabilities)
+/// first and score only the shortlist here. The expected-variance-reduction
+/// proxy sometimes suggested as a cheaper alternative is *not* cheaper: it
+/// needs the same hypothetical posteriors, so it shares the dominant cost.
+///
+/// # Errors
+///
+/// - `NotEnoughTeams` if fewer than two teams are supplied.
+/// - `EmptyTeam` if any team has no members.
+/// - `TooManyTeams` if the outcome space is too large to enumerate; see
+///   [`MAX_PREDICTED_TEAMS`](crate::MAX_PREDICTED_TEAMS).
+/// - `InvalidProbability` if `options.p_draw` is outside `[0.0, 1.0)`.
+/// - Anything [`Game::ranked`](crate::Game::ranked) returns for a hypothetical
+///   outcome.
+pub fn expected_information_gain<T: Time, D: Drift<T>>(
+    teams: &[&[Rating<T, D>]],
+    options: &GameOptions,
+) -> Result<f64, InferenceError> {
+    if teams.len() < 2 {
+        return Err(InferenceError::NotEnoughTeams { got: teams.len() });
+    }
+    if teams.len() > crate::MAX_PREDICTED_TEAMS {
+        return Err(InferenceError::TooManyTeams {
+            got: teams.len(),
+            max: crate::MAX_PREDICTED_TEAMS,
+        });
+    }
+    if !(0.0..1.0).contains(&options.p_draw) {
+        return Err(InferenceError::InvalidProbability {
+            value: options.p_draw,
+        });
+    }
+    for (idx, team) in teams.iter().enumerate() {
+        if team.is_empty() {
+            return Err(InferenceError::EmptyTeam { team: idx });
+        }
+    }
+
+    // Prediction runs on performances: skill inflated by each member's beta.
+    let performances: Vec<Gaussian> = teams
+        .iter()
+        .map(|team| {
+            team.iter()
+                .fold(crate::N00, |acc, rating| acc + rating.performance())
+        })
+        .collect();
+
+    // Draw margins per pair, derived from the teams' betas exactly as
+    // inference derives them, so the outcomes weighted here are the outcomes
+    // that would actually be fitted.
+    let beta_sq: Vec<f64> = teams
+        .iter()
+        .map(|team| team.iter().map(|r| r.beta().powi(2)).sum())
+        .collect();
+    let p_draw = options.p_draw;
+    let margins = predict::Margins::new(teams.len(), |i, j| {
+        if p_draw == 0.0 {
+            0.0
+        } else {
+            crate::compute_margin(p_draw, (beta_sq[i] + beta_sq[j]).sqrt())
+        }
+    });
+
+    let mut gain = 0.0;
+
+    for (ranks, probability) in predict::outcome_distribution(&performances, &margins) {
+        if probability <= NEGLIGIBLE {
+            continue;
+        }
+
+        let game = crate::Game::ranked(teams, Outcome::ranking(ranks), options)?;
+        let posteriors = game.posteriors();
+
+        // Beliefs factorise across competitors, so the joint divergence is the
+        // sum of the per-competitor ones.
+        let divergence: f64 = teams
+            .iter()
+            .zip(&posteriors)
+            .flat_map(|(team, posterior)| team.iter().zip(posterior))
+            .map(|(rating, &after)| kl_divergence(after, rating.prior()))
+            .sum();
+
+        gain += probability * divergence;
+    }
+
+    Ok(gain)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{BETA, ConstantDrift, GAMMA};
+
+    type R = Rating<i64, ConstantDrift>;
+
+    fn rating(mu: f64, sigma: f64) -> R {
+        R::new(Gaussian::from_ms(mu, sigma), BETA, ConstantDrift(GAMMA))
+    }
+
+    fn options(p_draw: f64) -> GameOptions {
+        GameOptions {
+            p_draw,
+            ..GameOptions::default()
+        }
+    }
+
+    fn eig(teams: &[&[R]], p_draw: f64) -> f64 {
+        expected_information_gain(teams, &options(p_draw)).unwrap()
+    }
+
+    /// The analytic ceiling. Information gain is the mutual information between
+    /// the outcome and the skills, so it cannot exceed the entropy of the
+    /// outcome variable — whatever the ratings. This is the check a subtly
+    /// wrong implementation fails while still returning plausible numbers: an
+    /// early prototype of this returned 4.77 nats from a sign error and passed
+    /// every monotonicity test.
+    #[test]
+    fn never_exceeds_the_entropy_of_the_outcome() {
+        let ceiling_two = std::f64::consts::LN_2;
+
+        for (a, b) in [
+            (rating(0.0, 6.0), rating(0.0, 6.0)),
+            (rating(0.0, 0.5), rating(0.0, 0.5)),
+            (rating(12.0, 6.0), rating(-12.0, 6.0)),
+            (rating(40.0, 1.0), rating(-40.0, 1.0)),
+            (rating(3.0, 6.0), rating(-2.0, 0.1)),
+            (rating(0.0, 25.0), rating(0.0, 25.0)),
+        ] {
+            let g = eig(&[&[a], &[b]], 0.0);
+            assert!(
+                g >= 0.0 && g <= ceiling_two,
+                "EIG {g} outside [0, ln 2] for mu=({}, {}) sigma=({}, {})",
+                a.prior().mu(),
+                b.prior().mu(),
+                a.prior().sigma(),
+                b.prior().sigma()
+            );
+        }
+    }
+
+    /// With draws enabled there are three outcomes, so the ceiling rises to
+    /// `ln 3` — and the two-outcome bound no longer applies.
+    #[test]
+    fn the_ceiling_follows_the_outcome_count() {
+        let ceiling_three = 3.0f64.ln();
+        for sigma in [0.5, 3.0, 6.0, 25.0] {
+            let g = eig(&[&[rating(0.0, sigma)], &[rating(0.0, sigma)]], 0.25);
+            assert!(
+                g >= 0.0 && g <= ceiling_three,
+                "EIG {g} outside [0, ln 3] at sigma {sigma}"
+            );
+        }
+    }
+
+    /// An even matchup between uncertain competitors is the informative one.
+    /// A hopelessly lopsided matchup teaches you almost nothing, because you
+    /// already know how it ends.
+    #[test]
+    fn an_even_matchup_beats_a_lopsided_one() {
+        let even = eig(&[&[rating(0.0, 6.0)], &[rating(0.0, 6.0)]], 0.0);
+        let lopsided = eig(&[&[rating(12.0, 6.0)], &[rating(-12.0, 6.0)]], 0.0);
+        assert!(
+            even > lopsided,
+            "even {even} should beat lopsided {lopsided}"
+        );
+    }
+
+    /// Certainty is the thing information gain is measuring the absence of:
+    /// the less you know, the more there is to learn.
+    #[test]
+    fn gain_falls_as_certainty_rises() {
+        let mut previous = f64::INFINITY;
+        for sigma in [12.0, 6.0, 3.0, 1.0, 0.5, 0.1] {
+            let g = eig(&[&[rating(0.0, sigma)], &[rating(0.0, sigma)]], 0.0);
+            assert!(
+                g < previous,
+                "sigma {sigma}: {g} did not fall below {previous}"
+            );
+            previous = g;
+        }
+        assert!(previous >= 0.0);
+    }
+
+    /// The heuristic this replaces is `quality * sigma_a^2 * sigma_b^2`. It is
+    /// not a monotone transform of information gain — it ranks a lopsided
+    /// matchup above a confident even one, and EIG ranks them the other way.
+    /// Pinning the disagreement down is what stops a future "simplification"
+    /// from quietly reverting to the heuristic.
+    #[test]
+    fn disagrees_with_the_quality_times_variance_heuristic() {
+        let heuristic = |a: &R, b: &R| {
+            crate::quality(&[&[a.prior()], &[b.prior()]], BETA)
+                * a.prior().sigma().powi(2)
+                * b.prior().sigma().powi(2)
+        };
+
+        let (confident_a, confident_b) = (rating(0.0, 0.5), rating(0.0, 0.5));
+        let (lopsided_a, lopsided_b) = (rating(12.0, 6.0), rating(-12.0, 6.0));
+
+        assert!(
+            heuristic(&lopsided_a, &lopsided_b) > heuristic(&confident_a, &confident_b),
+            "the heuristic should prefer the lopsided matchup"
+        );
+        assert!(
+            eig(&[&[confident_a], &[confident_b]], 0.0) > eig(&[&[lopsided_a], &[lopsided_b]], 0.0),
+            "information gain should prefer the even matchup"
+        );
+    }
+
+    #[test]
+    fn supports_more_than_two_teams() {
+        let teams: Vec<Vec<R>> = vec![
+            vec![rating(0.0, 6.0)],
+            vec![rating(0.0, 6.0)],
+            vec![rating(0.0, 6.0)],
+        ];
+        let refs: Vec<&[R]> = teams.iter().map(Vec::as_slice).collect();
+        let g = expected_information_gain(&refs, &options(0.0)).unwrap();
+        // Six distinguishable orderings with no draws.
+        assert!(
+            g > 0.0 && g <= 6.0f64.ln(),
+            "three-team EIG {g} out of range"
+        );
+    }
+
+    #[test]
+    fn multi_member_teams_are_supported() {
+        let a = [rating(0.0, 6.0), rating(1.0, 4.0)];
+        let b = [rating(0.0, 6.0)];
+        let g = expected_information_gain(&[&a, &b], &options(0.0)).unwrap();
+        assert!(g > 0.0 && g <= std::f64::consts::LN_2, "{g}");
+    }
+
+    #[test]
+    fn degenerate_shapes_are_errors() {
+        let a = [rating(0.0, 6.0)];
+        assert!(matches!(
+            expected_information_gain(&[&a], &options(0.0)),
+            Err(InferenceError::NotEnoughTeams { got: 1 })
+        ));
+        let empty: [R; 0] = [];
+        assert!(matches!(
+            expected_information_gain(&[&a, &empty], &options(0.0)),
+            Err(InferenceError::EmptyTeam { team: 1 })
+        ));
+        assert!(matches!(
+            expected_information_gain(&[&a, &a], &options(1.5)),
+            Err(InferenceError::InvalidProbability { .. })
+        ));
+    }
+
+    #[test]
+    fn kl_divergence_is_zero_for_identical_beliefs() {
+        let g = Gaussian::from_ms(3.0, 2.0);
+        assert!(kl_divergence(g, g).abs() < 1e-15);
+    }
+
+    #[test]
+    fn kl_divergence_is_non_negative_and_grows_with_separation() {
+        let prior = Gaussian::from_ms(0.0, 3.0);
+        let mut previous = 0.0;
+        for mu in [0.0, 0.5, 1.0, 2.0, 4.0] {
+            let d = kl_divergence(Gaussian::from_ms(mu, 3.0), prior);
+            assert!(d >= 0.0, "negative divergence at mu {mu}: {d}");
+            assert!(d >= previous, "not increasing at mu {mu}");
+            previous = d;
+        }
+    }
+}

src/competitor.rs

@@ -7,8 +7,8 @@ use crate::{
 
 /// Per-history, temporal state for someone competing.
 ///
-/// Renamed from `Agent` in T2; the former `.player` field is now
-/// `.rating` to match the `Player → Rating` rename.
+/// The mutable half of a competitor: `Rating` holds their static
+/// configuration, this holds what inference learns as it sweeps.
 #[derive(Debug)]
 pub struct Competitor<T: Time = i64, D: Drift<T> = ConstantDrift> {
     pub rating: Rating<T, D>,

src/convergence.rs

@@ -20,6 +20,37 @@ pub struct ConvergenceOptions {
     pub alpha: f64,
 }
 
+impl ConvergenceOptions {
+    /// Reject values that would make inference silently meaningless.
+    ///
+    /// `HistoryBuilder::convergence` asserts these eagerly, but the fields are
+    /// public and `GameOptions` carries a `ConvergenceOptions` — so a caller
+    /// can hand `Game::ranked` a set the builder never saw. In release the
+    /// engine's `debug_assert!`s are gone, and an `alpha` of zero leaves every
+    /// EP update unapplied: inference returns the priors, with every likelihood
+    /// uninformative and nothing to indicate anything went wrong.
+    ///
+    /// # Errors
+    ///
+    /// `InvalidParameter` if `alpha` is outside `(0.0, 1.0]` or `epsilon` is
+    /// negative. NaN fails both comparisons and is rejected.
+    pub(crate) fn validate(&self) -> Result<(), crate::InferenceError> {
+        if !(self.alpha > 0.0 && self.alpha <= 1.0) {
+            return Err(crate::InferenceError::InvalidParameter {
+                name: "alpha",
+                value: self.alpha,
+            });
+        }
+        if self.epsilon.is_nan() || self.epsilon < 0.0 {
+            return Err(crate::InferenceError::InvalidParameter {
+                name: "epsilon",
+                value: self.epsilon,
+            });
+        }
+        Ok(())
+    }
+}
+
 impl Default for ConvergenceOptions {
     fn default() -> Self {
         Self {
@@ -32,6 +63,9 @@ impl Default for ConvergenceOptions {
 
 /// Post-hoc summary of a `History::converge` call.
 #[derive(Clone, Debug)]
+#[must_use = "a ConvergenceReport carries `converged`, and a fit that stopped \
+              at `max_iter` is wrong by a little rather than loudly broken — \
+              check it, or bind it to `_` to say you have decided not to"]
 pub struct ConvergenceReport {
     pub iterations: usize,
     pub final_step: (f64, f64),

src/error.rs

@@ -1,5 +1,44 @@
 use std::fmt;
 
+/// How a prediction should treat a key the history has never seen.
+///
+/// Configured once per history via
+/// [`HistoryBuilder::unknown_keys`](crate::HistoryBuilder::unknown_keys).
+/// Neither known consumer wants this to vary between queries — one predicts
+/// thousands of candidate matchups in a loop, the other's headline feature is
+/// predicting a competitor nobody has faced — so it is a property of how you
+/// intend to use the model rather than an argument on five call sites.
+///
+/// # There is deliberately no `Skip`
+///
+/// Dropping an unknown member is the obvious third option and it is wrong. A
+/// team's performance is the *sum* of its members, so removing one removes its
+/// variance too: measured on a two-member team with one unknown, skipping gives
+/// a performance sigma of 2.37 where treating the member as unknown gives 6.53.
+/// An unknown competitor would make the model *more* certain, which is
+/// backwards. `Prior` is also the answer the model already gives for a
+/// competitor it knows about but has no evidence for, so it corresponds to a
+/// state the model can actually be in; skipping does not.
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
+#[non_exhaustive]
+pub enum UnknownKeys {
+    /// Reject the prediction with [`InferenceError::UnknownKey`].
+    ///
+    /// The default, and the right one when every key is expected to be known:
+    /// a team of strangers should not silently produce a confident-looking
+    /// answer.
+    #[default]
+    Reject,
+    /// Treat an unknown competitor as one sitting at the history's configured
+    /// prior.
+    ///
+    /// This is the honest Bayesian reading — a competitor you have never
+    /// observed is exactly the prior — and it makes "predict a matchup
+    /// involving someone new" a first-class question rather than something a
+    /// caller fakes with a neutral constant.
+    Prior,
+}
+
 #[derive(Debug, Clone, PartialEq)]
 #[non_exhaustive]
 pub enum InferenceError {
@@ -25,11 +64,6 @@ pub enum InferenceError {
     /// 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
@@ -38,8 +72,49 @@ pub enum InferenceError {
         context: &'static str,
         step: (f64, f64),
     },
-    /// Negative precision: a Gaussian with `pi < 0` slipped into an API call.
-    NegativePrecision { pi: f64 },
+    /// One batch declared two different values for the same competitor's
+    /// configuration.
+    ///
+    /// `prior` and `drift_scale` configure a competitor, not an event, so a
+    /// batch that sets one of them twice with different values has no
+    /// well-defined meaning: events within a batch are not ordered, so
+    /// "last one wins" would make the result depend on iteration order.
+    /// Declaring the same value repeatedly is fine and is the expected shape
+    /// when a competitor's configuration is a property of the domain.
+    ConflictingCompetitorConfig {
+        competitor: usize,
+        field: &'static str,
+    },
+    /// A prediction referenced a key the history has no skill for.
+    ///
+    /// Reported rather than skipped: dropping unknown keys turns a team of
+    /// strangers into a confident-looking probability about nobody.
+    ///
+    /// `key` is the offending key's `Debug` rendering. It is carried because
+    /// the indices alone are not actionable: a caller that logs
+    /// `UnknownKey { team: 0, member: 0 }` learns nothing about *which* of its
+    /// keys the history has not seen, and the natural handling — fall back to a
+    /// neutral value — turns the whole thing into a plausible constant.
+    UnknownKey {
+        team: usize,
+        member: usize,
+        key: String,
+    },
+    /// A prediction was given a team with no members.
+    EmptyTeam { team: usize },
+    /// A joint posterior was requested where one cannot be formed exactly.
+    JointUnavailable { reason: &'static str },
+    /// Fewer than two teams were supplied to a prediction.
+    NotEnoughTeams { got: usize },
+    /// The full outcome distribution was requested for too many teams.
+    ///
+    /// Each realisation sorts into exactly one (order, tie-pattern) event, so
+    /// the space holds `n! * 2^(n-1)` members — 1_920 at five teams, 23_040 at
+    /// six, 322_560 at seven. Past `max` this stops being something to
+    /// enumerate on a caller's behalf; ask for individual rankings with
+    /// `predict_ranking`, or for `predict_win_probabilities`, both of which
+    /// stay cheap at any team count.
+    TooManyTeams { got: usize, max: usize },
 }
 
 impl fmt::Display for InferenceError {
@@ -78,17 +153,35 @@ impl fmt::Display for InferenceError {
             Self::InvalidParameter { name, value } => {
                 write!(f, "{name} is invalid: {value}")
             }
-            Self::ConvergenceFailed {
-                last_step,
-                iterations,
-            } => {
+            Self::ConflictingCompetitorConfig { competitor, field } => {
                 write!(
                     f,
-                    "convergence failed after {iterations} iterations; last step = {last_step:?}"
+                    "competitor {competitor}: this batch sets {field} to two different values"
                 )
             }
-            Self::NegativePrecision { pi } => {
-                write!(f, "precision must be non-negative; got {pi}")
+            Self::UnknownKey { team, member, key } => {
+                write!(
+                    f,
+                    "team {team}, member {member}: no skill recorded for key {key} \
+                     (every key must already be known to the history; pre-filter \
+                     with `lookup` or `current_skill` if that is not guaranteed)"
+                )
+            }
+            Self::EmptyTeam { team } => {
+                write!(f, "team {team} has no members")
+            }
+            Self::JointUnavailable { reason } => {
+                write!(f, "no exact joint posterior is available: {reason}")
+            }
+            Self::NotEnoughTeams { got } => {
+                write!(f, "prediction needs at least 2 teams, got {got}")
+            }
+            Self::TooManyTeams { got, max } => {
+                write!(
+                    f,
+                    "the outcome distribution over {got} teams is too large to enumerate (limit {max}); \
+                     use predict_ranking or predict_win_probabilities instead"
+                )
             }
         }
     }

src/event.rs

@@ -1,8 +1,10 @@
 //! Typed event description for bulk ingestion.
 //!
-//! `Event<T, K>` is the new public event shape (spec Section 4). Replaces
-//! the nested `Vec<Vec<Vec<Index>>>`, `Vec<Vec<f64>>`, `Vec<Vec<Vec<f64>>>`
-//! that the old `add_events_with_prior` took.
+//! `Event<T, K>` is the public event shape taken by `History::add_events`. It
+//! is a typed front end, not a replacement: `add_events` flattens it into the
+//! nested `Vec<Vec<Vec<Index>>>` / `Vec<Vec<f64>>` / `Vec<Vec<Vec<f64>>>` that
+//! the internal `add_events_with_prior` chokepoint still takes, and which
+//! `record_winner` and `record_draw` also route through.
 
 use smallvec::SmallVec;
 
@@ -48,9 +50,17 @@ impl<K> Default for Team<K> {
 /// `weight` applies per event and defaults to 1.0.
 ///
 /// `prior` and `drift_scale` are **competitor configuration**, not per-event
-/// values: both are captured when the competitor is first created and ignored
-/// on every later appearance. Setting either on a key the history already knows
-/// has no effect.
+/// values. Setting either applies to the competitor for the whole history, not
+/// just to this event, and applies whenever it is supplied — including on a key
+/// the history already knows. Because configuration lives on the competitor and
+/// `converge` refits from competitor state, configuring one late still refits
+/// the whole history rather than taking effect only from that event onward.
+///
+/// Repeating the same value is inert, which is the expected shape when the
+/// configuration is a property of the domain. Supplying two *different* values
+/// for one competitor within a single batch is
+/// `InferenceError::ConflictingCompetitorConfig`: events in a batch have no
+/// order, so there would be no well-defined winner.
 #[derive(Clone, Debug)]
 pub struct Member<K> {
     pub key: K,

src/event_builder.rs

@@ -107,7 +107,10 @@ where
     /// Set explicit per-team continuous scores with a per-event noise override.
     ///
     /// `sigma` overrides `HistoryBuilder::score_sigma` for this event only.
-    /// Must be `> 0.0`; debug-asserts otherwise via `Outcome::scores_with_sigma`.
+    /// Must be `> 0.0`. Constructing the 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 from `commit` rather than a panic.
     pub fn scores_with_sigma<I: IntoIterator<Item = f64>>(mut self, scores: I, sigma: f64) -> Self {
         self.event.outcome = crate::Outcome::scores_with_sigma(scores, sigma);
         self

src/factor/margin.rs

@@ -1,8 +1,8 @@
 use crate::{
     N_INF,
-    factor::{Factor, VarId, VarStore},
+    factor::{VarId, VarStore},
     gaussian::Gaussian,
-    pdf,
+    ln_pdf,
 };
 
 /// Gaussian observation factor on a diff variable.
@@ -16,7 +16,7 @@ pub struct MarginFactor {
     pub m_obs: f64,
     pub sigma: f64,
     pub(crate) msg: Gaussian,
-    pub(crate) evidence_cached: Option<f64>,
+    pub(crate) log_evidence_cached: Option<f64>,
 }
 
 impl MarginFactor {
@@ -28,7 +28,7 @@ impl MarginFactor {
             m_obs,
             sigma,
             msg: N_INF,
-            evidence_cached: None,
+            log_evidence_cached: None,
         }
     }
 }
@@ -41,8 +41,8 @@ impl MarginFactor {
         let marginal = vars.get(self.diff);
         let cavity = marginal / self.msg;
 
-        if self.evidence_cached.is_none() {
-            self.evidence_cached = Some(cavity_evidence(cavity, self.m_obs, self.sigma));
+        if self.log_evidence_cached.is_none() {
+            self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.m_obs, self.sigma));
         }
 
         let new_msg = Gaussian::from_ms(self.m_obs, self.sigma);
@@ -55,23 +55,42 @@ impl MarginFactor {
     }
 }
 
-impl Factor for MarginFactor {
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
+/// Undamped wrappers, used by this module's tests. Inference drives these
+/// factors through `propagate_with_alpha` and reads the cached log evidence
+/// directly, so these are not on any production path.
+#[cfg(test)]
+impl MarginFactor {
+    pub(crate) fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
         self.propagate_with_alpha(vars, 1.0)
     }
 
-    fn log_evidence(&self, _vars: &VarStore) -> f64 {
-        self.evidence_cached.unwrap_or(1.0).ln()
+    pub(crate) fn log_evidence(&self) -> f64 {
+        self.log_evidence_cached.unwrap_or(0.0)
     }
 }
 
-/// 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();
+/// `ln` of the observed margin's density under the cavity.
+///
+/// Computed in log space rather than as `pdf(..).ln()`. The density underflows
+/// to zero past about 38 sigma of separation, and clamping that to
+/// `f64::MIN_POSITIVE` reported -708 nats however far out the observation
+/// actually was — 4292 nats adrift at 100 sigma, and unbounded beyond. A score
+/// far from what the model expected is exactly the observation a log-evidence
+/// figure exists to notice.
+fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 {
+    // `hypot`, not `sqrt(a^2 + b^2)`: squaring overflows to infinity above a
+    // sigma of ~1.3e154 and flushes to zero below ~1.5e-154, and `Gaussian`'s
+    // constructors are public so a caller can reach both.
+    let combined_sigma = cavity.sigma().hypot(sigma);
+    let value = ln_pdf(m_obs, cavity.mu(), combined_sigma);
 
-    pdf(m_obs, cavity.mu(), combined_sigma).max(f64::MIN_POSITIVE)
+    // A degenerate cavity (infinite sigma) is the only way to reach a
+    // non-finite result; fall back to the old floor rather than emit -inf.
+    if value.is_finite() {
+        value
+    } else {
+        f64::MIN_POSITIVE.ln()
+    }
 }
 
 #[cfg(test)]
@@ -113,16 +132,16 @@ mod tests {
         let mut vars = VarStore::new();
         let diff = vars.alloc(Gaussian::from_ms(0.0, 6.0));
         let mut f = MarginFactor::new(diff, 5.0, 1.0);
-        assert!(f.evidence_cached.is_none());
+        assert!(f.log_evidence_cached.is_none());
 
         f.propagate(&mut vars);
-        let z = f.evidence_cached.unwrap();
-        // pdf(5, 0, sqrt(37)) ≈ 0.046783
-        assert!((z - 0.04678300292616668).abs() < 1e-10);
+        let z = f.log_evidence_cached.unwrap();
+        // ln pdf(5, 0, sqrt(37)) = ln(0.046783...)
+        assert!((z.exp() - 0.04678300292616668).abs() < 1e-10);
 
         // Subsequent propagations don't change it.
         f.propagate(&mut vars);
-        assert_eq!(f.evidence_cached.unwrap(), z);
+        assert_eq!(f.log_evidence_cached.unwrap(), z);
     }
 
     #[test]
@@ -131,7 +150,7 @@ mod tests {
         let diff = vars.alloc(Gaussian::from_ms(0.0, 6.0));
         let mut f = MarginFactor::new(diff, 5.0, 1.0);
         f.propagate(&mut vars);
-        let logz = f.log_evidence(&vars);
+        let logz = f.log_evidence();
         assert!((logz - (-3.062235327364623)).abs() < 1e-10);
     }
 

src/factor/mod.rs

@@ -20,6 +20,8 @@ pub struct VarStore {
 }
 
 impl VarStore {
+    /// Test-only: inference allocates its store through `ScratchArena`.
+    #[cfg(test)]
     #[must_use]
     pub fn new() -> Self {
         Self::default()
@@ -29,16 +31,13 @@ impl VarStore {
         self.marginals.clear();
     }
 
+    /// Test-only, as `new`.
+    #[cfg(test)]
     #[must_use]
     pub fn len(&self) -> usize {
         self.marginals.len()
     }
 
-    #[must_use]
-    pub fn is_empty(&self) -> bool {
-        self.marginals.is_empty()
-    }
-
     pub fn alloc(&mut self, init: Gaussian) -> VarId {
         let id = VarId(self.marginals.len() as u32);
         self.marginals.push(init);
@@ -55,58 +54,7 @@ impl VarStore {
     }
 }
 
-/// A factor in the EP graph.
-///
-/// Factors hold their own outgoing messages and propagate them by reading
-/// connected variable marginals from a `VarStore` and writing back updated
-/// marginals.
-pub trait Factor: Send + Sync {
-    /// Update outgoing messages and write back to the var store.
-    ///
-    /// Returns the max delta `(|Δmu|, |Δsigma|)` across writes this
-    /// propagation. Used by the `Schedule` to detect convergence.
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64);
-
-    /// Optional log-evidence contribution. Default 0.0 (no contribution).
-    fn log_evidence(&self, _vars: &VarStore) -> f64 {
-        0.0
-    }
-}
-
-/// Enum dispatcher for the built-in factor types.
-///
-/// Using an enum instead of `Box<dyn Factor>` keeps factor data inline and
-/// avoids virtual-call overhead in the hot inference loop.
-#[derive(Debug)]
-pub enum BuiltinFactor {
-    TeamSum(team_sum::TeamSumFactor),
-    RankDiff(rank_diff::RankDiffFactor),
-    Trunc(trunc::TruncFactor),
-    Margin(margin::MarginFactor),
-}
-
-impl Factor for BuiltinFactor {
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
-        match self {
-            Self::TeamSum(f) => f.propagate(vars),
-            Self::RankDiff(f) => f.propagate(vars),
-            Self::Trunc(f) => f.propagate(vars),
-            Self::Margin(f) => f.propagate(vars),
-        }
-    }
-
-    fn log_evidence(&self, vars: &VarStore) -> f64 {
-        match self {
-            Self::Trunc(f) => f.log_evidence(vars),
-            Self::Margin(f) => f.log_evidence(vars),
-            Self::TeamSum(_) | Self::RankDiff(_) => 0.0,
-        }
-    }
-}
-
 pub mod margin;
-pub mod rank_diff;
-pub mod team_sum;
 pub mod trunc;
 
 #[cfg(test)]
@@ -153,20 +101,4 @@ mod tests {
         assert_eq!(store.len(), 0);
         assert_eq!(store.marginals.capacity(), cap);
     }
-
-    #[test]
-    fn builtin_factor_dispatches_to_margin() {
-        use super::margin::MarginFactor;
-        let mut vars = VarStore::new();
-        let diff = vars.alloc(Gaussian::from_ms(0.0, 6.0));
-        let mut f = BuiltinFactor::Margin(MarginFactor::new(diff, 5.0, 1.0));
-
-        f.propagate(&mut vars);
-
-        let result = vars.get(diff);
-        assert!((result.mu() - 4.864864864864865).abs() < 1e-12);
-
-        let logz = f.log_evidence(&vars);
-        assert!((logz - (-3.062235327364623)).abs() < 1e-10);
-    }
 }

src/factor/rank_diff.rs

@@ -1,95 +0,0 @@
-use crate::factor::{Factor, VarId, VarStore};
-
-/// Maintains the constraint `diff = team_a - team_b` between three vars.
-///
-/// On each propagation:
-/// - Reads marginals at `team_a` and `team_b` (which already incorporate any
-///   incoming messages from neighboring factors).
-/// - Computes `new_diff = team_a - team_b` (variance addition; see `Gaussian::Sub`).
-/// - Writes the new marginal to `diff`.
-/// - Returns the delta against the previous diff value.
-///
-/// This factor does NOT store an outgoing message; the diff variable is
-/// effectively replaced on each propagation. The `TruncFactor` on the same diff
-/// var holds the EP-divide message that produces the cavity.
-#[derive(Debug)]
-pub struct RankDiffFactor {
-    pub team_a: VarId,
-    pub team_b: VarId,
-    pub diff: VarId,
-}
-
-impl Factor for RankDiffFactor {
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
-        let a = vars.get(self.team_a);
-        let b = vars.get(self.team_b);
-        let new_diff = a - b;
-        let old = vars.get(self.diff);
-        vars.set(self.diff, new_diff);
-        old.delta(new_diff)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::{N_INF, gaussian::Gaussian};
-
-    #[test]
-    fn diff_of_two_known_gaussians() {
-        let mut vars = VarStore::new();
-        let team_a = vars.alloc(Gaussian::from_ms(25.0, 3.0));
-        let team_b = vars.alloc(Gaussian::from_ms(20.0, 4.0));
-        let diff = vars.alloc(N_INF);
-
-        let mut f = RankDiffFactor {
-            team_a,
-            team_b,
-            diff,
-        };
-        f.propagate(&mut vars);
-
-        let result = vars.get(diff);
-        // mu = 25 - 20 = 5; var = 9 + 16 = 25; sigma = 5
-        assert!((result.mu() - 5.0).abs() < 1e-12);
-        assert!((result.sigma() - 5.0).abs() < 1e-12);
-    }
-
-    #[test]
-    fn delta_zero_on_repeat() {
-        let mut vars = VarStore::new();
-        let team_a = vars.alloc(Gaussian::from_ms(10.0, 2.0));
-        let team_b = vars.alloc(Gaussian::from_ms(8.0, 1.0));
-        let diff = vars.alloc(N_INF);
-
-        let mut f = RankDiffFactor {
-            team_a,
-            team_b,
-            diff,
-        };
-        f.propagate(&mut vars);
-        let (dmu, dsig) = f.propagate(&mut vars);
-        assert!(dmu < 1e-12);
-        assert!(dsig < 1e-12);
-    }
-
-    #[test]
-    fn delta_reflects_team_change() {
-        let mut vars = VarStore::new();
-        let team_a = vars.alloc(Gaussian::from_ms(10.0, 1.0));
-        let team_b = vars.alloc(Gaussian::from_ms(0.0, 1.0));
-        let diff = vars.alloc(N_INF);
-
-        let mut f = RankDiffFactor {
-            team_a,
-            team_b,
-            diff,
-        };
-        f.propagate(&mut vars);
-
-        // change team_a, repropagate; delta should be positive
-        vars.set(team_a, Gaussian::from_ms(15.0, 1.0));
-        let (dmu, _dsig) = f.propagate(&mut vars);
-        assert!(dmu > 4.0, "expected ~5 delta, got {}", dmu);
-    }
-}

src/factor/team_sum.rs

@@ -1,98 +0,0 @@
-use crate::{
-    N00,
-    factor::{Factor, VarId, VarStore},
-    gaussian::Gaussian,
-};
-
-/// Computes the weighted sum of player performances into a team-perf var.
-///
-/// Inputs are pre-computed player performance Gaussians (i.e., rating priors
-/// already with beta² noise added via `Rating::performance()`). The factor
-/// runs once per game and writes the weighted sum to the output var.
-#[derive(Debug)]
-pub struct TeamSumFactor {
-    pub inputs: Vec<(Gaussian, f64)>,
-    pub out: VarId,
-}
-
-impl Factor for TeamSumFactor {
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
-        let perf = self.inputs.iter().fold(N00, |acc, (g, w)| acc + (*g * *w));
-        let old = vars.get(self.out);
-        vars.set(self.out, perf);
-        old.delta(perf)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::N_INF;
-
-    #[test]
-    fn single_player_unit_weight() {
-        let mut vars = VarStore::new();
-        let out = vars.alloc(N_INF);
-        let g = Gaussian::from_ms(25.0, 5.0);
-        let mut f = TeamSumFactor {
-            inputs: vec![(g, 1.0)],
-            out,
-        };
-
-        f.propagate(&mut vars);
-        let result = vars.get(out);
-        assert!((result.mu() - 25.0).abs() < 1e-12);
-        assert!((result.sigma() - 5.0).abs() < 1e-12);
-    }
-
-    #[test]
-    fn two_players_summed() {
-        let mut vars = VarStore::new();
-        let out = vars.alloc(N_INF);
-        let g1 = Gaussian::from_ms(20.0, 3.0);
-        let g2 = Gaussian::from_ms(30.0, 4.0);
-        let mut f = TeamSumFactor {
-            inputs: vec![(g1, 1.0), (g2, 1.0)],
-            out,
-        };
-
-        f.propagate(&mut vars);
-        let result = vars.get(out);
-        // sum: mu = 20 + 30 = 50, var = 9 + 16 = 25, sigma = 5
-        assert!((result.mu() - 50.0).abs() < 1e-12);
-        assert!((result.sigma() - 5.0).abs() < 1e-12);
-    }
-
-    #[test]
-    fn weighted_inputs() {
-        let mut vars = VarStore::new();
-        let out = vars.alloc(N_INF);
-        let g = Gaussian::from_ms(10.0, 2.0);
-        let mut f = TeamSumFactor {
-            inputs: vec![(g, 2.0)],
-            out,
-        };
-
-        f.propagate(&mut vars);
-        let result = vars.get(out);
-        // g * 2.0: mu = 10*2 = 20, sigma = 2*2 = 4
-        assert!((result.mu() - 20.0).abs() < 1e-12);
-        assert!((result.sigma() - 4.0).abs() < 1e-12);
-    }
-
-    #[test]
-    fn delta_is_zero_on_repeat_propagate() {
-        let mut vars = VarStore::new();
-        let out = vars.alloc(N_INF);
-        let g = Gaussian::from_ms(5.0, 1.0);
-        let mut f = TeamSumFactor {
-            inputs: vec![(g, 1.0)],
-            out,
-        };
-
-        f.propagate(&mut vars);
-        let (dmu, dsig) = f.propagate(&mut vars);
-        assert!(dmu < 1e-12, "expected ~0 delta on repeat, got {}", dmu);
-        assert!(dsig < 1e-12);
-    }
-}

src/factor/trunc.rs

@@ -1,7 +1,8 @@
 use crate::{
-    N_INF, approx, cdf,
-    factor::{Factor, VarId, VarStore},
+    N_INF, approx,
+    factor::{VarId, VarStore},
     gaussian::Gaussian,
+    ln_interval, ln_sf,
 };
 
 /// EP truncation factor on a diff variable.
@@ -18,7 +19,7 @@ pub struct TruncFactor {
     /// Outgoing message to the diff variable (initial: `N_INF`, the EP identity).
     pub(crate) msg: Gaussian,
     /// Cached evidence (linear, not log) computed from the cavity on first propagation.
-    pub(crate) evidence_cached: Option<f64>,
+    pub(crate) log_evidence_cached: Option<f64>,
 }
 
 impl TruncFactor {
@@ -29,7 +30,7 @@ impl TruncFactor {
             margin,
             tie,
             msg: N_INF,
-            evidence_cached: None,
+            log_evidence_cached: None,
         }
     }
 }
@@ -42,8 +43,8 @@ impl TruncFactor {
         let marginal = vars.get(self.diff);
         let cavity = marginal / self.msg;
 
-        if self.evidence_cached.is_none() {
-            self.evidence_cached = Some(cavity_evidence(cavity, self.margin, self.tie));
+        if self.log_evidence_cached.is_none() {
+            self.log_evidence_cached = Some(cavity_log_evidence(cavity, self.margin, self.tie));
         }
 
         let trunc = approx(cavity, self.margin, self.tie);
@@ -62,31 +63,40 @@ impl TruncFactor {
     }
 }
 
-impl Factor for TruncFactor {
-    fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
+/// Undamped wrappers, used by this module's tests. Inference drives these
+/// factors through `propagate_with_alpha` and reads the cached log evidence
+/// directly, so these are not on any production path.
+#[cfg(test)]
+impl TruncFactor {
+    pub(crate) fn propagate(&mut self, vars: &mut VarStore) -> (f64, f64) {
         self.propagate_with_alpha(vars, 1.0)
     }
-
-    fn log_evidence(&self, _vars: &VarStore) -> f64 {
-        self.evidence_cached.unwrap_or(1.0).ln()
-    }
 }
 
-/// P(diff > margin) for non-tie, P(|diff| < margin) for tie.
+/// `ln P(diff > margin)` for a win, `ln P(|diff| < margin)` for a 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 {
-    let raw = if tie {
-        cdf(margin, diff.mu(), diff.sigma()) - cdf(-margin, diff.mu(), diff.sigma())
+/// Computed in log space throughout. Two earlier shapes both lost the tail:
+/// `1 - cdf(..)` cancelled away every digit of an unlikely outcome, and even
+/// once that was fixed the linear probability underflows to zero past about 38
+/// sigma, where clamping reported -708 nats regardless of the truth. An upset
+/// is the observation a log-evidence figure exists to notice, so it has to stay
+/// exact precisely where it is smallest.
+fn cavity_log_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 {
+    let (mu, sigma) = (diff.mu(), diff.sigma());
+
+    let value = if tie {
+        ln_interval(-margin, margin, mu, sigma)
     } else {
-        1.0 - cdf(margin, diff.mu(), diff.sigma())
+        ln_sf(margin, mu, sigma)
     };
 
-    raw.clamp(f64::MIN_POSITIVE, 1.0)
+    // A degenerate cavity is the only route to a non-finite result; keep the
+    // old floor for it rather than letting -inf poison the whole history's sum.
+    if value.is_finite() {
+        value
+    } else {
+        f64::MIN_POSITIVE.ln()
+    }
 }
 
 #[cfg(test)]
@@ -117,19 +127,90 @@ mod tests {
         let diff = vars.alloc(Gaussian::from_ms(2.0, 3.0));
 
         let mut f = TruncFactor::new(diff, 0.0, false);
-        assert!(f.evidence_cached.is_none());
+        assert!(f.log_evidence_cached.is_none());
 
         f.propagate(&mut vars);
-        assert!(f.evidence_cached.is_some());
-        let first = f.evidence_cached.unwrap();
+        assert!(f.log_evidence_cached.is_some());
+        let first = f.log_evidence_cached.unwrap();
 
         // Evidence should be P(diff > 0) for diff ~ N(2, 9) ≈ 0.748
-        assert!(first > 0.7);
-        assert!(first < 0.8);
+        assert!(first.exp() > 0.7);
+        assert!(first.exp() < 0.8);
 
         // Subsequent propagations don't change it.
         f.propagate(&mut vars);
-        assert_eq!(f.evidence_cached.unwrap(), first);
+        assert_eq!(f.log_evidence_cached.unwrap(), first);
+    }
+
+    /// The defect this guards: `1 - cdf` collapsed to zero for a surprising
+    /// result, the clamp turned that into `f64::MIN_POSITIVE`, and
+    /// `log_evidence` reported ln of *that* — about -708 whatever the truth
+    /// was. An upset is the observation a model-comparison score exists to
+    /// notice, so it was wrong exactly where it mattered.
+    #[test]
+    fn evidence_of_an_upset_is_not_flattened_to_the_clamp_floor() {
+        // diff ~ N(-9, 1) with margin 0: the favoured side lost by nine sigma.
+        let evidence = cavity_log_evidence(Gaussian::from_ms(-9.0, 1.0), 0.0, false).exp();
+
+        assert!(
+            evidence > f64::MIN_POSITIVE,
+            "evidence collapsed onto the clamp floor: {evidence}"
+        );
+        // P(X > 0) for X ~ N(-9, 1) is the standard normal tail at 9 sigma.
+        assert!(
+            (evidence - 1.128_588e-19).abs() / 1.128_588e-19 < 1e-6,
+            "expected ~1.13e-19, got {evidence}"
+        );
+        assert!(
+            (evidence.ln() + 43.628).abs() < 1e-2,
+            "log evidence {} should be about -43.6, not -708",
+            evidence.ln()
+        );
+    }
+
+    /// Evidence must stay finite and positive however extreme the mismatch,
+    /// since `log_evidence` sums across the whole history and one `-inf` or
+    /// `NaN` poisons all of it.
+    ///
+    /// Finiteness alone is too weak a bar — the clamped version was finite too,
+    /// and wrong by hundreds of nats. `log_evidence_tracks_the_analytic_tail`
+    /// below is the assertion that actually holds this up.
+    #[test]
+    fn evidence_stays_positive_and_finite_at_any_separation() {
+        for mu in [-300.0f64, -50.0, -9.0, 0.0, 9.0, 50.0, 300.0] {
+            for tie in [false, true] {
+                let ln_e = cavity_log_evidence(Gaussian::from_ms(mu, 1.0), 1.0, tie);
+                assert!(
+                    ln_e.is_finite() && ln_e <= 0.0,
+                    "mu={mu} tie={tie}: log evidence {ln_e} is not a log-probability"
+                );
+            }
+        }
+    }
+
+    /// The clamp used to floor everything past ~38 sigma at `ln(MIN_POSITIVE)`
+    /// = -708, however far out the real observation was. In log space the
+    /// answer is a polynomial and stays exact: at 1000 sigma the truth is about
+    /// -500_000 nats, and -708 is not a rounding error.
+    #[test]
+    fn log_evidence_tracks_the_analytic_tail() {
+        for mu in [-40.0f64, -60.0, -100.0, -1000.0] {
+            // P(diff > 0) for diff ~ N(mu, 1), mu far below zero.
+            let got = cavity_log_evidence(Gaussian::from_ms(mu, 1.0), 0.0, false);
+
+            // ln Phi(mu) ~ -mu^2/2 - ln(-mu) - ln(sqrt(2 pi)) for mu << 0.
+            let z = -mu;
+            let approx = -0.5 * z * z - z.ln() - (2.0 * std::f64::consts::PI).sqrt().ln();
+
+            assert!(
+                got < f64::MIN_POSITIVE.ln(),
+                "mu={mu}: {got} is still stuck on the old clamp floor"
+            );
+            assert!(
+                (got - approx).abs() / approx.abs() < 1e-3,
+                "mu={mu}: got {got}, asymptotic expectation {approx}"
+            );
+        }
     }
 
     #[test]
@@ -141,7 +222,7 @@ mod tests {
         f.propagate(&mut vars);
 
         // For diff ~ N(0, 4), tie=true with margin=1: P(-1 < diff < 1) ≈ 0.383
-        let ev = f.evidence_cached.unwrap();
+        let ev = f.log_evidence_cached.unwrap().exp();
         assert!(ev > 0.35 && ev < 0.42);
     }
 

src/factors.rs

@@ -1,13 +0,0 @@
-//! Factor-graph public API.
-//!
-//! Power users can construct custom factor graphs via `Game::custom` (T2
-//! minimal; full ergonomics in T4) and drive them with custom `Schedule`
-//! implementations.
-
-pub use crate::{
-    factor::{
-        BuiltinFactor, Factor, VarId, VarStore, margin::MarginFactor, rank_diff::RankDiffFactor,
-        team_sum::TeamSumFactor, trunc::TruncFactor,
-    },
-    schedule::{EpsilonOrMax, Schedule, ScheduleReport},
-};

src/game.rs

@@ -46,8 +46,8 @@ impl DiffFactor {
     /// reaches.
     pub(crate) fn log_evidence(&self) -> f64 {
         match self {
-            Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0).ln(),
-            Self::Margin(f) => f.evidence_cached.unwrap_or(1.0).ln(),
+            Self::Trunc(f) => f.log_evidence_cached.unwrap_or(0.0),
+            Self::Margin(f) => f.log_evidence_cached.unwrap_or(0.0),
         }
     }
 
@@ -433,6 +433,9 @@ impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
 impl<T: Time, D: Drift<T>> Game<'_, T, D> {
     /// # Errors
     ///
+    /// - `InvalidParameter` if `options.convergence` is out of range — an
+    ///   `alpha` of zero would leave every EP update unapplied and silently
+    ///   return the priors.
     /// - `InvalidProbability` if `options.p_draw` is outside `[0.0, 1.0)`.
     /// - `MismatchedShape` if the outcome's rank count differs from `teams.len()`.
     /// - `WrongOutcomeKind` if `outcome` is not `Outcome::Ranked`.
@@ -444,6 +447,7 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
         outcome: crate::Outcome,
         options: &GameOptions,
     ) -> Result<OwnedGame<T, D>, crate::InferenceError> {
+        options.convergence.validate()?;
         if !(0.0..1.0).contains(&options.p_draw) {
             return Err(crate::InferenceError::InvalidProbability {
                 value: options.p_draw,
@@ -491,8 +495,8 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
 
     /// # Errors
     ///
-    /// - `InvalidParameter` if `options.score_sigma` is not strictly positive,
-    ///   or is NaN.
+    /// - `InvalidParameter` if `options.score_sigma` is not strictly positive
+    ///   or is NaN, or if `options.convergence` is out of range.
     /// - `MismatchedShape` if the outcome's score count differs from `teams.len()`.
     /// - `WrongOutcomeKind` if `outcome` is not `Outcome::Scored`.
     pub fn scored(
@@ -500,6 +504,7 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
         outcome: crate::Outcome,
         options: &GameOptions,
     ) -> Result<OwnedGame<T, D>, crate::InferenceError> {
+        options.convergence.validate()?;
         if options.score_sigma <= 0.0 || options.score_sigma.is_nan() {
             return Err(crate::InferenceError::InvalidParameter {
                 name: "score_sigma",
@@ -532,18 +537,20 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
         ))
     }
 
+    /// Convenience wrapper over [`Game::ranked`] for two single-player teams.
+    ///
     /// # Errors
     ///
-    /// Delegates to [`Game::ranked`] with default options, so it returns the
-    /// same errors — in practice `WrongOutcomeKind` for a non-ranked outcome,
-    /// or `TieWithoutDrawProbability` for a draw, since the default `p_draw`
-    /// applies rather than one you chose.
+    /// Delegates to [`Game::ranked`], so it returns the same errors — in
+    /// practice `WrongOutcomeKind` for a non-ranked outcome, or
+    /// `TieWithoutDrawProbability` for a draw when `options.p_draw` is zero.
     pub fn one_v_one(
         a: &Rating<T, D>,
         b: &Rating<T, D>,
         outcome: crate::Outcome,
+        options: &GameOptions,
     ) -> Result<(Gaussian, Gaussian), crate::InferenceError> {
-        let game = Self::ranked(&[&[*a], &[*b]], outcome, &GameOptions::default())?;
+        let game = Self::ranked(&[&[*a], &[*b]], outcome, options)?;
         let post = game.posteriors();
         Ok((post[0][0], post[1][0]))
     }
@@ -561,15 +568,6 @@ impl<T: Time, D: Drift<T>> Game<'_, T, D> {
         let team_refs: Vec<&[Rating<T, D>]> = teams.iter().map(|t| t.as_slice()).collect();
         Self::ranked(&team_refs, outcome, options)
     }
-
-    #[doc(hidden)]
-    pub fn custom<S: crate::factors::Schedule>(
-        factors: &mut [crate::factors::BuiltinFactor],
-        vars: &mut crate::factors::VarStore,
-        schedule: &S,
-    ) -> crate::factors::ScheduleReport {
-        schedule.run(factors, vars)
-    }
 }
 
 #[cfg(test)]
@@ -726,9 +724,15 @@ mod tests {
         let c = p[2][0];
 
         // T1 ULP shift: mu rounds to 25.0 (was 24.999999) under natural-parameter storage.
+        //
+        // The 1e-6-place values moved when `erfc_inv`'s sign error was fixed:
+        // this case runs at `p_draw = 0.5`, so it goes through `compute_margin`,
+        // and the margin is now 8.4e-8 from the exact quantile where it was
+        // 1.46e-7. Verified as movement *toward* analytic truth, not a
+        // regression — see `erfc_inv_matches_known_quantiles`.
         assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.092561), epsilon = 1e-6);
-        assert_ulps_eq!(b, Gaussian::from_ms(33.379314, 6.483575), epsilon = 1e-6);
-        assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483575), epsilon = 1e-6);
+        assert_ulps_eq!(b, Gaussian::from_ms(33.379315, 6.483576), epsilon = 1e-6);
+        assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483576), epsilon = 1e-6);
     }
 
     #[test]
@@ -1241,7 +1245,7 @@ mod tests {
         );
         assert_ulps_eq!(
             p[1][0],
-            Gaussian::from_ms(19.287197, 7.243465),
+            Gaussian::from_ms(19.287198285, 7.243465848),
             epsilon = 1e-6
         );
         assert_ulps_eq!(
@@ -1301,7 +1305,7 @@ mod tests {
 
         assert_ulps_eq!(
             p[0][0],
-            Gaussian::from_ms(31.674697, 7.501180),
+            Gaussian::from_ms(31.674698083, 7.501180037),
             epsilon = 1e-6
         );
         assert_ulps_eq!(

src/gaussian.rs

@@ -145,6 +145,45 @@ impl Gaussian {
         Self::from_mv(self.mu(), self.variance() + variance_delta)
     }
 
+    /// `P(X < x)` under this Gaussian.
+    ///
+    /// The question a stopping rule asks: *how sure am I that this competitor's
+    /// true skill is below the cutoff?* Expressing that as a probability keeps
+    /// its meaning as sigma changes, where a `mu + z * sigma` band silently
+    /// means different confidence at different uncertainties — which is exactly
+    /// the regime a stopping rule operates in.
+    ///
+    /// Accurate in the *lower* tail. For the upper tail use
+    /// [`Gaussian::probability_above`] rather than `1.0 - probability_below(x)`,
+    /// which cancels away every significant digit once the result is small.
+    ///
+    /// An improper Gaussian (non-positive precision) has no defined mean, so
+    /// this returns `0.5` — the same convention `mu()` and `sigma()` follow.
+    #[must_use]
+    pub fn probability_below(&self, x: f64) -> f64 {
+        if self.pi <= 0.0 {
+            return 0.5;
+        }
+        crate::cdf(x, self.mu(), self.sigma())
+    }
+
+    /// `P(X > x)` under this Gaussian.
+    ///
+    /// Computed as a survival function rather than `1 - cdf`, so it keeps full
+    /// relative precision in the upper tail: `1 - cdf` returns exactly zero
+    /// past about 8.3 sigma, where the true value is still 1e-19 and perfectly
+    /// representable. A stopping rule is evaluated precisely there — the
+    /// interesting cases are the ones near certainty.
+    ///
+    /// An improper Gaussian returns `0.5`, as [`Gaussian::probability_below`].
+    #[must_use]
+    pub fn probability_above(&self, x: f64) -> f64 {
+        if self.pi <= 0.0 {
+            return 0.5;
+        }
+        crate::sf(x, self.mu(), self.sigma())
+    }
+
     /// EP damping in natural-parameter space: `α·new + (1−α)·self`.
     ///
     /// Used by within-game inference to stabilise oscillating fixed-point
@@ -340,3 +379,68 @@ mod tests {
         assert!((damped.tau() - expected_tau).abs() < 1e-12);
     }
 }
+
+#[cfg(test)]
+mod tail_probability_tests {
+    use super::*;
+
+    #[test]
+    fn probability_below_matches_published_quantiles() {
+        let g = Gaussian::from_ms(0.0, 1.0);
+        for (x, expected) in [
+            (-1.959_963_984_540_054, 0.025),
+            (0.0, 0.5),
+            (1.281_551_565_544_6, 0.9),
+            (1.959_963_984_540_054, 0.975),
+        ] {
+            let got = g.probability_below(x);
+            assert!(
+                (got - expected).abs() < 1e-12,
+                "P(X < {x}) = {got}, expected {expected}"
+            );
+        }
+    }
+
+    #[test]
+    fn the_two_tails_partition_the_mass() {
+        let g = Gaussian::from_ms(3.0, 2.0);
+        for x in [-4.0f64, 0.0, 3.0, 7.5] {
+            let total = g.probability_below(x) + g.probability_above(x);
+            assert!((total - 1.0).abs() < 1e-15, "at {x}: {total}");
+        }
+    }
+
+    /// The reason `probability_above` exists rather than `1 - probability_below`.
+    #[test]
+    fn probability_above_keeps_precision_where_the_complement_collapses() {
+        let g = Gaussian::from_ms(0.0, 1.0);
+        for (x, expected) in [(9.0f64, 1.128_588e-19), (20.0, 2.753_624e-89)] {
+            let got = g.probability_above(x);
+            assert!(
+                (got - expected).abs() / expected < 1e-6,
+                "P(X > {x}) = {got}, expected ~{expected}"
+            );
+            assert_eq!(
+                1.0 - g.probability_below(x),
+                0.0,
+                "the complement should still collapse at {x}"
+            );
+        }
+    }
+
+    #[test]
+    fn a_scaled_gaussian_shifts_and_stretches() {
+        let g = Gaussian::from_ms(25.0, 6.0);
+        assert!((g.probability_below(25.0) - 0.5).abs() < 1e-15);
+        // One sigma either side of the mean.
+        assert!((g.probability_below(31.0) - 0.841_344_746_068_543).abs() < 1e-12);
+        assert!((g.probability_above(19.0) - 0.841_344_746_068_543).abs() < 1e-12);
+    }
+
+    #[test]
+    fn an_improper_gaussian_is_uninformative_rather_than_nan() {
+        let improper = Gaussian::from_ms(0.0, f64::INFINITY);
+        assert_eq!(improper.probability_below(5.0), 0.5);
+        assert_eq!(improper.probability_above(5.0), 0.5);
+    }
+}

src/joint.rs

@@ -0,0 +1,99 @@
+//! Posterior of a linear combination of competitors.
+//!
+//! Every accessor on `History` returns a per-competitor marginal, and almost
+//! nothing a consumer publishes is one competitor: "can we tell these two
+//! apart" is a difference, "what was this round worth" is a sum. Combining
+//! marginals means assuming the competitors are independent, and they are
+//! correlated through every event they share — which is the mechanism the model
+//! exists to exploit.
+//!
+//! Measured on a five-competitor round robin, the exact correlation is +0.857,
+//! so `sqrt(sa^2 + sb^2)` overstates the width of a difference by 2.6x.
+
+/// Solve `A z = b` for a symmetric positive-definite `A`, by Cholesky.
+///
+/// `a` is row-major and is consumed as scratch.
+///
+/// Returns `None` if the matrix is not positive-definite, which for a precision
+/// matrix means the model is improper — a competitor with no prior and no
+/// evidence.
+pub(crate) fn solve_spd(mut a: Vec<f64>, b: &[f64]) -> Option<Vec<f64>> {
+    let n = b.len();
+    debug_assert_eq!(a.len(), n * n);
+
+    // In-place Cholesky: A = L L^T, lower triangle.
+    for j in 0..n {
+        let mut d = a[j * n + j];
+        for k in 0..j {
+            d -= a[j * n + k] * a[j * n + k];
+        }
+        // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here too,
+        // and a negated comparison would let it through as "not positive".
+        if d.is_nan() || d <= 0.0 {
+            return None;
+        }
+        let d = d.sqrt();
+        a[j * n + j] = d;
+
+        for i in j + 1..n {
+            let mut s = a[i * n + j];
+            for k in 0..j {
+                s -= a[i * n + k] * a[j * n + k];
+            }
+            a[i * n + j] = s / d;
+        }
+    }
+
+    // Forward substitution, then back substitution.
+    let mut z = b.to_vec();
+    for i in 0..n {
+        let mut s = z[i];
+        for k in 0..i {
+            s -= a[i * n + k] * z[k];
+        }
+        z[i] = s / a[i * n + i];
+    }
+    for i in (0..n).rev() {
+        let mut s = z[i];
+        for k in i + 1..n {
+            s -= a[k * n + i] * z[k];
+        }
+        z[i] = s / a[i * n + i];
+    }
+
+    Some(z)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn solves_a_known_system() {
+        // [[4, 1], [1, 3]] z = [1, 2]  =>  z = [1/11, 7/11]
+        let a = vec![4.0, 1.0, 1.0, 3.0];
+        let z = solve_spd(a, &[1.0, 2.0]).unwrap();
+        assert!((z[0] - 1.0 / 11.0).abs() < 1e-12, "{z:?}");
+        assert!((z[1] - 7.0 / 11.0).abs() < 1e-12, "{z:?}");
+    }
+
+    #[test]
+    fn recovers_the_inverse_diagonal() {
+        // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is
+        // [0.75, 1.0, 0.75].
+        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() {
+            let mut e = vec![0.0; 3];
+            e[i] = 1.0;
+            let z = solve_spd(a.clone(), &e).unwrap();
+            assert!((z[i] - expected).abs() < 1e-12, "row {i}: {z:?}");
+        }
+    }
+
+    #[test]
+    fn rejects_a_non_positive_definite_matrix() {
+        // Singular: the second row is a multiple of the first.
+        let a = vec![1.0, 2.0, 2.0, 4.0];
+        assert!(solve_spd(a, &[1.0, 1.0]).is_none());
+    }
+}

src/lib.rs

@@ -86,6 +86,19 @@
 
 #![forbid(unsafe_code)]
 
+/// Compiles every `rust` block in `README.md` as a doctest.
+///
+/// The README is not the crate's front page — the module docs above are — so it
+/// is pulled in here rather than via a crate-level `#![doc = ...]`, purely so
+/// its examples are type-checked. Without this nothing compiled them, and they
+/// had drifted far enough that four blocks no longer built (#35). `cfg(doctest)`
+/// means this type exists only while collecting doctests.
+///
+/// Blocks that are illustrative rather than runnable are fenced as `text`.
+#[cfg(doctest)]
+#[doc = include_str!("../README.md")]
+pub struct ReadmeDoctests;
+
 use std::{
     cmp::Reverse,
     f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2},
@@ -97,6 +110,7 @@ pub(crate) mod arena;
 mod time;
 mod time_slice;
 pub use time_slice::{EventKind, TimeSlice};
+mod acquisition;
 mod color_group;
 mod competitor;
 mod convergence;
@@ -105,22 +119,24 @@ mod error;
 mod event;
 mod event_builder;
 pub(crate) mod factor;
-pub mod factors;
 mod game;
 pub mod gaussian;
 mod history;
+mod joint;
 mod key_table;
 mod matrix;
 mod observer;
 mod outcome;
+mod predict;
+pub(crate) mod quadrature;
 mod rating;
-pub(crate) mod schedule;
 pub mod storage;
 
+pub use acquisition::expected_information_gain;
 pub use competitor::Competitor;
 pub use convergence::{ConvergenceOptions, ConvergenceReport};
 pub use drift::{ConstantDrift, Drift};
-pub use error::InferenceError;
+pub use error::{InferenceError, UnknownKeys};
 pub use event::{Event, Member, Team};
 pub use event_builder::EventBuilder;
 pub use game::{Game, GameOptions, OwnedGame};
@@ -130,8 +146,8 @@ pub use key_table::KeyTable;
 use matrix::Matrix;
 pub use observer::{NullObserver, Observer};
 pub use outcome::Outcome;
+pub use predict::Prediction;
 pub use rating::Rating;
-pub use schedule::ScheduleReport;
 pub use time::{Time, Untimed};
 
 pub const BETA: f64 = 1.0;
@@ -140,9 +156,48 @@ pub const SIGMA: f64 = BETA * 6.0;
 pub const GAMMA: f64 = BETA * 0.03;
 pub const P_DRAW: f64 = 0.0;
 pub const EPSILON: f64 = 1e-6;
+/// Default cap on convergence sweeps.
+///
+/// **This is a floor, not a recommendation.** It is adequate for small
+/// histories and is quickly outgrown: a history of 400 events over 100
+/// competitors already stops here with a final step of ~7e-3 against the 1e-6
+/// default tolerance — four orders of magnitude short — and a dense joint model
+/// of ~2,000 nodes over ~3,300 events has been measured needing 76 to 161.
+///
+/// Overrunning it is not an error, and deliberately so: `converge` returns a
+/// [`ConvergenceReport`] whose `converged` flag says what happened. But a fit
+/// that stopped short is *wrong by a little*, which is the worst available
+/// failure — every rating is finite and ordered sensibly, and nothing in the
+/// numbers themselves says they were still moving. Read the report; the type is
+/// `#[must_use]` for that reason.
+///
+/// Raise it via [`ConvergenceOptions`]. Convergence cost is roughly linear in
+/// the cap, and for anything but a toy the extra sweeps are milliseconds.
 pub const ITERATIONS: usize = 30;
 
+/// Largest team count `History::predict_outcome` will enumerate.
+///
+/// The outcome space holds `n! * 2^(n-1)` events, so it grows factorially:
+/// 1_920 at five teams, 23_040 at six, 322_560 at seven. Six is where
+/// enumerating on a caller's behalf stops being reasonable.
+pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION;
+
 const SQRT_TAU: f64 = 2.5066282746310002;
+/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction.
+const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3;
+/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form.
+const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;
+/// How many window widths into the tail before a tie window is treated as a
+/// half-line. Beyond this the truncated mass is concentrated within `1/alpha`
+/// of the near edge, so the far edge contributes nothing measurable.
+const HALF_LINE_WINDOW: f64 = 10.0;
+/// Where `v - alpha` switches from subtraction to its asymptotic series.
+///
+/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the
+/// four-term series is good to ~1e-10 by here, so the two are at their closest
+/// agreement around this point. Below it the subtraction is exact; above it the
+/// series is.
+const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;
 
 pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
 pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
@@ -174,24 +229,47 @@ impl From<Index> for usize {
     }
 }
 
+/// Complementary error function.
+///
+/// # Why every transcendental in this crate goes through `libm`
+///
+/// IEEE 754 specifies the basic operations and `sqrt` exactly, but says nothing
+/// about `exp`, `log` or `erf`. `std`'s versions delegate to the *system* math
+/// library, so they differ between platforms: measured here, `f64::exp` and
+/// `libm::exp` disagree on 9.7% of inputs and `f64::ln` / `libm::log` on 5.0%,
+/// each by one ULP.
+///
+/// Inference is an iterative fixed point, so a one-ULP difference can change an
+/// iteration count and therefore the answer by more than one ULP. Routing every
+/// transcendental through `libm` makes a fit reproducible across platforms, not
+/// just across thread counts as `tests/determinism.rs` already checks.
+///
+/// **So: use `libm::exp` / `libm::log` in inference code, never `f64::exp` /
+/// `f64::ln`.** `sqrt` is exempt — IEEE specifies it exactly, so `f64::sqrt` is
+/// already portable. Test code may use whichever is clearer.
+///
+/// It costs nothing: `Batch::iteration` measured -2.7% [-5.7%, -0.3%] with the
+/// whole set swapped.
+///
+/// Delegates to `libm`, which is the Rust port of FDLIBM and accurate to about
+/// one ULP. This replaced a Numerical Recipes `erfcc` rational approximation
+/// whose documented bound was 1.2e-7 *relative* — measured at ~1e-7 across the
+/// whole range, and the binding accuracy constraint on the entire crate.
+///
+/// The swap is free. 98% of the arguments inference passes here have
+/// `|x| < 0.84375`, which is exactly where FDLIBM skips the exponential
+/// entirely, so the longer polynomial costs nothing on the distribution that
+/// actually occurs: `Batch::iteration` moved -1.6% [-4.7%, +0.9%], p = 0.31.
+///
+/// What it bought: `compute_margin` went from 8.4e-8 to 1.7e-16 against exact
+/// quantiles, `cdf(mu, mu, sigma)` is now exactly 0.5, and `sf + cdf` sums to
+/// one within a single ULP where it was 3e-8 out.
 fn erfc(x: f64) -> f64 {
-    let z = x.abs();
-    let t = 1.0 / (1.0 + z / 2.0);
-
-    let a = -0.82215223 + t * 0.17087277;
-    let b = 1.48851587 + t * a;
-    let c = -1.13520398 + t * b;
-    let d = 0.27886807 + t * c;
-    let e = -0.18628806 + t * d;
-    let f = 0.09678418 + t * e;
-    let g = 0.37409196 + t * f;
-    let h = 1.00002368 + t * g;
-
-    let r = t * (-z * z - 1.26551223 + t * h).exp();
-
-    if x >= 0.0 { r } else { 2.0 - r }
+    libm::erfc(x)
 }
 
+/// The previous Numerical Recipes `erfcc`, kept only so the timing test can
+/// compare both in one binary. Removed once the comparison is recorded.
 fn erfc_inv(mut y: f64) -> f64 {
     if y >= 2.0 {
         return f64::NEG_INFINITY;
@@ -207,14 +285,22 @@ fn erfc_inv(mut y: f64) -> f64 {
         y = 2.0 - y;
     }
 
-    let t = (-2.0 * (y / 2.0).ln()).sqrt();
+    let t = libm::sqrt(-2.0 * libm::log(y / 2.0));
 
-    let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
+    // The leading coefficient is NEGATIVE. `rational - t` is negative here, so
+    // a positive coefficient mirrors the starting point to `-x0` — the
+    // reflection of the root. Newton then has to cross the origin to get back,
+    // which a fixed iteration count does not manage: measured against the true
+    // value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error
+    // grew as y shrank until `compute_margin` stopped being monotone in
+    // `p_draw` altogether.
+    let mut x =
+        -FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);
 
     for _ in 0..3 {
         let err = erfc(x) - y;
 
-        x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err)
+        x += err / (FRAC_2_SQRT_PI * libm::exp(-(x * x)) - x * err)
     }
 
     if y < 1.0 { x } else { -x }
@@ -234,32 +320,216 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
     0.5 * erfc(z)
 }
 
+/// `P(X > x)` for `X ~ N(mu, sigma^2)`.
+///
+/// The survival function, computed directly rather than as `1 - cdf(..)`.
+///
+/// The two are algebraically identical and numerically are not. `cdf` returns
+/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels
+/// away every significant digit the tail had: measured against this function,
+/// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly
+/// zero beyond about 8.3 sigma — where the true value is still 1e-19 and
+/// perfectly representable. `erfc` holds *relative* accuracy all the way down
+/// to 1e-296, so the precision is there to keep; only the subtraction threw it
+/// away.
+///
+/// This matters most where evidence is smallest, which is exactly where an
+/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of
+/// whether the truth was -43 or -600.
+pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 {
+    0.5 * erfc((x - mu) / (sigma * SQRT_2))
+}
+
+/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`.
+///
+/// Exists so the exponential factor common to a Gaussian density and its tail
+/// integral can be cancelled *analytically* instead of being computed twice
+/// and divided. Both underflow to zero past about 26 sigma, and their ratio is
+/// then `0/0` — finite in the limit, `NaN` in floating point.
+fn erfcx(x: f64) -> f64 {
+    if x < 2.0 {
+        // Below the crossover neither factor is extreme: erfc is O(1) and
+        // exp(x^2) is at most e^4, so the direct product is exact enough and
+        // cheaper than the continued fraction.
+        libm::exp(x * x) * erfc(x)
+    } else {
+        // erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))),
+        // evaluated by backward recurrence. Converges quickly for x >= 2 and,
+        // unlike the product form, never touches an exponential.
+        let mut f = 0.0;
+        for n in (1..=60u32).rev() {
+            f = (f64::from(n) * 0.5) / (x + f);
+        }
+        FRAC_1_SQRT_PI / (x + f)
+    }
+}
+
+/// `ln` of the normal density at `x`.
+///
+/// The density itself underflows to zero past about 38 sigma, and `ln` of a
+/// clamped zero is -708 whatever the truth was. The log form is a polynomial:
+/// it stays exact at any separation, and the values it produces (-5001 nats at
+/// 100 sigma, -500001 at 1000) are perfectly representable.
+pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let z = (x - mu) / sigma;
+    -libm::log(SQRT_TAU * sigma) - 0.5 * z * z
+}
+
+/// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
+///
+/// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is
+/// factored out analytically via `erfcx`, so this never underflows — where
+/// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero.
+pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
+    let z = (x - mu) / sigma;
+
+    if z > 0.0 {
+        // ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y).
+        -std::f64::consts::LN_2 - 0.5 * z * z + libm::log(erfcx(z / SQRT_2))
+    } else {
+        // The mass here is at least a half; nothing to lose.
+        libm::log(sf(x, mu, sigma))
+    }
+}
+
+/// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`.
+///
+/// When the interval sits in a tail both endpoint probabilities underflow
+/// together, so their difference is taken in scaled form with the shared
+/// exponential factored out. When it straddles the mean nothing is small and
+/// the direct difference is exact.
+pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
+    let z_lo = (lo - mu) / sigma;
+    let z_hi = (hi - mu) / sigma;
+
+    if z_hi <= z_lo {
+        return f64::NEG_INFINITY;
+    }
+
+    // Fold a lower-tail interval onto the upper tail; the normal is symmetric.
+    let (near, far) = if z_lo >= 0.0 {
+        (z_lo, z_hi)
+    } else if z_hi <= 0.0 {
+        (-z_hi, -z_lo)
+    } else {
+        // Straddles the mean: the interval holds a non-negligible share of the
+        // mass, so neither endpoint is near enough to 1 to cancel.
+        return libm::log((cdf(hi, mu, sigma) - cdf(lo, mu, sigma)).max(f64::MIN_POSITIVE));
+    };
+
+    let (a, b) = (near / SQRT_2, far / SQRT_2);
+    // b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow.
+    let scale = libm::exp(a * a - b * b);
+    let bracket = erfcx(a) - scale * erfcx(b);
+
+    if bracket <= 0.0 {
+        return f64::NEG_INFINITY;
+    }
+
+    -std::f64::consts::LN_2 - a * a + libm::log(bracket)
+}
+
 fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
     let normalizer = (SQRT_TAU * sigma).powi(-1);
-    let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp();
+    let functional = libm::exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma));
 
     normalizer * functional
 }
 
+/// Truncated-Gaussian correction terms `(v, w)`.
+///
+/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose
+/// numerator and denominator underflow together in the tails, so both are
+/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled
+/// analytically rather than evaluated and divided out. Without that, a
+/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN`
+/// straight into the posterior.
+/// Truncation terms for a boundary `alpha` standard deviations into the upper
+/// tail, from the asymptotic expansion of the inverse Mills ratio.
+///
+/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by
+/// subtracting one from the other — the series computes the gap directly, and
+/// `w = v * gap` then never forms the difference of two large near-equal
+/// numbers. A far-tail *window* behaves like a half-line once it is more than a
+/// few multiples of its own width from the mean, so the tie branch shares this.
+fn half_line_truncation(alpha: f64) -> (f64, f64) {
+    let inv = alpha.recip();
+    let inv_sq = inv * inv;
+    let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
+    let v = alpha + gap;
+
+    (v, v * gap)
+}
+
 fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
     if !tie {
         let alpha = (margin - mu) / sigma;
 
-        let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
-        let w = v * (v + (-alpha));
+        // v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs
+        // the gap `v - alpha` as well as v itself. Far into the tail v tends to
+        // alpha, so that gap is a subtraction of two nearly equal numbers and
+        // loses every digit it has: at alpha = 1e6 it drove w above 1 and made
+        // `sqrt(1 - w)` NaN. Past the crossover the gap comes from its
+        // asymptotic series instead, which has no subtraction in it.
+        if alpha >= ASYMPTOTIC_MILLS_ALPHA {
+            return half_line_truncation(alpha);
+        }
 
-        (v, w)
+        let (v, gap) = if alpha > 0.0 {
+            // Both terms carry exp(-alpha^2 / 2); in scaled form it cancels
+            // and the result stays exact however far into the tail alpha sits.
+            let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2);
+            (v, v - alpha)
+        } else {
+            // Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing.
+            let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
+            (v, v - alpha)
+        };
+
+        (v, v * gap)
     } else {
+        // v is odd in mu and w is even, so fold to mu <= 0. Both truncation
+        // points then sit in the upper tail, where the scaled form applies.
+        let flipped = mu > 0.0;
+        let mu = if flipped { -mu } else { mu };
+
         let alpha = (-margin - mu) / sigma;
         let beta = (margin - mu) / sigma;
 
-        let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0))
-            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
-        let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0))
-            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
+        // `w` comes out of `v * v - u`, and both terms grow as alpha^2 while
+        // their difference stays O(1) — at alpha = 1e9 that subtraction had no
+        // digits left and returned w = -128, making `sqrt(1 - w)` nonsense.
+        // Once the window sits many of its own widths into the tail it is
+        // indistinguishable from a half-line, so the asymptotic covers it with
+        // no subtraction at all.
+        if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW {
+            let (v, w) = half_line_truncation(alpha);
+            return (if flipped { -v } else { v }, w);
+        }
+
+        let (v, u) = if alpha > 0.0 {
+            // beta > alpha > 0, so this ratio of exponentials is at most 1 and
+            // cannot overflow.
+            let scale = libm::exp(0.5 * (alpha * alpha - beta * beta));
+            let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
+
+            (
+                (1.0 - scale) / SQRT_TAU / denominator,
+                (alpha - beta * scale) / SQRT_TAU / denominator,
+            )
+        } else {
+            // The interval straddles the mean, so nothing here is small.
+            let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0);
+
+            (
+                (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator,
+                (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator,
+            )
+        };
+
         let w = -(u - v.powi(2));
 
-        (v, w)
+        (if flipped { -v } else { v }, w)
     }
 }
 
@@ -428,7 +698,7 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
     let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
     let s_arg = ata.determinant() / middle.determinant();
 
-    e_arg.exp() * s_arg.sqrt()
+    libm::exp(e_arg) * s_arg.sqrt()
 }
 
 #[cfg(test)]
@@ -442,6 +712,313 @@ mod tests {
         assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]);
     }
 
+    /// Upper-tail values of the standard normal, from published tables. The
+    /// point is not the digits — these are 7-digit table values — but that a
+    /// number comes back at all: `1 - cdf` returned exactly zero for every one
+    /// of these.
+    #[test]
+    fn survival_function_survives_the_far_tail() {
+        for (z, expected) in [
+            (9.0f64, 1.128_588e-19),
+            (12.0, 1.776_482e-33),
+            (20.0, 2.753_624e-89),
+            (37.0, 5.725_571e-300),
+        ] {
+            let got = sf(z, 0.0, 1.0);
+            assert!(got > 0.0, "sf({z}) collapsed to zero");
+            assert!(
+                (got - expected).abs() / expected < 1e-6, // published table values, 7 digits
+                "sf({z}) = {got}, expected ~{expected}"
+            );
+            assert_eq!(
+                1.0 - cdf(z, 0.0, 1.0),
+                0.0,
+                "the naive form should still be zero here"
+            );
+        }
+    }
+
+    /// Where no cancellation happens the two forms must agree exactly enough
+    /// that nothing else in the crate shifts.
+    #[test]
+    fn survival_function_matches_the_naive_form_where_that_form_works() {
+        for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] {
+            let naive = 1.0 - cdf(z, 0.0, 1.0);
+            let direct = sf(z, 0.0, 1.0);
+            assert!(
+                (naive - direct).abs() < 1e-15,
+                "z={z}: naive {naive} vs direct {direct}"
+            );
+        }
+    }
+
+    #[test]
+    fn survival_and_cdf_partition_the_mass() {
+        for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] {
+            let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0);
+            assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}");
+        }
+    }
+
+    /// `erfcx` switches formulation at x = 2; the two sides must meet.
+    #[test]
+    fn erfcx_is_continuous_across_its_crossover() {
+        for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] {
+            let direct = (x * x).exp() * erfc(x);
+            let scaled = erfcx(x);
+            assert!(
+                (direct - scaled).abs() / scaled < 1e-14,
+                "x={x}: direct {direct} vs erfcx {scaled}"
+            );
+        }
+    }
+
+    /// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly
+    /// where `exp(x^2)` overflows and `erfc(x)` underflows.
+    #[test]
+    fn erfcx_stays_finite_where_its_factors_do_not() {
+        for x in [27.0f64, 50.0, 1.0e3, 1.0e8] {
+            let scaled = erfcx(x);
+            assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}");
+            // Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)).
+            let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt());
+            assert!(
+                (scaled - asymptote).abs() / asymptote < 1e-2,
+                "erfcx({x}) = {scaled} strays from its asymptote {asymptote}"
+            );
+            assert!(
+                (x * x).exp().is_infinite(),
+                "x={x} should overflow the direct form"
+            );
+        }
+    }
+
+    /// Truncation must never produce a non-finite posterior. Before the scaled
+    /// formulation these returned NaN from `0 / 0` past about 39 sigma.
+    #[test]
+    fn truncation_stays_finite_arbitrarily_far_into_the_tail() {
+        for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] {
+            for tie in [false, true] {
+                let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
+                assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}");
+                assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}");
+                // sigma_trunc = sigma * sqrt(1 - w) must stay real.
+                assert!(
+                    (0.0..=1.0).contains(&w),
+                    "alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary"
+                );
+
+                let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
+                assert!(
+                    mu_t.is_finite() && sigma_t.is_finite(),
+                    "alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})"
+                );
+            }
+        }
+    }
+
+    /// The Mills gap switches from subtraction to series at alpha = 100. Both
+    /// are supposed to be right there; if they disagree, the crossover is in
+    /// the wrong place.
+    #[test]
+    fn the_mills_gap_series_meets_the_scaled_form() {
+        for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] {
+            let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha;
+            let inv = alpha.recip();
+            let inv_sq = inv * inv;
+            let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
+            assert!(
+                (scaled - series).abs() / series < 1e-9,
+                "alpha={alpha}: scaled {scaled} vs series {series}"
+            );
+        }
+    }
+
+    /// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and
+    /// w is even. Assert the symmetry the implementation relies on.
+    #[test]
+    fn tie_truncation_is_odd_in_v_and_even_in_w() {
+        for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] {
+            let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true);
+            let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true);
+            assert!(
+                (v_pos + v_neg).abs() < 1e-9,
+                "mu={mu}: v should be odd, got {v_pos} and {v_neg}"
+            );
+            assert!(
+                (w_pos - w_neg).abs() < 1e-9,
+                "mu={mu}: w should be even, got {w_pos} and {w_neg}"
+            );
+        }
+    }
+
+    /// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the
+    /// mirror image of the root. Three fixed iterations could not cross back,
+    /// so the error grew as the argument shrank: at `p_draw = 0.99` the margin
+    /// came out 0.503 where the answer is 2.576.
+    #[test]
+    fn erfc_inv_matches_known_quantiles() {
+        // sqrt(2) * erfc_inv(1 - p) is the standard normal quantile
+        // Phi^-1((1 + p) / 2).
+        for (p, exact) in [
+            (0.5f64, 0.674_489_750_196_081_7f64),
+            (0.9, 1.644_853_626_951_472_7),
+            (0.95, 1.959_963_984_540_054_2),
+            (0.99, 2.575_829_303_548_9),
+            (0.999, 3.290_526_731_491_896_4),
+        ] {
+            let got = SQRT_2 * erfc_inv(1.0 - p);
+            assert!(
+                (got - exact).abs() / exact < 1e-14,
+                "p={p}: got {got}, exact {exact}"
+            );
+        }
+    }
+
+    /// The draw margin must grow with the draw probability. It did not: it ran
+    /// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 ->
+    /// 0.999, which is not a rounding error but a broken function.
+    /// Deep in the tail the accuracy limit is the *caller's* argument, not this
+    /// function.
+    ///
+    /// `compute_margin(0.999999, ..)` computes `1.0 - p_draw`, and 0.999999 is
+    /// not representable: the subtraction cancels and leaves 2.9e-11 of
+    /// relative error in the argument before `erfc_inv` is even entered. Given
+    /// an exactly-representable argument the result is good to 1.8e-16, so this
+    /// is inherent to taking `p_draw` near one rather than something to fix
+    /// here. At `p_draw = 0.999` the whole path is still accurate to 4e-16.
+    ///
+    /// Worth pinning: measured against a 70-digit reference, `puruspe`'s
+    /// `inverfc` returns the identical wrong value for the identical reason,
+    /// which is what makes it clear the fault is upstream of both.
+    #[test]
+    fn erfc_inv_is_exact_given_an_exactly_representable_argument() {
+        // erfc(z / sqrt2) = 1e-6 exactly, so z = Phi^-1(0.9999995).
+        let got = SQRT_2 * erfc_inv(1e-6);
+        let exact = 4.891_638_475_698_59;
+        assert!(
+            (got - exact).abs() / exact < 1e-14,
+            "got {got}, exact {exact}"
+        );
+    }
+
+    #[test]
+    fn compute_margin_is_monotone_in_the_draw_probability() {
+        let mut previous = 0.0;
+        for p_draw in [
+            0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999,
+        ] {
+            let margin = compute_margin(p_draw, 1.0);
+            assert!(
+                margin > previous,
+                "p_draw={p_draw}: margin {margin} did not exceed {previous}"
+            );
+            previous = margin;
+        }
+    }
+
+    /// Round-tripping the margin back through the model's own CDF must recover
+    /// the draw probability it was built from.
+    #[test]
+    fn compute_margin_round_trips_through_the_cdf() {
+        for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] {
+            for sd in [0.5f64, 1.0, 5.892_557] {
+                let margin = compute_margin(p_draw, sd);
+                // P(|X| < margin) for X ~ N(0, sd^2).
+                let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd);
+                assert!(
+                    (recovered - p_draw).abs() < 1e-14,
+                    "p_draw={p_draw} sd={sd}: recovered {recovered}"
+                );
+            }
+        }
+    }
+
+    /// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where
+    /// the linear forms underflow. Past ~38 sigma the linear value is zero and
+    /// its log is whatever floor it was clamped to.
+    #[test]
+    fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() {
+        for z in [40.0f64, 60.0, 100.0, 1000.0] {
+            assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}");
+            assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}");
+
+            let lp = ln_pdf(z, 0.0, 1.0);
+            let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z;
+            assert!(
+                (lp - expected_lp).abs() < 1e-9,
+                "ln_pdf({z}) = {lp}, expected {expected_lp}"
+            );
+
+            let ls = ln_sf(z, 0.0, 1.0);
+            // ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z.
+            let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln();
+            assert!(
+                (ls - approx).abs() / approx.abs() < 1e-3,
+                "ln_sf({z}) = {ls}, asymptote {approx}"
+            );
+            assert!(
+                ls < f64::MIN_POSITIVE.ln(),
+                "ln_sf({z}) still on the clamp floor"
+            );
+        }
+    }
+
+    /// Where nothing underflows, the log helpers must agree with the direct
+    /// forms exactly enough that nothing else in the crate shifts.
+    #[test]
+    fn log_space_helpers_agree_with_the_linear_forms_in_range() {
+        for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] {
+            let lp = ln_pdf(z, 0.5, 2.0);
+            let direct_pdf = pdf(z, 0.5, 2.0);
+            assert!(
+                (lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf,
+                "ln_pdf at {z}: {} vs {direct_pdf}",
+                lp.exp()
+            );
+
+            let ls = ln_sf(z, 0.5, 2.0);
+            let direct = sf(z, 0.5, 2.0);
+            assert!(
+                (ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300),
+                "ln_sf at {z}: {} vs {direct}",
+                ls.exp()
+            );
+        }
+    }
+
+    #[test]
+    fn ln_interval_matches_the_direct_difference_when_nothing_is_small() {
+        for mu in [-2.0f64, 0.0, 0.5, 2.0] {
+            let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0);
+            let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp();
+            assert!(
+                (logged - direct).abs() <= 1e-13 * direct,
+                "mu={mu}: {logged} vs {direct}"
+            );
+        }
+    }
+
+    /// A window far out in the tail: both endpoints underflow together, so the
+    /// difference has to be taken in scaled form.
+    #[test]
+    fn ln_interval_survives_a_window_deep_in_the_tail() {
+        for mu in [-50.0f64, -100.0, -1000.0] {
+            let logged = ln_interval(-1.0, 1.0, mu, 1.0);
+            assert!(logged.is_finite(), "mu={mu}: {logged}");
+            assert!(
+                logged < f64::MIN_POSITIVE.ln(),
+                "mu={mu}: {logged} is stuck on the clamp floor"
+            );
+            // Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)).
+            let near = ln_sf(-1.0, mu, 1.0);
+            assert!(
+                (logged - near).abs() < 5.0,
+                "mu={mu}: {logged} strays from the near-edge tail {near}"
+            );
+        }
+    }
+
     #[test]
     fn test_quality() {
         let a = Gaussian::from_ms(25.0, 3.0);

src/observer.rs

@@ -14,13 +14,95 @@ pub trait Observer<T: Time>: Send + Sync {
     /// Called after each convergence iteration across the whole history.
     fn on_iteration_end(&self, _iter: usize, _max_step: (f64, f64)) {}
 
-    /// Called after each time slice is processed within an iteration.
-    fn on_batch_processed(&self, _time: &T, _slice_idx: usize, _n_events: usize) {}
+    /// Called after each time slice is swept within an iteration.
+    ///
+    /// A convergence iteration sweeps every slice twice — once travelling
+    /// backward through the history and once forward — so a multi-slice
+    /// history fires this twice per slice per iteration. A single-slice
+    /// history is swept once and fires once.
+    fn on_slice_processed(&self, _time: &T, _slice_idx: usize, _n_events: usize) {}
 
     /// Called once when convergence completes (or max iters is reached).
     fn on_converged(&self, _iters: usize, _final_step: (f64, f64), _converged: bool) {}
 }
 
+/// Shared and boxed observers forward to what they point at.
+///
+/// `History` takes its observer by value, so a caller who wants to *read* what
+/// an observer recorded has to keep a handle to it. Without these impls the
+/// natural spelling does not compile:
+///
+/// ```
+/// # use std::sync::{Arc, Mutex};
+/// # use trueskill_tt::{History, Observer};
+/// #[derive(Default)]
+/// struct Recorder {
+///     iterations: Mutex<Vec<usize>>,
+/// }
+///
+/// impl Observer<i64> for Recorder {
+///     fn on_iteration_end(&self, iter: usize, _step: (f64, f64)) {
+///         self.iterations.lock().unwrap().push(iter);
+///     }
+/// }
+///
+/// let recorder = Arc::new(Recorder::default());
+/// let mut h = History::builder().observer(Arc::clone(&recorder)).build();
+/// h.record_winner(&"a", &"b", 1).unwrap();
+/// h.converge().unwrap();
+///
+/// // The caller's handle sees what the history's copy recorded.
+/// assert!(!recorder.iterations.lock().unwrap().is_empty());
+/// ```
+///
+/// The alternative was for every observer to wrap each of its own fields in an
+/// `Arc` and derive `Clone` — one allocation and one lock per field, and a
+/// pattern each implementor had to rediscover.
+///
+/// `?Sized` is deliberate: it makes `Arc<dyn Observer<T>>` and
+/// `Box<dyn Observer<T>>` work, so observers can be chosen at runtime.
+impl<T: Time, O: Observer<T> + ?Sized> Observer<T> for std::sync::Arc<O> {
+    fn on_iteration_end(&self, iter: usize, max_step: (f64, f64)) {
+        (**self).on_iteration_end(iter, max_step);
+    }
+
+    fn on_slice_processed(&self, time: &T, slice_idx: usize, n_events: usize) {
+        (**self).on_slice_processed(time, slice_idx, n_events);
+    }
+
+    fn on_converged(&self, iters: usize, final_step: (f64, f64), converged: bool) {
+        (**self).on_converged(iters, final_step, converged);
+    }
+}
+
+impl<T: Time, O: Observer<T> + ?Sized> Observer<T> for Box<O> {
+    fn on_iteration_end(&self, iter: usize, max_step: (f64, f64)) {
+        (**self).on_iteration_end(iter, max_step);
+    }
+
+    fn on_slice_processed(&self, time: &T, slice_idx: usize, n_events: usize) {
+        (**self).on_slice_processed(time, slice_idx, n_events);
+    }
+
+    fn on_converged(&self, iters: usize, final_step: (f64, f64), converged: bool) {
+        (**self).on_converged(iters, final_step, converged);
+    }
+}
+
+impl<T: Time, O: Observer<T> + ?Sized> Observer<T> for &O {
+    fn on_iteration_end(&self, iter: usize, max_step: (f64, f64)) {
+        (**self).on_iteration_end(iter, max_step);
+    }
+
+    fn on_slice_processed(&self, time: &T, slice_idx: usize, n_events: usize) {
+        (**self).on_slice_processed(time, slice_idx, n_events);
+    }
+
+    fn on_converged(&self, iters: usize, final_step: (f64, f64), converged: bool) {
+        (**self).on_converged(iters, final_step, converged);
+    }
+}
+
 /// ZST no-op observer; the default when none is configured.
 #[derive(Copy, Clone, Debug, Default)]
 pub struct NullObserver;
@@ -35,6 +117,7 @@ mod tests {
     fn null_observer_compiles_for_i64() {
         let o = NullObserver;
         <NullObserver as Observer<i64>>::on_iteration_end(&o, 1, (0.0, 0.0));
+        <NullObserver as Observer<i64>>::on_slice_processed(&o, &7, 0, 3);
         <NullObserver as Observer<i64>>::on_converged(&o, 5, (1e-6, 1e-6), true);
     }
 

src/outcome.rs

@@ -34,12 +34,41 @@ impl Outcome {
     ///
     /// # Panics
     ///
-    /// Panics if `winner >= n`.
+    /// Panics if `winner >= n`. Use [`Outcome::try_winner`] when the index
+    /// comes from data rather than a literal.
+    ///
+    /// This is the one constructor here that validates, and deliberately so.
+    /// Its siblings build freely and let ingestion reject what it cannot use,
+    /// which works because a malformed rank vector stays recognisable. An
+    /// out-of-range winner does not: `winner(5, 2)` would produce ranks
+    /// `[1, 1]`, an all-tied draw that ingestion accepts without complaint when
+    /// `p_draw > 0`. Asking "team 5 won" and silently getting "everyone drew"
+    /// is exactly the class of quiet wrong answer this crate keeps removing, so
+    /// the check happens here where the mistake is.
     #[must_use]
     pub fn winner(winner: u32, n: u32) -> Self {
-        assert!(winner < n, "winner index {winner} out of range 0..{n}");
+        Self::try_winner(winner, n)
+            .unwrap_or_else(|_| panic!("winner index {winner} out of range 0..{n}"))
+    }
+
+    /// `n`-team outcome where team `winner` won, or an error if `winner` is not
+    /// a valid team index.
+    ///
+    /// The fallible form of [`Outcome::winner`], for when the index is computed
+    /// or parsed rather than written literally.
+    ///
+    /// # Errors
+    ///
+    /// `InvalidParameter` if `winner >= n`.
+    pub fn try_winner(winner: u32, n: u32) -> Result<Self, crate::InferenceError> {
+        if winner >= n {
+            return Err(crate::InferenceError::InvalidParameter {
+                name: "winner",
+                value: f64::from(winner),
+            });
+        }
         let ranks: SmallVec<[u32; 4]> = (0..n).map(|i| if i == winner { 0 } else { 1 }).collect();
-        Self::Ranked(ranks)
+        Ok(Self::Ranked(ranks))
     }
 
     /// All `n` teams tied.

src/predict.rs

@@ -0,0 +1,729 @@
+//! Outcome prediction: who wins, and how likely is a given finishing order.
+//!
+//! Prediction runs on *performances*, not skills. A competitor's skill is
+//! inflated by their performance noise `beta` before any comparison, which is
+//! what separates "how good are they" from "how will they do today".
+//!
+//! Two questions, two algorithms:
+//!
+//! - **Who finishes first.** Because performances are independent Gaussians,
+//!   the probability that team `i` beats every other team separates into a
+//!   *one-dimensional* integral — no multivariate orthant integral is
+//!   involved. [`quadrature::integrate`] evaluates it to near machine
+//!   precision for a few hundred `cdf` calls.
+//! - **A specific finishing order.** The factor graph only ever constrains
+//!   rank-*adjacent* teams (see `Game::run_chain`), so the joint probability
+//!   of a full order is a chain of local constraints rather than a general
+//!   orthant probability. That chain collapses into a sequential recursion:
+//!   one cumulative integral per adjacent pair, `O(teams * grid)` overall.
+//!
+//! Both are deterministic. A sampler would have been easier to write and
+//! would have made every `predict_*` call return a slightly different number,
+//! which is not a property a rating library should have.
+
+use crate::{Gaussian, quadrature};
+
+/// Teams beyond this count make the outcome enumeration impractical.
+///
+/// Each realisation sorts into exactly one (permutation, tie-pattern) event,
+/// so the space has `n! * 2^(n-1)` members: 24 at 3 teams, 192 at 4, 1_920 at
+/// 5, 23_040 at 6. The jump to 322_560 at 7 is where enumerating stops being
+/// a reasonable thing to do on a caller's behalf.
+pub(crate) const MAX_TEAMS_FOR_DISTRIBUTION: usize = 6;
+
+/// Relative tolerance for the first-place integrals.
+///
+/// The adaptive integrator reaches the exact two-team closed form to ~1e-15 at
+/// this tolerance, which is round-off for a probability. `cdf` is no longer the
+/// limit — it went to ~1 ULP when `erfc` moved to `libm` — so this is the
+/// integrator's own floor.
+const WIN_TOLERANCE: f64 = 1e-8;
+
+/// Nodes for the ranking grid, and the floor below which a grid is pointless.
+///
+/// The recursion converges as O(h^2), so this trades nodes against accuracy
+/// directly. Measured against the exact two-team closed form, 2_048 nodes leave
+/// ~1.2e-6 of discretisation error and 8_192 reach ~1e-7.
+///
+/// Unlike the adaptive path there is no approximation floor underneath this any
+/// more — `cdf` is accurate to ~1 ULP since `erfc` moved to `libm` — so the
+/// error here is purely the grid, and a caller who needs more can only get it
+/// by paying for more nodes. 8_192 is the accuracy/cost point chosen, not a
+/// point where refining stops helping.
+const MIN_GRID_POINTS: usize = 8_192;
+const MAX_GRID_POINTS: usize = 262_144;
+
+/// How many standard deviations of support the grid and integrals cover.
+///
+/// The normal density is below 1e-18 of its peak past nine sigma, far under
+/// the precision of everything else here.
+const SUPPORT_SIGMAS: f64 = 9.0;
+
+/// Standard normal CDF at `z`.
+fn phi(z: f64) -> f64 {
+    crate::cdf(z, 0.0, 1.0)
+}
+
+/// Normal density of `x` under `g`.
+fn density(g: Gaussian, x: f64) -> f64 {
+    let sigma = g.sigma();
+    let z = (x - g.mu()) / sigma;
+    libm::exp(-0.5 * z * z) / (sigma * (2.0 * std::f64::consts::PI).sqrt())
+}
+
+/// Per-pair draw margins.
+///
+/// The margin is *not* a single number for the whole game: inference derives
+/// it per rank-adjacent pair from those two teams' betas (`Game::likelihoods`).
+/// Prediction has to use the same per-pair values or it answers a question
+/// about a different model than the one that will actually be fitted.
+pub(crate) struct Margins {
+    n: usize,
+    values: Vec<f64>,
+}
+
+impl Margins {
+    /// Build from a per-pair margin function.
+    pub(crate) fn new<F: Fn(usize, usize) -> f64>(n: usize, f: F) -> Self {
+        let mut values = vec![0.0; n * n];
+        for i in 0..n {
+            for j in 0..n {
+                if i != j {
+                    values[i * n + j] = f(i, j);
+                }
+            }
+        }
+        Self { n, values }
+    }
+
+    fn get(&self, i: usize, j: usize) -> f64 {
+        self.values[i * self.n + j]
+    }
+
+    /// True when no pair can draw, so every tie has probability zero.
+    fn all_zero(&self) -> bool {
+        self.values.iter().all(|&v| v == 0.0)
+    }
+}
+
+/// `P(team i finishes strictly first)` for every team.
+///
+/// Strictly means beating each rival by more than that pair's draw margin, so
+/// with a non-zero margin these sum to less than one; the shortfall is the
+/// probability that the top place is shared.
+pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec<f64> {
+    (0..perf.len())
+        .map(|i| {
+            let (mu, sigma) = (perf[i].mu(), perf[i].sigma());
+            let (lo, hi) = (mu - SUPPORT_SIGMAS * sigma, mu + SUPPORT_SIGMAS * sigma);
+
+            // Each rival's CDF turns over near its own mean plus the margin.
+            // Seeding there is what keeps a rival with a tiny sigma — a step
+            // function in disguise — from being stepped over.
+            let mut seeds = Vec::with_capacity(3 * perf.len());
+            for (j, rival) in perf.iter().enumerate().filter(|&(j, _)| j != i) {
+                let centre = rival.mu() + margins.get(i, j);
+                seeds.extend_from_slice(&[centre - rival.sigma(), centre, centre + rival.sigma()]);
+            }
+
+            quadrature::integrate(
+                |x| {
+                    let d = density(perf[i], x);
+                    if d == 0.0 {
+                        return 0.0;
+                    }
+                    let beaten: f64 = (0..perf.len())
+                        .filter(|&j| j != i)
+                        .map(|j| phi((x - margins.get(i, j) - perf[j].mu()) / perf[j].sigma()))
+                        .product();
+                    d * beaten
+                },
+                lo,
+                hi,
+                &seeds,
+                WIN_TOLERANCE,
+            )
+        })
+        .collect()
+}
+
+/// Grid bounds and resolution covering every team's support.
+///
+/// Resolution is set by the *smallest* feature in play — the narrowest sigma,
+/// or a draw margin narrower still — because that is what the recursion has to
+/// resolve. A grid sized off the widest team would step over the narrow one.
+fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) {
+    let lo = perf
+        .iter()
+        .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma())
+        .fold(f64::INFINITY, f64::min);
+    let hi = perf
+        .iter()
+        .map(|g| g.mu() + SUPPORT_SIGMAS * g.sigma())
+        .fold(f64::NEG_INFINITY, f64::max);
+
+    let narrowest = perf
+        .iter()
+        .map(Gaussian::sigma)
+        .fold(f64::INFINITY, f64::min);
+    let smallest_margin = margins
+        .values
+        .iter()
+        .copied()
+        .filter(|&m| m > 0.0)
+        .fold(f64::INFINITY, f64::min);
+
+    let feature = narrowest.min(smallest_margin);
+    let wanted = if feature.is_finite() && feature > 0.0 {
+        ((hi - lo) / (feature / 12.0)).ceil()
+    } else {
+        MIN_GRID_POINTS as f64
+    };
+
+    let points = if wanted.is_finite() {
+        (wanted as usize).clamp(MIN_GRID_POINTS, MAX_GRID_POINTS)
+    } else {
+        MIN_GRID_POINTS
+    };
+
+    (lo, hi, points)
+}
+
+/// Densities of each team sampled on the shared grid.
+struct Sampled {
+    lo: f64,
+    step: f64,
+    points: usize,
+    density: Vec<Vec<f64>>,
+}
+
+impl Sampled {
+    fn new(perf: &[Gaussian], margins: &Margins) -> Self {
+        let (lo, hi, points) = grid_shape(perf, margins);
+        let step = (hi - lo) / (points - 1) as f64;
+        let density = perf
+            .iter()
+            .map(|&g| {
+                (0..points)
+                    .map(|i| density(g, lo + i as f64 * step))
+                    .collect()
+            })
+            .collect();
+        Self {
+            lo,
+            step,
+            points,
+            density,
+        }
+    }
+
+    fn node(&self, i: usize) -> f64 {
+        self.lo + i as f64 * self.step
+    }
+}
+
+/// `P(order[0] >= order[1] >= ... )` with the given adjacency pattern.
+///
+/// `tied[k]` says whether `order[k]` and `order[k + 1]` finish within that
+/// pair's draw margin. The recursion runs bottom-up: `carry` holds, for each
+/// grid node, the probability that everything *below* the current team holds
+/// given that team landed on that node. A strict gap reads a cumulative
+/// integral; a tie reads a window. Both are O(1) against one prefix array,
+/// so each level costs O(grid) and the whole order costs O(teams * grid).
+fn order_probability(margins: &Margins, sampled: &Sampled, order: &[usize], tied: &[bool]) -> f64 {
+    let mut carry = vec![1.0; sampled.points];
+
+    for k in (0..order.len() - 1).rev() {
+        let below = order[k + 1];
+        let above = order[k];
+        let margin = margins.get(above, below);
+
+        let integrand: Vec<f64> = (0..sampled.points)
+            .map(|i| sampled.density[below][i] * carry[i])
+            .collect();
+        let cumulative = quadrature::Grid::from_values(sampled.lo, sampled.step, integrand);
+
+        carry = (0..sampled.points)
+            .map(|i| {
+                let x = sampled.node(i);
+                if tied[k] {
+                    // Sorted order already implies `below <= above`, so the
+                    // tie window is one-sided: [x - margin, x].
+                    cumulative.integral_between(x - margin, x)
+                } else {
+                    cumulative.integral_to(x - margin)
+                }
+            })
+            .collect();
+    }
+
+    let top = order[0];
+    let integrand: Vec<f64> = (0..sampled.points)
+        .map(|i| sampled.density[top][i] * carry[i])
+        .collect();
+    quadrature::Grid::from_values(sampled.lo, sampled.step, integrand).total()
+}
+
+/// Dense ranks implied by a sorted order and its tie pattern.
+fn ranks_of(order: &[usize], tied: &[bool], n: usize) -> Vec<u32> {
+    let mut ranks = vec![0u32; n];
+    let mut rank = 0u32;
+    ranks[order[0]] = 0;
+    for k in 0..order.len() - 1 {
+        if !tied[k] {
+            rank += 1;
+        }
+        ranks[order[k + 1]] = rank;
+    }
+    ranks
+}
+
+/// Every (order, tie-pattern) event, or only the strict ones when no pair can
+/// draw — a tie then has probability exactly zero and is not worth integrating.
+fn events(n: usize, strict_only: bool) -> Vec<(Vec<usize>, Vec<bool>)> {
+    fn permute(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
+        if k == current.len() {
+            out.push(current.clone());
+            return;
+        }
+        for i in k..current.len() {
+            current.swap(k, i);
+            permute(current, k + 1, out);
+            current.swap(k, i);
+        }
+    }
+
+    let mut orders = Vec::new();
+    permute(&mut (0..n).collect(), 0, &mut orders);
+
+    let patterns: Vec<Vec<bool>> = if strict_only {
+        vec![vec![false; n - 1]]
+    } else {
+        (0..(1u32 << (n - 1)))
+            .map(|mask| (0..n - 1).map(|i| mask >> i & 1 == 1).collect())
+            .collect()
+    };
+
+    let mut out = Vec::with_capacity(orders.len() * patterns.len());
+    for order in orders {
+        for pattern in &patterns {
+            out.push((order.clone(), pattern.clone()));
+        }
+    }
+    out
+}
+
+/// The full distribution over finishing orders, aggregated by rank vector.
+///
+/// Orders that differ only *within* a tied group describe the same finishing
+/// order, so their probabilities are summed into one entry.
+pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<(Vec<u32>, f64)> {
+    let n = perf.len();
+    let sampled = Sampled::new(perf, margins);
+
+    let mut aggregated: Vec<(Vec<u32>, f64)> = Vec::new();
+    for (order, tied) in events(n, margins.all_zero()) {
+        let p = order_probability(margins, &sampled, &order, &tied);
+        let ranks = ranks_of(&order, &tied, n);
+        match aggregated.iter_mut().find(|(r, _)| *r == ranks) {
+            Some((_, acc)) => *acc += p,
+            None => aggregated.push((ranks, p)),
+        }
+    }
+
+    aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
+    aggregated
+}
+
+/// All permutations of `items`.
+fn permutations(items: &[usize]) -> Vec<Vec<usize>> {
+    fn go(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
+        if k == current.len() {
+            out.push(current.clone());
+            return;
+        }
+        for i in k..current.len() {
+            current.swap(k, i);
+            go(current, k + 1, out);
+            current.swap(k, i);
+        }
+    }
+    let mut out = Vec::new();
+    go(&mut items.to_vec(), 0, &mut out);
+    out
+}
+
+/// Every (order, tie-pattern) event consistent with a grouping by rank.
+///
+/// Teams sharing a rank may finish in any internal order, so this is the
+/// product of each group's permutations. Adjacencies inside a group are ties;
+/// the adjacency joining one group to the next is not.
+fn orders_for_groups(groups: &[Vec<usize>]) -> Vec<(Vec<usize>, Vec<bool>)> {
+    let per_group: Vec<Vec<Vec<usize>>> = groups.iter().map(|g| permutations(g)).collect();
+
+    let mut out = Vec::new();
+    let mut choice = vec![0usize; groups.len()];
+
+    loop {
+        let mut order = Vec::new();
+        let mut tied = Vec::new();
+        for (gi, group) in per_group.iter().enumerate() {
+            for (offset, &member) in group[choice[gi]].iter().enumerate() {
+                if !order.is_empty() {
+                    tied.push(offset != 0);
+                }
+                order.push(member);
+            }
+        }
+        out.push((order, tied));
+
+        let mut k = 0;
+        loop {
+            if k == choice.len() {
+                return out;
+            }
+            choice[k] += 1;
+            if choice[k] < per_group[k].len() {
+                break;
+            }
+            choice[k] = 0;
+            k += 1;
+        }
+    }
+}
+
+/// Probability of one specific rank vector.
+///
+/// Ties in `ranks` mean the tied teams may finish in any internal order, so
+/// this sums the orders consistent with the requested ranking rather than
+/// picking one.
+pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &[u32]) -> f64 {
+    let n = perf.len();
+    let sampled = Sampled::new(perf, margins);
+
+    let mut distinct: Vec<u32> = ranks.to_vec();
+    distinct.sort_unstable();
+    distinct.dedup();
+
+    let groups: Vec<Vec<usize>> = distinct
+        .iter()
+        .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect())
+        .collect();
+
+    orders_for_groups(&groups)
+        .iter()
+        .map(|(order, tied)| order_probability(margins, &sampled, order, tied))
+        .sum()
+}
+
+/// A distribution over the ways a contest could finish.
+///
+/// Each entry pairs a rank vector — the same shape [`crate::Outcome::ranking`]
+/// takes, with equal ranks meaning a tie — against its probability. Entries
+/// are ordered most likely first, and cover the whole outcome space, so the
+/// probabilities sum to one.
+///
+/// The rank vectors compose directly with inference: feeding one to
+/// `Game::ranked` asks "what would we believe if *this* happened", which is
+/// what an expected-information-gain calculation needs alongside the weight.
+#[derive(Clone, Debug, PartialEq)]
+pub struct Prediction {
+    outcomes: Vec<(Vec<u32>, f64)>,
+}
+
+impl Prediction {
+    pub(crate) fn new(outcomes: Vec<(Vec<u32>, f64)>) -> Self {
+        Self { outcomes }
+    }
+
+    /// Every possible finishing order and its probability, most likely first.
+    pub fn outcomes(&self) -> impl ExactSizeIterator<Item = (&[u32], f64)> {
+        self.outcomes.iter().map(|(r, p)| (r.as_slice(), *p))
+    }
+
+    /// The single most likely finishing order.
+    #[must_use]
+    pub fn most_likely(&self) -> Option<(&[u32], f64)> {
+        self.outcomes.first().map(|(r, p)| (r.as_slice(), *p))
+    }
+
+    /// Probability of one specific finishing order, or zero if it cannot occur.
+    #[must_use]
+    pub fn probability_of(&self, ranks: &[u32]) -> f64 {
+        self.outcomes
+            .iter()
+            .find(|(r, _)| r.as_slice() == ranks)
+            .map_or(0.0, |(_, p)| *p)
+    }
+
+    /// `P(team i finishes strictly first)`, for each team.
+    ///
+    /// Sums to less than one exactly when the top place can be shared; the
+    /// shortfall is [`Prediction::shared_first_place`].
+    #[must_use]
+    pub fn win_probabilities(&self) -> Vec<f64> {
+        let n = self.outcomes.first().map_or(0, |(r, _)| r.len());
+        let mut wins = vec![0.0; n];
+        for (ranks, p) in &self.outcomes {
+            let leaders = ranks.iter().filter(|&&r| r == 0).count();
+            if leaders == 1 {
+                let winner = ranks.iter().position(|&r| r == 0).expect("a rank-0 team");
+                wins[winner] += p;
+            }
+        }
+        wins
+    }
+
+    /// Probability that two or more teams share first place.
+    #[must_use]
+    pub fn shared_first_place(&self) -> f64 {
+        self.outcomes
+            .iter()
+            .filter(|(r, _)| r.iter().filter(|&&x| x == 0).count() > 1)
+            .map(|(_, p)| p)
+            .sum()
+    }
+
+    /// Total probability mass, which should be one.
+    ///
+    /// Exposed because it is a genuine check on the numerics rather than a
+    /// formality: the outcome space is exhaustive and disjoint by construction,
+    /// so any drift from one is integration error and nothing else.
+    #[must_use]
+    pub fn total(&self) -> f64 {
+        self.outcomes.iter().map(|(_, p)| p).sum()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn g(mu: f64, sigma: f64) -> Gaussian {
+        Gaussian::from_ms(mu, sigma)
+    }
+
+    fn flat(n: usize, eps: f64) -> Margins {
+        Margins::new(n, |_, _| eps)
+    }
+
+    /// Exact two-team result: `P(a first) = Phi((mu_a - mu_b - eps) / sd)`.
+    fn closed_form_two(a: Gaussian, b: Gaussian, eps: f64) -> (f64, f64) {
+        let sd = a.sigma().hypot(b.sigma());
+        (
+            phi((a.mu() - b.mu() - eps) / sd),
+            phi((b.mu() - a.mu() - eps) / sd),
+        )
+    }
+
+    #[test]
+    fn two_team_win_probabilities_match_the_closed_form() {
+        for (ma, sa, mb, sb, eps) in [
+            (0.0, 6.0, 0.0, 6.0, 0.0),
+            (3.0, 6.0, -2.0, 1.0, 0.0),
+            (0.0, 6.0, 0.0, 6.0, 2.0),
+            (3.0, 6.0, -2.0, 1.0, 1.5),
+            (40.0, 1.0, 0.0, 1.0, 0.0),
+        ] {
+            let perf = [g(ma, sa), g(mb, sb)];
+            let got = win_probabilities(&perf, &flat(2, eps));
+            let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
+            assert!(
+                (got[0] - wa).abs() < 1e-12 && (got[1] - wb).abs() < 1e-12,
+                "mu=({ma},{mb}) sigma=({sa},{sb}) eps={eps}: got {got:?}, want [{wa}, {wb}]"
+            );
+        }
+    }
+
+    /// The identity that a wrong-but-plausible implementation cannot fake:
+    /// with no draw margin, exactly one team finishes first.
+    #[test]
+    fn win_probabilities_sum_to_one_without_a_draw_margin() {
+        for perf in [
+            vec![g(0.0, 6.0), g(0.0, 6.0)],
+            vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)],
+            vec![
+                g(8.0, 2.0),
+                g(3.0, 6.0),
+                g(0.0, 1.0),
+                g(-3.0, 4.0),
+                g(-8.0, 6.0),
+            ],
+        ] {
+            let sum: f64 = win_probabilities(&perf, &flat(perf.len(), 0.0))
+                .iter()
+                .sum();
+            assert!(
+                (sum - 1.0).abs() < 1e-7,
+                "{} teams: sum = {sum}",
+                perf.len()
+            );
+        }
+    }
+
+    /// A rival with a tiny sigma is a step function in disguise. Fixed-node
+    /// quadrature steps over it and lands ~1e-2 out while still looking like a
+    /// probability; this is the case that rules that approach out.
+    #[test]
+    fn win_probabilities_survive_a_rival_with_a_tiny_sigma() {
+        let perf = [g(0.0, 0.001), g(0.5, 6.0), g(-0.5, 6.0)];
+        let got = win_probabilities(&perf, &flat(3, 0.0));
+        let sum: f64 = got.iter().sum();
+        assert!((sum - 1.0).abs() < 1e-6, "sum = {sum}, probs = {got:?}");
+    }
+
+    #[test]
+    fn a_stronger_team_is_more_likely_to_win() {
+        let perf = [g(10.0, 3.0), g(0.0, 3.0), g(-10.0, 3.0)];
+        let p = win_probabilities(&perf, &flat(3, 0.0));
+        assert!(p[0] > p[1] && p[1] > p[2], "not monotone: {p:?}");
+    }
+
+    #[test]
+    fn identical_teams_are_equally_likely_to_win() {
+        let perf = [g(1.0, 4.0), g(1.0, 4.0), g(1.0, 4.0)];
+        let p = win_probabilities(&perf, &flat(3, 0.0));
+        for probs in p.windows(2) {
+            assert!((probs[0] - probs[1]).abs() < 1e-9, "asymmetric: {p:?}");
+        }
+    }
+
+    /// Every realisation sorts into exactly one finishing order, so the whole
+    /// distribution must sum to one — with or without a draw margin.
+    #[test]
+    fn outcome_distribution_sums_to_one() {
+        for (perf, eps) in [
+            (vec![g(0.0, 6.0), g(0.0, 6.0)], 0.0),
+            (vec![g(0.0, 6.0), g(0.0, 6.0)], 2.0),
+            (vec![g(0.0, 6.0), g(0.0, 6.0), g(0.0, 6.0)], 0.0),
+            (vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)], 1.5),
+            (vec![g(0.0, 0.05), g(0.5, 6.0), g(-0.5, 6.0)], 1.0),
+            (
+                vec![g(6.0, 2.0), g(2.0, 6.0), g(-2.0, 1.0), g(-6.0, 4.0)],
+                1.0,
+            ),
+        ] {
+            let n = perf.len();
+            let dist = outcome_distribution(&perf, &flat(n, eps));
+            let sum: f64 = dist.iter().map(|(_, p)| p).sum();
+            assert!(
+                (sum - 1.0).abs() < 1e-6,
+                "{n} teams, eps={eps}: sum = {sum} over {} outcomes",
+                dist.len()
+            );
+            assert!(dist.iter().all(|(_, p)| *p >= 0.0), "negative probability");
+        }
+    }
+
+    /// With two teams the distribution is the exact win/draw/loss triple.
+    #[test]
+    fn two_team_distribution_matches_the_closed_form() {
+        let perf = [g(3.0, 6.0), g(-2.0, 1.0)];
+        let eps = 1.5;
+        let dist = outcome_distribution(&perf, &flat(2, eps));
+        let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
+
+        let find = |ranks: &[u32]| {
+            dist.iter()
+                .find(|(r, _)| r == ranks)
+                .map_or(0.0, |(_, p)| *p)
+        };
+
+        assert!(
+            (find(&[0, 1]) - wa).abs() < 1e-6,
+            "a wins: {}",
+            find(&[0, 1])
+        );
+        assert!(
+            (find(&[1, 0]) - wb).abs() < 1e-6,
+            "b wins: {}",
+            find(&[1, 0])
+        );
+        assert!(
+            (find(&[0, 0]) - (1.0 - wa - wb)).abs() < 1e-6,
+            "draw: {}",
+            find(&[0, 0])
+        );
+    }
+
+    /// Asking for one ranking must agree with that ranking's entry in the
+    /// full distribution — the two use different code paths to the same value.
+    #[test]
+    fn ranking_probability_agrees_with_the_distribution() {
+        let perf = [g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)];
+        let eps = 1.5;
+        let margins = flat(3, eps);
+        let dist = outcome_distribution(&perf, &margins);
+
+        for (ranks, expected) in &dist {
+            let direct = ranking_probability(&perf, &margins, ranks);
+            assert!(
+                (direct - expected).abs() < 1e-9,
+                "ranks {ranks:?}: direct {direct} vs distribution {expected}"
+            );
+        }
+    }
+
+    /// Tie mass is controlled by the draw margin. Only the *all-tied* outcome
+    /// is monotone in it: every one of its constraints is a window that widens
+    /// with the margin. A partially-tied outcome like `[0, 0, 1]` is not, and
+    /// must not be asserted to be — widening the margin makes its tie easier
+    /// but its "and the last team is strictly behind by more than the margin"
+    /// clause harder, so it peaks and then falls.
+    #[test]
+    fn all_tied_probability_grows_with_the_draw_margin() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
+        let mut previous = 0.0;
+        for eps in [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 24.0] {
+            let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]);
+            assert!(p >= previous, "eps={eps}: {p} < {previous}");
+            if eps == 0.0 {
+                assert!(p < 1e-12, "a tie needs a margin, got {p}");
+            }
+            previous = p;
+        }
+        assert!(
+            previous > 0.9,
+            "a very wide margin ties everyone: {previous}"
+        );
+    }
+
+    /// The converse, stated as the non-property it is: a partially-tied
+    /// outcome is non-monotone in the margin. Pinning this down stops a future
+    /// change from "fixing" it into monotonicity and quietly breaking the model.
+    #[test]
+    fn a_partially_tied_outcome_peaks_in_the_middle() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)];
+        let sweep: Vec<f64> = [0.5, 2.0, 4.0, 8.0, 16.0]
+            .iter()
+            .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]))
+            .collect();
+        let peak = sweep
+            .iter()
+            .enumerate()
+            .fold(
+                (0, 0.0),
+                |(bi, bv), (i, &v)| if v > bv { (i, v) } else { (bi, bv) },
+            )
+            .0;
+        assert!(
+            peak > 0 && peak < sweep.len() - 1,
+            "expected an interior peak: {sweep:?}"
+        );
+    }
+
+    /// With no draw margin a tie has probability exactly zero, and the
+    /// enumeration must not waste work pretending otherwise.
+    #[test]
+    fn ties_are_impossible_without_a_draw_margin() {
+        let perf = [g(0.0, 4.0), g(0.0, 4.0), g(0.0, 4.0)];
+        let dist = outcome_distribution(&perf, &flat(3, 0.0));
+        assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}");
+        assert!(dist.iter().all(|(r, _)| {
+            let mut seen = r.clone();
+            seen.sort_unstable();
+            seen.dedup();
+            seen.len() == r.len()
+        }));
+    }
+}

src/quadrature.rs

@@ -0,0 +1,322 @@
+//! Deterministic numerical integration for the prediction paths.
+//!
+//! Prediction asks two questions that have no closed form beyond two teams:
+//! "who finishes first" and "how likely is this exact finishing order". Both
+//! reduce to integrals over a single performance variable, so neither needs a
+//! sampler — and that matters, because a Monte Carlo predictor would make
+//! `predict_*` non-reproducible and would answer a slightly different question
+//! on every call.
+//!
+//! Two routines live here:
+//!
+//! - [`integrate`], adaptive Gauss-Kronrod G7-K15, for the first-place
+//!   marginals. It carries its own error estimate, so it can refine where the
+//!   integrand actually bends instead of guessing a node count up front.
+//! - [`Grid`], a uniform grid with trapezoid prefix sums, for the ranking
+//!   chain recursion, where each level needs the *running* integral of the
+//!   level below at arbitrary points rather than one definite integral.
+//!
+//! Fixed-node Gauss-Hermite is the obvious tool for the first of these and is
+//! a trap: the integrand is a product of normal CDFs, and when one team's
+//! sigma is much smaller than the integrating team's, that product turns into
+//! a near-step function narrower than the node spacing. The nodes step over
+//! it and the result is wrong by ~1e-2 while still looking like a probability.
+//! Adaptive refinement is what makes the small-sigma case safe.
+
+/// Kronrod 15-point abscissae, non-negative half, descending.
+const XGK: [f64; 8] = [
+    0.991_455_371_120_813,
+    0.949_107_912_342_759,
+    0.864_864_423_359_769,
+    0.741_531_185_599_394,
+    0.586_087_235_467_691,
+    0.405_845_151_377_397,
+    0.207_784_955_007_898,
+    0.0,
+];
+
+/// Kronrod 15-point weights, matching [`XGK`].
+const WGK: [f64; 8] = [
+    0.022_935_322_010_529,
+    0.063_092_092_629_979,
+    0.104_790_010_322_250,
+    0.140_653_259_715_525,
+    0.169_004_726_639_267,
+    0.190_350_578_064_785,
+    0.204_432_940_075_298,
+    0.209_482_141_084_728,
+];
+
+/// Gauss 7-point weights, applying to the odd-indexed [`XGK`] entries.
+const WG: [f64; 4] = [
+    0.129_484_966_168_870,
+    0.279_705_391_489_277,
+    0.381_830_050_505_119,
+    0.417_959_183_673_469,
+];
+
+/// Panels are bisected worst-first; this bounds the work on a pathological
+/// integrand rather than letting it spin.
+const MAX_SUBDIVISIONS: usize = 200;
+
+/// One G7-K15 panel over `[a, b]`: `(integral, absolute error estimate)`.
+///
+/// The error estimate is the gap between the embedded 7-point Gauss rule and
+/// the 15-point Kronrod extension. It is the only reason this is preferable
+/// to a fixed rule: it tells the caller *where* the integrand is hard.
+fn gk15<F: Fn(f64) -> f64>(f: &F, a: f64, b: f64) -> (f64, f64) {
+    let centre = 0.5 * (a + b);
+    let half = 0.5 * (b - a);
+
+    let mut kronrod = 0.0;
+    let mut gauss = 0.0;
+
+    for i in 0..8 {
+        let offset = XGK[i] * half;
+        // XGK[7] is the centre node and must not be counted twice.
+        let sum = if i == 7 {
+            f(centre)
+        } else {
+            f(centre - offset) + f(centre + offset)
+        };
+        kronrod += WGK[i] * sum;
+        if i % 2 == 1 {
+            gauss += WG[i / 2] * sum;
+        }
+    }
+
+    (kronrod * half, ((kronrod - gauss) * half).abs())
+}
+
+/// Adaptively integrate `f` over `[a, b]` to relative tolerance `tol`.
+///
+/// `seeds` are interior points where the integrand is known to bend sharply —
+/// for a product of normal CDFs, each rival's transition centre. Splitting
+/// there up front costs nothing and saves the adaptive loop from having to
+/// discover a step by bisection.
+///
+/// Returns the integral. The error estimate is consumed internally rather
+/// than returned: callers here integrate probability densities, where the
+/// meaningful check is the sum-to-one identity over a whole outcome space,
+/// not a per-integral residual.
+pub(crate) fn integrate<F: Fn(f64) -> f64>(f: F, a: f64, b: f64, seeds: &[f64], tol: f64) -> f64 {
+    // Explicit rather than `!(b > a)`: a NaN bound must fall through to zero
+    // rather than being read as a valid ordering.
+    if a.partial_cmp(&b) != Some(std::cmp::Ordering::Less) {
+        return 0.0;
+    }
+
+    let mut edges: Vec<f64> = Vec::with_capacity(seeds.len() + 2);
+    edges.push(a);
+    edges.push(b);
+    for &s in seeds {
+        if s > a && s < b {
+            edges.push(s);
+        }
+    }
+    edges.sort_by(|p, q| p.partial_cmp(q).expect("integration bounds are finite"));
+    edges.dedup();
+
+    // (lo, hi, integral, error)
+    let mut panels: Vec<(f64, f64, f64, f64)> = edges
+        .windows(2)
+        .map(|w| {
+            let (v, e) = gk15(&f, w[0], w[1]);
+            (w[0], w[1], v, e)
+        })
+        .collect();
+
+    for _ in 0..MAX_SUBDIVISIONS {
+        let total: f64 = panels.iter().map(|p| p.2).sum();
+        let error: f64 = panels.iter().map(|p| p.3).sum();
+
+        // Absolute floor as well as relative: these integrands are
+        // probabilities, so an absolute 1e-15 is already past the useful
+        // precision of the underlying `cdf`.
+        if error <= tol * total.abs().max(1e-12) || error < 1e-15 {
+            break;
+        }
+
+        let worst = panels
+            .iter()
+            .enumerate()
+            .fold((0usize, f64::NEG_INFINITY), |(bi, be), (i, p)| {
+                if p.3 > be { (i, p.3) } else { (bi, be) }
+            })
+            .0;
+
+        let (lo, hi, _, _) = panels[worst];
+        let mid = 0.5 * (lo + hi);
+        // Bisection has hit the floating-point floor; refining further would
+        // loop without reducing the error.
+        if !(mid > lo && mid < hi) {
+            break;
+        }
+
+        let (v1, e1) = gk15(&f, lo, mid);
+        let (v2, e2) = gk15(&f, mid, hi);
+        panels[worst] = (lo, mid, v1, e1);
+        panels.push((mid, hi, v2, e2));
+    }
+
+    panels.iter().map(|p| p.2).sum()
+}
+
+/// A uniform grid carrying trapezoid prefix sums of one integrand.
+///
+/// The ranking recursion needs, at every level, the running integral of the
+/// level below evaluated at arbitrary points — a cumulative integral, not a
+/// definite one. Prefix sums give that in O(1) per query after an O(G) build,
+/// which is what keeps a full ranking probability linear in the team count.
+pub(crate) struct Grid {
+    lo: f64,
+    step: f64,
+    /// Integrand sampled at each node.
+    values: Vec<f64>,
+    /// `prefix[i]` is the integral from `lo` to node `i`.
+    prefix: Vec<f64>,
+}
+
+impl Grid {
+    /// Build directly from already-sampled values.
+    ///
+    /// The ranking recursion evaluates every level on the same nodes, so the
+    /// per-team densities are sampled once and reused; re-evaluating `exp`
+    /// per level would dominate the cost.
+    pub(crate) fn from_values(lo: f64, step: f64, values: Vec<f64>) -> Self {
+        let mut prefix = vec![0.0; values.len()];
+        for i in 1..values.len() {
+            prefix[i] = prefix[i - 1] + 0.5 * step * (values[i - 1] + values[i]);
+        }
+        Self {
+            lo,
+            step,
+            values,
+            prefix,
+        }
+    }
+
+    /// Integral from the grid's lower bound up to `x`.
+    ///
+    /// Clamped at both ends: the caller sizes the grid to cover the whole
+    /// support, so a query outside it is asking for a tail that is zero (below)
+    /// or the whole mass (above).
+    pub(crate) fn integral_to(&self, x: f64) -> f64 {
+        let last = self.values.len() - 1;
+        if x <= self.lo {
+            return 0.0;
+        }
+        if x >= self.lo + last as f64 * self.step {
+            return self.prefix[last];
+        }
+
+        let scaled = (x - self.lo) / self.step;
+        let i = scaled.floor() as usize;
+        let frac = scaled - i as f64;
+
+        // Whole cells, plus the trapezoid over the partial cell. The integrand
+        // is linear within a cell under the trapezoid rule, so the partial
+        // piece is exact with respect to that same approximation.
+        self.prefix[i]
+            + frac
+                * self.step
+                * (self.values[i] + 0.5 * frac * (self.values[i + 1] - self.values[i]))
+    }
+
+    /// Integral over `[from, to]`.
+    pub(crate) fn integral_between(&self, from: f64, to: f64) -> f64 {
+        (self.integral_to(to) - self.integral_to(from)).max(0.0)
+    }
+
+    /// Total integral over the whole grid.
+    pub(crate) fn total(&self) -> f64 {
+        self.prefix[self.values.len() - 1]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    const TOL: f64 = 1e-10;
+
+    /// Sample `f` over `[lo, hi]` at `points` nodes.
+    fn sample<F: FnMut(f64) -> f64>(lo: f64, hi: f64, points: usize, mut f: F) -> Grid {
+        let step = (hi - lo) / (points - 1) as f64;
+        Grid::from_values(
+            lo,
+            step,
+            (0..points).map(|i| f(lo + i as f64 * step)).collect(),
+        )
+    }
+
+    #[test]
+    fn integrates_a_polynomial_exactly() {
+        // G7-K15 is exact for polynomials well past cubic, so a single panel
+        // should already be at round-off.
+        let v = integrate(|x| 3.0 * x * x + 2.0 * x + 1.0, 0.0, 2.0, &[], TOL);
+        assert!((v - 14.0).abs() < 1e-12, "got {v}");
+    }
+
+    #[test]
+    fn integrates_a_gaussian_density_to_one() {
+        let f = |x: f64| (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt();
+        let v = integrate(f, -10.0, 10.0, &[], TOL);
+        assert!((v - 1.0).abs() < 1e-12, "got {v}");
+    }
+
+    #[test]
+    fn resolves_a_step_far_narrower_than_the_initial_panel() {
+        // The failure mode that rules out fixed-node quadrature: a transition
+        // 1e-4 wide inside a range of 20. A fixed rule steps over it.
+        let f = |x: f64| if x < 0.5 { 0.0 } else { 1.0 };
+        let v = integrate(f, -10.0, 10.0, &[0.5], TOL);
+        assert!((v - 9.5).abs() < 1e-6, "got {v}");
+    }
+
+    #[test]
+    fn seeds_do_not_change_the_value_of_a_smooth_integrand() {
+        let f = |x: f64| (-0.5 * x * x).exp();
+        let plain = integrate(f, -8.0, 8.0, &[], TOL);
+        let seeded = integrate(f, -8.0, 8.0, &[-3.0, 0.25, 5.5], TOL);
+        assert!((plain - seeded).abs() < 1e-12, "{plain} vs {seeded}");
+    }
+
+    #[test]
+    fn empty_or_inverted_range_integrates_to_zero() {
+        assert_eq!(integrate(|_| 1.0, 1.0, 1.0, &[], TOL), 0.0);
+        assert_eq!(integrate(|_| 1.0, 2.0, 1.0, &[], TOL), 0.0);
+    }
+
+    #[test]
+    fn grid_prefix_matches_a_known_cumulative_integral() {
+        // f(x) = x over [0, 4]; integral to x is x^2/2.
+        let g = sample(0.0, 4.0, 4001, |x| x);
+        for probe in [0.0, 0.5, 1.0, 2.5, 3.75, 4.0] {
+            let want = probe * probe / 2.0;
+            let got = g.integral_to(probe);
+            assert!(
+                (got - want).abs() < 1e-9,
+                "at {probe}: got {got}, want {want}"
+            );
+        }
+        assert!((g.total() - 8.0).abs() < 1e-9);
+    }
+
+    #[test]
+    fn grid_between_is_the_difference_of_two_prefixes() {
+        let g = sample(-5.0, 5.0, 8001, |x| (-0.5 * x * x).exp());
+        let whole = g.integral_between(-5.0, 5.0);
+        let split = g.integral_between(-5.0, 0.3) + g.integral_between(0.3, 5.0);
+        assert!((whole - split).abs() < 1e-12, "{whole} vs {split}");
+    }
+
+    #[test]
+    fn grid_clamps_queries_outside_its_support() {
+        let g = sample(0.0, 1.0, 101, |_| 1.0);
+        assert_eq!(g.integral_to(-3.0), 0.0);
+        assert!((g.integral_to(9.0) - 1.0).abs() < 1e-12);
+        // Reversed bounds must not produce negative probability mass.
+        assert_eq!(g.integral_between(0.8, 0.2), 0.0);
+    }
+}

src/rating.rs

@@ -9,8 +9,8 @@ use crate::{
 
 /// Static rating configuration: prior skill, performance noise `beta`, drift.
 ///
-/// Renamed from `Player` in T2; `Rating` better describes the data
-/// (a configuration) vs. a person (who's a `Competitor` with state).
+/// A configuration rather than a person: the per-history temporal state
+/// (messages, last appearance) lives on `Competitor`.
 #[derive(Clone, Copy, Debug)]
 pub struct Rating<T: Time = i64, D: Drift<T> = ConstantDrift> {
     pub(crate) prior: Gaussian,

src/schedule.rs

@@ -1,152 +0,0 @@
-//! Schedule trait and built-in implementations.
-//!
-//! A schedule drives factor propagation to convergence. The default
-//! `EpsilonOrMax` performs one `TeamSum` sweep (setup) then alternating
-//! forward/backward sweeps over the iterating factors until the max
-//! delta drops below epsilon or `max` iterations is reached.
-
-use crate::factor::{BuiltinFactor, Factor, VarStore};
-
-/// Result returned by a `Schedule::run` call.
-#[derive(Debug, Clone, Copy)]
-pub struct ScheduleReport {
-    pub iterations: usize,
-    pub final_step: (f64, f64),
-    pub converged: bool,
-}
-
-/// Drives factor propagation to convergence.
-pub trait Schedule: Send + Sync {
-    fn run(&self, factors: &mut [BuiltinFactor], vars: &mut VarStore) -> ScheduleReport;
-}
-
-/// Default schedule: sweep forward then backward until step ≤ eps or iter == max.
-///
-/// Matches the existing `Game::likelihoods` loop bit-for-bit when given the
-/// same factor layout (`TeamSums` first, then alternating RankDiff/Trunc pairs).
-#[derive(Debug, Clone, Copy)]
-pub struct EpsilonOrMax {
-    pub eps: f64,
-    pub max: usize,
-}
-
-impl Default for EpsilonOrMax {
-    fn default() -> Self {
-        // 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,
-        }
-    }
-}
-
-impl Schedule for EpsilonOrMax {
-    fn run(&self, factors: &mut [BuiltinFactor], vars: &mut VarStore) -> ScheduleReport {
-        // Partition: leading run of TeamSum factors run exactly once (setup).
-        let n_setup = factors
-            .iter()
-            .position(|f| !matches!(f, BuiltinFactor::TeamSum(_)))
-            .unwrap_or(factors.len());
-
-        for f in factors[..n_setup].iter_mut() {
-            f.propagate(vars);
-        }
-
-        let mut iterations = 0;
-        // 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);
-
-                // Forward sweep over iterating factors.
-                for f in factors[n_setup..].iter_mut() {
-                    let d = f.propagate(vars);
-                    step.0 = step.0.max(d.0);
-                    step.1 = step.1.max(d.1);
-                }
-
-                // Backward sweep.
-                for f in factors[n_setup..].iter_mut().rev() {
-                    let d = f.propagate(vars);
-                    step.0 = step.0.max(d.0);
-                    step.1 = step.1.max(d.1);
-                }
-
-                iterations += 1;
-                final_step = step;
-
-                if step.0 <= self.eps && step.1 <= self.eps {
-                    converged = true;
-                    break;
-                }
-            }
-        }
-
-        ScheduleReport {
-            iterations,
-            final_step,
-            converged,
-        }
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::{N_INF, factor::team_sum::TeamSumFactor, gaussian::Gaussian};
-
-    #[test]
-    fn schedule_runs_setup_factors_once() {
-        // Single TeamSum factor; schedule should propagate it exactly once and report 0 iterations.
-        let mut vars = VarStore::new();
-        let out = vars.alloc(N_INF);
-        let mut factors = vec![BuiltinFactor::TeamSum(TeamSumFactor {
-            inputs: vec![(Gaussian::from_ms(5.0, 1.0), 1.0)],
-            out,
-        })];
-        let schedule = EpsilonOrMax::default();
-        let report = schedule.run(&mut factors, &mut vars);
-        assert_eq!(report.iterations, 0);
-        // The team-perf var should hold the sum.
-        let result = vars.get(out);
-        assert!((result.mu() - 5.0).abs() < 1e-12);
-    }
-
-    #[test]
-    fn report_marks_converged_when_no_iterating_factors() {
-        // 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 {
-            inputs: vec![(Gaussian::from_ms(0.0, 1.0), 1.0)],
-            out,
-        })];
-        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/time_slice.rs

@@ -809,6 +809,93 @@ pub(crate) fn compute_elapsed<T: Time>(last: Option<&T>, current: &T) -> i64 {
     elapsed.max(0)
 }
 
+impl<T: Time> TimeSlice<T> {
+    /// Precision matrix of the joint posterior over this slice's competitors.
+    ///
+    /// Message passing produces per-competitor marginals and throws the
+    /// correlation away — `Item::likelihood` is already the projection of an
+    /// event's factor down onto one competitor. So the joint has to be rebuilt
+    /// from the factor structure rather than recovered from the messages.
+    ///
+    /// Usefully, a precision matrix depends only on *structure* — who played
+    /// whom, with what weights and what observation noise — and not on the
+    /// observed outcomes. The means are already exact (Gaussian belief
+    /// propagation gets those right even with cycles), so only the second
+    /// moment needs rebuilding.
+    ///
+    /// Returns the competitor order and the dense matrix in row-major order.
+    /// Only scored events contribute their factors exactly; see the caller.
+    pub(crate) fn joint_precision<D: Drift<T>>(
+        &self,
+        agents: &CompetitorStore<T, D>,
+    ) -> (Vec<Index>, Vec<f64>) {
+        let order: Vec<Index> = self.skills.keys().collect();
+        let n = order.len();
+        let mut row_of: HashMap<Index, usize> = HashMap::with_capacity(n);
+        for (r, idx) in order.iter().enumerate() {
+            row_of.insert(*idx, r);
+        }
+
+        let mut lambda = vec![0.0; n * n];
+
+        // Everything outside this slice enters as each competitor's forward and
+        // backward messages, which message passing treats as independent.
+        for (r, idx) in order.iter().enumerate() {
+            let skill = self.skills.get(*idx).expect("slice key has a skill");
+            lambda[r * n + r] += (skill.forward * skill.backward).pi();
+        }
+
+        for event in &self.events {
+            let EventKind::Scored { score_sigma } = event.kind else {
+                continue;
+            };
+
+            // Teams best-first, matching the diff chain inference builds.
+            let mut order_idx: Vec<usize> = (0..event.teams.len()).collect();
+            order_idx.sort_by(|&a, &b| {
+                event.teams[b]
+                    .output
+                    .partial_cmp(&event.teams[a].output)
+                    .unwrap_or(std::cmp::Ordering::Equal)
+            });
+
+            for pair in order_idx.windows(2) {
+                let (hi, lo) = (pair[0], pair[1]);
+
+                // Contrast vector, and the observation noise that sits on top
+                // of the skills: per-member performance noise plus the score
+                // noise itself.
+                let mut contrast: HashMap<usize, f64> = HashMap::new();
+                let mut noise = score_sigma * score_sigma;
+
+                for (team, sign) in [(hi, 1.0), (lo, -1.0)] {
+                    for (m, item) in event.teams[team].items.iter().enumerate() {
+                        let w = event.weights[team][m];
+                        let beta = agents[item.agent].rating.beta;
+                        noise += w * w * beta * beta;
+                        *contrast.entry(row_of[&item.agent]).or_insert(0.0) += sign * w;
+                    }
+                }
+
+                for (&i, &ci) in &contrast {
+                    for (&j, &cj) in &contrast {
+                        lambda[i * n + j] += ci * cj / noise;
+                    }
+                }
+            }
+        }
+
+        (order, lambda)
+    }
+
+    /// True when every event here is scored, so `joint_precision` is exact.
+    pub(crate) fn all_scored(&self) -> bool {
+        self.events
+            .iter()
+            .all(|e| matches!(e.kind, EventKind::Scored { .. }))
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use approx::assert_ulps_eq;

tests/additive_model.rs

@@ -0,0 +1,142 @@
+//! What an additive model does to uncertainty, and why "add the marginals" is
+//! unsafe in one direction and merely wasteful in the other.
+//!
+//! Structurally this is the shape a joint player/layout model takes: every
+//! observation measures a *sum* of nodes against a reference, so the data pins
+//! differences and leaves the overall level to the prior. That is the classic
+//! rating-scale indeterminacy, not a defect.
+//!
+//! The consequence for a consumer is that combining marginals is wrong in
+//! opposite directions depending on the combination, which is worth pinning
+//! because the unsafe direction is not the one you would guess:
+//!
+//! - **Differences** (`a - b`): the shared level cancels, so the exact width is
+//!   small — and adding marginals lands within a couple of percent of it here,
+//!   because the loopy underestimate offsets the ignored correlation.
+//! - **Sums** (`a + b`): the shared level does *not* cancel, so the exact width
+//!   is large, and adding marginals is roughly five times too narrow. That is
+//!   overconfident, and it is the direction that publishes a claim the data
+//!   does not support.
+
+use smallvec::smallvec;
+use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team};
+
+#[test]
+fn additive_structure_makes_sums_wide_and_differences_tight() {
+    // Structurally like ustat: every round is (player + hole) measured against
+    // a fixed reference. Only SUMS are pinned by the data; the split between
+    // player and hole is pinned only by the prior.
+    let players = ["p0", "p1", "p2"];
+    let holes = ["h0", "h1"];
+
+    let mut h: History<i64, _, _, &'static str> = History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.0))
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-12,
+            alpha: 1.0,
+        })
+        .build();
+
+    let mut seed = 3u64;
+    let mut rnd = move || {
+        seed ^= seed << 13;
+        seed ^= seed >> 7;
+        seed ^= seed << 17;
+        seed
+    };
+    // true skills, so we know what the data encodes
+    let truth_p = [2.0, 0.0, -2.0];
+    let truth_h = [1.0, -1.0];
+
+    let mut events = Vec::new();
+    for _ in 0..60 {
+        let p = (rnd() as usize) % 3;
+        let q = (rnd() as usize) % 2;
+        let noise = ((rnd() % 1000) as f64 / 1000.0 - 0.5) * 2.0;
+        let score = truth_p[p] + truth_h[q] + noise;
+        events.push(Event {
+            time: 1,
+            teams: smallvec![
+                Team::with_members([Member::new(players[p]), Member::new(holes[q])]),
+                Team::with_members([Member::new("reference")]),
+            ],
+            outcome: Outcome::scores([score, 0.0]),
+        });
+    }
+    h.add_events(events).unwrap();
+    let r = h.converge().unwrap();
+    assert!(r.converged, "{:?}", r.final_step);
+
+    println!("\n== marginals (what current_skill reports) ==");
+    for k in players.iter().chain(holes.iter()) {
+        let g = h.current_skill(k).unwrap();
+        println!("  {k}: mu {:>8.4}  sigma {:>8.4}", g.mu(), g.sigma());
+    }
+
+    println!("\n== the same nodes via posterior_of (exact marginal) ==");
+    for k in players.iter().chain(holes.iter()) {
+        let g = h.posterior_of(&[(k, 1.0)]).unwrap();
+        println!("  {k}: mu {:>8.4}  sigma {:>8.4}", g.mu(), g.sigma());
+    }
+
+    println!("\n== combinations the data actually pins ==");
+    for (label, terms) in [
+        ("p0 + h0 (a round)", vec![(&"p0", 1.0), (&"h0", 1.0)]),
+        (
+            "p0 - p1 (rank two players)",
+            vec![(&"p0", 1.0), (&"p1", -1.0)],
+        ),
+        ("p0 - p2", vec![(&"p0", 1.0), (&"p2", -1.0)]),
+        (
+            "h0 - h1 (rank two holes)",
+            vec![(&"h0", 1.0), (&"h1", -1.0)],
+        ),
+    ] {
+        let joint = h.posterior_of(&terms).unwrap();
+        // what a consumer gets today by adding marginals
+        let naive: f64 = terms
+            .iter()
+            .map(|(k, c)| c * c * h.current_skill(*k).unwrap().sigma().powi(2))
+            .sum::<f64>()
+            .sqrt();
+        println!(
+            "  {label:<28} exact sigma {:>7.4}   adding marginals {:>7.4}   {:>5.2}x over",
+            joint.sigma(),
+            naive,
+            naive / joint.sigma()
+        );
+
+        let ratio = naive / joint.sigma();
+        if label.contains('+') {
+            assert!(
+                ratio < 0.5,
+                "{label}: adding marginals should be badly OVERconfident for a \
+                 sum, got {ratio:.3}x"
+            );
+        } else {
+            assert!(
+                (0.8..1.25).contains(&ratio),
+                "{label}: adding marginals happens to be close for a difference, \
+                 got {ratio:.3}x"
+            );
+        }
+    }
+
+    // A single node in an additive model is weakly identified: its exact
+    // posterior is far wider than message passing reports, because the level it
+    // shares with its partners is pinned only by the prior.
+    for k in players.iter().chain(holes.iter()) {
+        let bp = h.current_skill(k).unwrap().sigma();
+        let exact = h.posterior_of(&[(k, 1.0)]).unwrap().sigma();
+        assert!(
+            exact > 3.0 * bp,
+            "{k}: exact marginal {exact} should be much wider than the reported \
+             {bp} in an additive model"
+        );
+    }
+}

tests/api_shape.rs

@@ -65,7 +65,7 @@ fn add_events_draw() {
         outcome: Outcome::draw(2),
     }];
     h.add_events(events).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 }
 
 #[test]
@@ -123,7 +123,7 @@ fn fluent_event_builder_winner_convenience() {
         .winner(0)
         .commit()
         .unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 }
 
 #[test]
@@ -141,7 +141,7 @@ fn fluent_event_builder_draw() {
         .draw()
         .commit()
         .unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 }
 
 #[test]
@@ -155,7 +155,7 @@ fn current_skill_and_learning_curve() {
         .build();
     h.record_winner(&"a", &"b", 1).unwrap();
     h.record_winner(&"a", &"b", 2).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let a = h.current_skill(&"a").unwrap();
     assert!(a.mu() > 25.0);
@@ -201,9 +201,9 @@ fn predict_quality_two_teams() {
         .p_draw(0.0)
         .build();
     h.record_winner(&"a", &"b", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
-    let q = h.predict_quality(&[&[&"a"], &[&"b"]]);
+    let q = h.predict_quality(&[&[&"a"], &[&"b"]]).unwrap();
     assert!(q > 0.0 && q <= 1.0);
 }
 
@@ -217,12 +217,16 @@ fn predict_outcome_two_teams_sums_to_one() {
         .p_draw(0.0)
         .build();
     h.record_winner(&"a", &"b", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
-    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]);
-    assert_eq!(p.len(), 2);
-    assert!((p[0] + p[1] - 1.0).abs() < 1e-9);
-    assert!(p[0] > p[1]);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+    let wins = p.win_probabilities();
+    assert_eq!(wins.len(), 2);
+    // With p_draw == 0 there is no draw outcome, so the two win
+    // probabilities are the whole space.
+    assert!((p.total() - 1.0).abs() < 1e-9, "total = {}", p.total());
+    assert!((wins[0] + wins[1] - 1.0).abs() < 1e-9);
+    assert!(wins[0] > wins[1]);
 }
 
 #[test]
@@ -241,7 +245,7 @@ fn fluent_event_builder_scores() {
         .scores([12.0, 4.0])
         .commit()
         .unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let a = h.current_skill(&"alice").unwrap();
     let b = h.current_skill(&"bob").unwrap();

tests/competitor_config.rs

@@ -0,0 +1,222 @@
+//! `Member::with_prior` / `with_drift_scale` — competitor configuration.
+//!
+//! Both were previously consumed only on the branch that *creates* a
+//! competitor, so configuration supplied for a key the history already knew was
+//! dropped with no error. `with_prior` had no coverage in this directory at
+//! all, which is how that survived.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConvergenceOptions, Event, Gaussian, History, InferenceError, Member, Outcome, Team,
+};
+
+const CONVERGENCE: ConvergenceOptions = ConvergenceOptions {
+    max_iter: 2_000,
+    epsilon: 1e-12,
+    alpha: 1.0,
+};
+
+fn history() -> History {
+    History::builder()
+        .mu(25.0)
+        .sigma(25.0 / 3.0)
+        .beta(25.0 / 6.0)
+        .p_draw(0.0)
+        .convergence(CONVERGENCE)
+        .build()
+}
+
+/// One event, optionally configuring `a`.
+fn bout(
+    a: &'static str,
+    b: &'static str,
+    time: i64,
+    prior: Option<Gaussian>,
+    scale: Option<f64>,
+) -> Event<i64, &'static str> {
+    let mut member = Member::new(a);
+    if let Some(p) = prior {
+        member = member.with_prior(p);
+    }
+    if let Some(s) = scale {
+        member = member.with_drift_scale(s);
+    }
+
+    Event {
+        time,
+        teams: smallvec![
+            Team::with_members([member]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::winner(0, 2),
+    }
+}
+
+fn skill_of(h: &History, key: &str) -> Gaussian {
+    h.current_skill(&key).expect("key in history")
+}
+
+/// Baseline: the mechanism works at all on a competitor's first appearance.
+#[test]
+fn a_prior_applies_to_a_new_competitor() {
+    let seeded = Gaussian::from_ms(40.0, 1.0);
+
+    let mut with = history();
+    with.add_events(vec![bout("a", "b", 0, Some(seeded), None)])
+        .unwrap();
+    let _ = with.converge().unwrap();
+
+    let mut without = history();
+    without
+        .add_events(vec![bout("a", "b", 0, None, None)])
+        .unwrap();
+    let _ = without.converge().unwrap();
+
+    assert!(
+        (skill_of(&with, "a").mu() - skill_of(&without, "a").mu()).abs() > 1.0,
+        "a seeded prior should move the fit"
+    );
+}
+
+/// The defect in #10: a prior supplied for a competitor the history already
+/// knows was silently discarded, and the caller got output computed from the
+/// default prior with no indication anything had been dropped.
+#[test]
+fn a_prior_applies_to_a_competitor_the_history_already_knows() {
+    let seeded = Gaussian::from_ms(40.0, 1.0);
+
+    let mut late = history();
+    late.add_events(vec![bout("a", "b", 0, None, None)])
+        .unwrap();
+    // "a" now exists. Configuring it here used to do nothing whatsoever.
+    late.add_events(vec![bout("a", "b", 1, Some(seeded), None)])
+        .unwrap();
+    let _ = late.converge().unwrap();
+
+    let mut never = history();
+    never
+        .add_events(vec![
+            bout("a", "b", 0, None, None),
+            bout("a", "b", 1, None, None),
+        ])
+        .unwrap();
+    let _ = never.converge().unwrap();
+
+    assert!(
+        (skill_of(&late, "a").mu() - skill_of(&never, "a").mu()).abs() > 1.0,
+        "a late prior must not be silently dropped: {} vs {}",
+        skill_of(&late, "a").mu(),
+        skill_of(&never, "a").mu()
+    );
+}
+
+/// Configuration is competitor-scoped, not event-scoped, and `converge` refits
+/// from competitor state — so seeding late reaches the same fit as seeding from
+/// the start. This is the documented scope, asserted rather than assumed.
+#[test]
+fn a_prior_is_whole_history_scoped_not_per_event() {
+    let seeded = Gaussian::from_ms(40.0, 1.0);
+
+    let mut late = history();
+    late.add_events(vec![bout("a", "b", 0, None, None)])
+        .unwrap();
+    late.add_events(vec![bout("a", "b", 1, Some(seeded), None)])
+        .unwrap();
+    let _ = late.converge().unwrap();
+
+    let mut early = history();
+    early
+        .add_events(vec![
+            bout("a", "b", 0, Some(seeded), None),
+            bout("a", "b", 1, Some(seeded), None),
+        ])
+        .unwrap();
+    let _ = early.converge().unwrap();
+
+    let (l, e) = (skill_of(&late, "a"), skill_of(&early, "a"));
+    assert!(
+        (l.mu() - e.mu()).abs() < 1e-9 && (l.sigma() - e.sigma()).abs() < 1e-9,
+        "late seeding should refit the whole history: {l:?} vs {e:?}"
+    );
+}
+
+#[test]
+fn repeating_the_same_prior_is_inert() {
+    let seeded = Gaussian::from_ms(40.0, 1.0);
+
+    let mut once = history();
+    once.add_events(vec![
+        bout("a", "b", 0, Some(seeded), None),
+        bout("a", "b", 1, None, None),
+    ])
+    .unwrap();
+    let _ = once.converge().unwrap();
+
+    let mut every_time = history();
+    every_time
+        .add_events(vec![
+            bout("a", "b", 0, Some(seeded), None),
+            bout("a", "b", 1, Some(seeded), None),
+        ])
+        .unwrap();
+    let _ = every_time.converge().unwrap();
+
+    let (o, e) = (skill_of(&once, "a"), skill_of(&every_time, "a"));
+    assert!(
+        (o.mu() - e.mu()).abs() < 1e-12 && (o.sigma() - e.sigma()).abs() < 1e-12,
+        "declaring the same prior repeatedly changed the fit: {o:?} vs {e:?}"
+    );
+}
+
+/// Events within a batch have no order, so two different values for one
+/// competitor have no well-defined winner. Rejecting is what keeps the answer
+/// independent of iteration order.
+#[test]
+fn a_batch_declaring_two_different_priors_is_rejected() {
+    let mut h = history();
+    let err = h
+        .add_events(vec![
+            bout("a", "b", 0, Some(Gaussian::from_ms(40.0, 1.0)), None),
+            bout("a", "b", 1, Some(Gaussian::from_ms(10.0, 1.0)), None),
+        ])
+        .expect_err("two different priors for one competitor in one batch");
+
+    assert!(
+        matches!(
+            err,
+            InferenceError::ConflictingCompetitorConfig { field: "prior", .. }
+        ),
+        "got {err:?}"
+    );
+}
+
+/// A member setting only `drift_scale` must not also assert the default prior,
+/// or it would silently undo a prior seeded earlier. This is why the collected
+/// configuration tracks each field separately rather than a merged `Rating`.
+#[test]
+fn setting_one_field_late_leaves_the_other_alone() {
+    let seeded = Gaussian::from_ms(40.0, 1.0);
+
+    let mut h = history();
+    h.add_events(vec![bout("a", "b", 0, Some(seeded), None)])
+        .unwrap();
+    // Only the scale this time — the prior above must survive.
+    h.add_events(vec![bout("a", "b", 1, None, Some(0.5))])
+        .unwrap();
+    let _ = h.converge().unwrap();
+
+    let mut both_upfront = history();
+    both_upfront
+        .add_events(vec![
+            bout("a", "b", 0, Some(seeded), Some(0.5)),
+            bout("a", "b", 1, None, None),
+        ])
+        .unwrap();
+    let _ = both_upfront.converge().unwrap();
+
+    let (a, b) = (skill_of(&h, "a"), skill_of(&both_upfront, "a"));
+    assert!(
+        (a.mu() - b.mu()).abs() < 1e-9 && (a.sigma() - b.sigma()).abs() < 1e-9,
+        "setting drift_scale late clobbered the earlier prior: {a:?} vs {b:?}"
+    );
+}

tests/degenerate_inputs.rs

@@ -351,7 +351,7 @@ fn zero_weight_does_not_produce_a_non_finite_posterior() {
         .commit()
         .expect("a zero weight is accepted today; update this test if that changes");
 
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     assert_curve_finite(&h, &["a", "b"], "zero weight");
 }
@@ -368,7 +368,7 @@ fn negative_weight_does_not_produce_a_non_finite_posterior() {
         .commit()
         .expect("a negative weight is accepted today; update this test if that changes");
 
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     assert_curve_finite(&h, &["a", "b"], "negative weight");
 }
@@ -389,7 +389,7 @@ fn out_of_order_timestamps_converge_to_the_same_answer() {
             h.record_winner(&"a", &"b", time).unwrap();
         }
 
-        h.converge().unwrap();
+        let _ = h.converge().unwrap();
         h
     }
 
@@ -416,7 +416,7 @@ fn extreme_beta_and_sigma_stay_finite() {
 
         h.record_winner(&"a", &"b", 1).unwrap();
         h.record_winner(&"a", &"b", 2).unwrap();
-        h.converge().unwrap();
+        let _ = h.converge().unwrap();
 
         assert_curve_finite(&h, &["a", "b"], &format!("beta={beta} sigma={sigma}"));
     }

tests/determinism.rs

@@ -47,7 +47,7 @@ fn build_and_converge(seed: u64) -> Vec<(i64, trueskill_tt::Gaussian)> {
         });
     }
     h.add_events(events).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
     // Sample one competitor's curve for the comparison.
     h.learning_curve("p0")
 }

tests/drift_scale.rs

@@ -58,7 +58,7 @@ fn fit(events: Vec<Event<i64, &'static str>>, gamma: f64) -> Fit {
         .build();
 
     h.add_events(events).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
     h
 }
 
@@ -341,13 +341,20 @@ fn zero_scale_pins_a_competitor_in_the_filtered_pass() {
     );
 }
 
-/// `drift_scale` is competitor configuration captured at first appearance, the
-/// same as `prior` — a later `with_drift_scale` on a key the history already
-/// knows is ignored. This guards that decision rather than driving it: the
-/// behaviour falls out of where the capture happens, and the point of the test
-/// is that moving the capture would be a visible break, not a silent one.
+/// `drift_scale` is competitor configuration, and configuration supplied for a
+/// competitor the history already knows is now *applied* rather than dropped.
+///
+/// This test previously asserted the opposite. It was written as a deliberate
+/// change-detector — "moving the capture would be a visible break, not a silent
+/// one" — and that is exactly what happened: the capture moved, and the
+/// assertion inverted rather than being deleted.
+///
+/// Because configuration lives on the competitor and `converge` refits from
+/// competitor state, a late pin applies to the *whole* history, not just to
+/// events after it. So a scale set on the second batch must reach the same fit
+/// as one set from the very first event.
 #[test]
-fn drift_scale_is_ignored_after_first_appearance() {
+fn drift_scale_applies_when_set_after_first_appearance() {
     let mut late = History::builder()
         .mu(25.0)
         .sigma(25.0 / 3.0)
@@ -368,7 +375,7 @@ fn drift_scale_is_ignored_after_first_appearance() {
     }])
     .unwrap();
 
-    // Second batch asks for a pin. Too late: the competitor already exists.
+    // Second batch asks for a pin. No longer too late.
     late.add_events(vec![Event {
         time: 1000,
         teams: smallvec![
@@ -378,25 +385,115 @@ fn drift_scale_is_ignored_after_first_appearance() {
         outcome: Outcome::winner(1, 2),
     }])
     .unwrap();
-    late.converge().unwrap();
+    let _ = late.converge().unwrap();
 
-    let ignored = curve(&late, "anchor");
-    let drifting = curve(&fit(distant_pair(None), 25.0 / 300.0), "anchor");
+    let applied = curve(&late, "anchor");
+    let pinned_from_the_start = curve(&fit(distant_pair(Some(0.0)), 25.0 / 300.0), "anchor");
+    let never_pinned = curve(&fit(distant_pair(None), 25.0 / 300.0), "anchor");
 
-    for ((t_l, g_l), (t_r, g_r)) in ignored.iter().zip(drifting.iter()) {
+    for ((t_l, g_l), (t_r, g_r)) in applied.iter().zip(pinned_from_the_start.iter()) {
         assert_eq!(t_l, t_r);
         assert!(
             (g_l.sigma() - g_r.sigma()).abs() < 1e-9,
-            "a scale set after first appearance must be ignored, leaving the fit \
-             identical to one that never set it: t={t_l}, {} vs {}",
+            "a late pin should refit the whole history: t={t_l}, {} vs {}",
             g_l.sigma(),
             g_r.sigma()
         );
     }
 
-    let pinned = curve(&fit(distant_pair(Some(0.0)), 25.0 / 300.0), "anchor");
+    // And it must actually have done something.
     assert!(
-        (ignored[1].1.sigma() - pinned[1].1.sigma()).abs() > 1e-6,
-        "sanity: the pinned fit must actually differ, or the assertion above is vacuous"
+        applied
+            .iter()
+            .zip(never_pinned.iter())
+            .any(|((_, a), (_, b))| (a.sigma() - b.sigma()).abs() > 1e-9),
+        "the pin had no effect at all — the silent drop is back"
+    );
+}
+
+/// Re-declaring the same configuration must be inert. This is the shape a
+/// caller gets when the configuration is a property of the domain — "layouts
+/// are static" — so every ingestion path repeats it on every event.
+///
+/// Both histories see exactly the same events; only how many times the scale
+/// is declared differs.
+#[test]
+fn repeating_the_same_configuration_changes_nothing() {
+    let events = |declare_every_time: bool| {
+        let anchor = |first: bool| {
+            if first || declare_every_time {
+                Member::new("anchor").with_drift_scale(0.0)
+            } else {
+                Member::new("anchor")
+            }
+        };
+        vec![
+            Event {
+                time: 0,
+                teams: smallvec![
+                    Team::with_members([anchor(true)]),
+                    Team::with_members([Member::new("player")]),
+                ],
+                outcome: Outcome::winner(0, 2),
+            },
+            Event {
+                time: 1000,
+                teams: smallvec![
+                    Team::with_members([anchor(false)]),
+                    Team::with_members([Member::new("player")]),
+                ],
+                outcome: Outcome::winner(1, 2),
+            },
+        ]
+    };
+
+    let once = curve(&fit(events(false), 25.0 / 300.0), "anchor");
+    let every_time = curve(&fit(events(true), 25.0 / 300.0), "anchor");
+
+    for ((t_l, a), (t_r, b)) in once.iter().zip(every_time.iter()) {
+        assert_eq!(t_l, t_r);
+        assert!(
+            (a.sigma() - b.sigma()).abs() < 1e-12,
+            "t={t_l}: declaring the same scale repeatedly changed the fit, {} vs {}",
+            a.sigma(),
+            b.sigma()
+        );
+    }
+}
+
+#[test]
+fn a_batch_that_contradicts_itself_is_rejected() {
+    let mut h = History::builder().convergence(CONVERGENCE).build();
+
+    let err = h
+        .add_events(vec![
+            Event {
+                time: 0,
+                teams: smallvec![
+                    Team::with_members([Member::new("anchor").with_drift_scale(0.0)]),
+                    Team::with_members([Member::new("player")]),
+                ],
+                outcome: Outcome::winner(0, 2),
+            },
+            Event {
+                time: 1,
+                teams: smallvec![
+                    Team::with_members([Member::new("anchor").with_drift_scale(1.0)]),
+                    Team::with_members([Member::new("player")]),
+                ],
+                outcome: Outcome::winner(0, 2),
+            },
+        ])
+        .expect_err("two different scales for one competitor in one batch");
+
+    assert!(
+        matches!(
+            err,
+            InferenceError::ConflictingCompetitorConfig {
+                field: "drift_scale",
+                ..
+            }
+        ),
+        "got {err:?}"
     );
 }

tests/equivalence.rs

@@ -19,7 +19,8 @@ fn ts_rating(mu: f64, sigma: f64, beta: f64, gamma: f64) -> R {
 fn game_1v1_golden_matches_historical() {
     let a = ts_rating(25.0, 25.0 / 3.0, 25.0 / 6.0, 25.0 / 300.0);
     let b = ts_rating(25.0, 25.0 / 3.0, 25.0 / 6.0, 25.0 / 300.0);
-    let (a_post, b_post) = Game::<i64, _>::one_v_one(&a, &b, Outcome::winner(0, 2)).unwrap();
+    let (a_post, b_post) =
+        Game::<i64, _>::one_v_one(&a, &b, Outcome::winner(0, 2), &GameOptions::default()).unwrap();
     // Historical golden from pre-T2 test_1vs1 (team 0 wins):
     assert_ulps_eq!(
         a_post,

tests/filtered.rs

@@ -47,7 +47,7 @@ fn tight() -> ConvergenceOptions {
 fn filtered_evidence_sits_between_coin_flip_and_batch() {
     let mut history = repeated_winner(5);
 
-    history.converge().unwrap();
+    let _ = history.converge().unwrap();
 
     let coin_flip = 5.0 * 0.5f64.ln();
     let batch = history.log_evidence();
@@ -71,7 +71,7 @@ fn filtered_evidence_sits_between_coin_flip_and_batch() {
 fn filtered_first_point_is_less_certain_than_smoothed() {
     let mut history = repeated_winner(12);
 
-    history.converge().unwrap();
+    let _ = history.converge().unwrap();
 
     let smoothed = history.learning_curve("a");
     let filtered = history.filtered_learning_curve("a");
@@ -121,7 +121,7 @@ fn filtered_first_point_is_less_certain_than_smoothed() {
 fn filtered_curves_plural_agrees_with_singular() {
     let mut history = repeated_winner(4);
 
-    history.converge().unwrap();
+    let _ = history.converge().unwrap();
 
     let curves = history.filtered_learning_curves();
 
@@ -180,7 +180,7 @@ fn single_slice_filtered_matches_smoothed() {
         ])
         .unwrap();
 
-    history.converge().unwrap();
+    let _ = history.converge().unwrap();
 
     let smoothed = history.learning_curve("a");
     let filtered = history.filtered_learning_curve("a");
@@ -223,13 +223,13 @@ fn filtered_curves_do_not_depend_on_ingestion_order() {
 
     let mut batched = History::builder().convergence(tight()).build();
     batched.add_events(all.clone()).unwrap();
-    batched.converge().unwrap();
+    let _ = batched.converge().unwrap();
 
     let mut incremental = History::builder().convergence(tight()).build();
     for event in all {
         incremental.add_events([event]).unwrap();
     }
-    incremental.converge().unwrap();
+    let _ = incremental.converge().unwrap();
 
     let from_batched = batched.filtered_learning_curve("a");
     let from_incremental = incremental.filtered_learning_curve("a");

tests/game.rs

@@ -32,7 +32,8 @@ fn game_ranked_1v1_golden() {
 fn game_one_v_one_shortcut() {
     let a = default_rating();
     let b = default_rating();
-    let (a_post, b_post) = Game::<i64, _>::one_v_one(&a, &b, Outcome::winner(0, 2)).unwrap();
+    let (a_post, b_post) =
+        Game::<i64, _>::one_v_one(&a, &b, Outcome::winner(0, 2), &GameOptions::default()).unwrap();
     assert!(a_post.mu() > 25.0);
     assert!(b_post.mu() < 25.0);
 }
@@ -95,3 +96,45 @@ fn game_log_evidence_is_finite() {
     assert!(g.log_evidence().is_finite());
     assert!(g.log_evidence() < 0.0);
 }
+
+/// `one_v_one` used to hardcode `GameOptions::default()`, so a 1v1 could
+/// never set `p_draw` and a drawn 1v1 was unreachable through it.
+#[test]
+fn one_v_one_honours_the_draw_probability_it_is_given() {
+    let a = default_rating();
+    let b = default_rating();
+
+    // Default options still reject a draw, because the default p_draw is zero.
+    let err = Game::<i64, _>::one_v_one(&a, &b, Outcome::draw(2), &GameOptions::default())
+        .expect_err("a draw needs a positive p_draw");
+    assert!(matches!(
+        err,
+        InferenceError::TieWithoutDrawProbability { .. }
+    ));
+
+    // With a draw probability supplied it succeeds — which was impossible
+    // before the signature took options.
+    let options = GameOptions {
+        p_draw: 0.25,
+        ..GameOptions::default()
+    };
+    let (a_post, b_post) = Game::<i64, _>::one_v_one(&a, &b, Outcome::draw(2), &options)
+        .expect("a draw is representable once p_draw is positive");
+
+    // A symmetric draw leaves the means alone and sharpens both sides.
+    assert!((a_post.mu() - b_post.mu()).abs() < 1e-9);
+    assert!(a_post.sigma() < 25.0 / 3.0);
+}
+
+/// Convergence options reach the 1v1 path too, not just `p_draw`.
+#[test]
+fn one_v_one_honours_convergence_options() {
+    let a = default_rating();
+    let b = default_rating();
+    let options = GameOptions {
+        convergence: ConvergenceOptions::default(),
+        ..GameOptions::default()
+    };
+    let (a_post, _) = Game::<i64, _>::one_v_one(&a, &b, Outcome::winner(0, 2), &options).unwrap();
+    assert!(a_post.mu() > 25.0);
+}

tests/ingestion_equivalence.rs

@@ -30,6 +30,22 @@ fn event(a: &str, b: &str, time: i64) -> Event<i64, String> {
     }
 }
 
+/// Like [`event`], but `a` carries competitor configuration.
+///
+/// `prior` and `drift_scale` configure the competitor rather than the event, so
+/// they are the part of ingestion most exposed to order: they are consumed once,
+/// where the competitor's state is written.
+fn configured_event(a: &str, b: &str, time: i64, scale: f64) -> Event<i64, String> {
+    Event {
+        time,
+        teams: smallvec![
+            Team::with_members([Member::new(a.to_string()).with_drift_scale(scale)]),
+            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();
@@ -145,3 +161,65 @@ fn back_dated_event_matches_batched() {
     let incremental = converged_skills(events, false);
     assert_same(&batched, &incremental, "back-dated event");
 }
+
+/// The invariant this file protects was only ever checked for *unconfigured*
+/// competitors — every helper above built members with `Member::new`.
+///
+/// Configuration is the part most exposed to ordering, because it is consumed
+/// once at the point the competitor's state is written rather than replayed per
+/// event. These cover it.
+#[test]
+fn configured_competitors_are_order_independent() {
+    let events = vec![
+        configured_event("a", "b", 0, 0.0),
+        configured_event("a", "c", 1, 0.0),
+        configured_event("a", "b", 2, 0.0),
+        event("b", "c", 3),
+    ];
+
+    assert_same(
+        &converged_skills(events.clone(), true),
+        &converged_skills(events, false),
+        "configuration repeated on every appearance",
+    );
+}
+
+/// Configuration supplied only on a *later* event is the case that used to be
+/// silently dropped. It must now reach the same fit either way it is ingested.
+#[test]
+fn late_configuration_is_order_independent() {
+    let events = vec![
+        event("a", "b", 0),
+        configured_event("a", "c", 1, 0.0),
+        event("a", "b", 2),
+    ];
+
+    assert_same(
+        &converged_skills(events.clone(), true),
+        &converged_skills(events, false),
+        "configuration supplied after first appearance",
+    );
+}
+
+/// And it must actually be doing something — an implementation that dropped
+/// configuration entirely would pass both tests above.
+#[test]
+fn configuration_changes_the_fit_however_it_is_ingested() {
+    let configured = vec![
+        event("a", "b", 0),
+        configured_event("a", "c", 1, 0.0),
+        event("a", "b", 2),
+    ];
+    let plain = vec![event("a", "b", 0), event("a", "c", 1), event("a", "b", 2)];
+
+    for batched in [true, false] {
+        let with = converged_skills(configured.clone(), batched);
+        let without = converged_skills(plain.clone(), batched);
+        assert!(
+            with.iter()
+                .zip(&without)
+                .any(|((_, x), (_, y))| (x.sigma() - y.sigma()).abs() > 1e-9),
+            "batched={batched}: configuration had no effect, so the order tests are vacuous"
+        );
+    }
+}

tests/large_history_converges_finite.rs

@@ -46,7 +46,7 @@ fn nan_after_fit(players: usize) -> usize {
         let (w, l) = if rng.coin() { (a, b) } else { (b, a) };
         h.record_winner(&ids[w], &ids[l], 0).unwrap();
     }
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     ids.iter()
         .filter(|id| {

tests/marginal_calibration.rs

@@ -0,0 +1,367 @@
+//! Calibration of the crate's marginals against the EXACT posterior.
+//!
+//! A scored history is linear-Gaussian — `MarginFactor` encodes
+//! `score_a - score_b ~ N(perf_a - perf_b, score_sigma^2)` — so the true joint
+//! posterior has a closed form and the crate can be checked against ground
+//! truth rather than against intuition. That is not possible for ranked
+//! outcomes, whose truncation likelihood EP genuinely approximates.
+//!
+//! Two things are pinned here, and one is deliberately only recorded.
+//!
+//! **Pinned: on a tree the crate is exact**, means and variances both. Message
+//! passing has no approximation to make when the factor graph has no cycles, so
+//! any drift here would be a real defect.
+//!
+//! **Pinned: means are exact even with cycles.** This is the standard result
+//! for Gaussian belief propagation (Weiss & Freeman 2001) and it is what makes
+//! ratings trustworthy.
+//!
+//! **Recorded, not asserted: with cycles, marginal variances are too narrow.**
+//! Measured on the round-robin fixture below, the crate reports sigma 1.430
+//! where the exact posterior is 2.851 — a ratio of 0.502. That is the known
+//! behaviour of loopy Gaussian BP, not a bug in this crate, and it is left
+//! unasserted because fixing it is exactly what #46 proposes.
+//!
+//! Why that matters for a consumer, and why #46 cannot be implemented as "add
+//! a covariance accessor": the exact correlation between two nodes here is
+//! +0.857, so a consumer computing `sqrt(sa^2 + sb^2)` for a difference
+//! overstates its width. But the too-narrow marginals partially cancel that,
+//! leaving 1.327x rather than 2.646x. Adding true correlations to these
+//! marginals without also correcting them would give 0.765 against a true
+//! 1.524 — *overconfident*, which is the unsafe direction.
+
+use smallvec::smallvec;
+use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team};
+
+const N: usize = 5;
+const MU0: f64 = 0.0;
+const SIGMA0: f64 = 6.0;
+const BETA: f64 = 1.0;
+const SCORE_SIGMA: f64 = 2.0;
+
+/// A STAR: every event touches c0, so the node-event graph is a tree and
+/// Gaussian BP is exact. Any discrepancy here is not caused by loops.
+fn tree_fixture() -> Vec<(usize, usize, f64)> {
+    vec![(0, 1, 3.0), (0, 2, 5.0), (0, 3, 4.0), (0, 4, 6.0)]
+}
+
+/// (winner, loser, score_diff)
+fn fixture() -> Vec<(usize, usize, f64)> {
+    vec![
+        (0, 1, 3.0),
+        (0, 2, 5.0),
+        (1, 2, 2.0),
+        (3, 4, 1.0),
+        (0, 3, 4.0),
+        (1, 4, 2.5),
+        (2, 3, 0.5),
+        (0, 4, 6.0),
+        (1, 3, 1.5),
+        (2, 4, 3.0),
+    ]
+}
+
+/// Invert a small symmetric positive-definite matrix by Gauss-Jordan.
+fn inverse(mut a: Vec<Vec<f64>>) -> Vec<Vec<f64>> {
+    let n = a.len();
+    let mut inv: Vec<Vec<f64>> = (0..n)
+        .map(|i| (0..n).map(|j| if i == j { 1.0 } else { 0.0 }).collect())
+        .collect();
+    for col in 0..n {
+        // partial pivot
+        let mut piv = col;
+        for r in col + 1..n {
+            if a[r][col].abs() > a[piv][col].abs() {
+                piv = r;
+            }
+        }
+        a.swap(col, piv);
+        inv.swap(col, piv);
+        let d = a[col][col];
+        for j in 0..n {
+            a[col][j] /= d;
+            inv[col][j] /= d;
+        }
+        for r in 0..n {
+            if r == col {
+                continue;
+            }
+            let f = a[r][col];
+            for j in 0..n {
+                a[r][j] -= f * a[col][j];
+                inv[r][j] -= f * inv[col][j];
+            }
+        }
+    }
+    inv
+}
+
+/// The exact posterior of a linear-Gaussian model:
+/// precision = prior precision + sum of a_k a_k^T / v_k.
+fn exact_for(obs: &[(usize, usize, f64)]) -> (Vec<f64>, Vec<Vec<f64>>) {
+    let mut lambda = vec![vec![0.0; N]; N];
+    let mut eta = [0.0; N];
+    for (i, row) in lambda.iter_mut().enumerate() {
+        row[i] = 1.0 / (SIGMA0 * SIGMA0);
+        eta[i] = MU0 / (SIGMA0 * SIGMA0);
+    }
+
+    // Each 1v1 observation: d ~ N(x_a - x_b, score_sigma^2 + 2 beta^2)
+    let v = SCORE_SIGMA * SCORE_SIGMA + 2.0 * BETA * BETA;
+    for &(a, b, d) in obs {
+        let mut vec_a = [0.0; N];
+        vec_a[a] = 1.0;
+        vec_a[b] = -1.0;
+        for i in 0..N {
+            for j in 0..N {
+                lambda[i][j] += vec_a[i] * vec_a[j] / v;
+            }
+            eta[i] += vec_a[i] * d / v;
+        }
+    }
+
+    let cov = inverse(lambda);
+    let mean: Vec<f64> = (0..N)
+        .map(|i| (0..N).map(|j| cov[i][j] * eta[j]).sum())
+        .collect();
+    (mean, cov)
+}
+
+fn key(i: usize) -> &'static str {
+    ["c0", "c1", "c2", "c3", "c4"][i]
+}
+
+/// Returns (worst mean error, worst sd ratio).
+fn fitted(
+    obs: &[(usize, usize, f64)],
+) -> History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str> {
+    let mut h: History<i64, _, _, &'static str> = History::builder()
+        .mu(MU0)
+        .sigma(SIGMA0)
+        .beta(BETA)
+        .score_sigma(SCORE_SIGMA)
+        .drift(ConstantDrift(0.0))
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build();
+
+    let events: Vec<Event<i64, &'static str>> = obs
+        .iter()
+        .copied()
+        .map(|(a, b, d)| Event {
+            time: 1,
+            teams: smallvec![
+                Team::with_members([Member::new(key(a))]),
+                Team::with_members([Member::new(key(b))]),
+            ],
+            outcome: Outcome::scores([d, 0.0]),
+        })
+        .collect();
+    h.add_events(events).unwrap();
+    let report = h.converge().unwrap();
+    assert!(
+        report.converged,
+        "fixture must converge: {:?}",
+        report.final_step
+    );
+
+    h
+}
+
+/// Returns (worst mean error, worst sd ratio gap).
+fn run(name: &str, obs: Vec<(usize, usize, f64)>) -> (f64, f64) {
+    println!("\n########## {name} ##########");
+    let h = fitted(&obs);
+    let (mean, cov) = exact_for(&obs);
+
+    println!("\n== marginals: crate vs the exact linear-Gaussian posterior ==");
+    println!(
+        "{:>4} {:>12} {:>12} {:>12} {:>12} {:>8}",
+        "node", "crate mu", "exact mu", "crate sd", "exact sd", "sd ratio"
+    );
+    for i in 0..N {
+        let g = h.current_skill(&key(i)).unwrap();
+        let exact_sd = cov[i][i].sqrt();
+        println!(
+            "{:>4} {:>12.6} {:>12.6} {:>12.6} {:>12.6} {:>8.3}",
+            key(i),
+            g.mu(),
+            mean[i],
+            g.sigma(),
+            exact_sd,
+            g.sigma() / exact_sd
+        );
+    }
+
+    let mut worst_mean = 0.0f64;
+    let mut worst_ratio_gap = 0.0f64;
+    for i in 0..N {
+        let g = h.current_skill(&key(i)).unwrap();
+        worst_mean = worst_mean.max((g.mu() - mean[i]).abs());
+        worst_ratio_gap = worst_ratio_gap.max((g.sigma() / cov[i][i].sqrt() - 1.0).abs());
+    }
+
+    println!("\n== what a consumer actually computes for a DIFFERENCE ==");
+    println!(
+        "{:>8} {:>12} {:>14} {:>14} {:>12}",
+        "pair", "exact", "naive(exact)", "naive(crate)", "crate err"
+    );
+    for i in 0..N {
+        for j in i + 1..N {
+            if i != 0 && j != 1 {
+                continue;
+            }
+            let gi = h.current_skill(&key(i)).unwrap();
+            let gj = h.current_skill(&key(j)).unwrap();
+            let exact_sd = (cov[i][i] + cov[j][j] - 2.0 * cov[i][j]).sqrt();
+            let naive_exact = (cov[i][i] + cov[j][j]).sqrt();
+            let naive_crate = (gi.sigma().powi(2) + gj.sigma().powi(2)).sqrt();
+            let corr = cov[i][j] / (cov[i][i].sqrt() * cov[j][j].sqrt());
+            println!(
+                "{:>8} {:>12.6} {:>14.6} {:>14.6} {:>11.3}x   (corr {corr:.4})",
+                format!("{}-{}", key(i), key(j)),
+                exact_sd,
+                naive_exact,
+                naive_crate,
+                naive_crate / exact_sd
+            );
+        }
+    }
+
+    (worst_mean, worst_ratio_gap)
+}
+
+/// With no cycles there is nothing for message passing to approximate.
+#[test]
+fn on_a_tree_the_marginals_are_exact() {
+    let (mean_err, sd_gap) = run("TREE (star: no loops, BP is exact)", tree_fixture());
+    assert!(
+        mean_err < 1e-9,
+        "tree means should be exact, worst error {mean_err}"
+    );
+    assert!(
+        sd_gap < 1e-9,
+        "tree sigmas should be exact, worst ratio gap {sd_gap}"
+    );
+}
+
+/// With cycles the means stay exact — the property ratings depend on — while
+/// the variances do not. The variance gap is measured and reported rather than
+/// asserted; see the module docs.
+#[test]
+fn with_cycles_the_means_stay_exact_but_the_variances_shrink() {
+    let (mean_err, sd_gap) = run("LOOPY (round robin)", fixture());
+    assert!(
+        mean_err < 1e-9,
+        "loopy means must still be exact, worst error {mean_err}"
+    );
+    assert!(
+        sd_gap > 0.1,
+        "the loopy variance gap is the premise of #46; if it has closed, that \
+         issue and these docs need revisiting (worst ratio gap {sd_gap})"
+    );
+}
+
+/// The point of #46: `posterior_of` must reproduce the exact joint, including
+/// the correlation that marginals cannot express.
+#[test]
+fn posterior_of_matches_the_exact_joint() {
+    for (name, obs) in [("tree", tree_fixture()), ("loopy", fixture())] {
+        let h = fitted(&obs);
+        let (_, cov) = exact_for(&obs);
+
+        println!("\n== posterior_of vs exact ({name}) ==");
+        println!(
+            "{:>12} {:>14} {:>14} {:>10}",
+            "functional", "posterior_of", "exact", "ratio"
+        );
+
+        for (i, j) in [(0usize, 1usize), (0, 2), (1, 3), (2, 4)] {
+            let got = h
+                .posterior_of(&[(&key(i), 1.0), (&key(j), -1.0)])
+                .expect("scored slice should have a joint");
+            let exact_sd = (cov[i][i] + cov[j][j] - 2.0 * cov[i][j]).sqrt();
+            println!(
+                "{:>12} {:>14.6} {:>14.6} {:>10.4}",
+                format!("{}-{}", key(i), key(j)),
+                got.sigma(),
+                exact_sd,
+                got.sigma() / exact_sd
+            );
+            assert!(
+                (got.sigma() - exact_sd).abs() / exact_sd < 1e-9,
+                "{name} {}-{}: posterior_of gave {} where the exact joint is {exact_sd}",
+                key(i),
+                key(j),
+                got.sigma()
+            );
+        }
+
+        // A single competitor: this is where the loopy marginal was 2x narrow.
+        for (i, row) in cov.iter().enumerate() {
+            let got = h.posterior_of(&[(&key(i), 1.0)]).unwrap();
+            let exact_sd = row[i].sqrt();
+            assert!(
+                (got.sigma() - exact_sd).abs() / exact_sd < 1e-9,
+                "{name} {}: posterior_of gave {} where exact is {exact_sd}",
+                key(i),
+                got.sigma()
+            );
+        }
+        println!("  single-competitor marginals also exact");
+    }
+}
+
+/// Cost of the dense solve as the slice grows. Recorded, not asserted.
+#[test]
+#[ignore = "timing probe, run explicitly"]
+fn cost_scaling() {
+    use std::time::Instant;
+    for n in [50usize, 100, 200, 400, 800] {
+        let names: Vec<String> = (0..n).map(|i| format!("c{i}")).collect();
+        let mut h: History<i64, _, _, String> = History::builder_with_key()
+            .score_sigma(2.0)
+            .drift(ConstantDrift(0.0))
+            .convergence(ConvergenceOptions {
+                max_iter: 200,
+                epsilon: 1e-8,
+                alpha: 1.0,
+            })
+            .build();
+        let mut seed = 5u64;
+        let mut rnd = move || {
+            seed ^= seed << 13;
+            seed ^= seed >> 7;
+            seed ^= seed << 17;
+            seed
+        };
+        let events: Vec<Event<i64, String>> = (0..n * 4)
+            .map(|_| {
+                let a = (rnd() as usize) % n;
+                let mut b = (rnd() as usize) % n;
+                if b == a {
+                    b = (b + 1) % n;
+                }
+                Event {
+                    time: 1,
+                    teams: smallvec![
+                        Team::with_members([Member::new(names[a].clone())]),
+                        Team::with_members([Member::new(names[b].clone())]),
+                    ],
+                    outcome: Outcome::scores([1.0, 0.0]),
+                }
+            })
+            .collect();
+        h.add_events(events).unwrap();
+        let _ = h.converge().unwrap();
+
+        let t = Instant::now();
+        let g = h
+            .posterior_of(&[(&names[0], 1.0), (&names[1], -1.0)])
+            .unwrap();
+        println!("  n={n:>4}: {:>10.2?}  sigma {:.6}", t.elapsed(), g.sigma());
+    }
+}

tests/observer.rs

@@ -0,0 +1,161 @@
+//! `Observer` callbacks must actually fire.
+//!
+//! `on_slice_processed` (formerly `on_batch_processed`) was declared on the
+//! trait and never called from anywhere, so implementors wired up a callback
+//! that could not run. These tests exist so that cannot silently recur.
+
+use std::sync::{Arc, Mutex};
+
+use trueskill_tt::{History, Observer};
+
+/// Plain fields. `Arc<O>` implements `Observer`, so the caller shares the
+/// observer itself rather than wrapping each field in its own `Arc`.
+#[derive(Default)]
+struct Recorder {
+    iterations: Mutex<Vec<usize>>,
+    slices: Mutex<Vec<(i64, usize, usize)>>,
+    converged: Mutex<Vec<(usize, bool)>>,
+}
+
+impl Observer<i64> for Recorder {
+    fn on_iteration_end(&self, iter: usize, _max_step: (f64, f64)) {
+        self.iterations.lock().unwrap().push(iter);
+    }
+
+    fn on_slice_processed(&self, time: &i64, slice_idx: usize, n_events: usize) {
+        self.slices
+            .lock()
+            .unwrap()
+            .push((*time, slice_idx, n_events));
+    }
+
+    fn on_converged(&self, iters: usize, _final_step: (f64, f64), converged: bool) {
+        self.converged.lock().unwrap().push((iters, converged));
+    }
+}
+
+#[test]
+fn every_observer_callback_fires() {
+    let recorder = Arc::new(Recorder::default());
+    let mut h = History::builder().observer(Arc::clone(&recorder)).build();
+
+    h.record_winner(&"a", &"b", 1).unwrap();
+    h.record_winner(&"b", &"c", 2).unwrap();
+    h.record_winner(&"c", &"a", 3).unwrap();
+    let _ = h.converge().unwrap();
+
+    assert!(
+        !recorder.iterations.lock().unwrap().is_empty(),
+        "on_iteration_end never fired"
+    );
+    assert!(
+        !recorder.converged.lock().unwrap().is_empty(),
+        "on_converged never fired"
+    );
+    assert!(
+        !recorder.slices.lock().unwrap().is_empty(),
+        "on_slice_processed never fired — the defect this test exists for"
+    );
+}
+
+#[test]
+fn slice_callbacks_report_the_slice_they_swept() {
+    let recorder = Arc::new(Recorder::default());
+    let mut h = History::builder().observer(Arc::clone(&recorder)).build();
+
+    h.record_winner(&"a", &"b", 10).unwrap();
+    h.record_winner(&"a", &"b", 20).unwrap();
+    let _ = h.converge().unwrap();
+
+    let slices = recorder.slices.lock().unwrap();
+
+    // Only the times actually in the history, and each with its own events.
+    for &(time, idx, events) in slices.iter() {
+        assert!(time == 10 || time == 20, "unexpected slice time {time}");
+        assert!(idx < 2, "slice index {idx} out of range");
+        assert_eq!(events, 1, "each slice holds exactly one event");
+    }
+
+    // Both slices must be reported, not just one end of the sweep.
+    assert!(
+        slices.iter().any(|&(t, ..)| t == 10),
+        "slice 10 never reported"
+    );
+    assert!(
+        slices.iter().any(|&(t, ..)| t == 20),
+        "slice 20 never reported"
+    );
+}
+
+#[test]
+fn a_single_slice_history_still_reports_its_sweep() {
+    let recorder = Arc::new(Recorder::default());
+    let mut h = History::builder().observer(Arc::clone(&recorder)).build();
+
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let _ = h.converge().unwrap();
+
+    let slices = recorder.slices.lock().unwrap();
+    assert!(
+        !slices.is_empty(),
+        "the single-slice path must report its sweep too"
+    );
+    assert!(slices.iter().all(|&(t, idx, _)| t == 1 && idx == 0));
+}
+
+/// The gap #40 closed: without `impl Observer for Arc<O>`, an observer that
+/// accumulates anything had to wrap every field in its own `Arc` and derive
+/// `Clone`, because `History` consumes the observer and never hands it back.
+#[test]
+fn a_shared_observer_reaches_the_callers_handle() {
+    let recorder = Arc::new(Recorder::default());
+    let mut h = History::builder().observer(Arc::clone(&recorder)).build();
+
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let _ = h.converge().unwrap();
+
+    assert!(!recorder.iterations.lock().unwrap().is_empty());
+    assert!(!recorder.slices.lock().unwrap().is_empty());
+    assert!(!recorder.converged.lock().unwrap().is_empty());
+}
+
+/// `?Sized` on the blanket impls means the observer can be chosen at runtime.
+#[test]
+fn a_trait_object_observer_works() {
+    let boxed: Box<dyn Observer<i64>> = Box::new(Recorder::default());
+    let mut h = History::builder().observer(boxed).build();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let _ = h.converge().unwrap();
+
+    let shared: Arc<dyn Observer<i64>> = Arc::new(Recorder::default());
+    let mut h = History::builder().observer(Arc::clone(&shared)).build();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let _ = h.converge().unwrap();
+}
+
+/// A non-shared observer can be reclaimed after convergence instead.
+#[test]
+fn into_observer_returns_the_accumulated_state() {
+    let mut h = History::builder().observer(Recorder::default()).build();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    let _ = h.converge().unwrap();
+
+    // Readable in place...
+    assert!(!h.observer().iterations.lock().unwrap().is_empty());
+
+    // ...and reclaimable by value.
+    let recorder = h.into_observer();
+    assert!(!recorder.slices.lock().unwrap().is_empty());
+}
+
+/// Borrowing works too, for an observer that outlives the history.
+#[test]
+fn a_borrowed_observer_works() {
+    let recorder = Recorder::default();
+    {
+        let mut h = History::builder().observer(&recorder).build();
+        h.record_winner(&"a", &"b", 1).unwrap();
+        let _ = h.converge().unwrap();
+    }
+    assert!(!recorder.iterations.lock().unwrap().is_empty());
+}

tests/predict_margin.rs

@@ -0,0 +1,153 @@
+//! `predict_margin`: the predictive distribution of a scored matchup.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, InferenceError, Member, Outcome, Team,
+    UnknownKeys,
+};
+
+fn builder(
+    policy: UnknownKeys,
+) -> History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str> {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.0))
+        .unknown_keys(policy)
+        .convergence(ConvergenceOptions {
+            max_iter: 5_000,
+            epsilon: 1e-12,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+fn round(a: &'static str, b: &'static str, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: 1,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+/// A history where "veteran" and "regular" are well observed and "novice"
+/// appears once.
+fn fitted(
+    policy: UnknownKeys,
+) -> History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str> {
+    let mut h = builder(policy);
+    let mut events: Vec<_> = (0..40)
+        .map(|t| round("veteran", "regular", 10.0 + f64::from(t % 3), 5.0))
+        .collect();
+    events.push(round("veteran", "novice", 10.0, 6.0));
+    h.add_events(events).unwrap();
+    let _ = h.converge().unwrap();
+    h
+}
+
+/// The property #48 exists for: the interval must widen when the model knows
+/// less. Their hand-fitted noise law quoted the same sigma for a competitor
+/// with forty rounds and one with none.
+#[test]
+fn the_interval_widens_as_the_model_knows_less() {
+    let h = fitted(UnknownKeys::Prior);
+
+    let well_known = h
+        .predict_margin(&[&[&"veteran"], &[&"regular"]])
+        .unwrap()
+        .sigma();
+    let thin = h
+        .predict_margin(&[&[&"veteran"], &[&"novice"]])
+        .unwrap()
+        .sigma();
+    let unseen = h
+        .predict_margin(&[&[&"veteran"], &[&"stranger"]])
+        .unwrap()
+        .sigma();
+
+    assert!(
+        well_known < thin && thin < unseen,
+        "margin width should grow as evidence thins: {well_known} < {thin} < {unseen}"
+    );
+}
+
+/// #48's second requirement: an unseen competitor is a legitimate question, not
+/// an error, and the answer should come from the prior rather than be faked.
+#[test]
+fn an_unseen_competitor_is_answered_from_the_prior() {
+    let h = fitted(UnknownKeys::Prior);
+    let g = h.predict_margin(&[&[&"nobody"], &[&"no_one"]]).unwrap();
+
+    // Two unknowns: the gap is centred on zero and carries both priors plus
+    // both performance noises plus the observation noise.
+    assert!(g.mu().abs() < 1e-9, "mu {}", g.mu());
+    let expected = (2.0 * 36.0 + 2.0 * 1.0 + 4.0f64).sqrt();
+    assert!(
+        (g.sigma() - expected).abs() < 1e-9,
+        "sigma {} vs expected {expected}",
+        g.sigma()
+    );
+}
+
+#[test]
+fn reject_still_rejects() {
+    let h = fitted(UnknownKeys::Reject);
+    assert!(matches!(
+        h.predict_margin(&[&[&"veteran"], &[&"stranger"]]),
+        Err(InferenceError::UnknownKey { .. })
+    ));
+}
+
+/// The margin is the *difference*, so it must be antisymmetric in the teams.
+#[test]
+fn swapping_the_teams_negates_the_margin() {
+    let h = fitted(UnknownKeys::Prior);
+    let forward = h.predict_margin(&[&[&"veteran"], &[&"regular"]]).unwrap();
+    let reverse = h.predict_margin(&[&[&"regular"], &[&"veteran"]]).unwrap();
+
+    assert!((forward.mu() + reverse.mu()).abs() < 1e-9);
+    assert!((forward.sigma() - reverse.sigma()).abs() < 1e-12);
+}
+
+/// The predictive interval must be wider than the skill gap alone: it also
+/// carries per-event performance noise and the observation noise.
+#[test]
+fn the_predictive_interval_exceeds_the_skill_uncertainty() {
+    let h = fitted(UnknownKeys::Prior);
+    let skill_gap = h
+        .posterior_of(&[(&"veteran", 1.0), (&"regular", -1.0)])
+        .unwrap();
+    let predictive = h.predict_margin(&[&[&"veteran"], &[&"regular"]]).unwrap();
+
+    assert!(
+        (predictive.mu() - skill_gap.mu()).abs() < 1e-12,
+        "means agree"
+    );
+    // beta^2 twice plus score_sigma^2 = 2 + 4.
+    let expected = (skill_gap.sigma().powi(2) + 6.0).sqrt();
+    assert!((predictive.sigma() - expected).abs() < 1e-12);
+    assert!(predictive.sigma() > skill_gap.sigma());
+}
+
+#[test]
+fn shape_errors_are_reported() {
+    let h = fitted(UnknownKeys::Prior);
+    assert!(matches!(
+        h.predict_margin(&[&[&"veteran"]]),
+        Err(InferenceError::MismatchedShape {
+            expected: 2,
+            got: 1,
+            ..
+        })
+    ));
+    let empty: [&&str; 0] = [];
+    assert!(matches!(
+        h.predict_margin(&[&[&"veteran"], &empty]),
+        Err(InferenceError::EmptyTeam { team: 1 })
+    ));
+}

tests/prediction.rs

@@ -0,0 +1,417 @@
+//! Prediction API: N-team outcomes, draw mass, and the error paths that used
+//! to be panics or silent wrong answers.
+
+use trueskill_tt::{History, InferenceError, MAX_PREDICTED_TEAMS};
+
+fn history_with(names: &[&'static str], p_draw: f64) -> History {
+    let mut h = History::builder().p_draw(p_draw).build();
+    // Give every competitor a recorded skill by playing a small round robin.
+    for pair in names.windows(2) {
+        h.record_winner(&pair[0], &pair[1], 1).unwrap();
+    }
+    let _ = h.converge().unwrap();
+    h
+}
+
+#[test]
+fn unknown_keys_are_reported_not_silently_dropped() {
+    let h = history_with(&["a", "b"], 0.0);
+
+    let err = h
+        .predict_outcome(&[&[&"a"], &[&"ghost"]])
+        .expect_err("an unknown key must not yield a confident prediction");
+    assert_eq!(
+        err,
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"ghost\"".to_owned(),
+        }
+    );
+
+    // Every prediction entry point, not just one.
+    assert!(
+        h.predict_win_probabilities(&[&[&"a"], &[&"ghost"]])
+            .is_err()
+    );
+    assert!(h.predict_quality(&[&[&"a"], &[&"ghost"]]).is_err());
+    assert!(h.predict_ranking(&[&[&"a"], &[&"ghost"]], &[0, 1]).is_err());
+}
+
+#[test]
+fn an_entirely_unknown_team_is_an_error() {
+    let h = history_with(&["a", "b"], 0.0);
+    let err = h.predict_outcome(&[&[&"a"], &[&"x", &"y"]]).unwrap_err();
+    assert_eq!(
+        err,
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"x\"".to_owned(),
+        }
+    );
+}
+
+#[test]
+fn degenerate_team_shapes_are_errors_rather_than_panics() {
+    let h = history_with(&["a", "b"], 0.0);
+
+    assert_eq!(
+        h.predict_outcome(&[&[&"a"]]).unwrap_err(),
+        InferenceError::NotEnoughTeams { got: 1 }
+    );
+    assert_eq!(
+        h.predict_outcome(&[]).unwrap_err(),
+        InferenceError::NotEnoughTeams { got: 0 }
+    );
+    assert_eq!(
+        h.predict_outcome(&[&[&"a"], &[]]).unwrap_err(),
+        InferenceError::EmptyTeam { team: 1 }
+    );
+}
+
+#[test]
+fn more_than_two_teams_no_longer_panics() {
+    let h = history_with(&["a", "b", "c"], 0.0);
+    let p = h
+        .predict_outcome(&[&[&"a"], &[&"b"], &[&"c"]])
+        .expect("three teams must be supported");
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+    // Three teams, no draws possible: exactly the six strict orderings.
+    assert_eq!(p.outcomes().len(), 6);
+}
+
+#[test]
+fn the_outcome_space_is_capped_rather_than_hanging() {
+    let names: Vec<&'static str> = vec!["a", "b", "c", "d", "e", "f", "g", "h"];
+    let h = history_with(&names, 0.0);
+
+    let teams: Vec<&[&&'static str]> = Vec::new();
+    let _ = teams;
+
+    let too_many: Vec<Vec<&&str>> = names.iter().map(|n| vec![n]).collect();
+    let refs: Vec<&[&&str]> = too_many.iter().map(Vec::as_slice).collect();
+
+    let err = h.predict_outcome(&refs).unwrap_err();
+    assert_eq!(
+        err,
+        InferenceError::TooManyTeams {
+            got: 8,
+            max: MAX_PREDICTED_TEAMS
+        }
+    );
+
+    // The cheap paths stay available at any size.
+    let wins = h.predict_win_probabilities(&refs).unwrap();
+    assert_eq!(wins.len(), 8);
+    assert!(
+        (wins.iter().sum::<f64>() - 1.0).abs() < 1e-6,
+        "win probabilities must still sum to one: {wins:?}"
+    );
+}
+
+/// The defect that made every draw-enabled prediction wrong: `[p, 1 - p]`
+/// allocated no mass to a draw even with `p_draw > 0`.
+#[test]
+fn a_draw_carries_probability_mass_when_p_draw_is_positive() {
+    let h = history_with(&["a", "b"], 0.25);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+
+    let draw = p.probability_of(&[0, 0]);
+    assert!(draw > 0.0, "a draw-enabled model must give draws mass");
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+
+    let wins = p.win_probabilities();
+    assert!(
+        (wins.iter().sum::<f64>() + draw - 1.0).abs() < 1e-6,
+        "wins {wins:?} plus draw {draw} must be the whole space"
+    );
+    assert!(
+        (p.shared_first_place() - draw).abs() < 1e-12,
+        "a two-team draw is a shared first place"
+    );
+}
+
+#[test]
+fn a_zero_draw_probability_admits_no_ties() {
+    let h = history_with(&["a", "b"], 0.0);
+    let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+    assert_eq!(p.probability_of(&[0, 0]), 0.0);
+    assert!(p.shared_first_place() < 1e-12);
+}
+
+/// The two routes to a win probability run through entirely different
+/// algorithms — adaptive quadrature versus the enumerated chain recursion —
+/// so agreement between them is a real cross-check, not a tautology.
+#[test]
+fn the_cheap_and_exhaustive_paths_agree() {
+    for p_draw in [0.0, 0.1] {
+        let h = history_with(&["a", "b", "c"], p_draw);
+        let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]];
+
+        let cheap = h.predict_win_probabilities(teams).unwrap();
+        let exhaustive = h.predict_outcome(teams).unwrap().win_probabilities();
+
+        for (i, (a, b)) in cheap.iter().zip(&exhaustive).enumerate() {
+            assert!(
+                (a - b).abs() < 1e-6,
+                "p_draw={p_draw} team {i}: quadrature {a} vs enumeration {b}"
+            );
+        }
+    }
+}
+
+#[test]
+fn predict_ranking_agrees_with_the_distribution() {
+    let h = history_with(&["a", "b", "c"], 0.1);
+    let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]];
+    let dist = h.predict_outcome(teams).unwrap();
+
+    for (ranks, expected) in dist.outcomes() {
+        let direct = h.predict_ranking(teams, ranks).unwrap();
+        assert!(
+            (direct - expected).abs() < 1e-9,
+            "ranks {ranks:?}: {direct} vs {expected}"
+        );
+    }
+}
+
+#[test]
+fn predict_ranking_checks_its_shape() {
+    let h = history_with(&["a", "b"], 0.0);
+    let err = h
+        .predict_ranking(&[&[&"a"], &[&"b"]], &[0, 1, 2])
+        .unwrap_err();
+    assert!(matches!(
+        err,
+        InferenceError::MismatchedShape {
+            expected: 2,
+            got: 3,
+            ..
+        }
+    ));
+}
+
+#[test]
+fn the_stronger_competitor_is_favoured() {
+    let mut h = History::builder().build();
+    for t in 1..=10 {
+        h.record_winner(&"strong", &"weak", t).unwrap();
+    }
+    let _ = h.converge().unwrap();
+
+    let p = h.predict_outcome(&[&[&"strong"], &[&"weak"]]).unwrap();
+    let (best, _) = p.most_likely().expect("a most likely outcome");
+    assert_eq!(best, &[0, 1], "the winner should be favoured");
+
+    let wins = p.win_probabilities();
+    assert!(wins[0] > wins[1], "{wins:?}");
+}
+
+/// Unequal team sizes change the draw margin, because inference derives it
+/// from the teams' betas. Prediction has to follow, or it describes a
+/// different model than the one that will be fitted.
+#[test]
+fn team_size_affects_the_prediction() {
+    let mut h = History::builder().p_draw(0.2).build();
+    h.event(1)
+        .team(["a", "b"])
+        .team(["c"])
+        .winner(0)
+        .commit()
+        .unwrap();
+    let _ = h.converge().unwrap();
+
+    let p = h.predict_outcome(&[&[&"a", &"b"], &[&"c"]]).unwrap();
+    assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total());
+    assert!(p.probability_of(&[0, 0]) > 0.0);
+}
+
+// ---------------------------------------------------------------------------
+// Expected information gain
+// ---------------------------------------------------------------------------
+
+/// The whole point of #39: "which comparison should I run next?" is a
+/// different question from "who will win?" or "is this fair?".
+#[test]
+fn information_gain_prefers_the_uncertain_pairing() {
+    let mut h = History::builder().build();
+
+    // "known" and "rival" have played a lot; "newcomer" has played once.
+    for t in 1..=15 {
+        h.record_winner(&"known", &"rival", t).unwrap();
+        h.record_winner(&"rival", &"known", t + 100).unwrap();
+    }
+    h.record_winner(&"known", &"newcomer", 500).unwrap();
+    let _ = h.converge().unwrap();
+
+    let settled = h
+        .expected_information_gain(&[&[&"known"], &[&"rival"]])
+        .unwrap();
+    let unknown = h
+        .expected_information_gain(&[&[&"known"], &[&"newcomer"]])
+        .unwrap();
+
+    assert!(
+        unknown > settled,
+        "pairing against the newcomer should teach more: {unknown} vs {settled}"
+    );
+}
+
+/// The analytic ceiling, through the `History` entry point rather than the
+/// standalone one.
+#[test]
+fn information_gain_respects_the_entropy_ceiling() {
+    let h = history_with(&["a", "b", "c"], 0.0);
+
+    let two = h.expected_information_gain(&[&[&"a"], &[&"b"]]).unwrap();
+    assert!(
+        (0.0..=std::f64::consts::LN_2).contains(&two),
+        "two-team EIG {two} outside [0, ln 2]"
+    );
+
+    let three = h
+        .expected_information_gain(&[&[&"a"], &[&"b"], &[&"c"]])
+        .unwrap();
+    assert!(
+        (0.0..=6.0f64.ln()).contains(&three),
+        "three-team EIG {three} outside [0, ln 6]"
+    );
+}
+
+#[test]
+fn information_gain_reports_unknown_keys() {
+    let h = history_with(&["a", "b"], 0.0);
+    assert_eq!(
+        h.expected_information_gain(&[&[&"a"], &[&"ghost"]])
+            .unwrap_err(),
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"ghost\"".to_owned(),
+        }
+    );
+}
+
+/// A draw-enabled history has three outcomes to weigh rather than two, so the
+/// draw branch must actually be reachable through this path.
+#[test]
+fn information_gain_accounts_for_draws() {
+    let with_draws = history_with(&["a", "b"], 0.25);
+    let g = with_draws
+        .expected_information_gain(&[&[&"a"], &[&"b"]])
+        .unwrap();
+    assert!(g > 0.0 && g <= 3.0f64.ln(), "{g}");
+
+    // The draw outcome carries mass, so it is genuinely being weighed.
+    let dist = with_draws.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap();
+    assert!(dist.probability_of(&[0, 0]) > 0.0);
+}
+
+/// The defect that cost a consumer a day: `UnknownKey { team: 0, member: 0 }`
+/// says nothing about *which* key is unknown, so the natural handling — log it,
+/// fall back to a neutral value — converts a total miss into a plausible
+/// constant. The key has to be in the error, and in its `Display`.
+#[test]
+fn unknown_key_names_the_key_it_could_not_find() {
+    let h = history_with(&["a", "b"], 0.0);
+    let err = h.predict_outcome(&[&[&"a"], &[&"never_seen"]]).unwrap_err();
+
+    match &err {
+        InferenceError::UnknownKey { key, .. } => {
+            assert!(
+                key.contains("never_seen"),
+                "the error should name the key, got {key}"
+            );
+        }
+        other => panic!("expected UnknownKey, got {other:?}"),
+    }
+
+    let rendered = err.to_string();
+    assert!(
+        rendered.contains("never_seen"),
+        "Display should name the key: {rendered}"
+    );
+    assert!(
+        rendered.contains("pre-filter"),
+        "Display should say what to do about it: {rendered}"
+    );
+}
+
+// ---------------------------------------------------------------------------
+// UnknownKeys policy
+// ---------------------------------------------------------------------------
+
+fn history_with_policy(names: &[&'static str], policy: trueskill_tt::UnknownKeys) -> History {
+    let mut h = History::builder().unknown_keys(policy).build();
+    for pair in names.windows(2) {
+        h.record_winner(&pair[0], &pair[1], 1).unwrap();
+    }
+    let _ = h.converge().unwrap();
+    h
+}
+
+#[test]
+fn reject_is_the_default() {
+    let h = history_with(&["a", "b"], 0.0);
+    assert!(matches!(
+        h.predict_outcome(&[&[&"a"], &[&"ghost"]]),
+        Err(InferenceError::UnknownKey { .. })
+    ));
+}
+
+#[test]
+fn prior_answers_instead_of_erroring() {
+    let h = history_with_policy(&["a", "b"], trueskill_tt::UnknownKeys::Prior);
+    let p = h
+        .predict_outcome(&[&[&"a"], &[&"ghost"]])
+        .expect("Prior should answer rather than reject");
+    assert!((p.total() - 1.0).abs() < 1e-6);
+}
+
+/// Two competitors the model has never seen are genuinely a coin flip. The
+/// point is that this is now *derived* rather than a constant a caller
+/// substitutes after swallowing an error.
+#[test]
+fn two_unknown_competitors_are_an_honest_coin_flip() {
+    let h = history_with_policy(&["a", "b"], trueskill_tt::UnknownKeys::Prior);
+    let wins = h
+        .predict_win_probabilities(&[&[&"nobody"], &[&"no_one"]])
+        .unwrap();
+    assert!((wins[0] - 0.5).abs() < 1e-9, "{wins:?}");
+    assert!((wins[1] - 0.5).abs() < 1e-9, "{wins:?}");
+}
+
+/// The property that rules out a `Skip` mode: an unknown member must make a
+/// team *less* certain, never more. Skipping would drop the member's variance
+/// from the sum and narrow the team, which is backwards.
+#[test]
+fn an_unknown_member_widens_its_team_rather_than_narrowing_it() {
+    let h = history_with_policy(&["a", "b", "c"], trueskill_tt::UnknownKeys::Prior);
+
+    // "a" alone against "b" — then "a" plus an unknown partner against "b".
+    let solo = h.predict_win_probabilities(&[&[&"a"], &[&"b"]]).unwrap();
+    let with_unknown = h
+        .predict_win_probabilities(&[&[&"a", &"stranger"], &[&"b"]])
+        .unwrap();
+
+    // Adding an unknown partner pulls the outcome toward even, because the
+    // team's performance spread grew.
+    assert!(
+        (with_unknown[0] - 0.5).abs() < (solo[0] - 0.5).abs(),
+        "an unknown partner should make the result less certain: solo {solo:?}, \
+         with unknown {with_unknown:?}"
+    );
+}
+
+#[test]
+fn prior_reaches_every_prediction_entry_point() {
+    let h = history_with_policy(&["a", "b"], trueskill_tt::UnknownKeys::Prior);
+    let teams: &[&[&&str]] = &[&[&"a"], &[&"ghost"]];
+
+    assert!(h.predict_quality(teams).is_ok());
+    assert!(h.predict_win_probabilities(teams).is_ok());
+    assert!(h.predict_outcome(teams).is_ok());
+    assert!(h.predict_ranking(teams, &[0, 1]).is_ok());
+    assert!(h.expected_information_gain(teams).is_ok());
+}

tests/properties.proptest-regressions

@@ -0,0 +1,7 @@
+# Seeds for failure cases proptest has generated in the past. It is
+# automatically read and these particular cases re-run before any
+# novel cases are generated.
+#
+# It is recommended to check this file in to source control so that
+# everyone who runs the test benefits from these saved cases.
+cc 8859be600e638573980f78622b8fcd8b4553ca34a9a041746c417f7e4293f89c # shrinks to games = [(0, 1), (0, 1), (4, 6), (2, 0), (0, 1), (0, 6), (0, 1), (2, 0), (6, 4), (0, 1), (0, 2), (6, 4), (1, 0), (4, 0), (0, 2)]

tests/properties.rs

@@ -26,7 +26,10 @@ const KEYS: [&str; 8] = ["a", "b", "c", "d", "e", "f", "g", "h"];
 fn history_from(games: &[(usize, usize)]) -> History {
     let mut h = History::builder()
         .convergence(ConvergenceOptions {
-            max_iter: 200,
+            // 200 was not enough: the batched side stopped at the cap with a
+            // step of 3.4e-9, so this test was comparing two truncated fits and
+            // attributing the gap to ingestion order.
+            max_iter: 20_000,
             epsilon: 1e-10,
             ..ConvergenceOptions::default()
         })
@@ -61,7 +64,7 @@ proptest! {
     fn converged_posteriors_are_always_finite(games in pairs()) {
         let mut h = history_from(&games);
 
-        h.converge().unwrap();
+        let _ = h.converge().unwrap();
 
         for key in KEYS {
             for (time, g) in h.learning_curve(key) {
@@ -79,7 +82,7 @@ proptest! {
     fn log_evidence_is_a_finite_log_probability(games in pairs()) {
         let mut h = history_from(&games);
 
-        h.converge().unwrap();
+        let _ = h.converge().unwrap();
 
         let batch = h.log_evidence();
         let filtered = h.filtered_log_evidence();
@@ -98,7 +101,7 @@ proptest! {
 
         let before = h.filtered_log_evidence();
 
-        h.converge().unwrap();
+        let _ = h.converge().unwrap();
 
         let after = h.filtered_log_evidence();
 
@@ -114,14 +117,21 @@ proptest! {
     fn ingestion_order_does_not_change_the_answer(games in pairs()) {
         let batched = {
             let mut h = history_from(&games);
-            h.converge().unwrap();
+            let report = h.converge().unwrap();
+            prop_assert!(
+                report.converged,
+                "batched side stopped at {} iterations with step {:?}; comparing \
+                 two fits that have not converged measures truncation, not order",
+                report.iterations,
+                report.final_step
+            );
             h
         };
 
         let incremental = {
             let mut h = History::builder()
                 .convergence(ConvergenceOptions {
-                    max_iter: 200,
+                    max_iter: 20_000,
                     epsilon: 1e-10,
                     ..ConvergenceOptions::default()
                 })
@@ -139,7 +149,13 @@ proptest! {
                 .unwrap();
             }
 
-            h.converge().unwrap();
+            let report = h.converge().unwrap();
+            prop_assert!(
+                report.converged,
+                "incremental side stopped at {} iterations with step {:?}",
+                report.iterations,
+                report.final_step
+            );
             h
         };
 

tests/quality.rs

@@ -108,12 +108,59 @@ fn history_predict_quality_supports_three_teams() {
     let mut h = History::default();
     h.record_winner(&"a", &"b", 1).unwrap();
     h.record_winner(&"b", &"c", 2).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
-    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]);
+    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]).unwrap();
     assert!(
         q.is_finite(),
         "3-team predict_quality must be finite, got {q}"
     );
     assert!((0.0..=1.0).contains(&q), "out of range: {q}");
 }
+
+/// `quality()` for N identical teams has a closed form, which pins the N-group
+/// determinant path across the whole range rather than at a single golden.
+///
+/// For two identical single-player teams the standard result is
+/// `sqrt(2b^2 / (2b^2 + s1^2 + s2^2))`. With the conventional parameters
+/// (`sigma = 25/3`, `beta = 25/6`) that ratio is exactly `1/5`, and the N-group
+/// generalisation is `(1/5)^((n-1)/2)` — one factor per adjacent pair.
+///
+/// The n=3 and n=5 values this produces (0.200 and 0.040) are also what the
+/// `trueskill` Python package returns for the same configuration, so this
+/// doubles as the cross-implementation check the README asked for.
+#[test]
+fn quality_of_identical_teams_follows_its_closed_form() {
+    let g = Gaussian::from_ms(25.0, 25.0 / 3.0);
+    let beta = 25.0 / 6.0;
+
+    for n in 2..=10usize {
+        let groups: Vec<Vec<Gaussian>> = (0..n).map(|_| vec![g]).collect();
+        let refs: Vec<&[Gaussian]> = groups.iter().map(Vec::as_slice).collect();
+
+        let got = quality(&refs, beta);
+        let expected = 0.2f64.powf((n - 1) as f64 / 2.0);
+
+        assert!(
+            (got - expected).abs() / expected < 1e-9,
+            "n={n}: quality {got}, closed form {expected}"
+        );
+    }
+}
+
+/// Spot-check against the two values the `trueskill` Python package is known
+/// to produce for this configuration, stated as literals so a future change to
+/// the closed-form reasoning above cannot quietly take these with it.
+#[test]
+fn quality_matches_the_reference_implementation() {
+    let g = Gaussian::from_ms(25.0, 25.0 / 3.0);
+    let beta = 25.0 / 6.0;
+
+    let three: Vec<Vec<Gaussian>> = (0..3).map(|_| vec![g]).collect();
+    let refs: Vec<&[Gaussian]> = three.iter().map(Vec::as_slice).collect();
+    assert!((quality(&refs, beta) - 0.200).abs() < 1e-9);
+
+    let five: Vec<Vec<Gaussian>> = (0..5).map(|_| vec![g]).collect();
+    let refs: Vec<&[Gaussian]> = five.iter().map(Vec::as_slice).collect();
+    assert!((quality(&refs, beta) - 0.040).abs() < 1e-9);
+}

tests/reconvergence_equivalence.rs

@@ -0,0 +1,165 @@
+//! Converging, appending, and converging again must reach the same fixed point
+//! as converging once over the whole event set.
+//!
+//! `tests/ingestion_equivalence.rs` covers a different question: it varies how
+//! events are *batched* but converges only at the end. This file converges
+//! between batches, which is the path a caller takes when it fits, serves for a
+//! while, then ingests more.
+//!
+//! The property matters beyond ergonomics. It says `converge` reaches a fixed
+//! point determined by the events, ratings and configuration alone — not by the
+//! message state it started from. That is what makes a restored snapshot safe:
+//! an inexact one cannot corrupt the answer, only cost an extra sweep. See #45.
+
+use smallvec::smallvec;
+use trueskill_tt::{ConvergenceOptions, Event, Gaussian, History, Member, Outcome, Team};
+
+fn tight() -> ConvergenceOptions {
+    ConvergenceOptions {
+        max_iter: 5_000,
+        epsilon: 1e-12,
+        alpha: 1.0,
+    }
+}
+
+fn ev(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),
+    }
+}
+
+/// Ingest each chunk in turn, converging fully after every one.
+fn fit_in_chunks(chunks: Vec<Events>) -> Vec<(String, Gaussian)> {
+    let mut h: History<i64, _, _, String> =
+        History::builder_with_key().convergence(tight()).build();
+
+    for chunk in chunks {
+        h.add_events(chunk).unwrap();
+        let report = h.converge().unwrap();
+        assert!(
+            report.converged,
+            "a chunk failed to converge, so any comparison would be measuring \
+             truncation rather than the fixed point; final step {:?}",
+            report.final_step
+        );
+    }
+
+    let mut skills: Vec<(String, Gaussian)> = h
+        .learning_curves()
+        .into_iter()
+        .map(|(k, curve)| (k, curve.last().unwrap().1))
+        .collect();
+    skills.sort_by(|a, b| a.0.cmp(&b.0));
+    skills
+}
+
+fn assert_same(a: &[(String, Gaussian)], b: &[(String, Gaussian)], what: &str) {
+    assert_eq!(a.len(), b.len(), "{what}: competitor count differs");
+    for ((ka, ga), (kb, gb)) in a.iter().zip(b) {
+        assert_eq!(ka, kb, "{what}: key order differs");
+        // Measured: 6.2e-13 for a later append, 8.9e-11 for an interleaved one.
+        // The bar is well clear of both but far under anything that would let a
+        // genuine divergence through.
+        assert!(
+            (ga.mu() - gb.mu()).abs() < 1e-8 && (ga.sigma() - gb.sigma()).abs() < 1e-8,
+            "{what}: {ka} differs — one-shot mu={} sigma={}, chunked mu={} sigma={}",
+            ga.mu(),
+            ga.sigma(),
+            gb.mu(),
+            gb.sigma()
+        );
+    }
+}
+
+type Events = Vec<Event<i64, String>>;
+
+/// Two chunks of events: the first at times 0..20, the second at 100..120.
+fn fixture() -> (Events, Events) {
+    let names = ["a", "b", "c", "d", "e"];
+    let mut seed = 7u64;
+    let mut rnd = move || {
+        seed ^= seed << 13;
+        seed ^= seed >> 7;
+        seed ^= seed << 17;
+        seed
+    };
+
+    let (mut early, mut late) = (Vec::new(), Vec::new());
+    for t in 0..40i64 {
+        let i = (rnd() % 5) as usize;
+        let mut j = (rnd() % 5) as usize;
+        if j == i {
+            j = (j + 1) % 5;
+        }
+        if t < 20 {
+            early.push(ev(names[i], names[j], t));
+        } else {
+            late.push(ev(names[i], names[j], 100 + t));
+        }
+    }
+    (early, late)
+}
+
+/// The ordinary case: new events are strictly later than everything fitted.
+#[test]
+fn appending_later_events_matches_a_single_fit() {
+    let (early, late) = fixture();
+    let all: Vec<_> = early.iter().cloned().chain(late.iter().cloned()).collect();
+
+    assert_same(
+        &fit_in_chunks(vec![all]),
+        &fit_in_chunks(vec![early, late]),
+        "append strictly later",
+    );
+}
+
+/// The case the design question suspected might be weaker: appended events
+/// interleave with slices that are already fitted, so the append legitimately
+/// revises the past. It is not weaker — Through Time revises the past on every
+/// converge regardless, so there is nothing special about doing it in two steps.
+#[test]
+fn appending_interleaved_events_matches_a_single_fit() {
+    let (early, late) = fixture();
+    let all: Vec<_> = early.iter().cloned().chain(late.iter().cloned()).collect();
+
+    // Split by parity so the second chunk is back-dated into the first's range.
+    let first: Vec<_> = all.iter().step_by(2).cloned().collect();
+    let second: Vec<_> = all.iter().skip(1).step_by(2).cloned().collect();
+    let together: Vec<_> = first
+        .iter()
+        .cloned()
+        .chain(second.iter().cloned())
+        .collect();
+
+    assert_same(
+        &fit_in_chunks(vec![together]),
+        &fit_in_chunks(vec![first, second]),
+        "append interleaved",
+    );
+}
+
+/// Converging an already-converged history is a no-op, which is what makes a
+/// restored snapshot worth having: the work is skipped rather than redone.
+#[test]
+fn re_converging_an_unchanged_history_costs_one_iteration() {
+    let (early, late) = fixture();
+    let all: Vec<_> = early.into_iter().chain(late).collect();
+
+    let mut h: History<i64, _, _, String> =
+        History::builder_with_key().convergence(tight()).build();
+    h.add_events(all).unwrap();
+    let first = h.converge().unwrap();
+    assert!(first.converged);
+
+    let again = h.converge().unwrap();
+    assert_eq!(
+        again.iterations, 1,
+        "a converged history should settle immediately, not re-grind"
+    );
+    assert!(again.converged);
+}

tests/record_winner.rs

@@ -15,7 +15,7 @@ fn record_winner_builds_history() {
         .build();
 
     h.record_winner(&"alice", &"bob", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let a_idx = h.lookup(&"alice").unwrap();
     let b_idx = h.lookup(&"bob").unwrap();
@@ -48,7 +48,7 @@ fn record_draw_with_p_draw_set() {
         .build();
 
     h.record_draw(&"alice", &"bob", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     assert!(h.lookup(&"alice").is_some());
     assert!(h.lookup(&"bob").is_some());

tests/validation.rs

@@ -0,0 +1,186 @@
+//! Input validation must hold in **release**, where `debug_assert!` is gone.
+//!
+//! The engine guards itself with `debug_assert!`, which documents invariants
+//! but vanishes in the profile users actually ship. Anything reachable from the
+//! public API has to be rejected with an `InferenceError` instead, at the
+//! boundary, rather than becoming NaN or an out-of-bounds panic deep inside
+//! `run_chain`.
+//!
+//! `GameOptions` and `ConvergenceOptions` both have public fields, so the
+//! eager asserts on `HistoryBuilder` do not cover the `Game` constructors —
+//! a caller can build the options struct directly.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, Game, GameOptions, Gaussian, History, InferenceError,
+    Member, Outcome, Rating, Team,
+};
+
+type R = Rating<i64, ConstantDrift>;
+
+fn rating() -> R {
+    R::new(
+        Gaussian::from_ms(25.0, 25.0 / 3.0),
+        25.0 / 6.0,
+        ConstantDrift(0.0),
+    )
+}
+
+fn options_with_alpha(alpha: f64) -> GameOptions {
+    GameOptions {
+        convergence: ConvergenceOptions {
+            alpha,
+            ..ConvergenceOptions::default()
+        },
+        ..GameOptions::default()
+    }
+}
+
+/// `alpha == 0.0` leaves every EP update unapplied, so inference silently
+/// returns the priors — the worst possible failure, since the output looks
+/// entirely reasonable.
+#[test]
+fn ranked_rejects_a_zero_damping_factor() {
+    let (a, b) = (rating(), rating());
+    let err = Game::<i64, _>::ranked(
+        &[&[a], &[b]],
+        Outcome::winner(0, 2),
+        &options_with_alpha(0.0),
+    )
+    .expect_err("alpha = 0 must be rejected");
+    assert!(
+        matches!(err, InferenceError::InvalidParameter { name: "alpha", .. }),
+        "got {err:?}"
+    );
+}
+
+#[test]
+fn ranked_rejects_an_out_of_range_damping_factor() {
+    let (a, b) = (rating(), rating());
+    for alpha in [-0.5, 1.5, f64::NAN] {
+        let err = Game::<i64, _>::ranked(
+            &[&[a], &[b]],
+            Outcome::winner(0, 2),
+            &options_with_alpha(alpha),
+        )
+        .expect_err("alpha out of (0, 1] must be rejected");
+        assert!(
+            matches!(err, InferenceError::InvalidParameter { name: "alpha", .. }),
+            "alpha={alpha}: got {err:?}"
+        );
+    }
+}
+
+#[test]
+fn scored_rejects_a_bad_damping_factor() {
+    let (a, b) = (rating(), rating());
+    let err = Game::<i64, _>::scored(
+        &[&[a], &[b]],
+        Outcome::scores([21.0, 9.0]),
+        &options_with_alpha(0.0),
+    )
+    .expect_err("alpha = 0 must be rejected");
+    assert!(
+        matches!(err, InferenceError::InvalidParameter { name: "alpha", .. }),
+        "got {err:?}"
+    );
+}
+
+/// Already covered by `Game::ranked`, asserted here so the release-mode
+/// guarantee is stated in one place.
+#[test]
+fn ranked_rejects_an_out_of_range_draw_probability() {
+    let (a, b) = (rating(), rating());
+    for p_draw in [-0.5, 1.0, 1.5] {
+        let options = GameOptions {
+            p_draw,
+            ..GameOptions::default()
+        };
+        assert!(
+            Game::<i64, _>::ranked(&[&[a], &[b]], Outcome::winner(0, 2), &options).is_err(),
+            "p_draw={p_draw} must be rejected"
+        );
+    }
+}
+
+#[test]
+fn scored_rejects_a_non_positive_noise() {
+    let (a, b) = (rating(), rating());
+    for score_sigma in [0.0, -1.0, f64::NAN] {
+        let options = GameOptions {
+            score_sigma,
+            ..GameOptions::default()
+        };
+        assert!(
+            Game::<i64, _>::scored(&[&[a], &[b]], Outcome::scores([21.0, 9.0]), &options).is_err(),
+            "score_sigma={score_sigma} must be rejected"
+        );
+    }
+}
+
+/// A tie with no draw probability makes the truncation margin zero and the
+/// two-sided update evaluate 0/0. Ingestion must refuse it.
+#[test]
+fn ingestion_rejects_a_tie_without_a_draw_probability() {
+    let mut h = History::builder().p_draw(0.0).build();
+    let err = h
+        .add_events(vec![Event {
+            time: 0,
+            teams: smallvec![
+                Team::with_members([Member::new("a")]),
+                Team::with_members([Member::new("b")]),
+            ],
+            outcome: Outcome::draw(2),
+        }])
+        .expect_err("a tie with p_draw = 0 must be rejected");
+    assert!(
+        matches!(err, InferenceError::TieWithoutDrawProbability { .. }),
+        "got {err:?}"
+    );
+}
+
+/// `Outcome::scores_with_sigma` documents that a non-positive sigma is
+/// accepted at construction and rejected at ingestion.
+#[test]
+fn ingestion_rejects_a_non_positive_per_event_score_sigma() {
+    for sigma in [0.0, -1.0, f64::NAN] {
+        let mut h = History::builder().build();
+        let err = h
+            .add_events(vec![Event {
+                time: 0,
+                teams: smallvec![
+                    Team::with_members([Member::new("a")]),
+                    Team::with_members([Member::new("b")]),
+                ],
+                outcome: Outcome::scores_with_sigma([21.0, 9.0], sigma),
+            }])
+            .expect_err("a non-positive per-event sigma must be rejected");
+        assert!(
+            matches!(err, InferenceError::InvalidParameter { .. }),
+            "sigma={sigma}: got {err:?}"
+        );
+    }
+}
+
+/// Per-team weights must match that team's membership. The top-level length
+/// checks in ingestion do not cover the inner dimension.
+#[test]
+fn ingestion_rejects_weights_that_do_not_match_their_team() {
+    let mut h = History::builder().build();
+    let mut team = Team::with_members([Member::new("a"), Member::new("b")]);
+    team.members[0].weight = 1.0;
+
+    let err = h
+        .event(0)
+        .team(["a", "b"])
+        .team(["c"])
+        // Three weights for a two-member team.
+        .weights([1.0, 1.0, 1.0])
+        .winner(0)
+        .commit()
+        .expect_err("a weight/member length mismatch must be rejected");
+    assert!(
+        matches!(err, InferenceError::MismatchedShape { .. }),
+        "got {err:?}"
+    );
+}

tests/variance_reduction.rs

@@ -0,0 +1,156 @@
+//! `expected_variance_reduction`: which matchup best sharpens a given question.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, InferenceError, Member, Outcome, Team,
+    UnknownKeys,
+};
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn round(a: &'static str, b: &'static str, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: 1,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+fn base() -> Vec<Event<i64, &'static str>> {
+    vec![
+        round("a", "b", 5.0, 2.0),
+        round("a", "c", 6.0, 1.0),
+        round("b", "c", 4.0, 3.0),
+        round("c", "d", 2.0, 1.0),
+        round("a", "d", 7.0, 2.0),
+    ]
+}
+
+fn fit(extra: Option<Event<i64, &'static str>>, policy: UnknownKeys) -> H {
+    let mut h: History<i64, _, _, &'static str> = History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.0))
+        .unknown_keys(policy)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build();
+    let mut ev = base();
+    if let Some(e) = extra {
+        ev.push(e);
+    }
+    h.add_events(ev).unwrap();
+    let _ = h.converge().unwrap();
+    h
+}
+
+/// The closed form must equal what actually happens if the matchup is played.
+/// This is the assertion that makes the whole call trustworthy: a wrong
+/// acquisition function returns plausible numbers and quietly picks worse
+/// matchups forever.
+#[test]
+fn the_closed_form_matches_an_actual_refit() {
+    let h = fit(None, UnknownKeys::Reject);
+    let target: Vec<(&&str, f64)> = vec![(&"a", 1.0), (&"b", -1.0)];
+    let before = h.posterior_of(&target).unwrap().sigma().powi(2);
+
+    for (x, y) in [("a", "b"), ("c", "d"), ("a", "c"), ("b", "d")] {
+        let predicted = h
+            .expected_variance_reduction(&[&[&x], &[&y]], &target)
+            .unwrap();
+
+        let after = fit(Some(round(x, y, 3.0, 1.0)), UnknownKeys::Reject);
+        let actual = before - after.posterior_of(&target).unwrap().sigma().powi(2);
+
+        assert!(
+            (predicted - actual).abs() / actual.abs() < 1e-9,
+            "{x} vs {y}: predicted {predicted}, actual {actual}"
+        );
+    }
+}
+
+/// The reduction cannot depend on the score, because for a Gaussian likelihood
+/// the posterior variance update is data-independent. This is why the call
+/// needs no expectation despite its name.
+#[test]
+fn the_outcome_does_not_change_the_reduction() {
+    let target: Vec<(&&str, f64)> = vec![(&"a", 1.0), (&"b", -1.0)];
+    let h = fit(None, UnknownKeys::Reject);
+    let before = h.posterior_of(&target).unwrap().sigma().powi(2);
+
+    let mut seen = Vec::new();
+    for (sa, sb) in [(3.0, 1.0), (100.0, -50.0), (0.0, 0.0)] {
+        let after = fit(Some(round("c", "d", sa, sb)), UnknownKeys::Reject);
+        seen.push(before - after.posterior_of(&target).unwrap().sigma().powi(2));
+    }
+    for w in seen.windows(2) {
+        assert!(
+            (w[0] - w[1]).abs() < 1e-12,
+            "variance reduction moved with the observed score: {seen:?}"
+        );
+    }
+}
+
+/// The point of the call: it must rank candidate matchups usefully. Playing the
+/// pair you are trying to separate helps most; an unrelated pair helps least.
+#[test]
+fn it_ranks_candidates_by_how_much_they_answer_the_question() {
+    let h = fit(None, UnknownKeys::Reject);
+    let target: Vec<(&&str, f64)> = vec![(&"a", 1.0), (&"b", -1.0)];
+
+    let direct = h
+        .expected_variance_reduction(&[&[&"a"], &[&"b"]], &target)
+        .unwrap();
+    let unrelated = h
+        .expected_variance_reduction(&[&[&"c"], &[&"d"]], &target)
+        .unwrap();
+
+    assert!(direct > 0.0 && unrelated > 0.0);
+    assert!(
+        direct > 5.0 * unrelated,
+        "playing the target pair should dominate: {direct} vs {unrelated}"
+    );
+}
+
+/// A matchup between two competitors nobody has seen still teaches something
+/// about them, but nothing about a target that does not involve them.
+#[test]
+fn an_unrelated_unseen_matchup_teaches_nothing_about_the_target() {
+    let h = fit(None, UnknownKeys::Prior);
+    let target: Vec<(&&str, f64)> = vec![(&"a", 1.0), (&"b", -1.0)];
+
+    let reduction = h
+        .expected_variance_reduction(&[&[&"stranger"], &[&"nobody"]], &target)
+        .unwrap();
+    assert!(
+        reduction.abs() < 1e-12,
+        "an unseen pair shares nothing with the target: {reduction}"
+    );
+}
+
+#[test]
+fn shape_errors_are_reported() {
+    let h = fit(None, UnknownKeys::Reject);
+    let target: Vec<(&&str, f64)> = vec![(&"a", 1.0), (&"b", -1.0)];
+
+    assert!(matches!(
+        h.expected_variance_reduction(&[&[&"a"]], &target),
+        Err(InferenceError::MismatchedShape {
+            expected: 2,
+            got: 1,
+            ..
+        })
+    ));
+    assert!(matches!(
+        h.expected_variance_reduction(&[&[&"a"], &[&"ghost"]], &target),
+        Err(InferenceError::UnknownKey { .. })
+    ));
+}