use std::cmp::Ordering;

use crate::{
    N_INF, N00,
    arena::ScratchArena,
    compute_margin,
    drift::Drift,
    factor::{VarId, margin::MarginFactor, trunc::TruncFactor},
    gaussian::Gaussian,
    rating::Rating,
    time::Time,
    tuple_gt, tuple_max,
};

/// Per-adjacent-pair link factor in the game's diff chain.
///
/// `Trunc` is used for `Outcome::Ranked` (rank-based truncation).
/// `Margin` is used for `Outcome::Scored` (Gaussian observation on the diff).
#[derive(Debug)]
pub(crate) enum DiffFactor {
    Trunc(TruncFactor),
    Margin(MarginFactor),
}

impl DiffFactor {
    pub(crate) fn diff(&self) -> VarId {
        match self {
            Self::Trunc(f) => f.diff,
            Self::Margin(f) => f.diff,
        }
    }

    pub(crate) fn msg(&self) -> Gaussian {
        match self {
            Self::Trunc(f) => f.msg,
            Self::Margin(f) => f.msg,
        }
    }

    pub(crate) fn evidence(&self) -> f64 {
        match self {
            Self::Trunc(f) => f.evidence_cached.unwrap_or(1.0),
            Self::Margin(f) => f.evidence_cached.unwrap_or(1.0),
        }
    }

    pub(crate) fn propagate(
        &mut self,
        vars: &mut crate::factor::VarStore,
        alpha: f64,
    ) -> (f64, f64) {
        match self {
            Self::Trunc(f) => f.propagate_with_alpha(vars, alpha),
            Self::Margin(f) => f.propagate_with_alpha(vars, alpha),
        }
    }
}

/// Per-game inference options.
///
/// `p_draw` and `convergence` apply to ranked outcomes (`Game::ranked`).
/// `score_sigma` applies only to scored outcomes (`Game::scored`); it controls
/// how much the engine trusts the observed score margin (smaller σ = more trust).
#[derive(Clone, Copy, Debug)]
pub struct GameOptions {
    pub p_draw: f64,
    pub score_sigma: f64,
    pub convergence: crate::ConvergenceOptions,
}

impl Default for GameOptions {
    fn default() -> Self {
        Self {
            p_draw: crate::P_DRAW,
            score_sigma: 1.0,
            convergence: crate::ConvergenceOptions::default(),
        }
    }
}

/// Owned variant of `Game` returned by public constructors.
///
/// Unlike `Game<'a, T, D>` (which borrows its result/weights slices from
/// History's internal state), `OwnedGame<T, D>` owns its inputs so it can
/// be returned freely from public constructors.
#[derive(Debug)]
#[allow(dead_code)]
pub struct OwnedGame<T: Time, D: Drift<T>> {
    teams: Vec<Vec<Rating<T, D>>>,
    result: Vec<f64>,
    weights: Vec<Vec<f64>>,
    p_draw: f64,
    pub(crate) convergence: crate::ConvergenceOptions,
    pub(crate) likelihoods: Vec<Vec<Gaussian>>,
    pub(crate) evidence: f64,
}

impl<T: Time, D: Drift<T>> OwnedGame<T, D> {
    pub(crate) fn new(
        teams: Vec<Vec<Rating<T, D>>>,
        result: Vec<f64>,
        weights: Vec<Vec<f64>>,
        p_draw: f64,
        convergence: crate::ConvergenceOptions,
    ) -> Self {
        let mut arena = ScratchArena::new();
        let g = Game::ranked_with_arena(
            teams.clone(),
            &result,
            &weights,
            p_draw,
            convergence,
            &mut arena,
        );
        let likelihoods = g.likelihoods;
        let evidence = g.evidence;
        Self {
            teams,
            result,
            weights,
            p_draw,
            convergence,
            likelihoods,
            evidence,
        }
    }

    pub(crate) fn new_scored(
        teams: Vec<Vec<Rating<T, D>>>,
        scores: Vec<f64>,
        weights: Vec<Vec<f64>>,
        score_sigma: f64,
        convergence: crate::ConvergenceOptions,
    ) -> Self {
        let mut arena = ScratchArena::new();
        let g = Game::scored_with_arena(
            teams.clone(),
            &scores,
            &weights,
            score_sigma,
            convergence,
            &mut arena,
        );
        let likelihoods = g.likelihoods;
        let evidence = g.evidence;
        Self {
            teams,
            result: scores,
            weights,
            p_draw: 0.0,
            convergence,
            likelihoods,
            evidence,
        }
    }

    pub fn posteriors(&self) -> Vec<Vec<Gaussian>> {
        self.likelihoods
            .iter()
            .zip(self.teams.iter())
            .map(|(l, t)| l.iter().zip(t.iter()).map(|(&l, r)| l * r.prior).collect())
            .collect()
    }

    pub fn log_evidence(&self) -> f64 {
        self.evidence.ln()
    }
}

#[derive(Debug)]
pub struct Game<'a, T: Time = i64, D: Drift<T> = crate::drift::ConstantDrift> {
    teams: Vec<Vec<Rating<T, D>>>,
    result: &'a [f64],
    weights: &'a [Vec<f64>],
    p_draw: f64,
    pub(crate) convergence: crate::ConvergenceOptions,
    pub(crate) likelihoods: Vec<Vec<Gaussian>>,
    pub(crate) evidence: f64,
}

impl<'a, T: Time, D: Drift<T>> Game<'a, T, D> {
    pub(crate) fn ranked_with_arena(
        teams: Vec<Vec<Rating<T, D>>>,
        result: &'a [f64],
        weights: &'a [Vec<f64>],
        p_draw: f64,
        convergence: crate::ConvergenceOptions,
        arena: &mut ScratchArena,
    ) -> Self {
        debug_assert!(
            result.len() == teams.len(),
            "result must have the same length as teams"
        );
        debug_assert!(
            weights
                .iter()
                .zip(teams.iter())
                .all(|(w, t)| w.len() == t.len()),
            "weights must have the same dimensions as teams"
        );
        debug_assert!(
            (0.0..1.0).contains(&p_draw),
            "draw probability must be >= 0.0 and < 1.0"
        );
        debug_assert!(
            p_draw > 0.0 || {
                let mut r = result.to_vec();
                r.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap());
                r.windows(2).all(|w| w[0] != w[1])
            },
            "draw must be > 0.0 if there are teams with draw"
        );
        debug_assert!(
            convergence.alpha > 0.0 && convergence.alpha <= 1.0,
            "convergence alpha must be in (0.0, 1.0]"
        );

        let mut this = Self {
            teams,
            result,
            weights,
            p_draw,
            convergence,
            likelihoods: Vec::new(),
            evidence: 0.0,
        };

        this.likelihoods(arena);
        this
    }

    pub(crate) fn scored_with_arena(
        teams: Vec<Vec<Rating<T, D>>>,
        scores: &'a [f64],
        weights: &'a [Vec<f64>],
        score_sigma: f64,
        convergence: crate::ConvergenceOptions,
        arena: &mut ScratchArena,
    ) -> Self {
        debug_assert!(
            scores.len() == teams.len(),
            "scores must have the same length as teams"
        );
        debug_assert!(
            weights
                .iter()
                .zip(teams.iter())
                .all(|(w, t)| w.len() == t.len()),
            "weights must have the same dimensions as teams"
        );
        debug_assert!(score_sigma > 0.0, "score_sigma must be positive");
        debug_assert!(
            convergence.alpha > 0.0 && convergence.alpha <= 1.0,
            "convergence alpha must be in (0.0, 1.0]"
        );

        let mut this = Self {
            teams,
            result: scores,
            weights,
            p_draw: 0.0,
            convergence,
            likelihoods: Vec::new(),
            evidence: 0.0,
        };

        this.likelihoods_scored(arena, score_sigma);
        this
    }

    fn run_chain<F>(&self, arena: &mut ScratchArena, mut make_link: F) -> (f64, Vec<Vec<Gaussian>>)
    where
        F: FnMut(usize, &[usize], &mut crate::factor::VarStore) -> DiffFactor,
    {
        arena.reset();

        let alpha = self.convergence.alpha;
        let epsilon = self.convergence.epsilon;
        let max_iter = self.convergence.max_iter;

        let n_teams = self.teams.len();

        arena.sort_buf.extend(0..n_teams);
        arena.sort_buf.sort_by(|&i, &j| {
            self.result[j]
                .partial_cmp(&self.result[i])
                .unwrap_or(Ordering::Equal)
        });

        arena.team_prior.extend(arena.sort_buf.iter().map(|&t| {
            self.teams[t]
                .iter()
                .zip(self.weights[t].iter())
                .fold(N00, |p, (player, &w)| p + (player.performance() * w))
        }));

        let n_diffs = n_teams.saturating_sub(1);

        let mut links: Vec<DiffFactor> = (0..n_diffs)
            .map(|i| make_link(i, &arena.sort_buf, &mut arena.vars))
            .collect();

        arena.lhood_lose.resize(n_teams, N_INF);
        arena.lhood_win.resize(n_teams, N_INF);

        let mut step = (f64::INFINITY, f64::INFINITY);
        let mut iter = 0;

        while tuple_gt(step, epsilon) && iter < max_iter {
            step = (0.0_f64, 0.0_f64);

            for (e, lf) in links[..n_diffs.saturating_sub(1)].iter_mut().enumerate() {
                let pw = arena.team_prior[e] * arena.lhood_lose[e];
                let pl = arena.team_prior[e + 1] * arena.lhood_win[e + 1];
                let raw = pw - pl;
                arena.vars.set(lf.diff(), raw * lf.msg());
                let d = lf.propagate(&mut arena.vars, alpha);
                step = tuple_max(step, d);

                let new_ll = pw - lf.msg();
                step = tuple_max(step, arena.lhood_lose[e + 1].delta(new_ll));
                arena.lhood_lose[e + 1] = new_ll;
            }

            for (rev_i, lf) in links[1..].iter_mut().rev().enumerate() {
                let e = n_diffs - 1 - rev_i;
                let pw = arena.team_prior[e] * arena.lhood_lose[e];
                let pl = arena.team_prior[e + 1] * arena.lhood_win[e + 1];
                let raw = pw - pl;
                arena.vars.set(lf.diff(), raw * lf.msg());
                let d = lf.propagate(&mut arena.vars, alpha);
                step = tuple_max(step, d);

                let new_lw = pl + lf.msg();
                step = tuple_max(step, arena.lhood_win[e].delta(new_lw));
                arena.lhood_win[e] = new_lw;
            }

            iter += 1;
        }

        // Special case: exactly 1 diff (2-team game); loop body was empty.
        if n_diffs == 1 {
            let raw = (arena.team_prior[0] * arena.lhood_lose[0])
                - (arena.team_prior[1] * arena.lhood_win[1]);
            arena.vars.set(links[0].diff(), raw * links[0].msg());
            links[0].propagate(&mut arena.vars, alpha);
        }

        // Boundary updates: close the chain at both ends.
        if n_diffs > 0 {
            let pl1 = arena.team_prior[1] * arena.lhood_win[1];
            arena.lhood_win[0] = pl1 + links[0].msg();
            let pw_last = arena.team_prior[n_teams - 2] * arena.lhood_lose[n_teams - 2];
            arena.lhood_lose[n_teams - 1] = pw_last - links[n_diffs - 1].msg();
        }

        let evidence: f64 = links.iter().map(|l| l.evidence()).product();

        // Inverse permutation: inv_buf[orig_i] = sorted_i.
        arena.inv_buf.resize(n_teams, 0);
        for (si, &orig_i) in arena.sort_buf.iter().enumerate() {
            arena.inv_buf[orig_i] = si;
        }

        let likelihoods = self
            .teams
            .iter()
            .zip(self.weights.iter())
            .enumerate()
            .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));
                players
                    .iter()
                    .zip(weights.iter())
                    .map(|(player, &w)| {
                        ((m - performance.exclude(player.performance() * w)) * (1.0 / w))
                            .forget(player.beta.powi(2))
                    })
                    .collect::<Vec<_>>()
            })
            .collect::<Vec<_>>();

        (evidence, likelihoods)
    }

    fn likelihoods(&mut self, arena: &mut ScratchArena) {
        let (evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
            let tie = self.result[sort_buf[i]] == self.result[sort_buf[i + 1]];
            let margin = if self.p_draw == 0.0 {
                0.0
            } else {
                let a: f64 = self.teams[sort_buf[i]].iter().map(|p| p.beta.powi(2)).sum();
                let b: f64 = self.teams[sort_buf[i + 1]]
                    .iter()
                    .map(|p| p.beta.powi(2))
                    .sum();
                compute_margin(self.p_draw, (a + b).sqrt())
            };
            let vid = vars.alloc(N_INF);
            DiffFactor::Trunc(TruncFactor::new(vid, margin, tie))
        });
        self.evidence = evidence;
        self.likelihoods = likelihoods;
    }

    fn likelihoods_scored(&mut self, arena: &mut ScratchArena, score_sigma: f64) {
        let (evidence, likelihoods) = self.run_chain(arena, |i, sort_buf, vars| {
            let m_obs = self.result[sort_buf[i]] - self.result[sort_buf[i + 1]];
            let vid = vars.alloc(N_INF);
            DiffFactor::Margin(MarginFactor::new(vid, m_obs, score_sigma))
        });
        self.evidence = evidence;
        self.likelihoods = likelihoods;
    }

    pub fn posteriors(&self) -> Vec<Vec<Gaussian>> {
        self.likelihoods
            .iter()
            .zip(self.teams.iter())
            .map(|(l, t)| {
                l.iter()
                    .zip(t.iter())
                    .map(|(&l, p)| l * p.prior)
                    .collect::<Vec<_>>()
            })
            .collect::<Vec<_>>()
    }

    pub fn log_evidence(&self) -> f64 {
        self.evidence.ln()
    }
}

impl<T: Time, D: Drift<T>> Game<'_, T, D> {
    pub fn ranked(
        teams: &[&[Rating<T, D>]],
        outcome: crate::Outcome,
        options: &GameOptions,
    ) -> Result<OwnedGame<T, D>, crate::InferenceError> {
        if !(0.0..1.0).contains(&options.p_draw) {
            return Err(crate::InferenceError::InvalidProbability {
                value: options.p_draw,
            });
        }
        if outcome.team_count() != teams.len() {
            return Err(crate::InferenceError::MismatchedShape {
                kind: "outcome ranks vs teams",
                expected: teams.len(),
                got: outcome.team_count(),
            });
        }

        let ranks = outcome
            .as_ranks()
            .ok_or(crate::InferenceError::MismatchedShape {
                kind: "Game::ranked requires Outcome::Ranked",
                expected: 0,
                got: 0,
            })?;
        let max_rank = ranks.iter().copied().max().unwrap_or(0) as f64;
        let result: Vec<f64> = ranks.iter().map(|&r| max_rank - r as f64).collect();
        let teams_owned: Vec<Vec<Rating<T, D>>> = teams.iter().map(|t| t.to_vec()).collect();
        let weights: Vec<Vec<f64>> = teams.iter().map(|t| vec![1.0; t.len()]).collect();

        Ok(OwnedGame::new(
            teams_owned,
            result,
            weights,
            options.p_draw,
            options.convergence,
        ))
    }

    pub fn scored(
        teams: &[&[Rating<T, D>]],
        outcome: crate::Outcome,
        options: &GameOptions,
    ) -> Result<OwnedGame<T, D>, crate::InferenceError> {
        if options.score_sigma <= 0.0 || options.score_sigma.is_nan() {
            return Err(crate::InferenceError::InvalidParameter {
                name: "score_sigma",
                value: options.score_sigma,
            });
        }
        if outcome.team_count() != teams.len() {
            return Err(crate::InferenceError::MismatchedShape {
                kind: "outcome scores vs teams",
                expected: teams.len(),
                got: outcome.team_count(),
            });
        }
        let scores = outcome
            .as_scores()
            .ok_or(crate::InferenceError::MismatchedShape {
                kind: "Game::scored requires Outcome::Scored",
                expected: 0,
                got: 0,
            })?
            .to_vec();
        let teams_owned: Vec<Vec<Rating<T, D>>> = teams.iter().map(|t| t.to_vec()).collect();
        let weights: Vec<Vec<f64>> = teams.iter().map(|t| vec![1.0; t.len()]).collect();
        Ok(OwnedGame::new_scored(
            teams_owned,
            scores,
            weights,
            options.score_sigma,
            options.convergence,
        ))
    }

    pub fn one_v_one(
        a: &Rating<T, D>,
        b: &Rating<T, D>,
        outcome: crate::Outcome,
    ) -> Result<(Gaussian, Gaussian), crate::InferenceError> {
        let game = Self::ranked(&[&[*a], &[*b]], outcome, &GameOptions::default())?;
        let post = game.posteriors();
        Ok((post[0][0], post[1][0]))
    }

    pub fn free_for_all(
        players: &[&Rating<T, D>],
        outcome: crate::Outcome,
        options: &GameOptions,
    ) -> Result<OwnedGame<T, D>, crate::InferenceError> {
        let teams: Vec<Vec<Rating<T, D>>> = players.iter().map(|p| vec![**p]).collect();
        let team_refs: Vec<&[Rating<T, D>]> = teams.iter().map(|t| t.as_slice()).collect();
        Self::ranked(&team_refs, outcome, options)
    }

    #[doc(hidden)]
    pub fn custom<S: crate::factors::Schedule>(
        factors: &mut [crate::factors::BuiltinFactor],
        vars: &mut crate::factors::VarStore,
        schedule: &S,
    ) -> crate::factors::ScheduleReport {
        schedule.run(factors, vars)
    }
}

#[cfg(test)]
mod tests {
    use ::approx::assert_ulps_eq;

    use super::*;
    use crate::{ConstantDrift, GAMMA, Gaussian, N_INF, Rating, arena::ScratchArena};

    type R = Rating<i64, ConstantDrift>;

    #[test]
    fn test_1vs1() {
        let t_a = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_b = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );

        let w = [vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b]],
            &[0.0, 1.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];

        assert_ulps_eq!(a, Gaussian::from_ms(20.794779, 7.194481), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(29.205220, 7.194481), epsilon = 1e-6);

        let t_a = R::new(
            Gaussian::from_ms(29.0, 1.0),
            25.0 / 6.0,
            ConstantDrift(GAMMA),
        );
        let t_b = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(GAMMA),
        );

        let w = [vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b]],
            &[0.0, 1.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];

        assert_ulps_eq!(a, Gaussian::from_ms(28.896475, 0.996604), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(32.189211, 6.062063), epsilon = 1e-6);

        let t_a = R::new(Gaussian::from_ms(1.139, 0.531), 1.0, ConstantDrift(0.2125));
        let t_b = R::new(Gaussian::from_ms(15.568, 0.51), 1.0, ConstantDrift(0.2125));

        let w = [vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b]],
            &[0.0, 1.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );

        assert_eq!(g.likelihoods[0][0], N_INF);
        assert_eq!(g.likelihoods[1][0], N_INF);
    }

    #[test]
    fn test_1vs1vs1() {
        let teams = vec![
            vec![R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            )],
            vec![R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            )],
            vec![R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            )],
        ];

        let w = [vec![1.0], vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            teams.clone(),
            &[1.0, 2.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];

        assert_ulps_eq!(a, Gaussian::from_ms(25.000000, 6.238469), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(31.311358, 6.698818), epsilon = 1e-6);

        let w = [vec![1.0], vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            teams.clone(),
            &[2.0, 1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];

        assert_ulps_eq!(a, Gaussian::from_ms(31.311358, 6.698818), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(25.000000, 6.238469), epsilon = 1e-6);

        let w = [vec![1.0], vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            teams,
            &[1.0, 2.0, 0.0],
            &w,
            0.5,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];
        let c = p[2][0];

        // T1 ULP shift: mu rounds to 25.0 (was 24.999999) under natural-parameter storage.
        assert_ulps_eq!(a, Gaussian::from_ms(25.0, 6.092561), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(33.379314, 6.483575), epsilon = 1e-6);
        assert_ulps_eq!(c, Gaussian::from_ms(16.620685, 6.483575), epsilon = 1e-6);
    }

    #[test]
    fn test_1vs1_draw() {
        let t_a = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_b = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );

        let w = [vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b]],
            &[0.0, 0.0],
            &w,
            0.25,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        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);

        let t_a = R::new(
            Gaussian::from_ms(25.0, 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_b = R::new(
            Gaussian::from_ms(29.0, 2.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );

        let w = [vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b]],
            &[0.0, 0.0],
            &w,
            0.25,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];

        assert_ulps_eq!(a, Gaussian::from_ms(25.736001, 2.709956), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(28.672888, 1.916471), epsilon = 1e-6);
    }

    #[test]
    fn test_1vs1vs1_draw() {
        let t_a = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_b = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_c = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );

        let w = [vec![1.0], vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b], vec![t_c]],
            &[0.0, 0.0, 0.0],
            &w,
            0.25,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];
        let c = p[2][0];

        // Goldens updated for natural-parameter storage: mu rounds to 25.0 (was 24.999999),
        // sigma shifts by ~3e-7 ULPs (within 1e-6 of original). Both bounded differences.
        assert_ulps_eq!(a, Gaussian::from_ms(25.0, 5.729069), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(25.0, 5.707424), epsilon = 1e-6);
        assert_ulps_eq!(c, Gaussian::from_ms(25.0, 5.729069), epsilon = 1e-6);

        let t_a = R::new(
            Gaussian::from_ms(25.0, 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_b = R::new(
            Gaussian::from_ms(25.0, 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let t_c = R::new(
            Gaussian::from_ms(29.0, 2.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );

        let w = [vec![1.0], vec![1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![vec![t_a], vec![t_b], vec![t_c]],
            &[0.0, 0.0, 0.0],
            &w,
            0.25,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        let a = p[0][0];
        let b = p[1][0];
        let c = p[2][0];

        assert_ulps_eq!(a, Gaussian::from_ms(25.488507, 2.638208), epsilon = 1e-6);
        assert_ulps_eq!(b, Gaussian::from_ms(25.510671, 2.628751), epsilon = 1e-6);
        assert_ulps_eq!(c, Gaussian::from_ms(28.555920, 1.885689), epsilon = 1e-6);
    }

    #[test]
    fn test_2vs1vs2_mixed() {
        let t_a = vec![
            R::new(
                Gaussian::from_ms(12.0, 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            ),
            R::new(
                Gaussian::from_ms(18.0, 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            ),
        ];
        let t_b = vec![R::new(
            Gaussian::from_ms(30.0, 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        )];
        let t_c = vec![
            R::new(
                Gaussian::from_ms(14.0, 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            ),
            R::new(
                Gaussian::from_ms(16., 3.0),
                25.0 / 6.0,
                ConstantDrift(25.0 / 300.0),
            ),
        ];

        let w = [vec![1.0, 1.0], vec![1.0], vec![1.0, 1.0]];
        let g = Game::ranked_with_arena(
            vec![t_a, t_b, t_c],
            &[1.0, 0.0, 0.0],
            &w,
            0.25,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(p[0][0], Gaussian::from_ms(13.051, 2.864), epsilon = 1e-3);
        assert_ulps_eq!(p[0][1], Gaussian::from_ms(19.051, 2.864), epsilon = 1e-3);
        assert_ulps_eq!(p[1][0], Gaussian::from_ms(29.292, 2.764), epsilon = 1e-3);
        assert_ulps_eq!(p[2][0], Gaussian::from_ms(13.658, 2.813), epsilon = 1e-3);
        assert_ulps_eq!(p[2][1], Gaussian::from_ms(15.658, 2.813), epsilon = 1e-3);
    }

    #[test]
    fn test_1vs1_weighted() {
        let w_a = vec![1.0];
        let w_b = vec![2.0];

        let t_a = vec![R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(0.0),
        )];
        let t_b = vec![R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(0.0),
        )];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a.clone(), t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(30.625173, 7.765472),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(13.749653, 5.733840),
            epsilon = 1e-6
        );

        let w_a = vec![1.0];
        let w_b = vec![0.7];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a.clone(), t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(27.630080, 7.206676),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(23.158943, 7.801628),
            epsilon = 1e-6
        );

        let w_a = vec![1.6];
        let w_b = vec![0.7];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a, t_b],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(26.142438, 7.573088),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(24.500183, 8.193278),
            epsilon = 1e-6
        );

        let w_a = vec![1.0];
        let w_b = vec![0.0];

        let t_a = vec![R::new(Gaussian::from_ms(2.0, 6.0), 1.0, ConstantDrift(0.0))];
        let t_b = vec![R::new(Gaussian::from_ms(2.0, 6.0), 1.0, ConstantDrift(0.0))];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a, t_b],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(5.557067, 4.052826),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(2.000000, 6.000000),
            epsilon = 1e-6
        );

        let w_a = vec![1.0];
        let w_b = vec![-1.0];

        let t_a = vec![R::new(Gaussian::from_ms(2.0, 6.0), 1.0, ConstantDrift(0.0))];
        let t_b = vec![R::new(Gaussian::from_ms(2.0, 6.0), 1.0, ConstantDrift(0.0))];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a, t_b],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(p[0][0], p[1][0], epsilon = 1e-6);
    }

    #[test]
    fn diff_factor_dispatch_trunc_and_margin() {
        use super::DiffFactor;
        use crate::factor::{VarStore, margin::MarginFactor, trunc::TruncFactor};

        let mut vars = VarStore::new();
        let dt = vars.alloc(Gaussian::from_ms(0.0, 6.0));
        let dm = vars.alloc(Gaussian::from_ms(0.0, 6.0));

        let mut t = DiffFactor::Trunc(TruncFactor::new(dt, 0.0, false));
        let mut m = DiffFactor::Margin(MarginFactor::new(dm, 5.0, 1.0));

        let _ = t.propagate(&mut vars, 1.0);
        let _ = m.propagate(&mut vars, 1.0);

        // Smoke: both diffs got written; their msgs are non-N_INF.
        assert!(t.msg().pi() > 0.0);
        assert!(m.msg().pi() > 0.0);
        assert_eq!(t.diff(), dt);
        assert_eq!(m.diff(), dm);
    }

    #[test]
    fn scored_path_sharper_when_margin_is_large() {
        let prior = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let teams = vec![vec![prior], vec![prior]];
        let result = vec![10.0, 0.0]; // a beat b by 10
        let weights = [vec![1.0], vec![1.0]];
        let mut arena = ScratchArena::new();
        let g = Game::scored_with_arena(
            teams,
            &result,
            &weights,
            1.0,
            crate::ConvergenceOptions::default(),
            &mut arena,
        );
        let p = g.posteriors();
        let a = p[0][0];
        let b = p[1][0];
        assert!(
            a.mu() > b.mu(),
            "expected team a posterior mu > team b; got {} vs {}",
            a.mu(),
            b.mu()
        );

        // Tighter score_sigma should produce a stronger update.
        let mut arena2 = ScratchArena::new();
        let g_tight = Game::scored_with_arena(
            vec![vec![prior], vec![prior]],
            &result,
            &weights,
            0.1,
            crate::ConvergenceOptions::default(),
            &mut arena2,
        );
        let p_tight = g_tight.posteriors();
        let a_tight = p_tight[0][0];
        assert!(
            a_tight.mu() > a.mu(),
            "expected tighter sigma to push posterior further; {} vs {}",
            a_tight.mu(),
            a.mu()
        );
    }

    #[test]
    fn game_scored_public_ctor() {
        use crate::Outcome;
        let prior = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let opts = GameOptions {
            score_sigma: 1.0,
            ..GameOptions::default()
        };
        let g = Game::scored(&[&[prior], &[prior]], Outcome::scores([8.0, 2.0]), &opts).unwrap();
        let p = g.posteriors();
        assert!(p[0][0].mu() > p[1][0].mu());
    }

    #[test]
    fn game_scored_rejects_ranked_outcome() {
        let prior = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let err = Game::scored(
            &[&[prior], &[prior]],
            crate::Outcome::winner(0, 2),
            &GameOptions::default(),
        )
        .unwrap_err();
        assert!(matches!(err, crate::InferenceError::MismatchedShape { .. }));
    }

    #[test]
    fn game_scored_rejects_zero_score_sigma() {
        let prior = R::new(
            Gaussian::from_ms(25.0, 25.0 / 3.0),
            25.0 / 6.0,
            ConstantDrift(25.0 / 300.0),
        );
        let opts = GameOptions {
            score_sigma: 0.0,
            ..GameOptions::default()
        };
        let err = Game::scored(
            &[&[prior], &[prior]],
            crate::Outcome::scores([1.0, 0.0]),
            &opts,
        )
        .unwrap_err();
        assert!(matches!(
            err,
            crate::InferenceError::InvalidParameter {
                name: "score_sigma",
                ..
            }
        ));
    }

    #[test]
    fn test_2vs2_weighted() {
        let t_a = vec![
            R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(0.0),
            ),
            R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(0.0),
            ),
        ];
        let w_a = vec![0.4, 0.8];

        let t_b = vec![
            R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(0.0),
            ),
            R::new(
                Gaussian::from_ms(25.0, 25.0 / 3.0),
                25.0 / 6.0,
                ConstantDrift(0.0),
            ),
        ];
        let w_b = vec![0.9, 0.6];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a.clone(), t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(27.539023, 8.129639),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[0][1],
            Gaussian::from_ms(30.078046, 7.485372),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(19.287197, 7.243465),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][1],
            Gaussian::from_ms(21.191465, 7.867608),
            epsilon = 1e-6
        );

        let w_a = vec![1.3, 1.5];
        let w_b = vec![0.7, 0.4];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a.clone(), t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(25.190190, 8.220511),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[0][1],
            Gaussian::from_ms(25.219450, 8.182783),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(24.897589, 8.300779),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][1],
            Gaussian::from_ms(24.941479, 8.322717),
            epsilon = 1e-6
        );

        let w_a = vec![1.6, 0.2];
        let w_b = vec![0.7, 2.4];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a.clone(), t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(
            p[0][0],
            Gaussian::from_ms(31.674697, 7.501180),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[0][1],
            Gaussian::from_ms(25.834337, 8.320970),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][0],
            Gaussian::from_ms(22.079819, 8.180607),
            epsilon = 1e-6
        );
        assert_ulps_eq!(
            p[1][1],
            Gaussian::from_ms(14.987953, 6.308469),
            epsilon = 1e-6
        );

        let w = [vec![1.0, 1.0], vec![1.0]];
        let g = Game::ranked_with_arena(
            vec![
                t_a.clone(),
                vec![R::new(
                    Gaussian::from_ms(25.0, 25.0 / 3.0),
                    25.0 / 6.0,
                    ConstantDrift(0.0),
                )],
            ],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let post_2vs1 = g.posteriors();

        let w_a = vec![1.0, 1.0];
        let w_b = vec![1.0, 0.0];

        let w = [w_a, w_b];
        let g = Game::ranked_with_arena(
            vec![t_a, t_b.clone()],
            &[1.0, 0.0],
            &w,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut ScratchArena::new(),
        );
        let p = g.posteriors();

        assert_ulps_eq!(p[0][0], post_2vs1[0][0], epsilon = 1e-6);
        assert_ulps_eq!(p[0][1], post_2vs1[0][1], epsilon = 1e-6);
        assert_ulps_eq!(p[1][0], post_2vs1[1][0], epsilon = 1e-6);
        assert_ulps_eq!(p[1][1], t_b[1].prior, epsilon = 1e-6);
    }

    #[test]
    fn run_chain_honours_max_iter_in_convergence_options() {
        let players: Vec<R> = (0..4).map(|_| R::default()).collect();
        let teams: Vec<Vec<_>> = players.iter().map(|p| vec![*p]).collect();
        let result = vec![3.0, 2.0, 1.0, 0.0];
        let weights = vec![vec![1.0]; 4];

        // Capped at 1 iteration: cannot fully propagate down a 4-team chain.
        let mut arena = ScratchArena::new();
        let g_capped = Game::ranked_with_arena(
            teams.clone(),
            &result,
            &weights,
            0.0,
            crate::ConvergenceOptions {
                max_iter: 1,
                ..crate::ConvergenceOptions::default()
            },
            &mut arena,
        );
        let posteriors_capped = g_capped.posteriors();

        // Same inputs, plenty of iterations: fully converged.
        let mut arena = ScratchArena::new();
        let g_full = Game::ranked_with_arena(
            teams,
            &result,
            &weights,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut arena,
        );
        let posteriors_full = g_full.posteriors();

        // The two posteriors should differ — capped did not converge.
        let mut max_diff: f64 = 0.0;
        for (team_capped, team_full) in posteriors_capped.iter().zip(posteriors_full.iter()) {
            for (g_capped, g_full) in team_capped.iter().zip(team_full.iter()) {
                max_diff = max_diff.max((g_capped.mu() - g_full.mu()).abs());
                max_diff = max_diff.max((g_capped.sigma() - g_full.sigma()).abs());
            }
        }
        assert!(
            max_diff > 1e-6,
            "max_iter=1 should differ from full convergence; max_diff={max_diff}"
        );
    }

    #[test]
    fn run_chain_with_damping_converges_to_same_posterior() {
        let players: Vec<R> = (0..4).map(|_| R::default()).collect();
        let teams: Vec<Vec<_>> = players.iter().map(|p| vec![*p]).collect();
        let result = vec![3.0, 2.0, 1.0, 0.0];
        let weights = vec![vec![1.0]; 4];

        let mut arena = ScratchArena::new();
        let g_undamped = Game::ranked_with_arena(
            teams.clone(),
            &result,
            &weights,
            0.0,
            crate::ConvergenceOptions::default(),
            &mut arena,
        );
        let posteriors_undamped = g_undamped.posteriors();

        // alpha=0.5 with extra iterations: should reach the same fixed point.
        let mut arena = ScratchArena::new();
        let g_damped = Game::ranked_with_arena(
            teams,
            &result,
            &weights,
            0.0,
            crate::ConvergenceOptions {
                alpha: 0.5,
                max_iter: 100,
                ..crate::ConvergenceOptions::default()
            },
            &mut arena,
        );
        let posteriors_damped = g_damped.posteriors();

        let mut max_diff: f64 = 0.0;
        for (team_u, team_d) in posteriors_undamped.iter().zip(posteriors_damped.iter()) {
            for (g_u, g_d) in team_u.iter().zip(team_d.iter()) {
                max_diff = max_diff.max((g_u.mu() - g_d.mu()).abs());
                max_diff = max_diff.max((g_u.sigma() - g_d.sigma()).abs());
            }
        }
        assert!(
            max_diff < 1e-4,
            "α=0.5 should reach the same fixed point as α=1.0; max_diff={max_diff}"
        );
    }
}
