"""Exact permutation checks and figures for the strong convergence post.

Run from any directory with Python 3, NumPy, and Matplotlib installed.
All finite expectation comparisons use integers and fractions. The figures
illustrate explicit formulas, not simulated spectral or rare-event data.
"""

from fractions import Fraction
from itertools import permutations
from math import factorial, floor, sqrt
from pathlib import Path


def compose(left, right):
    return tuple(left[right[i]] for i in range(len(left)))


def inverse(perm):
    result = [0] * len(perm)
    for i, target in enumerate(perm):
        result[target] = i
    return tuple(result)


def fixed_points(perm):
    return sum(i == target for i, target in enumerate(perm))


def exact_checks():
    checks = 0
    for n in range(2, 8):
        sums = [0] * 8
        for perm in permutations(range(n)):
            power = tuple(range(n))
            for k in range(1, 9):
                power = compose(perm, power)
                sums[k - 1] += fixed_points(power)
        for k, total in enumerate(sums, start=1):
            expected = sum(k % length == 0 for length in range(1, n + 1))
            assert Fraction(total, factorial(n)) == expected
            checks += 1
    print(f"Single-permutation divisor identities: {checks} exact checks passed.")
    print("N | pairs | expected full fixed points | expected restricted trace / N")
    for n in range(2, 6):
        perms = list(permutations(range(n)))
        inverses = {perm: inverse(perm) for perm in perms}
        total = 0
        for left in perms:
            for right in perms:
                commutator = compose(
                    compose(compose(left, right), inverses[left]), inverses[right]
                )
                total += fixed_points(commutator)
        count = factorial(n) ** 2
        expectation = Fraction(total, count)
        restricted = (expectation - 1) / n
        assert expectation == Fraction(n, n - 1)
        assert restricted == Fraction(1, n * (n - 1))
        print(f"{n} | {count} | {expectation} | {restricted}")


def make_figures():
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    import numpy as np

    out = Path(__file__).resolve().parents[2] / "images" / "2026 Fall"
    out.mkdir(parents=True, exist_ok=True)
    plt.rcParams.update(
        {
            "font.family": "DejaVu Sans",
            "font.size": 11,
            "axes.spines.top": False,
            "axes.spines.right": False,
            "axes.labelcolor": "#313244",
            "text.color": "#313244",
            "xtick.color": "#313244",
            "ytick.color": "#313244",
            "svg.fonttype": "path",
            "svg.hashsalt": "strong-convergence-2405-16026",
            "figure.facecolor": "#eff1f5",
            "axes.facecolor": "#eff1f5",
        }
    )
    metadata = {"Date": None, "Creator": "strong_convergence_checks.py"}

    # Smooth step exp(-1/t)/(exp(-1/t)+exp(-1/(1-t))), with flat endpoints.
    def smooth_step(t):
        result = np.zeros_like(t)
        result[t >= 1] = 1
        inside = (t > 0) & (t < 1)
        left = np.exp(-1 / t[inside])
        right = np.exp(-1 / (1 - t[inside]))
        result[inside] = left / (left + right)
        return result

    radius = 2 * sqrt(3)
    cap = 4.0
    epsilon = 0.36
    fig, ax = plt.subplots(figsize=(8.2, 3.7), layout="constrained")
    x = np.linspace(-cap, cap, 3000)
    detector = smooth_step((np.abs(x) - radius - epsilon / 2) / (epsilon / 2))
    ax.axvspan(-radius, radius, color="#7287fd", alpha=0.18)
    ax.plot(x, detector, color="#8839ef", linewidth=2.6)
    ax.axvspan(radius + epsilon / 2, radius + epsilon, color="#9ca0b0", alpha=0.25)
    ax.axvspan(-radius - epsilon, -radius - epsilon / 2, color="#9ca0b0", alpha=0.25)
    ax.axvline(radius, color="#7287fd", linestyle=":", linewidth=1.3)
    ax.axvline(-radius, color="#7287fd", linestyle=":", linewidth=1.3)
    ax.text(0, 0.58, "Limiting spectral interval\nThe detector is identically zero here", ha="center")
    ax.annotate(
        "Outlier counted with weight one",
        xy=(radius + epsilon + 0.05, 1),
        xytext=(0.0, 1.25),
        ha="center",
        arrowprops={"arrowstyle": "->", "color": "#313244"},
    )
    ax.set(xlim=(-4.08, 4.08), ylim=(-0.08, 1.48), ylabel=r"Detector $\chi(x)$", xlabel="Spectral parameter x")
    ax.set_xticks([-radius, 0, radius], [r"$-R$", "0", r"$R$"])
    ax.set_yticks([0, 1])
    ax.set_title(r"Smooth detection of outliers   $R=2\sqrt{3}$, $K=4$", pad=14)
    path = out / "strong-convergence-spectral-detector.svg"
    fig.savefig(path, metadata=metadata)
    fig.savefig("/tmp/strong-convergence-spectral-detector.png", dpi=150)
    plt.close(fig)

    d = 6
    b = 2 * d - 1
    m_star = floor((sqrt(b) + 1) / 2)

    def rho(m):
        return 2 * sqrt(b) if 2 * m - 1 <= sqrt(b) else (2 * m - 1) + b / (2 * m - 1)

    fig, ax = plt.subplots(figsize=(8.2, 4.3), layout="constrained")
    colors = ["#1e66f5", "#8839ef", "#d20f39", "#179299"]
    for m, color in zip(range(m_star, d), colors):
        left, right = rho(m), rho(m + 1)
        ax.plot([left, right], [m, m], color=color, linewidth=3)
        ax.plot([left, right], [m, m], linestyle="none", marker="o", markersize=7,
                markerfacecolor="#eff1f5", markeredgecolor=color, markeredgewidth=1.8)
        ax.text((left + right) / 2, m + 0.18, rf"$N^{{-{m}}}$", ha="center", color=color)
    ticks = [rho(m) for m in range(m_star, d + 1)]
    labels = [r"$2\sqrt{11}$", "7.2", r"$60/7$", r"$92/9$", "12"]
    ax.set_xticks(ticks, labels)
    ax.set_yticks(range(m_star, d))
    ax.set(xlim=(ticks[0] - 0.25, 12.2), ylim=(1.6, 5.65),
           xlabel="Fixed threshold x", ylabel="Asymptotic power of 1/N")
    ax.grid(axis="y", alpha=0.2)
    ax.set_title("Rare-event staircase for degree 12\nTail probabilities on open threshold intervals", pad=12)
    fig.savefig(out / "strong-convergence-staircase.svg", metadata=metadata)
    fig.savefig("/tmp/strong-convergence-staircase.png", dpi=150)
    plt.close(fig)
    for name in ("strong-convergence-spectral-detector.svg", "strong-convergence-staircase.svg"):
        svg = out / name
        svg.write_text(
            "\n".join(line.rstrip() for line in svg.read_text().splitlines()) + "\n"
        )
    print(f"Figures written to {out}")


if __name__ == "__main__":
    exact_checks()
    make_figures()