"""Illustrate the baseline series in the Lorentz counterexamples.

Run with Python, NumPy, Matplotlib, and mpmath, for example:
  uv run --no-project --with numpy --with matplotlib --with mpmath \
    python 'content/codes/2026 Fall/lorentz_baseline_cusps.py'

This is an illustration, not an exact certificate. Sixty terms give the
analytic truncation bound rho**60/(1-rho). The overview uses double-precision
angular iteration; the contact deficits use 70-digit mpmath arithmetic.
Neither plot is used to accept an identity, sign, or supremum assertion.
"""

from pathlib import Path

import matplotlib

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


TERMS = 60
CONTENT = Path(__file__).resolve().parents[2]
OUTPUT = CONTENT / "images" / "2026 Fall" / "lorentz-baseline-cusps.svg"


def baseline(theta, rho):
    phase = np.asarray(theta).copy()
    value = np.zeros_like(phase)
    coefficient = 1.0
    for _ in range(TERMS):
        value += coefficient * np.cos(phase)
        phase = np.remainder(3 * phase, 2 * np.pi)
        coefficient *= -rho
    return value


def normalized_deficit(displacement, numerator, denominator):
    """Evaluate the deficit from the exact maximum with guard precision."""
    rho = mp.mpf(numerator) / denominator
    theta = mp.pi / 4 + mp.mpf(str(displacement))
    phase = theta
    value = mp.mpf(0)
    coefficient = mp.mpf(1)
    for _ in range(TERMS):
        value += coefficient * mp.cos(phase)
        phase = mp.fmod(3 * phase, 2 * mp.pi)
        coefficient *= -rho
    maximum = (1 / mp.sqrt(2)) / (1 - rho)
    return float((maximum - value) / maximum)


def main():
    mp.mp.dps = 70
    plt.rcParams.update(
        {
            "font.size": 15,
            "axes.titlesize": 17,
            "axes.labelsize": 16,
            "legend.fontsize": 14,
            "font.family": "DejaVu Sans",
            "figure.facecolor": "#eff1f5",
            "axes.facecolor": "#eff1f5",
            "text.color": "#4c4f69",
            "axes.labelcolor": "#4c4f69",
            "xtick.color": "#4c4f69",
            "ytick.color": "#4c4f69",
            "axes.edgecolor": "#8c8fa1",
            "svg.fonttype": "path",
            "svg.hashsalt": "lorentz-baseline-cusps-v1",
        }
    )
    fig, axes = plt.subplots(2, 1, figsize=(6.4, 8.6), constrained_layout=True)
    theta = np.linspace(np.pi, 0, 9001)
    x = np.cos(theta)
    u = 1 / np.sqrt(2)
    series = [(1, 3, "#1e66f5"), (3, 5, "#d20f39")]
    for numerator, denominator, color in series:
        rho = numerator / denominator
        maximum = u / (1 - rho)
        axes[0].plot(
            x,
            baseline(theta, rho) / maximum,
            lw=1.25,
            color=color,
            label=rf"$\rho={numerator}/{denominator}$",
        )
    axes[0].scatter([-u, u], [-1, 1], color="#4c4f69", s=28, zorder=5)
    axes[0].axhline(0, color="#9ca0b0", lw=0.6)
    axes[0].set(
        title="Baseline extrema at a common two-cycle",
        xlabel="$x$",
        ylabel=r"$G_\rho(x)/G_\rho(u)$",
        xlim=(-1.02, 1.02),
        ylim=(-1.12, 1.12),
    )
    axes[0].legend(loc="lower right", framealpha=0.95)

    t = np.geomspace(1e-5, 1e-2, 240)
    for numerator, denominator, color in series:
        deficits = np.array(
            [normalized_deficit(v, numerator, denominator) for v in t]
        )
        if not np.all(deficits > 0):
            raise RuntimeError("Positive deficit expected in illustration domain")
        axes[1].loglog(t, deficits, color=color, lw=1.7)
    alpha = np.log(5 / 3) / np.log(3)
    axes[1].loglog(t, 0.6 * t, "--", color="#1e66f5", lw=1)
    axes[1].loglog(t, 0.30 * t**alpha, "--", color="#d20f39", lw=1)
    axes[1].text(2.3e-4, 5.0e-5, "slope 1", color="#1e66f5", rotation=20)
    axes[1].text(2.0e-5, 0.035, r"slope $\alpha\approx0.465$", color="#d20f39")
    axes[1].set(
        title="Deficit near the maximum",
        xlabel=r"Angular displacement $t>0$",
        ylabel=r"$1-G_\rho(\cos(\pi/4+t))/G_\rho(u)$",
        xlim=(1e-5, 1e-2),
    )
    axes[1].grid(which="major", color="#ccd0da", lw=0.6)
    for ax in axes:
        ax.spines[["top", "right"]].set_visible(False)

    OUTPUT.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(OUTPUT, metadata={"Date": None})
    plt.close(fig)
    OUTPUT.write_text("\n".join(line.rstrip() for line in OUTPUT.read_text().splitlines()) + "\n")
    print(f"Wrote {OUTPUT}")
    print(f"Largest analytic truncation bound: {(3 / 5) ** TERMS / (1 - 3 / 5):.4g}")


if __name__ == "__main__":
    main()