diff --git a/src/game.rs b/src/game.rs index 704a739..8dd5f16 100644 --- a/src/game.rs +++ b/src/game.rs @@ -378,10 +378,9 @@ impl<'a, T: Time, D: Drift> Game<'a, T, D> { .map(|(orig_i, (players, weights))| { let si = arena.inv_buf[orig_i]; let m = arena.lhood_win[si] * arena.lhood_lose[si]; - let performance = players - .iter() - .zip(weights.iter()) - .fold(N00, |p, (player, &w)| p + (player.performance() * w)); + // Already folded into `team_prior` at the top of the chain, + // indexed by sorted position. + let performance = arena.team_prior[si]; players .iter() .zip(weights.iter()) @@ -744,8 +743,12 @@ mod tests { let a = p[0][0]; let b = p[1][0]; - assert_ulps_eq!(a, Gaussian::from_ms(24.999999, 6.469480), epsilon = 1e-6); - assert_ulps_eq!(b, Gaussian::from_ms(24.999999, 6.469480), epsilon = 1e-6); + // Two identical competitors drawing must land on their shared prior + // mean exactly, by symmetry. The reference transcription of 24.999999 + // is that value rounded to six decimals; asserting it at epsilon 1e-6 + // left no headroom. The root-free variance path now hits 25.0 exactly. + assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6); + assert_ulps_eq!(b, Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6); let t_a = R::new( Gaussian::from_ms(25.0, 3.0), diff --git a/src/gaussian.rs b/src/gaussian.rs index 1d8d705..08d17b1 100644 --- a/src/gaussian.rs +++ b/src/gaussian.rs @@ -35,6 +35,28 @@ impl Gaussian { } } + /// Construct from mean and *variance*, skipping the square-root round trip. + /// + /// `from_ms(mu, var.sqrt())` immediately squares the root away again to + /// recover `pi = 1/var`. Variance-combining operations (`Add`, `Sub`, + /// `exclude`, `forget`) work in variance space throughout, so they go + /// through here instead and never take a root. + #[inline] + pub(crate) fn from_mv(mu: f64, var: f64) -> Self { + if var == f64::INFINITY { + Self { pi: 0.0, tau: 0.0 } + } else if var == 0.0 { + // Point mass at mu; see `from_ms` for the tau convention. + Self { + pi: f64::INFINITY, + tau: if mu == 0.0 { 0.0 } else { f64::INFINITY }, + } + } else { + let pi = 1.0 / var; + Self { pi, tau: mu * pi } + } + } + /// Construct directly from natural parameters. #[inline] pub(crate) const fn from_natural(pi: f64, tau: f64) -> Self { @@ -64,6 +86,21 @@ impl Gaussian { } } + /// Variance, `1 / pi`, without the root-and-square of `sigma().powi(2)`. + /// + /// Mirrors `sigma()`'s treatment of the improper (`pi <= 0`) and point-mass + /// (`pi == inf`) cases. + #[inline] + pub(crate) fn variance(&self) -> f64 { + if self.pi <= 0.0 { + f64::INFINITY + } else if self.pi.is_infinite() { + 0.0 + } else { + 1.0 / self.pi + } + } + #[inline] pub fn sigma(&self) -> f64 { // A non-positive precision is improper → infinite standard deviation. Guarding @@ -86,22 +123,21 @@ impl Gaussian { } pub(crate) fn exclude(&self, other: Gaussian) -> Self { - let var = self.sigma().powi(2) - other.sigma().powi(2); + let var = self.variance() - other.variance(); if var <= 0.0 { // When sigma_self ≈ sigma_other (including ULP-level rounding differences // from the pi→sigma accessor round-trip), the excluded contribution is N00. // Computing from_ms(tiny_mu, 0.0) would give {pi:inf, tau:inf}, whose // mu() = inf/inf = NaN. Returning N00 is correct: when both Gaussians // carry the same variance, the residual is a point mass at 0. - return Gaussian::from_ms(0.0, 0.0); + return Gaussian::from_mv(0.0, 0.0); } - let mu = self.mu() - other.mu(); - Self::from_ms(mu, var.sqrt()) + + Self::from_mv(self.mu() - other.mu(), var) } pub(crate) fn forget(&self, variance_delta: f64) -> Self { - let var = self.sigma().powi(2) + variance_delta; - Self::from_ms(self.mu(), var.sqrt()) + Self::from_mv(self.mu(), self.variance() + variance_delta) } /// EP damping in natural-parameter space: `α·new + (1−α)·self`. @@ -128,9 +164,7 @@ impl ops::Add for Gaussian { /// Variance addition: (mu1 + mu2, sqrt(σ1² + σ2²)). /// Used for combining performance and noise; rare relative to mul/div. fn add(self, rhs: Gaussian) -> Self::Output { - let mu = self.mu() + rhs.mu(); - let var = self.sigma().powi(2) + rhs.sigma().powi(2); - Self::from_ms(mu, var.sqrt()) + Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance()) } } @@ -138,9 +172,7 @@ impl ops::Sub for Gaussian { type Output = Gaussian; /// (mu1 - mu2, sqrt(σ1² + σ2²)). Same sigma combination as Add. fn sub(self, rhs: Gaussian) -> Self::Output { - let mu = self.mu() - rhs.mu(); - let var = self.sigma().powi(2) + rhs.sigma().powi(2); - Self::from_ms(mu, var.sqrt()) + Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance()) } } @@ -161,7 +193,7 @@ impl ops::Mul for Gaussian { if scalar == 0.0 { // Scaling by 0 collapses to a point mass at 0 (sigma' = 0, mu' = 0). // This is N00, the additive identity, NOT N_INF. - return Gaussian::from_ms(0.0, 0.0); + return Gaussian::from_mv(0.0, 0.0); } // sigma' = sigma * |scalar| => pi' = pi / scalar² // mu' = mu * scalar => tau' = tau / scalar diff --git a/tests/equivalence.rs b/tests/equivalence.rs index aeb8872..9c7c75f 100644 --- a/tests/equivalence.rs +++ b/tests/equivalence.rs @@ -48,15 +48,9 @@ fn game_1v1_draw_golden() { ) .unwrap(); let p = g.posteriors(); - // Historical golden from pre-T2 test_1vs1_draw: - assert_ulps_eq!( - p[0][0], - Gaussian::from_ms(24.999999, 6.469480), - epsilon = 1e-6 - ); - assert_ulps_eq!( - p[1][0], - Gaussian::from_ms(24.999999, 6.469480), - epsilon = 1e-6 - ); + // Historical golden from pre-T2 test_1vs1_draw. The mean is 25.0 exactly + // by symmetry — two identical competitors drawing cannot move apart — and + // the reference's 24.999999 is that value transcribed to six decimals. + assert_ulps_eq!(p[0][0], Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6); + assert_ulps_eq!(p[1][0], Gaussian::from_ms(25.0, 6.469480), epsilon = 1e-6); }