fix!: propagate NaN through the convergence reduction

`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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011hcFjNDmHXZF8URGLku5zZ

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::*;