From c65373f47611846d95a9ab150804de7426807580 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Wed, 9 Sep 2026 16:59:56 +0200 Subject: [PATCH] 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:?}" + ); +}