#!/usr/bin/env python3
"""Exact verifier for the sharp degree-thirteen Karlin-Wronskian theorem."""

from __future__ import annotations

from fractions import Fraction
from itertools import combinations, permutations
from math import comb, factorial

import numpy as np
import sympy as sp


EXPECTED_MONOMIAL_COEFFICIENTS = {
    (4, 4, 4): 108,
    (3, 3, 3, 3): -288,
    (4, 3, 3, 2): 96,
    (4, 4, 2, 2): -24,
    (3, 3, 2, 2, 2): 96,
    (4, 2, 2, 2, 2): 72,
    (2, 2, 2, 2, 2, 2): 720,
}
EXPECTED_JOHNSON_TABLE = [
    [1, 27, 108, 84],
    [1, 15, 12, -28],
    [1, 5, -13, 7],
    [1, -3, 3, -1],
]
EXPECTED_L_EIGENVALUES = [
    Fraction(0),
    Fraction(0),
    Fraction(0),
    Fraction(90, 7),
]
EXPECTED_WITNESS_DETERMINANT = -527889869252540784
EXPECTED_BOUNDARY_MINOR = -16384


if not __debug__:
    raise SystemExit("verify.py requires assertions")


def determinant_polynomial() -> sp.Expr:
    q = sp.symbols("q1:7")
    derivatives = [sp.Integer(1)]
    for degree in range(6):
        value = 0
        for index in range(degree + 1):
            cumulant = (-1) ** index * sp.factorial(index) * q[index]
            value += (
                sp.binomial(degree, index)
                * derivatives[degree - index]
                * cumulant
            )
        derivatives.append(sp.expand(value))

    matrix = sp.Matrix(
        [
            [derivatives[3 + column - row] for column in range(4)]
            for row in range(4)
        ]
    )
    determinant = sp.Poly(sp.expand(matrix.det()), *q)
    assert all(coefficient % 12 == 0 for coefficient in determinant.coeffs())
    return sp.expand(determinant.as_expr() / 12)


def monomial_coefficients(polynomial: sp.Expr) -> dict[tuple[int, ...], int]:
    q = sp.symbols("q1:7")
    roots = sp.symbols("alpha0:6")
    substituted = polynomial.subs(
        {
            q[degree - 1]: sum(root**degree for root in roots)
            for degree in range(2, 7)
        }
    )
    expanded = sp.Poly(sp.expand(substituted), *roots)
    coefficients: dict[tuple[int, ...], int] = {}
    for exponents, coefficient in expanded.terms():
        partition = tuple(
            sorted((exponent for exponent in exponents if exponent), reverse=True)
        )
        integer = int(coefficient)
        previous = coefficients.setdefault(partition, integer)
        assert previous == integer
    return coefficients


def verify_am_gm_orbit() -> None:
    orbit = sorted(set(permutations((4, 3, 3, 2))))
    assert len(orbit) == 12
    exponent_totals = [
        sum(exponents[index] for exponents in orbit)
        for index in range(4)
    ]
    assert exponent_totals == [36, 36, 36, 36]


def johnson_eigenvalue(n: int, distance: int, level: int) -> int:
    value = 0
    for index in range(distance + 1):
        value += (
            (-1) ** (distance - index)
            * comb(3 - index, distance - index)
            * comb(3 - level, index)
            * comb(n - 3 + index - level, index)
        )
    return value


def verify_johnson_certificate() -> None:
    table = [
        [johnson_eigenvalue(12, distance, level) for distance in range(4)]
        for level in range(4)
    ]
    assert table == EXPECTED_JOHNSON_TABLE

    coefficients = [
        Fraction(9),
        Fraction(-1),
        Fraction(1, 4),
        Fraction(-3, 28),
    ]
    eigenvalues = [
        sum(
            coefficient * eigenvalue
            for coefficient, eigenvalue in zip(coefficients, row, strict=True)
        )
        for row in table
    ]
    assert eigenvalues == EXPECTED_L_EIGENVALUES

    kernel_values = [
        Fraction(9),
        Fraction(-1),
        Fraction(1),
        Fraction(3),
    ]
    remainder = [
        kernel - certificate
        for kernel, certificate in zip(kernel_values, coefficients, strict=True)
    ]
    assert remainder == [
        Fraction(0),
        Fraction(0),
        Fraction(3, 4),
        Fraction(87, 28),
    ]


def verify_incidence_projection() -> None:
    triples = list(combinations(range(12), 3))
    pairs = list(combinations(range(12), 2))
    incidence = np.zeros((len(pairs), len(triples)), dtype=np.int64)
    for pair_index, pair in enumerate(pairs):
        pair_set = set(pair)
        for triple_index, triple in enumerate(triples):
            incidence[pair_index, triple_index] = int(
                pair_set.issubset(triple)
            )

    moment = incidence.T @ incidence
    identity = np.eye(len(triples), dtype=np.int64)
    minimal_polynomial = (
        moment
        @ (moment - 8 * identity)
        @ (moment - 18 * identity)
        @ (moment - 30 * identity)
    )
    assert not np.any(minimal_polynomial)

    scaled_certificate = np.empty_like(moment)
    for left_index, left in enumerate(triples):
        left_set = set(left)
        for right_index, right in enumerate(triples):
            scaled_certificate[left_index, right_index] = {
                3: 252,
                2: -28,
                1: 7,
                0: -3,
            }[len(left_set.intersection(right))]

    # The certificate annihilates the pair-incidence space.
    assert not np.any(scaled_certificate @ incidence.T)

    projection_polynomial = (
        (8 * identity - moment)
        @ (18 * identity - moment)
        @ (30 * identity - moment)
    )
    assert np.array_equal(projection_polynomial, 12 * scaled_certificate)


def verify_scope_boundary() -> None:
    """Show why the theorem must not be called full-PF4 minimality."""
    matrix = sp.Matrix(
        [
            [16, 128, 0, 0],
            [1, 16, 128, 0],
            [0, 1, 16, 128],
            [0, 0, 1, 16],
        ]
    )
    assert int(matrix.det()) == EXPECTED_BOUNDARY_MINOR


def centered_witness_determinant() -> int:
    coefficients = [
        sum(
            comb(11, degree - offset) * [1, 324, 26244][offset]
            for offset in range(3)
            if 0 <= degree - offset <= 11
        )
        for degree in range(14)
    ]
    derivatives = [
        factorial(degree) * coefficients[degree]
        for degree in range(7)
    ]
    matrix = sp.Matrix(
        [
            [derivatives[3 + column - row] for column in range(4)]
            for row in range(4)
        ]
    )
    return int(matrix.det())


def main() -> None:
    polynomial = determinant_polynomial()
    assert monomial_coefficients(polynomial) == EXPECTED_MONOMIAL_COEFFICIENTS
    verify_am_gm_orbit()
    verify_johnson_certificate()
    verify_incidence_projection()
    verify_scope_boundary()
    assert centered_witness_determinant() == EXPECTED_WITNESS_DETERMINANT
    print("monomial expansion: PASS")
    print("AM-GM orbit count: PASS")
    print("Johnson L eigenvalues:", EXPECTED_L_EIGENVALUES)
    print("220x220 incidence projection identity: PASS")
    print("entrywise nonnegative remainder: PASS")
    print("full-PF4 scope boundary: PASS")
    print("degree-thirteen negative witness: PASS")


if __name__ == "__main__":
    main()
