//! `TrueSkill` Through Time — Bayesian skill rating over a time axis.
//!
//! Where plain `TrueSkill` gives each competitor one running estimate, `TrueSkill`
//! Through Time treats a whole history as a single model and infers skill *at
//! every point in time*. Evidence flows both directions: a result today
//! sharpens the estimate of who someone was last year, so early estimates stop
//! being frozen guesses and comparisons across eras become meaningful.
//!
//! This is a Rust port of
//! [TrueSkillThroughTime.py](https://github.com/glandfried/TrueSkillThroughTime.py).
//!
//! # Getting started
//!
//! Record results, converge, then read off skills:
//!
//! ```
//! use trueskill_tt::History;
//!
//! let mut history = History::default();
//!
//! history.record_winner(&"alice", &"bob", 1)?;
//! history.record_winner(&"bob", &"carol", 2)?;
//! history.record_winner(&"alice", &"carol", 3)?;
//!
//! let report = history.converge()?;
//! assert!(report.converged);
//!
//! let alice = history.current_skill("alice").unwrap();
//! assert!(alice.mu() > 0.0, "alice won every game she played");
//! # Ok::<(), trueskill_tt::InferenceError>(())
//! ```
//!
//! Teams, weights, explicit rankings and continuous scores go through the
//! fluent event builder:
//!
//! ```
//! use trueskill_tt::History;
//!
//! let mut history = History::builder().p_draw(0.1).build();
//!
//! history
//!     .event(1)
//!     .team(["alice", "bob"])
//!     .team(["carol", "dave"])
//!     .ranking([0, 1])
//!     .commit()?;
//!
//! history.converge()?;
//! # Ok::<(), trueskill_tt::InferenceError>(())
//! ```
//!
//! # Draws need a draw probability
//!
//! A `p_draw` of zero asserts that draws cannot happen, so a tied result has
//! no representable likelihood and is rejected:
//!
//! ```
//! use trueskill_tt::{History, InferenceError};
//!
//! let mut history = History::default(); // p_draw defaults to 0.0
//! let err = history.record_draw(&"alice", &"bob", 1).unwrap_err();
//! assert!(matches!(err, InferenceError::TieWithoutDrawProbability { .. }));
//! ```
//!
//! This also applies to [`Outcome::winner`] for three or more teams, which
//! ties every loser. Configure a positive `p_draw` for those.
//!
//! # Core types
//!
//! - [`History`] — the top-level container: ingests events, runs
//!   forward/backward message passing, and answers queries.
//! - [`Gaussian`] — the probability type, stored in natural parameters
//!   (`pi = 1/sigma²`, `tau = mu/sigma²`) so message passing is add/subtract.
//! - [`Game`] — one match in isolation, for scoring a hypothetical without a
//!   history.
//! - [`Outcome`] — how a match ended: ranks, or continuous scores.
//! - [`Rating`] — a competitor's static configuration (prior, `beta`, drift).
//!
//! # Feature flags
//!
//! - `approx` — implements [`approx`](https://docs.rs/approx) equality traits
//!   for [`Gaussian`]. Useful in tests.
//! - `rayon` — parallelises the within-slice sweep and the per-slice passes of
//!   `learning_curves`/`log_evidence`. Opt-in; results stay bit-identical
//!   regardless of worker count.

#![forbid(unsafe_code)]

use std::{
    cmp::Reverse,
    f64::consts::{FRAC_1_SQRT_2, FRAC_2_SQRT_PI, SQRT_2},
};

#[cfg(feature = "approx")]
mod approx;
pub(crate) mod arena;
mod time;
mod time_slice;
pub use time_slice::{EventKind, TimeSlice};
mod color_group;
mod competitor;
mod convergence;
pub mod drift;
mod error;
mod event;
mod event_builder;
pub(crate) mod factor;
pub mod factors;
mod game;
pub mod gaussian;
mod history;
mod key_table;
mod matrix;
mod observer;
mod outcome;
mod rating;
pub(crate) mod schedule;
pub mod storage;

pub use competitor::Competitor;
pub use convergence::{ConvergenceOptions, ConvergenceReport};
pub use drift::{ConstantDrift, Drift};
pub use error::InferenceError;
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 key_table::KeyTable;
use matrix::Matrix;
pub use observer::{NullObserver, Observer};
pub use outcome::Outcome;
pub use rating::Rating;
pub use schedule::ScheduleReport;
pub use time::{Time, Untimed};

pub const BETA: f64 = 1.0;
pub const MU: f64 = 0.0;
pub const SIGMA: f64 = BETA * 6.0;
pub const GAMMA: f64 = BETA * 0.03;
pub const P_DRAW: f64 = 0.0;
pub const EPSILON: f64 = 1e-6;
pub const ITERATIONS: usize = 30;

const SQRT_TAU: f64 = 2.5066282746310002;

pub const N01: Gaussian = Gaussian::from_ms(0.0, 1.0);
pub const N00: Gaussian = Gaussian::from_ms(0.0, 0.0);
pub const N_INF: Gaussian = Gaussian::from_ms(0.0, f64::INFINITY);

#[derive(Copy, Clone, Default, PartialEq, PartialOrd, Eq, Ord, Hash, Debug)]
pub struct Index(usize);

impl Index {
    /// The underlying slot number.
    ///
    /// Indices are dense and assigned in interning order, so this is usable as
    /// a key into a caller-side side table.
    #[must_use]
    pub fn get(self) -> usize {
        self.0
    }
}

impl From<usize> for Index {
    fn from(ix: usize) -> Self {
        Self(ix)
    }
}

impl From<Index> for usize {
    fn from(idx: Index) -> Self {
        idx.0
    }
}

fn erfc(x: f64) -> f64 {
    let z = x.abs();
    let t = 1.0 / (1.0 + z / 2.0);

    let a = -0.82215223 + t * 0.17087277;
    let b = 1.48851587 + t * a;
    let c = -1.13520398 + t * b;
    let d = 0.27886807 + t * c;
    let e = -0.18628806 + t * d;
    let f = 0.09678418 + t * e;
    let g = 0.37409196 + t * f;
    let h = 1.00002368 + t * g;

    let r = t * (-z * z - 1.26551223 + t * h).exp();

    if x >= 0.0 { r } else { 2.0 - r }
}

fn erfc_inv(mut y: f64) -> f64 {
    if y >= 2.0 {
        return f64::NEG_INFINITY;
    }

    debug_assert!(y >= 0.0, "y must be nonnegative");

    if y == 0.0 {
        return f64::INFINITY;
    }

    if y >= 1.0 {
        y = 2.0 - y;
    }

    let t = (-2.0 * (y / 2.0).ln()).sqrt();

    let mut x = FRAC_1_SQRT_2 * ((2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t);

    for _ in 0..3 {
        let err = erfc(x) - y;

        x += err / (FRAC_2_SQRT_PI * (-(x.powi(2))).exp() - x * err)
    }

    if y < 1.0 { x } else { -x }
}

fn ppf(p: f64, mu: f64, sigma: f64) -> f64 {
    mu - sigma * SQRT_2 * erfc_inv(2.0 * p)
}

fn compute_margin(p_draw: f64, sd: f64) -> f64 {
    ppf(0.5 - p_draw / 2.0, 0.0, sd).abs()
}

pub(crate) fn cdf(x: f64, mu: f64, sigma: f64) -> f64 {
    let z = -(x - mu) / (sigma * SQRT_2);

    0.5 * erfc(z)
}

fn pdf(x: f64, mu: f64, sigma: f64) -> f64 {
    let normalizer = (SQRT_TAU * sigma).powi(-1);
    let functional = (-((x - mu).powi(2)) / (2.0 * sigma.powi(2))).exp();

    normalizer * functional
}

fn v_w(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
    if !tie {
        let alpha = (margin - mu) / sigma;

        let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
        let w = v * (v + (-alpha));

        (v, w)
    } else {
        let alpha = (-margin - mu) / sigma;
        let beta = (margin - mu) / sigma;

        let v = (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0))
            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
        let u = (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0))
            / (cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0));
        let w = -(u - v.powi(2));

        (v, w)
    }
}

fn trunc(mu: f64, sigma: f64, margin: f64, tie: bool) -> (f64, f64) {
    let (v, w) = v_w(mu, sigma, margin, tie);

    let mu_trunc = mu + sigma * v;
    let sigma_trunc = sigma * (1.0 - w).sqrt();

    (mu_trunc, sigma_trunc)
}

pub(crate) fn approx(n: Gaussian, margin: f64, tie: bool) -> Gaussian {
    let (mu, sigma) = trunc(n.mu(), n.sigma(), margin, tie);

    Gaussian::from_ms(mu, sigma)
}

pub(crate) fn tuple_max(v1: (f64, f64), v2: (f64, f64)) -> (f64, f64) {
    (
        if v1.0 > v2.0 { v1.0 } else { v2.0 },
        if v1.1 > v2.1 { v1.1 } else { v2.1 },
    )
}

pub(crate) fn tuple_gt(t: (f64, f64), e: f64) -> bool {
    t.0 > e || t.1 > e
}

/// Whether a convergence step is finite in both components.
///
/// A NaN step means EP broke down numerically. Because every comparison
/// against NaN is false, `tuple_gt` reads NaN as "below epsilon" — so
/// convergence checks must test finiteness explicitly rather than inferring
/// success from `!tuple_gt(..)`.
pub(crate) fn step_is_finite(t: (f64, f64)) -> bool {
    t.0.is_finite() && t.1.is_finite()
}

/// Whether a step counts as converged: finite *and* within `epsilon`.
pub(crate) fn step_converged(t: (f64, f64), epsilon: f64) -> bool {
    step_is_finite(t) && !tuple_gt(t, epsilon)
}

/// Indices of the first pair of teams sharing a rank, if any.
///
/// A tie is only representable when the draw probability is positive: with
/// `p_draw == 0.0` the truncation margin collapses to zero and the two-sided
/// tie update evaluates `0/0`. Callers use this to reject such events before
/// they reach inference.
pub(crate) fn first_tied_pair(ranks: &[u32]) -> Option<(usize, usize)> {
    for (i, a) in ranks.iter().enumerate() {
        for (j, b) in ranks.iter().enumerate().skip(i + 1) {
            if a == b {
                return Some((i, j));
            }
        }
    }

    None
}

/// As `first_tied_pair`, but over the engine's internal `f64` outputs.
///
/// Ranks reach the engine already converted to descending `f64` outputs, and
/// `Game` decides a tie by exact equality of those values — so this mirrors
/// the comparison inference itself performs.
pub(crate) fn first_tied_output(outputs: &[f64]) -> Option<(usize, usize)> {
    for (i, a) in outputs.iter().enumerate() {
        for (j, b) in outputs.iter().enumerate().skip(i + 1) {
            if a == b {
                return Some((i, j));
            }
        }
    }

    None
}

pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
    let mut x: Vec<(usize, T)> = xs.iter().enumerate().map(|(i, &t)| (i, t)).collect();

    if reverse {
        x.sort_by_key(|&(_, t)| Reverse(t));
    } else {
        x.sort_by_key(|&(_, t)| t);
    }

    x.into_iter().map(|(i, _)| i).collect()
}

/// Calculates the match quality of the given rating groups. A result is the draw probability in the association
///
/// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a
/// perfectly balanced match.
///
/// # Panics
///
/// Panics if fewer than two rating groups are supplied, or if any group is
/// empty — match quality is a property of a contest between at least two
/// non-empty sides.
#[must_use]
pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
    assert!(
        rating_groups.len() >= 2,
        "quality() requires at least 2 rating groups, got {}",
        rating_groups.len()
    );
    assert!(
        rating_groups.iter().all(|group| !group.is_empty()),
        "quality() requires every rating group to be non-empty"
    );

    let flatten_ratings = rating_groups
        .iter()
        .flat_map(|group| group.iter())
        .collect::<Vec<_>>();

    let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice();

    let length = flatten_ratings.len();

    let mut mean_matrix = Matrix::new(length, 1);

    for (i, rating) in flatten_ratings.iter().enumerate() {
        mean_matrix[(i, 0)] = rating.mu();
    }

    let mut variance_matrix = Matrix::new(length, length);

    for (i, rating) in flatten_ratings.iter().enumerate() {
        variance_matrix[(i, i)] = rating.sigma().powi(2);
    }

    let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length);

    // Row `row` contrasts group `row` (+weight) against group `row + 1`
    // (-weight). `t` is the column where the current group's players start;
    // the negative block begins immediately after it.
    let mut t = 0;

    for (row, group) in rating_groups.windows(2).enumerate() {
        let current = group[0];
        let next = group[1];

        for n in t..t + current.len() {
            rotated_a_matrix[(row, n)] = flatten_weights[n];
        }

        t += current.len();

        for n in t..t + next.len() {
            rotated_a_matrix[(row, n)] = -flatten_weights[n];
        }
    }

    let a_matrix = rotated_a_matrix.transpose();

    let ata = beta.powi(2) * &rotated_a_matrix * &a_matrix;
    let atsa = &rotated_a_matrix * &variance_matrix * &a_matrix;

    let start = mean_matrix.transpose() * &a_matrix;
    let middle = &ata + &atsa;
    let end = &rotated_a_matrix * &mean_matrix;

    let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
    let s_arg = ata.determinant() / middle.determinant();

    e_arg.exp() * s_arg.sqrt()
}

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

    use super::*;

    #[test]
    fn test_sort_time() {
        assert_eq!(sort_time(&[0i64, 1, 2, 0], true), vec![2, 1, 0, 3]);
    }

    #[test]
    fn test_quality() {
        let a = Gaussian::from_ms(25.0, 3.0);
        let b = Gaussian::from_ms(25.0, 3.0);

        let q = quality(&[&[a], &[b]], 25.0 / 3.0 / 2.0);

        assert_ulps_eq!(q, 0.8115343414514944, epsilon = 1e-6)
    }
}