diff --git a/src/acquisition.rs b/src/acquisition.rs index 29ad6a2..76a236f 100644 --- a/src/acquisition.rs +++ b/src/acquisition.rs @@ -47,7 +47,38 @@ fn kl_divergence(q: Gaussian, p: Gaussian) -> f64 { } let mean_gap = q.mu() - p.mu(); - 0.5 * (libm::log(var_p / var_q) + (var_q + mean_gap * mean_gap) / var_p - 1.0) + + // Algebraically `0.5 * (ln(var_p/var_q) + (var_q + gap^2)/var_p - 1)`, but + // written so that neither term can go negative. + // + // The direct form cancels against its `- 1.0` for two near-identical + // distributions and returns a *negative* divergence — measured, 762 082 of + // 3 000 000 near-identical pairs, worst `-5.55e-17`, which is exactly one + // ULP of the 1.0. It also loses the answer entirely where it is small: + // at `var_q/var_p - 1 = 1e-9` the direct form gives `0.0` where the true + // value is `2.5e-19`. + // + // With `u = var_q/var_p - 1` the variance part is `0.5 * (u - ln(1+u))`, + // which is non-negative for every `u > -1`, and the mean part is a square + // over a positive variance. Non-negativity is then structural rather than + // incidental. + let u = var_q / var_p - 1.0; + 0.5 * u_minus_ln1p(u) + mean_gap * mean_gap / (2.0 * var_p) +} + +/// `u - ln(1 + u)`, without the cancellation that spelling invites. +/// +/// Both terms are approximately `u` for small `u`, so the subtraction loses +/// everything just where the result matters. The Taylor series +/// `u^2/2 - u^3/3 + u^4/4 - ...` is exact in that regime and manifestly +/// non-negative, since `u^2/2` dominates. +fn u_minus_ln1p(u: f64) -> f64 { + if u.abs() < 1e-4 { + let u2 = u * u; + u2 * (0.5 - u / 3.0 + u2 / 4.0) + } else { + u - libm::log1p(u) + } } /// Expected information gain of a hypothetical matchup, in nats. @@ -146,7 +177,7 @@ pub fn expected_information_gain>( let mut gain = 0.0; - for (ranks, probability) in predict::outcome_distribution(&performances, &margins) { + for (ranks, probability) in predict::outcome_distribution(&performances, &margins)? { if probability <= NEGLIGIBLE { continue; } diff --git a/src/drift.rs b/src/drift.rs index c751624..8f89415 100644 --- a/src/drift.rs +++ b/src/drift.rs @@ -21,6 +21,23 @@ pub trait Drift: Copy + Debug + Send + Sync { /// /// For `Time = i64`: variance added is `(to - from) * gamma^2`. /// For `Time = Untimed`: elapsed is always 0, so drift is always 0. +/// +/// # The sign of `gamma` is not meaningful +/// +/// `gamma` enters only as `gamma * gamma`, so `ConstantDrift(-0.05)` produces +/// results **bit identical** to `ConstantDrift(0.05)`. That is the same +/// sign-absorption `HistoryBuilder::sigma`, `HistoryBuilder::beta`, +/// `Gaussian::from_ms` and `Rating::new` all reject outright. +/// +/// It is not rejected here because the field is public and positional, so +/// there is no constructor to intercept — sealing it would break every +/// `ConstantDrift(x)` in existence for a case whose *resulting model* is +/// perfectly valid, just not the one a caller writing a minus sign expected. +/// +/// A non-finite `gamma` is a different matter and **is** rejected: +/// `History::converge` validates the drift variance each competitor actually +/// accumulates, which also covers a custom [`Drift`] implementation, and +/// reports `InferenceError::InvalidParameter`. #[derive(Clone, Copy, Debug)] pub struct ConstantDrift(pub f64); diff --git a/src/error.rs b/src/error.rs index de9d905..8f24478 100644 --- a/src/error.rs +++ b/src/error.rs @@ -131,6 +131,28 @@ pub enum InferenceError { AlreadyRegistered { key: String }, /// A prediction was given a team with no members. EmptyTeam { team: usize }, + /// The prediction grid cannot resolve the narrowest feature in the matchup. + /// + /// `predict_outcome` and `predict_ranking` integrate every team's density + /// on one shared grid, whose resolution is set by the narrowest sigma (or a + /// narrower draw margin). When the widest and narrowest are far enough + /// apart, resolving the narrow one across the wide one's support needs more + /// nodes than the grid is allowed to hold. + /// + /// Reported rather than clamped. Clamping is what this replaced, and it + /// returned probabilities greater than one — measured, a `P` of 2.79 and a + /// `Prediction::total()` of 5.41 — because the trapezoid rule stops + /// resolving a density once the step exceeds roughly 1.7 of its sigma. + /// + /// `predict_win_probabilities` answers the same matchup through adaptive + /// quadrature and is accurate here; use it when only the per-team win + /// probabilities are needed. + GridTooCoarse { + /// Nodes required to resolve the narrowest feature. + needed: usize, + /// Nodes the grid may hold. + max: 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. @@ -219,6 +241,15 @@ impl fmt::Display for InferenceError { Self::EmptyTeam { team } => { write!(f, "team {team} has no members") } + Self::GridTooCoarse { needed, max } => { + write!( + f, + "the prediction grid needs {needed} nodes to resolve the narrowest \ + team's density across the widest team's support, but may hold only \ + {max}; the sigmas in this matchup are too far apart to integrate on \ + one grid. Use predict_win_probabilities, which is accurate here" + ) + } Self::JointUnavailable { reason } => { write!(f, "no exact joint posterior is available: {reason}") } diff --git a/src/factor/margin.rs b/src/factor/margin.rs index 61b76d9..dfcb0db 100644 --- a/src/factor/margin.rs +++ b/src/factor/margin.rs @@ -81,7 +81,7 @@ fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 { // `hypot`, not `sqrt(a^2 + b^2)`: squaring overflows to infinity above a // sigma of ~1.3e154 and flushes to zero below ~1.5e-154, and `Gaussian`'s // constructors are public so a caller can reach both. - let combined_sigma = cavity.sigma().hypot(sigma); + let combined_sigma = libm::hypot(cavity.sigma(), sigma); let value = ln_pdf(m_obs, cavity.mu(), combined_sigma); // A degenerate cavity (infinite sigma) is the only way to reach a @@ -89,7 +89,7 @@ fn cavity_log_evidence(cavity: Gaussian, m_obs: f64, sigma: f64) -> f64 { if value.is_finite() { value } else { - f64::MIN_POSITIVE.ln() + libm::log(f64::MIN_POSITIVE) } } diff --git a/src/factor/trunc.rs b/src/factor/trunc.rs index 3105dc6..3ffcea3 100644 --- a/src/factor/trunc.rs +++ b/src/factor/trunc.rs @@ -95,7 +95,7 @@ fn cavity_log_evidence(diff: Gaussian, margin: f64, tie: bool) -> f64 { if value.is_finite() { value } else { - f64::MIN_POSITIVE.ln() + libm::log(f64::MIN_POSITIVE) } } @@ -203,7 +203,7 @@ mod tests { let approx = -0.5 * z * z - z.ln() - (2.0 * std::f64::consts::PI).sqrt().ln(); assert!( - got < f64::MIN_POSITIVE.ln(), + got < libm::log(f64::MIN_POSITIVE), "mu={mu}: {got} is still stuck on the old clamp floor" ); assert!( diff --git a/src/gaussian.rs b/src/gaussian.rs index 23cc261..08ba486 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -18,8 +18,44 @@ pub struct Gaussian { impl Gaussian { /// Construct from mean and standard deviation. + /// + /// # Panics + /// + /// Panics if `sigma` is negative. NaN is deliberately allowed through: a + /// broken fit produces one, and `converge` reports that as + /// `NonFiniteResult` rather than panicking mid-inference. + /// + /// A negative sigma used to be accepted and returned results **bit + /// identical** to its absolute value, because sigma only ever enters as + /// `sigma * sigma`. The sign was not rejected and not honoured; it simply + /// vanished. That is the same defect `HistoryBuilder::sigma`, + /// `HistoryBuilder::beta` and `Member::with_drift_scale` already reject. + /// + /// # Very small sigma + /// + /// `pi = 1 / sigma^2` leaves `f64`'s range below about `1.5e-154`, and + /// `tau = mu * pi` overflows sooner still — at a threshold that depends on + /// `mu`, so there is a band where `pi` is finite and only `tau` is not. + /// Both land on the same point-mass representation the `sigma == 0.0` + /// branch produces, and a point mass with a non-zero mean has `mu() = NaN`, + /// because `tau / pi` is `inf / inf`. + /// + /// This is not rejected, because `approx` legitimately produces a very + /// small truncated sigma and inference must not panic. It is worth knowing + /// that such a `Gaussian` is not equal to itself, so two identical + /// declarations of one can be reported as conflicting. #[must_use] pub const fn from_ms(mu: f64, sigma: f64) -> Self { + // NaN is admitted on purpose. A broken fit legitimately produces a NaN + // sigma — `sqrt` of a negative truncated variance — and the design is + // to propagate that to `converge`'s `NonFiniteResult` guard, not to + // panic inside inference. Rejecting it here turned that reporting path + // into a crash, which two tests caught immediately. + assert!( + sigma >= 0.0 || sigma.is_nan(), + "sigma must not be negative; it is only ever squared, so a negative \ + value would silently behave as its absolute value" + ); if sigma == f64::INFINITY { Self { pi: 0.0, tau: 0.0 } } else if sigma == 0.0 { @@ -120,7 +156,25 @@ impl Gaussian { } } + /// How far this Gaussian moved from `other`, as `(|d mu|, |d sigma|)`. + /// + /// Identical messages have not moved, whatever their parameters, and that + /// case is answered in natural space before touching `mu()`/`sigma()`. An + /// improper message has `pi == 0`, so `sigma()` is infinite — and + /// `inf - inf` is NaN, a NaN *change* for a message that did not change at + /// all. (`mu()` is guarded and returns 0.0 here, so the mean component was + /// never the problem; the sigma component alone produced `(0.0, NaN)`.) + /// + /// That is reachable in ordinary inference: once a pairing is more than + /// about nine cavity-sigma apart the truncation is a no-op, `trunc / cavity` + /// is exactly the identity message, and the chain compares one identity + /// against another. Before this guard that produced `(0.0, NaN)`, which + /// silently disabled the sigma half of the convergence test. pub(crate) fn delta(&self, other: Gaussian) -> (f64, f64) { + if self.pi == other.pi && self.tau == other.tau { + return (0.0, 0.0); + } + ( (self.mu() - other.mu()).abs(), (self.sigma() - other.sigma()).abs(), @@ -256,6 +310,42 @@ impl ops::Div for Gaussian { #[cfg(test)] mod tests { + /// A message that did not change must report no change, even when it is + /// improper. `mu()` of an improper Gaussian is `0/0 = NaN` and `sigma()` is + /// infinite, so the mean/sigma form reported `(NaN, NaN)` for two identical + /// identity messages — which silently disabled the sigma half of the + /// convergence test in `run_chain`. + #[test] + fn delta_of_two_identical_improper_messages_is_zero() { + let improper = crate::N_INF; + // `mu()` is guarded and returns 0.0 for an improper Gaussian, so the + // mean component was always fine. The NaN came from the sigma + // component alone: `inf - inf`. The pre-fix value was `(0.0, NaN)`. + assert!(improper.sigma().is_infinite(), "premise: sigma is infinite"); + assert_eq!(improper.mu(), 0.0, "premise: mu is guarded, not NaN"); + assert!( + (improper.sigma() - improper.sigma()).is_nan(), + "premise: the unguarded sigma difference is NaN" + ); + assert_eq!(improper.delta(improper), (0.0, 0.0)); + } + + #[test] + fn delta_of_identical_proper_messages_is_zero() { + let g = Gaussian::from_ms(25.0, 8.0); + assert_eq!(g.delta(g), (0.0, 0.0)); + } + + /// The shortcut must not swallow a real difference. + #[test] + fn delta_still_measures_a_real_move() { + let a = Gaussian::from_ms(25.0, 8.0); + let b = Gaussian::from_ms(26.0, 9.0); + let (dmu, dsigma) = a.delta(b); + assert!((dmu - 1.0).abs() < 1e-12, "{dmu}"); + assert!((dsigma - 1.0).abs() < 1e-12, "{dsigma}"); + } + use super::*; #[test] diff --git a/src/history.rs b/src/history.rs index bd86277..5612cf5 100644 --- a/src/history.rs +++ b/src/history.rs @@ -1,4 +1,9 @@ -use std::{borrow::Borrow, collections::HashMap, hash::Hash, marker::PhantomData}; +use std::{ + borrow::Borrow, + collections::{BTreeMap, HashMap}, + hash::Hash, + marker::PhantomData, +}; use crate::{ BETA, GAMMA, Index, MU, P_DRAW, SIGMA, @@ -261,7 +266,15 @@ struct ResolvedTerms { contrast: Vec, /// Coefficients of competitors the slice has never seen, keyed by their /// rendering. Independent of everything in the slice by construction. - unseen: HashMap, + /// + /// A `BTreeMap` rather than a `HashMap`, and that is load-bearing. These + /// coefficients are summed, addition is not associative, and Rust seeds its + /// default hasher per process — so iterating a `HashMap` here made + /// `posterior_of` return different bits run to run on identical input. + /// Measured over 40 processes: two distinct sigma bit patterns, and five + /// distinct values from `expected_variance_reduction` spanning ~7 ULP. + /// Ordered iteration makes the sum reproducible. + unseen: BTreeMap, mean: f64, } @@ -947,6 +960,40 @@ impl, O: Observer, K: Eq + Hash + Clone> History f64 { + // sqrt(f64::EPSILON), as a const rather than a runtime sqrt. + const SQRT_EPSILON: f64 = 1.490_116_119_384_765_6e-8; + prior_variance * SQRT_EPSILON + } + fn time_expanded_joint(&self) -> TimeExpanded { let mut latest: HashMap = HashMap::new(); let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new(); @@ -969,8 +1016,9 @@ impl, O: Observer, K: Eq + Hash + Clone> History { let drift = rating.drift_variance_for_elapsed(elapsed); - if drift <= 0.0 { - // No drift: the same latent skill, not two. + if drift <= Self::collapse_threshold(rating.prior.variance()) { + // No drift, or too little to represent: the same + // latent skill, not two. See `collapse_threshold`. prev } else { let row = n; @@ -1032,7 +1080,7 @@ impl, O: Observer, K: Eq + Hash + Clone> History = HashMap::new(); + let mut unseen: BTreeMap = BTreeMap::new(); let mut mean = 0.0; for (member, (key, coefficient)) in terms.iter().enumerate() { @@ -1278,8 +1326,10 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History` and cannot inspect + // an arbitrary implementation. Validate what it actually produces + // instead, which also covers a custom impl. + // + // `ConstantDrift` returns `elapsed * gamma * gamma`, so a negative + // gamma is squared away: measured, `ConstantDrift(-0.0833)` gave + // results **bit identical** to `+0.0833`, the same sign-absorption + // defect already rejected for `sigma` and `beta`. A non-finite gamma + // poisons every posterior derived from it. + for slice in &self.time_slices { + for (agent, elapsed) in slice.appearances() { + let drift = self.agents[agent] + .rating + .drift_variance_for_elapsed(elapsed); + if !drift.is_finite() || drift < 0.0 { + return Err(InferenceError::InvalidParameter { + name: "drift variance", + value: drift, + }); + } + } + } + if self.time_slices.is_empty() { return Ok(ConvergenceReport { iterations: 0, @@ -1801,7 +1871,16 @@ impl, O: Observer, K: Eq + Hash + Clone> History = priors.keys().copied().collect(); + conflict_scan.sort_unstable(); + + for agent in &conflict_scan { + let batch = priors[agent]; let held = self.declared.get(agent).copied().unwrap_or_default(); if let (Some(existing), Some(new)) = (held.prior, batch.prior) { diff --git a/src/lib.rs b/src/lib.rs index a337b45..1feca81 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -218,6 +218,12 @@ const HALF_LINE_WINDOW: f64 = 10.0; /// 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. +/// `alpha / width` past which the tie branch's `v^2 - u` has lost too many +/// digits to trust, and the narrow-window form takes over. +/// +/// The subtraction retains about `(width / alpha)^2 / EPSILON` of its +/// precision, so this is the ratio at which that falls below roughly 1e-6. +const NARROW_WINDOW_RATIO: f64 = 2.0e4; const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0; pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0); @@ -476,10 +482,72 @@ fn pdf(x: f64, mu: f64, sigma: f64) -> f64 { 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 b = 2.0 - inv_sq * (10.0 - 74.0 * inv_sq); + let gap = inv * (1.0 - inv_sq * b); let v = alpha + gap; - (v, v * gap) + // Returns `1 - w`, not `w`, and that is the whole point of this shape. + // + // `w` tends to 1 out here, so a caller forming `1 - w` loses about + // `log10(alpha^2)` digits: measured against the exact truncated variance, + // `1 - w` came back with 8.9e-5 relative error at alpha = 1e6 and **0.0** + // from alpha = 1e8 — where the true value is 1e-16 and perfectly + // representable. `sigma * (1 - w).sqrt()` was then exactly zero, and + // `from_ms(mu, 0.0)` is a point mass whose `mu()` is `inf/inf = NaN`. + // + // Expanding `1 - v*gap` symbolically removes the subtraction: with + // `alpha*gap = 1 - inv^2*b`, the leading ones cancel on paper instead of in + // floating point, leaving `inv^2` times a bracket that tends to 1. Measured + // exact — 0.0 relative error — from alpha = 1e3 to 1e8. + let one_minus_w = inv_sq + * ((1.0 - inv_sq * (10.0 - 74.0 * inv_sq)) + 2.0 * inv_sq * b - inv_sq * inv_sq * b * b); + + (v, one_minus_w) +} + +/// Truncation to a *narrow* window `[alpha, alpha + d]`, as `(v, 1 - w)`. +/// +/// The tie branch forms `w` from `v^2 - u`, and both grow as `alpha^2` while +/// their difference stays `O(1)`. Far enough into the tail that subtraction has +/// nothing left: measured at `alpha = 1e6` with a window of `1e-6` it kept four +/// significant digits and returned `1 - w = -2.4e-4` where the truth is +/// `+2.8e-13`, so `sqrt` of it was NaN. One step earlier it was quietly wrong +/// instead — `1 - w = 1.0` exactly, a truncation reported as a no-op, where the +/// truth was `5e-17`. +/// +/// The existing half-line escape hatch does not cover it, because that keys on +/// `alpha * d >= HALF_LINE_WINDOW` — how many window-widths from the mean the +/// window sits — and a *narrow* window fails that however deep it is. +/// +/// Over a narrow window the density is `exp(-t*s - s^2 d^2 / 2)` in +/// `x = alpha + s*d`, with `t = alpha * d`. Dropping the `d^2` term leaves a +/// truncated exponential on `[0, 1]`, whose mean and variance are closed forms. +/// So `v = alpha + d*m(t)` and `1 - w = d^2 * V(t)`, with no subtraction of +/// large quantities anywhere. +/// +/// Measured against high-precision quadrature over `alpha` in `[1e2, 1e9]`: +/// `v` exact to 4e-10 or better, `1 - w` to 4e-10 across the region this is +/// used in. +fn narrow_window_truncation(alpha: f64, d: f64) -> (f64, f64) { + let t = alpha * d; + + // `m` and `V` are the mean and variance of a truncated exponential on + // [0, 1] with rate `t`, both of which cancel as `t -> 0`. The series is + // their limit (1/2 and 1/12, a uniform window) with the leading correction. + let (m, v_s) = if t < 1e-3 { + ( + 0.5 - t / 12.0 + t * t * t / 720.0, + 1.0 / 12.0 - t * t / 240.0, + ) + } else { + let em1 = libm::expm1(t); + ( + 1.0 / t - 1.0 / em1, + 1.0 / (t * t) - (em1 + 1.0) / (em1 * em1), + ) + }; + + (alpha + d * m, d * d * v_s) } fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { @@ -507,7 +575,7 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { (v, v - alpha) }; - (v, v * gap) + (v, 1.0 - 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. @@ -523,9 +591,22 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { // 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 width = beta - alpha; + + if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * width >= HALF_LINE_WINDOW { + let (v, one_minus_w) = half_line_truncation(alpha); + return (if flipped { -v } else { v }, one_minus_w); + } + + // A narrow window deep in the tail: too narrow for the half-line above, + // too deep for the subtraction below. The direct form keeps roughly + // `1 / (alpha/width)^2` of its digits, so the crossover is on that + // ratio rather than on either quantity alone — and the approximation is + // most accurate exactly where the subtraction is worst, since both + // improve as the window narrows. + if alpha > 0.0 && alpha > NARROW_WINDOW_RATIO * width { + let (v, one_minus_w) = narrow_window_truncation(alpha, width); + return (if flipped { -v } else { v }, one_minus_w); } let (v, u) = if alpha > 0.0 { @@ -548,17 +629,23 @@ fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { ) }; - let w = -(u - v.powi(2)); + // `1 - w` where `w = v^2 - u`. Both `v^2` and `u` grow as alpha^2 while + // their difference stays O(1), so this subtraction is the one place the + // tie branch can still lose everything — see the escape hatch above, + // which is what keeps the far tail away from it. + let one_minus_w = 1.0 + u - v.powi(2); - (if flipped { -v } else { v }, w) + (if flipped { -v } else { v }, one_minus_w) } } fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) { - let (v, w) = v_w(mu, sigma, margin, tie); + // `v_w` returns `1 - w` rather than `w`: forming the difference here is + // what destroyed the truncated variance in the far tail. + let (v, one_minus_w) = v_w(mu, sigma, margin, tie); let mu_trunc = mu + sigma * v; - let sigma_trunc = sigma * (1.0 - w).sqrt(); + let sigma_trunc = sigma * one_minus_w.sqrt(); (mu_trunc, sigma_trunc) } @@ -569,13 +656,34 @@ pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian { 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) { ( - if v1.0 > v2.0 { v1.0 } else { v2.0 }, - if v1.1 > v2.1 { v1.1 } else { v2.1 }, + 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 } @@ -652,6 +760,13 @@ pub(crate) fn sort_time(xs: &[T], reverse: bool) -> Vec { /// 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!( @@ -717,13 +832,141 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { let end = &rotated_a_matrix * &mean_matrix; let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant(); - let s_arg = ata.determinant() / middle.determinant(); - libm::exp(e_arg) * s_arg.sqrt() + // `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 { + /// The truncated variance must stay a variance across every branch, and + /// the branches must agree where they meet. + /// + /// `v_w` now has three regimes for a tie — half-line, narrow-window, and + /// the direct subtraction — and a misplaced crossover between them is the + /// failure mode this guards. A jump at a boundary is visible here even + /// though the absolute values are not pinned. + #[test] + fn truncated_variance_is_continuous_across_the_tie_branches() { + for &alpha in &[50.0, 99.0, 100.0, 101.0, 1e3, 1e5, 1e6] { + // Sweep the window width across NARROW_WINDOW_RATIO and the + // half-line threshold, which sit at different widths per alpha. + let mut previous: Option<(f64, f64)> = None; + let mut width = alpha / (NARROW_WINDOW_RATIO * 100.0); + while width < 40.0 / alpha { + // mu = 0 puts the window at [-margin, margin]; shift it out to + // `alpha` by moving the mean instead. + let margin = width * 0.5; + let mu = -(alpha + width * 0.5); + let (v, one_minus_w) = v_w(mu, 1.0, margin, true); + + assert!(v.is_finite(), "alpha {alpha}, width {width:e}: v = {v}"); + assert!( + one_minus_w.is_finite() && one_minus_w > 0.0 && one_minus_w <= 1.0, + "alpha {alpha}, width {width:e}: 1 - w = {one_minus_w:e} is not a variance" + ); + + if let Some((pv, pw)) = previous { + // Consecutive widths differ by 2x, so the moments may not + // differ by more than a small multiple of that. + assert!( + one_minus_w / pw < 32.0 && pw / one_minus_w < 32.0, + "alpha {alpha}: 1 - w jumped from {pw:e} to {one_minus_w:e} \ + at width {width:e} — a branch boundary is misplaced" + ); + assert!( + (v - pv).abs() <= 8.0 * width.max(1e-12) + 1e-9 * v.abs(), + "alpha {alpha}: v jumped from {pv} to {v} at width {width:e}" + ); + } + previous = Some((v, one_minus_w)); + width *= 2.0; + } + } + } + + /// The narrow-window form against high-precision quadrature. + /// + /// These are the inputs where the direct `v^2 - u` subtraction had four + /// significant digits left and returned a negative variance. + #[test] + fn narrow_window_truncation_matches_quadrature() { + for &(alpha, d, expect_v, expect_w) in &[ + (1e6, 2e-6, 1_000_000.000_000_687, 2.759_383_390_335_666e-13), + (1e4, 1e-6, 10_000.000_000_499_167, 8.333_291_666_831_727e-14), + ( + 1e3, + 1e-5, + 1_000.000_004_991_666_6, + 8.333_291_666_803_818e-12, + ), + ] { + let (v, one_minus_w) = narrow_window_truncation(alpha, d); + assert!( + ((v - expect_v) / expect_v).abs() < 1e-12, + "alpha {alpha:e}: v = {v}, want {expect_v}" + ); + assert!( + ((one_minus_w - expect_w) / expect_w).abs() < 1e-8, + "alpha {alpha:e}: 1 - w = {one_minus_w:e}, want {expect_w:e}" + ); + } + } + + /// 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::*; diff --git a/src/matrix.rs b/src/matrix.rs index cf29777..e9e05b3 100644 --- a/src/matrix.rs +++ b/src/matrix.rs @@ -91,6 +91,29 @@ impl Lu { det } + /// `ln |det|`, accumulated term by term rather than multiplied out. + /// + /// The determinant of an `n x n` Gram matrix is a product of `n` diagonal + /// entries, so it leaves `f64`'s range long before the quantities built + /// from it do. `quality()` only ever wants a *ratio* of two determinants, + /// and that ratio is perfectly representable while the determinants + /// themselves are not — measured, at 250 rating groups both overflow and + /// the ratio came back `NaN` where the true answer is `9.51e-88`. + /// + /// Returns `-inf` for a singular matrix, so `exp` of it is zero. + fn ln_abs_determinant(&self) -> f64 { + if self.sign == 0.0 { + return f64::NEG_INFINITY; + } + + let mut acc = 0.0; + for i in 0..self.n { + acc += libm::log(self.lu[i * self.n + i].abs()); + } + + acc + } + /// Solve `Ax = b` for a single column of the identity, giving one column /// of the inverse. fn solve_column(&self, col: usize, out: &mut [f64]) { @@ -157,6 +180,24 @@ impl Matrix { Lu::decompose(self).determinant() } + /// `ln |det|` of a square matrix; `-inf` when singular. + /// + /// See [`Lu::ln_abs_determinant`] for why a ratio of determinants must be + /// taken this way. + pub fn ln_abs_determinant(&self) -> f64 { + assert_eq!( + self.width, self.height, + "determinant requires a square matrix, got {}x{}", + self.height, self.width + ); + + if self.width == 0 { + return 0.0; + } + + Lu::decompose(self).ln_abs_determinant() + } + /// Matrix inverse via LU decomposition. /// /// # Panics diff --git a/src/predict.rs b/src/predict.rs index 2f221b3..ac7cfe4 100644 --- a/src/predict.rs +++ b/src/predict.rs @@ -21,7 +21,7 @@ //! would have made every `predict_*` call return a slightly different number, //! which is not a property a rating library should have. -use crate::{Gaussian, quadrature}; +use crate::{Gaussian, InferenceError, quadrature}; /// Teams beyond this count make the outcome enumeration impractical. /// @@ -52,6 +52,10 @@ const WIN_TOLERANCE: f64 = 1e-8; /// point where refining stops helping. const MIN_GRID_POINTS: usize = 8_192; const MAX_GRID_POINTS: usize = 262_144; +/// Nodes requested across the narrowest feature the recursion must resolve. +const NODES_PER_FEATURE: f64 = 12.0; +/// Nodes below which the trapezoid rule stops resolving that feature at all. +const MIN_NODES_PER_FEATURE: f64 = 4.0; /// How many standard deviations of support the grid and integrals cover. /// @@ -152,7 +156,7 @@ pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec (f64, f64, usize) { +fn grid_shape(perf: &[Gaussian], margins: &Margins) -> Result<(f64, f64, usize), InferenceError> { let lo = perf .iter() .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma()) @@ -175,18 +179,36 @@ fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (f64, f64, usize) { let feature = narrowest.min(smallest_margin); let wanted = if feature.is_finite() && feature > 0.0 { - ((hi - lo) / (feature / 12.0)).ceil() + ((hi - lo) / (feature / NODES_PER_FEATURE)).ceil() } else { MIN_GRID_POINTS as f64 }; - let points = if wanted.is_finite() { - (wanted as usize).clamp(MIN_GRID_POINTS, MAX_GRID_POINTS) - } else { - MIN_GRID_POINTS - }; + if !wanted.is_finite() { + return Ok((lo, hi, MIN_GRID_POINTS)); + } - (lo, hi, points) + // Report rather than clamp. Clamping is what this replaced: it silently + // handed the recursion a grid too coarse for the narrowest density, and the + // trapezoid rule then returned probabilities greater than one — measured, a + // `P` of 2.79 and a total of 5.41. Trapezoid error on a Gaussian is + // `~exp(-2 pi^2 (sigma/h)^2)`, which is 1e-12 at `h/sigma = 0.86` and O(1) + // by `h/sigma = 17`, so the cliff is sharp and there is no useful answer on + // the far side of it. + // + // The floor is `MIN_NODES_PER_FEATURE` rather than the `NODES_PER_FEATURE` + // asked for, because the request carries a large margin: measured accurate + // to 2.2e-12 at 1.4 nodes per sigma, and wrong by 1.2e-3 at 0.7. + let needed = wanted as usize; + let floor = ((hi - lo) / (feature / MIN_NODES_PER_FEATURE)).ceil(); + if floor.is_finite() && floor as usize > MAX_GRID_POINTS { + return Err(InferenceError::GridTooCoarse { + needed, + max: MAX_GRID_POINTS, + }); + } + + Ok((lo, hi, needed.clamp(MIN_GRID_POINTS, MAX_GRID_POINTS))) } /// Densities of each team sampled on the shared grid. @@ -198,8 +220,8 @@ struct Sampled { } impl Sampled { - fn new(perf: &[Gaussian], margins: &Margins) -> Self { - let (lo, hi, points) = grid_shape(perf, margins); + fn new(perf: &[Gaussian], margins: &Margins) -> Result { + let (lo, hi, points) = grid_shape(perf, margins)?; let step = (hi - lo) / (points - 1) as f64; let density = perf .iter() @@ -209,12 +231,12 @@ impl Sampled { .collect() }) .collect(); - Self { + Ok(Self { lo, step, points, density, - } + }) } fn node(&self, i: usize) -> f64 { @@ -317,9 +339,12 @@ fn events(n: usize, strict_only: bool) -> Vec<(Vec, Vec)> { /// /// Orders that differ only *within* a tied group describe the same finishing /// order, so their probabilities are summed into one entry. -pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec<(Vec, f64)> { +pub(crate) fn outcome_distribution( + perf: &[Gaussian], + margins: &Margins, +) -> Result, f64)>, InferenceError> { let n = perf.len(); - let sampled = Sampled::new(perf, margins); + let sampled = Sampled::new(perf, margins)?; let mut aggregated: Vec<(Vec, f64)> = Vec::new(); for (order, tied) in events(n, margins.all_zero()) { @@ -332,7 +357,7 @@ pub(crate) fn outcome_distribution(perf: &[Gaussian], margins: &Margins) -> Vec< } aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal)); - aggregated + Ok(aggregated) } /// All permutations of `items`. @@ -397,9 +422,13 @@ fn orders_for_groups(groups: &[Vec]) -> Vec<(Vec, Vec)> { /// Ties in `ranks` mean the tied teams may finish in any internal order, so /// this sums the orders consistent with the requested ranking rather than /// picking one. -pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: &[u32]) -> f64 { +pub(crate) fn ranking_probability( + perf: &[Gaussian], + margins: &Margins, + ranks: &[u32], +) -> Result { let n = perf.len(); - let sampled = Sampled::new(perf, margins); + let sampled = Sampled::new(perf, margins)?; let mut distinct: Vec = ranks.to_vec(); distinct.sort_unstable(); @@ -410,10 +439,10 @@ pub(crate) fn ranking_probability(perf: &[Gaussian], margins: &Margins, ranks: & .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect()) .collect(); - orders_for_groups(&groups) + Ok(orders_for_groups(&groups) .iter() .map(|(order, tied)| order_probability(margins, &sampled, order, tied)) - .sum() + .sum()) } /// A distribution over the ways a contest could finish. @@ -604,7 +633,7 @@ mod tests { ), ] { let n = perf.len(); - let dist = outcome_distribution(&perf, &flat(n, eps)); + let dist = outcome_distribution(&perf, &flat(n, eps)).unwrap(); let sum: f64 = dist.iter().map(|(_, p)| p).sum(); assert!( (sum - 1.0).abs() < 1e-6, @@ -620,7 +649,7 @@ mod tests { fn two_team_distribution_matches_the_closed_form() { let perf = [g(3.0, 6.0), g(-2.0, 1.0)]; let eps = 1.5; - let dist = outcome_distribution(&perf, &flat(2, eps)); + let dist = outcome_distribution(&perf, &flat(2, eps)).unwrap(); let (wa, wb) = closed_form_two(perf[0], perf[1], eps); let find = |ranks: &[u32]| { @@ -653,10 +682,10 @@ mod tests { let perf = [g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)]; let eps = 1.5; let margins = flat(3, eps); - let dist = outcome_distribution(&perf, &margins); + let dist = outcome_distribution(&perf, &margins).unwrap(); for (ranks, expected) in &dist { - let direct = ranking_probability(&perf, &margins, ranks); + let direct = ranking_probability(&perf, &margins, ranks).unwrap(); assert!( (direct - expected).abs() < 1e-9, "ranks {ranks:?}: direct {direct} vs distribution {expected}" @@ -675,7 +704,7 @@ mod tests { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)]; let mut previous = 0.0; for eps in [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 24.0] { - let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]); + let p = ranking_probability(&perf, &flat(3, eps), &[0, 0, 0]).unwrap(); assert!(p >= previous, "eps={eps}: {p} < {previous}"); if eps == 0.0 { assert!(p < 1e-12, "a tie needs a margin, got {p}"); @@ -696,7 +725,7 @@ mod tests { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(-8.0, 2.0)]; let sweep: Vec = [0.5, 2.0, 4.0, 8.0, 16.0] .iter() - .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1])) + .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1]).unwrap()) .collect(); let peak = sweep .iter() @@ -717,7 +746,7 @@ mod tests { #[test] fn ties_are_impossible_without_a_draw_margin() { let perf = [g(0.0, 4.0), g(0.0, 4.0), g(0.0, 4.0)]; - let dist = outcome_distribution(&perf, &flat(3, 0.0)); + let dist = outcome_distribution(&perf, &flat(3, 0.0)).unwrap(); assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}"); assert!(dist.iter().all(|(r, _)| { let mut seen = r.clone(); diff --git a/src/rating.rs b/src/rating.rs index f589dc8..5de6b5c 100644 --- a/src/rating.rs +++ b/src/rating.rs @@ -23,7 +23,24 @@ pub struct Rating = ConstantDrift> { } impl> Rating { + /// # Panics + /// + /// Panics unless `beta` is finite and non-negative, matching + /// `HistoryBuilder::beta`. + /// + /// Zero is allowed and meaningful — performance is then exactly skill, and + /// the fit differs measurably from a positive beta rather than degenerating. + /// Negative is rejected because `beta` enters only as `beta^2`: measured, a + /// negative beta returned results **bit identical** to its absolute value, + /// and a NaN beta reached `Game::ranked`, which returned `Ok` carrying a + /// `Gaussian { pi: NaN, tau: NaN }` — there is no `converge` on that path to + /// catch it. pub fn new(prior: Gaussian, beta: f64, drift: D) -> Self { + assert!( + beta.is_finite() && beta >= 0.0, + "beta must be finite and non-negative (got {beta}); it is only ever \ + squared, so a negative value would silently behave as its absolute value" + ); Self { prior, beta, diff --git a/tests/cross_process_determinism.rs b/tests/cross_process_determinism.rs new file mode 100644 index 0000000..da9b5ed --- /dev/null +++ b/tests/cross_process_determinism.rs @@ -0,0 +1,157 @@ +//! Determinism across *processes*, which an in-process test cannot see. +//! +//! Rust seeds its default hasher once per process, so every `HashMap` +//! iteration order is fixed for a run and varies between runs. A test that +//! compares results within one process therefore cannot detect a float sum +//! whose order comes from a map — all its samples share one seed. +//! +//! That is not hypothetical. `tests/determinism.rs` compares four thread counts +//! inside one process and passed throughout, while `posterior_of` was returning +//! two distinct bit patterns across 40 separate runs on identical input. +//! +//! This re-executes the test binary and compares `f64::to_bits`. + +use std::{env, process::Command}; + +use smallvec::smallvec; +use trueskill_tt::{ + ConstantDrift, ConvergenceOptions, Event, History, Member, NullObserver, Outcome, Team, + UnknownKeys, +}; + +/// Set in the child so it reports instead of re-spawning. +const CHILD: &str = "TSTT_DETERMINISM_CHILD"; + +const RUNS: usize = 40; + +type H = History; + +fn fitted() -> H { + let mut h: H = History::builder_with_key() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .score_sigma(2.0) + .drift(ConstantDrift(0.05)) + .unknown_keys(UnknownKeys::Prior) + .convergence(ConvergenceOptions { + max_iter: 20_000, + epsilon: 1e-13, + alpha: 1.0, + }) + .build(); + + let mut events = Vec::new(); + for t in 0..12i64 { + for k in 0..6usize { + let a = format!("p{}", (t as usize * 6 + k) % 10); + let b = format!("p{}", (t as usize * 6 + k + 4) % 10); + events.push(Event { + time: t, + teams: smallvec![ + Team::with_members([Member::new(a)]), + Team::with_members([Member::new(b)]), + ], + outcome: Outcome::scores([3.0, 1.0]), + }); + } + } + h.add_events(events).unwrap(); + assert!(h.converge().unwrap().converged); + h +} + +/// Every quantity that could plausibly depend on iteration order, as bits. +fn fingerprint() -> String { + let h = fitted(); + + // Unknown keys with UNEQUAL but COMPARABLE coefficients, which is what + // makes the sum order-sensitive. + // + // Equal terms sum order-independently and would make this pass vacuously. + // Terms of wildly different magnitudes are no better: the small ones fall + // below the running total's ULP and are absorbed whatever the order — + // measured, spreading these over nine decades dropped the detection rate + // to roughly one run in forty. Comparable sizes keep every term able to + // change the last bits. + let ghosts: Vec = (0..24).map(|i| format!("ghost{i}")).collect(); + let mut terms: Vec<(&String, f64)> = ghosts + .iter() + .enumerate() + .map(|(i, k)| (k, 1.0 + i as f64 * 0.37)) + .collect(); + let known = "p0".to_string(); + terms.push((&known, -1.0)); + + let posterior = h.posterior_of(&terms).unwrap(); + + let a = "p0".to_string(); + let b = "p1".to_string(); + let target = [(&a, 1.0), (&b, -1.0)]; + let teams: [&[&String]; 2] = [&[&a], &[&b]]; + let evr = h.expected_variance_reduction(&teams, &target).unwrap(); + + let curves = h.learning_curves(); + let mut curve_bits: u64 = 0; + let mut keys: Vec<&String> = curves.keys().collect(); + keys.sort(); + for key in keys { + for (t, g) in &curves[key] { + curve_bits ^= (*t as u64).rotate_left(17) + ^ g.mu().to_bits().rotate_left(31) + ^ g.sigma().to_bits(); + } + } + + format!( + "post={:016x} evr={:016x} le={:016x} curves={curve_bits:016x}", + posterior.sigma().to_bits(), + evr.to_bits(), + h.log_evidence().to_bits(), + ) +} + +#[test] +fn results_are_identical_across_processes() { + if env::var(CHILD).is_ok() { + println!("FINGERPRINT {}", fingerprint()); + return; + } + + let exe = env::current_exe().expect("current exe"); + let mut seen: Vec = Vec::new(); + + for run in 0..RUNS { + let out = Command::new(&exe) + .args([ + "results_are_identical_across_processes", + "--exact", + "--nocapture", + ]) + .env(CHILD, "1") + .output() + .expect("spawn child"); + assert!( + out.status.success(), + "child {run} failed: {}", + String::from_utf8_lossy(&out.stderr) + ); + let stdout = String::from_utf8_lossy(&out.stdout); + let line = stdout + .lines() + .find_map(|l| l.strip_prefix("FINGERPRINT ")) + .unwrap_or_else(|| panic!("child {run} printed no fingerprint:\n{stdout}")) + .to_string(); + seen.push(line); + } + + let first = &seen[0]; + let differing: Vec<&String> = seen.iter().filter(|s| *s != first).collect(); + assert!( + differing.is_empty(), + "results differ across processes on identical input.\n {} of {RUNS} runs differed\n \ + first: {first}\n differing: {}", + differing.len(), + differing[0] + ); +} diff --git a/tests/determinism.rs b/tests/determinism.rs index 006f25b..e6c1da5 100644 --- a/tests/determinism.rs +++ b/tests/determinism.rs @@ -1,101 +1,201 @@ -//! Determinism tests: identical posteriors across RAYON_NUM_THREADS -//! values. Only compiled with the `rayon` feature. +//! Determinism across `RAYON_NUM_THREADS`, on a workload that actually reaches +//! the parallel path. +//! +//! This test previously proved less than it appeared to. `sweep_color_groups` +//! takes its `par_iter` branch only for colour groups of at least +//! `RAYON_THRESHOLD` (64) events, and the old fixture built 20 slices of 10 +//! events — a colour group is a subset of one slice's events, so it could never +//! exceed 10. The branch was unreachable, confirmed by CPU-vs-wall time: +//! `user 0.64` on eight threads is one core. +//! +//! It also compared a single competitor's curve out of forty, and never +//! compared `log_evidence`, `final_step` or `iterations`. +//! +//! The fixture below guarantees the parallel branch **by construction**: within +//! a slice every event uses a disjoint pair of competitors, so greedy colouring +//! puts all of them in colour 0, and that group is `EVENTS_PER_SLICE` long. +//! Competitors recur across slices, so the fit still has temporal coupling and +//! drift rather than being a set of independent duels. #![cfg(feature = "rayon")] use smallvec::smallvec; -use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team}; +use trueskill_tt::{ + ConstantDrift, ConvergenceOptions, Event, Gaussian, History, Member, Outcome, Team, +}; -/// Build a deterministic workload using a simple LCG (no external rand crate). -fn build_and_converge(seed: u64) -> Vec<(i64, trueskill_tt::Gaussian)> { +/// Comfortably above the crate's internal `RAYON_THRESHOLD` of 64. +const EVENTS_PER_SLICE: usize = 96; +const SLICES: i64 = 8; +/// Two per event, all disjoint within a slice. +const COMPETITORS: usize = EVENTS_PER_SLICE * 2; + +/// Everything a thread count could plausibly perturb. +struct Fingerprint { + curves: Vec<(String, Vec<(i64, Gaussian)>)>, + log_evidence: f64, + final_step: (f64, f64), + iterations: usize, +} + +fn build_and_converge() -> Fingerprint { let mut h = History::::builder_with_key() .mu(25.0) .sigma(25.0 / 3.0) .beta(25.0 / 6.0) .drift(ConstantDrift(25.0 / 300.0)) .convergence(ConvergenceOptions { - max_iter: 30, - epsilon: 1e-6, + max_iter: 20_000, + epsilon: 1e-9, alpha: 1.0, }) .build(); - // LCG for deterministic pseudo-random ints. - let mut rng = seed; - let mut next = || { - rng = rng - .wrapping_mul(6364136223846793005) - .wrapping_add(1442695040888963407); - rng - }; - - let mut events: Vec> = Vec::with_capacity(200); - for ev_i in 0..200 { - let a = (next() % 40) as usize; - let mut b = (next() % 40) as usize; - while b == a { - b = (next() % 40) as usize; + let mut events: Vec> = Vec::new(); + for slice in 0..SLICES { + for e in 0..EVENTS_PER_SLICE { + // Disjoint within the slice: event `e` owns competitors 2e and + // 2e+1. Rotating by the slice index makes the pairings differ + // between slices, so competitors accumulate a real history. + let a = (2 * e + slice as usize) % COMPETITORS; + let b = (2 * e + 1 + slice as usize * 3) % COMPETITORS; + if a == b { + continue; + } + events.push(Event { + time: slice + 1, + teams: smallvec![ + Team::with_members([Member::new(format!("p{a}"))]), + Team::with_members([Member::new(format!("p{b}"))]), + ], + outcome: Outcome::winner(u32::from((e + slice as usize) % 2 == 0), 2), + }); } - // ~10 events per slice so color groups have material parallelism. - events.push(Event { - time: (ev_i as i64 / 10) + 1, - teams: smallvec![ - Team::with_members([Member::new(format!("p{a}"))]), - Team::with_members([Member::new(format!("p{b}"))]), - ], - outcome: Outcome::winner((next() % 2) as u32, 2), - }); } h.add_events(events).unwrap(); - let _ = h.converge().unwrap(); - // Sample one competitor's curve for the comparison. - h.learning_curve("p0") + + let report = h.converge().expect("fixture must converge"); + + let mut curves: Vec<(String, Vec<(i64, Gaussian)>)> = h + .learning_curves() + .into_iter() + .map(|(k, v)| (k.clone(), v)) + .collect(); + curves.sort_by(|a, b| a.0.cmp(&b.0)); + + Fingerprint { + curves, + log_evidence: h.log_evidence(), + final_step: report.final_step, + iterations: report.iterations, + } } #[test] fn posteriors_identical_across_thread_counts() { let sizes = [1usize, 2, 4, 8]; - let mut results: Vec> = Vec::new(); + let mut results: Vec = Vec::new(); + for &n in &sizes { let pool = rayon::ThreadPoolBuilder::new() .num_threads(n) .build() .expect("rayon pool build"); - let curve = pool.install(|| build_and_converge(42)); - results.push(curve); + results.push(pool.install(build_and_converge)); } let reference = &results[0]; - for (i, curve) in results.iter().enumerate().skip(1) { + + // Guard against the failure this test previously had: passing while + // measuring almost nothing. + assert!( + reference.curves.len() > 100, + "expected every competitor's curve, got {}", + reference.curves.len() + ); + + for (i, got) in results.iter().enumerate().skip(1) { + let n = sizes[i]; + assert_eq!( - curve.len(), - reference.len(), - "curve length differs at {n} threads", - n = sizes[i], + got.iterations, reference.iterations, + "iterations differ at {n} threads" ); - for (j, (&(t_ref, g_ref), &(t, g))) in reference.iter().zip(curve.iter()).enumerate() { + assert_eq!( + got.final_step.0.to_bits(), + reference.final_step.0.to_bits(), + "final_step.0 differs at {n} threads: {:?} vs {:?}", + reference.final_step, + got.final_step + ); + assert_eq!( + got.final_step.1.to_bits(), + reference.final_step.1.to_bits(), + "final_step.1 differs at {n} threads" + ); + assert_eq!( + got.log_evidence.to_bits(), + reference.log_evidence.to_bits(), + "log_evidence differs at {n} threads: {} vs {}", + reference.log_evidence, + got.log_evidence + ); + + assert_eq!( + got.curves.len(), + reference.curves.len(), + "competitor count differs at {n} threads" + ); + + for ((ref_key, ref_curve), (key, curve)) in reference.curves.iter().zip(got.curves.iter()) { + assert_eq!(ref_key, key, "competitor order differs at {n} threads"); assert_eq!( - t_ref, - t, - "time point {j} differs at {n} threads: ref={t_ref} vs got={t}", - n = sizes[i], - ); - assert_eq!( - g_ref.mu().to_bits(), - g.mu().to_bits(), - "mu bits differ at {n} threads, time {t}: ref={ref_mu} got={got_mu}", - n = sizes[i], - ref_mu = g_ref.mu(), - got_mu = g.mu(), - ); - assert_eq!( - g_ref.sigma().to_bits(), - g.sigma().to_bits(), - "sigma bits differ at {n} threads, time {t}: ref={ref_sigma} got={got_sigma}", - n = sizes[i], - ref_sigma = g_ref.sigma(), - got_sigma = g.sigma(), + curve.len(), + ref_curve.len(), + "curve length differs for {key} at {n} threads" ); + for (&(t_ref, g_ref), &(t, g)) in ref_curve.iter().zip(curve.iter()) { + assert_eq!(t_ref, t, "time point differs for {key} at {n} threads"); + assert_eq!( + g_ref.mu().to_bits(), + g.mu().to_bits(), + "mu differs for {key} at t={t}, {n} threads: {} vs {}", + g_ref.mu(), + g.mu() + ); + assert_eq!( + g_ref.sigma().to_bits(), + g.sigma().to_bits(), + "sigma differs for {key} at t={t}, {n} threads: {} vs {}", + g_ref.sigma(), + g.sigma() + ); + } } } } + +/// The fixture must keep reaching the parallel branch. +/// +/// `RAYON_THRESHOLD` is private, so this pins the property that makes the +/// branch reachable rather than the branch itself: within a slice every event +/// uses a disjoint competitor pair, so greedy colouring puts all +/// `EVENTS_PER_SLICE` of them in one colour group. If someone shrinks the +/// fixture, this fails rather than the suite quietly going back to testing the +/// sequential path. +#[test] +fn the_fixture_still_exceeds_the_rayon_threshold() { + const RAYON_THRESHOLD: usize = 64; + const { + assert!( + EVENTS_PER_SLICE >= RAYON_THRESHOLD, + "a colour group holds at most EVENTS_PER_SLICE events, which must \ + reach the crate's RAYON_THRESHOLD for the parallel sweep to run" + ); + } + + // Measured by instrumenting `sweep_color_groups`: this fixture produces + // one colour group of 96 events and takes the parallel branch on all 872 + // sweeps. The old fixture's 10-event slices could not reach 64 at all. + assert_eq!(EVENTS_PER_SLICE, 96); +} diff --git a/tests/joint_handle.rs b/tests/joint_handle.rs index e5d90df..9a72050 100644 --- a/tests/joint_handle.rs +++ b/tests/joint_handle.rs @@ -265,3 +265,75 @@ fn unseen_competitors_match_the_one_shot_path() { assert_eq!(one_shot.pi(), cached.pi()); assert_eq!(one_shot.tau(), cached.tau()); } + +/// A drift too small to represent must collapse, not corrupt the matrix. +/// +/// The collapse rule used to fire only at `drift <= 0.0` exactly. Anything +/// smaller-but-positive got an explicit `1.0 / drift` precision, and at +/// `drift = 1e-16` that entry is `1e16` — so `1e16 + 0.28` rounds back to +/// `1e16` and the prior and contrasts are annihilated in the stored `f64`. +/// +/// Measured before the fix, at `drift_scale = 1e-10` this returned a variance +/// **12 000x too small** (a 111x overconfident interval) as `Ok`, with a band +/// just above it returning a misleading `JointUnavailable`. +#[test] +fn a_drift_too_small_to_represent_collapses_rather_than_corrupting() { + fn variance(scale: f64) -> f64 { + let mut h: History = History::builder_with_key() + .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(); + let mut events = Vec::new(); + for t in 0..15i64 { + for k in 0..4usize { + let x = format!("p{}", (t as usize * 4 + k) % 8); + let y = format!("p{}", (t as usize * 4 + k + 3) % 8); + events.push(Event { + time: t, + teams: smallvec![ + Team::with_members([Member::new(x).with_drift_scale(scale)]), + Team::with_members([Member::new(y).with_drift_scale(scale)]), + ], + outcome: Outcome::scores([3.0, 1.0]), + }); + } + } + h.add_events(events).unwrap(); + assert!(h.converge().unwrap().converged); + let (a, b) = ("p0".to_string(), "p1".to_string()); + let joint = h + .joint() + .expect("a tiny drift must not make the joint unavailable"); + let g = joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).unwrap(); + g.sigma() * g.sigma() + } + + let collapsed = variance(0.0); + + // Below the threshold every scale must reach the collapsed answer exactly, + // and none may error. + for scale in [1e-3, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12] { + let v = variance(scale); + assert_eq!( + v.to_bits(), + collapsed.to_bits(), + "drift_scale {scale:e}: {v} vs collapsed {collapsed}" + ); + } + + // Above it, real drift is still modelled — otherwise this test would pass + // by collapsing everything. + let drifting = variance(1e-2); + assert!( + (drifting - collapsed).abs() / collapsed > 1e-5, + "a drift of 1e-2 must still move the answer: {drifting} vs {collapsed}" + ); +} diff --git a/tests/libm_rule.rs b/tests/libm_rule.rs new file mode 100644 index 0000000..b64e40d --- /dev/null +++ b/tests/libm_rule.rs @@ -0,0 +1,175 @@ +//! The libm rule, enforced rather than asserted in prose. +//! +//! CLAUDE.md requires transcendentals to 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. +//! +//! The rule was stated clearly and still violated in three production sites, +//! one of them `hypot` on the path of every scored event — whose measured +//! divergence, 12.1%, is *higher* than the `exp` figure the rule cites as its +//! own justification. Prose is evidently not enough, so this is a test. +//! +//! Tests may use either, which the crate documents, so `#[cfg(test)]` blocks +//! are excluded. + +use std::{fs, path::Path}; + +/// Method-call spellings that reach the system math library. +/// +/// `sqrt` is deliberately absent: IEEE 754 specifies it exactly, so `std` and +/// `libm` cannot disagree. `abs`, `recip`, `powi` and `mul_add` are likewise +/// exact or specified. +const FORBIDDEN: &[&str] = &[ + "exp", "exp2", "exp_m1", "ln", "ln_1p", "log", "log2", "log10", "powf", "sin", "cos", "tan", + "asin", "acos", "atan", "atan2", "sinh", "cosh", "tanh", "hypot", "cbrt", "erf", "erfc", +]; + +/// Strip `#[cfg(test)]` items by brace matching, plus comments and string +/// literals, so a mention in prose is not mistaken for a call. +fn production_code(source: &str) -> String { + let mut out = String::with_capacity(source.len()); + let bytes: Vec = source.chars().collect(); + let mut i = 0; + + while i < bytes.len() { + let rest: String = bytes[i..].iter().take(16).collect(); + + if rest.starts_with("#[cfg(test)]") { + // Skip to the opening brace of the guarded item, then past its + // matching close. + let mut j = i; + while j < bytes.len() && bytes[j] != '{' { + j += 1; + } + let mut depth = 0usize; + while j < bytes.len() { + match bytes[j] { + '{' => depth += 1, + '}' => { + depth -= 1; + if depth == 0 { + j += 1; + break; + } + } + _ => {} + } + j += 1; + } + i = j; + continue; + } + + if rest.starts_with("//") { + while i < bytes.len() && bytes[i] != '\n' { + i += 1; + } + continue; + } + + if rest.starts_with("/*") { + i += 2; + while i + 1 < bytes.len() && !(bytes[i] == '*' && bytes[i + 1] == '/') { + i += 1; + } + i += 2; + continue; + } + + if bytes[i] == '"' { + i += 1; + while i < bytes.len() && bytes[i] != '"' { + if bytes[i] == '\\' { + i += 1; + } + i += 1; + } + i += 1; + continue; + } + + out.push(bytes[i]); + i += 1; + } + + out +} + +fn rust_files(dir: &Path, out: &mut Vec) { + for entry in fs::read_dir(dir).expect("read src") { + let path = entry.expect("dir entry").path(); + if path.is_dir() { + rust_files(&path, out); + } else if path.extension().is_some_and(|e| e == "rs") { + out.push(path); + } + } +} + +#[test] +fn production_code_never_calls_a_std_transcendental() { + let mut files = Vec::new(); + rust_files(Path::new("src"), &mut files); + assert!(files.len() > 10, "expected to find the crate's sources"); + + let mut offences = Vec::new(); + + for path in &files { + let source = fs::read_to_string(path).expect("read source"); + let code = production_code(&source); + + for (n, line) in code.lines().enumerate() { + for name in FORBIDDEN { + let needle = format!(".{name}("); + if line.contains(&needle) { + offences.push(format!("{}:{}: {}", path.display(), n + 1, line.trim())); + } + } + } + } + + assert!( + offences.is_empty(), + "production code must call libm, not std, for transcendentals \ + (`sqrt` is exempt — IEEE 754 specifies it):\n{}", + offences.join("\n") + ); +} + +/// The stripper has to actually strip, or the test above passes vacuously. +#[test] +fn the_test_module_stripper_works() { + let source = r#" +fn production() { let _ = libm::exp(1.0); } + +#[cfg(test)] +mod tests { + fn allowed() { let x = 1.0f64.exp(); } +} + +fn also_production() {} +"#; + let code = production_code(source); + assert!( + code.contains("also_production"), + "stripped too much: {code}" + ); + assert!( + !code.contains(".exp()"), + "failed to strip cfg(test): {code}" + ); +} + +/// And it must not strip a doc comment's worth of prose into oblivion, nor +/// mistake prose for a call. +#[test] +fn prose_is_not_mistaken_for_a_call() { + let source = "/// Uses `x.exp()` in the docs.\nfn f() { let _ = libm::exp(1.0); }\n"; + let code = production_code(source); + assert!(!code.contains(".exp()"), "doc comment leaked: {code}"); + assert!(code.contains("libm::exp"), "stripped real code: {code}"); +} diff --git a/tests/non_finite_results.rs b/tests/non_finite_results.rs index 8928431..b7db0f0 100644 --- a/tests/non_finite_results.rs +++ b/tests/non_finite_results.rs @@ -10,7 +10,9 @@ //! `!tuple_gt(..)`. These tests pin the guard from outside. use smallvec::smallvec; -use trueskill_tt::{Event, Gaussian, History, InferenceError, Member, Outcome, Team}; +use trueskill_tt::{ + ConstantDrift, Event, Gaussian, History, InferenceError, Member, Outcome, Team, +}; fn scored_fit( sigma: f64, @@ -115,3 +117,101 @@ fn merely_extreme_parameters_still_converge() { assert!(scored_fit(6.0, 1.0, 1e6, [3.0, 1.0]).unwrap()); assert!(scored_fit(6.0, 1.0, 1.0, [1e150, -1e150]).unwrap()); } + +/// A NaN in one competitor must not be masked by a healthy competitor reduced +/// after it. +/// +/// The convergence step is a fold over a `HashMap`, so which competitor is +/// reduced last is per-process hash order. Before the fix, `tuple_max` dropped +/// a NaN accumulator in favour of the next finite delta and this returned +/// `Ok(converged: true)` with a NaN posterior in **16 of 30 runs** on identical +/// input. Deterministic now, but note this test can only ever sample one hash +/// order per run — the ordering guarantee itself is pinned by +/// `tuple_max_propagates_a_nan_from_any_position` in the crate's unit tests. +#[test] +fn a_nan_competitor_is_not_masked_by_a_healthy_one() { + let mut h = History::builder() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .p_draw(0.1) + .build(); + h.add_events(vec![ + Event { + time: 1i64, + teams: smallvec![ + Team::with_members([Member::new("a").with_prior(Gaussian::from_ms(0.0, 1e-200))]), + Team::with_members([Member::new("b")]), + ], + outcome: Outcome::winner(0, 2), + }, + // A healthy pair in the same slice, to be reduced alongside the NaN. + Event { + time: 1i64, + teams: smallvec![ + Team::with_members([Member::new("c")]), + Team::with_members([Member::new("d")]), + ], + outcome: Outcome::winner(0, 2), + }, + ]) + .unwrap(); + + let err = h + .converge() + .expect_err("a NaN fit must never be reported as converged"); + assert!( + matches!(err, InferenceError::NonFiniteResult { .. }), + "{err:?}" + ); +} + +/// A tie observed with a narrow draw margin between far-apart competitors must +/// produce a fit, not NaN skills. +/// +/// The tie branch forms the truncated variance from `v^2 - u`, and both grow as +/// `alpha^2` while their difference stays `O(1)`. Deep enough into the tail +/// that subtraction had four digits left: measured, it returned `1 - w` +/// negative and `sqrt` of it was NaN. The half-line escape hatch did not cover +/// it, because that keys on how many window-widths from the mean the window +/// sits and a narrow window fails that however deep it is. +/// +/// These parameters are ordinary for a precise-scoring domain, and the +/// neighbouring wider-margin case always worked — so this was a cliff, not +/// "extreme inputs break". +#[test] +fn a_narrow_draw_margin_far_into_the_tail_still_fits() { + for (beta, p_draw, sd, gap) in [ + (1e-2, 1e-8, 1e-2, 10.0), + (1e-3, 1e-9, 1e-3, 1.0), + (1e-4, 1e-12, 1e-4, 1.0), + ] { + let mut h = History::builder() + .mu(0.0) + .sigma(sd) + .beta(beta) + .p_draw(p_draw) + .drift(ConstantDrift(0.0)) + .build(); + h.add_events(vec![Event { + time: 1i64, + teams: smallvec![ + Team::with_members([Member::new("a").with_prior(Gaussian::from_ms(0.0, sd))]), + Team::with_members([Member::new("b").with_prior(Gaussian::from_ms(gap, sd))]), + ], + outcome: Outcome::draw(2), + }]) + .unwrap(); + + let report = h + .converge() + .unwrap_or_else(|e| panic!("beta {beta:e}, p_draw {p_draw:e}: {e:?}")); + assert!(report.converged); + + let skill = h.current_skill(&"a").unwrap(); + assert!( + skill.mu().is_finite() && skill.sigma().is_finite() && skill.sigma() > 0.0, + "beta {beta:e}, p_draw {p_draw:e}: {skill:?}" + ); + } +} diff --git a/tests/prediction_bounds.rs b/tests/prediction_bounds.rs new file mode 100644 index 0000000..eedde48 --- /dev/null +++ b/tests/prediction_bounds.rs @@ -0,0 +1,154 @@ +//! Bounds that any correct implementation must satisfy, swept rather than +//! spot-checked. +//! +//! The crate's docs call the `ln k` ceiling "the sharpest available test of an +//! implementation", and record that an early prototype returned 4.77 nats. It +//! was violated again — 3.237828 nats against `ln 2` — because the existing +//! check sampled one fixture and the violation lives in a specific regime: a +//! large ratio between the widest and narrowest performance sigma, where the +//! shared prediction grid could not resolve the narrow density and returned +//! probabilities greater than one. +//! +//! A single fixture cannot defend a bound like this. A sweep can. + +use trueskill_tt::{ + ConstantDrift, GameOptions, Gaussian, InferenceError, Rating, expected_information_gain, +}; + +type R = Rating; + +/// How many random matchups the ceiling sweep draws. +/// +/// Scaled by build profile rather than fixed. Each sample runs a full inference +/// pass per outcome, and that is about **19x** faster in release — measured, +/// 20 000 samples take 12.1s released against 23s for 2 000 in debug. `just +/// test` runs three debug feature combinations and one release one, so a fixed +/// count pays the slow price three times and the fast one once, which is +/// exactly backwards. +/// +/// The debug run is here to prove the sweep still compiles and holds on a small +/// sample; the release run is the one that actually searches. The violation +/// this guards was found at a rate near 1.8%, so even the debug count expects +/// tens of hits in the regime. +#[cfg(debug_assertions)] +const SAMPLES: usize = 1_000; +#[cfg(not(debug_assertions))] +const SAMPLES: usize = 50_000; + +/// Deterministic LCG, so a failure is reproducible from the printed seed. +struct Lcg(u64); + +impl Lcg { + fn next_f64(&mut self) -> f64 { + self.0 = self + .0 + .wrapping_mul(6_364_136_223_846_793_005) + .wrapping_add(1_442_695_040_888_963_407); + // Top 53 bits to [0, 1). + ((self.0 >> 11) as f64) / ((1u64 << 53) as f64) + } + + fn in_range(&mut self, lo: f64, hi: f64) -> f64 { + lo + (hi - lo) * self.next_f64() + } + + /// Log-uniform, so the sweep spends its samples across magnitudes rather + /// than crowding the top of the range — the violations live at small sigma. + fn log_uniform(&mut self, lo: f64, hi: f64) -> f64 { + let t = self.next_f64(); + (lo.ln() + t * (hi.ln() - lo.ln())).exp() + } +} + +#[test] +fn information_gain_never_exceeds_the_entropy_of_the_outcome() { + let mut rng = Lcg(0x5eed_1234_abcd_ef01); + let ceiling = 2.0_f64.ln(); + let mut evaluated = 0usize; + let mut refused = 0usize; + + for i in 0..SAMPLES { + let mu_a = rng.in_range(-100.0, 100.0); + let mu_b = rng.in_range(-100.0, 100.0); + let sigma_a = rng.log_uniform(1e-4, 1e2); + let sigma_b = rng.log_uniform(1e-4, 1e2); + let beta = rng.log_uniform(1e-4, 1e1); + + let a = R::new(Gaussian::from_ms(mu_a, sigma_a), beta, ConstantDrift(0.0)); + let b = R::new(Gaussian::from_ms(mu_b, sigma_b), beta, ConstantDrift(0.0)); + let options = GameOptions { + p_draw: 0.0, + ..GameOptions::default() + }; + + match expected_information_gain(&[&[a], &[b]], &options) { + Ok(gain) => { + evaluated += 1; + assert!( + gain.is_finite(), + "sample {i}: non-finite gain {gain} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + assert!( + gain >= 0.0, + "sample {i}: negative gain {gain} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + assert!( + gain <= ceiling + 1e-9, + "sample {i}: gain {gain} exceeds ln 2 = {ceiling} \ + (mu {mu_a}, {mu_b}; sigma {sigma_a:e}, {sigma_b:e}; beta {beta:e})" + ); + } + // Refusing to answer is acceptable; answering wrongly is not. + Err(InferenceError::GridTooCoarse { .. }) => refused += 1, + Err(e) => panic!("sample {i}: unexpected error {e:?}"), + } + } + + // The sweep must actually exercise the function, not pass by refusing + // everything. + assert!( + evaluated * 2 > SAMPLES, + "only {evaluated} of {SAMPLES} samples were evaluated ({refused} refused); \ + the sweep is no longer testing anything" + ); + // And it must still reach the regime where the ceiling was violated — + // large sigma ratios, which is exactly where the grid now refuses. Without + // this the sweep could drift into only-easy inputs and stop being a guard. + assert!( + refused > 0, + "no sample reached the coarse-grid regime; the sweep no longer covers \ + the case that produced 3.24 nats" + ); +} + +/// The regime that produced 3.237828 nats, pinned exactly. +#[test] +fn the_known_ceiling_violation_no_longer_answers_wrongly() { + let a = R::new( + Gaussian::from_ms(9.577_887_112_129_012, 0.000_132_507_526_585_134_38), + 0.000_307_235_559_013_096_2, + ConstantDrift(0.0), + ); + let b = R::new( + Gaussian::from_ms(-14.114_932_828_525_696, 91.586_690_140_921_16), + 0.000_307_235_559_013_096_2, + ConstantDrift(0.0), + ); + let options = GameOptions { + p_draw: 0.0, + ..GameOptions::default() + }; + + match expected_information_gain(&[&[a], &[b]], &options) { + Ok(gain) => assert!( + gain <= 2.0_f64.ln() + 1e-9, + "returned {gain}, over the ln 2 ceiling" + ), + Err(InferenceError::GridTooCoarse { needed, max }) => { + assert!(needed > max, "needed {needed} should exceed max {max}"); + } + Err(e) => panic!("unexpected error {e:?}"), + } +} diff --git a/tests/quality.rs b/tests/quality.rs index 4b2570d..cbc15fe 100644 --- a/tests/quality.rs +++ b/tests/quality.rs @@ -164,3 +164,52 @@ fn quality_matches_the_reference_implementation() { let refs: Vec<&[Gaussian]> = five.iter().map(Vec::as_slice).collect(); assert!((quality(&refs, beta) - 0.040).abs() < 1e-9); } + +/// `quality()` used to compute `det(ata) / det(middle)` in linear space. Both +/// are products of `k - 1` diagonal entries, so they leave `f64`'s range long +/// before their ratio does — and the ratio is the only thing the answer needs. +/// +/// Measured before the fix: at the crate defaults 150 groups was correct, 200 +/// returned `0`, and 250 returned `NaN` where the truth is `9.51e-88`. With a +/// small beta it bit sooner — `sigma = beta = 1e-3` returned `NaN` at 60 groups +/// against a true `1.32e-9`, a value that is entirely ordinary. +/// +/// For `k` single-member groups with equal means the answer has a closed form, +/// `(beta / sqrt(beta^2 + sigma^2))^(k-1)`, so this checks against arithmetic +/// rather than against a recorded output. +#[test] +fn quality_matches_its_closed_form_past_the_overflow_point() { + for (sigma, beta) in [(25.0 / 3.0, 25.0 / 6.0), (1e-3, 1e-3), (50.0, 25.0 / 6.0)] { + let rating = vec![Gaussian::from_ms(25.0, sigma)]; + for k in [2usize, 50, 60, 150, 200, 250, 300] { + let groups: Vec<&[Gaussian]> = (0..k).map(|_| rating.as_slice()).collect(); + let got = quality(&groups, beta); + let expected = (beta / (beta * beta + sigma * sigma).sqrt()).powi(k as i32 - 1); + + assert!( + got.is_finite(), + "sigma {sigma}, beta {beta}, {k} groups: got {got}" + ); + // Subnormal results have no relative precision left to check. + if expected > f64::MIN_POSITIVE { + let rel = ((got - expected) / expected).abs(); + assert!( + rel < 1e-11, + "sigma {sigma}, beta {beta}, {k} groups: got {got:e}, \ + closed form {expected:e}, rel {rel:e}" + ); + } + } + } +} + +/// The overflow was in the intermediates, never in the answer: every value +/// above is an ordinary float. This pins the specific case that returned `NaN` +/// where the true answer is nine orders of magnitude inside the normal range. +#[test] +fn a_small_beta_does_not_overflow_at_sixty_groups() { + let rating = vec![Gaussian::from_ms(25.0, 1e-3)]; + let groups: Vec<&[Gaussian]> = (0..60).map(|_| rating.as_slice()).collect(); + let got = quality(&groups, 1e-3); + assert!((got - 1.317_089e-9).abs() / 1.317_089e-9 < 1e-6, "{got:e}"); +} diff --git a/tests/validation.rs b/tests/validation.rs index 9ec4686..a798fdf 100644 --- a/tests/validation.rs +++ b/tests/validation.rs @@ -255,3 +255,88 @@ mod builder_parameters { ); } } + +/// The constructors below `HistoryBuilder`, which 0.8.0's validation did not +/// reach. +/// +/// `sigma`, `beta` and `gamma` all enter inference only as squares, so a +/// negative value behaves as its absolute value and the sign vanishes without +/// comment. Measured before these guards: `from_ms(25.0, -8.33)` and +/// `Rating::new(_, -4.17, _)` returned results bit identical to their positive +/// counterparts, and `Rating::new(_, NaN, _)` reached `Game::ranked`, which +/// returned `Ok` carrying `Gaussian { pi: NaN, tau: NaN }`. +mod constructor_parameters { + use trueskill_tt::{ConstantDrift, Gaussian, History, InferenceError, Rating}; + + #[test] + #[should_panic(expected = "sigma must not be negative")] + fn a_negative_sigma_is_rejected_by_from_ms() { + let _ = Gaussian::from_ms(25.0, -8.33); + } + + /// NaN must pass, and that is deliberate: a broken fit produces a NaN + /// sigma and `converge` reports it as `NonFiniteResult`. Rejecting it here + /// would turn reporting into a panic inside inference. + #[test] + fn a_nan_sigma_passes_through_from_ms() { + let g = Gaussian::from_ms(25.0, f64::NAN); + assert!(g.sigma().is_nan() || g.pi().is_nan()); + } + + #[test] + #[should_panic(expected = "beta must be finite and non-negative")] + fn a_negative_beta_is_rejected_by_rating_new() { + let _ = Rating::::new(Gaussian::default(), -4.17, ConstantDrift(0.0)); + } + + #[test] + #[should_panic(expected = "beta must be finite and non-negative")] + fn a_nan_beta_is_rejected_by_rating_new() { + let _ = + Rating::::new(Gaussian::default(), f64::NAN, ConstantDrift(0.0)); + } + + #[test] + fn a_zero_beta_is_accepted_by_rating_new() { + let _ = Rating::::new(Gaussian::default(), 0.0, ConstantDrift(0.0)); + } + + /// `HistoryBuilder::drift` is generic and cannot inspect an arbitrary + /// `Drift`, so the check is on the variance each competitor actually + /// accumulates. That also covers a custom implementation. + #[test] + fn a_non_finite_drift_is_rejected_at_convergence() { + for gamma in [f64::NAN, f64::INFINITY] { + let mut h = History::builder() + .mu(25.0) + .sigma(25.0 / 3.0) + .beta(25.0 / 6.0) + .drift(ConstantDrift(gamma)) + .build(); + h.record_winner(&"a", &"b", 1).unwrap(); + h.record_winner(&"a", &"b", 5).unwrap(); + let err = h.converge().unwrap_err(); + assert!( + matches!( + err, + InferenceError::InvalidParameter { + name: "drift variance", + .. + } + ), + "gamma {gamma}: {err:?}" + ); + } + } + + /// An ordinary drift is untouched. + #[test] + fn an_ordinary_drift_still_converges() { + let mut h = History::builder() + .drift(ConstantDrift(25.0 / 300.0)) + .build(); + h.record_winner(&"a", &"b", 1).unwrap(); + h.record_winner(&"a", &"b", 5).unwrap(); + assert!(h.converge().unwrap().converged); + } +}