//! Outcome prediction: who wins, and how likely is a given finishing order.
//!
//! Prediction runs on *performances*, not skills. A competitor's skill is
//! inflated by their performance noise `beta` before any comparison, which is
//! what separates "how good are they" from "how will they do today".
//!
//! Two questions, two algorithms:
//!
//! - **Who finishes first.** Because performances are independent Gaussians,
//!   the probability that team `i` beats every other team separates into a
//!   *one-dimensional* integral — no multivariate orthant integral is
//!   involved. [`quadrature::integrate`] evaluates it to near machine
//!   precision for a few hundred `cdf` calls.
//! - **A specific finishing order.** The factor graph only ever constrains
//!   rank-*adjacent* teams (see `Game::run_chain`), so the joint probability
//!   of a full order is a chain of local constraints rather than a general
//!   orthant probability. That chain collapses into a sequential recursion:
//!   one cumulative integral per adjacent pair, `O(teams * grid)` overall.
//!
//! Both are deterministic. A sampler would have been easier to write and
//! would have made every `predict_*` call return a slightly different number,
//! which is not a property a rating library should have.

use crate::{Gaussian, InferenceError, quadrature};

/// Teams beyond this count make the outcome enumeration impractical.
///
/// Each realisation sorts into exactly one (permutation, tie-pattern) event,
/// so the space has `n! * 2^(n-1)` members: 24 at 3 teams, 192 at 4, 1_920 at
/// 5, 23_040 at 6. The jump to 322_560 at 7 is where enumerating stops being
/// a reasonable thing to do on a caller's behalf.
pub(crate) const MAX_TEAMS_FOR_DISTRIBUTION: usize = 6;

/// Relative tolerance for the first-place integrals.
///
/// The adaptive integrator reaches the exact two-team closed form to ~1e-15 at
/// this tolerance, which is round-off for a probability. `cdf` is no longer the
/// limit — it went to ~1 ULP when `erfc` moved to `libm` — so this is the
/// integrator's own floor.
const WIN_TOLERANCE: f64 = 1e-8;

/// Nodes for the ranking grid, and the floor below which a grid is pointless.
///
/// The recursion converges as O(h^2), so this trades nodes against accuracy
/// directly. Measured against the exact two-team closed form, 2_048 nodes leave
/// ~1.2e-6 of discretisation error and 8_192 reach ~1e-7.
///
/// Unlike the adaptive path there is no approximation floor underneath this any
/// more — `cdf` is accurate to ~1 ULP since `erfc` moved to `libm` — so the
/// error here is purely the grid, and a caller who needs more can only get it
/// by paying for more nodes. 8_192 is the accuracy/cost point chosen, not a
/// 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.
///
/// The normal density is below 1e-18 of its peak past nine sigma, far under
/// the precision of everything else here.
const SUPPORT_SIGMAS: f64 = 9.0;

/// Standard normal CDF at `z`.
fn phi(z: f64) -> f64 {
    crate::cdf(z, 0.0, 1.0)
}

/// Normal density of `x` under `g`.
fn density(g: Gaussian, x: f64) -> f64 {
    let sigma = g.sigma();
    let z = (x - g.mu()) / sigma;
    libm::exp(-0.5 * z * z) / (sigma * (2.0 * std::f64::consts::PI).sqrt())
}

/// Per-pair draw margins.
///
/// The margin is *not* a single number for the whole game: inference derives
/// it per rank-adjacent pair from those two teams' betas (`Game::likelihoods`).
/// Prediction has to use the same per-pair values or it answers a question
/// about a different model than the one that will actually be fitted.
pub(crate) struct Margins {
    n: usize,
    values: Vec<f64>,
}

impl Margins {
    /// Build from a per-pair margin function.
    pub(crate) fn new<F: Fn(usize, usize) -> f64>(n: usize, f: F) -> Self {
        let mut values = vec![0.0; n * n];
        for i in 0..n {
            for j in 0..n {
                if i != j {
                    values[i * n + j] = f(i, j);
                }
            }
        }
        Self { n, values }
    }

    fn get(&self, i: usize, j: usize) -> f64 {
        self.values[i * self.n + j]
    }

    /// True when no pair can draw, so every tie has probability zero.
    fn all_zero(&self) -> bool {
        self.values.iter().all(|&v| v == 0.0)
    }
}

/// `P(team i finishes strictly first)` for every team.
///
/// Strictly means beating each rival by more than that pair's draw margin, so
/// with a non-zero margin these sum to less than one; the shortfall is the
/// probability that the top place is shared.
pub(crate) fn win_probabilities(perf: &[Gaussian], margins: &Margins) -> Vec<f64> {
    (0..perf.len())
        .map(|i| {
            let (mu, sigma) = (perf[i].mu(), perf[i].sigma());
            let (lo, hi) = (mu - SUPPORT_SIGMAS * sigma, mu + SUPPORT_SIGMAS * sigma);

            // Each rival's CDF turns over near its own mean plus the margin.
            // Seeding there is what keeps a rival with a tiny sigma — a step
            // function in disguise — from being stepped over.
            let mut seeds = Vec::with_capacity(3 * perf.len());
            for (j, rival) in perf.iter().enumerate().filter(|&(j, _)| j != i) {
                let centre = rival.mu() + margins.get(i, j);
                seeds.extend_from_slice(&[centre - rival.sigma(), centre, centre + rival.sigma()]);
            }

            quadrature::integrate(
                |x| {
                    let d = density(perf[i], x);
                    if d == 0.0 {
                        return 0.0;
                    }
                    let beaten: f64 = (0..perf.len())
                        .filter(|&j| j != i)
                        .map(|j| phi((x - margins.get(i, j) - perf[j].mu()) / perf[j].sigma()))
                        .product();
                    d * beaten
                },
                lo,
                hi,
                &seeds,
                WIN_TOLERANCE,
            )
        })
        .collect()
}

/// Grid bounds and resolution covering every team's support.
///
/// 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) -> Result<(f64, f64, usize), InferenceError> {
    let lo = perf
        .iter()
        .map(|g| g.mu() - SUPPORT_SIGMAS * g.sigma())
        .fold(f64::INFINITY, f64::min);
    let hi = perf
        .iter()
        .map(|g| g.mu() + SUPPORT_SIGMAS * g.sigma())
        .fold(f64::NEG_INFINITY, f64::max);

    let narrowest = perf
        .iter()
        .map(Gaussian::sigma)
        .fold(f64::INFINITY, f64::min);
    let smallest_margin = margins
        .values
        .iter()
        .copied()
        .filter(|&m| m > 0.0)
        .fold(f64::INFINITY, f64::min);

    let feature = narrowest.min(smallest_margin);
    let wanted = if feature.is_finite() && feature > 0.0 {
        ((hi - lo) / (feature / NODES_PER_FEATURE)).ceil()
    } else {
        MIN_GRID_POINTS as f64
    };

    if !wanted.is_finite() {
        return Ok((lo, hi, MIN_GRID_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.
struct Sampled {
    lo: f64,
    step: f64,
    points: usize,
    density: Vec<Vec<f64>>,
}

impl Sampled {
    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()
            .map(|&g| {
                (0..points)
                    .map(|i| density(g, lo + i as f64 * step))
                    .collect()
            })
            .collect();
        Ok(Self {
            lo,
            step,
            points,
            density,
        })
    }

    fn node(&self, i: usize) -> f64 {
        self.lo + i as f64 * self.step
    }
}

/// `P(order[0] >= order[1] >= ... )` with the given adjacency pattern.
///
/// `tied[k]` says whether `order[k]` and `order[k + 1]` finish within that
/// pair's draw margin. The recursion runs bottom-up: `carry` holds, for each
/// grid node, the probability that everything *below* the current team holds
/// given that team landed on that node. A strict gap reads a cumulative
/// integral; a tie reads a window. Both are O(1) against one prefix array,
/// so each level costs O(grid) and the whole order costs O(teams * grid).
fn order_probability(margins: &Margins, sampled: &Sampled, order: &[usize], tied: &[bool]) -> f64 {
    let mut carry = vec![1.0; sampled.points];

    for k in (0..order.len() - 1).rev() {
        let below = order[k + 1];
        let above = order[k];
        let margin = margins.get(above, below);

        let integrand: Vec<f64> = (0..sampled.points)
            .map(|i| sampled.density[below][i] * carry[i])
            .collect();
        let cumulative = quadrature::Grid::from_values(sampled.lo, sampled.step, integrand);

        carry = (0..sampled.points)
            .map(|i| {
                let x = sampled.node(i);
                if tied[k] {
                    // Sorted order already implies `below <= above`, so the
                    // tie window is one-sided: [x - margin, x].
                    cumulative.integral_between(x - margin, x)
                } else {
                    cumulative.integral_to(x - margin)
                }
            })
            .collect();
    }

    let top = order[0];
    let integrand: Vec<f64> = (0..sampled.points)
        .map(|i| sampled.density[top][i] * carry[i])
        .collect();
    quadrature::Grid::from_values(sampled.lo, sampled.step, integrand).total()
}

/// Dense ranks implied by a sorted order and its tie pattern.
fn ranks_of(order: &[usize], tied: &[bool], n: usize) -> Vec<u32> {
    let mut ranks = vec![0u32; n];
    let mut rank = 0u32;
    ranks[order[0]] = 0;
    for k in 0..order.len() - 1 {
        if !tied[k] {
            rank += 1;
        }
        ranks[order[k + 1]] = rank;
    }
    ranks
}

/// Every (order, tie-pattern) event, or only the strict ones when no pair can
/// draw — a tie then has probability exactly zero and is not worth integrating.
fn events(n: usize, strict_only: bool) -> Vec<(Vec<usize>, Vec<bool>)> {
    fn permute(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
        if k == current.len() {
            out.push(current.clone());
            return;
        }
        for i in k..current.len() {
            current.swap(k, i);
            permute(current, k + 1, out);
            current.swap(k, i);
        }
    }

    let mut orders = Vec::new();
    permute(&mut (0..n).collect(), 0, &mut orders);

    let patterns: Vec<Vec<bool>> = if strict_only {
        vec![vec![false; n - 1]]
    } else {
        (0..(1u32 << (n - 1)))
            .map(|mask| (0..n - 1).map(|i| mask >> i & 1 == 1).collect())
            .collect()
    };

    let mut out = Vec::with_capacity(orders.len() * patterns.len());
    for order in orders {
        for pattern in &patterns {
            out.push((order.clone(), pattern.clone()));
        }
    }
    out
}

/// The full distribution over finishing orders, aggregated by rank vector.
///
/// 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,
) -> Result<Vec<(Vec<u32>, f64)>, InferenceError> {
    let n = perf.len();
    let sampled = Sampled::new(perf, margins)?;

    let mut aggregated: Vec<(Vec<u32>, f64)> = Vec::new();
    for (order, tied) in events(n, margins.all_zero()) {
        let p = order_probability(margins, &sampled, &order, &tied);
        let ranks = ranks_of(&order, &tied, n);
        match aggregated.iter_mut().find(|(r, _)| *r == ranks) {
            Some((_, acc)) => *acc += p,
            None => aggregated.push((ranks, p)),
        }
    }

    aggregated.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
    Ok(aggregated)
}

/// All permutations of `items`.
fn permutations(items: &[usize]) -> Vec<Vec<usize>> {
    fn go(current: &mut Vec<usize>, k: usize, out: &mut Vec<Vec<usize>>) {
        if k == current.len() {
            out.push(current.clone());
            return;
        }
        for i in k..current.len() {
            current.swap(k, i);
            go(current, k + 1, out);
            current.swap(k, i);
        }
    }
    let mut out = Vec::new();
    go(&mut items.to_vec(), 0, &mut out);
    out
}

/// Every (order, tie-pattern) event consistent with a grouping by rank.
///
/// Teams sharing a rank may finish in any internal order, so this is the
/// product of each group's permutations. Adjacencies inside a group are ties;
/// the adjacency joining one group to the next is not.
fn orders_for_groups(groups: &[Vec<usize>]) -> Vec<(Vec<usize>, Vec<bool>)> {
    let per_group: Vec<Vec<Vec<usize>>> = groups.iter().map(|g| permutations(g)).collect();

    let mut out = Vec::new();
    let mut choice = vec![0usize; groups.len()];

    loop {
        let mut order = Vec::new();
        let mut tied = Vec::new();
        for (gi, group) in per_group.iter().enumerate() {
            for (offset, &member) in group[choice[gi]].iter().enumerate() {
                if !order.is_empty() {
                    tied.push(offset != 0);
                }
                order.push(member);
            }
        }
        out.push((order, tied));

        let mut k = 0;
        loop {
            if k == choice.len() {
                return out;
            }
            choice[k] += 1;
            if choice[k] < per_group[k].len() {
                break;
            }
            choice[k] = 0;
            k += 1;
        }
    }
}

/// Probability of one specific rank vector.
///
/// 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],
) -> Result<f64, InferenceError> {
    let n = perf.len();
    let sampled = Sampled::new(perf, margins)?;

    let mut distinct: Vec<u32> = ranks.to_vec();
    distinct.sort_unstable();
    distinct.dedup();

    let groups: Vec<Vec<usize>> = distinct
        .iter()
        .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect())
        .collect();

    Ok(orders_for_groups(&groups)
        .iter()
        .map(|(order, tied)| order_probability(margins, &sampled, order, tied))
        .sum())
}

/// A distribution over the ways a contest could finish.
///
/// Each entry pairs a rank vector — the same shape [`crate::Outcome::ranking`]
/// takes, with equal ranks meaning a tie — against its probability. Entries
/// are ordered most likely first, and cover the whole outcome space, so the
/// probabilities sum to one.
///
/// The rank vectors compose directly with inference: feeding one to
/// `Game::ranked` asks "what would we believe if *this* happened", which is
/// what an expected-information-gain calculation needs alongside the weight.
#[derive(Clone, Debug, PartialEq)]
#[must_use]
pub struct Prediction {
    outcomes: Vec<(Vec<u32>, f64)>,
}

impl Prediction {
    pub(crate) fn new(outcomes: Vec<(Vec<u32>, f64)>) -> Self {
        Self { outcomes }
    }

    /// Every possible finishing order and its probability, most likely first.
    #[must_use]
    pub fn outcomes(&self) -> impl ExactSizeIterator<Item = (&[u32], f64)> {
        self.outcomes.iter().map(|(r, p)| (r.as_slice(), *p))
    }

    /// The single most likely finishing order.
    #[must_use]
    pub fn most_likely(&self) -> Option<(&[u32], f64)> {
        self.outcomes.first().map(|(r, p)| (r.as_slice(), *p))
    }

    /// Probability of one specific finishing order, or zero if it cannot occur.
    #[must_use]
    pub fn probability_of(&self, ranks: &[u32]) -> f64 {
        self.outcomes
            .iter()
            .find(|(r, _)| r.as_slice() == ranks)
            .map_or(0.0, |(_, p)| *p)
    }

    /// `P(team i finishes strictly first)`, for each team.
    ///
    /// Sums to less than one exactly when the top place can be shared; the
    /// shortfall is [`Prediction::shared_first_place`].
    #[must_use]
    pub fn win_probabilities(&self) -> Vec<f64> {
        let n = self.outcomes.first().map_or(0, |(r, _)| r.len());
        let mut wins = vec![0.0; n];
        for (ranks, p) in &self.outcomes {
            let leaders = ranks.iter().filter(|&&r| r == 0).count();
            if leaders == 1 {
                let winner = ranks.iter().position(|&r| r == 0).expect("a rank-0 team");
                wins[winner] += p;
            }
        }
        wins
    }

    /// Probability that two or more teams share first place.
    #[must_use]
    pub fn shared_first_place(&self) -> f64 {
        self.outcomes
            .iter()
            .filter(|(r, _)| r.iter().filter(|&&x| x == 0).count() > 1)
            .map(|(_, p)| p)
            .sum()
    }

    /// Total probability mass, which should be one.
    ///
    /// Exposed because it is a genuine check on the numerics rather than a
    /// formality: the outcome space is exhaustive and disjoint by construction,
    /// so any drift from one is integration error and nothing else.
    #[must_use]
    pub fn total(&self) -> f64 {
        self.outcomes.iter().map(|(_, p)| p).sum()
    }
}

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

    fn g(mu: f64, sigma: f64) -> Gaussian {
        Gaussian::from_ms(mu, sigma)
    }

    fn flat(n: usize, eps: f64) -> Margins {
        Margins::new(n, |_, _| eps)
    }

    /// Exact two-team result: `P(a first) = Phi((mu_a - mu_b - eps) / sd)`.
    fn closed_form_two(a: Gaussian, b: Gaussian, eps: f64) -> (f64, f64) {
        let sd = a.sigma().hypot(b.sigma());
        (
            phi((a.mu() - b.mu() - eps) / sd),
            phi((b.mu() - a.mu() - eps) / sd),
        )
    }

    #[test]
    fn two_team_win_probabilities_match_the_closed_form() {
        for (ma, sa, mb, sb, eps) in [
            (0.0, 6.0, 0.0, 6.0, 0.0),
            (3.0, 6.0, -2.0, 1.0, 0.0),
            (0.0, 6.0, 0.0, 6.0, 2.0),
            (3.0, 6.0, -2.0, 1.0, 1.5),
            (40.0, 1.0, 0.0, 1.0, 0.0),
        ] {
            let perf = [g(ma, sa), g(mb, sb)];
            let got = win_probabilities(&perf, &flat(2, eps));
            let (wa, wb) = closed_form_two(perf[0], perf[1], eps);
            assert!(
                (got[0] - wa).abs() < 1e-12 && (got[1] - wb).abs() < 1e-12,
                "mu=({ma},{mb}) sigma=({sa},{sb}) eps={eps}: got {got:?}, want [{wa}, {wb}]"
            );
        }
    }

    /// The identity that a wrong-but-plausible implementation cannot fake:
    /// with no draw margin, exactly one team finishes first.
    #[test]
    fn win_probabilities_sum_to_one_without_a_draw_margin() {
        for perf in [
            vec![g(0.0, 6.0), g(0.0, 6.0)],
            vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)],
            vec![
                g(8.0, 2.0),
                g(3.0, 6.0),
                g(0.0, 1.0),
                g(-3.0, 4.0),
                g(-8.0, 6.0),
            ],
        ] {
            let sum: f64 = win_probabilities(&perf, &flat(perf.len(), 0.0))
                .iter()
                .sum();
            assert!(
                (sum - 1.0).abs() < 1e-7,
                "{} teams: sum = {sum}",
                perf.len()
            );
        }
    }

    /// A rival with a tiny sigma is a step function in disguise. Fixed-node
    /// quadrature steps over it and lands ~1e-2 out while still looking like a
    /// probability; this is the case that rules that approach out.
    #[test]
    fn win_probabilities_survive_a_rival_with_a_tiny_sigma() {
        let perf = [g(0.0, 0.001), g(0.5, 6.0), g(-0.5, 6.0)];
        let got = win_probabilities(&perf, &flat(3, 0.0));
        let sum: f64 = got.iter().sum();
        assert!((sum - 1.0).abs() < 1e-6, "sum = {sum}, probs = {got:?}");
    }

    #[test]
    fn a_stronger_team_is_more_likely_to_win() {
        let perf = [g(10.0, 3.0), g(0.0, 3.0), g(-10.0, 3.0)];
        let p = win_probabilities(&perf, &flat(3, 0.0));
        assert!(p[0] > p[1] && p[1] > p[2], "not monotone: {p:?}");
    }

    #[test]
    fn identical_teams_are_equally_likely_to_win() {
        let perf = [g(1.0, 4.0), g(1.0, 4.0), g(1.0, 4.0)];
        let p = win_probabilities(&perf, &flat(3, 0.0));
        for probs in p.windows(2) {
            assert!((probs[0] - probs[1]).abs() < 1e-9, "asymmetric: {p:?}");
        }
    }

    /// Every realisation sorts into exactly one finishing order, so the whole
    /// distribution must sum to one — with or without a draw margin.
    #[test]
    fn outcome_distribution_sums_to_one() {
        for (perf, eps) in [
            (vec![g(0.0, 6.0), g(0.0, 6.0)], 0.0),
            (vec![g(0.0, 6.0), g(0.0, 6.0)], 2.0),
            (vec![g(0.0, 6.0), g(0.0, 6.0), g(0.0, 6.0)], 0.0),
            (vec![g(5.0, 6.0), g(0.0, 3.0), g(-5.0, 1.0)], 1.5),
            (vec![g(0.0, 0.05), g(0.5, 6.0), g(-0.5, 6.0)], 1.0),
            (
                vec![g(6.0, 2.0), g(2.0, 6.0), g(-2.0, 1.0), g(-6.0, 4.0)],
                1.0,
            ),
        ] {
            let n = perf.len();
            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,
                "{n} teams, eps={eps}: sum = {sum} over {} outcomes",
                dist.len()
            );
            assert!(dist.iter().all(|(_, p)| *p >= 0.0), "negative probability");
        }
    }

    /// With two teams the distribution is the exact win/draw/loss triple.
    #[test]
    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)).unwrap();
        let (wa, wb) = closed_form_two(perf[0], perf[1], eps);

        let find = |ranks: &[u32]| {
            dist.iter()
                .find(|(r, _)| r == ranks)
                .map_or(0.0, |(_, p)| *p)
        };

        assert!(
            (find(&[0, 1]) - wa).abs() < 1e-6,
            "a wins: {}",
            find(&[0, 1])
        );
        assert!(
            (find(&[1, 0]) - wb).abs() < 1e-6,
            "b wins: {}",
            find(&[1, 0])
        );
        assert!(
            (find(&[0, 0]) - (1.0 - wa - wb)).abs() < 1e-6,
            "draw: {}",
            find(&[0, 0])
        );
    }

    /// Asking for one ranking must agree with that ranking's entry in the
    /// full distribution — the two use different code paths to the same value.
    #[test]
    fn ranking_probability_agrees_with_the_distribution() {
        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).unwrap();

        for (ranks, expected) in &dist {
            let direct = ranking_probability(&perf, &margins, ranks).unwrap();
            assert!(
                (direct - expected).abs() < 1e-9,
                "ranks {ranks:?}: direct {direct} vs distribution {expected}"
            );
        }
    }

    /// Tie mass is controlled by the draw margin. Only the *all-tied* outcome
    /// is monotone in it: every one of its constraints is a window that widens
    /// with the margin. A partially-tied outcome like `[0, 0, 1]` is not, and
    /// must not be asserted to be — widening the margin makes its tie easier
    /// but its "and the last team is strictly behind by more than the margin"
    /// clause harder, so it peaks and then falls.
    #[test]
    fn all_tied_probability_grows_with_the_draw_margin() {
        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]).unwrap();
            assert!(p >= previous, "eps={eps}: {p} < {previous}");
            if eps == 0.0 {
                assert!(p < 1e-12, "a tie needs a margin, got {p}");
            }
            previous = p;
        }
        assert!(
            previous > 0.9,
            "a very wide margin ties everyone: {previous}"
        );
    }

    /// The converse, stated as the non-property it is: a partially-tied
    /// outcome is non-monotone in the margin. Pinning this down stops a future
    /// change from "fixing" it into monotonicity and quietly breaking the model.
    #[test]
    fn a_partially_tied_outcome_peaks_in_the_middle() {
        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]).unwrap())
            .collect();
        let peak = sweep
            .iter()
            .enumerate()
            .fold(
                (0, 0.0),
                |(bi, bv), (i, &v)| if v > bv { (i, v) } else { (bi, bv) },
            )
            .0;
        assert!(
            peak > 0 && peak < sweep.len() - 1,
            "expected an interior peak: {sweep:?}"
        );
    }

    /// With no draw margin a tie has probability exactly zero, and the
    /// enumeration must not waste work pretending otherwise.
    #[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)).unwrap();
        assert_eq!(dist.len(), 6, "expected only the 6 strict orders: {dist:?}");
        assert!(dist.iter().all(|(r, _)| {
            let mut seen = r.clone();
            seen.sort_unstable();
            seen.dedup();
            seen.len() == r.len()
        }));
    }
}