trueskill-tt/src/quadrature.rs

//! Deterministic numerical integration for the prediction paths.
//!
//! Prediction asks two questions that have no closed form beyond two teams:
//! "who finishes first" and "how likely is this exact finishing order". Both
//! reduce to integrals over a single performance variable, so neither needs a
//! sampler — and that matters, because a Monte Carlo predictor would make
//! `predict_*` non-reproducible and would answer a slightly different question
//! on every call.
//!
//! Two routines live here:
//!
//! - [`integrate`], adaptive Gauss-Kronrod G7-K15, for the first-place
//!   marginals. It carries its own error estimate, so it can refine where the
//!   integrand actually bends instead of guessing a node count up front.
//! - [`Grid`], a uniform grid with trapezoid prefix sums, for the ranking
//!   chain recursion, where each level needs the *running* integral of the
//!   level below at arbitrary points rather than one definite integral.
//!
//! Fixed-node Gauss-Hermite is the obvious tool for the first of these and is
//! a trap: the integrand is a product of normal CDFs, and when one team's
//! sigma is much smaller than the integrating team's, that product turns into
//! a near-step function narrower than the node spacing. The nodes step over
//! it and the result is wrong by ~1e-2 while still looking like a probability.
//! Adaptive refinement is what makes the small-sigma case safe.

/// Kronrod 15-point abscissae, non-negative half, descending.
const XGK: [f64; 8] = [
    0.991_455_371_120_813,
    0.949_107_912_342_759,
    0.864_864_423_359_769,
    0.741_531_185_599_394,
    0.586_087_235_467_691,
    0.405_845_151_377_397,
    0.207_784_955_007_898,
    0.0,
];

/// Kronrod 15-point weights, matching [`XGK`].
const WGK: [f64; 8] = [
    0.022_935_322_010_529,
    0.063_092_092_629_979,
    0.104_790_010_322_250,
    0.140_653_259_715_525,
    0.169_004_726_639_267,
    0.190_350_578_064_785,
    0.204_432_940_075_298,
    0.209_482_141_084_728,
];

/// Gauss 7-point weights, applying to the odd-indexed [`XGK`] entries.
const WG: [f64; 4] = [
    0.129_484_966_168_870,
    0.279_705_391_489_277,
    0.381_830_050_505_119,
    0.417_959_183_673_469,
];

/// Panels are bisected worst-first; this bounds the work on a pathological
/// integrand rather than letting it spin.
const MAX_SUBDIVISIONS: usize = 200;

/// One G7-K15 panel over `[a, b]`: `(integral, absolute error estimate)`.
///
/// The error estimate is the gap between the embedded 7-point Gauss rule and
/// the 15-point Kronrod extension. It is the only reason this is preferable
/// to a fixed rule: it tells the caller *where* the integrand is hard.
fn gk15<F: Fn(f64) -> f64>(f: &F, a: f64, b: f64) -> (f64, f64) {
    let centre = 0.5 * (a + b);
    let half = 0.5 * (b - a);

    let mut kronrod = 0.0;
    let mut gauss = 0.0;

    for i in 0..8 {
        let offset = XGK[i] * half;
        // XGK[7] is the centre node and must not be counted twice.
        let sum = if i == 7 {
            f(centre)
        } else {
            f(centre - offset) + f(centre + offset)
        };
        kronrod += WGK[i] * sum;
        if i % 2 == 1 {
            gauss += WG[i / 2] * sum;
        }
    }

    (kronrod * half, ((kronrod - gauss) * half).abs())
}

/// Adaptively integrate `f` over `[a, b]` to relative tolerance `tol`.
///
/// `seeds` are interior points where the integrand is known to bend sharply —
/// for a product of normal CDFs, each rival's transition centre. Splitting
/// there up front costs nothing and saves the adaptive loop from having to
/// discover a step by bisection.
///
/// Returns the integral. The error estimate is consumed internally rather
/// than returned: callers here integrate probability densities, where the
/// meaningful check is the sum-to-one identity over a whole outcome space,
/// not a per-integral residual.
pub(crate) fn integrate<F: Fn(f64) -> f64>(f: F, a: f64, b: f64, seeds: &[f64], tol: f64) -> f64 {
    // Explicit rather than `!(b > a)`: a NaN bound must fall through to zero
    // rather than being read as a valid ordering.
    if a.partial_cmp(&b) != Some(std::cmp::Ordering::Less) {
        return 0.0;
    }

    let mut edges: Vec<f64> = Vec::with_capacity(seeds.len() + 2);
    edges.push(a);
    edges.push(b);
    for &s in seeds {
        if s > a && s < b {
            edges.push(s);
        }
    }
    edges.sort_by(|p, q| p.partial_cmp(q).expect("integration bounds are finite"));
    edges.dedup();

    // (lo, hi, integral, error)
    let mut panels: Vec<(f64, f64, f64, f64)> = edges
        .windows(2)
        .map(|w| {
            let (v, e) = gk15(&f, w[0], w[1]);
            (w[0], w[1], v, e)
        })
        .collect();

    for _ in 0..MAX_SUBDIVISIONS {
        let total: f64 = panels.iter().map(|p| p.2).sum();
        let error: f64 = panels.iter().map(|p| p.3).sum();

        // Absolute floor as well as relative: these integrands are
        // probabilities, so an absolute 1e-15 is already past the useful
        // precision of the underlying `cdf`.
        if error <= tol * total.abs().max(1e-12) || error < 1e-15 {
            break;
        }

        let worst = panels
            .iter()
            .enumerate()
            .fold((0usize, f64::NEG_INFINITY), |(bi, be), (i, p)| {
                if p.3 > be { (i, p.3) } else { (bi, be) }
            })
            .0;

        let (lo, hi, _, _) = panels[worst];
        let mid = 0.5 * (lo + hi);
        // Bisection has hit the floating-point floor; refining further would
        // loop without reducing the error.
        if !(mid > lo && mid < hi) {
            break;
        }

        let (v1, e1) = gk15(&f, lo, mid);
        let (v2, e2) = gk15(&f, mid, hi);
        panels[worst] = (lo, mid, v1, e1);
        panels.push((mid, hi, v2, e2));
    }

    panels.iter().map(|p| p.2).sum()
}

/// A uniform grid carrying trapezoid prefix sums of one integrand.
///
/// The ranking recursion needs, at every level, the running integral of the
/// level below evaluated at arbitrary points — a cumulative integral, not a
/// definite one. Prefix sums give that in O(1) per query after an O(G) build,
/// which is what keeps a full ranking probability linear in the team count.
pub(crate) struct Grid {
    lo: f64,
    step: f64,
    /// Integrand sampled at each node.
    values: Vec<f64>,
    /// `prefix[i]` is the integral from `lo` to node `i`.
    prefix: Vec<f64>,
}

impl Grid {
    /// Build directly from already-sampled values.
    ///
    /// The ranking recursion evaluates every level on the same nodes, so the
    /// per-team densities are sampled once and reused; re-evaluating `exp`
    /// per level would dominate the cost.
    pub(crate) fn from_values(lo: f64, step: f64, values: Vec<f64>) -> Self {
        let mut prefix = vec![0.0; values.len()];
        for i in 1..values.len() {
            prefix[i] = prefix[i - 1] + 0.5 * step * (values[i - 1] + values[i]);
        }
        Self {
            lo,
            step,
            values,
            prefix,
        }
    }

    /// Integral from the grid's lower bound up to `x`.
    ///
    /// Clamped at both ends: the caller sizes the grid to cover the whole
    /// support, so a query outside it is asking for a tail that is zero (below)
    /// or the whole mass (above).
    pub(crate) fn integral_to(&self, x: f64) -> f64 {
        let last = self.values.len() - 1;
        if x <= self.lo {
            return 0.0;
        }
        if x >= self.lo + last as f64 * self.step {
            return self.prefix[last];
        }

        let scaled = (x - self.lo) / self.step;
        let i = scaled.floor() as usize;
        let frac = scaled - i as f64;

        // Whole cells, plus the trapezoid over the partial cell. The integrand
        // is linear within a cell under the trapezoid rule, so the partial
        // piece is exact with respect to that same approximation.
        self.prefix[i]
            + frac
                * self.step
                * (self.values[i] + 0.5 * frac * (self.values[i + 1] - self.values[i]))
    }

    /// Integral over `[from, to]`.
    pub(crate) fn integral_between(&self, from: f64, to: f64) -> f64 {
        (self.integral_to(to) - self.integral_to(from)).max(0.0)
    }

    /// Total integral over the whole grid.
    pub(crate) fn total(&self) -> f64 {
        self.prefix[self.values.len() - 1]
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    const TOL: f64 = 1e-10;

    /// Sample `f` over `[lo, hi]` at `points` nodes.
    fn sample<F: FnMut(f64) -> f64>(lo: f64, hi: f64, points: usize, mut f: F) -> Grid {
        let step = (hi - lo) / (points - 1) as f64;
        Grid::from_values(
            lo,
            step,
            (0..points).map(|i| f(lo + i as f64 * step)).collect(),
        )
    }

    #[test]
    fn integrates_a_polynomial_exactly() {
        // G7-K15 is exact for polynomials well past cubic, so a single panel
        // should already be at round-off.
        let v = integrate(|x| 3.0 * x * x + 2.0 * x + 1.0, 0.0, 2.0, &[], TOL);
        assert!((v - 14.0).abs() < 1e-12, "got {v}");
    }

    #[test]
    fn integrates_a_gaussian_density_to_one() {
        let f = |x: f64| (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt();
        let v = integrate(f, -10.0, 10.0, &[], TOL);
        assert!((v - 1.0).abs() < 1e-12, "got {v}");
    }

    #[test]
    fn resolves_a_step_far_narrower_than_the_initial_panel() {
        // The failure mode that rules out fixed-node quadrature: a transition
        // 1e-4 wide inside a range of 20. A fixed rule steps over it.
        let f = |x: f64| if x < 0.5 { 0.0 } else { 1.0 };
        let v = integrate(f, -10.0, 10.0, &[0.5], TOL);
        assert!((v - 9.5).abs() < 1e-6, "got {v}");
    }

    #[test]
    fn seeds_do_not_change_the_value_of_a_smooth_integrand() {
        let f = |x: f64| (-0.5 * x * x).exp();
        let plain = integrate(f, -8.0, 8.0, &[], TOL);
        let seeded = integrate(f, -8.0, 8.0, &[-3.0, 0.25, 5.5], TOL);
        assert!((plain - seeded).abs() < 1e-12, "{plain} vs {seeded}");
    }

    #[test]
    fn empty_or_inverted_range_integrates_to_zero() {
        assert_eq!(integrate(|_| 1.0, 1.0, 1.0, &[], TOL), 0.0);
        assert_eq!(integrate(|_| 1.0, 2.0, 1.0, &[], TOL), 0.0);
    }

    #[test]
    fn grid_prefix_matches_a_known_cumulative_integral() {
        // f(x) = x over [0, 4]; integral to x is x^2/2.
        let g = sample(0.0, 4.0, 4001, |x| x);
        for probe in [0.0, 0.5, 1.0, 2.5, 3.75, 4.0] {
            let want = probe * probe / 2.0;
            let got = g.integral_to(probe);
            assert!(
                (got - want).abs() < 1e-9,
                "at {probe}: got {got}, want {want}"
            );
        }
        assert!((g.total() - 8.0).abs() < 1e-9);
    }

    #[test]
    fn grid_between_is_the_difference_of_two_prefixes() {
        let g = sample(-5.0, 5.0, 8001, |x| (-0.5 * x * x).exp());
        let whole = g.integral_between(-5.0, 5.0);
        let split = g.integral_between(-5.0, 0.3) + g.integral_between(0.3, 5.0);
        assert!((whole - split).abs() < 1e-12, "{whole} vs {split}");
    }

    #[test]
    fn grid_clamps_queries_outside_its_support() {
        let g = sample(0.0, 1.0, 101, |_| 1.0);
        assert_eq!(g.integral_to(-3.0), 0.0);
        assert!((g.integral_to(9.0) - 1.0).abs() < 1e-12);
        // Reversed bounds must not produce negative probability mass.
        assert_eq!(g.integral_between(0.8, 0.2), 0.0);
    }
}