diff --git a/src/error.rs b/src/error.rs index 520e58d..86c6618 100644 --- a/src/error.rs +++ b/src/error.rs @@ -40,6 +40,24 @@ pub enum InferenceError { }, /// Negative precision: a Gaussian with `pi < 0` slipped into an API call. NegativePrecision { pi: f64 }, + /// A prediction referenced a key the history has no skill for. + /// + /// Reported rather than skipped: dropping unknown keys turns a team of + /// strangers into a confident-looking probability about nobody. + UnknownKey { team: usize, member: usize }, + /// A prediction was given a team with no members. + EmptyTeam { team: usize }, + /// Fewer than two teams were supplied to a prediction. + NotEnoughTeams { got: usize }, + /// The full outcome distribution was requested for too many teams. + /// + /// Each realisation sorts into exactly one (order, tie-pattern) event, so + /// the space holds `n! * 2^(n-1)` members — 1_920 at five teams, 23_040 at + /// six, 322_560 at seven. Past `max` this stops being something to + /// enumerate on a caller's behalf; ask for individual rankings with + /// `predict_ranking`, or for `predict_win_probabilities`, both of which + /// stay cheap at any team count. + TooManyTeams { got: usize, max: usize }, } impl fmt::Display for InferenceError { @@ -90,6 +108,25 @@ impl fmt::Display for InferenceError { Self::NegativePrecision { pi } => { write!(f, "precision must be non-negative; got {pi}") } + Self::UnknownKey { team, member } => { + write!( + f, + "team {team}, member {member}: no skill recorded for this key" + ) + } + Self::EmptyTeam { team } => { + write!(f, "team {team} has no members") + } + Self::NotEnoughTeams { got } => { + write!(f, "prediction needs at least 2 teams, got {got}") + } + Self::TooManyTeams { got, max } => { + write!( + f, + "the outcome distribution over {got} teams is too large to enumerate (limit {max}); \ + use predict_ranking or predict_win_probabilities instead" + ) + } } } } diff --git a/src/history.rs b/src/history.rs index 2f34251..1a29bbf 100644 --- a/src/history.rs +++ b/src/history.rs @@ -9,6 +9,7 @@ use crate::{ gaussian::Gaussian, key_table::KeyTable, observer::{NullObserver, Observer}, + predict::Prediction, rating::Rating, sort_time, storage::CompetitorStore, @@ -522,61 +523,213 @@ impl, O: Observer, K: Eq + Hash + Clone> History f64 { - let groups: Vec> = teams - .iter() - .map(|team| { - team.iter() - .filter_map(|k| self.keys.get(*k)) - .filter_map(|idx| { - self.time_slices - .iter() - .rev() - .find_map(|ts| ts.skills.get(idx).map(|s| s.posterior())) - }) - .collect() - }) - .collect(); - let group_refs: Vec<&[Gaussian]> = groups.iter().map(|g| g.as_slice()).collect(); - crate::quality(&group_refs, self.beta) + /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`. Unknown keys are + /// reported rather than dropped — silently skipping them would turn a team + /// of strangers into a confident-looking prediction about nobody, which is + /// the failure this replaced. + fn performances(&self, teams: &[&[&K]]) -> Result<(Vec, Vec), InferenceError> { + if teams.len() < 2 { + return Err(InferenceError::NotEnoughTeams { got: teams.len() }); + } + + let mut performances = Vec::with_capacity(teams.len()); + let mut sizes = Vec::with_capacity(teams.len()); + + for (team_idx, team) in teams.iter().enumerate() { + if team.is_empty() { + return Err(InferenceError::EmptyTeam { team: team_idx }); + } + + let mut total = crate::N00; + for (member_idx, key) in team.iter().enumerate() { + let unknown = InferenceError::UnknownKey { + team: team_idx, + member: member_idx, + }; + let index = self.keys.get(*key).ok_or(unknown.clone())?; + let skill = self + .time_slices + .iter() + .rev() + .find_map(|ts| ts.skills.get(index).map(|s| s.posterior())) + .ok_or(unknown)?; + total = total + skill.forget(self.beta.powi(2)); + } + + performances.push(total); + sizes.push(team.len()); + } + + Ok((performances, sizes)) } - /// 2-team win probability: returns `[P(team0 wins), P(team1 wins)]`. + /// Draw margins per team pair. /// - /// Only two teams are supported. + /// Inference derives the margin per rank-adjacent pair from those two + /// teams' betas (`Game::likelihoods`), so prediction must too — a single + /// game-wide margin would describe a different model than the one that + /// will actually be fitted. + fn margins(&self, sizes: &[usize]) -> crate::predict::Margins { + let beta_sq = self.beta.powi(2); + let p_draw = self.p_draw; + + crate::predict::Margins::new(sizes.len(), |i, j| { + if p_draw == 0.0 { + 0.0 + } else { + let sd = ((sizes[i] + sizes[j]) as f64 * beta_sq).sqrt(); + crate::compute_margin(p_draw, sd) + } + }) + } + + /// Draw-probability quality metric for the given teams (key slices). /// - /// # Panics + /// Values range roughly `[0, 1]`; 1 == perfectly matched. Supports any + /// number of teams. /// - /// Panics if `teams.len() != 2`. - pub fn predict_outcome(&self, teams: &[&[&K]]) -> Vec { - assert_eq!(teams.len(), 2, "predict_outcome supports exactly 2 teams"); - let gather = |team: &[&K]| -> Gaussian { - team.iter() - .filter_map(|k| self.keys.get(*k)) - .filter_map(|idx| { + /// Note this answers "is this matchup *fair*", which is not the same as + /// "is this matchup *informative*" — the two coincide for two evenly + /// matched teams and diverge elsewhere. + /// + /// # Errors + /// + /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`. + pub fn predict_quality(&self, teams: &[&[&K]]) -> Result { + let mut groups: Vec> = Vec::with_capacity(teams.len()); + + for (team_idx, team) in teams.iter().enumerate() { + if team.is_empty() { + return Err(InferenceError::EmptyTeam { team: team_idx }); + } + let mut members = Vec::with_capacity(team.len()); + for (member_idx, key) in team.iter().enumerate() { + let unknown = InferenceError::UnknownKey { + team: team_idx, + member: member_idx, + }; + let index = self.keys.get(*key).ok_or(unknown.clone())?; + members.push( self.time_slices .iter() .rev() - .find_map(|ts| ts.skills.get(idx).map(|s| s.posterior())) - }) - .fold(crate::N00, |acc, g| acc + g.forget(self.beta.powi(2))) - }; - let a = gather(teams[0]); - let b = gather(teams[1]); - let diff = a - b; - let p_a = 1.0 - crate::cdf(0.0, diff.mu(), diff.sigma()); - vec![p_a, 1.0 - p_a] + .find_map(|ts| ts.skills.get(index).map(|s| s.posterior())) + .ok_or(unknown)?, + ); + } + groups.push(members); + } + + if groups.len() < 2 { + return Err(InferenceError::NotEnoughTeams { got: groups.len() }); + } + + let group_refs: Vec<&[Gaussian]> = groups.iter().map(Vec::as_slice).collect(); + Ok(crate::quality(&group_refs, self.beta)) + } + + /// `P(team i finishes strictly first)`, for every team. + /// + /// Supports any number of teams. Because performances are independent + /// Gaussians, this separates into a one-dimensional integral per team — + /// no multivariate orthant probability is involved — and is evaluated by + /// adaptive quadrature to within the precision of the underlying normal + /// CDF (~1e-8). + /// + /// With a zero `p_draw` these sum to one. With a positive `p_draw` the + /// shortfall is the probability that the top place is shared. + /// + /// Cheap at any team count: cost grows as the square of the team count, + /// not factorially. Prefer this to [`History::predict_outcome`] when you + /// only need to know who wins. + /// + /// # Errors + /// + /// `NotEnoughTeams`, `EmptyTeam`, or `UnknownKey`. + pub fn predict_win_probabilities(&self, teams: &[&[&K]]) -> Result, InferenceError> { + let (performances, sizes) = self.performances(teams)?; + Ok(crate::predict::win_probabilities( + &performances, + &self.margins(&sizes), + )) + } + + /// The full distribution over finishing orders. + /// + /// Every entry is a rank vector — the shape [`crate::Outcome::ranking`] + /// takes, equal ranks meaning a tie — paired with its probability. The + /// entries are exhaustive and disjoint, so they sum to one; that identity + /// is the strongest available check on the numerics and is worth asserting + /// in tests via [`Prediction::total`]. + /// + /// Accounts for `p_draw`: with a positive draw probability, tied outcomes + /// carry real mass rather than being silently omitted. + /// + /// # Cost + /// + /// This enumerates the outcome space, which holds `n! * 2^(n-1)` events — + /// 24 at three teams, 192 at four, 1_920 at five, 23_040 at six. Each + /// costs one `O(teams * grid)` pass, so this is milliseconds at three or + /// four teams and seconds at six. Above + /// [`MAX_TEAMS_FOR_DISTRIBUTION`](crate::MAX_PREDICTED_TEAMS) it returns + /// `TooManyTeams` rather than hanging. When you need one specific ordering + /// use [`History::predict_ranking`], and when you only need the winner use + /// [`History::predict_win_probabilities`]; both stay cheap at any size. + /// + /// # Errors + /// + /// `NotEnoughTeams`, `EmptyTeam`, `UnknownKey`, or `TooManyTeams`. + pub fn predict_outcome(&self, teams: &[&[&K]]) -> Result { + if teams.len() > crate::MAX_PREDICTED_TEAMS { + return Err(InferenceError::TooManyTeams { + got: teams.len(), + max: crate::MAX_PREDICTED_TEAMS, + }); + } + + let (performances, sizes) = self.performances(teams)?; + Ok(Prediction::new(crate::predict::outcome_distribution( + &performances, + &self.margins(&sizes), + ))) + } + + /// Probability of one specific finishing order. + /// + /// `ranks` follows [`crate::Outcome::ranking`]: lower is better, and equal + /// values mean those teams tied. Teams sharing a rank may finish in any + /// internal order, so this sums over those orders rather than picking one. + /// + /// Unlike [`History::predict_outcome`] this does not enumerate the outcome + /// space, so it stays cheap at any team count — use it when you know which + /// orderings you care about. + /// + /// # Errors + /// + /// `NotEnoughTeams`, `EmptyTeam`, `UnknownKey`, or `MismatchedShape` if + /// `ranks` does not have one entry per team. + pub fn predict_ranking(&self, teams: &[&[&K]], ranks: &[u32]) -> Result { + if ranks.len() != teams.len() { + return Err(InferenceError::MismatchedShape { + kind: "ranks vs teams", + expected: teams.len(), + got: ranks.len(), + }); + } + + let (performances, sizes) = self.performances(teams)?; + Ok(crate::predict::ranking_probability( + &performances, + &self.margins(&sizes), + ranks, + )) } /// Run the full forward+backward convergence loop and return a summary. diff --git a/src/lib.rs b/src/lib.rs index 6ddfaac..db2c9d9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -126,6 +126,8 @@ mod key_table; mod matrix; mod observer; mod outcome; +mod predict; +pub(crate) mod quadrature; mod rating; pub(crate) mod schedule; pub mod storage; @@ -143,6 +145,7 @@ pub use key_table::KeyTable; use matrix::Matrix; pub use observer::{NullObserver, Observer}; pub use outcome::Outcome; +pub use predict::Prediction; pub use rating::Rating; pub use schedule::ScheduleReport; pub use time::{Time, Untimed}; @@ -155,6 +158,13 @@ pub const P_DRAW: f64 = 0.0; pub const EPSILON: f64 = 1e-6; pub const ITERATIONS: usize = 30; +/// Largest team count `History::predict_outcome` will enumerate. +/// +/// The outcome space holds `n! * 2^(n-1)` events, so it grows factorially: +/// 1_920 at five teams, 23_040 at six, 322_560 at seven. Six is where +/// enumerating on a caller's behalf stops being reasonable. +pub const MAX_PREDICTED_TEAMS: usize = predict::MAX_TEAMS_FOR_DISTRIBUTION; + const SQRT_TAU: f64 = 2.5066282746310002; pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0); diff --git a/src/predict.rs b/src/predict.rs new file mode 100644 index 0000000..ae4a258 --- /dev/null +++ b/src/predict.rs @@ -0,0 +1,723 @@ +//! 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, 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. +/// +/// Tightening past this buys nothing: the underlying `cdf` is a rational +/// approximation with fractional error ~1.2e-7, which contributes ~6e-9 to a +/// finished probability and dominates any further quadrature refinement. +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). Measured against the exact two-team +/// closed form, 2_048 nodes leave ~1.2e-6 of discretisation error while 8_192 +/// reach ~1e-7 — at which point the residual is the `cdf` rational +/// approximation (~2.4e-8), not the grid, and refining further buys nothing. +const MIN_GRID_POINTS: usize = 8_192; +const MAX_GRID_POINTS: usize = 262_144; + +/// 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; + (-0.5 * z * z).exp() / (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, +} + +impl Margins { + /// Build from a per-pair margin function. + pub(crate) fn new 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 { + (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) -> (f64, f64, usize) { + 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 / 12.0)).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 + }; + + (lo, hi, points) +} + +/// Densities of each team sampled on the shared grid. +struct Sampled { + lo: f64, + step: f64, + points: usize, + density: Vec>, +} + +impl Sampled { + fn new(perf: &[Gaussian], margins: &Margins) -> Self { + 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(); + 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 = (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 = (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 { + 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, Vec)> { + fn permute(current: &mut Vec, k: usize, out: &mut Vec>) { + 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> = 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) -> Vec<(Vec, f64)> { + let n = perf.len(); + let sampled = Sampled::new(perf, margins); + + let mut aggregated: Vec<(Vec, 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)); + aggregated +} + +/// All permutations of `items`. +fn permutations(items: &[usize]) -> Vec> { + fn go(current: &mut Vec, k: usize, out: &mut Vec>) { + 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]) -> Vec<(Vec, Vec)> { + let per_group: Vec>> = 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]) -> f64 { + let n = perf.len(); + let sampled = Sampled::new(perf, margins); + + let mut distinct: Vec = ranks.to_vec(); + distinct.sort_unstable(); + distinct.dedup(); + + let groups: Vec> = distinct + .iter() + .map(|&r| (0..n).filter(|&i| ranks[i] == r).collect()) + .collect(); + + 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)] +pub struct Prediction { + outcomes: Vec<(Vec, f64)>, +} + +impl Prediction { + pub(crate) fn new(outcomes: Vec<(Vec, f64)>) -> Self { + Self { outcomes } + } + + /// Every possible finishing order and its probability, most likely first. + pub fn outcomes(&self) -> impl ExactSizeIterator { + 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 { + 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().powi(2) + b.sigma().powi(2)).sqrt(); + ( + 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-7 && (got[1] - wb).abs() < 1e-7, + "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)); + 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)); + 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); + + for (ranks, expected) in &dist { + let direct = ranking_probability(&perf, &margins, ranks); + 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]); + 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 = [0.5, 2.0, 4.0, 8.0, 16.0] + .iter() + .map(|&eps| ranking_probability(&perf, &flat(3, eps), &[0, 0, 1])) + .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)); + 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() + })); + } +} diff --git a/src/quadrature.rs b/src/quadrature.rs new file mode 100644 index 0000000..33bcce7 --- /dev/null +++ b/src/quadrature.rs @@ -0,0 +1,322 @@ +//! 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 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 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 = 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, + /// `prefix[i]` is the integral from `lo` to node `i`. + prefix: Vec, +} + +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) -> 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 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); + } +} diff --git a/tests/api_shape.rs b/tests/api_shape.rs index 2ac8971..0331be5 100644 --- a/tests/api_shape.rs +++ b/tests/api_shape.rs @@ -203,7 +203,7 @@ fn predict_quality_two_teams() { h.record_winner(&"a", &"b", 1).unwrap(); h.converge().unwrap(); - let q = h.predict_quality(&[&[&"a"], &[&"b"]]); + let q = h.predict_quality(&[&[&"a"], &[&"b"]]).unwrap(); assert!(q > 0.0 && q <= 1.0); } @@ -219,10 +219,14 @@ fn predict_outcome_two_teams_sums_to_one() { h.record_winner(&"a", &"b", 1).unwrap(); h.converge().unwrap(); - let p = h.predict_outcome(&[&[&"a"], &[&"b"]]); - assert_eq!(p.len(), 2); - assert!((p[0] + p[1] - 1.0).abs() < 1e-9); - assert!(p[0] > p[1]); + let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap(); + let wins = p.win_probabilities(); + assert_eq!(wins.len(), 2); + // With p_draw == 0 there is no draw outcome, so the two win + // probabilities are the whole space. + assert!((p.total() - 1.0).abs() < 1e-9, "total = {}", p.total()); + assert!((wins[0] + wins[1] - 1.0).abs() < 1e-9); + assert!(wins[0] > wins[1]); } #[test] diff --git a/tests/prediction.rs b/tests/prediction.rs new file mode 100644 index 0000000..bee0be8 --- /dev/null +++ b/tests/prediction.rs @@ -0,0 +1,214 @@ +//! Prediction API: N-team outcomes, draw mass, and the error paths that used +//! to be panics or silent wrong answers. + +use trueskill_tt::{History, InferenceError, MAX_PREDICTED_TEAMS}; + +fn history_with(names: &[&'static str], p_draw: f64) -> History { + let mut h = History::builder().p_draw(p_draw).build(); + // Give every competitor a recorded skill by playing a small round robin. + for pair in names.windows(2) { + h.record_winner(&pair[0], &pair[1], 1).unwrap(); + } + h.converge().unwrap(); + h +} + +#[test] +fn unknown_keys_are_reported_not_silently_dropped() { + let h = history_with(&["a", "b"], 0.0); + + let err = h + .predict_outcome(&[&[&"a"], &[&"ghost"]]) + .expect_err("an unknown key must not yield a confident prediction"); + assert_eq!(err, InferenceError::UnknownKey { team: 1, member: 0 }); + + // Every prediction entry point, not just one. + assert!( + h.predict_win_probabilities(&[&[&"a"], &[&"ghost"]]) + .is_err() + ); + assert!(h.predict_quality(&[&[&"a"], &[&"ghost"]]).is_err()); + assert!(h.predict_ranking(&[&[&"a"], &[&"ghost"]], &[0, 1]).is_err()); +} + +#[test] +fn an_entirely_unknown_team_is_an_error() { + let h = history_with(&["a", "b"], 0.0); + let err = h.predict_outcome(&[&[&"a"], &[&"x", &"y"]]).unwrap_err(); + assert_eq!(err, InferenceError::UnknownKey { team: 1, member: 0 }); +} + +#[test] +fn degenerate_team_shapes_are_errors_rather_than_panics() { + let h = history_with(&["a", "b"], 0.0); + + assert_eq!( + h.predict_outcome(&[&[&"a"]]).unwrap_err(), + InferenceError::NotEnoughTeams { got: 1 } + ); + assert_eq!( + h.predict_outcome(&[]).unwrap_err(), + InferenceError::NotEnoughTeams { got: 0 } + ); + assert_eq!( + h.predict_outcome(&[&[&"a"], &[]]).unwrap_err(), + InferenceError::EmptyTeam { team: 1 } + ); +} + +#[test] +fn more_than_two_teams_no_longer_panics() { + let h = history_with(&["a", "b", "c"], 0.0); + let p = h + .predict_outcome(&[&[&"a"], &[&"b"], &[&"c"]]) + .expect("three teams must be supported"); + assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total()); + // Three teams, no draws possible: exactly the six strict orderings. + assert_eq!(p.outcomes().len(), 6); +} + +#[test] +fn the_outcome_space_is_capped_rather_than_hanging() { + let names: Vec<&'static str> = vec!["a", "b", "c", "d", "e", "f", "g", "h"]; + let h = history_with(&names, 0.0); + + let teams: Vec<&[&&'static str]> = Vec::new(); + let _ = teams; + + let too_many: Vec> = names.iter().map(|n| vec![n]).collect(); + let refs: Vec<&[&&str]> = too_many.iter().map(Vec::as_slice).collect(); + + let err = h.predict_outcome(&refs).unwrap_err(); + assert_eq!( + err, + InferenceError::TooManyTeams { + got: 8, + max: MAX_PREDICTED_TEAMS + } + ); + + // The cheap paths stay available at any size. + let wins = h.predict_win_probabilities(&refs).unwrap(); + assert_eq!(wins.len(), 8); + assert!( + (wins.iter().sum::() - 1.0).abs() < 1e-6, + "win probabilities must still sum to one: {wins:?}" + ); +} + +/// The defect that made every draw-enabled prediction wrong: `[p, 1 - p]` +/// allocated no mass to a draw even with `p_draw > 0`. +#[test] +fn a_draw_carries_probability_mass_when_p_draw_is_positive() { + let h = history_with(&["a", "b"], 0.25); + let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap(); + + let draw = p.probability_of(&[0, 0]); + assert!(draw > 0.0, "a draw-enabled model must give draws mass"); + assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total()); + + let wins = p.win_probabilities(); + assert!( + (wins.iter().sum::() + draw - 1.0).abs() < 1e-6, + "wins {wins:?} plus draw {draw} must be the whole space" + ); + assert!( + (p.shared_first_place() - draw).abs() < 1e-12, + "a two-team draw is a shared first place" + ); +} + +#[test] +fn a_zero_draw_probability_admits_no_ties() { + let h = history_with(&["a", "b"], 0.0); + let p = h.predict_outcome(&[&[&"a"], &[&"b"]]).unwrap(); + assert_eq!(p.probability_of(&[0, 0]), 0.0); + assert!(p.shared_first_place() < 1e-12); +} + +/// The two routes to a win probability run through entirely different +/// algorithms — adaptive quadrature versus the enumerated chain recursion — +/// so agreement between them is a real cross-check, not a tautology. +#[test] +fn the_cheap_and_exhaustive_paths_agree() { + for p_draw in [0.0, 0.1] { + let h = history_with(&["a", "b", "c"], p_draw); + let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]]; + + let cheap = h.predict_win_probabilities(teams).unwrap(); + let exhaustive = h.predict_outcome(teams).unwrap().win_probabilities(); + + for (i, (a, b)) in cheap.iter().zip(&exhaustive).enumerate() { + assert!( + (a - b).abs() < 1e-6, + "p_draw={p_draw} team {i}: quadrature {a} vs enumeration {b}" + ); + } + } +} + +#[test] +fn predict_ranking_agrees_with_the_distribution() { + let h = history_with(&["a", "b", "c"], 0.1); + let teams: &[&[&&str]] = &[&[&"a"], &[&"b"], &[&"c"]]; + let dist = h.predict_outcome(teams).unwrap(); + + for (ranks, expected) in dist.outcomes() { + let direct = h.predict_ranking(teams, ranks).unwrap(); + assert!( + (direct - expected).abs() < 1e-9, + "ranks {ranks:?}: {direct} vs {expected}" + ); + } +} + +#[test] +fn predict_ranking_checks_its_shape() { + let h = history_with(&["a", "b"], 0.0); + let err = h + .predict_ranking(&[&[&"a"], &[&"b"]], &[0, 1, 2]) + .unwrap_err(); + assert!(matches!( + err, + InferenceError::MismatchedShape { + expected: 2, + got: 3, + .. + } + )); +} + +#[test] +fn the_stronger_competitor_is_favoured() { + let mut h = History::builder().build(); + for t in 1..=10 { + h.record_winner(&"strong", &"weak", t).unwrap(); + } + h.converge().unwrap(); + + let p = h.predict_outcome(&[&[&"strong"], &[&"weak"]]).unwrap(); + let (best, _) = p.most_likely().expect("a most likely outcome"); + assert_eq!(best, &[0, 1], "the winner should be favoured"); + + let wins = p.win_probabilities(); + assert!(wins[0] > wins[1], "{wins:?}"); +} + +/// Unequal team sizes change the draw margin, because inference derives it +/// from the teams' betas. Prediction has to follow, or it describes a +/// different model than the one that will be fitted. +#[test] +fn team_size_affects_the_prediction() { + let mut h = History::builder().p_draw(0.2).build(); + h.event(1) + .team(["a", "b"]) + .team(["c"]) + .winner(0) + .commit() + .unwrap(); + h.converge().unwrap(); + + let p = h.predict_outcome(&[&[&"a", &"b"], &[&"c"]]).unwrap(); + assert!((p.total() - 1.0).abs() < 1e-6, "total = {}", p.total()); + assert!(p.probability_of(&[0, 0]) > 0.0); +} diff --git a/tests/quality.rs b/tests/quality.rs index fe7c2c4..3dcea0f 100644 --- a/tests/quality.rs +++ b/tests/quality.rs @@ -110,7 +110,7 @@ fn history_predict_quality_supports_three_teams() { h.record_winner(&"b", &"c", 2).unwrap(); h.converge().unwrap(); - let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]); + let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]).unwrap(); assert!( q.is_finite(), "3-team predict_quality must be finite, got {q}"