From 83bdb84152048950c1606e13fad2f3e9e271b5c0 Mon Sep 17 00:00:00 2001
From: Anders Olsson <anders.e.olsson@gmail.com>
Date: Wed, 9 Sep 2026 17:03:00 +0200
Subject: [PATCH] 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) <noreply@anthropic.com>
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<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
     ///
     /// Returns the matrix, each competitor's row at its latest appearance, and
     /// the row at each `(competitor, slice)` for time-addressed queries.
+    /// Drift below which two consecutive appearances are one latent variable.
+    ///
+    /// The rule used to be `drift <= 0.0` exactly, and everything above it got an
+    /// explicit `1.0 / drift` precision. That is a representation 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, `drift_scale =
+    /// 1e-10` returned a posterior variance **12 000x too small** — a 111x
+    /// overconfident interval — as `Ok`, and the band just above it returned a
+    /// misleading `JointUnavailable`.
+    ///
+    /// Solved exactly in high precision the same system is perfectly well
+    /// conditioned: it converges smoothly onto the collapsed value and is flat from
+    /// `1e-16` down to `1e-40`. So this is a representation problem, not a
+    /// conditioning one, and scaling cannot fix it — symmetric (Jacobi)
+    /// equilibration was measured **30x worse**, because the information is already
+    /// gone from the assembled matrix by the time a solver sees it.
+    ///
+    /// The threshold balances the two errors that trade off here. Ignoring a real
+    /// drift costs roughly `drift / V`; representing one costs roughly
+    /// `EPSILON * V / drift`, since `1 / drift` swamps the other precisions in the
+    /// row. They cross at `drift ~ V * sqrt(EPSILON)`, which is what this returns.
+    /// `V` is the competitor's own prior variance, so the threshold follows the
+    /// scale each competitor is actually measured on.
+    ///
+    /// Ordinary drift is far above this and is unaffected: the crate's default
+    /// `gamma = 25/300` accumulates `0.0069` per unit time against a threshold of
+    /// `1.0e-6` at the default prior.
+    fn collapse_threshold(prior_variance: f64) -> 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<Index, (usize, usize)> = HashMap::new();
         let mut at_slice: HashMap<(Index, usize), usize> = HashMap::new();
@@ -969,8 +1003,9 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
                     }
                     Some(&prev) => {
                         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<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
 
         let cholesky = crate::joint::Cholesky::factor(lambda, width).ok_or(
             InferenceError::JointUnavailable {
-                reason: "the precision matrix is not positive-definite, which means \
-                         a competitor has neither a proper prior nor any evidence",
+                reason: "the precision matrix is not positive-definite; the usual \
+                         cause is a competitor with neither a proper prior nor any \
+                         evidence, but an extreme prior or drift can also make the \
+                         assembled matrix indefinite in floating point",
             },
         )?;
 
diff --git a/tests/joint_handle.rs b/tests/joint_handle.rs
index e5d90df..9a72050 100644
--- a/tests/joint_handle.rs
+++ b/tests/joint_handle.rs
@@ -265,3 +265,75 @@ fn unseen_competitors_match_the_one_shot_path() {
     assert_eq!(one_shot.pi(), cached.pi());
     assert_eq!(one_shot.tau(), cached.tau());
 }
+
+/// A drift too small to represent must collapse, not corrupt the matrix.
+///
+/// The collapse rule used to fire only at `drift <= 0.0` exactly. Anything
+/// smaller-but-positive got an explicit `1.0 / drift` precision, and at
+/// `drift = 1e-16` that entry is `1e16` — so `1e16 + 0.28` rounds back to
+/// `1e16` and the prior and contrasts are annihilated in the stored `f64`.
+///
+/// Measured before the fix, at `drift_scale = 1e-10` this returned a variance
+/// **12 000x too small** (a 111x overconfident interval) as `Ok`, with a band
+/// just above it returning a misleading `JointUnavailable`.
+#[test]
+fn a_drift_too_small_to_represent_collapses_rather_than_corrupting() {
+    fn variance(scale: f64) -> f64 {
+        let mut h: History<i64, ConstantDrift, _, String> = 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}"
+    );
+}