Merge branch 'feat/convergence-strictness'

Make a short fit an error, raise the default iteration cap, validate the
remaining HistoryBuilder parameters, add History::register and
History::rating, reject competitor config conflicts across batches, and
document what the joint's cost scales in.

Closes #50

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

CHANGELOG.md

@@ -2,12 +2,88 @@
 
 All notable changes to this project will be documented in this file.
 
+## 0.7.0 - 2026-09-08
+
+### Features
+
+- feat: factorise the joint once with History::joint
+
+### Other (unconventional)
+
+- Merge branch 'feat/joint-handle'
+
+## 0.6.0 - 2026-09-08
+
+### Breaking Changes
+
+- fix!: make the joint span slices, not just the latest one
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.6.0
+
+## 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
+
+### Miscellaneous Tasks
+
+- chore: Release trueskill-tt version 0.5.0
+
+### 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

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):
@@ -66,11 +80,11 @@ History  →  TimeSlice[]  →  Event[]  →  Item[]
   `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, `pub(crate)`) and `CompetitorStore`
@@ -97,6 +111,12 @@ History  →  TimeSlice[]  →  Event[]  →  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.4.1"
+version = "0.7.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"
 
@@ -78,3 +79,7 @@ debug = true
 
 [profile.dev]
 debug = true
+
+[[bench]]
+name = "joint"
+harness = false

README.md

@@ -134,14 +134,20 @@ h.add_events(vec![Event {
 h.converge().unwrap();
 ```
 
-Like `with_prior`, the scale is **competitor configuration captured at first
-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`.
+Like `with_prior`, the scale is **competitor configuration, not a per-event
+value**: it applies to the competitor for the whole history, and it applies
+whenever it is supplied — including on a key the history already knows.
+Configuring one late still refits the whole history rather than taking effect
+only from that event onward, because `converge` refits from competitor state.
+Repeating the same value is inert; supplying two *different* values for one
+competitor within a single batch is `InferenceError::ConflictingCompetitorConfig`,
+since events in a batch have no order. The scale 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.
+The fluent `EventBuilder` reaches this too: `.team([...])` is the common case
+and leaves both unset, while `.members([...])` takes `Member` values directly,
+so `h.event(t).members([Member::new("layout_7").with_drift_scale(0.0)])` is
+equivalent to the typed shape above.
 
 ## Scored outcomes
 
@@ -195,8 +201,47 @@ stay available at any size:
   quadratic in team count.
 - `predict_ranking(teams, ranks)` — one specific finishing order.
 
-Unknown keys are an error, not a silent omission: a team the history has never
-seen cannot produce a confident-looking probability.
+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
 

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/joint.rs

@@ -0,0 +1,71 @@
+//! Cost of the joint posterior: factorising versus querying.
+//!
+//! The split is the whole point of `History::joint`. Factorising is `O(n^3)` in
+//! the history's appearances and depends only on the fit; a query is `O(n^2)`
+//! and depends only on the question. `posterior_of_one_shot` pays both every
+//! time, `joint_query` pays only the second.
+
+use criterion::{Criterion, criterion_group, criterion_main};
+use smallvec::smallvec;
+use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team};
+
+/// 30 slices of 8 duels: 480 appearances over 100 competitors.
+fn fitted() -> History<i64, ConstantDrift, trueskill_tt::NullObserver, String> {
+    let mut h: History<i64, ConstantDrift, _, String> = History::builder_with_key()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.05))
+        .convergence(ConvergenceOptions {
+            max_iter: 30,
+            epsilon: 1e-10,
+            alpha: 1.0,
+        })
+        .build();
+
+    let mut events: Vec<Event<i64, String>> = Vec::new();
+    let mut k = 0usize;
+    for t in 0..30i64 {
+        for _ in 0..8 {
+            k += 1;
+            events.push(Event {
+                time: t,
+                teams: smallvec![
+                    Team::with_members([Member::new(format!("p{}", k % 100))]),
+                    Team::with_members([Member::new(format!("p{}", (k + 37) % 100))]),
+                ],
+                outcome: Outcome::scores([
+                    (k as f64 * 0.3).sin().abs() * 20.0,
+                    (k as f64 * 0.3).cos().abs() * 20.0,
+                ]),
+            });
+        }
+    }
+    h.add_events(events).unwrap();
+    let _ = h.converge().unwrap();
+    h
+}
+
+fn bench_joint(c: &mut Criterion) {
+    let h = fitted();
+    let a = "p0".to_string();
+    let b = "p1".to_string();
+    let terms = [(&a, 1.0), (&b, -1.0)];
+
+    c.bench_function("joint_factorise_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(h.joint().unwrap().variables()));
+    });
+
+    c.bench_function("posterior_of_one_shot_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(h.posterior_of(&terms).unwrap()));
+    });
+
+    let joint = h.joint().unwrap();
+    c.bench_function("joint_query_480_appearances", |bencher| {
+        bencher.iter(|| std::hint::black_box(joint.posterior_of(&terms).unwrap()));
+    });
+}
+
+criterion_group!(benches, bench_joint);
+criterion_main!(benches);

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();
         });
     });
 }

docs/superpowers/specs/2026-04-23-trueskill-engine-redesign-design.md

@@ -500,6 +500,26 @@ All public traits (`Time`, `Drift`, `Observer`, `Factor`, `Schedule`) require `S
 
 `rayon` as default-on feature; with `default-features = false`, parallel paths fall back to sequential iterators behind `cfg(feature = "rayon")`.
 
+> **Not implemented. Deliberate deviation, decided 2026-09-08 (issue #5).**
+>
+> `rayon` ships **opt-in**: `Cargo.toml` has no `default = [...]` key. The
+> measured speedups are 1.0x on realistic workloads and 1.3x on a pathological
+> one (issue #4), because typical slices hold too few events to amortize
+> rayon's task-spawn overhead. Default-on would hand every downstream user a
+> thread pool and a dependency for approximately no gain.
+>
+> This section made the trade conditional on cross-slice dirty-bit skipping
+> landing and changing the parallel story. It did not land: #4 was closed on
+> 2026-08-27 by removing the inert `ConvergenceReport::slices_skipped` field
+> rather than by implementing the mechanism, so the re-measurement this was
+> waiting on will not arrive.
+>
+> The "Trade-offs" note below also cited an `unsafe` concurrent-write path
+> through `SkillStore` as a cost of default-on. That cost does not exist: the
+> crate is `#![forbid(unsafe_code)]`, and the compute/apply split on the
+> internal `Event` is what lets a color group run in parallel without it. The
+> case for opt-in rests on the measurements alone.
+
 ### Expected speedup ballpark
 
 For 1000 players, 60 events/slice × 1000 slices, 30 convergence iterations:
@@ -521,7 +541,7 @@ These are pre-implementation estimates. Each tier validates with criterion.
 - Color-group parallelism requires up-front graph coloring at ingestion. Cost: linear in events, run once per `add_events`. Cheap.
 - Default = asynchronous EP (preserves current semantics). Synchronous opt-in only.
 - Cross-slice sweep stays sequential; no speculative parallel sweeps.
-- Rayon default-on but feature-gated.
+- Rayon default-on but feature-gated. **Superseded — shipped opt-in; see the deviation note in Section 6.**
 
 ### Open question
 

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),

src/acquisition.rs

@@ -47,7 +47,7 @@ fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 {
     }
 
     let mean_gap = q.mu() - p.mu();
-    0.5 * ((var_p / var_q).ln() + (var_q + mean_gap * mean_gap) / var_p - 1.0)
+    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.

src/convergence.rs

@@ -62,7 +62,17 @@ impl Default for ConvergenceOptions {
 }
 
 /// Post-hoc summary of a `History::converge` call.
+///
+/// From [`History::converge`](crate::History::converge) this always describes a
+/// converged fit — stopping at `max_iter` is
+/// [`InferenceError::NotConverged`](crate::InferenceError::NotConverged) there.
+/// From [`History::converge_partial`](crate::History::converge_partial) it may
+/// not be, and `converged` is what says so.
 #[derive(Clone, Debug)]
+#[must_use = "from `converge_partial` this may describe a fit that stopped at \
+              `max_iter`, which is wrong by a little rather than loudly \
+              broken — check `converged`, 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,6 +64,24 @@ 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) },
+    /// The convergence sweep hit `max_iter` with the step still above
+    /// `epsilon`.
+    ///
+    /// A fit that stops short is wrong by a little, which is the worst
+    /// available failure: every rating is finite, the ordering looks sensible,
+    /// and nothing in the numbers says they were still moving. Reported rather
+    /// than returned as a flag on an `Ok`, because a flag has to be checked
+    /// and `let _ = h.converge()` is the natural way not to.
+    ///
+    /// Either the history needs more iterations — raise `max_iter` — or it is
+    /// oscillating rather than converging, in which case `alpha < 1.0` damps
+    /// the within-game EP loop. [`History::converge_partial`](crate::History::converge_partial)
+    /// returns the short fit instead when that is genuinely what is wanted.
+    NotConverged {
+        iterations: usize,
+        final_step: (f64, f64),
+        epsilon: f64,
+    },
     /// Inference produced a non-finite value (NaN or infinity).
     ///
     /// Indicates numerical breakdown; the resulting skills are meaningless
@@ -50,9 +107,32 @@ pub enum InferenceError {
     ///
     /// Reported rather than skipped: dropping unknown keys turns a team of
     /// strangers into a confident-looking probability about nobody.
-    UnknownKey { team: usize, member: usize },
+    ///
+    /// `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,
+    },
+    /// `History::register` was called for a competitor that already exists.
+    ///
+    /// Registration states a competitor's configuration before anything has
+    /// been observed about them, so a competitor that already exists has
+    /// already been configured — by an earlier `register`, or by an event that
+    /// created them. Silently overwriting would reintroduce exactly the
+    /// order-dependence registration exists to remove.
+    ///
+    /// To change an existing competitor's configuration, supply it on an event
+    /// through `Member`; that refits the whole history.
+    AlreadyRegistered { 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.
@@ -93,6 +173,18 @@ impl fmt::Display for InferenceError {
                     teams.0, teams.1
                 )
             }
+            Self::NotConverged {
+                iterations,
+                final_step,
+                epsilon,
+            } => {
+                write!(
+                    f,
+                    "did not converge in {iterations} iterations: final step {final_step:?} \
+                     is still above epsilon {epsilon}; raise max_iter, or damp with \
+                     alpha < 1.0 if it is oscillating"
+                )
+            }
             Self::NonFiniteResult { context, step } => {
                 write!(
                     f,
@@ -108,15 +200,28 @@ impl fmt::Display for InferenceError {
                     "competitor {competitor}: this batch sets {field} to two different values"
                 )
             }
-            Self::UnknownKey { team, member } => {
+            Self::UnknownKey { team, member, key } => {
                 write!(
                     f,
-                    "team {team}, member {member}: no skill recorded for this key"
+                    "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::AlreadyRegistered { key } => {
+                write!(
+                    f,
+                    "competitor {key} is already registered; registration states \
+                     configuration before anything is observed, so re-registering \
+                     would silently overwrite it"
                 )
             }
             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}")
             }

src/event.rs

@@ -88,7 +88,9 @@ impl<K> Member<K> {
 
     /// Set this competitor's starting skill estimate.
     ///
-    /// Captured at the competitor's first appearance; see the type docs.
+    /// Competitor configuration, not a per-event value: it applies for the
+    /// whole history and applies whenever it is supplied, including on a key
+    /// the history already knows. See the type docs.
     pub fn with_prior(mut self, prior: Gaussian) -> Self {
         self.prior = Some(prior);
         self
@@ -104,7 +106,8 @@ impl<K> Member<K> {
     /// shares a scale with moving competitors but should not itself move: a bot
     /// at a known strength, a rating floor, a course difficulty.
     ///
-    /// Captured at the competitor's first appearance; see the type docs.
+    /// Applies for the whole history and whenever it is supplied, including on
+    /// a key the history already knows; see the type docs.
     /// Must be finite and non-negative, or ingestion fails with
     /// [`InferenceError::InvalidParameter`](crate::InferenceError::InvalidParameter).
     pub fn with_drift_scale(mut self, scale: f64) -> Self {

src/event_builder.rs

@@ -50,6 +50,8 @@ where
     }
 
     /// Add a team by its member keys (weight 1.0 each, no prior overrides).
+    ///
+    /// Use [`EventBuilder::members`] to set `prior` or `drift_scale`.
     pub fn team<I: IntoIterator<Item = K>>(mut self, keys: I) -> Self {
         let members: SmallVec<[Member<K>; 4]> = keys.into_iter().map(Member::new).collect();
         self.event.teams.push(Team { members });
@@ -57,6 +59,40 @@ where
         self
     }
 
+    /// Add a team from fully-specified [`Member`] values.
+    ///
+    /// [`EventBuilder::team`] is the common case and builds members with
+    /// `Member::new`, which leaves `prior` and `drift_scale` unset. This is the
+    /// escape hatch for when they matter:
+    ///
+    /// ```
+    /// # use trueskill_tt::{Gaussian, History, Member};
+    /// # let mut h = History::builder().build();
+    /// h.event(0)
+    ///     .team(["player"])
+    ///     .members([Member::new("layout_7")
+    ///         .with_drift_scale(0.0)
+    ///         .with_prior(Gaussian::from_ms(0.0, 1.0))])
+    ///     .ranking([0, 1])
+    ///     .commit()?;
+    /// # Ok::<(), trueskill_tt::InferenceError>(())
+    /// ```
+    ///
+    /// One method rather than a `priors` and a `drift_scales` setter beside
+    /// `weights`: those would have to grow a parallel array — and a parallel
+    /// length check — every time `Member` gains a field, and each one would be
+    /// a new way to get the lengths wrong. `Member`'s own builder already
+    /// expresses all of it.
+    ///
+    /// `prior` and `drift_scale` are competitor configuration rather than
+    /// per-event values; see [`Member`] for what that means for a key the
+    /// history already knows.
+    pub fn members<I: IntoIterator<Item = Member<K>>>(mut self, members: I) -> Self {
+        self.event.teams.push(Team::with_members(members));
+        self.current_team_idx = Some(self.event.teams.len() - 1);
+        self
+    }
+
     /// Set per-member weights for the most recently added team.
     ///
     /// A length mismatch is recorded and returned by [`EventBuilder::commit`]

src/factor/margin.rs

@@ -1,6 +1,6 @@
 use crate::{
     N_INF,
-    factor::{Factor, VarId, VarStore},
+    factor::{VarId, VarStore},
     gaussian::Gaussian,
     ln_pdf,
 };
@@ -55,12 +55,16 @@ 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 {
+    pub(crate) fn log_evidence(&self) -> f64 {
         self.log_evidence_cached.unwrap_or(0.0)
     }
 }
@@ -74,7 +78,10 @@ impl Factor for MarginFactor {
 /// 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 {
-    let combined_sigma = (cavity.sigma().powi(2) + sigma.powi(2)).sqrt();
+    // `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);
 
     // A degenerate cavity (infinite sigma) is the only way to reach a
@@ -143,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,6 +1,6 @@
 use crate::{
     N_INF, approx,
-    factor::{Factor, VarId, VarStore},
+    factor::{VarId, VarStore},
     gaussian::Gaussian,
     ln_interval, ln_sf,
 };
@@ -63,14 +63,14 @@ 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.log_evidence_cached.unwrap_or(0.0)
-    }
 }
 
 /// `ln P(diff > margin)` for a win, `ln P(|diff| < margin)` for a tie.

src/game.rs

@@ -568,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::graph::Schedule>(
-        factors: &mut [crate::graph::BuiltinFactor],
-        vars: &mut crate::graph::VarStore,
-        schedule: &S,
-    ) -> crate::graph::ScheduleReport {
-        schedule.run(factors, vars)
-    }
 }
 
 #[cfg(test)]
@@ -1254,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!(
@@ -1314,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/graph.rs

@@ -1,20 +0,0 @@
-//! Factor-graph public API.
-//!
-//! Named `graph` rather than `factors` because the private implementation
-//! module beside it is `factor`: two module paths differing by one character,
-//! one public and one not, was a standing invitation to import the wrong one.
-//!
-//! The factor types, `VarStore` and the `Schedule` trait are public so custom
-//! schedules can be written against them.
-//!
-//! Building a factor graph by hand goes through `Game::custom`, which is
-//! deliberately `#[doc(hidden)]`: it works, but its signature is not yet
-//! considered stable API and so is not listed in these docs.
-
-pub use crate::{
-    factor::{
-        BuiltinFactor, Factor, VarId, VarStore, margin::MarginFactor, rank_diff::RankDiffFactor,
-        team_sum::TeamSumFactor, trunc::TruncFactor,
-    },
-    schedule::{EpsilonOrMax, Schedule, ScheduleReport},
-};

src/joint.rs

@@ -0,0 +1,152 @@
+//! Cholesky factorisation of a joint precision matrix.
+//!
+//! Every question the joint answers is a *bilinear form* in the precision
+//! matrix's inverse — the variance of a contrast is `c^T L^-1 c`, and the
+//! covariance of two contrasts is `c^T L^-1 a`. None of them wants `L^-1 c`
+//! itself, which is what makes the shape here worth stating explicitly.
+//!
+//! Writing the precision as `A = L L^T`,
+//!
+//! ```text
+//! c^T A^-1 a = c^T L^-T L^-1 a = (L^-1 c) . (L^-1 a)
+//! ```
+//!
+//! so a single forward substitution per contrast answers everything, and the
+//! back substitution a general solve would do is wasted work. That halves the
+//! cost of a query, and it removes a failure mode: a variance computed as
+//! `c . (A^-1 c)` is a difference of products that can round to a small
+//! negative number, where the same quantity as `|L^-1 c|^2` is a sum of
+//! squares and cannot.
+//!
+//! Factorising is `O(n^3)` and whitening is `O(n^2)`, so the split also
+//! matters structurally: the expensive half depends only on the fit, and is
+//! shared across every query a [`Joint`](crate::Joint) answers.
+
+/// A factorised symmetric positive-definite matrix, reusable across queries.
+pub(crate) struct Cholesky {
+    /// Lower triangle of `L`, row-major `n * n`. The upper triangle is
+    /// leftover scratch from the factorisation and is never read.
+    l: Vec<f64>,
+    n: usize,
+}
+
+impl Cholesky {
+    /// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`.
+    ///
+    /// `a` 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
+    /// neither a proper prior nor any evidence.
+    pub(crate) fn factor(mut a: Vec<f64>, n: usize) -> Option<Self> {
+        debug_assert_eq!(a.len(), n * n);
+
+        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;
+            }
+        }
+
+        Some(Self { l: a, n })
+    }
+
+    /// Whiten a contrast: `y = L^-1 b`.
+    ///
+    /// The point of the result is the dot product, not the vector: for two
+    /// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs.
+    pub(crate) fn whiten(&self, b: &[f64]) -> Vec<f64> {
+        debug_assert_eq!(b.len(), self.n);
+        let n = self.n;
+        let mut y = b.to_vec();
+        for i in 0..n {
+            // Folded from `y[i]` rather than summed and subtracted once, so the
+            // accumulation order matches a plain substitution loop exactly.
+            let row = &self.l[i * n..i * n + i];
+            let s = row
+                .iter()
+                .zip(&y[..i])
+                .fold(y[i], |acc, (l, v)| acc - l * v);
+            y[i] = s / self.l[i * n + i];
+        }
+        y
+    }
+}
+
+/// `b^T A^-1 b'`, given the two whitened contrasts.
+pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 {
+    y.iter().zip(y_prime).map(|(a, b)| a * b).sum()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    /// `[[4, 1], [1, 3]] z = [1, 2]` has `z = [1/11, 7/11]`, so the quadratic
+    /// form `b^T A^-1 b` is `1 * 1/11 + 2 * 7/11 = 15/11`.
+    #[test]
+    fn reproduces_a_known_quadratic_form() {
+        let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap();
+        let y = c.whiten(&[1.0, 2.0]);
+        assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12);
+    }
+
+    /// Whitening `e_i` recovers the inverse's diagonal, which is the variance
+    /// of a single variable.
+    #[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];
+        let c = Cholesky::factor(a, 3).unwrap();
+        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 y = c.whiten(&e);
+            assert!((bilinear(&y, &y) - expected).abs() < 1e-12, "row {i}");
+        }
+    }
+
+    /// The off-diagonal bilinear form is symmetric and matches the inverse.
+    #[test]
+    fn recovers_an_off_diagonal_covariance() {
+        // Same A; (A^-1)_{0,1} = 0.5.
+        let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0];
+        let c = Cholesky::factor(a, 3).unwrap();
+        let y0 = c.whiten(&[1.0, 0.0, 0.0]);
+        let y1 = c.whiten(&[0.0, 1.0, 0.0]);
+        assert!((bilinear(&y0, &y1) - 0.5).abs() < 1e-12);
+        assert!((bilinear(&y1, &y0) - 0.5).abs() < 1e-12);
+    }
+
+    /// A variance can never come out negative, because it is a sum of squares.
+    #[test]
+    fn a_quadratic_form_is_never_negative() {
+        let a = vec![1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12];
+        let c = Cholesky::factor(a, 2).unwrap();
+        let y = c.whiten(&[1.0, -1.0]);
+        assert!(bilinear(&y, &y) >= 0.0);
+    }
+
+    #[test]
+    fn rejects_a_non_positive_definite_matrix() {
+        // Singular: the second row is a multiple of the first.
+        assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none());
+    }
+}

src/lib.rs

@@ -121,8 +121,8 @@ mod event_builder;
 pub(crate) mod factor;
 mod game;
 pub mod gaussian;
-pub mod graph;
 mod history;
+mod joint;
 mod key_table;
 mod matrix;
 mod observer;
@@ -130,26 +130,24 @@ 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};
 pub use gaussian::Gaussian;
-pub use history::{History, HistoryBuilder};
+pub use history::{History, HistoryBuilder, Joint};
 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;
@@ -158,7 +156,45 @@ 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;
-pub const ITERATIONS: usize = 30;
+/// Default cap on convergence sweeps.
+///
+/// **A runaway guard, not a budget.** The sweep exits as soon as the step falls
+/// below `epsilon`, so the cap is never reached by a history that converges and
+/// raising it costs nothing. Measured on a history that needs four sweeps:
+///
+/// ```text
+/// max_iter      30: 4 iterations, 129.9 us
+/// max_iter 100_000: 4 iterations, 131.9 us
+/// ```
+///
+/// This was `30` until it was measured, and 30 truncated ordinary healthy
+/// histories: 160 events over 100 competitors already needs 42. Because a short
+/// fit is finite and sensibly ordered, that was invisible.
+///
+/// # Why it is not scaled to the history
+///
+/// The obvious improvement — pick the cap from the node or event count — does
+/// not work, because iteration count is driven by how *loopy* the graph is
+/// rather than how big it is. At a fixed 320 events over 40 slices, varying
+/// only the number of competitors sharing them:
+///
+/// ```text
+/// competitors    appearances each    iterations
+///           3                 213         2_789
+///          10                  64         1_068
+///          50                12.8           206
+///         100                 6.4            90
+///         400                 1.6             2
+/// ```
+///
+/// Three orders of magnitude apart on identical event and slice counts. Any
+/// formula in those two numbers would be badly wrong on some real shape, so the
+/// cap is a single value set high enough that reaching it means the fit is
+/// oscillating rather than merely large.
+///
+/// Reaching it is [`InferenceError::NotConverged`]. See
+/// [`History::converge`](crate::History::converge).
+pub const ITERATIONS: usize = 10_000;
 
 /// Largest team count `History::predict_outcome` will enumerate.
 ///
@@ -214,24 +250,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;
@@ -247,7 +306,7 @@ 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));
 
     // The leading coefficient is NEGATIVE. `rational - t` is negative here, so
     // a positive coefficient mirrors the starting point to `-x0` — the
@@ -262,7 +321,7 @@ fn erfc_inv(mut y: f64) -> f64 {
     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 }
@@ -291,7 +350,7 @@ pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
 /// 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` itself holds ~1e-7 *relative* accuracy down
+/// 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.
 ///
@@ -313,7 +372,7 @@ fn erfcx(x: f64) -> f64 {
         // 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.
-        (x * x).exp() * erfc(x)
+        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,
@@ -334,7 +393,7 @@ fn erfcx(x: f64) -> f64 {
 /// 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;
-    -(SQRT_TAU * sigma).ln() - 0.5 * z * z
+    -libm::log(SQRT_TAU * sigma) - 0.5 * z * z
 }
 
 /// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
@@ -347,10 +406,10 @@ pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
 
     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 + erfcx(z / SQRT_2).ln()
+        -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.
-        sf(x, mu, sigma).ln()
+        libm::log(sf(x, mu, sigma))
     }
 }
 
@@ -376,26 +435,24 @@ pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
     } 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 (cdf(hi, mu, sigma) - cdf(lo, mu, sigma))
-            .max(f64::MIN_POSITIVE)
-            .ln();
+        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 = (a * a - b * b).exp();
+    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 + bracket.ln()
+    -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
 }
@@ -474,7 +531,7 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
         let (v, u) = if alpha > 0.0 {
             // beta > alpha > 0, so this ratio of exponentials is at most 1 and
             // cannot overflow.
-            let scale = (0.5 * (alpha * alpha - beta * beta)).exp();
+            let scale = libm::exp(0.5 * (alpha * alpha - beta * beta));
             let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));
 
             (
@@ -662,7 +719,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)]
@@ -677,9 +734,9 @@ mod tests {
     }
 
     /// Upper-tail values of the standard normal, from published tables. The
-    /// point is not the digits — `erfc` only carries ~1e-7 relative — but that
-    /// a number comes back at all: `1 - cdf` returned exactly zero for every
-    /// one of these.
+    /// 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 [
@@ -691,7 +748,7 @@ mod tests {
             let got = sf(z, 0.0, 1.0);
             assert!(got > 0.0, "sf({z}) collapsed to zero");
             assert!(
-                (got - expected).abs() / expected < 1e-6,
+                (got - expected).abs() / expected < 1e-6, // published table values, 7 digits
                 "sf({z}) = {got}, expected ~{expected}"
             );
             assert_eq!(
@@ -709,11 +766,8 @@ mod tests {
         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);
-            // Bounded by `erfc`'s own ~1e-7 relative error, not by the
-            // subtraction: the two forms evaluate `erfc` at different points
-            // and the approximation is not exactly antisymmetric.
             assert!(
-                (naive - direct).abs() < 1e-6,
+                (naive - direct).abs() < 1e-15,
                 "z={z}: naive {naive} vs direct {direct}"
             );
         }
@@ -723,9 +777,7 @@ mod tests {
     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);
-            // `erfc(z) + erfc(-z) == 2` only to the accuracy of the
-            // approximation, which is ~1e-7 relative.
-            assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}");
+            assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}");
         }
     }
 
@@ -736,7 +788,7 @@ mod tests {
             let direct = (x * x).exp() * erfc(x);
             let scaled = erfcx(x);
             assert!(
-                (direct - scaled).abs() / scaled < 1e-6,
+                (direct - scaled).abs() / scaled < 1e-14,
                 "x={x}: direct {direct} vs erfcx {scaled}"
             );
         }
@@ -838,7 +890,7 @@ mod tests {
         ] {
             let got = SQRT_2 * erfc_inv(1.0 - p);
             assert!(
-                (got - exact).abs() / exact < 1e-6,
+                (got - exact).abs() / exact < 1e-14,
                 "p={p}: got {got}, exact {exact}"
             );
         }
@@ -847,6 +899,30 @@ mod tests {
     /// 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;
@@ -872,7 +948,7 @@ mod tests {
                 // 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-6,
+                    (recovered - p_draw).abs() < 1e-14,
                     "p_draw={p_draw} sd={sd}: recovered {recovered}"
                 );
             }
@@ -922,14 +998,10 @@ mod tests {
                 lp.exp()
             );
 
-            // Bounded by `erfc`'s ~1e-7, not tighter: above x = 2 `erfcx` uses
-            // a continued fraction accurate to ~1e-15, so the log path is the
-            // *more* accurate of the two and they part company at `erfc`'s
-            // error rather than at round-off.
             let ls = ln_sf(z, 0.5, 2.0);
             let direct = sf(z, 0.5, 2.0);
             assert!(
-                (ls.exp() - direct).abs() <= 1e-6 * direct.max(1e-300),
+                (ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300),
                 "ln_sf at {z}: {} vs {direct}",
                 ls.exp()
             );
@@ -941,11 +1013,8 @@ mod tests {
         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();
-            // See `log_space_helpers_agree_with_the_linear_forms_in_range`:
-            // the gap here is `erfc`'s own error, and the log path is the more
-            // accurate side of it.
             assert!(
-                (logged - direct).abs() <= 1e-6 * direct,
+                (logged - direct).abs() <= 1e-13 * direct,
                 "mu={mu}: {logged} vs {direct}"
             );
         }

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

@@ -33,17 +33,23 @@ pub(crate) const MAX_TEAMS_FOR_DISTRIBUTION: usize = 6;
 
 /// Relative tolerance for the first-place integrals.
 ///
-/// Tightening past this buys nothing: the underlying `cdf` is a rational
-/// approximation with fractional error ~1.2e-7, which contributes ~6e-9 to a
-/// finished probability and dominates any further quadrature refinement.
+/// 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). Measured against the exact two-team
-/// closed form, 2_048 nodes leave ~1.2e-6 of discretisation error while 8_192
-/// reach ~1e-7 — at which point the residual is the `cdf` rational
-/// approximation (~2.4e-8), not the grid, and refining further buys nothing.
+/// 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;
 
@@ -62,7 +68,7 @@ fn phi(z: f64) -> f64 {
 fn density(g: Gaussian, x: f64) -> f64 {
     let sigma = g.sigma();
     let z = (x - g.mu()) / sigma;
-    (-0.5 * z * z).exp() / (sigma * (2.0 * std::f64::consts::PI).sqrt())
+    libm::exp(-0.5 * z * z) / (sigma * (2.0 * std::f64::consts::PI).sqrt())
 }
 
 /// Per-pair draw margins.
@@ -503,7 +509,7 @@ mod tests {
 
     /// 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().powi(2) + b.sigma().powi(2)).sqrt();
+        let sd = a.sigma().hypot(b.sigma());
         (
             phi((a.mu() - b.mu() - eps) / sd),
             phi((b.mu() - a.mu() - eps) / sd),
@@ -523,7 +529,7 @@ mod tests {
             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-7 && (got[1] - wb).abs() < 1e-7,
+                (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}]"
             );
         }

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,78 @@ pub(crate) fn compute_elapsed<T: Time>(last: Option<&T>, current: &T) -> i64 {
     elapsed.max(0)
 }
 
+impl<T: Time> TimeSlice<T> {
+    /// This slice's scored event factors, as contrasts over competitors.
+    ///
+    /// Message passing produces per-competitor marginals and throws the
+    /// correlation away — `Item::likelihood` is already the projection of an
+    /// event's factor onto one competitor. So a 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, so only the second
+    /// moment needs rebuilding.
+    ///
+    /// Each entry is a contrast and the observation variance that sits on it.
+    /// Ranked events contribute nothing: their truncation factors are EP
+    /// approximations that inference does not retain.
+    pub(crate) fn scored_contrasts<D: Drift<T>>(
+        &self,
+        agents: &CompetitorStore<T, D>,
+    ) -> Vec<(Vec<(Index, f64)>, f64)> {
+        let mut out = Vec::new();
+
+        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: Vec<usize> = (0..event.teams.len()).collect();
+            order.sort_by(|&a, &b| {
+                event.teams[b]
+                    .output
+                    .partial_cmp(&event.teams[a].output)
+                    .unwrap_or(std::cmp::Ordering::Equal)
+            });
+
+            for pair in order.windows(2) {
+                let (hi, lo) = (pair[0], pair[1]);
+                let mut contrast: Vec<(Index, f64)> = Vec::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];
+                        noise += w * w * agents[item.agent].rating.beta.powi(2);
+                        contrast.push((item.agent, sign * w));
+                    }
+                }
+
+                out.push((contrast, noise));
+            }
+        }
+
+        out
+    }
+
+    /// True when every event here is scored, so the joint is exact.
+    pub(crate) fn all_scored(&self) -> bool {
+        self.events
+            .iter()
+            .all(|e| matches!(e.kind, EventKind::Scored { .. }))
+    }
+
+    /// The competitors appearing in this slice, with the elapsed count since
+    /// each one's previous appearance.
+    pub(crate) fn appearances(&self) -> impl Iterator<Item = (Index, i64)> + '_ {
+        self.skills
+            .keys()
+            .map(|idx| (idx, self.skills.get(idx).expect("slice key").elapsed))
+    }
+}
+
 #[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,7 +201,7 @@ 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"]]).unwrap();
     assert!(q > 0.0 && q <= 1.0);
@@ -217,7 +217,7 @@ 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"]]).unwrap();
     let wins = p.win_probabilities();
@@ -245,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

@@ -64,13 +64,13 @@ fn a_prior_applies_to_a_new_competitor() {
     let mut with = history();
     with.add_events(vec![bout("a", "b", 0, Some(seeded), None)])
         .unwrap();
-    with.converge().unwrap();
+    let _ = with.converge().unwrap();
 
     let mut without = history();
     without
         .add_events(vec![bout("a", "b", 0, None, None)])
         .unwrap();
-    without.converge().unwrap();
+    let _ = without.converge().unwrap();
 
     assert!(
         (skill_of(&with, "a").mu() - skill_of(&without, "a").mu()).abs() > 1.0,
@@ -91,7 +91,7 @@ fn a_prior_applies_to_a_competitor_the_history_already_knows() {
     // "a" now exists. Configuring it here used to do nothing whatsoever.
     late.add_events(vec![bout("a", "b", 1, Some(seeded), None)])
         .unwrap();
-    late.converge().unwrap();
+    let _ = late.converge().unwrap();
 
     let mut never = history();
     never
@@ -100,7 +100,7 @@ fn a_prior_applies_to_a_competitor_the_history_already_knows() {
             bout("a", "b", 1, None, None),
         ])
         .unwrap();
-    never.converge().unwrap();
+    let _ = never.converge().unwrap();
 
     assert!(
         (skill_of(&late, "a").mu() - skill_of(&never, "a").mu()).abs() > 1.0,
@@ -122,7 +122,7 @@ fn a_prior_is_whole_history_scoped_not_per_event() {
         .unwrap();
     late.add_events(vec![bout("a", "b", 1, Some(seeded), None)])
         .unwrap();
-    late.converge().unwrap();
+    let _ = late.converge().unwrap();
 
     let mut early = history();
     early
@@ -131,7 +131,7 @@ fn a_prior_is_whole_history_scoped_not_per_event() {
             bout("a", "b", 1, Some(seeded), None),
         ])
         .unwrap();
-    early.converge().unwrap();
+    let _ = early.converge().unwrap();
 
     let (l, e) = (skill_of(&late, "a"), skill_of(&early, "a"));
     assert!(
@@ -150,7 +150,7 @@ fn repeating_the_same_prior_is_inert() {
         bout("a", "b", 1, None, None),
     ])
     .unwrap();
-    once.converge().unwrap();
+    let _ = once.converge().unwrap();
 
     let mut every_time = history();
     every_time
@@ -159,7 +159,7 @@ fn repeating_the_same_prior_is_inert() {
             bout("a", "b", 1, Some(seeded), None),
         ])
         .unwrap();
-    every_time.converge().unwrap();
+    let _ = every_time.converge().unwrap();
 
     let (o, e) = (skill_of(&once, "a"), skill_of(&every_time, "a"));
     assert!(
@@ -203,7 +203,7 @@ fn setting_one_field_late_leaves_the_other_alone() {
     // Only the scale this time — the prior above must survive.
     h.add_events(vec![bout("a", "b", 1, None, Some(0.5))])
         .unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let mut both_upfront = history();
     both_upfront
@@ -212,7 +212,7 @@ fn setting_one_field_late_leaves_the_other_alone() {
             bout("a", "b", 1, None, None),
         ])
         .unwrap();
-    both_upfront.converge().unwrap();
+    let _ = both_upfront.converge().unwrap();
 
     let (a, b) = (skill_of(&h, "a"), skill_of(&both_upfront, "a"));
     assert!(

tests/convergence_strictness.rs

@@ -0,0 +1,152 @@
+//! Stopping short of convergence is an error, not a flag on a success.
+//!
+//! A fit that hits `max_iter` is wrong by a little: every rating is finite,
+//! the ordering looks sensible, and nothing in the numbers says they were
+//! still moving. When that was `Ok` with `converged: false`, detecting it was
+//! opt-in and `let _ = h.converge()` was the natural way to opt out — which is
+//! how a real defect once hid in this crate's own suite.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, InferenceError, Member, Outcome, Team,
+};
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn duel(a: &'static str, b: &'static str, t: i64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([3.0, 1.0]),
+    }
+}
+
+fn capped(max_iter: usize) -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.5))
+        .convergence(ConvergenceOptions {
+            max_iter,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+fn fill(h: &mut H) {
+    h.add_events((1..=6).map(|t| duel("a", "b", t)).collect::<Vec<_>>())
+        .unwrap();
+}
+
+#[test]
+fn hitting_the_cap_is_an_error() {
+    let mut h = capped(1);
+    fill(&mut h);
+    let err = h.converge().unwrap_err();
+    match err {
+        InferenceError::NotConverged {
+            iterations,
+            final_step,
+            epsilon,
+        } => {
+            assert_eq!(iterations, 1);
+            assert!(
+                final_step.0 > epsilon || final_step.1 > epsilon,
+                "{final_step:?}"
+            );
+        }
+        other => panic!("expected NotConverged, got {other:?}"),
+    }
+}
+
+/// The message has to name what to do about it, since the fit looks fine.
+#[test]
+fn the_error_says_how_to_fix_it() {
+    let mut h = capped(1);
+    fill(&mut h);
+    let text = h.converge().unwrap_err().to_string();
+    assert!(text.contains("did not converge in 1 iterations"), "{text}");
+    assert!(text.contains("max_iter"), "{text}");
+    assert!(text.contains("alpha"), "{text}");
+}
+
+/// The escape hatch: a deliberately capped fit is still reachable.
+#[test]
+fn converge_partial_returns_the_short_fit() {
+    let mut h = capped(1);
+    fill(&mut h);
+    let report = h.converge_partial().unwrap();
+    assert_eq!(report.iterations, 1);
+    assert!(!report.converged);
+    assert!(h.current_skill(&"a").is_some());
+}
+
+/// Both agree when the fit does converge, so the strict path costs nothing.
+#[test]
+fn the_two_agree_on_a_converged_fit() {
+    let mut strict = capped(20_000);
+    fill(&mut strict);
+    let a = strict.converge().unwrap();
+
+    let mut partial = capped(20_000);
+    fill(&mut partial);
+    let b = partial.converge_partial().unwrap();
+
+    assert!(a.converged && b.converged);
+    assert_eq!(a.iterations, b.iterations);
+    assert_eq!(a.final_step, b.final_step);
+}
+
+/// The default cap must be high enough that an ordinary history clears it.
+/// At the old value of 30 this history stopped short and said nothing.
+#[test]
+fn the_default_cap_clears_an_ordinary_history() {
+    let mut h: History<i64, ConstantDrift, _, String> = History::builder_with_key()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.05))
+        .build();
+
+    let mut events = Vec::new();
+    for t in 0..20i64 {
+        for j in 0..8usize {
+            let k = (t as usize) * 8 + j;
+            events.push(Event {
+                time: t,
+                teams: smallvec![
+                    Team::with_members([Member::new(format!("p{}", k % 100))]),
+                    Team::with_members([Member::new(format!("p{}", (k + 37) % 100))]),
+                ],
+                outcome: Outcome::scores([3.0, 1.0]),
+            });
+        }
+    }
+    h.add_events(events).unwrap();
+
+    let report = h
+        .converge()
+        .expect("an ordinary history must converge by default");
+    assert!(
+        report.iterations > 30,
+        "needed {} sweeps",
+        report.iterations
+    );
+    assert!(report.iterations < trueskill_tt::ITERATIONS);
+}
+
+/// An empty history converges trivially rather than erroring.
+#[test]
+fn an_empty_history_converges() {
+    let mut h = capped(1);
+    let report = h.converge().unwrap();
+    assert!(report.converged);
+    assert_eq!(report.iterations, 0);
+}

tests/degenerate_inputs.rs

@@ -170,8 +170,9 @@ fn event_builder_rejects_a_weights_length_mismatch() {
 fn event_builder_weights_mismatch_leaves_the_history_untouched() {
     let mut h = History::default();
 
-    // Two teams, so ingestion would otherwise succeed — a one-team event is
-    // rejected for an unrelated reason and would pass this vacuously.
+    // Two teams, so ingestion would otherwise succeed. A one-team event is
+    // rejected as `NotEnoughTeams` before the weights are ever examined, so
+    // building this with one team would pass vacuously.
     let _ = h
         .event(1)
         .team(["a"])
@@ -351,7 +352,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 +369,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 +390,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 +417,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
 }
 
@@ -385,7 +385,7 @@ fn drift_scale_applies_when_set_after_first_appearance() {
         outcome: Outcome::winner(1, 2),
     }])
     .unwrap();
-    late.converge().unwrap();
+    let _ = late.converge().unwrap();
 
     let applied = curve(&late, "anchor");
     let pinned_from_the_start = curve(&fit(distant_pair(Some(0.0)), 25.0 / 300.0), "anchor");

tests/event_builder_members.rs

@@ -0,0 +1,193 @@
+//! `EventBuilder::members` must reach exactly what the typed path reaches.
+//!
+//! Before this existed, `EventBuilder` could set weights and nothing else, so
+//! `prior` and `drift_scale` were expressible only through `Event`/`Team`/
+//! `Member` + `add_events`. Which ingestion route a competitor arrived through
+//! decided whether it could be configured at all.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, Gaussian, History, InferenceError, Member, Outcome,
+    Team,
+};
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn history() -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.5))
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+const PRIOR: Gaussian = Gaussian::from_ms(3.0, 1.5);
+
+/// The contract that makes the escape hatch worth having: same configuration,
+/// same fit, bit for bit.
+#[test]
+fn members_matches_the_typed_path_exactly() {
+    let mut typed = history();
+    typed
+        .add_events(vec![Event {
+            time: 1,
+            teams: smallvec![
+                Team::with_members([Member::new("player")]),
+                Team::with_members([Member::new("layout_7")
+                    .with_drift_scale(0.0)
+                    .with_prior(PRIOR)]),
+            ],
+            outcome: Outcome::scores([5.0, 2.0]),
+        }])
+        .unwrap();
+    assert!(typed.converge().unwrap().converged);
+
+    let mut fluent = history();
+    fluent
+        .event(1)
+        .team(["player"])
+        .members([Member::new("layout_7")
+            .with_drift_scale(0.0)
+            .with_prior(PRIOR)])
+        .scores([5.0, 2.0])
+        .commit()
+        .unwrap();
+    assert!(fluent.converge().unwrap().converged);
+
+    for key in ["player", "layout_7"] {
+        let a = typed.current_skill(&key).unwrap();
+        let b = fluent.current_skill(&key).unwrap();
+        assert_eq!(a.pi(), b.pi(), "{key} pi");
+        assert_eq!(a.tau(), b.tau(), "{key} tau");
+    }
+}
+
+/// The configuration has to actually take effect, not merely round-trip: a
+/// competitor pinned with `drift_scale = 0.0` must not move across slices,
+/// where an unpinned one does.
+///
+/// The comparison is against a control rather than against a fixed epsilon.
+/// Pinned marginals are not bit-identical across slices — each slice combines
+/// its own forward and backward messages, so the arithmetic order differs and
+/// the last bit moves. What "pinned" promises is that no drift variance
+/// accumulates, and the control is what makes that measurable.
+#[test]
+fn a_drift_scale_set_through_members_is_applied() {
+    fn spread(h: &H, key: &'static str) -> f64 {
+        let curve = h.learning_curve(&key);
+        assert!(curve.len() >= 2, "{key}: expected several appearances");
+        let (lo, hi) = curve.iter().fold((f64::MAX, f64::MIN), |(lo, hi), (_, g)| {
+            (lo.min(g.sigma()), hi.max(g.sigma()))
+        });
+        (hi - lo) / hi
+    }
+
+    let mut h = history();
+    for t in 1..=4 {
+        h.event(t)
+            .team(["player"])
+            .members([Member::new("pinned").with_drift_scale(0.0)])
+            .scores([5.0, 2.0])
+            .commit()
+            .unwrap();
+        // Same shape, no pinning: the control.
+        h.event(t)
+            .team(["rival"])
+            .team(["drifting"])
+            .scores([5.0, 2.0])
+            .commit()
+            .unwrap();
+    }
+    assert!(h.converge().unwrap().converged);
+
+    let pinned = spread(&h, "pinned");
+    let drifting = spread(&h, "drifting");
+    assert!(pinned < 1e-9, "pinned competitor moved: {pinned:e}");
+    assert!(
+        drifting > 1e-3,
+        "control did not move, so the test proves nothing: {drifting:e}"
+    );
+}
+
+/// `weights` still applies to a team added through `members`, and still
+/// records a mismatch rather than partially applying it.
+#[test]
+fn weights_still_guards_a_members_team() {
+    let mut h = history();
+    let err = h
+        .event(1)
+        .team(["a"])
+        .members([Member::new("b"), Member::new("c")])
+        .weights([1.0])
+        .winner(0)
+        .commit()
+        .unwrap_err();
+    assert!(
+        matches!(
+            err,
+            InferenceError::MismatchedShape {
+                kind: "weights",
+                expected: 2,
+                got: 1
+            }
+        ),
+        "{err:?}"
+    );
+    assert!(h.current_skill(&"b").is_none(), "nothing may reach history");
+}
+
+/// An invalid `drift_scale` surfaces from `commit`, not from a panic and not
+/// silently.
+#[test]
+fn an_invalid_drift_scale_surfaces_from_commit() {
+    for bad in [-1.0, f64::NAN, f64::INFINITY] {
+        let mut h = history();
+        let err = h
+            .event(1)
+            .team(["a"])
+            .members([Member::new("b").with_drift_scale(bad)])
+            .winner(0)
+            .commit()
+            .unwrap_err();
+        assert!(
+            matches!(
+                err,
+                InferenceError::InvalidParameter {
+                    name: "drift_scale",
+                    ..
+                }
+            ),
+            "{bad}: {err:?}"
+        );
+        assert!(h.current_skill(&"b").is_none(), "{bad} reached the history");
+    }
+}
+
+/// `members` and `team` compose in either order.
+#[test]
+fn members_and_team_interleave() {
+    let mut h = history();
+    h.event(1)
+        .members([Member::new("a").with_prior(PRIOR)])
+        .team(["b"])
+        .scores([3.0, 1.0])
+        .commit()
+        .unwrap();
+    h.event(2)
+        .team(["b"])
+        .members([Member::new("c").with_prior(PRIOR)])
+        .scores([2.0, 4.0])
+        .commit()
+        .unwrap();
+    assert!(h.converge().unwrap().converged);
+    for key in ["a", "b", "c"] {
+        assert!(h.current_skill(&key).is_some(), "{key} missing");
+    }
+}

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/ingestion_shape.rs

@@ -0,0 +1,147 @@
+//! Malformed events must be rejected at the ingestion boundary.
+//!
+//! Every case here was reachable from safe public API in a release build. Two
+//! of them are the two shapes this crate's defects keep taking: a panic from
+//! deep inside inference, and a finite, plausible-looking posterior computed
+//! from an event that should never have been accepted.
+//!
+//! `InferenceError::NotEnoughTeams` and `EmptyTeam` already existed when these
+//! were found — they were checked on the prediction paths and nowhere else, so
+//! ingestion could still manufacture the states they describe.
+
+use smallvec::smallvec;
+use trueskill_tt::{Event, History, InferenceError, Member, Outcome, Team};
+
+type Ev = Event<i64, &'static str>;
+
+fn history() -> History<i64, trueskill_tt::ConstantDrift, trueskill_tt::NullObserver, &'static str>
+{
+    History::builder().score_sigma(1.0).build()
+}
+
+fn teams(names: &[&[&'static str]]) -> smallvec::SmallVec<[Team<&'static str>; 4]> {
+    names
+        .iter()
+        .map(|team| Team::with_members(team.iter().map(|k| Member::new(*k))))
+        .collect()
+}
+
+/// The regression this file exists for: `run_chain` builds one diff link per
+/// adjacent pair of teams, so a one-team event left it indexing `links[1..]`
+/// on an empty vector and panicked — in release, from `History::add_events`.
+#[test]
+fn a_one_team_event_is_an_error_not_a_panic() {
+    let mut h = history();
+    let err = h
+        .add_events(vec![Ev {
+            time: 1,
+            teams: teams(&[&["a"]]),
+            outcome: Outcome::winner(0, 1),
+        }])
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::NotEnoughTeams { got: 1 }),
+        "{err:?}"
+    );
+}
+
+#[test]
+fn a_zero_team_event_is_an_error() {
+    let mut h = history();
+    let err = h
+        .add_events(vec![Ev {
+            time: 1,
+            teams: smallvec![],
+            outcome: Outcome::ranking([]),
+        }])
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::NotEnoughTeams { got: 0 }),
+        "{err:?}"
+    );
+}
+
+/// The quiet half. An empty team contributes no performance, so before this
+/// was rejected the event converged and handed back a finite posterior for its
+/// opponent — a plausible constant computed from nothing.
+#[test]
+fn an_empty_team_is_an_error_rather_than_a_free_win() {
+    let mut h = history();
+    let err = h
+        .add_events(vec![Ev {
+            time: 1,
+            teams: teams(&[&[], &["b"]]),
+            outcome: Outcome::winner(0, 2),
+        }])
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::EmptyTeam { team: 0 }),
+        "{err:?}"
+    );
+    // Nothing was recorded, so the history is still empty.
+    assert!(h.current_skill(&"b").is_none());
+}
+
+#[test]
+fn an_empty_team_is_reported_by_position() {
+    let mut h = history();
+    let err = h
+        .add_events(vec![Ev {
+            time: 1,
+            teams: teams(&[&["a"], &[]]),
+            outcome: Outcome::winner(0, 2),
+        }])
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::EmptyTeam { team: 1 }),
+        "{err:?}"
+    );
+}
+
+/// A NaN score used to ingest cleanly. `converge` reported `NonFiniteResult`,
+/// but a caller who read `current_skill` first was handed `tau: NaN` with
+/// nothing to say so.
+#[test]
+fn a_non_finite_score_is_rejected_at_ingestion() {
+    for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
+        let mut h = history();
+        let err = h
+            .add_events(vec![Ev {
+                time: 1,
+                teams: teams(&[&["a"], &["b"]]),
+                outcome: Outcome::scores([bad, 0.0]),
+            }])
+            .unwrap_err();
+        assert!(
+            matches!(err, InferenceError::InvalidParameter { name: "score", .. }),
+            "{bad}: {err:?}"
+        );
+        assert!(h.current_skill(&"a").is_none(), "{bad} was recorded anyway");
+    }
+}
+
+/// The fluent builder routes through the same chokepoint, so it inherits the
+/// checks rather than needing its own.
+#[test]
+fn the_event_builder_inherits_the_shape_checks() {
+    let mut h = history();
+    let err = h.event(1).team(["a"]).winner(0).commit().unwrap_err();
+    assert!(
+        matches!(err, InferenceError::NotEnoughTeams { got: 1 }),
+        "{err:?}"
+    );
+}
+
+/// A well-formed event is untouched by any of this.
+#[test]
+fn a_well_formed_event_still_ingests() {
+    let mut h = history();
+    h.add_events(vec![Ev {
+        time: 1,
+        teams: teams(&[&["a"], &["b"]]),
+        outcome: Outcome::scores([3.0, 1.0]),
+    }])
+    .unwrap();
+    assert!(h.converge().unwrap().converged);
+    assert!(h.current_skill(&"a").unwrap().mu() > h.current_skill(&"b").unwrap().mu());
+}

tests/joint_handle.rs

@@ -0,0 +1,267 @@
+//! `History::joint` factorises once and answers many questions.
+//!
+//! The contract that matters is *identity*: a `Joint` must return exactly what
+//! the one-shot call returns, bit for bit. A faster path that quietly disagreed
+//! with the slow one would be worse than no fast path — a caller would get
+//! different numbers depending on how many questions they happened to ask.
+
+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 duel(a: &'static str, b: &'static str, t: i64, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+fn ranked(a: &'static str, b: &'static str, t: i64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::winner(0, 2),
+    }
+}
+
+fn history(unknown: UnknownKeys) -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.5))
+        .unknown_keys(unknown)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+/// Several slices, competitors with different last appearances, so `latest`
+/// and `at_slice` both have work to do.
+fn fitted(unknown: UnknownKeys) -> H {
+    let mut h = history(unknown);
+    h.add_events(vec![
+        duel("a", "b", 1, 5.0, 2.0),
+        duel("c", "d", 1, 3.0, 3.5),
+        duel("a", "c", 2, 6.0, 1.0),
+        duel("b", "d", 3, 4.0, 3.0),
+        duel("a", "d", 4, 7.0, 2.0),
+        duel("b", "c", 5, 2.0, 4.0),
+    ])
+    .unwrap();
+    let report = h.converge().unwrap();
+    assert!(report.converged, "fixture must converge");
+    h
+}
+
+const PAIRS: [(&str, &str); 6] = [
+    ("a", "b"),
+    ("a", "c"),
+    ("a", "d"),
+    ("b", "c"),
+    ("b", "d"),
+    ("c", "d"),
+];
+
+#[test]
+fn a_joint_answers_exactly_what_the_one_shot_call_does() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+
+    for (a, b) in PAIRS {
+        let terms = [(&a, 1.0), (&b, -1.0)];
+        let one_shot = h.posterior_of(&terms).unwrap();
+        let cached = joint.posterior_of(&terms).unwrap();
+        assert_eq!(one_shot.pi(), cached.pi(), "{a} - {b}");
+        assert_eq!(one_shot.tau(), cached.tau(), "{a} - {b}");
+    }
+}
+
+#[test]
+fn a_joint_agrees_at_a_pinned_time_too() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+
+    for time in 1..=5 {
+        for (a, b) in PAIRS {
+            let terms = [(&a, 1.0), (&b, -1.0)];
+            let one_shot = h.posterior_of_at(time, &terms);
+            let cached = joint.posterior_of_at(time, &terms);
+            match (one_shot, cached) {
+                (Ok(x), Ok(y)) => {
+                    assert_eq!(x.pi(), y.pi(), "t={time} {a} - {b}");
+                    assert_eq!(x.tau(), y.tau(), "t={time} {a} - {b}");
+                }
+                (Err(x), Err(y)) => assert_eq!(x, y, "t={time} {a} - {b}"),
+                (x, y) => panic!("t={time} {a} - {b}: disagreed on success: {x:?} vs {y:?}"),
+            }
+        }
+    }
+}
+
+#[test]
+fn a_joint_scores_candidate_matchups_identically() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    let (a, b) = ("a", "b");
+    let target = [(&a, 1.0), (&b, -1.0)];
+
+    for (x, y) in PAIRS {
+        let teams: [&[&&str]; 2] = [&[&x], &[&y]];
+        let one_shot = h.expected_variance_reduction(&teams, &target).unwrap();
+        let cached = joint.expected_variance_reduction(&teams, &target).unwrap();
+        assert_eq!(one_shot, cached, "{x} vs {y}");
+    }
+}
+
+/// The whole point: a competitor appears once per slice, so the joint is over
+/// appearances rather than competitors, and a caller sizing a batch needs to
+/// know which.
+#[test]
+fn variables_counts_appearances_not_competitors() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    // Four competitors, twelve appearances across five slices, all with
+    // positive drift between them, so no two collapse.
+    assert_eq!(joint.variables(), 12);
+}
+
+/// How much the collapse is worth, which is the part a caller has to plan
+/// around: a drift-free competitor contributes **one** variable however long
+/// the history, so the same events at `gamma = 0` and `gamma > 0` differ by
+/// roughly the slice count in problem size — and by its cube in solve time.
+///
+/// Reported by a consumer as an 8x difference in solve time on a ~2,000-node,
+/// 76-slice model (787 ms career against 6,214 ms drifting). This pins the
+/// mechanism behind that so a change to the collapse rule cannot quietly
+/// remove it.
+#[test]
+fn drift_free_competitors_shrink_the_joint_by_the_slice_count() {
+    fn variables(gamma: f64) -> usize {
+        let mut h = History::builder()
+            .mu(0.0)
+            .sigma(6.0)
+            .beta(1.0)
+            .score_sigma(2.0)
+            .drift(ConstantDrift(gamma))
+            .convergence(ConvergenceOptions {
+                max_iter: 20_000,
+                epsilon: 1e-13,
+                alpha: 1.0,
+            })
+            .build();
+        h.add_events(
+            (1..=10)
+                .map(|t| duel("a", "b", t, 5.0, 2.0))
+                .collect::<Vec<_>>(),
+        )
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h.joint().unwrap().variables()
+    }
+
+    let drifting = variables(0.5);
+    let career = variables(0.0);
+
+    // Two competitors over ten slices: twenty appearances, or two variables.
+    assert_eq!(drifting, 20);
+    assert_eq!(career, 2);
+    assert_eq!(
+        drifting / career,
+        10,
+        "collapse should track the slice count"
+    );
+}
+
+/// With `drift = 0` consecutive appearances are the same latent variable, so
+/// the joint is smaller than the appearance count.
+#[test]
+fn pinned_competitors_collapse_consecutive_appearances() {
+    let mut h = 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-13,
+            alpha: 1.0,
+        })
+        .build();
+    h.add_events(vec![
+        duel("a", "b", 1, 5.0, 2.0),
+        duel("a", "b", 2, 4.0, 3.0),
+        duel("a", "b", 3, 6.0, 1.0),
+    ])
+    .unwrap();
+    assert!(h.converge().unwrap().converged);
+    assert_eq!(h.joint().unwrap().variables(), 2);
+}
+
+#[test]
+fn a_ranked_history_has_no_exact_joint() {
+    let mut h = history(UnknownKeys::Reject);
+    h.add_events(vec![duel("a", "b", 1, 5.0, 2.0), ranked("a", "b", 2)])
+        .unwrap();
+    let _ = h.converge().unwrap();
+    assert!(matches!(
+        h.joint().unwrap_err(),
+        InferenceError::JointUnavailable { .. }
+    ));
+}
+
+#[test]
+fn an_empty_history_has_no_joint() {
+    let h = history(UnknownKeys::Reject);
+    assert!(matches!(
+        h.joint().unwrap_err(),
+        InferenceError::JointUnavailable { .. }
+    ));
+}
+
+/// Unknown keys are decided per query, not when the joint is factorised — the
+/// factorisation does not depend on the question.
+#[test]
+fn unknown_keys_are_rejected_per_query() {
+    let h = fitted(UnknownKeys::Reject);
+    let joint = h.joint().unwrap();
+    let (a, z) = ("a", "nobody");
+    assert!(matches!(
+        joint.posterior_of(&[(&a, 1.0), (&z, -1.0)]).unwrap_err(),
+        InferenceError::UnknownKey { .. }
+    ));
+    // The handle is still usable afterwards.
+    let b = "b";
+    assert!(joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).is_ok());
+}
+
+/// Under `Prior`, an unseen competitor is independent of everything in the
+/// history, and the cached path must add the same prior variance the one-shot
+/// path does.
+#[test]
+fn unseen_competitors_match_the_one_shot_path() {
+    let h = fitted(UnknownKeys::Prior);
+    let joint = h.joint().unwrap();
+    let (a, z) = ("a", "nobody");
+    let terms = [(&a, 1.0), (&z, -1.0)];
+    let one_shot = h.posterior_of(&terms).unwrap();
+    let cached = joint.posterior_of(&terms).unwrap();
+    assert_eq!(one_shot.pi(), cached.pi());
+    assert_eq!(one_shot.tau(), cached.tau());
+}

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

@@ -42,7 +42,7 @@ fn every_observer_callback_fires() {
     h.record_winner(&"a", &"b", 1).unwrap();
     h.record_winner(&"b", &"c", 2).unwrap();
     h.record_winner(&"c", &"a", 3).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     assert!(
         !recorder.iterations.lock().unwrap().is_empty(),
@@ -65,7 +65,7 @@ fn slice_callbacks_report_the_slice_they_swept() {
 
     h.record_winner(&"a", &"b", 10).unwrap();
     h.record_winner(&"a", &"b", 20).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let slices = recorder.slices.lock().unwrap();
 
@@ -93,7 +93,7 @@ fn a_single_slice_history_still_reports_its_sweep() {
     let mut h = History::builder().observer(Arc::clone(&recorder)).build();
 
     h.record_winner(&"a", &"b", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let slices = recorder.slices.lock().unwrap();
     assert!(
@@ -112,7 +112,7 @@ fn a_shared_observer_reaches_the_callers_handle() {
     let mut h = History::builder().observer(Arc::clone(&recorder)).build();
 
     h.record_winner(&"a", &"b", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     assert!(!recorder.iterations.lock().unwrap().is_empty());
     assert!(!recorder.slices.lock().unwrap().is_empty());
@@ -125,12 +125,12 @@ 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();
-    h.converge().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();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 }
 
 /// A non-shared observer can be reclaimed after convergence instead.
@@ -138,7 +138,7 @@ fn a_trait_object_observer_works() {
 fn into_observer_returns_the_accumulated_state() {
     let mut h = History::builder().observer(Recorder::default()).build();
     h.record_winner(&"a", &"b", 1).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     // Readable in place...
     assert!(!h.observer().iterations.lock().unwrap().is_empty());
@@ -155,7 +155,7 @@ fn a_borrowed_observer_works() {
     {
         let mut h = History::builder().observer(&recorder).build();
         h.record_winner(&"a", &"b", 1).unwrap();
-        h.converge().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

@@ -9,7 +9,7 @@ fn history_with(names: &[&'static str], p_draw: f64) -> History {
     for pair in names.windows(2) {
         h.record_winner(&pair[0], &pair[1], 1).unwrap();
     }
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
     h
 }
 
@@ -20,7 +20,14 @@ fn unknown_keys_are_reported_not_silently_dropped() {
     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 });
+    assert_eq!(
+        err,
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"ghost\"".to_owned(),
+        }
+    );
 
     // Every prediction entry point, not just one.
     assert!(
@@ -35,7 +42,14 @@ fn unknown_keys_are_reported_not_silently_dropped() {
 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 });
+    assert_eq!(
+        err,
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"x\"".to_owned(),
+        }
+    );
 }
 
 #[test]
@@ -184,7 +198,7 @@ fn the_stronger_competitor_is_favoured() {
     for t in 1..=10 {
         h.record_winner(&"strong", &"weak", t).unwrap();
     }
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let p = h.predict_outcome(&[&[&"strong"], &[&"weak"]]).unwrap();
     let (best, _) = p.most_likely().expect("a most likely outcome");
@@ -206,7 +220,7 @@ fn team_size_affects_the_prediction() {
         .winner(0)
         .commit()
         .unwrap();
-    h.converge().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());
@@ -229,7 +243,7 @@ fn information_gain_prefers_the_uncertain_pairing() {
         h.record_winner(&"rival", &"known", t + 100).unwrap();
     }
     h.record_winner(&"known", &"newcomer", 500).unwrap();
-    h.converge().unwrap();
+    let _ = h.converge().unwrap();
 
     let settled = h
         .expected_information_gain(&[&[&"known"], &[&"rival"]])
@@ -271,7 +285,11 @@ fn information_gain_reports_unknown_keys() {
     assert_eq!(
         h.expected_information_gain(&[&[&"a"], &[&"ghost"]])
             .unwrap_err(),
-        InferenceError::UnknownKey { team: 1, member: 0 }
+        InferenceError::UnknownKey {
+            team: 1,
+            member: 0,
+            key: "\"ghost\"".to_owned(),
+        }
     );
 }
 
@@ -289,3 +307,111 @@ fn information_gain_accounts_for_draws() {
     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,7 +108,7 @@ 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"]]).unwrap();
     assert!(

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/registration.rs

@@ -0,0 +1,338 @@
+//! Configuring a competitor before anything is observed about them.
+//!
+//! The configuration a competitor needs is usually a property of the domain —
+//! "every layout is static" — not of whichever event happens to mention them
+//! first. Stating it per-event meant every ingestion path had to remember it,
+//! and two of the four paths could not state it at all.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, Gaussian, History, InferenceError, Member, Outcome,
+    Team,
+};
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+const PINNED: Gaussian = Gaussian::from_ms(2.0, 0.5);
+
+fn history() -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(6.0)
+        .beta(1.0)
+        .score_sigma(2.0)
+        .drift(ConstantDrift(0.5))
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+fn duel(
+    a: &'static str,
+    b: &'static str,
+    t: i64,
+    m: Option<Member<&'static str>>,
+) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([m.unwrap_or_else(|| Member::new(b))]),
+        ],
+        outcome: Outcome::scores([5.0, 2.0]),
+    }
+}
+
+fn skills(h: &H) -> Vec<(&'static str, Gaussian)> {
+    ["player", "layout"]
+        .into_iter()
+        .map(|k| (k, h.current_skill(&k).unwrap()))
+        .collect()
+}
+
+/// The headline contract.
+#[test]
+fn registering_matches_configuring_on_the_first_event() {
+    let configured = {
+        let mut h = history();
+        h.add_events(vec![
+            duel(
+                "player",
+                "layout",
+                1,
+                Some(
+                    Member::new("layout")
+                        .with_drift_scale(0.0)
+                        .with_prior(PINNED),
+                ),
+            ),
+            duel("player", "layout", 2, None),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h
+    };
+
+    let registered = {
+        let mut h = history();
+        h.register(
+            Member::new("layout")
+                .with_drift_scale(0.0)
+                .with_prior(PINNED),
+        )
+        .unwrap();
+        h.add_events(vec![
+            duel("player", "layout", 1, None),
+            duel("player", "layout", 2, None),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h
+    };
+
+    for ((k, a), (_, b)) in skills(&configured).into_iter().zip(skills(&registered)) {
+        assert_eq!(a.pi(), b.pi(), "{k} pi");
+        assert_eq!(a.tau(), b.tau(), "{k} tau");
+    }
+}
+
+/// The case `EventBuilder` and the typed path cannot reach: a competitor whose
+/// first appearance arrives through the two-argument convenience route.
+#[test]
+fn registration_reaches_a_competitor_first_seen_through_record_winner() {
+    let mut h = history();
+    h.register(
+        Member::new("layout")
+            .with_drift_scale(0.0)
+            .with_prior(PINNED),
+    )
+    .unwrap();
+    h.record_winner(&"player", &"layout", 1).unwrap();
+    h.record_winner(&"player", &"layout", 2).unwrap();
+    let _ = h.converge().unwrap();
+
+    let rating = h.rating(&"layout").unwrap();
+    assert_eq!(rating.drift_scale(), 0.0);
+    assert_eq!(rating.prior().mu(), PINNED.mu());
+
+    // Pinned means pinned: no drift across the two slices.
+    let curve = h.learning_curve(&"layout");
+    assert!(curve.len() >= 2);
+    let widest = curve
+        .iter()
+        .map(|(_, g)| g.sigma())
+        .fold(f64::MIN, f64::max);
+    let narrowest = curve
+        .iter()
+        .map(|(_, g)| g.sigma())
+        .fold(f64::MAX, f64::min);
+    assert!(
+        (widest - narrowest) / widest < 1e-9,
+        "{narrowest} .. {widest}"
+    );
+}
+
+#[test]
+fn registering_a_known_competitor_is_an_error() {
+    let mut h = history();
+    h.record_winner(&"player", &"layout", 1).unwrap();
+    let err = h.register(Member::new("layout")).unwrap_err();
+    assert!(
+        matches!(err, InferenceError::AlreadyRegistered { .. }),
+        "{err:?}"
+    );
+}
+
+#[test]
+fn registering_twice_is_an_error() {
+    let mut h = history();
+    h.register(Member::new("layout").with_drift_scale(0.0))
+        .unwrap();
+    let err = h
+        .register(Member::new("layout").with_drift_scale(1.0))
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::AlreadyRegistered { .. }),
+        "{err:?}"
+    );
+    // The first registration stands.
+    assert_eq!(h.rating(&"layout").unwrap().drift_scale(), 0.0);
+}
+
+/// `weight` is per-event and meaningless here, so it is rejected rather than
+/// dropped — dropping it silently is the defect class this whole area keeps
+/// producing.
+#[test]
+fn a_weight_on_a_registration_is_rejected() {
+    let mut h = history();
+    let err = h
+        .register(Member::new("layout").with_weight(0.5))
+        .unwrap_err();
+    assert!(
+        matches!(err, InferenceError::InvalidParameter { name: "weight", .. }),
+        "{err:?}"
+    );
+}
+
+#[test]
+fn an_invalid_drift_scale_on_a_registration_is_rejected() {
+    for bad in [-1.0, f64::NAN, f64::INFINITY] {
+        let mut h = history();
+        let err = h
+            .register(Member::new("layout").with_drift_scale(bad))
+            .unwrap_err();
+        assert!(
+            matches!(
+                err,
+                InferenceError::InvalidParameter {
+                    name: "drift_scale",
+                    ..
+                }
+            ),
+            "{bad}: {err:?}"
+        );
+    }
+}
+
+/// Registration makes the fit independent of the order events arrive in,
+/// which is what the per-event shape could not guarantee.
+#[test]
+fn registration_makes_the_fit_order_independent() {
+    let build = |reversed: bool| {
+        let mut h = history();
+        h.register(
+            Member::new("layout")
+                .with_drift_scale(0.0)
+                .with_prior(PINNED),
+        )
+        .unwrap();
+        let mut events = vec![
+            duel("player", "layout", 1, None),
+            duel("player", "layout", 2, None),
+            duel("player", "layout", 3, None),
+        ];
+        if reversed {
+            events.reverse();
+        }
+        h.add_events(events).unwrap();
+        let _ = h.converge().unwrap();
+        h
+    };
+
+    let forward = build(false);
+    let backward = build(true);
+    for ((k, a), (_, b)) in skills(&forward).into_iter().zip(skills(&backward)) {
+        assert_eq!(a.pi(), b.pi(), "{k} pi");
+        assert_eq!(a.tau(), b.tau(), "{k} tau");
+    }
+}
+
+/// `rating` is the read-back that made a configuration mistake detectable from
+/// outside the crate at all. Every other accessor reports what inference
+/// inferred; this reports what it was told.
+#[test]
+fn rating_reads_back_what_was_stored() {
+    let mut h = history();
+    assert!(h.rating(&"nobody").is_none());
+
+    h.register(
+        Member::new("layout")
+            .with_drift_scale(0.25)
+            .with_prior(PINNED),
+    )
+    .unwrap();
+    let r = h.rating(&"layout").unwrap();
+    assert_eq!(r.drift_scale(), 0.25);
+    assert_eq!(r.prior().pi(), PINNED.pi());
+    assert_eq!(r.prior().tau(), PINNED.tau());
+
+    // A competitor created by an event reports the history defaults.
+    h.record_winner(&"player", &"layout", 1).unwrap();
+    assert_eq!(h.rating(&"player").unwrap().drift_scale(), 1.0);
+}
+
+/// The decision this issue turned on: two different values for one competitor
+/// are an error whether they arrive in one batch or two.
+///
+/// Last-write-wins across batches cut against the invariant
+/// `tests/ingestion_equivalence.rs` protects — the same contradictory events
+/// errored when batched and succeeded, order-dependently, one at a time.
+mod conflicting_configuration {
+    use super::*;
+
+    fn seed(scale: f64) -> Event<i64, &'static str> {
+        duel(
+            "player",
+            "layout",
+            1,
+            Some(Member::new("layout").with_drift_scale(scale)),
+        )
+    }
+
+    #[test]
+    fn within_one_batch_is_an_error() {
+        let mut h = history();
+        let err = h.add_events(vec![seed(0.0), seed(1.0)]).unwrap_err();
+        assert!(
+            matches!(
+                err,
+                InferenceError::ConflictingCompetitorConfig {
+                    field: "drift_scale",
+                    ..
+                }
+            ),
+            "{err:?}"
+        );
+    }
+
+    #[test]
+    fn across_two_batches_is_also_an_error() {
+        let mut h = history();
+        h.add_events(vec![seed(0.0)]).unwrap();
+        let err = h.add_events(vec![seed(1.0)]).unwrap_err();
+        assert!(
+            matches!(
+                err,
+                InferenceError::ConflictingCompetitorConfig {
+                    field: "drift_scale",
+                    ..
+                }
+            ),
+            "{err:?}"
+        );
+        // Rejected before anything mutates: the first declaration stands.
+        assert_eq!(h.rating(&"layout").unwrap().drift_scale(), 0.0);
+    }
+
+    /// Repeating the *same* value stays inert, which is the expected shape
+    /// when the configuration is a property of the domain.
+    #[test]
+    fn repeating_the_same_value_is_inert() {
+        let mut h = history();
+        h.add_events(vec![seed(0.0)]).unwrap();
+        h.add_events(vec![seed(0.0)]).unwrap();
+        assert_eq!(h.rating(&"layout").unwrap().drift_scale(), 0.0);
+    }
+
+    /// A registration and a later event that agree are fine; one that
+    /// disagrees is the same error.
+    #[test]
+    fn a_registration_conflicts_with_a_later_event() {
+        let mut h = history();
+        h.register(Member::new("layout").with_drift_scale(0.0))
+            .unwrap();
+        h.add_events(vec![seed(0.0)]).unwrap();
+
+        let mut h2 = history();
+        h2.register(Member::new("layout").with_drift_scale(0.0))
+            .unwrap();
+        let err = h2.add_events(vec![seed(1.0)]).unwrap_err();
+        assert!(
+            matches!(err, InferenceError::ConflictingCompetitorConfig { .. }),
+            "{err:?}"
+        );
+    }
+}

tests/time_expanded_joint.rs

@@ -0,0 +1,296 @@
+//! The joint must span slices, because Through Time reads each competitor at
+//! their own last appearance.
+//!
+//! The exact posterior of a multi-slice scored history is still Gaussian: the
+//! prior, the drift between appearances, and the scored likelihoods are all
+//! Gaussian. So it can be written out by hand and compared against, which is
+//! the check a single-slice fixture cannot make.
+
+use smallvec::smallvec;
+use trueskill_tt::{
+    ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team, UnknownKeys,
+};
+
+const SIGMA0: f64 = 6.0;
+const BETA: f64 = 1.0;
+const SCORE_SIGMA: f64 = 2.0;
+const GAMMA: f64 = 0.5;
+
+type H = History<i64, ConstantDrift, trueskill_tt::NullObserver, &'static str>;
+
+fn history(gamma: f64) -> H {
+    History::builder()
+        .mu(0.0)
+        .sigma(SIGMA0)
+        .beta(BETA)
+        .score_sigma(SCORE_SIGMA)
+        .drift(ConstantDrift(gamma))
+        .unknown_keys(UnknownKeys::Reject)
+        .convergence(ConvergenceOptions {
+            max_iter: 20_000,
+            epsilon: 1e-13,
+            alpha: 1.0,
+        })
+        .build()
+}
+
+fn duel(a: &'static str, b: &'static str, t: i64, sa: f64, sb: f64) -> Event<i64, &'static str> {
+    Event {
+        time: t,
+        teams: smallvec![
+            Team::with_members([Member::new(a)]),
+            Team::with_members([Member::new(b)]),
+        ],
+        outcome: Outcome::scores([sa, sb]),
+    }
+}
+
+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| f64::from(u8::from(i == j))).collect())
+        .collect();
+    for col in 0..n {
+        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
+}
+
+/// Two competitors, two slices ten units apart, one duel in each.
+///
+/// The exact precision is written out explicitly here rather than obtained
+/// from the crate, so this is an independent check rather than a restatement.
+/// Variables are `[a0, b0, a1, b1]`.
+#[test]
+fn a_two_slice_joint_matches_the_exact_posterior() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 10, 4.0, 3.0),
+    ])
+    .unwrap();
+    let report = h.converge().unwrap();
+    assert!(report.converged, "{:?}", report.final_step);
+
+    let prior_prec = 1.0 / (SIGMA0 * SIGMA0);
+    let drift_prec = 1.0 / (10.0 * GAMMA * GAMMA);
+    let obs_prec = 1.0 / (SCORE_SIGMA * SCORE_SIGMA + 2.0 * BETA * BETA);
+
+    let mut lambda = vec![vec![0.0; 4]; 4];
+    // priors on the first appearances
+    lambda[0][0] += prior_prec;
+    lambda[1][1] += prior_prec;
+    // drift a0-a1 and b0-b1
+    for (p, q) in [(0usize, 2usize), (1, 3)] {
+        lambda[p][p] += drift_prec;
+        lambda[q][q] += drift_prec;
+        lambda[p][q] -= drift_prec;
+        lambda[q][p] -= drift_prec;
+    }
+    // one duel per slice: contrast (+1, -1) on that slice's variables
+    for (p, q) in [(0usize, 1usize), (2, 3)] {
+        lambda[p][p] += obs_prec;
+        lambda[q][q] += obs_prec;
+        lambda[p][q] -= obs_prec;
+        lambda[q][p] -= obs_prec;
+    }
+    let cov = inverse(lambda);
+
+    // The crate reads each competitor at their latest appearance: a1, b1.
+    let exact_gap = (cov[2][2] + cov[3][3] - 2.0 * cov[2][3]).sqrt();
+    let got = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    assert!(
+        (got.sigma() - exact_gap).abs() / exact_gap < 1e-9,
+        "difference: got {} exact {exact_gap}",
+        got.sigma()
+    );
+
+    let exact_single = cov[2][2].sqrt();
+    let got_single = h.posterior_of(&[(&"a", 1.0)]).unwrap();
+    assert!(
+        (got_single.sigma() - exact_single).abs() / exact_single < 1e-9,
+        "single node: got {} exact {exact_single}",
+        got_single.sigma()
+    );
+}
+
+/// The case that motivated this: competitors read at *different* slices, with
+/// the last slice holding only one of them. Under the old latest-slice joint
+/// this was `UnknownKey`.
+#[test]
+fn competitors_last_seen_in_different_slices_are_comparable() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "c", 10, 4.0, 3.0),
+        // the final slice holds one duel that does not involve b at all
+        duel("a", "c", 20, 6.0, 1.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    // b last appeared at time 0; a and c at time 20. All three must resolve.
+    for (x, y) in [("a", "b"), ("b", "c"), ("a", "c")] {
+        let g = h
+            .posterior_of(&[(&x, 1.0), (&y, -1.0)])
+            .unwrap_or_else(|e| panic!("{x} - {y} should resolve across slices: {e}"));
+        assert!(g.sigma() > 0.0 && g.sigma().is_finite());
+    }
+}
+
+/// The mean must agree with what message passing reports, which is exact even
+/// with cycles. Only the second moment needs the joint.
+#[test]
+fn means_agree_with_the_marginals() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("b", "c", 5, 3.0, 1.0),
+        duel("a", "c", 10, 4.0, 2.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    for k in ["a", "b", "c"] {
+        let marginal = h.current_skill(&k).unwrap().mu();
+        let joint = h.posterior_of(&[(&k, 1.0)]).unwrap().mu();
+        assert!(
+            (marginal - joint).abs() < 1e-9,
+            "{k}: marginal {marginal}, joint {joint}"
+        );
+    }
+}
+
+/// With zero drift a competitor has one latent skill however many slices it
+/// appears in, so spreading the same events over time must not change the
+/// answer. This exercises the appearance-merging path.
+#[test]
+fn zero_drift_makes_slice_layout_irrelevant() {
+    let spread = {
+        let mut h = history(0.0);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "b", 10, 4.0, 3.0),
+            duel("a", "b", 20, 6.0, 1.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap()
+    };
+    let together = {
+        let mut h = history(0.0);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "b", 0, 4.0, 3.0),
+            duel("a", "b", 0, 6.0, 1.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+        h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap()
+    };
+
+    assert!(
+        (spread.sigma() - together.sigma()).abs() < 1e-9,
+        "zero drift: spread {} vs together {}",
+        spread.sigma(),
+        together.sigma()
+    );
+}
+
+/// More drift means less is carried forward from old evidence, so a comparison
+/// against a competitor last seen long ago must widen.
+#[test]
+fn drift_widens_a_comparison_across_time() {
+    let mut previous = 0.0;
+    for gamma in [0.0f64, 0.1, 0.5, 2.0] {
+        let mut h = history(gamma);
+        h.add_events(vec![
+            duel("a", "b", 0, 5.0, 2.0),
+            duel("a", "c", 100, 4.0, 3.0),
+        ])
+        .unwrap();
+        let _ = h.converge().unwrap();
+
+        // b was last seen at time 0; a at time 100.
+        let g = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+        assert!(
+            g.sigma() > previous,
+            "gamma={gamma}: sigma {} did not exceed {previous}",
+            g.sigma()
+        );
+        previous = g.sigma();
+    }
+}
+
+/// `posterior_of_at` pins the reading to a moment, where `posterior_of` takes
+/// each competitor wherever they were last seen.
+#[test]
+fn posterior_of_at_reads_as_of_a_time() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 10, 4.0, 3.0),
+        duel("a", "b", 20, 6.0, 1.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    let early = h.posterior_of_at(0, &[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    let late = h.posterior_of_at(20, &[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+    let latest = h.posterior_of(&[(&"a", 1.0), (&"b", -1.0)]).unwrap();
+
+    // Asking as of the final slice is the same as asking for the latest.
+    assert!((late.mu() - latest.mu()).abs() < 1e-9);
+    assert!((late.sigma() - latest.sigma()).abs() < 1e-9);
+
+    // Reading at time 0 is a different quantity, and the smoothed estimate
+    // there is informed by everything that came after.
+    assert!(
+        (early.mu() - late.mu()).abs() > 1e-6,
+        "as-of-0 and as-of-20 should differ: {} vs {}",
+        early.mu(),
+        late.mu()
+    );
+
+    // A time before any event has nothing to read.
+    assert!(h.posterior_of_at(-1, &[(&"a", 1.0)]).is_err());
+}
+
+/// Times between slices resolve to the latest appearance at or before them.
+#[test]
+fn a_time_between_slices_reads_the_previous_appearance() {
+    let mut h = history(GAMMA);
+    h.add_events(vec![
+        duel("a", "b", 0, 5.0, 2.0),
+        duel("a", "b", 100, 4.0, 3.0),
+    ])
+    .unwrap();
+    let _ = h.converge().unwrap();
+
+    let at_zero = h.posterior_of_at(0, &[(&"a", 1.0)]).unwrap();
+    let between = h.posterior_of_at(50, &[(&"a", 1.0)]).unwrap();
+    assert!((at_zero.mu() - between.mu()).abs() < 1e-12);
+    assert!((at_zero.sigma() - between.sigma()).abs() < 1e-12);
+}

tests/validation.rs

@@ -184,3 +184,74 @@ fn ingestion_rejects_weights_that_do_not_match_their_team() {
         "got {err:?}"
     );
 }
+
+/// `mu`, `sigma` and `beta` were the last unvalidated setters on
+/// `HistoryBuilder`, next to `p_draw`, `score_sigma` and `convergence`, which
+/// all assert eagerly.
+///
+/// Two of the rejected values are the quiet kind. A negative `sigma` or `beta`
+/// enters inference only as its square, so it produced bit-identical results
+/// to the positive value — the sign was dropped without comment.
+mod builder_parameters {
+    use trueskill_tt::History;
+
+    #[test]
+    #[should_panic(expected = "mu must be finite")]
+    fn a_non_finite_mu_is_rejected() {
+        let _ = History::builder().mu(f64::NAN);
+    }
+
+    #[test]
+    #[should_panic(expected = "sigma must be finite and positive")]
+    fn a_zero_sigma_is_rejected() {
+        let _ = History::builder().sigma(0.0);
+    }
+
+    #[test]
+    #[should_panic(expected = "sigma must be finite and positive")]
+    fn a_negative_sigma_is_rejected() {
+        let _ = History::builder().sigma(-8.33);
+    }
+
+    #[test]
+    #[should_panic(expected = "sigma must be finite and positive")]
+    fn an_infinite_sigma_is_rejected() {
+        let _ = History::builder().sigma(f64::INFINITY);
+    }
+
+    #[test]
+    #[should_panic(expected = "beta must be finite and non-negative")]
+    fn a_negative_beta_is_rejected() {
+        let _ = History::builder().beta(-4.17);
+    }
+
+    #[test]
+    #[should_panic(expected = "beta must be finite and non-negative")]
+    fn a_non_finite_beta_is_rejected() {
+        let _ = History::builder().beta(f64::NAN);
+    }
+
+    /// Zero beta is deliberately allowed: performance is then exactly skill.
+    /// It has to reach a different fit than a positive beta, or "allowed"
+    /// would just mean "not checked".
+    #[test]
+    fn a_zero_beta_is_allowed_and_changes_the_fit() {
+        let fit = |beta: f64| {
+            let mut h = History::builder()
+                .mu(25.0)
+                .sigma(25.0 / 3.0)
+                .beta(beta)
+                .build();
+            h.record_winner(&"a", &"b", 1).unwrap();
+            let _ = h.converge().unwrap();
+            h.current_skill(&"a").unwrap()
+        };
+        let zero = fit(0.0);
+        let positive = fit(25.0 / 6.0);
+        assert!(zero.pi().is_finite() && zero.pi() > 0.0);
+        assert!(
+            (zero.pi() - positive.pi()).abs() > 1e-6,
+            "zero beta must not merely be ignored: {zero:?} vs {positive:?}"
+        );
+    }
+}

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 { .. })
+    ));
+}