diff --git a/Cargo.toml b/Cargo.toml index f8efdb7..713c421 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -79,3 +79,7 @@ debug = true [profile.dev] debug = true + +[[bench]] +name = "joint" +harness = false diff --git a/benches/joint.rs b/benches/joint.rs new file mode 100644 index 0000000..eb5aff1 --- /dev/null +++ b/benches/joint.rs @@ -0,0 +1,71 @@ +//! Cost of the joint posterior: factorising versus querying. +//! +//! The split is the whole point of `History::joint`. Factorising is `O(n^3)` in +//! the history's appearances and depends only on the fit; a query is `O(n^2)` +//! and depends only on the question. `posterior_of_one_shot` pays both every +//! time, `joint_query` pays only the second. + +use criterion::{Criterion, criterion_group, criterion_main}; +use smallvec::smallvec; +use trueskill_tt::{ConstantDrift, ConvergenceOptions, Event, History, Member, Outcome, Team}; + +/// 30 slices of 8 duels: 480 appearances over 100 competitors. +fn fitted() -> History { + let mut h: History = History::builder_with_key() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .score_sigma(2.0) + .drift(ConstantDrift(0.05)) + .convergence(ConvergenceOptions { + max_iter: 30, + epsilon: 1e-10, + alpha: 1.0, + }) + .build(); + + let mut events: Vec> = Vec::new(); + let mut k = 0usize; + for t in 0..30i64 { + for _ in 0..8 { + k += 1; + events.push(Event { + time: t, + teams: smallvec![ + Team::with_members([Member::new(format!("p{}", k % 100))]), + Team::with_members([Member::new(format!("p{}", (k + 37) % 100))]), + ], + outcome: Outcome::scores([ + (k as f64 * 0.3).sin().abs() * 20.0, + (k as f64 * 0.3).cos().abs() * 20.0, + ]), + }); + } + } + h.add_events(events).unwrap(); + let _ = h.converge().unwrap(); + h +} + +fn bench_joint(c: &mut Criterion) { + let h = fitted(); + let a = "p0".to_string(); + let b = "p1".to_string(); + let terms = [(&a, 1.0), (&b, -1.0)]; + + c.bench_function("joint_factorise_480_appearances", |bencher| { + bencher.iter(|| std::hint::black_box(h.joint().unwrap().variables())); + }); + + c.bench_function("posterior_of_one_shot_480_appearances", |bencher| { + bencher.iter(|| std::hint::black_box(h.posterior_of(&terms).unwrap())); + }); + + let joint = h.joint().unwrap(); + c.bench_function("joint_query_480_appearances", |bencher| { + bencher.iter(|| std::hint::black_box(joint.posterior_of(&terms).unwrap())); + }); +} + +criterion_group!(benches, bench_joint); +criterion_main!(benches); diff --git a/src/history.rs b/src/history.rs index f5af948..989138a 100644 --- a/src/history.rs +++ b/src/history.rs @@ -950,6 +950,14 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History() - + unseen.values().map(|c| c * c * prior_var).sum::(); - - Ok(Gaussian::from_mv(mean, variance)) + self.joint()?.posterior_of(terms) } /// Posterior of a linear combination, read as of `time`. @@ -1018,6 +987,9 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History = HashMap::new(); - for (slice_idx, slice) in self.time_slices.iter().enumerate() { - if slice.time > time { - break; - } - for (agent, _) in slice.appearances() { - if let Some(row) = at_slice.get(&(agent, slice_idx)) { - as_of.insert(agent, (*row, slice_idx)); - } - } - } - - let ResolvedTerms { - contrast, - unseen, - mean, - } = self.resolve_terms(terms, width, |index| as_of.get(&index).copied())?; - - let z = - crate::joint::solve_spd(lambda, &contrast).ok_or(InferenceError::JointUnavailable { - reason: "the precision matrix is not positive-definite", - })?; - - let prior_var = self.sigma * self.sigma; - let variance: f64 = contrast.iter().zip(&z).map(|(c, z)| c * z).sum::() - + unseen.values().map(|c| c * c * prior_var).sum::(); - - Ok(Gaussian::from_mv(mean, variance)) + self.joint()?.posterior_of_at(time, terms) } /// How much observing this matchup would shrink the variance of `target`. @@ -1089,6 +1014,11 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> History(()) + /// ``` + /// + /// The handle borrows the history, so the borrow checker enforces what a + /// cache would otherwise have to invalidate: no events can be added and no + /// refit can run while it is alive. Drop it to release the factorisation, + /// which is `n^2` floats and is the largest thing this crate allocates. + /// + /// # Errors + /// + /// `JointUnavailable` if the history is empty, contains ranked events, or + /// yields a precision matrix that is not positive-definite. + pub fn joint(&self) -> Result, InferenceError> { if self.time_slices.is_empty() { return Err(InferenceError::JointUnavailable { reason: "the history has no events", @@ -1138,70 +1110,28 @@ impl, O: Observer, K: Eq + Hash + Clone> History = Vec::new(); - let mut noise = self.score_sigma * self.score_sigma; - for (team_idx, team) in teams.iter().enumerate() { - if team.is_empty() { - return Err(InferenceError::EmptyTeam { team: team_idx }); - } - let sign = if team_idx == 0 { 1.0 } else { -1.0 }; - for key in team.iter() { - matchup.push((*key, sign)); - let beta = self - .keys - .get(*key) - .map_or(self.beta, |index| self.agents[index].rating.beta); - noise += beta * beta; - } - } - let TimeExpanded { lambda, latest, + at_slice, width, - .. } = self.time_expanded_joint(); - let target = self.resolve_terms(target, width, |i| latest.get(&i).copied())?; - let matchup = self.resolve_terms(&matchup, width, |i| latest.get(&i).copied())?; - let (target_contrast, target_unseen) = (target.contrast, target.unseen); - let (matchup_contrast, matchup_unseen) = (matchup.contrast, matchup.unseen); - - // One solve: z = L^-1 a serves both inner products, since - // c^T L^-1 a = c^T z and a^T L^-1 a = a^T z. - let z = crate::joint::solve_spd(lambda, &matchup_contrast).ok_or( + let cholesky = crate::joint::Cholesky::factor(lambda, width).ok_or( InferenceError::JointUnavailable { - reason: "the precision matrix is not positive-definite", + reason: "the precision matrix is not positive-definite, which means \ + a competitor has neither a proper prior nor any evidence", }, )?; - let prior_var = self.sigma * self.sigma; - // Competitors outside the history are independent, so they contribute - // only where the same key appears in both functionals. - let cross_unseen: f64 = target_unseen - .iter() - .map(|(k, tc)| tc * matchup_unseen.get(k).copied().unwrap_or(0.0) * prior_var) - .sum(); - let self_unseen: f64 = matchup_unseen.values().map(|c| c * c * prior_var).sum(); - - let cross: f64 = target_contrast - .iter() - .zip(&z) - .map(|(c, z)| c * z) - .sum::() - + cross_unseen; - let matchup_var: f64 = matchup_contrast - .iter() - .zip(&z) - .map(|(a, z)| a * z) - .sum::() - + self_unseen; - - Ok(cross * cross / (noise + matchup_var)) + Ok(Joint { + history: self, + cholesky, + latest, + at_slice, + width, + }) } - /// Predictive distribution of the score margin between two teams. /// /// Answers "what will the gap be, and how wide is that interval" for a @@ -1980,6 +1910,195 @@ impl, O: Observer, K: Eq + Hash + Clone> History, O: Observer, K: Eq + Hash + Clone> { + history: &'h History, + cholesky: crate::joint::Cholesky, + /// `(row, slice)` of each competitor's latest appearance. + latest: HashMap, + /// Row of each `(competitor, slice)` appearance. + at_slice: HashMap<(Index, usize), usize>, + /// Side length of the precision matrix. + width: usize, +} + +/// Deliberately does not print the factorisation, which is `n^2` floats and +/// would make a `{:?}` of a large joint unreadable and slow. +impl, O: Observer, K: Eq + Hash + Clone> std::fmt::Debug + for Joint<'_, T, D, O, K> +{ + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + f.debug_struct("Joint") + .field("variables", &self.width) + .finish_non_exhaustive() + } +} + +impl, O: Observer, K: Eq + Hash + Clone> Joint<'_, T, D, O, K> { + /// Number of variables in the joint: the history's appearances, after + /// collapsing consecutive pairs a competitor does not drift between. + /// + /// This is what the cost scales in, and it is not the competitor count — a + /// competitor contributes one variable per slice it appears in. Worth + /// checking before asking for a joint over a long history. + #[must_use] + pub fn variables(&self) -> usize { + self.width + } + + /// Turn a resolved functional into its posterior. + /// + /// The variance is `|L^-1 c|^2` over the competitors the history knows, + /// plus an independent prior variance for each competitor it does not — + /// unseen competitors are uncorrelated with everything by construction. + fn distribution(&self, resolved: &ResolvedTerms) -> Gaussian { + let y = self.cholesky.whiten(&resolved.contrast); + let prior_var = self.history.sigma * self.history.sigma; + let variance = crate::joint::bilinear(&y, &y) + + resolved + .unseen + .values() + .map(|c| c * c * prior_var) + .sum::(); + Gaussian::from_mv(resolved.mean, variance) + } + + /// Posterior of a linear combination of competitors' skills. + /// + /// Identical to [`History::posterior_of`], including which appearance each + /// competitor is read at, without re-paying the factorisation. + /// + /// # Errors + /// + /// `UnknownKey` for a competitor the history has never seen. + pub fn posterior_of(&self, terms: &[(&K, f64)]) -> Result + where + K: std::fmt::Debug, + { + let resolved = self + .history + .resolve_terms(terms, self.width, |index| self.latest.get(&index).copied())?; + Ok(self.distribution(&resolved)) + } + + /// Posterior of a linear combination, read as of `time`. + /// + /// Identical to [`History::posterior_of_at`] without re-paying the + /// factorisation. + /// + /// # Errors + /// + /// `UnknownKey` for a competitor with no appearance at or before `time`. + pub fn posterior_of_at(&self, time: T, terms: &[(&K, f64)]) -> Result + where + K: std::fmt::Debug, + { + let as_of = self.rows_as_of(time); + let resolved = self + .history + .resolve_terms(terms, self.width, |index| as_of.get(&index).copied())?; + Ok(self.distribution(&resolved)) + } + + /// Latest appearance at or before `time`, per competitor. + fn rows_as_of(&self, time: T) -> HashMap { + let mut as_of: HashMap = HashMap::new(); + for (slice_idx, slice) in self.history.time_slices.iter().enumerate() { + if slice.time > time { + break; + } + for (agent, _) in slice.appearances() { + if let Some(row) = self.at_slice.get(&(agent, slice_idx)) { + as_of.insert(agent, (*row, slice_idx)); + } + } + } + as_of + } + + /// How much observing this matchup would shrink the variance of `target`. + /// + /// Identical to [`History::expected_variance_reduction`] without re-paying + /// the factorisation, which is the shape this call is normally used in: + /// one target, a field of candidate matchups, one unchanged fit. + /// + /// # Errors + /// + /// `MismatchedShape` unless exactly two teams are supplied, `EmptyTeam` for + /// an empty one, and `UnknownKey` for an unseen competitor. + pub fn expected_variance_reduction( + &self, + teams: &[&[&K]], + target: &[(&K, f64)], + ) -> Result + where + K: std::fmt::Debug, + { + if teams.len() != 2 { + return Err(InferenceError::MismatchedShape { + kind: "expected_variance_reduction takes exactly 2 teams", + expected: 2, + got: teams.len(), + }); + } + + // The candidate matchup, expressed as the same kind of linear + // functional as the target. + let mut matchup: Vec<(&K, f64)> = Vec::new(); + let mut noise = self.history.score_sigma * self.history.score_sigma; + for (team_idx, team) in teams.iter().enumerate() { + if team.is_empty() { + return Err(InferenceError::EmptyTeam { team: team_idx }); + } + let sign = if team_idx == 0 { 1.0 } else { -1.0 }; + for key in team.iter() { + matchup.push((*key, sign)); + let beta = self + .history + .keys + .get(*key) + .map_or(self.history.beta, |index| { + self.history.agents[index].rating.beta + }); + noise += beta * beta; + } + } + + let row_for = |index: Index| self.latest.get(&index).copied(); + let target = self.history.resolve_terms(target, self.width, row_for)?; + let matchup = self.history.resolve_terms(&matchup, self.width, row_for)?; + + let y_target = self.cholesky.whiten(&target.contrast); + let y_matchup = self.cholesky.whiten(&matchup.contrast); + + let prior_var = self.history.sigma * self.history.sigma; + // Competitors outside the history are independent, so they contribute + // only where the same key appears in both functionals. + let cross_unseen: f64 = target + .unseen + .iter() + .map(|(k, tc)| tc * matchup.unseen.get(k).copied().unwrap_or(0.0) * prior_var) + .sum(); + let self_unseen: f64 = matchup.unseen.values().map(|c| c * c * prior_var).sum(); + + let cross = crate::joint::bilinear(&y_target, &y_matchup) + cross_unseen; + let matchup_var = crate::joint::bilinear(&y_matchup, &y_matchup) + self_unseen; + + Ok(cross * cross / (noise + matchup_var)) + } +} + #[cfg(test)] mod tests { use approx::assert_ulps_eq; diff --git a/src/joint.rs b/src/joint.rs index 7746126..ae4c266 100644 --- a/src/joint.rs +++ b/src/joint.rs @@ -1,99 +1,152 @@ -//! Posterior of a linear combination of competitors. +//! Cholesky factorisation of a joint precision matrix. //! -//! Every accessor on `History` returns a per-competitor marginal, and almost -//! nothing a consumer publishes is one competitor: "can we tell these two -//! apart" is a difference, "what was this round worth" is a sum. Combining -//! marginals means assuming the competitors are independent, and they are -//! correlated through every event they share — which is the mechanism the model -//! exists to exploit. +//! Every question the joint answers is a *bilinear form* in the precision +//! matrix's inverse — the variance of a contrast is `c^T L^-1 c`, and the +//! covariance of two contrasts is `c^T L^-1 a`. None of them wants `L^-1 c` +//! itself, which is what makes the shape here worth stating explicitly. //! -//! Measured on a five-competitor round robin, the exact correlation is +0.857, -//! so `sqrt(sa^2 + sb^2)` overstates the width of a difference by 2.6x. +//! Writing the precision as `A = L L^T`, +//! +//! ```text +//! c^T A^-1 a = c^T L^-T L^-1 a = (L^-1 c) . (L^-1 a) +//! ``` +//! +//! so a single forward substitution per contrast answers everything, and the +//! back substitution a general solve would do is wasted work. That halves the +//! cost of a query, and it removes a failure mode: a variance computed as +//! `c . (A^-1 c)` is a difference of products that can round to a small +//! negative number, where the same quantity as `|L^-1 c|^2` is a sum of +//! squares and cannot. +//! +//! Factorising is `O(n^3)` and whitening is `O(n^2)`, so the split also +//! matters structurally: the expensive half depends only on the fit, and is +//! shared across every query a [`Joint`](crate::Joint) answers. -/// Solve `A z = b` for a symmetric positive-definite `A`, by Cholesky. -/// -/// `a` is row-major and is consumed as scratch. -/// -/// Returns `None` if the matrix is not positive-definite, which for a precision -/// matrix means the model is improper — a competitor with no prior and no -/// evidence. -pub(crate) fn solve_spd(mut a: Vec, b: &[f64]) -> Option> { - let n = b.len(); - debug_assert_eq!(a.len(), n * n); +/// A factorised symmetric positive-definite matrix, reusable across queries. +pub(crate) struct Cholesky { + /// Lower triangle of `L`, row-major `n * n`. The upper triangle is + /// leftover scratch from the factorisation and is never read. + l: Vec, + n: usize, +} - // In-place Cholesky: A = L L^T, lower triangle. - for j in 0..n { - let mut d = a[j * n + j]; - for k in 0..j { - d -= a[j * n + k] * a[j * n + k]; - } - // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here too, - // and a negated comparison would let it through as "not positive". - if d.is_nan() || d <= 0.0 { - return None; - } - let d = d.sqrt(); - a[j * n + j] = d; +impl Cholesky { + /// Factorise `a` (row-major, `n * n`, symmetric) into `L L^T`. + /// + /// `a` is consumed as scratch. + /// + /// Returns `None` if the matrix is not positive-definite, which for a + /// precision matrix means the model is improper — a competitor with + /// neither a proper prior nor any evidence. + pub(crate) fn factor(mut a: Vec, n: usize) -> Option { + debug_assert_eq!(a.len(), n * n); - for i in j + 1..n { - let mut s = a[i * n + j]; + for j in 0..n { + let mut d = a[j * n + j]; for k in 0..j { - s -= a[i * n + k] * a[j * n + k]; + d -= a[j * n + k] * a[j * n + k]; + } + // Explicit rather than `!(d > 0.0)`: a NaN pivot must fail here + // too, and a negated comparison would let it through as "not + // positive". + if d.is_nan() || d <= 0.0 { + return None; + } + let d = d.sqrt(); + a[j * n + j] = d; + + for i in j + 1..n { + let mut s = a[i * n + j]; + for k in 0..j { + s -= a[i * n + k] * a[j * n + k]; + } + a[i * n + j] = s / d; } - a[i * n + j] = s / d; } + + Some(Self { l: a, n }) } - // Forward substitution, then back substitution. - let mut z = b.to_vec(); - for i in 0..n { - let mut s = z[i]; - for k in 0..i { - s -= a[i * n + k] * z[k]; + /// Whiten a contrast: `y = L^-1 b`. + /// + /// The point of the result is the dot product, not the vector: for two + /// contrasts `b` and `b'`, `y . y'` is `b^T A^-1 b'`. See the module docs. + pub(crate) fn whiten(&self, b: &[f64]) -> Vec { + debug_assert_eq!(b.len(), self.n); + let n = self.n; + let mut y = b.to_vec(); + for i in 0..n { + // Folded from `y[i]` rather than summed and subtracted once, so the + // accumulation order matches a plain substitution loop exactly. + let row = &self.l[i * n..i * n + i]; + let s = row + .iter() + .zip(&y[..i]) + .fold(y[i], |acc, (l, v)| acc - l * v); + y[i] = s / self.l[i * n + i]; } - z[i] = s / a[i * n + i]; - } - for i in (0..n).rev() { - let mut s = z[i]; - for k in i + 1..n { - s -= a[k * n + i] * z[k]; - } - z[i] = s / a[i * n + i]; + y } +} - Some(z) +/// `b^T A^-1 b'`, given the two whitened contrasts. +pub(crate) fn bilinear(y: &[f64], y_prime: &[f64]) -> f64 { + y.iter().zip(y_prime).map(|(a, b)| a * b).sum() } #[cfg(test)] mod tests { use super::*; + /// `[[4, 1], [1, 3]] z = [1, 2]` has `z = [1/11, 7/11]`, so the quadratic + /// form `b^T A^-1 b` is `1 * 1/11 + 2 * 7/11 = 15/11`. #[test] - fn solves_a_known_system() { - // [[4, 1], [1, 3]] z = [1, 2] => z = [1/11, 7/11] - let a = vec![4.0, 1.0, 1.0, 3.0]; - let z = solve_spd(a, &[1.0, 2.0]).unwrap(); - assert!((z[0] - 1.0 / 11.0).abs() < 1e-12, "{z:?}"); - assert!((z[1] - 7.0 / 11.0).abs() < 1e-12, "{z:?}"); + fn reproduces_a_known_quadratic_form() { + let c = Cholesky::factor(vec![4.0, 1.0, 1.0, 3.0], 2).unwrap(); + let y = c.whiten(&[1.0, 2.0]); + assert!((bilinear(&y, &y) - 15.0 / 11.0).abs() < 1e-12); } + /// Whitening `e_i` recovers the inverse's diagonal, which is the variance + /// of a single variable. #[test] fn recovers_the_inverse_diagonal() { // A = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]; inverse diagonal is // [0.75, 1.0, 0.75]. let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; + let c = Cholesky::factor(a, 3).unwrap(); for (i, expected) in [0.75, 1.0, 0.75].into_iter().enumerate() { let mut e = vec![0.0; 3]; e[i] = 1.0; - let z = solve_spd(a.clone(), &e).unwrap(); - assert!((z[i] - expected).abs() < 1e-12, "row {i}: {z:?}"); + let y = c.whiten(&e); + assert!((bilinear(&y, &y) - expected).abs() < 1e-12, "row {i}"); } } + /// The off-diagonal bilinear form is symmetric and matches the inverse. + #[test] + fn recovers_an_off_diagonal_covariance() { + // Same A; (A^-1)_{0,1} = 0.5. + let a = vec![2.0, -1.0, 0.0, -1.0, 2.0, -1.0, 0.0, -1.0, 2.0]; + let c = Cholesky::factor(a, 3).unwrap(); + let y0 = c.whiten(&[1.0, 0.0, 0.0]); + let y1 = c.whiten(&[0.0, 1.0, 0.0]); + assert!((bilinear(&y0, &y1) - 0.5).abs() < 1e-12); + assert!((bilinear(&y1, &y0) - 0.5).abs() < 1e-12); + } + + /// A variance can never come out negative, because it is a sum of squares. + #[test] + fn a_quadratic_form_is_never_negative() { + let a = vec![1e12, 1e12 - 1.0, 1e12 - 1.0, 1e12]; + let c = Cholesky::factor(a, 2).unwrap(); + let y = c.whiten(&[1.0, -1.0]); + assert!(bilinear(&y, &y) >= 0.0); + } + #[test] fn rejects_a_non_positive_definite_matrix() { // Singular: the second row is a multiple of the first. - let a = vec![1.0, 2.0, 2.0, 4.0]; - assert!(solve_spd(a, &[1.0, 1.0]).is_none()); + assert!(Cholesky::factor(vec![1.0, 2.0, 2.0, 4.0], 2).is_none()); } } diff --git a/src/lib.rs b/src/lib.rs index f58c41b..a918605 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -141,7 +141,7 @@ pub use event::{Event, Member, Team}; pub use event_builder::EventBuilder; pub use game::{Game, GameOptions, OwnedGame}; pub use gaussian::Gaussian; -pub use history::{History, HistoryBuilder}; +pub use history::{History, HistoryBuilder, Joint}; pub use key_table::KeyTable; use matrix::Matrix; pub use observer::{NullObserver, Observer}; diff --git a/tests/joint_handle.rs b/tests/joint_handle.rs new file mode 100644 index 0000000..f93ff77 --- /dev/null +++ b/tests/joint_handle.rs @@ -0,0 +1,220 @@ +//! `History::joint` factorises once and answers many questions. +//! +//! The contract that matters is *identity*: a `Joint` must return exactly what +//! the one-shot call returns, bit for bit. A faster path that quietly disagreed +//! with the slow one would be worse than no fast path — a caller would get +//! different numbers depending on how many questions they happened to ask. + +use smallvec::smallvec; +use trueskill_tt::{ + ConstantDrift, ConvergenceOptions, Event, History, InferenceError, Member, Outcome, Team, + UnknownKeys, +}; + +type H = History; + +fn duel(a: &'static str, b: &'static str, t: i64, sa: f64, sb: f64) -> Event { + Event { + time: t, + teams: smallvec![ + Team::with_members([Member::new(a)]), + Team::with_members([Member::new(b)]), + ], + outcome: Outcome::scores([sa, sb]), + } +} + +fn ranked(a: &'static str, b: &'static str, t: i64) -> Event { + Event { + time: t, + teams: smallvec![ + Team::with_members([Member::new(a)]), + Team::with_members([Member::new(b)]), + ], + outcome: Outcome::winner(0, 2), + } +} + +fn history(unknown: UnknownKeys) -> H { + History::builder() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .score_sigma(2.0) + .drift(ConstantDrift(0.5)) + .unknown_keys(unknown) + .convergence(ConvergenceOptions { + max_iter: 20_000, + epsilon: 1e-13, + alpha: 1.0, + }) + .build() +} + +/// Several slices, competitors with different last appearances, so `latest` +/// and `at_slice` both have work to do. +fn fitted(unknown: UnknownKeys) -> H { + let mut h = history(unknown); + h.add_events(vec![ + duel("a", "b", 1, 5.0, 2.0), + duel("c", "d", 1, 3.0, 3.5), + duel("a", "c", 2, 6.0, 1.0), + duel("b", "d", 3, 4.0, 3.0), + duel("a", "d", 4, 7.0, 2.0), + duel("b", "c", 5, 2.0, 4.0), + ]) + .unwrap(); + let report = h.converge().unwrap(); + assert!(report.converged, "fixture must converge"); + h +} + +const PAIRS: [(&str, &str); 6] = [ + ("a", "b"), + ("a", "c"), + ("a", "d"), + ("b", "c"), + ("b", "d"), + ("c", "d"), +]; + +#[test] +fn a_joint_answers_exactly_what_the_one_shot_call_does() { + let h = fitted(UnknownKeys::Reject); + let joint = h.joint().unwrap(); + + for (a, b) in PAIRS { + let terms = [(&a, 1.0), (&b, -1.0)]; + let one_shot = h.posterior_of(&terms).unwrap(); + let cached = joint.posterior_of(&terms).unwrap(); + assert_eq!(one_shot.pi(), cached.pi(), "{a} - {b}"); + assert_eq!(one_shot.tau(), cached.tau(), "{a} - {b}"); + } +} + +#[test] +fn a_joint_agrees_at_a_pinned_time_too() { + let h = fitted(UnknownKeys::Reject); + let joint = h.joint().unwrap(); + + for time in 1..=5 { + for (a, b) in PAIRS { + let terms = [(&a, 1.0), (&b, -1.0)]; + let one_shot = h.posterior_of_at(time, &terms); + let cached = joint.posterior_of_at(time, &terms); + match (one_shot, cached) { + (Ok(x), Ok(y)) => { + assert_eq!(x.pi(), y.pi(), "t={time} {a} - {b}"); + assert_eq!(x.tau(), y.tau(), "t={time} {a} - {b}"); + } + (Err(x), Err(y)) => assert_eq!(x, y, "t={time} {a} - {b}"), + (x, y) => panic!("t={time} {a} - {b}: disagreed on success: {x:?} vs {y:?}"), + } + } + } +} + +#[test] +fn a_joint_scores_candidate_matchups_identically() { + let h = fitted(UnknownKeys::Reject); + let joint = h.joint().unwrap(); + let (a, b) = ("a", "b"); + let target = [(&a, 1.0), (&b, -1.0)]; + + for (x, y) in PAIRS { + let teams: [&[&&str]; 2] = [&[&x], &[&y]]; + let one_shot = h.expected_variance_reduction(&teams, &target).unwrap(); + let cached = joint.expected_variance_reduction(&teams, &target).unwrap(); + assert_eq!(one_shot, cached, "{x} vs {y}"); + } +} + +/// The whole point: a competitor appears once per slice, so the joint is over +/// appearances rather than competitors, and a caller sizing a batch needs to +/// know which. +#[test] +fn variables_counts_appearances_not_competitors() { + let h = fitted(UnknownKeys::Reject); + let joint = h.joint().unwrap(); + // Four competitors, twelve appearances across five slices, all with + // positive drift between them, so no two collapse. + assert_eq!(joint.variables(), 12); +} + +/// With `drift = 0` consecutive appearances are the same latent variable, so +/// the joint is smaller than the appearance count. +#[test] +fn pinned_competitors_collapse_consecutive_appearances() { + let mut h = History::builder() + .mu(0.0) + .sigma(6.0) + .beta(1.0) + .score_sigma(2.0) + .drift(ConstantDrift(0.0)) + .convergence(ConvergenceOptions { + max_iter: 20_000, + epsilon: 1e-13, + alpha: 1.0, + }) + .build(); + h.add_events(vec![ + duel("a", "b", 1, 5.0, 2.0), + duel("a", "b", 2, 4.0, 3.0), + duel("a", "b", 3, 6.0, 1.0), + ]) + .unwrap(); + assert!(h.converge().unwrap().converged); + assert_eq!(h.joint().unwrap().variables(), 2); +} + +#[test] +fn a_ranked_history_has_no_exact_joint() { + let mut h = history(UnknownKeys::Reject); + h.add_events(vec![duel("a", "b", 1, 5.0, 2.0), ranked("a", "b", 2)]) + .unwrap(); + let _ = h.converge().unwrap(); + assert!(matches!( + h.joint().unwrap_err(), + InferenceError::JointUnavailable { .. } + )); +} + +#[test] +fn an_empty_history_has_no_joint() { + let h = history(UnknownKeys::Reject); + assert!(matches!( + h.joint().unwrap_err(), + InferenceError::JointUnavailable { .. } + )); +} + +/// Unknown keys are decided per query, not when the joint is factorised — the +/// factorisation does not depend on the question. +#[test] +fn unknown_keys_are_rejected_per_query() { + let h = fitted(UnknownKeys::Reject); + let joint = h.joint().unwrap(); + let (a, z) = ("a", "nobody"); + assert!(matches!( + joint.posterior_of(&[(&a, 1.0), (&z, -1.0)]).unwrap_err(), + InferenceError::UnknownKey { .. } + )); + // The handle is still usable afterwards. + let b = "b"; + assert!(joint.posterior_of(&[(&a, 1.0), (&b, -1.0)]).is_ok()); +} + +/// Under `Prior`, an unseen competitor is independent of everything in the +/// history, and the cached path must add the same prior variance the one-shot +/// path does. +#[test] +fn unseen_competitors_match_the_one_shot_path() { + let h = fitted(UnknownKeys::Prior); + let joint = h.joint().unwrap(); + let (a, z) = ("a", "nobody"); + let terms = [(&a, 1.0), (&z, -1.0)]; + let one_shot = h.posterior_of(&terms).unwrap(); + let cached = joint.posterior_of(&terms).unwrap(); + assert_eq!(one_shot.pi(), cached.pi()); + assert_eq!(one_shot.tau(), cached.tau()); +}