//! `TrueSkill` Through Time — Bayesian skill rating over a time axis. //! //! Where plain `TrueSkill` gives each competitor one running estimate, `TrueSkill` //! Through Time treats a whole history as a single model and infers skill *at //! every point in time*. Evidence flows both directions: a result today //! sharpens the estimate of who someone was last year, so early estimates stop //! being frozen guesses and comparisons across eras become meaningful. //! //! This is a Rust port of //! [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py). //! //! # Getting started //! //! Record results, converge, then read off skills: //! //! ``` //! use trueskill_tt::History; //! //! let mut history = History::default(); //! //! history.record_winner(&"alice", &"bob", 1)?; //! history.record_winner(&"bob", &"carol", 2)?; //! history.record_winner(&"alice", &"carol", 3)?; //! //! let report = history.converge()?; //! assert!(report.converged); //! //! let alice = history.current_skill("alice").unwrap(); //! assert!(alice.mu() > 0.0, "alice won every game she played"); //! # Ok::<(), trueskill_tt::InferenceError>(()) //! ``` //! //! Teams, weights, explicit rankings and continuous scores go through the //! fluent event builder: //! //! ``` //! use trueskill_tt::History; //! //! let mut history = History::builder().p_draw(0.1).build(); //! //! history //! .event(1) //! .team(["alice", "bob"]) //! .team(["carol", "dave"]) //! .ranking([0, 1]) //! .commit()?; //! //! history.converge()?; //! # Ok::<(), trueskill_tt::InferenceError>(()) //! ``` //! //! # Draws need a draw probability //! //! A `p_draw` of zero asserts that draws cannot happen, so a tied result has //! no representable likelihood and is rejected: //! //! ``` //! use trueskill_tt::{History, InferenceError}; //! //! let mut history = History::default(); // p_draw defaults to 0.0 //! let err = history.record_draw(&"alice", &"bob", 1).unwrap_err(); //! assert!(matches!(err, InferenceError::TieWithoutDrawProbability { .. })); //! ``` //! //! This also applies to [`Outcome::winner`] for three or more teams, which //! ties every loser. Configure a positive `p_draw` for those. //! //! # Core types //! //! - [`History`] — the top-level container: ingests events, runs //! forward/backward message passing, and answers queries. //! - [`Gaussian`] — the probability type, stored in natural parameters //! (`pi = 1/sigma²`, `tau = mu/sigma²`) so message passing is add/subtract. //! - [`Game`] — one match in isolation, for scoring a hypothetical without a //! history. //! - [`Outcome`] — how a match ended: ranks, or continuous scores. //! - [`Rating`] — a competitor's static configuration (prior, `beta`, drift). //! //! # Feature flags //! //! - `approx` — implements [`approx`](https://docs.rs/approx) equality traits //! for [`Gaussian`]. Useful in tests. //! - `rayon` — parallelises the within-slice sweep and the per-slice passes of //! `learning_curves`/`log_evidence`. Opt-in; results stay bit-identical //! regardless of worker count. #![forbid(unsafe_code)] /// Compiles every `rust` block in `README.md` as a doctest. /// /// The README is not the crate's front page — the module docs above are — so it /// is pulled in here rather than via a crate-level `#![doc = ...]`, purely so /// its examples are type-checked. Without this nothing compiled them, and they /// had drifted far enough that four blocks no longer built (#35). `cfg(doctest)` /// means this type exists only while collecting doctests. /// /// Blocks that are illustrative rather than runnable are fenced as `text`. #[cfg(doctest)] #[doc = include_str!("../README.md")] pub struct ReadmeDoctests; use std::{ cmp::Reverse, f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2}, }; #[cfg(feature = "approx")] mod approx; pub(crate) mod arena; mod time; mod time_slice; pub use time_slice::{EventKind, TimeSlice}; mod acquisition; mod color_group; mod competitor; mod convergence; pub mod drift; mod error; mod event; mod event_builder; pub(crate) mod factor; mod game; pub mod gaussian; mod history; mod joint; mod key_table; mod matrix; mod observer; mod outcome; mod predict; pub(crate) mod quadrature; mod rating; pub 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, 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, 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 time::{Time, Untimed}; pub const BETA: f64 = 1.0; pub const MU: f64 = 0.0; pub const SIGMA: f64 = BETA * 6.0; pub const GAMMA: f64 = BETA * 0.03; pub const P_DRAW: f64 = 0.0; pub const EPSILON: f64 = 1e-6; /// Default cap on convergence sweeps. /// /// **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. /// /// The outcome space holds `n! * 2^(n-1)` events, so it grows factorially: /// 1_920 at five teams, 23_040 at six, 322_560 at seven. Six is where /// enumerating on a caller's behalf stops being reasonable. pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION; const SQRT_TAU: f64 = 2.5066282746310002; /// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction. const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3; /// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form. const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4; /// How many window widths into the tail before a tie window is treated as a /// half-line. Beyond this the truncated mass is concentrated within `1/alpha` /// of the near edge, so the far edge contributes nothing measurable. const HALF_LINE_WINDOW: f64 = 10.0; /// Where `v - alpha` switches from subtraction to its asymptotic series. /// /// The subtraction loses roughly `eps * alpha^2` of relative precision, and the /// four-term series is good to ~1e-10 by here, so the two are at their closest /// agreement around this point. Below it the subtraction is exact; above it the /// series is. const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0; pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0); pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0); pub const N_INF: Gaussian = Gaussian::from_ms(0.0, f64::INFINITY); #[derive(Copy, Clone, Default, PartialEq, PartialOrd, Eq, Ord, Hash, Debug)] pub struct Index(usize); impl Index { /// The underlying slot number. /// /// Indices are dense and assigned in interning order, so this is usable as /// a key into a caller-side side table. #[must_use] pub fn get(self) -> usize { self.0 } } impl From for Index { fn from(ix: usize) -> Self { Self(ix) } } impl From for usize { fn from(idx: Index) -> Self { idx.0 } } /// 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 { 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; } debug_assert!(y >= 0.0, "y must be nonnegative"); if y == 0.0 { return f64::INFINITY; } if y >= 1.0 { y = 2.0 - y; } 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 // reflection of the root. Newton then has to cross the origin to get back, // which a fixed iteration count does not manage: measured against the true // value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error // grew as y shrank until `compute_margin` stopped being monotone in // `p_draw` altogether. let mut x = -FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t); for _ in 0..3 { let err = erfc(x) - y; x += err / (FRAC_2_SQRT_PI * libm::exp(-(x * x)) - x * err) } if y < 1.0 { x } else { -x } } fn ppf(p: f64, mu: f64, sigma: f64) -> f64 { mu - sigma * SQRT_2 * erfc_inv(2.0 * p) } fn compute_margin(p_draw: f64, sd: f64) -> f64 { ppf(0.5 - p_draw / 2.0, 0.0, sd).abs() } pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 { let z = -(x - mu) / (sigma * SQRT_2); 0.5 * erfc(z) } /// `P(X > x)` for `X ~ N(mu, sigma^2)`. /// /// The survival function, computed directly rather than as `1 - cdf(..)`. /// /// The two are algebraically identical and numerically are not. `cdf` returns /// a value approaching 1 for an upper tail, so subtracting it from 1 cancels /// away every significant digit the tail had: measured against this function, /// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly /// zero beyond about 8.3 sigma — where the true value is still 1e-19 and /// perfectly representable. `erfc` holds *relative* accuracy all the way down /// to 1e-296, so the precision is there to keep; only the subtraction threw it /// away. /// /// This matters most where evidence is smallest, which is exactly where an /// upset makes it interesting: `ln` of a clamped zero is -708 regardless of /// whether the truth was -43 or -600. pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 { 0.5 * erfc((x - mu) / (sigma * SQRT_2)) } /// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`. /// /// Exists so the exponential factor common to a Gaussian density and its tail /// integral can be cancelled *analytically* instead of being computed twice /// and divided. Both underflow to zero past about 26 sigma, and their ratio is /// then `0/0` — finite in the limit, `NaN` in floating point. fn erfcx(x: f64) -> f64 { if x < 2.0 { // Below the crossover neither factor is extreme: erfc is O(1) and // exp(x^2) is at most e^4, so the direct product is exact enough and // cheaper than the continued fraction. libm::exp(x * x) * erfc(x) } else { // erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))), // evaluated by backward recurrence. Converges quickly for x >= 2 and, // unlike the product form, never touches an exponential. let mut f = 0.0; for n in (1..=60u32).rev() { f = (f64::from(n) * 0.5) / (x + f); } FRAC_1_SQRT_PI / (x + f) } } /// `ln` of the normal density at `x`. /// /// The density itself underflows to zero past about 38 sigma, and `ln` of a /// clamped zero is -708 whatever the truth was. The log form is a polynomial: /// it stays exact at any separation, and the values it produces (-5001 nats at /// 100 sigma, -500001 at 1000) are perfectly representable. pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 { let z = (x - mu) / sigma; -libm::log(SQRT_TAU * sigma) - 0.5 * z * z } /// `ln P(X > x)` for `X ~ N(mu, sigma^2)`. /// /// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is /// factored out analytically via `erfcx`, so this never underflows — where /// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero. pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 { let z = (x - mu) / sigma; if z > 0.0 { // ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y). -std::f64::consts::LN_2 - 0.5 * z * z + libm::log(erfcx(z / SQRT_2)) } else { // The mass here is at least a half; nothing to lose. libm::log(sf(x, mu, sigma)) } } /// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`. /// /// When the interval sits in a tail both endpoint probabilities underflow /// together, so their difference is taken in scaled form with the shared /// exponential factored out. When it straddles the mean nothing is small and /// the direct difference is exact. pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 { let z_lo = (lo - mu) / sigma; let z_hi = (hi - mu) / sigma; if z_hi <= z_lo { return f64::NEG_INFINITY; } // Fold a lower-tail interval onto the upper tail; the normal is symmetric. let (near, far) = if z_lo >= 0.0 { (z_lo, z_hi) } else if z_hi <= 0.0 { (-z_hi, -z_lo) } else { // Straddles the mean: the interval holds a non-negligible share of the // mass, so neither endpoint is near enough to 1 to cancel. return libm::log((cdf(hi, mu, sigma) - cdf(lo, mu, sigma)).max(f64::MIN_POSITIVE)); }; let (a, b) = (near / SQRT_2, far / SQRT_2); // b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow. let scale = libm::exp(a * a - b * b); let bracket = erfcx(a) - scale * erfcx(b); if bracket <= 0.0 { return f64::NEG_INFINITY; } -std::f64::consts::LN_2 - a * a + libm::log(bracket) } fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { let normalizer = (SQRT_TAU * sigma).powi(-1); let functional = libm::exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma)); normalizer * functional } /// Truncated-Gaussian correction terms `(v, w)`. /// /// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose /// numerator and denominator underflow together in the tails, so both are /// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled /// analytically rather than evaluated and divided out. Without that, a /// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN` /// straight into the posterior. /// Truncation terms for a boundary `alpha` standard deviations into the upper /// tail, from the asymptotic expansion of the inverse Mills ratio. /// /// `v` tends to `alpha` out here, so the gap between them cannot be obtained by /// subtracting one from the other — the series computes the gap directly, and /// `w = v * gap` then never forms the difference of two large near-equal /// numbers. A far-tail *window* behaves like a half-line once it is more than a /// few multiples of its own width from the mean, so the tie branch shares this. fn half_line_truncation(alpha: f64) -> (f64, f64) { let inv = alpha.recip(); let inv_sq = inv * inv; let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq))); let v = alpha + gap; (v, v * gap) } fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { if !tie { let alpha = (margin - mu) / sigma; // v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs // the gap `v - alpha` as well as v itself. Far into the tail v tends to // alpha, so that gap is a subtraction of two nearly equal numbers and // loses every digit it has: at alpha = 1e6 it drove w above 1 and made // `sqrt(1 - w)` NaN. Past the crossover the gap comes from its // asymptotic series instead, which has no subtraction in it. if alpha >= ASYMPTOTIC_MILLS_ALPHA { return half_line_truncation(alpha); } let (v, gap) = if alpha > 0.0 { // Both terms carry exp(-alpha^2 / 2); in scaled form it cancels // and the result stays exact however far into the tail alpha sits. let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2); (v, v - alpha) } else { // Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing. let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0); (v, v - alpha) }; (v, v * gap) } else { // v is odd in mu and w is even, so fold to mu <= 0. Both truncation // points then sit in the upper tail, where the scaled form applies. let flipped = mu > 0.0; let mu = if flipped { -mu } else { mu }; let alpha = (-margin - mu) / sigma; let beta = (margin - mu) / sigma; // `w` comes out of `v * v - u`, and both terms grow as alpha^2 while // their difference stays O(1) — at alpha = 1e9 that subtraction had no // digits left and returned w = -128, making `sqrt(1 - w)` nonsense. // Once the window sits many of its own widths into the tail it is // indistinguishable from a half-line, so the asymptotic covers it with // no subtraction at all. if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW { let (v, w) = half_line_truncation(alpha); return (if flipped { -v } else { v }, w); } let (v, u) = if alpha > 0.0 { // beta > alpha > 0, so this ratio of exponentials is at most 1 and // cannot overflow. let scale = libm::exp(0.5 * (alpha * alpha - beta * beta)); let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2)); ( (1.0 - scale) / SQRT_TAU / denominator, (alpha - beta * scale) / SQRT_TAU / denominator, ) } else { // The interval straddles the mean, so nothing here is small. let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0); ( (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator, (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator, ) }; let w = -(u - v.powi(2)); (if flipped { -v } else { v }, w) } } fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { let (v, w) = v_w(mu, sigma, margin, tie); let mu_trunc = mu + sigma * v; let sigma_trunc = sigma * (1.0 - w).sqrt(); (mu_trunc, sigma_trunc) } pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian { let (mu, sigma) = trunc(n.mu(), n.sigma(), margin, tie); Gaussian::from_ms(mu, sigma) } /// Componentwise maximum that **propagates** NaN rather than dropping it. /// /// Every caller folds this as `tuple_max(accumulator, new)`. A plain `>` /// comparison is false against NaN, so a NaN accumulator would be replaced by /// the next finite delta and the breakdown would vanish — leaving `step_is_finite` /// to pass on a fit that is already NaN. Because the fold runs over a `HashMap`, /// whether that happened depended on per-process hash order: measured, a NaN fit /// was reported as `converged: true` in 16 of 30 runs on identical input. /// /// `f64::max` is not a substitute: it also ignores NaN by design, which is the /// same defect wearing a standard-library name. pub(crate) fn tuple_max(v1: (f64, f64), v2: (f64, f64)) -> (f64, f64) { ( max_propagating_nan(v1.0, v2.0), max_propagating_nan(v1.1, v2.1), ) } fn max_propagating_nan(a: f64, b: f64) -> f64 { if a.is_nan() || b.is_nan() { f64::NAN } else if a > b { a } else { b } } pub(crate) fn tuple_gt(t: (f64, f64), e: f64) -> bool { t.0 > e || t.1 > e } /// Whether a convergence step is finite in both components. /// /// A NaN step means EP broke down numerically. Because every comparison /// against NaN is false, `tuple_gt` reads NaN as "below epsilon" — so /// convergence checks must test finiteness explicitly rather than inferring /// success from `!tuple_gt(..)`. pub(crate) fn step_is_finite(t: (f64, f64)) -> bool { t.0.is_finite() && t.1.is_finite() } /// Whether a step counts as converged: finite *and* within `epsilon`. pub(crate) fn step_converged(t: (f64, f64), epsilon: f64) -> bool { step_is_finite(t) && !tuple_gt(t, epsilon) } /// Indices of the first pair of teams sharing a rank, if any. /// /// A tie is only representable when the draw probability is positive: with /// `p_draw == 0.0` the truncation margin collapses to zero and the two-sided /// tie update evaluates `0/0`. Callers use this to reject such events before /// they reach inference. pub(crate) fn first_tied_pair(ranks: &[u32]) -> Option<(usize, usize)> { for (i, a) in ranks.iter().enumerate() { for (j, b) in ranks.iter().enumerate().skip(i + 1) { if a == b { return Some((i, j)); } } } None } /// As `first_tied_pair`, but over the engine's internal `f64` outputs. /// /// Ranks reach the engine already converted to descending `f64` outputs, and /// `Game` decides a tie by exact equality of those values — so this mirrors /// the comparison inference itself performs. pub(crate) fn first_tied_output(outputs: &[f64]) -> Option<(usize, usize)> { for (i, a) in outputs.iter().enumerate() { for (j, b) in outputs.iter().enumerate().skip(i + 1) { if a == b { return Some((i, j)); } } } None } pub(crate) fn sort_time(xs: &[T], reverse: bool) -> Vec { let mut x: Vec<(usize, T)> = xs.iter().enumerate().map(|(i, &t)| (i, t)).collect(); if reverse { x.sort_by_key(|&(_, t)| Reverse(t)); } else { x.sort_by_key(|&(_, t)| t); } x.into_iter().map(|(i, _)| i).collect() } /// Calculates the match quality of the given rating groups. A result is the draw probability in the association /// /// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a /// perfectly balanced match. /// /// # Panics /// /// Panics if fewer than two rating groups are supplied, or if any group is /// empty — match quality is a property of a contest between at least two /// non-empty sides. /// /// Also panics with "cannot invert a singular matrix" when every rating has /// zero sigma *and* `beta` is zero. Nothing is then uncertain, so there is no /// distribution to take the quality of; `Gaussian::from_ms(mu, 0.0)` is a point /// mass and its `mu()` is not even well defined. Documented rather than /// converted, because the input has no meaningful answer rather than an /// awkward one. #[must_use] pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { assert!( rating_groups.len() >= 2, "quality() requires at least 2 rating groups, got {}", rating_groups.len() ); assert!( rating_groups.iter().all(|group| !group.is_empty()), "quality() requires every rating group to be non-empty" ); let flatten_ratings = rating_groups .iter() .flat_map(|group| group.iter()) .collect::>(); let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice(); let length = flatten_ratings.len(); let mut mean_matrix = Matrix::new(length, 1); for (i, rating) in flatten_ratings.iter().enumerate() { mean_matrix[(i, 0)] = rating.mu(); } let mut variance_matrix = Matrix::new(length, length); for (i, rating) in flatten_ratings.iter().enumerate() { variance_matrix[(i, i)] = rating.sigma().powi(2); } let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length); // Row `row` contrasts group `row` (+weight) against group `row + 1` // (-weight). `t` is the column where the current group's players start; // the negative block begins immediately after it. let mut t = 0; for (row, group) in rating_groups.windows(2).enumerate() { let current = group[0]; let next = group[1]; for n in t..t + current.len() { rotated_a_matrix[(row, n)] = flatten_weights[n]; } t += current.len(); for n in t..t + next.len() { rotated_a_matrix[(row, n)] = -flatten_weights[n]; } } let a_matrix = rotated_a_matrix.transpose(); let ata = beta.powi(2) * &rotated_a_matrix * &a_matrix; let atsa = &rotated_a_matrix * &variance_matrix * &a_matrix; let start = mean_matrix.transpose() * &a_matrix; let middle = &ata + &atsa; let end = &rotated_a_matrix * &mean_matrix; let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant(); // `sqrt(det(ata) / det(middle))`, taken in log space. Both determinants are // products of `k - 1` diagonal entries, so they leave `f64`'s range long // before their ratio does: measured at the crate defaults, 150 groups was // correct at `8.45e-53`, 200 returned `0`, and 250 returned `NaN` where the // true value is `9.51e-88`. With a small beta it is sharper still — at // `sigma = beta = 1e-3`, 60 groups returned `NaN` against a true `1.32e-9`. // // The ratio is what the answer needs and it is representable throughout, so // the intermediates are the only thing that ever overflowed. let ln_s_arg = ata.ln_abs_determinant() - middle.ln_abs_determinant(); libm::exp(e_arg + 0.5 * ln_s_arg) } #[cfg(test)] mod tests { /// A NaN must survive the fold from ANY position, not only the last. /// /// The fold runs over a `HashMap`, so "last" is per-process hash order. The /// end-to-end symptom was a NaN fit reported as `converged: true` in 16 of /// 30 runs on identical input; these three cases are the deterministic form /// of that, so a regression cannot hide behind a lucky seed. #[test] fn tuple_max_propagates_a_nan_from_any_position() { let nan = (f64::NAN, f64::NAN); let small = (1e-9, 1e-9); let big = (1e-3, 1e-3); // NaN last. let step = tuple_max(tuple_max(big, small), nan); assert!(!step_is_finite(step), "NaN last: {step:?}"); // NaN middle. let step = tuple_max(tuple_max(big, nan), small); assert!(!step_is_finite(step), "NaN middle: {step:?}"); // NaN first — the case a plain `>` comparison drops. let step = tuple_max(tuple_max(nan, big), small); assert!(!step_is_finite(step), "NaN first: {step:?}"); } /// `f64::max` would pass the test above's first two cases and fail the /// third, so pin that it is not what we use. #[test] fn tuple_max_is_not_f64_max() { assert!( f64::max(f64::NAN, 1.0) == 1.0, "premise: f64::max drops NaN" ); let (a, _) = tuple_max((f64::NAN, 0.0), (1.0, 0.0)); assert!(a.is_nan(), "tuple_max must not drop what f64::max drops"); } /// Ordinary values are unaffected. #[test] fn tuple_max_still_takes_the_larger_component() { assert_eq!(tuple_max((1.0, 5.0), (3.0, 2.0)), (3.0, 5.0)); assert_eq!(tuple_max((3.0, 2.0), (1.0, 5.0)), (3.0, 5.0)); } use ::approx::assert_ulps_eq; use super::*; #[test] fn test_sort_time() { assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]); } /// Upper-tail values of the standard normal, from published tables. The /// point is not the digits — these are 7-digit table values — but that a /// number comes back at all: `1 - cdf` returned exactly zero for every one /// of these. #[test] fn survival_function_survives_the_far_tail() { for (z, expected) in [ (9.0f64, 1.128_588e-19), (12.0, 1.776_482e-33), (20.0, 2.753_624e-89), (37.0, 5.725_571e-300), ] { let got = sf(z, 0.0, 1.0); assert!(got > 0.0, "sf({z}) collapsed to zero"); assert!( (got - expected).abs() / expected < 1e-6, // published table values, 7 digits "sf({z}) = {got}, expected ~{expected}" ); assert_eq!( 1.0 - cdf(z, 0.0, 1.0), 0.0, "the naive form should still be zero here" ); } } /// Where no cancellation happens the two forms must agree exactly enough /// that nothing else in the crate shifts. #[test] fn survival_function_matches_the_naive_form_where_that_form_works() { for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] { let naive = 1.0 - cdf(z, 0.0, 1.0); let direct = sf(z, 0.0, 1.0); assert!( (naive - direct).abs() < 1e-15, "z={z}: naive {naive} vs direct {direct}" ); } } #[test] fn survival_and_cdf_partition_the_mass() { for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] { let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0); assert!((total - 1.0).abs() < 1e-15, "z={z}: {total}"); } } /// `erfcx` switches formulation at x = 2; the two sides must meet. #[test] fn erfcx_is_continuous_across_its_crossover() { for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] { let direct = (x * x).exp() * erfc(x); let scaled = erfcx(x); assert!( (direct - scaled).abs() / scaled < 1e-14, "x={x}: direct {direct} vs erfcx {scaled}" ); } } /// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly /// where `exp(x^2)` overflows and `erfc(x)` underflows. #[test] fn erfcx_stays_finite_where_its_factors_do_not() { for x in [27.0f64, 50.0, 1.0e3, 1.0e8] { let scaled = erfcx(x); assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}"); // Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)). let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt()); assert!( (scaled - asymptote).abs() / asymptote < 1e-2, "erfcx({x}) = {scaled} strays from its asymptote {asymptote}" ); assert!( (x * x).exp().is_infinite(), "x={x} should overflow the direct form" ); } } /// Truncation must never produce a non-finite posterior. Before the scaled /// formulation these returned NaN from `0 / 0` past about 39 sigma. #[test] fn truncation_stays_finite_arbitrarily_far_into_the_tail() { for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] { for tie in [false, true] { let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie); assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}"); assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}"); // sigma_trunc = sigma * sqrt(1 - w) must stay real. assert!( (0.0..=1.0).contains(&w), "alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary" ); let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie); assert!( mu_t.is_finite() && sigma_t.is_finite(), "alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})" ); } } } /// The Mills gap switches from subtraction to series at alpha = 100. Both /// are supposed to be right there; if they disagree, the crossover is in /// the wrong place. #[test] fn the_mills_gap_series_meets_the_scaled_form() { for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] { let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha; let inv = alpha.recip(); let inv_sq = inv * inv; let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq))); assert!( (scaled - series).abs() / series < 1e-9, "alpha={alpha}: scaled {scaled} vs series {series}" ); } } /// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and /// w is even. Assert the symmetry the implementation relies on. #[test] fn tie_truncation_is_odd_in_v_and_even_in_w() { for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] { let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true); let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true); assert!( (v_pos + v_neg).abs() < 1e-9, "mu={mu}: v should be odd, got {v_pos} and {v_neg}" ); assert!( (w_pos - w_neg).abs() < 1e-9, "mu={mu}: w should be even, got {w_pos} and {w_neg}" ); } } /// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the /// mirror image of the root. Three fixed iterations could not cross back, /// so the error grew as the argument shrank: at `p_draw = 0.99` the margin /// came out 0.503 where the answer is 2.576. #[test] fn erfc_inv_matches_known_quantiles() { // sqrt(2) * erfc_inv(1 - p) is the standard normal quantile // Phi^-1((1 + p) / 2). for (p, exact) in [ (0.5f64, 0.674_489_750_196_081_7f64), (0.9, 1.644_853_626_951_472_7), (0.95, 1.959_963_984_540_054_2), (0.99, 2.575_829_303_548_9), (0.999, 3.290_526_731_491_896_4), ] { let got = SQRT_2 * erfc_inv(1.0 - p); assert!( (got - exact).abs() / exact < 1e-14, "p={p}: got {got}, exact {exact}" ); } } /// The draw margin must grow with the draw probability. It did not: it ran /// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 -> /// 0.999, which is not a rounding error but a broken function. /// Deep in the tail the accuracy limit is the *caller's* argument, not this /// function. /// /// `compute_margin(0.999999, ..)` computes `1.0 - p_draw`, and 0.999999 is /// not representable: the subtraction cancels and leaves 2.9e-11 of /// relative error in the argument before `erfc_inv` is even entered. Given /// an exactly-representable argument the result is good to 1.8e-16, so this /// is inherent to taking `p_draw` near one rather than something to fix /// here. At `p_draw = 0.999` the whole path is still accurate to 4e-16. /// /// Worth pinning: measured against a 70-digit reference, `puruspe`'s /// `inverfc` returns the identical wrong value for the identical reason, /// which is what makes it clear the fault is upstream of both. #[test] fn erfc_inv_is_exact_given_an_exactly_representable_argument() { // erfc(z / sqrt2) = 1e-6 exactly, so z = Phi^-1(0.9999995). let got = SQRT_2 * erfc_inv(1e-6); let exact = 4.891_638_475_698_59; assert!( (got - exact).abs() / exact < 1e-14, "got {got}, exact {exact}" ); } #[test] fn compute_margin_is_monotone_in_the_draw_probability() { let mut previous = 0.0; for p_draw in [ 0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999, ] { let margin = compute_margin(p_draw, 1.0); assert!( margin > previous, "p_draw={p_draw}: margin {margin} did not exceed {previous}" ); previous = margin; } } /// Round-tripping the margin back through the model's own CDF must recover /// the draw probability it was built from. #[test] fn compute_margin_round_trips_through_the_cdf() { for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] { for sd in [0.5f64, 1.0, 5.892_557] { let margin = compute_margin(p_draw, sd); // P(|X| < margin) for X ~ N(0, sd^2). let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd); assert!( (recovered - p_draw).abs() < 1e-14, "p_draw={p_draw} sd={sd}: recovered {recovered}" ); } } } /// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where /// the linear forms underflow. Past ~38 sigma the linear value is zero and /// its log is whatever floor it was clamped to. #[test] fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() { for z in [40.0f64, 60.0, 100.0, 1000.0] { assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}"); assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}"); let lp = ln_pdf(z, 0.0, 1.0); let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z; assert!( (lp - expected_lp).abs() < 1e-9, "ln_pdf({z}) = {lp}, expected {expected_lp}" ); let ls = ln_sf(z, 0.0, 1.0); // ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z. let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln(); assert!( (ls - approx).abs() / approx.abs() < 1e-3, "ln_sf({z}) = {ls}, asymptote {approx}" ); assert!( ls < f64::MIN_POSITIVE.ln(), "ln_sf({z}) still on the clamp floor" ); } } /// Where nothing underflows, the log helpers must agree with the direct /// forms exactly enough that nothing else in the crate shifts. #[test] fn log_space_helpers_agree_with_the_linear_forms_in_range() { for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] { let lp = ln_pdf(z, 0.5, 2.0); let direct_pdf = pdf(z, 0.5, 2.0); assert!( (lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf, "ln_pdf at {z}: {} vs {direct_pdf}", lp.exp() ); let ls = ln_sf(z, 0.5, 2.0); let direct = sf(z, 0.5, 2.0); assert!( (ls.exp() - direct).abs() <= 1e-13 * direct.max(1e-300), "ln_sf at {z}: {} vs {direct}", ls.exp() ); } } #[test] fn ln_interval_matches_the_direct_difference_when_nothing_is_small() { for mu in [-2.0f64, 0.0, 0.5, 2.0] { let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0); let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp(); assert!( (logged - direct).abs() <= 1e-13 * direct, "mu={mu}: {logged} vs {direct}" ); } } /// A window far out in the tail: both endpoints underflow together, so the /// difference has to be taken in scaled form. #[test] fn ln_interval_survives_a_window_deep_in_the_tail() { for mu in [-50.0f64, -100.0, -1000.0] { let logged = ln_interval(-1.0, 1.0, mu, 1.0); assert!(logged.is_finite(), "mu={mu}: {logged}"); assert!( logged < f64::MIN_POSITIVE.ln(), "mu={mu}: {logged} is stuck on the clamp floor" ); // Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)). let near = ln_sf(-1.0, mu, 1.0); assert!( (logged - near).abs() < 5.0, "mu={mu}: {logged} strays from the near-edge tail {near}" ); } } #[test] fn test_quality() { let a = Gaussian::from_ms(25.0, 3.0); let b = Gaussian::from_ms(25.0, 3.0); let q = quality(&[&[a], &[b]], 25.0 / 3.0 / 2.0); assert_ulps_eq!(q, 0.8115343414514944, epsilon = 1e-6) } }