//! Minimal dense matrix used by `quality()`.
//!
//! `determinant` and `inverse` go through one LU decomposition with partial
//! pivoting — O(n³) and numerically stable. The previous implementation
//! expanded cofactors recursively (O(n!), allocating a `Vec` per minor) and
//! only implemented `inverse` for the 1×1 case, which limited `quality()` to
//! exactly two rating groups.

use std::ops;

#[derive(Clone, Debug)]
pub struct Matrix {
    data: Box<[f64]>,
    height: usize,
    width: usize,
}

/// LU decomposition with partial pivoting: `PA = LU`, stored compactly.
///
/// `lu` holds `L` below the diagonal (unit diagonal implied) and `U` on and
/// above it. `sign` is the determinant sign contributed by row swaps, or 0.0
/// when the matrix is singular.
struct Lu {
    lu: Vec<f64>,
    perm: Vec<usize>,
    n: usize,
    sign: f64,
}

impl Lu {
    fn decompose(m: &Matrix) -> Self {
        debug_assert_eq!(m.width, m.height, "LU requires a square matrix");

        let n = m.width;
        let mut lu = m.data.to_vec();
        let mut perm: Vec<usize> = (0..n).collect();
        let mut sign = 1.0;

        for col in 0..n {
            // Partial pivot: take the largest-magnitude candidate to limit
            // growth of round-off in the elimination below.
            let mut pivot_row = col;
            let mut pivot_max = lu[col * n + col].abs();

            for row in (col + 1)..n {
                let candidate = lu[row * n + col].abs();
                if candidate > pivot_max {
                    pivot_max = candidate;
                    pivot_row = row;
                }
            }

            if pivot_max == 0.0 {
                sign = 0.0;
                continue;
            }

            if pivot_row != col {
                for k in 0..n {
                    lu.swap(col * n + k, pivot_row * n + k);
                }
                perm.swap(col, pivot_row);
                sign = -sign;
            }

            let pivot = lu[col * n + col];

            for row in (col + 1)..n {
                let factor = lu[row * n + col] / pivot;
                lu[row * n + col] = factor;

                for k in (col + 1)..n {
                    lu[row * n + k] -= factor * lu[col * n + k];
                }
            }
        }

        Self { lu, perm, n, sign }
    }

    fn determinant(&self) -> f64 {
        if self.sign == 0.0 {
            return 0.0;
        }

        let mut det = self.sign;
        for i in 0..self.n {
            det *= self.lu[i * self.n + i];
        }

        det
    }

    /// `ln |det|`, accumulated term by term rather than multiplied out.
    ///
    /// The determinant of an `n x n` Gram matrix is a product of `n` diagonal
    /// entries, so it leaves `f64`'s range long before the quantities built
    /// from it do. `quality()` only ever wants a *ratio* of two determinants,
    /// and that ratio is perfectly representable while the determinants
    /// themselves are not — measured, at 250 rating groups both overflow and
    /// the ratio came back `NaN` where the true answer is `9.51e-88`.
    ///
    /// Returns `-inf` for a singular matrix, so `exp` of it is zero.
    fn ln_abs_determinant(&self) -> f64 {
        if self.sign == 0.0 {
            return f64::NEG_INFINITY;
        }

        let mut acc = 0.0;
        for i in 0..self.n {
            acc += libm::log(self.lu[i * self.n + i].abs());
        }

        acc
    }

    /// Solve `Ax = b` for a single column of the identity, giving one column
    /// of the inverse.
    fn solve_column(&self, col: usize, out: &mut [f64]) {
        let n = self.n;

        // Forward substitution through L, applying the row permutation.
        for i in 0..n {
            let mut sum = if self.perm[i] == col { 1.0 } else { 0.0 };
            for (k, &solved) in out.iter().enumerate().take(i) {
                sum -= self.lu[i * n + k] * solved;
            }
            out[i] = sum;
        }

        // Back substitution through U.
        for i in (0..n).rev() {
            let mut sum = out[i];
            for (k, &solved) in out.iter().enumerate().skip(i + 1) {
                sum -= self.lu[i * n + k] * solved;
            }
            out[i] = sum / self.lu[i * n + i];
        }
    }
}

impl Matrix {
    pub(crate) fn new(height: usize, width: usize) -> Matrix {
        Matrix {
            data: vec![0.0; height * width].into_boxed_slice(),
            height,
            width,
        }
    }

    pub(crate) fn transpose(&self) -> Matrix {
        let mut matrix = Matrix::new(self.width, self.height);

        for c in 0..self.width {
            for r in 0..self.height {
                matrix[(c, r)] = self[(r, c)];
            }
        }

        matrix
    }

    /// Determinant of a square matrix. The 0×0 determinant is 1 by convention
    /// (the empty product).
    ///
    /// # Panics
    ///
    /// Panics if the matrix is not square.
    pub(crate) fn determinant(&self) -> f64 {
        assert_eq!(
            self.width, self.height,
            "determinant requires a square matrix, got {}x{}",
            self.height, self.width
        );

        if self.width == 0 {
            return 1.0;
        }

        Lu::decompose(self).determinant()
    }

    /// `ln |det|` of a square matrix; `-inf` when singular.
    ///
    /// See [`Lu::ln_abs_determinant`] for why a ratio of determinants must be
    /// taken this way.
    pub(crate) fn ln_abs_determinant(&self) -> f64 {
        assert_eq!(
            self.width, self.height,
            "determinant requires a square matrix, got {}x{}",
            self.height, self.width
        );

        if self.width == 0 {
            return 0.0;
        }

        Lu::decompose(self).ln_abs_determinant()
    }

    /// Matrix inverse via LU decomposition.
    ///
    /// # Panics
    ///
    /// Panics if the matrix is not square or is singular.
    pub(crate) fn inverse(&self) -> Matrix {
        assert_eq!(
            self.width, self.height,
            "inverse requires a square matrix, got {}x{}",
            self.height, self.width
        );

        let n = self.width;
        let mut inverse = Matrix::new(n, n);

        if n == 0 {
            return inverse;
        }

        let lu = Lu::decompose(self);
        assert!(lu.sign != 0.0, "cannot invert a singular matrix");

        let mut column = vec![0.0; n];

        for c in 0..n {
            lu.solve_column(c, &mut column);

            for (r, &value) in column.iter().enumerate() {
                inverse[(r, c)] = value;
            }
        }

        inverse
    }
}

impl ops::Index<(usize, usize)> for Matrix {
    type Output = f64;

    fn index(&self, pos: (usize, usize)) -> &Self::Output {
        debug_assert!(
            pos.0 < self.height && pos.1 < self.width,
            "index ({}, {}) out of bounds for {}x{} matrix",
            pos.0,
            pos.1,
            self.height,
            self.width
        );

        &self.data[(self.width * pos.0) + pos.1]
    }
}

impl ops::IndexMut<(usize, usize)> for Matrix {
    fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output {
        debug_assert!(
            pos.0 < self.height && pos.1 < self.width,
            "index ({}, {}) out of bounds for {}x{} matrix",
            pos.0,
            pos.1,
            self.height,
            self.width
        );

        &mut self.data[(self.width * pos.0) + pos.1]
    }
}

fn multiply(lhs: &Matrix, rhs: &Matrix) -> Matrix {
    assert_eq!(
        lhs.width, rhs.height,
        "cannot multiply {}x{} by {}x{}",
        lhs.height, lhs.width, rhs.height, rhs.width
    );

    let mut matrix = Matrix::new(lhs.height, rhs.width);

    for r in 0..matrix.height {
        for c in 0..matrix.width {
            let mut value = 0.0;

            for x in 0..lhs.width {
                value += lhs[(r, x)] * rhs[(x, c)];
            }

            matrix[(r, c)] = value;
        }
    }

    matrix
}

impl ops::Mul<&Matrix> for f64 {
    type Output = Matrix;

    fn mul(self, rhs: &Matrix) -> Matrix {
        let mut matrix = Matrix::new(rhs.height, rhs.width);

        for r in 0..rhs.height {
            for c in 0..rhs.width {
                matrix[(r, c)] = self * rhs[(r, c)];
            }
        }

        matrix
    }
}

impl ops::Mul<&Matrix> for Matrix {
    type Output = Matrix;

    fn mul(self, rhs: &Matrix) -> Matrix {
        multiply(&self, rhs)
    }
}

impl ops::Mul<&Matrix> for &Matrix {
    type Output = Matrix;

    fn mul(self, rhs: &Matrix) -> Matrix {
        multiply(self, rhs)
    }
}

impl ops::Add<&Matrix> for &Matrix {
    type Output = Matrix;

    fn add(self, rhs: &Matrix) -> Matrix {
        assert!(
            self.height == rhs.height && self.width == rhs.width,
            "cannot add {}x{} to {}x{}",
            self.height,
            self.width,
            rhs.height,
            rhs.width
        );

        let mut matrix = Matrix::new(self.height, self.width);

        for r in 0..matrix.height {
            for c in 0..matrix.width {
                matrix[(r, c)] = self[(r, c)] + rhs[(r, c)];
            }
        }

        matrix
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn from_rows(rows: &[&[f64]]) -> Matrix {
        let mut m = Matrix::new(rows.len(), rows[0].len());
        for (r, row) in rows.iter().enumerate() {
            for (c, &v) in row.iter().enumerate() {
                m[(r, c)] = v;
            }
        }
        m
    }

    #[test]
    fn determinant_1x1() {
        assert!((from_rows(&[&[3.0]]).determinant() - 3.0).abs() < 1e-12);
    }

    #[test]
    fn determinant_2x2() {
        let m = from_rows(&[&[1.0, 2.0], &[3.0, 4.0]]);
        assert!((m.determinant() - (-2.0)).abs() < 1e-12);
    }

    #[test]
    fn determinant_3x3() {
        let m = from_rows(&[&[6.0, 1.0, 1.0], &[4.0, -2.0, 5.0], &[2.0, 8.0, 7.0]]);
        assert!((m.determinant() - (-306.0)).abs() < 1e-10);
    }

    #[test]
    fn determinant_requires_no_pivot_at_origin() {
        // A zero in the top-left forces a row swap; the sign must follow.
        let m = from_rows(&[&[0.0, 1.0], &[1.0, 0.0]]);
        assert!((m.determinant() - (-1.0)).abs() < 1e-12);
    }

    #[test]
    fn determinant_of_singular_is_zero() {
        let m = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]);
        assert!(m.determinant().abs() < 1e-12);
    }

    #[test]
    fn inverse_1x1() {
        let inv = from_rows(&[&[4.0]]).inverse();
        assert!((inv[(0, 0)] - 0.25).abs() < 1e-12);
    }

    #[test]
    fn inverse_times_original_is_identity() {
        for rows in [
            vec![vec![1.0, 2.0], vec![3.0, 4.0]],
            vec![
                vec![6.0, 1.0, 1.0],
                vec![4.0, -2.0, 5.0],
                vec![2.0, 8.0, 7.0],
            ],
            vec![
                vec![2.0, 0.0, 1.0, 3.0],
                vec![1.0, 5.0, 2.0, 0.0],
                vec![0.0, 1.0, 4.0, 1.0],
                vec![3.0, 2.0, 0.0, 6.0],
            ],
        ] {
            let refs: Vec<&[f64]> = rows.iter().map(|r| r.as_slice()).collect();
            let m = from_rows(&refs);
            let product = &m * &m.inverse();

            for r in 0..product.height {
                for c in 0..product.width {
                    let expected = if r == c { 1.0 } else { 0.0 };
                    assert!(
                        (product[(r, c)] - expected).abs() < 1e-9,
                        "({r},{c}) = {} expected {expected}",
                        product[(r, c)]
                    );
                }
            }
        }
    }

    #[test]
    #[should_panic(expected = "singular")]
    fn inverse_of_singular_panics() {
        let _ = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]).inverse();
    }

    #[test]
    fn empty_determinant_is_one() {
        assert!((Matrix::new(0, 0).determinant() - 1.0).abs() < 1e-12);
    }

    #[test]
    fn transpose_round_trips() {
        let m = from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
        let t = m.transpose();
        assert_eq!((t.height, t.width), (3, 2));
        assert_eq!(t.transpose()[(1, 2)], m[(1, 2)]);
    }
}