From c65373f47611846d95a9ab150804de7426807580 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 16:59:56 +0200 Subject: [PATCH 01/10] fix!: propagate NaN through the convergence reduction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `tuple_max` compared with a plain `>`, which is false against NaN, so a NaN accumulator was replaced by the next finite delta. The fold runs over `TimeSlice::posteriors()`, a HashMap, so whether a NaN survived to `step` depended on per-process hash order. Measured before, four competitors in one slice with one pathological pair, same binary and input, 30 separate processes: 16 Ok converged=true, iterations=1, a = Gaussian { pi: NaN, tau: NaN } 14 Err NonFiniteResult After: 30/30 Err. A coin flip on whether a NaN fit was reported as an error or as a successful, converged fit — inside the guard whose entire purpose is "NaN is never convergence". `f64::max` would not have fixed it. It also ignores NaN by design, which is the same defect wearing a standard-library name, and a test pins that we do not use it. `Gaussian::delta` had to be fixed FIRST, and that ordering is the whole subtlety. Two identical improper messages produced `(0.0, NaN)` — not from `mu()`, which is guarded and returns 0.0, but from `inf - inf` in the sigma component. That NaN is reachable in ordinary healthy inference: once a pairing is more than about nine cavity-sigma apart the truncation is a no-op and the chain compares one identity message against another. Propagating NaN without fixing `delta` would therefore have turned correct fits into NonFiniteResult errors. `delta` now answers the identical-message case in natural space before touching the accessors. My first version of the `delta` test asserted `mu()` was NaN. It is not; the accessor guards `pi <= 0.0`. The test caught my own wrong premise, and the doc comment is corrected to match. BREAKING CHANGE: a fit that produced NaN in a non-final reduction position previously returned `Ok` with `converged: true` and a NaN posterior; it now returns `Err(NonFiniteResult)`. That was always the documented intent. Closes #58 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/gaussian.rs | 54 +++++++++++++++++++++++++++++ src/lib.rs | 69 +++++++++++++++++++++++++++++++++++-- tests/non_finite_results.rs | 48 ++++++++++++++++++++++++++ 3 files changed, 169 insertions(+), 2 deletions(-) diff --git a/src/gaussian.rs b/src/gaussian.rs index 23cc261..c1d76be 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -120,7 +120,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 +274,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/lib.rs b/src/lib.rs index a337b45..10d8396 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -569,13 +569,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 } @@ -724,6 +745,50 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { #[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::*; diff --git a/tests/non_finite_results.rs b/tests/non_finite_results.rs index 8928431..2434172 100644 --- a/tests/non_finite_results.rs +++ b/tests/non_finite_results.rs @@ -115,3 +115,51 @@ 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:?}" + ); +} From 83bdb84152048950c1606e13fad2f3e9e271b5c0 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:03:00 +0200 Subject: [PATCH 02/10] fix!: collapse a drift too small to represent, on a relative threshold MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `time_expanded_joint` collapsed consecutive appearances only at `drift <= 0.0` exactly. Anything smaller-but-positive got an explicit `1.0 / drift` precision, which the matrix cannot hold: at `drift = 1e-16` the entry is `1e16`, and `1e16 + 0.28` rounds back to `1e16`, so the prior and the event contrasts are annihilated in the stored f64 before the factorisation ever runs. Measured, 8 competitors over 15 slices: drift_scale before after 1e-6 1.3e-3 relative error exact 1e-7..1e-9 Err(JointUnavailable) exact 1e-10 12 200x TOO SMALL, as Ok exact At 1e-10 the caller was handed sigma = 0.0055 where the truth is 0.6108 — a 111x overconfident interval, returned as a success. This is representation, not conditioning. Solved in 200-digit precision the same system converges smoothly onto the collapsed value and is flat from 1e-16 to 1e-40, so the quantity is perfectly well conditioned. That also rules out the obvious fix: symmetric (Jacobi) equilibration measured 30x WORSE, because the information is gone from the assembled matrix before any solver sees it. The fix has to be at assembly. The threshold balances the two errors that trade off. Ignoring a real drift costs about `drift / V`; representing one costs about `EPSILON * V / drift`. They cross at `V * sqrt(EPSILON)`, scaled to each competitor's own prior variance. Ordinary drift is far above it and unaffected — the default gamma accumulates 0.0069 per unit time against a threshold of 1.0e-6 — and the test asserts both halves: everything below the threshold reaches the collapsed answer bit-identically, and a drift of 1e-2 still moves it, so the test cannot pass by collapsing everything. Also corrects the `JointUnavailable` message, which asserted "a competitor has neither a proper prior nor any evidence" for a fixture where every competitor had both. BREAKING CHANGE: a drift variance below `prior_variance * sqrt(EPSILON)` now collapses two appearances into one latent variable. Affected fits previously returned a badly wrong variance or an error. Closes #57 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/history.rs | 45 ++++++++++++++++++++++++--- tests/joint_handle.rs | 72 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 113 insertions(+), 4 deletions(-) diff --git a/src/history.rs b/src/history.rs index bd86277..6aaeb9e 100644 --- a/src/history.rs +++ b/src/history.rs @@ -947,6 +947,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 +1003,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; @@ -1278,8 +1313,10 @@ impl, O: Observer, K: Eq + Hash + Clone> History 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}" + ); +} From bbc7705c75b1736c9a653ab40c2a0a8348cb6e3f Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:21:00 +0200 Subject: [PATCH 03/10] fix!: report an unresolvable prediction grid instead of clamping MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `grid_shape` asked for 12 nodes across the narrowest feature and then clamped to MAX_GRID_POINTS with no detection that the request was not met. Past `step/sigma ~ 1.7` the trapezoid rule stops resolving the density, and the result is unbounded: sigma_a step/sig_a P(a first) exact total 2.0e-3 0.86 0.515953 0.515953 1.000000 1.0e-3 1.72 0.517185 0.515953 1.002388 1.0e-4 17.17 2.791336 0.515953 5.410065 A probability of 2.79. Reachable through `predict_outcome` with a pinned reference competitor — a documented pattern — where `predict_outcome` and `predict_win_probabilities` disagreed 44x and `predict_outcome` was the wrong one. There is no useful answer on the far side of that cliff, so this reports `GridTooCoarse` rather than guessing, and the message points at `predict_win_probabilities`, which answers the same matchup through adaptive quadrature and is accurate there to 1e-13. The floor is 4 nodes per feature rather than the 12 requested, because the request carries margin: measured accurate to 2.2e-12 at 1.4 nodes per sigma and wrong by 1.2e-3 at 0.7. This also fixes the `ln k` ceiling violation. `expected_information_gain` weights `probability * divergence`, so probabilities of 3.97 and 2.62 made it return 3.237828 nats against `ln 2 = 0.693147` — 4.67x over. The crate's docs call that ceiling its sharpest test and record a prototype once returning 4.77 nats; it was live again by a different route. The new sweep then caught a second, independent defect: `kl_divergence` returned NEGATIVE values, worst -5.55e-17, exactly one ULP of its `- 1.0`. Rewritten as `0.5*(u - ln1p(u)) + gap^2/(2*var_p)` with `u = var_q/var_p - 1`, so both terms are non-negative by construction. It is also more accurate where it matters: at `u = 1e-9` the old form returned 0.0 where the true value is 2.5e-19, and well-conditioned cases are unchanged. tests/prediction_bounds.rs sweeps rather than spot-checks, because a single fixture cannot defend a bound like this — the previous check passed throughout. It asserts the sweep still reaches the coarse-grid regime, so it cannot quietly stop testing the case it was written for. BREAKING CHANGE: `predict_outcome`, `predict_ranking` and `expected_information_gain` return `GridTooCoarse` for matchups whose performance sigmas are too far apart to integrate on one grid. They previously returned wrong answers, including probabilities above 1. Closes #55, closes #56 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/acquisition.rs | 35 +++++++++- src/error.rs | 31 +++++++++ src/history.rs | 8 +-- src/predict.rs | 83 ++++++++++++++-------- tests/prediction_bounds.rs | 136 +++++++++++++++++++++++++++++++++++++ 5 files changed, 258 insertions(+), 35 deletions(-) create mode 100644 tests/prediction_bounds.rs 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/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/history.rs b/src/history.rs index 6aaeb9e..dafc48b 100644 --- a/src/history.rs +++ b/src/history.rs @@ -1539,7 +1539,7 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History 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/tests/prediction_bounds.rs b/tests/prediction_bounds.rs new file mode 100644 index 0000000..06e5b19 --- /dev/null +++ b/tests/prediction_bounds.rs @@ -0,0 +1,136 @@ +//! 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; + +/// 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..2_000 { + 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 > 1_000, + "only {evaluated} of 2000 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:?}"), + } +} From 31cf0998b073fe983c905d3d1f180ea72d2c6fc8 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:27:09 +0200 Subject: [PATCH 04/10] test: scale the ceiling sweep by build profile MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Each sample runs a full inference pass per outcome, and that is about 19x faster in release: 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 sample count pays the slow price three times and the fast one once — exactly backwards. Scaling by `cfg!(debug_assertions)` puts the search where it is cheap: debug 1 000 samples 11.7s release 50 000 samples 31.6s Across the whole `just test` that is 67s against 70s before, for 25x the samples. The debug run proves the sweep compiles and holds; the release run is the one that actually searches. Not moving the suite to release-only, which was the alternative considered. `debug_assert!` is compiled out in release, and this crate documents that as load-bearing — several defects have hidden there — so dropping the debug runs would trade one class of coverage for another rather than adding any. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- tests/prediction_bounds.rs | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/tests/prediction_bounds.rs b/tests/prediction_bounds.rs index 06e5b19..eedde48 100644 --- a/tests/prediction_bounds.rs +++ b/tests/prediction_bounds.rs @@ -17,6 +17,24 @@ use trueskill_tt::{ 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); @@ -49,7 +67,7 @@ fn information_gain_never_exceeds_the_entropy_of_the_outcome() { let mut evaluated = 0usize; let mut refused = 0usize; - for i in 0..2_000 { + 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); @@ -91,8 +109,8 @@ fn information_gain_never_exceeds_the_entropy_of_the_outcome() { // The sweep must actually exercise the function, not pass by refusing // everything. assert!( - evaluated > 1_000, - "only {evaluated} of 2000 samples were evaluated ({refused} refused); \ + 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 — From f1219036b337a6992567573d49229eb67d58058e Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:33:00 +0200 Subject: [PATCH 05/10] fix: take quality's determinant ratio in log space MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `quality()` computed `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 at the crate defaults: 150 groups correct at 8.45e-53, 200 returned 0, 250 returned NaN where the truth is 9.51e-88. With a small beta it bit far sooner: at sigma = beta = 1e-3, 60 groups returned NaN against a true 1.32e-9 — a value nine orders of magnitude inside the normal range. Neither `quality()` nor `History::predict_quality` caps the group count, unlike `predict_outcome`, so those are supported calls. `Lu::ln_abs_determinant` accumulates `ln|diagonal|` instead of multiplying, and the call site becomes `exp(e_arg + 0.5 * ln_ratio)`. Verified against the closed form `(beta / sqrt(beta^2 + sigma^2))^(k-1)` rather than against recorded output, across three parameter sets and group counts to 300: every case now agrees to 1e-11 or better, including 9.88e-324 at 300 groups, which is subnormal. Also documents the remaining panic: every rating at zero sigma with a zero beta makes `middle` singular and `inverse()` panics. Documented rather than converted — nothing is uncertain there, so there is no distribution to take the quality of. Closes #59 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/lib.rs | 21 +++++++++++++++++++-- src/matrix.rs | 41 ++++++++++++++++++++++++++++++++++++++++ tests/quality.rs | 49 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 10d8396..88deb5a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -673,6 +673,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!( @@ -738,9 +745,19 @@ 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)] 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/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}"); +} From 6139061740adc1edc35ed9dce267e34a1a3aeb8f Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:42:16 +0200 Subject: [PATCH 06/10] fix: keep the truncated variance representable in the far tail MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `v_w` returned `w` and let `trunc` form `1 - w`. `w` tends to 1 out in the tail, so that subtraction lost about log10(alpha^2) digits — and the quantity it was destroying is perfectly representable. Two separate cancellations, fixed separately. The non-tie half: `half_line_truncation` now returns `1 - w` computed symbolically rather than as `1 - v*gap`. With `alpha*gap = 1 - inv^2*b` the leading ones cancel on paper instead of in floating point. Measured against the exact truncated variance: alpha before after 1e6 8.9e-5 rel 0.0 rel (exact) 1e8 returns 0.0 0.0 rel (exact) At 1e8 the old form gave `sigma_trunc = 0`, and `from_ms(mu, 0.0)` is a point mass whose `mu()` is inf/inf = NaN. `beta(1e-8).sigma(1e-8)` with priors 1000 apart went from Err + NaN skills to a finite fit. The tie half is a different subtraction — `w = v^2 - u`, where both grow as alpha^2 while their difference stays O(1). The existing escape hatch could not cover it: it keys on `alpha * width >= HALF_LINE_WINDOW`, how many window-widths from the mean the window sits, and a NARROW window fails that however deep it is. Measured at alpha 1e6 with a 1e-6 window it kept four digits and returned `1 - w = -2.4e-4` where the truth is +2.8e-13. 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. Over a narrow window the density is a truncated exponential in `s = (x - alpha)/width`, whose mean and variance are closed forms, so `v = alpha + width*m(t)` and `1 - w = width^2 * V(t)` with no large subtraction at all. Validated against high-precision quadrature: v exact to 4e-10, `1 - w` to 4e-10 across the region it is used in. The crossover is on `alpha / width` rather than on either alone, because that ratio is what says how many digits the subtraction has left — and the approximation is most accurate exactly where the subtraction is worst, since both improve as the window narrows. Defaults are bit-identical (pi 0.02398318151216503 before and after). Tests: the three reproductions from the issue, the narrow-window form against pinned quadrature values, and a continuity sweep across all three tie branches — a misplaced crossover is the real risk here, and a jump at a boundary is visible even without pinning absolute values. Closes #60 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/lib.rs | 181 ++++++++++++++++++++++++++++++++++-- tests/non_finite_results.rs | 54 ++++++++++- 2 files changed, 224 insertions(+), 11 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 88deb5a..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) } @@ -762,6 +849,80 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 { #[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 diff --git a/tests/non_finite_results.rs b/tests/non_finite_results.rs index 2434172..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, @@ -163,3 +165,53 @@ fn a_nan_competitor_is_not_masked_by_a_healthy_one() { "{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:?}" + ); + } +} From ab23476aaf51ad8ea4878d462580ad12c6fd6b7b Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:48:41 +0200 Subject: [PATCH 07/10] fix!: validate the constructors below HistoryBuilder MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 0.8.0 closed the sign-absorption defect at `HistoryBuilder::mu/sigma/beta` and at both ingestion paths. It was still open one layer down, in the constructors those paths call. Measured, all bit identical to their positive counterparts: Gaussian::from_ms(25.0, -8.33) == from_ms(25.0, +8.33) Rating::new(_, -4.17, _) == Rating::new(_, +4.17, _) ConstantDrift(-0.0833) == ConstantDrift(+0.0833) sigma, beta and gamma enter only as squares, so the sign vanished without comment. Worst of the set: `Rating::new(_, NaN, _)` reached `Game::ranked` which returned **Ok** carrying `Gaussian { pi: NaN, tau: NaN }` — no `converge` on that path to catch it. `from_ms` and `Rating::new` now reject. `ConstantDrift` cannot: the field is public and positional, so there is no constructor to intercept, and sealing it would break every `ConstantDrift(x)` for a case whose resulting model is perfectly valid. Documented instead. Its non-finite half IS rejected — `converge` validates the drift variance each competitor accumulates, which also covers a custom `Drift` impl. Two things the tests caught that I had wrong: NaN sigma must PASS `from_ms`. My first version rejected it, and two existing tests went red immediately: a broken fit legitimately produces a NaN sigma from `sqrt` of a negative truncated variance, and the design is to propagate that to `NonFiniteResult`. Rejecting it turned the reporting path into a panic inside inference. Written as `sigma >= 0.0 || sigma.is_nan()` so the intent is explicit rather than hidden in a negated comparison. Very small sigma is also not rejected, and that is deliberate: `approx` produces small truncated sigmas legitimately. `pi = 1/sigma^2` leaves f64's range below ~1.5e-154 and `tau = mu*pi` overflows sooner, 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 existing point-mass representation. Documented, including that such a Gaussian is not equal to itself and can make two identical declarations report as conflicting. BREAKING CHANGE: `Gaussian::from_ms` panics on a negative sigma, and `Rating::new` panics unless beta is finite and non-negative. `converge` returns `InvalidParameter` for a non-finite drift variance. Closes #61 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/drift.rs | 17 +++++++++ src/gaussian.rs | 36 +++++++++++++++++++ src/history.rs | 24 +++++++++++++ src/rating.rs | 17 +++++++++ tests/validation.rs | 85 +++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 179 insertions(+) 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/gaussian.rs b/src/gaussian.rs index c1d76be..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 { diff --git a/src/history.rs b/src/history.rs index dafc48b..fa1bf70 100644 --- a/src/history.rs +++ b/src/history.rs @@ -1644,6 +1644,30 @@ impl, 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, 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/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); + } +} From 305f8229644cc973a2ddff86c3ad62216fb59794 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 17:56:23 +0200 Subject: [PATCH 08/10] fix: route the last three transcendentals through libm, and enforce it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CLAUDE.md requires transcendentals to go through libm rather than std, because std delegates to the system math library and the two disagree by one ULP often enough to change an iteration count in a fixed point. Three production sites did not: factor/margin.rs:84 cavity.sigma().hypot(sigma) 12.136% of 1e6 inputs factor/margin.rs:92 f64::MIN_POSITIVE.ln() 3.437% factor/trunc.rs:98 f64::MIN_POSITIVE.ln() same `hypot` is the material one: it is on the path of every scored event, and its divergence rate is HIGHER than the 9.7% the rule cites for `exp` as its own justification. The two `ln` calls happen to agree bit-for-bit on this host, which is exactly the platform dependence the rule exists to remove. The `hypot` choice itself was right and stays — the comment above it explains why, and it is measured: naive sqrt(a^2 + b^2) overflows to inf at 1e200 and flushes to zero at 1e-200 where hypot does neither. Only the implementation moves. tests/libm_rule.rs enforces it. The rule was stated plainly in CLAUDE.md and still violated three times, so prose is evidently not sufficient. The test strips `#[cfg(test)]` items by brace matching, plus comments and string literals so prose is not mistaken for a call, then scans for std method spellings. `sqrt` is exempt: IEEE 754 specifies it, so std and libm cannot disagree. Confirmed non-vacuous by reintroducing the `hypot` violation and watching it fail with the offending line, then pass again on restore. Two further tests pin the stripper itself, since a stripper that removed everything would make the guard pass on anything. Closes #63 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/factor/margin.rs | 4 +- src/factor/trunc.rs | 4 +- tests/libm_rule.rs | 175 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 179 insertions(+), 4 deletions(-) create mode 100644 tests/libm_rule.rs 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/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}"); +} From 7aa7fb62ddf7ef8214041baa88698009753bbac7 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 18:02:42 +0200 Subject: [PATCH 09/10] fix: make posterior_of reproducible across processes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `ResolvedTerms::unseen` was a `HashMap` and three float reductions iterated it. Addition is not associative and Rust seeds its default hasher per process, so `posterior_of` returned 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 about 7 ULP. A `BTreeMap` fixes it by construction. 40/40 identical after, 24/16 before. The cross-batch conflict scan had the same cause with a different symptom. It returns on the FIRST conflict, so hash order decided WHICH competitor the error blamed — 15 different competitors named across 40 runs on identical input. The error fired every time; only its content was a lottery, which sends a reader after the wrong key. Now scanned in sorted order. Magnitude was 1-7 ULP throughout, so no decision changes. The cost was reproducibility: a golden test over these would flake at a low rate, which is the worst kind of CI failure to diagnose. tests/cross_process_determinism.rs re-executes the test binary and compares bits, because an in-process test CANNOT see this — every sample in one process shares one hasher seed. That is not hypothetical: tests/determinism.rs compares four thread counts inside one process and passed throughout while this was live. Tuning that fixture took a measurement. Coefficients spread over nine decades detected the bug in roughly one run in forty, because the small terms fall below the running total's ULP and are absorbed whatever the order. Comparable magnitudes keep every term able to change the last bits: 5 of 5 attempts detected it, with 3 to 38 of 40 runs differing. Verified non-vacuous by reverting the BTreeMap and watching it fail. Closes #62 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- src/history.rs | 30 +++++- tests/cross_process_determinism.rs | 157 +++++++++++++++++++++++++++++ 2 files changed, 183 insertions(+), 4 deletions(-) create mode 100644 tests/cross_process_determinism.rs diff --git a/src/history.rs b/src/history.rs index fa1bf70..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, } @@ -1067,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() { @@ -1858,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/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] + ); +} From c69a397d80a52843fbc12b974dd562b1688181b7 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 18:08:06 +0200 Subject: [PATCH 10/10] test: make the determinism test exercise the parallel sweep MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit It 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 a colour group is a subset of ONE slice's events. The fixture built 20 slices of 10, so the branch was unreachable — the test named the parallel path and ran the sequential one. It also compared one competitor's curve out of forty, and never compared log_evidence, final_step or iterations. The new fixture reaches the branch by construction: within a slice every event uses a disjoint competitor pair, so greedy colouring puts all 96 in colour 0. Competitors recur across slices, so the fit keeps temporal coupling and drift rather than degenerating into independent duels. Verified by instrumenting `sweep_color_groups`: 872 sweeps, one colour group of 96 each, parallel branch taken all 872 times. Worth recording how that verification went, because I nearly drew the opposite conclusion. My first two instrumented runs printed nothing and I read that as "the branch is still unreachable" — but `cargo test` captures stderr without `--nocapture`, so the probe was invisible, not absent. An instrument that cannot report is indistinguishable from a negative result. Now compares every competitor's curve plus log_evidence, final_step and iterations, and asserts the curve count so it cannot silently go back to measuring almost nothing. A companion test pins EVENTS_PER_SLICE against the threshold, so shrinking the fixture fails loudly rather than quietly returning the suite to the sequential path. Cross-process coverage is separate, in tests/cross_process_determinism.rs (#62) — an in-process test cannot see hasher-order effects at all, since every sample shares one seed. Closes #64 Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ --- tests/determinism.rs | 226 +++++++++++++++++++++++++++++++------------ 1 file changed, 163 insertions(+), 63 deletions(-) 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); +}