//! `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)]

/// Compiles every `rust` block in `README.md` as a doctest.
///
/// The README is not the crate's front page — the module docs above are — so it
/// is pulled in here rather than via a crate-level `#![doc = ...]`, purely so
/// its examples are type-checked. Without this nothing compiled them, and they
/// had drifted far enough that four blocks no longer built (#35). `cfg(doctest)`
/// means this type exists only while collecting doctests.
///
/// Blocks that are illustrative rather than runnable are fenced as `text`.
#[cfg(doctest)]
#[doc = include_str!("../README.md")]
pub struct ReadmeDoctests;

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 acquisition;
mod color_group;
mod competitor;
mod convergence;
pub mod drift;
mod error;
mod event;
mod event_builder;
pub(crate) mod factor;
mod game;
pub mod gaussian;
pub mod graph;
mod history;
mod key_table;
mod matrix;
mod observer;
mod outcome;
mod predict;
pub(crate) mod quadrature;
mod rating;
pub(crate) mod schedule;
pub mod storage;

pub use acquisition::expected_information_gain;
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 predict::Prediction;
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;

/// 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;
/// `1 / sqrt(pi)`, the leading factor of the `erfcx` continued fraction.
const FRAC_1_SQRT_PI: f64 = 0.564_189_583_547_756_3;
/// `sqrt(2 / pi)`, the numerator of the inverse Mills ratio in scaled form.
const SQRT_2_OVER_PI: f64 = 0.797_884_560_802_865_4;
/// How many window widths into the tail before a tie window is treated as a
/// half-line. Beyond this the truncated mass is concentrated within `1/alpha`
/// of the near edge, so the far edge contributes nothing measurable.
const HALF_LINE_WINDOW: f64 = 10.0;
/// Where `v - alpha` switches from subtraction to its asymptotic series.
///
/// The subtraction loses roughly `eps * alpha^2` of relative precision, and the
/// four-term series is good to ~1e-10 by here, so the two are at their closest
/// agreement around this point. Below it the subtraction is exact; above it the
/// series is.
const ASYMPTOTIC_MILLS_ALPHA: f64 = 100.0;

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

    // The leading coefficient is NEGATIVE. `rational - t` is negative here, so
    // a positive coefficient mirrors the starting point to `-x0` — the
    // reflection of the root. Newton then has to cross the origin to get back,
    // which a fixed iteration count does not manage: measured against the true
    // value, `erfc_inv(0.1)` returned 1.044 instead of 1.16309, and the error
    // grew as y shrank until `compute_margin` stopped being monotone in
    // `p_draw` altogether.
    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)
}

/// `P(X > x)` for `X ~ N(mu, sigma^2)`.
///
/// The survival function, computed directly rather than as `1 - cdf(..)`.
///
/// The two are algebraically identical and numerically are not. `cdf` returns
/// a value approaching 1 for an upper tail, so subtracting it from 1 cancels
/// away every significant digit the tail had: measured against this function,
/// `1 - cdf` carries 7% error by four sigma past the mean and returns exactly
/// zero beyond about 8.3 sigma — where the true value is still 1e-19 and
/// perfectly representable. `erfc` itself holds ~1e-7 *relative* accuracy down
/// to 1e-296, so the precision is there to keep; only the subtraction threw it
/// away.
///
/// This matters most where evidence is smallest, which is exactly where an
/// upset makes it interesting: `ln` of a clamped zero is -708 regardless of
/// whether the truth was -43 or -600.
pub(crate) fn sf(x: f64, mu: f64, sigma: f64) -> f64 {
    0.5 * erfc((x - mu) / (sigma * SQRT_2))
}

/// `e^(x^2) * erfc(x)`, the scaled complementary error function, for `x >= 0`.
///
/// Exists so the exponential factor common to a Gaussian density and its tail
/// integral can be cancelled *analytically* instead of being computed twice
/// and divided. Both underflow to zero past about 26 sigma, and their ratio is
/// then `0/0` — finite in the limit, `NaN` in floating point.
fn erfcx(x: f64) -> f64 {
    if x < 2.0 {
        // Below the crossover neither factor is extreme: erfc is O(1) and
        // exp(x^2) is at most e^4, so the direct product is exact enough and
        // cheaper than the continued fraction.
        (x * x).exp() * erfc(x)
    } else {
        // erfcx(x) = 1/sqrt(pi) * 1/(x + (1/2)/(x + 1/(x + (3/2)/(x + ...)))),
        // evaluated by backward recurrence. Converges quickly for x >= 2 and,
        // unlike the product form, never touches an exponential.
        let mut f = 0.0;
        for n in (1..=60u32).rev() {
            f = (f64::from(n) * 0.5) / (x + f);
        }
        FRAC_1_SQRT_PI / (x + f)
    }
}

/// `ln` of the normal density at `x`.
///
/// The density itself underflows to zero past about 38 sigma, and `ln` of a
/// clamped zero is -708 whatever the truth was. The log form is a polynomial:
/// it stays exact at any separation, and the values it produces (-5001 nats at
/// 100 sigma, -500001 at 1000) are perfectly representable.
pub(crate) fn ln_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
    let z = (x - mu) / sigma;
    -(SQRT_TAU * sigma).ln() - 0.5 * z * z
}

/// `ln P(X > x)` for `X ~ N(mu, sigma^2)`.
///
/// In the upper tail the `exp(-z^2 / 2)` common to the tail integral is
/// factored out analytically via `erfcx`, so this never underflows — where
/// `sf(..).ln()` bottoms out at -708 once `erfc` itself reaches zero.
pub(crate) fn ln_sf(x: f64, mu: f64, sigma: f64) -> f64 {
    let z = (x - mu) / sigma;

    if z > 0.0 {
        // ln(0.5 * erfc(z/sqrt2)) with erfc(y) = exp(-y^2) * erfcx(y).
        -std::f64::consts::LN_2 - 0.5 * z * z + erfcx(z / SQRT_2).ln()
    } else {
        // The mass here is at least a half; nothing to lose.
        sf(x, mu, sigma).ln()
    }
}

/// `ln P(lo < X < hi)` for `X ~ N(mu, sigma^2)`.
///
/// When the interval sits in a tail both endpoint probabilities underflow
/// together, so their difference is taken in scaled form with the shared
/// exponential factored out. When it straddles the mean nothing is small and
/// the direct difference is exact.
pub(crate) fn ln_interval(lo: f64, hi: f64, mu: f64, sigma: f64) -> f64 {
    let z_lo = (lo - mu) / sigma;
    let z_hi = (hi - mu) / sigma;

    if z_hi <= z_lo {
        return f64::NEG_INFINITY;
    }

    // Fold a lower-tail interval onto the upper tail; the normal is symmetric.
    let (near, far) = if z_lo >= 0.0 {
        (z_lo, z_hi)
    } else if z_hi <= 0.0 {
        (-z_hi, -z_lo)
    } else {
        // Straddles the mean: the interval holds a non-negligible share of the
        // mass, so neither endpoint is near enough to 1 to cancel.
        return (cdf(hi, mu, sigma) - cdf(lo, mu, sigma))
            .max(f64::MIN_POSITIVE)
            .ln();
    };

    let (a, b) = (near / SQRT_2, far / SQRT_2);
    // b > a >= 0, so this ratio of exponentials is at most 1 and cannot overflow.
    let scale = (a * a - b * b).exp();
    let bracket = erfcx(a) - scale * erfcx(b);

    if bracket <= 0.0 {
        return f64::NEG_INFINITY;
    }

    -std::f64::consts::LN_2 - a * a + bracket.ln()
}

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
}

/// Truncated-Gaussian correction terms `(v, w)`.
///
/// `v` shifts the mean and `w` shrinks the variance. Both are ratios whose
/// numerator and denominator underflow together in the tails, so both are
/// computed in scaled form there: the shared `exp(-alpha^2 / 2)` is cancelled
/// analytically rather than evaluated and divided out. Without that, a
/// truncation point beyond about 39 sigma produced `0 / 0` and put `NaN`
/// straight into the posterior.
/// Truncation terms for a boundary `alpha` standard deviations into the upper
/// tail, from the asymptotic expansion of the inverse Mills ratio.
///
/// `v` tends to `alpha` out here, so the gap between them cannot be obtained by
/// subtracting one from the other — the series computes the gap directly, and
/// `w = v * gap` then never forms the difference of two large near-equal
/// numbers. A far-tail *window* behaves like a half-line once it is more than a
/// few multiples of its own width from the mean, so the tie branch shares this.
fn half_line_truncation(alpha: f64) -> (f64, f64) {
    let inv = alpha.recip();
    let inv_sq = inv * inv;
    let gap = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
    let v = alpha + gap;

    (v, v * gap)
}

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

        // v is the inverse Mills ratio, phi(alpha) / Phi(-alpha), and w needs
        // the gap `v - alpha` as well as v itself. Far into the tail v tends to
        // alpha, so that gap is a subtraction of two nearly equal numbers and
        // loses every digit it has: at alpha = 1e6 it drove w above 1 and made
        // `sqrt(1 - w)` NaN. Past the crossover the gap comes from its
        // asymptotic series instead, which has no subtraction in it.
        if alpha >= ASYMPTOTIC_MILLS_ALPHA {
            return half_line_truncation(alpha);
        }

        let (v, gap) = if alpha > 0.0 {
            // Both terms carry exp(-alpha^2 / 2); in scaled form it cancels
            // and the result stays exact however far into the tail alpha sits.
            let v = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2);
            (v, v - alpha)
        } else {
            // Phi(-alpha) >= 1/2 here, so the direct ratio loses nothing.
            let v = pdf(-alpha, 0.0, 1.0) / cdf(-alpha, 0.0, 1.0);
            (v, v - alpha)
        };

        (v, v * gap)
    } else {
        // v is odd in mu and w is even, so fold to mu <= 0. Both truncation
        // points then sit in the upper tail, where the scaled form applies.
        let flipped = mu > 0.0;
        let mu = if flipped { -mu } else { mu };

        let alpha = (-margin - mu) / sigma;
        let beta = (margin - mu) / sigma;

        // `w` comes out of `v * v - u`, and both terms grow as alpha^2 while
        // their difference stays O(1) — at alpha = 1e9 that subtraction had no
        // digits left and returned w = -128, making `sqrt(1 - w)` nonsense.
        // Once the window sits many of its own widths into the tail it is
        // indistinguishable from a half-line, so the asymptotic covers it with
        // no subtraction at all.
        if alpha >= ASYMPTOTIC_MILLS_ALPHA && alpha * (beta - alpha) >= HALF_LINE_WINDOW {
            let (v, w) = half_line_truncation(alpha);
            return (if flipped { -v } else { v }, w);
        }

        let (v, u) = if alpha > 0.0 {
            // beta > alpha > 0, so this ratio of exponentials is at most 1 and
            // cannot overflow.
            let scale = (0.5 * (alpha * alpha - beta * beta)).exp();
            let denominator = 0.5 * (erfcx(alpha / SQRT_2) - scale * erfcx(beta / SQRT_2));

            (
                (1.0 - scale) / SQRT_TAU / denominator,
                (alpha - beta * scale) / SQRT_TAU / denominator,
            )
        } else {
            // The interval straddles the mean, so nothing here is small.
            let denominator = cdf(beta, 0.0, 1.0) - cdf(alpha, 0.0, 1.0);

            (
                (pdf(alpha, 0.0, 1.0) - pdf(beta, 0.0, 1.0)) / denominator,
                (alpha * pdf(alpha, 0.0, 1.0) - beta * pdf(beta, 0.0, 1.0)) / denominator,
            )
        };

        let w = -(u - v.powi(2));

        (if flipped { -v } else { 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]);
    }

    /// Upper-tail values of the standard normal, from published tables. The
    /// point is not the digits — `erfc` only carries ~1e-7 relative — but that
    /// a number comes back at all: `1 - cdf` returned exactly zero for every
    /// one of these.
    #[test]
    fn survival_function_survives_the_far_tail() {
        for (z, expected) in [
            (9.0f64, 1.128_588e-19),
            (12.0, 1.776_482e-33),
            (20.0, 2.753_624e-89),
            (37.0, 5.725_571e-300),
        ] {
            let got = sf(z, 0.0, 1.0);
            assert!(got > 0.0, "sf({z}) collapsed to zero");
            assert!(
                (got - expected).abs() / expected < 1e-6,
                "sf({z}) = {got}, expected ~{expected}"
            );
            assert_eq!(
                1.0 - cdf(z, 0.0, 1.0),
                0.0,
                "the naive form should still be zero here"
            );
        }
    }

    /// Where no cancellation happens the two forms must agree exactly enough
    /// that nothing else in the crate shifts.
    #[test]
    fn survival_function_matches_the_naive_form_where_that_form_works() {
        for z in [-4.0f64, -1.0, 0.0, 0.5, 1.0, 2.0, 3.0, 4.0] {
            let naive = 1.0 - cdf(z, 0.0, 1.0);
            let direct = sf(z, 0.0, 1.0);
            // Bounded by `erfc`'s own ~1e-7 relative error, not by the
            // subtraction: the two forms evaluate `erfc` at different points
            // and the approximation is not exactly antisymmetric.
            assert!(
                (naive - direct).abs() < 1e-6,
                "z={z}: naive {naive} vs direct {direct}"
            );
        }
    }

    #[test]
    fn survival_and_cdf_partition_the_mass() {
        for z in [-3.0f64, -0.5, 0.0, 1.0, 2.5] {
            let total = sf(z, 1.0, 2.0) + cdf(z, 1.0, 2.0);
            // `erfc(z) + erfc(-z) == 2` only to the accuracy of the
            // approximation, which is ~1e-7 relative.
            assert!((total - 1.0).abs() < 1e-6, "z={z}: {total}");
        }
    }

    /// `erfcx` switches formulation at x = 2; the two sides must meet.
    #[test]
    fn erfcx_is_continuous_across_its_crossover() {
        for x in [1.90f64, 1.99, 1.999, 2.0, 2.001, 2.01, 2.10] {
            let direct = (x * x).exp() * erfc(x);
            let scaled = erfcx(x);
            assert!(
                (direct - scaled).abs() / scaled < 1e-6,
                "x={x}: direct {direct} vs erfcx {scaled}"
            );
        }
    }

    /// The whole reason `erfcx` exists: it stays finite and O(1/x) exactly
    /// where `exp(x^2)` overflows and `erfc(x)` underflows.
    #[test]
    fn erfcx_stays_finite_where_its_factors_do_not() {
        for x in [27.0f64, 50.0, 1.0e3, 1.0e8] {
            let scaled = erfcx(x);
            assert!(scaled.is_finite() && scaled > 0.0, "erfcx({x}) = {scaled}");
            // Asymptotically erfcx(x) -> 1 / (x * sqrt(pi)).
            let asymptote = 1.0 / (x * std::f64::consts::PI.sqrt());
            assert!(
                (scaled - asymptote).abs() / asymptote < 1e-2,
                "erfcx({x}) = {scaled} strays from its asymptote {asymptote}"
            );
            assert!(
                (x * x).exp().is_infinite(),
                "x={x} should overflow the direct form"
            );
        }
    }

    /// Truncation must never produce a non-finite posterior. Before the scaled
    /// formulation these returned NaN from `0 / 0` past about 39 sigma.
    #[test]
    fn truncation_stays_finite_arbitrarily_far_into_the_tail() {
        for alpha in [0.0f64, 8.0, 38.0, 40.0, 100.0, 1.0e3, 1.0e6, 1.0e9, 1.0e15] {
            for tie in [false, true] {
                let (v, w) = v_w(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
                assert!(v.is_finite(), "alpha={alpha} tie={tie}: v = {v}");
                assert!(w.is_finite(), "alpha={alpha} tie={tie}: w = {w}");
                // sigma_trunc = sigma * sqrt(1 - w) must stay real.
                assert!(
                    (0.0..=1.0).contains(&w),
                    "alpha={alpha} tie={tie}: w = {w} leaves sqrt(1 - w) imaginary"
                );

                let (mu_t, sigma_t) = trunc(-alpha, 1.0, if tie { 1.0 } else { 0.0 }, tie);
                assert!(
                    mu_t.is_finite() && sigma_t.is_finite(),
                    "alpha={alpha} tie={tie}: trunc = ({mu_t}, {sigma_t})"
                );
            }
        }
    }

    /// The Mills gap switches from subtraction to series at alpha = 100. Both
    /// are supposed to be right there; if they disagree, the crossover is in
    /// the wrong place.
    #[test]
    fn the_mills_gap_series_meets_the_scaled_form() {
        for alpha in [50.0f64, 99.0, 100.0, 101.0, 200.0] {
            let scaled = SQRT_2_OVER_PI / erfcx(alpha / SQRT_2) - alpha;
            let inv = alpha.recip();
            let inv_sq = inv * inv;
            let series = inv * (1.0 - inv_sq * (2.0 - inv_sq * (10.0 - 74.0 * inv_sq)));
            assert!(
                (scaled - series).abs() / series < 1e-9,
                "alpha={alpha}: scaled {scaled} vs series {series}"
            );
        }
    }

    /// Folding the tie branch to `mu <= 0` is only valid if v is odd in mu and
    /// w is even. Assert the symmetry the implementation relies on.
    #[test]
    fn tie_truncation_is_odd_in_v_and_even_in_w() {
        for mu in [0.5f64, 3.0, 20.0, 40.0, 100.0, 1.0e3] {
            let (v_pos, w_pos) = v_w(mu, 1.0, 1.0, true);
            let (v_neg, w_neg) = v_w(-mu, 1.0, 1.0, true);
            assert!(
                (v_pos + v_neg).abs() < 1e-9,
                "mu={mu}: v should be odd, got {v_pos} and {v_neg}"
            );
            assert!(
                (w_pos - w_neg).abs() < 1e-9,
                "mu={mu}: w should be even, got {w_pos} and {w_neg}"
            );
        }
    }

    /// `erfc_inv`'s initial guess had the wrong sign, putting Newton on the
    /// mirror image of the root. Three fixed iterations could not cross back,
    /// so the error grew as the argument shrank: at `p_draw = 0.99` the margin
    /// came out 0.503 where the answer is 2.576.
    #[test]
    fn erfc_inv_matches_known_quantiles() {
        // sqrt(2) * erfc_inv(1 - p) is the standard normal quantile
        // Phi^-1((1 + p) / 2).
        for (p, exact) in [
            (0.5f64, 0.674_489_750_196_081_7f64),
            (0.9, 1.644_853_626_951_472_7),
            (0.95, 1.959_963_984_540_054_2),
            (0.99, 2.575_829_303_548_9),
            (0.999, 3.290_526_731_491_896_4),
        ] {
            let got = SQRT_2 * erfc_inv(1.0 - p);
            assert!(
                (got - exact).abs() / exact < 1e-6,
                "p={p}: got {got}, exact {exact}"
            );
        }
    }

    /// The draw margin must grow with the draw probability. It did not: it ran
    /// 0.674 -> 1.476 -> 0.503 -> 0.982 as `p_draw` went 0.5 -> 0.9 -> 0.99 ->
    /// 0.999, which is not a rounding error but a broken function.
    #[test]
    fn compute_margin_is_monotone_in_the_draw_probability() {
        let mut previous = 0.0;
        for p_draw in [
            0.001f64, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999,
        ] {
            let margin = compute_margin(p_draw, 1.0);
            assert!(
                margin > previous,
                "p_draw={p_draw}: margin {margin} did not exceed {previous}"
            );
            previous = margin;
        }
    }

    /// Round-tripping the margin back through the model's own CDF must recover
    /// the draw probability it was built from.
    #[test]
    fn compute_margin_round_trips_through_the_cdf() {
        for p_draw in [0.001f64, 0.1, 0.5, 0.9, 0.99, 0.999] {
            for sd in [0.5f64, 1.0, 5.892_557] {
                let margin = compute_margin(p_draw, sd);
                // P(|X| < margin) for X ~ N(0, sd^2).
                let recovered = 1.0 - 2.0 * cdf(-margin, 0.0, sd);
                assert!(
                    (recovered - p_draw).abs() < 1e-6,
                    "p_draw={p_draw} sd={sd}: recovered {recovered}"
                );
            }
        }
    }

    /// `ln_pdf`, `ln_sf` and `ln_interval` exist so evidence stays exact where
    /// the linear forms underflow. Past ~38 sigma the linear value is zero and
    /// its log is whatever floor it was clamped to.
    #[test]
    fn log_space_helpers_stay_exact_where_the_linear_forms_underflow() {
        for z in [40.0f64, 60.0, 100.0, 1000.0] {
            assert_eq!(pdf(z, 0.0, 1.0), 0.0, "pdf should underflow at {z}");
            assert_eq!(sf(z, 0.0, 1.0), 0.0, "sf should underflow at {z}");

            let lp = ln_pdf(z, 0.0, 1.0);
            let expected_lp = -(SQRT_TAU).ln() - 0.5 * z * z;
            assert!(
                (lp - expected_lp).abs() < 1e-9,
                "ln_pdf({z}) = {lp}, expected {expected_lp}"
            );

            let ls = ln_sf(z, 0.0, 1.0);
            // ln Phi(-z) ~ -z^2/2 - ln(z) - ln(sqrt(2 pi)) for large z.
            let approx = -0.5 * z * z - z.ln() - SQRT_TAU.ln();
            assert!(
                (ls - approx).abs() / approx.abs() < 1e-3,
                "ln_sf({z}) = {ls}, asymptote {approx}"
            );
            assert!(
                ls < f64::MIN_POSITIVE.ln(),
                "ln_sf({z}) still on the clamp floor"
            );
        }
    }

    /// Where nothing underflows, the log helpers must agree with the direct
    /// forms exactly enough that nothing else in the crate shifts.
    #[test]
    fn log_space_helpers_agree_with_the_linear_forms_in_range() {
        for z in [-3.0f64, -1.0, 0.0, 1.0, 2.0, 5.0, 10.0, 20.0] {
            let lp = ln_pdf(z, 0.5, 2.0);
            let direct_pdf = pdf(z, 0.5, 2.0);
            assert!(
                (lp.exp() - direct_pdf).abs() <= 1e-12 * direct_pdf,
                "ln_pdf at {z}: {} vs {direct_pdf}",
                lp.exp()
            );

            // Bounded by `erfc`'s ~1e-7, not tighter: above x = 2 `erfcx` uses
            // a continued fraction accurate to ~1e-15, so the log path is the
            // *more* accurate of the two and they part company at `erfc`'s
            // error rather than at round-off.
            let ls = ln_sf(z, 0.5, 2.0);
            let direct = sf(z, 0.5, 2.0);
            assert!(
                (ls.exp() - direct).abs() <= 1e-6 * direct.max(1e-300),
                "ln_sf at {z}: {} vs {direct}",
                ls.exp()
            );
        }
    }

    #[test]
    fn ln_interval_matches_the_direct_difference_when_nothing_is_small() {
        for mu in [-2.0f64, 0.0, 0.5, 2.0] {
            let direct = cdf(1.0, mu, 1.0) - cdf(-1.0, mu, 1.0);
            let logged = ln_interval(-1.0, 1.0, mu, 1.0).exp();
            // See `log_space_helpers_agree_with_the_linear_forms_in_range`:
            // the gap here is `erfc`'s own error, and the log path is the more
            // accurate side of it.
            assert!(
                (logged - direct).abs() <= 1e-6 * direct,
                "mu={mu}: {logged} vs {direct}"
            );
        }
    }

    /// A window far out in the tail: both endpoints underflow together, so the
    /// difference has to be taken in scaled form.
    #[test]
    fn ln_interval_survives_a_window_deep_in_the_tail() {
        for mu in [-50.0f64, -100.0, -1000.0] {
            let logged = ln_interval(-1.0, 1.0, mu, 1.0);
            assert!(logged.is_finite(), "mu={mu}: {logged}");
            assert!(
                logged < f64::MIN_POSITIVE.ln(),
                "mu={mu}: {logged} is stuck on the clamp floor"
            );
            // Dominated by the near edge: ln P ~ ln Phi(-(|mu| - 1)).
            let near = ln_sf(-1.0, mu, 1.0);
            assert!(
                (logged - near).abs() < 5.0,
                "mu={mu}: {logged} strays from the near-edge tail {near}"
            );
        }
    }

    #[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)
    }
}