fix!: report an unresolvable prediction grid instead of clamping

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

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<T: Time, D: Drift<T>>(
 
     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;
         }

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}")
             }

src/history.rs

@@ -1539,7 +1539,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         Ok(Prediction::new(crate::predict::outcome_distribution(
             &performances,
             &self.margins(&sizes),
-        )))
+        )?))
     }
 
     /// Probability of one specific finishing order.
@@ -1579,11 +1579,7 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
         }
 
         let (performances, sizes) = self.performances(teams)?;
-        Ok(crate::predict::ranking_probability(
-            &performances,
-            &self.margins(&sizes),
-            ranks,
-        ))
+        crate::predict::ranking_probability(&performances, &self.margins(&sizes), ranks)
     }
 
     /// Run the full forward+backward convergence loop to a fixed point.

src/predict.rs

@@ -21,7 +21,7 @@
 //! would have made every `predict_*` call return a slightly different number,
 //! which is not a property a rating library should have.
 
-use crate::{Gaussian, quadrature};
+use crate::{Gaussian, InferenceError, quadrature};
 
 /// Teams beyond this count make the outcome enumeration impractical.
 ///
@@ -52,6 +52,10 @@ const WIN_TOLERANCE: f64 = 1e-8;
 /// point where refining stops helping.
 const MIN_GRID_POINTS: usize = 8_192;
 const MAX_GRID_POINTS: usize = 262_144;
+/// Nodes requested across the narrowest feature the recursion must resolve.
+const NODES_PER_FEATURE: f64 = 12.0;
+/// Nodes below which the trapezoid rule stops resolving that feature at all.
+const MIN_NODES_PER_FEATURE: f64 = 4.0;
 
 /// How many standard deviations of support the grid and integrals cover.
 ///
@@ -152,7 +156,7 @@ pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec<f64
 /// Resolution is set by the *smallest* feature in play — the narrowest sigma,
 /// or a draw margin narrower still — because that is what the recursion has to
 /// resolve. A grid sized off the widest team would step over the narrow one.
-fn grid_shape(perf: &[Gaussian], margins: &Margins) -> (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<Self, InferenceError> {
+        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<usize>, Vec<bool>)> {
 ///
 /// 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<u32>, f64)> {
+pub(crate) fn outcome_distribution(
+    perf: &[Gaussian],
+    margins: &Margins,
+) -> Result<Vec<(Vec<u32>, f64)>, InferenceError> {
     let n = perf.len();
-    let sampled = Sampled::new(perf, margins);
+    let sampled = Sampled::new(perf, margins)?;
 
     let mut aggregated: Vec<(Vec<u32>, 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<usize>]) -> Vec<(Vec<usize>, Vec<bool>)> {
 /// 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<f64, InferenceError> {
     let n = perf.len();
-    let sampled = Sampled::new(perf, margins);
+    let sampled = Sampled::new(perf, margins)?;
 
     let mut distinct: Vec<u32> = 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<f64> = [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();