#!/usr/bin/env python3
"""Y-Theory toroidal Mind phase model.

This module keeps three layers separate:
1. triadic phase coordinates for a local relation;
2. N delayed recursive depth phases and their circular closure;
3. a toroidal projection formed by the collective local phase and a whole-return phase.

The torus is not drawn arbitrarily.  It follows from the product of two phase
circles, S^1 x S^1, and can also be obtained by stereographic projection of a
Clifford torus in S^3.
"""
from __future__ import annotations

from dataclasses import asdict, dataclass
from pathlib import Path
from typing import Literal
import argparse
import csv
import json
import math

import numpy as np

TAU = 2.0 * math.pi
PHI = (1.0 + math.sqrt(5.0)) / 2.0


@dataclass(slots=True)
class ModelParams:
    duration: float = 80.0
    dt: float = 0.02
    layers: int = 7
    # Neutral recursive preset: phi is measured, not injected.
    omega_local: float = 1.20
    omega_whole: float = 0.90
    omega_spread: float = 0.005
    coupling_depth: float = 0.20
    coupling_down: float = 0.08
    coupling_up: float = 0.06
    delay_depth: float = 0.12
    delay_down: float = 0.32
    delay_up: float = 0.42
    residue_decay: float = 0.65
    residue_gain: float = 0.045
    major_radius: float = 3.0
    minor_radius: float = 1.0
    weights: Literal["equal", "fibonacci", "reverse_fibonacci"] = "fibonacci"

    def validate(self) -> None:
        if self.layers < 1:
            raise ValueError("layers must be >= 1")
        if self.dt <= 0 or self.duration <= self.dt:
            raise ValueError("duration and dt must be positive, with duration > dt")
        if not (self.major_radius > self.minor_radius > 0):
            raise ValueError("require major_radius > minor_radius > 0")
        for name in ("delay_depth", "delay_down", "delay_up", "residue_decay"):
            if getattr(self, name) < 0:
                raise ValueError(f"{name} must be non-negative")


@dataclass(slots=True)
class Simulation:
    params: ModelParams
    time: np.ndarray
    layer_phase: np.ndarray
    whole_phase: np.ndarray
    layer_residue: np.ndarray
    whole_residue: np.ndarray
    collective_phase: np.ndarray
    coherence: np.ndarray


def normalized_weights(layers: int, mode: str) -> np.ndarray:
    if layers < 1:
        raise ValueError("layers must be >= 1")
    if mode == "equal":
        w = np.ones(layers, dtype=float)
    else:
        fib = [1.0, 1.0]
        while len(fib) < layers:
            fib.append(fib[-1] + fib[-2])
        w = np.asarray(fib[:layers], dtype=float)
        if mode == "reverse_fibonacci":
            w = w[::-1]
        elif mode != "fibonacci":
            raise ValueError(f"unknown weight mode: {mode}")
    return w / w.sum()


def triad_from_phase(angle: np.ndarray | float) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Normalized triadic phase chart.

    I + O = 1 and (O-I)^2 + R^2 = 1.
    This is a coordinate representation of one bipolar relational cycle, not a
    claim that every raw inside/outside amplitude in nature sums to one.
    """
    a = np.asarray(angle)
    polarity = np.cos(a)
    relation = np.sin(a)
    inside = 0.5 * (1.0 - polarity)
    outside = 0.5 * (1.0 + polarity)
    return inside, outside, relation


def circular_closure(phases: np.ndarray, weights: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    z = np.sum(weights[None, :] * np.exp(1j * phases), axis=1)
    coherence = np.abs(z)
    safe = np.where(coherence > 1e-12, z / coherence, 1.0 + 0.0j)
    phase = np.unwrap(np.angle(safe))
    return phase, coherence


def _delay_index(i: int, delay: float, dt: float) -> int:
    """Return the sampled history index for a non-negative delay.

    A zero delay means the current Euler state at index ``i``. Positive delays
    are rounded to the nearest sampled history index and clamped at the initial
    condition. This makes the UI label ``delay = 0`` semantically exact rather
    than silently inserting one time step.
    """
    # Half-up rounding matches JavaScript Math.round for non-negative delays.
    steps = int(math.floor(delay / dt + 0.5))
    return max(0, i - steps)


def initial_layer_phases(layers: int) -> np.ndarray:
    """Shared initial-state rule used by the Python and browser runtimes."""
    if layers < 1:
        raise ValueError("layers must be >= 1")
    return np.linspace(0.0, -0.7, layers)


def simulate(params: ModelParams) -> Simulation:
    params.validate()
    n = int(round(params.duration / params.dt)) + 1
    t = np.arange(n, dtype=float) * params.dt
    layers = params.layers

    theta = np.zeros((n, layers), dtype=float)
    phi = np.zeros(n, dtype=float)
    h = np.zeros((n, layers), dtype=float)
    h_whole = np.zeros(n, dtype=float)

    # The stagger is an initial history, not a hardcoded final geometry.
    theta[0, :] = initial_layer_phases(layers)
    phi[0] = 0.25
    omegas = params.omega_local + params.omega_spread * np.linspace(-1.0, 1.0, layers)

    for i in range(n - 1):
        jd = _delay_index(i, params.delay_depth, params.dt)
        jdown = _delay_index(i, params.delay_down, params.dt)
        jup = _delay_index(i, params.delay_up, params.dt)

        mismatch = np.empty(layers, dtype=float)
        mismatch[0] = math.sin(phi[jdown] - theta[i, 0])
        for layer in range(1, layers):
            mismatch[layer] = math.sin(theta[jd, layer - 1] - theta[i, layer])

        h[i + 1, :] = h[i, :] + params.dt * (
            -params.residue_decay * h[i, :] + mismatch
        )

        dtheta = np.empty(layers, dtype=float)
        dtheta[0] = (
            omegas[0]
            + params.coupling_down * mismatch[0]
            + params.residue_gain * h[i, 0]
        )
        for layer in range(1, layers):
            dtheta[layer] = (
                omegas[layer]
                + params.coupling_depth * mismatch[layer]
                + params.residue_gain * h[i, layer]
            )

        whole_mismatch = math.sin(theta[jup, -1] - phi[i])
        h_whole[i + 1] = h_whole[i] + params.dt * (
            -params.residue_decay * h_whole[i] + whole_mismatch
        )
        dphi = (
            params.omega_whole
            + params.coupling_up * whole_mismatch
            + params.residue_gain * h_whole[i]
        )

        theta[i + 1, :] = theta[i, :] + params.dt * dtheta
        phi[i + 1] = phi[i] + params.dt * dphi

    weights = normalized_weights(layers, params.weights)
    collective, coherence = circular_closure(theta, weights)
    return Simulation(params, t, theta, phi, h, h_whole, collective, coherence)


def torus_standard(theta: np.ndarray, phi: np.ndarray, major: float, minor: float) -> np.ndarray:
    theta = np.asarray(theta)
    phi = np.asarray(phi)
    x = (major + minor * np.cos(theta)) * np.cos(phi)
    y = (major + minor * np.cos(theta)) * np.sin(phi)
    z = minor * np.sin(theta)
    return np.stack((x, y, z), axis=-1)


def torus_from_triads(
    local_phase: np.ndarray,
    whole_phase: np.ndarray,
    major: float,
    minor: float,
) -> np.ndarray:
    i_l, o_l, r_l = triad_from_phase(local_phase)
    i_h, o_h, r_h = triad_from_phase(whole_phase)
    p_l = o_l - i_l
    p_h = o_h - i_h
    x = (major + minor * p_l) * p_h
    y = (major + minor * p_l) * r_h
    z = minor * r_l
    return np.stack((x, y, z), axis=-1)


def triads_from_torus(points: np.ndarray, major: float, minor: float) -> dict[str, np.ndarray]:
    """Invert a ring-torus point into local and whole normalized triads.

    Requires major > minor > 0 so rho never vanishes on the ring torus.
    Returns component arrays with the same leading shape as ``points``.
    """
    if not (major > minor > 0):
        raise ValueError("require major > minor > 0 for an everywhere-invertible ring-torus chart")
    pts = np.asarray(points, dtype=float)
    if pts.shape[-1] != 3:
        raise ValueError("points must have final dimension 3")
    x, y, z = pts[..., 0], pts[..., 1], pts[..., 2]
    rho = np.sqrt(x * x + y * y)
    if np.any(rho <= 0):
        raise ValueError("rho must remain positive for whole-phase recovery")
    p_l = (rho - major) / minor
    r_l = z / minor
    p_h = x / rho
    r_h = y / rho
    return {
        "local_inside": 0.5 * (1.0 - p_l),
        "local_outside": 0.5 * (1.0 + p_l),
        "local_relation": r_l,
        "whole_inside": 0.5 * (1.0 - p_h),
        "whole_outside": 0.5 * (1.0 + p_h),
        "whole_relation": r_h,
        "local_polarity": p_l,
        "whole_polarity": p_h,
        "local_phase": np.arctan2(r_l, p_l),
        "whole_phase": np.arctan2(r_h, p_h),
    }


def clifford_4d(local_phase: np.ndarray, whole_phase: np.ndarray, a: float = 1 / math.sqrt(2)) -> np.ndarray:
    if not 0 < a < 1:
        raise ValueError("a must lie in (0, 1)")
    b = math.sqrt(1.0 - a * a)
    q1 = a * np.cos(whole_phase)
    q2 = a * np.sin(whole_phase)
    q3 = b * np.cos(local_phase)
    q4 = b * np.sin(local_phase)
    return np.stack((q1, q2, q3, q4), axis=-1)


def rotation4d(angle_13: float = 0.0, angle_24: float = 0.0, angle_14: float = 0.0) -> np.ndarray:
    M = np.eye(4)
    for i, j, angle in ((0, 2, angle_13), (1, 3, angle_24), (0, 3, angle_14)):
        c, s = math.cos(angle), math.sin(angle)
        R = np.eye(4)
        R[i, i] = R[j, j] = c
        R[i, j] = -s
        R[j, i] = s
        M = R @ M
    return M


def stereographic(q: np.ndarray, rotation: np.ndarray | None = None, eps: float = 1e-8) -> np.ndarray:
    q = np.asarray(q, dtype=float)
    if rotation is not None:
        q = np.einsum("ij,...j->...i", rotation, q)
    den = 1.0 - q[..., 3]
    den = np.where(np.abs(den) < eps, np.sign(den) * eps + (den == 0) * eps, den)
    return q[..., :3] / den[..., None]


def phase_distance(a1: np.ndarray, a2: np.ndarray) -> np.ndarray:
    return np.abs(np.angle(np.exp(1j * (a1 - a2))))


def validate(sim: Simulation) -> dict:
    p = sim.params
    I, O, R = triad_from_phase(sim.collective_phase)
    triad_sum_error = float(np.max(np.abs(I + O - 1.0)))
    triad_circle_error = float(np.max(np.abs((O - I) ** 2 + R**2 - 1.0)))

    xyz_direct = torus_standard(sim.collective_phase, sim.whole_phase, p.major_radius, p.minor_radius)
    xyz_triad = torus_from_triads(sim.collective_phase, sim.whole_phase, p.major_radius, p.minor_radius)
    embedding_error = float(np.max(np.abs(xyz_direct - xyz_triad)))
    radial = np.sqrt(xyz_direct[:, 0] ** 2 + xyz_direct[:, 1] ** 2)
    torus_error = float(np.max(np.abs((radial - p.major_radius) ** 2 + xyz_direct[:, 2] ** 2 - p.minor_radius**2)))

    inv = triads_from_torus(xyz_direct, p.major_radius, p.minor_radius)
    i_l, o_l, r_l = triad_from_phase(sim.collective_phase)
    i_h, o_h, r_h = triad_from_phase(sim.whole_phase)
    inverse_triad_error = float(np.max(np.abs(np.column_stack((
        inv["local_inside"] - i_l,
        inv["local_outside"] - o_l,
        inv["local_relation"] - r_l,
        inv["whole_inside"] - i_h,
        inv["whole_outside"] - o_h,
        inv["whole_relation"] - r_h,
    )))))
    xyz_roundtrip = torus_from_triads(
        inv["local_phase"], inv["whole_phase"], p.major_radius, p.minor_radius
    )
    inverse_coordinate_error = float(np.max(np.abs(xyz_roundtrip - xyz_direct)))

    q = clifford_4d(sim.collective_phase, sim.whole_phase)
    s3_error = float(np.max(np.abs(np.sum(q * q, axis=1) - 1.0)))
    proj = stereographic(q)
    a = 1 / math.sqrt(2)
    b = math.sqrt(1 - a * a)
    Rproj, rproj = 1 / a, b / a
    rho = np.sqrt(proj[:, 0] ** 2 + proj[:, 1] ** 2)
    stereo_torus_error = float(np.max(np.abs((rho - Rproj) ** 2 + proj[:, 2] ** 2 - rproj**2)))

    # Approximate winding ratio after discarding transient.
    start = max(1, int(0.15 * len(sim.time)))
    whole_delta = sim.whole_phase[-1] - sim.whole_phase[start]
    local_delta = sim.collective_phase[-1] - sim.collective_phase[start]
    winding_ratio = float(local_delta / whole_delta) if abs(whole_delta) > 1e-12 else math.nan
    phase_lock_value = float(abs(np.mean(np.exp(1j * (sim.collective_phase[start:] - sim.whole_phase[start:])))))

    # Search for a near-return in torus phase while history/residue differs.
    phase_pair = np.column_stack((np.mod(sim.collective_phase, TAU), np.mod(sim.whole_phase, TAU)))
    target = phase_pair[-1]
    earlier = phase_pair[: max(1, len(phase_pair) - int(5.0 / p.dt))]
    dlocal = phase_distance(earlier[:, 0], target[0])
    dwhole = phase_distance(earlier[:, 1], target[1])
    dist = np.hypot(dlocal, dwhole)
    idx = int(np.argmin(dist))
    residue_distance = float(np.linalg.norm(sim.layer_residue[-1] - sim.layer_residue[idx]))
    full_state_phase_distance = float(dist[idx])

    return {
        "triad_sum_max_error": triad_sum_error,
        "triad_circle_max_error": triad_circle_error,
        "triad_to_standard_embedding_max_error": embedding_error,
        "standard_torus_implicit_max_error": torus_error,
        "torus_to_triad_inverse_max_error": inverse_triad_error,
        "torus_inverse_coordinate_roundtrip_max_error": inverse_coordinate_error,
        "clifford_S3_norm_max_error": s3_error,
        "stereographic_torus_implicit_max_error": stereo_torus_error,
        "measured_winding_ratio": winding_ratio,
        "target_golden_ratio": PHI,
        "phase_lock_value_0_to_1": phase_lock_value,
        "near_return_phase_distance": full_state_phase_distance,
        "near_return_residue_distance": residue_distance,
        "near_return_earlier_time": float(sim.time[idx]),
        "near_return_later_time": float(sim.time[-1]),
        "coherence_min": float(sim.coherence.min()),
        "coherence_mean": float(sim.coherence.mean()),
        "coherence_max": float(sim.coherence.max()),
    }


def write_csv(sim: Simulation, path: Path) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    xyz = torus_standard(sim.collective_phase, sim.whole_phase, sim.params.major_radius, sim.params.minor_radius)
    I_l, O_l, R_l = triad_from_phase(sim.collective_phase)
    I_h, O_h, R_h = triad_from_phase(sim.whole_phase)
    with path.open("w", newline="", encoding="utf-8") as f:
        writer = csv.writer(f)
        header = [
            "time", "collective_local_phase", "whole_phase", "coherence",
            "local_inside", "local_outside", "local_relation",
            "whole_inside", "whole_outside", "whole_relation",
            "x", "y", "z",
        ]
        header += [f"layer_{i+1}_phase" for i in range(sim.params.layers)]
        header += [f"layer_{i+1}_residue" for i in range(sim.params.layers)]
        writer.writerow(header)
        for n in range(len(sim.time)):
            row = [
                sim.time[n], sim.collective_phase[n], sim.whole_phase[n], sim.coherence[n],
                I_l[n], O_l[n], R_l[n], I_h[n], O_h[n], R_h[n],
                xyz[n, 0], xyz[n, 1], xyz[n, 2],
            ]
            row += list(sim.layer_phase[n])
            row += list(sim.layer_residue[n])
            writer.writerow(row)


def project_3d(points: np.ndarray, yaw: float = 0.7, pitch: float = 0.45, scale: float = 150.0, cx: float = 400.0, cy: float = 300.0) -> np.ndarray:
    pts = np.asarray(points, dtype=float)
    cyaw, syaw = math.cos(yaw), math.sin(yaw)
    cp, sp = math.cos(pitch), math.sin(pitch)
    Ry = np.array([[cyaw, -syaw, 0], [syaw, cyaw, 0], [0, 0, 1]], dtype=float)
    Rx = np.array([[1, 0, 0], [0, cp, -sp], [0, sp, cp]], dtype=float)
    rot = pts @ (Rx @ Ry).T
    # Orthographic projection preserves the torus invariant in model space.
    out = np.column_stack((cx + scale * rot[:, 0] / 4.5, cy - scale * rot[:, 1] / 4.5, rot[:, 2]))
    return out


def write_svg(sim: Simulation, path: Path, frame_index: int | None = None) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    idx = len(sim.time) - 1 if frame_index is None else int(frame_index)
    idx = max(0, min(idx, len(sim.time) - 1))
    width, height = 1000, 760
    lines: list[str] = []
    lines.append(f'<svg xmlns="http://www.w3.org/2000/svg" width="{width}" height="{height}" viewBox="0 0 {width} {height}">')
    lines.append('<rect width="100%" height="100%" fill="white"/>')
    lines.append('<text x="40" y="48" font-family="sans-serif" font-size="28">Y-Theory Toroidal Mind — equation-derived projection</text>')
    lines.append('<text x="40" y="78" font-family="sans-serif" font-size="15">Each surface address is a pair of triadic phases: local depth-cycle × whole return-cycle.</text>')

    # Mesh curves generated from the equation, not hand-drawn.
    curves = []
    for local in np.linspace(0, TAU, 15, endpoint=False):
        whole = np.linspace(0, TAU, 120)
        pts = torus_standard(np.full_like(whole, local), whole, sim.params.major_radius, sim.params.minor_radius)
        curves.append(pts)
    for whole in np.linspace(0, TAU, 24, endpoint=False):
        local = np.linspace(0, TAU, 90)
        pts = torus_standard(local, np.full_like(local, whole), sim.params.major_radius, sim.params.minor_radius)
        curves.append(pts)
    for pts in curves:
        p2 = project_3d(pts, cx=500, cy=390, scale=330)
        coords = " ".join(f"{x:.2f},{y:.2f}" for x, y, _ in p2)
        lines.append(f'<polyline points="{coords}" fill="none" stroke="currentColor" stroke-opacity="0.18" stroke-width="0.8"/>')

    start = max(0, idx - 900)
    traj = torus_standard(sim.collective_phase[start:idx+1], sim.whole_phase[start:idx+1], sim.params.major_radius, sim.params.minor_radius)
    p2 = project_3d(traj, cx=500, cy=390, scale=330)
    coords = " ".join(f"{x:.2f},{y:.2f}" for x, y, _ in p2)
    lines.append(f'<polyline points="{coords}" fill="none" stroke="currentColor" stroke-width="2.4"/>')

    layer_pts = torus_standard(sim.layer_phase[idx], np.full(sim.params.layers, sim.whole_phase[idx]), sim.params.major_radius, sim.params.minor_radius)
    layer2 = project_3d(layer_pts, cx=500, cy=390, scale=330)
    for i, (x, y, _) in enumerate(layer2):
        lines.append(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="7" fill="white" stroke="currentColor" stroke-width="2"/>')
        lines.append(f'<text x="{x+9:.2f}" y="{y-9:.2f}" font-family="sans-serif" font-size="13">L{i+1}</text>')

    current = torus_standard(np.array([sim.collective_phase[idx]]), np.array([sim.whole_phase[idx]]), sim.params.major_radius, sim.params.minor_radius)
    cur2 = project_3d(current, cx=500, cy=390, scale=330)[0]
    lines.append(f'<circle cx="{cur2[0]:.2f}" cy="{cur2[1]:.2f}" r="10" fill="currentColor"/>')

    I, O, Rel = triad_from_phase(sim.collective_phase[idx])
    IH, OH, RH = triad_from_phase(sim.whole_phase[idx])
    lines.append('<rect x="40" y="620" width="920" height="105" fill="none" stroke="currentColor" stroke-opacity="0.3"/>')
    text = [
        f't = {sim.time[idx]:.2f}',
        f'Local triad: I={float(I):.3f}, O={float(O):.3f}, R={float(Rel):.3f}',
        f'Whole triad: I={float(IH):.3f}, O={float(OH):.3f}, R={float(RH):.3f}',
        f'Coherence of {sim.params.layers}-layer closure = {sim.coherence[idx]:.3f}',
    ]
    for j, s in enumerate(text):
        lines.append(f'<text x="60" y="{650+j*22}" font-family="monospace" font-size="15">{s}</text>')
    lines.append('</svg>')
    path.write_text("\n".join(lines), encoding="utf-8")


def render_mp4(sim: Simulation, path: Path, seconds: float = 9.0, fps: int = 24) -> None:
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from matplotlib.animation import FuncAnimation, FFMpegWriter

    path.parent.mkdir(parents=True, exist_ok=True)
    frames = max(48, int(seconds * fps))
    indices = np.linspace(max(0, len(sim.time) - int(35 / sim.params.dt)), len(sim.time) - 1, frames).astype(int)

    u = np.linspace(0, TAU, 42)
    v = np.linspace(0, TAU, 22)
    U, V = np.meshgrid(u, v)
    surf = torus_standard(V, U, sim.params.major_radius, sim.params.minor_radius)

    fig = plt.figure(figsize=(12.8, 7.2))
    ax = fig.add_subplot(111, projection="3d")
    ax.plot_wireframe(surf[..., 0], surf[..., 1], surf[..., 2])
    ax.set_box_aspect((1, 1, 0.45))
    ax.set_xlim(-(sim.params.major_radius + sim.params.minor_radius + 0.3), sim.params.major_radius + sim.params.minor_radius + 0.3)
    ax.set_ylim(-(sim.params.major_radius + sim.params.minor_radius + 0.3), sim.params.major_radius + sim.params.minor_radius + 0.3)
    ax.set_zlim(-sim.params.minor_radius - 0.3, sim.params.minor_radius + 0.3)
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")
    ax.set_title(
        f"Y-Theory candidate model: {sim.params.layers} recursive relations "
        "carried through whole-level return"
    )

    trajectory, = ax.plot([], [], [])
    current, = ax.plot([], [], [], marker="o", linestyle="None")
    layer_lines = [ax.plot([], [], [], marker="o", linestyle="None")[0] for _ in range(sim.params.layers)]
    annotation = ax.text2D(0.02, 0.97, "", transform=ax.transAxes, va="top", family="monospace")

    def update(frame_no: int):
        idx = indices[frame_no]
        start = max(0, idx - int(18 / sim.params.dt))
        pts = torus_standard(sim.collective_phase[start:idx+1], sim.whole_phase[start:idx+1], sim.params.major_radius, sim.params.minor_radius)
        trajectory.set_data(pts[:, 0], pts[:, 1])
        trajectory.set_3d_properties(pts[:, 2])
        p = pts[-1]
        current.set_data([p[0]], [p[1]])
        current.set_3d_properties([p[2]])
        lpts = torus_standard(sim.layer_phase[idx], np.full(sim.params.layers, sim.whole_phase[idx]), sim.params.major_radius, sim.params.minor_radius)
        for j, line in enumerate(layer_lines):
            line.set_data([lpts[j, 0]], [lpts[j, 1]])
            line.set_3d_properties([lpts[j, 2]])
        I, O, Rel = triad_from_phase(sim.collective_phase[idx])
        annotation.set_text(
            f"t={sim.time[idx]:6.2f}\n"
            f"local triad I/O/R={float(I):.3f}/{float(O):.3f}/{float(Rel):+.3f}\n"
            f"whole phase={np.mod(sim.whole_phase[idx], TAU):.3f}\n"
            f"{sim.params.layers}-layer coherence={sim.coherence[idx]:.3f}"
        )
        ax.view_init(elev=28, azim=35 + 360 * frame_no / frames)
        return [trajectory, current, annotation, *layer_lines]

    anim = FuncAnimation(fig, update, frames=frames, interval=1000 / fps, blit=False)
    anim.save(path, writer=FFMpegWriter(fps=fps, bitrate=2400))
    plt.close(fig)


def preset(name: str) -> ModelParams:
    if name == "neutral":
        return ModelParams(
            duration=50.0,
            dt=0.02,
            omega_local=1.0,
            omega_whole=1.0,
            omega_spread=0.0,
            coupling_depth=0.0,
            coupling_down=0.0,
            coupling_up=0.0,
            residue_gain=0.0,
            weights="equal",
        )
    if name == "rational":
        return ModelParams(
            duration=8 * math.pi,
            dt=0.015,
            omega_local=1.5,
            omega_whole=1.0,
            omega_spread=0.0,
            coupling_depth=0.0,
            coupling_down=0.0,
            coupling_up=0.0,
            residue_gain=0.0,
            weights="equal",
        )
    if name == "golden":
        return ModelParams(
            duration=80.0,
            dt=0.02,
            omega_local=PHI,
            omega_whole=1.0,
            omega_spread=0.0,
            coupling_depth=0.0,
            coupling_down=0.0,
            coupling_up=0.0,
            residue_gain=0.0,
            weights="equal",
        )
    if name == "recursive":
        # N delayed depth relations close through the whole return.
        # No golden ratio is supplied to this preset.
        return ModelParams()
    raise ValueError(f"unknown preset: {name}")


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--preset", choices=["neutral", "rational", "golden", "recursive"], default="recursive")
    parser.add_argument("--out", type=Path, default=Path("exports"))
    parser.add_argument("--mp4", action="store_true")
    args = parser.parse_args()

    p = preset(args.preset)
    sim = simulate(p)
    args.out.mkdir(parents=True, exist_ok=True)
    write_csv(sim, args.out / f"{args.preset}_trajectory.csv")
    write_svg(sim, args.out / f"{args.preset}_projection.svg")
    report = {"params": asdict(p), "validation": validate(sim)}
    (args.out / f"{args.preset}_validation.json").write_text(json.dumps(report, indent=2), encoding="utf-8")
    if args.mp4:
        render_mp4(sim, args.out / f"{args.preset}_loop.mp4")
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()