"""Independent log-amplitude/Riccati mode-transport reference for AP1-M1.

This module is deliberately restricted to short, nonphysical diagnostic
background prefixes rooted at the chronology-safe 128-shell junction.  It
validates a representation that remains finite through tachyonic growth and a
moving physical-momentum chart whose overlap state is transported rather than
reinitialized.  It does not implement the long AP1-M1 background, UV-tail
entry, backreaction, a physical response kernel, or observables.

Chronology is explicit: a non-tachyonic pilot is evaluated first, a separate
tolerance manifest is then frozen, and only afterwards may disjoint momentum
shells be evaluated on the held-out tachyonic prefixes.
"""
from __future__ import annotations

from dataclasses import asdict, dataclass
from datetime import datetime, timezone
from hashlib import sha256
import json
import math
from pathlib import Path
from typing import Any

import numpy as np
from scipy.integrate import solve_ivp

from ap1_m1_background_multiscale_preflight import _diagnostic_rhs, load_seed
from ap1_r2c_high_precision_modes import propagate_endpoint_decimal
from chi_background_closure import ChiParameters, ChiTrajectory, PV_J, portal_terms
from planck2018_neutrino_closure import Planck2018Pilot


EXPECTED_AUTHORITIES = {
    "AP1/APEIRON_AP1_M1_JUNCTION_128_SEED_CHECKPOINT_LATEST.json": (
        "cf946bfcde4c10a36d5eee99a18fe1a8745483e111304d52f4cc17b3fdf103b4"
    ),
    "AP1/APEIRON_AP1_M1_JUNCTION_128_SEED_LATEST.npz": (
        "644126925fcff32df6f881ee604787751dbe48829072671bf541552bd65dc729"
    ),
    "AP1/APEIRON_AP1_M1_BACKGROUND_MULTISCALE_PREFLIGHT_LATEST.json": (
        "c58ee0221e284dc9a5a892cae68b31138b4f3b9c7634f2570880e4f1086d0d60"
    ),
    "AP1/CODE/ap1_m1_background_multiscale_preflight.py": (
        "1482b548881af15983b3887434b99713fe0547f3a1e23536e001b3e679b43488"
    ),
    "AP1/CODE/ap1_r2c_high_precision_modes.py": (
        "3b76dc1bba339fc7576381661b7a0099501231154f45cf9bb2eba3db16f0d484"
    ),
    "AP1/CODE/chi_background_closure.py": (
        "1097888c72fbe9152e1905a2a3ab37c15df9ca681639571b87d76de6ddf072a7"
    ),
}


FREEZE_POLICY = {
    "comparison_multiplier": 64.0,
    "comparison_floor": 5.0e-5,
    "time_grid_multiplier": 64.0,
    "time_grid_floor": 5.0e-5,
    "rounding_multiplier": 64.0,
    "rounding_floor": 1.0e-13,
    "riccati_refinement_multiplier": 64.0,
    "riccati_refinement_floor": 1.0e-10,
    "chart_overlap_multiplier": 128.0,
    "chart_overlap_floor": 1.0e-9,
    "initial_seed_match_cap": 1.0e-12,
    "wronskian_relative_cap": 1.0e-11,
}


class LogTransportError(RuntimeError):
    """A prerequisite or fail-closed transport condition was violated."""


@dataclass(frozen=True)
class ReferenceConfig:
    pilot_endpoints_relative_N: tuple[float, ...] = (0.01, 0.02, 0.03)
    heldout_endpoints_relative_N: tuple[float, ...] = (0.065, 0.08, 0.10)
    pilot_mode_indices: tuple[int, ...] = (0, 16, 48, 80, 112, 127)
    heldout_mode_indices: tuple[int, ...] = (1, 8, 32, 64, 96, 126)
    coarse_step_N: float = 2.5e-5
    fine_step_N: float = 1.25e-5
    diagnostic_rtol: float = 2.0e-11
    riccati_loose_rtol: float = 5.0e-10
    riccati_fine_rtol: float = 2.0e-11
    decimal_low_digits: int = 60
    decimal_high_digits: int = 80

    def validate(self) -> None:
        for endpoints in (
            self.pilot_endpoints_relative_N,
            self.heldout_endpoints_relative_N,
        ):
            if tuple(sorted(set(endpoints))) != endpoints:
                raise ValueError("endpoints must be unique and strictly increasing")
            if not endpoints or endpoints[0] <= 0.0 or endpoints[-1] > 0.125:
                raise ValueError("only short positive diagnostic prefixes are allowed")
            for endpoint in endpoints:
                intervals = endpoint / self.fine_step_N
                if abs(intervals - round(intervals)) > 1.0e-10:
                    raise ValueError("fine step must exactly tile every endpoint")
        if set(self.pilot_mode_indices) & set(self.heldout_mode_indices):
            raise ValueError("pilot and held-out momentum shells must be disjoint")
        if min(self.pilot_mode_indices + self.heldout_mode_indices) < 0:
            raise ValueError("momentum indices must be non-negative")
        if max(self.pilot_mode_indices + self.heldout_mode_indices) >= 128:
            raise ValueError("momentum indices exceed the chronology-safe seed")
        if not 0.0 < self.fine_step_N < self.coarse_step_N <= 5.0e-5:
            raise ValueError("bounded nested diagnostic steps required")
        ratio = self.coarse_step_N / self.fine_step_N
        if abs(ratio - round(ratio)) > 1.0e-12:
            raise ValueError("coarse and fine steps must be nested")
        if not 0.0 < self.riccati_fine_rtol < self.riccati_loose_rtol <= 1.0e-8:
            raise ValueError("ordered Riccati tolerances required")
        if not 0.0 < self.diagnostic_rtol <= 1.0e-9:
            raise ValueError("bounded diagnostic tolerance required")
        if not 32 <= self.decimal_low_digits < self.decimal_high_digits:
            raise ValueError("ordered independent Decimal precisions required")


def file_sha256(path: Path) -> str:
    digest = sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            digest.update(block)
    return digest.hexdigest()


def load_inputs(apeiron_root: Path) -> tuple[dict, dict[str, np.ndarray], dict[str, str]]:
    observed = {
        name: file_sha256(apeiron_root / name) for name in EXPECTED_AUTHORITIES
    }
    if observed != EXPECTED_AUTHORITIES:
        changed = [
            name for name, expected in EXPECTED_AUTHORITIES.items()
            if observed.get(name) != expected
        ]
        raise LogTransportError(f"authority drift: {changed}")
    preflight = json.loads(
        (apeiron_root / "AP1/APEIRON_AP1_M1_BACKGROUND_MULTISCALE_PREFLIGHT_LATEST.json")
        .read_text(encoding="utf-8")
    )
    if not preflight.get("all_preflight_gates_pass"):
        raise LogTransportError("G26 multiscale architecture preflight is not PASS")
    if preflight.get("new_AP1_M1_background_runs") != 0:
        raise LogTransportError("G26 authority already claims a background run")
    seed_report, arrays, _ = load_seed(apeiron_root)
    return seed_report, arrays, observed


def _intervals(endpoint: float, step: float) -> int:
    value = endpoint / step
    rounded = int(round(value))
    if rounded < 2 or abs(value - rounded) > 1.0e-10:
        raise ValueError("step does not tile endpoint")
    return rounded


def short_diagnostic_trajectory(
    arrays: dict[str, np.ndarray],
    endpoint_relative_N: float,
    step_N: float,
    rtol: float,
) -> ChiTrajectory:
    """Return one short fixed-source diagnostic prefix; never a background."""

    intervals = _intervals(endpoint_relative_N, step_N)
    relative_N = np.linspace(0.0, endpoint_relative_N, intervals + 1)
    state6 = np.asarray(arrays["m1_state_physical_N"], dtype=float)
    initial = state6[:5].copy()
    N_seed = float(state6[5])
    quantum_initial = np.asarray(arrays["quantum_rho_pressure_chi2"], dtype=float)
    pilot = Planck2018Pilot()
    chi = ChiParameters()
    solution = solve_ivp(
        lambda n, y: _diagnostic_rhs(
            n, y, N_seed, quantum_initial, pilot, chi
        ),
        (0.0, endpoint_relative_N),
        initial,
        t_eval=relative_N,
        method="DOP853",
        rtol=rtol,
        atol=np.array([2.0e-13, 2.0e-17, 2.0e-13, 2.0e-17, 2.0e-17]),
        max_step=min(1.0e-3, endpoint_relative_N / 32.0),
    )
    if not solution.success or solution.y.shape != (5, intervals + 1):
        raise LogTransportError(f"short diagnostic envelope failed: {solution.message}")
    rows = solution.y.T
    Hdot = np.empty(intervals + 1)
    for index, (n_relative, row) in enumerate(zip(relative_N, rows)):
        Hdot[index] = (
            _diagnostic_rhs(
                float(n_relative), row, N_seed, quantum_initial, pilot, chi
            )[4]
            * row[4]
        )
    # The chronology-safe capsule owns the junction derivative.  The
    # fixed-source diagnostic RHS is only used immediately after the junction;
    # replacing the stored initial Hdot would silently reinitialize the mode
    # frequency and defeat the seed transport test.
    Hdot[0] = float(np.asarray(arrays["Hdot_Mpl2"], dtype=float)[0])
    return ChiTrajectory(
        N=N_seed + relative_N,
        H=rows[:, 4],
        Hdot=Hdot,
        sigma=rows[:, 0],
        theta=rows[:, 2],
    ).validated()


def prefix_trajectory(trajectory: ChiTrajectory, endpoint_relative_N: float) -> ChiTrajectory:
    tr = trajectory.validated()
    target = float(tr.N[0] + endpoint_relative_N)
    index = int(np.argmin(np.abs(tr.N - target)))
    if abs(float(tr.N[index] - target)) > 2.0e-13:
        raise LogTransportError("requested prefix endpoint is not on the grid")
    sl = slice(0, index + 1)
    return ChiTrajectory(
        N=tr.N[sl], H=tr.H[sl], Hdot=tr.Hdot[sl],
        sigma=tr.sigma[sl], theta=tr.theta[sl],
    ).validated()


def seed_comoving_momenta(
    arrays: dict[str, np.ndarray], indices: tuple[int, ...]
) -> tuple[np.ndarray, np.ndarray]:
    """Map seed-local physical labels to present-normalized comoving momenta."""

    q_seed = np.asarray(arrays["k_comoving_Mpl"], dtype=float)[list(indices)]
    N_seed = float(np.asarray(arrays["m1_state_physical_N"], dtype=float)[5])
    k_comoving = q_seed * math.exp(N_seed)
    recovered = k_comoving * math.exp(-N_seed)
    if not np.allclose(recovered, q_seed, rtol=4.0e-15, atol=0.0):
        raise LogTransportError("physical/comoving seed momentum map lost precision")
    return q_seed, k_comoving


def seed_mode_state(
    arrays: dict[str, np.ndarray], indices: tuple[int, ...]
) -> dict[str, np.ndarray]:
    selection = np.asarray(indices, dtype=int)
    u = (
        np.asarray(arrays["u_real"], dtype=float)[:, selection]
        + 1j * np.asarray(arrays["u_imag"], dtype=float)[:, selection]
    )
    v = (
        np.asarray(arrays["v_real"], dtype=float)[:, selection]
        + 1j * np.asarray(arrays["v_imag"], dtype=float)[:, selection]
    )
    target = np.asarray(arrays["wronskian_target"], dtype=float)[:, selection]
    if np.any(np.abs(u) == 0.0) or np.any(target <= 0.0):
        raise LogTransportError("nonzero seed modes and positive Wronskians required")
    ratio = v / u
    observed = -2.0 * np.imag(np.conj(u) * v)
    if not np.allclose(observed, target, rtol=1.0e-12, atol=1.0e-14):
        raise LogTransportError("selected seed Wronskian mismatch")
    if not np.all(np.imag(ratio) < 0.0):
        raise LogTransportError("negative canonical Riccati imaginary part required")
    return {
        "log_amplitude": np.log(np.abs(u)),
        "riccati_real": np.real(ratio),
        "log_abs_riccati_imag": np.log(-np.imag(ratio)),
        "phase": np.angle(u),
        "wronskian": target,
    }


def _omega2_independent(
    N_value: float,
    trajectory: ChiTrajectory,
    k_comoving: np.ndarray,
) -> np.ndarray:
    tr = trajectory
    sigma = float(np.interp(N_value, tr.N, tr.sigma))
    theta = float(np.interp(N_value, tr.N, tr.theta))
    H = float(np.interp(N_value, tr.N, tr.H))
    Hdot = float(np.interp(N_value, tr.N, tr.Hdot))
    p = ChiParameters()
    mass2 = float(portal_terms(sigma, theta, p)[0])
    muR2 = (p.muR_over_Lambda * p.Lambda) ** 2
    physical2 = k_comoving * k_comoving * math.exp(-2.0 * N_value)
    return (
        physical2[None, :]
        + mass2
        + PV_J[:, None] * muR2
        - 2.25 * H * H
        - 1.5 * Hdot
    )


def transport_log_riccati(
    trajectory: ChiTrajectory,
    k_comoving: np.ndarray,
    initial: dict[str, np.ndarray],
    rtol: float,
) -> dict[str, Any]:
    """Transport u through L=log|u| and R=v/u without amplitude overflow."""

    tr = trajectory.validated()
    k = np.asarray(k_comoving, dtype=float)
    shape = (4, len(k))
    fields = {
        key: np.asarray(initial[key], dtype=float)
        for key in (
            "log_amplitude",
            "riccati_real",
            "log_abs_riccati_imag",
            "phase",
            "wronskian",
        )
    }
    if any(value.shape != shape for value in fields.values()):
        raise ValueError(f"initial log transport arrays must have shape {shape}")
    if not np.all(fields["wronskian"] > 0.0):
        raise ValueError("positive Wronskian targets required")
    if not 0.0 < rtol <= 1.0e-8:
        raise ValueError("bounded log-Riccati tolerance required")
    count = int(np.prod(shape))
    log_w = np.log(fields["wronskian"].reshape(-1))
    y0 = np.concatenate(
        [
            fields["log_amplitude"].reshape(-1),
            fields["riccati_real"].reshape(-1),
            fields["log_abs_riccati_imag"].reshape(-1),
            fields["phase"].reshape(-1),
        ]
    )

    def rhs(N_value: float, state: np.ndarray) -> np.ndarray:
        L = state[:count]
        x = state[count:2 * count]
        log_abs_beta = state[2 * count:3 * count]
        H = float(np.interp(N_value, tr.N, tr.H))
        omega2 = _omega2_independent(N_value, tr, k).reshape(-1)
        if np.any(log_abs_beta > 0.5 * math.log(np.finfo(float).max)):
            raise LogTransportError("Riccati imaginary component exceeded Float64")
        beta = -np.exp(np.maximum(log_abs_beta, math.log(np.nextafter(0.0, 1.0))))
        beta2 = np.exp(
            np.maximum(2.0 * log_abs_beta, math.log(np.nextafter(0.0, 1.0)))
        )
        return np.concatenate(
            [
                x / H,
                (-omega2 - x * x + beta2) / H,
                -2.0 * x / H,
                beta / H,
            ]
        )

    solution = solve_ivp(
        rhs,
        (float(tr.N[0]), float(tr.N[-1])),
        y0,
        method="DOP853",
        rtol=rtol,
        atol=np.full(4 * count, min(1.0e-12, rtol * 1.0e-2)),
        max_step=min(2.5e-4, float(tr.N[-1] - tr.N[0]) / 32.0),
    )
    if not solution.success or not np.all(np.isfinite(solution.y[:, -1])):
        raise LogTransportError(f"log-Riccati transport failed: {solution.message}")
    endpoint = solution.y[:, -1]
    L = endpoint[:count].reshape(shape)
    x = endpoint[count:2 * count].reshape(shape)
    log_abs_beta = endpoint[2 * count:3 * count].reshape(shape)
    phase = endpoint[3 * count:].reshape(shape)
    beta = -np.exp(
        np.maximum(log_abs_beta, math.log(np.nextafter(0.0, 1.0)))
    )
    reconstructed_relative_error = 0.0
    if float(np.max(L)) < 300.0 and float(np.min(L)) > -300.0:
        u = np.exp(L) * np.exp(1j * phase)
        v = (x + 1j * beta) * u
        observed = -2.0 * np.imag(np.conj(u) * v)
        reconstructed_relative_error = float(
            np.max(np.abs(observed / fields["wronskian"] - 1.0))
        )
    identity = math.log(2.0) + log_abs_beta + 2.0 * L - np.log(
        fields["wronskian"]
    )
    log_wronskian_relative_error = float(
        np.max(np.abs(np.expm1(np.clip(identity, -1.0, 1.0))))
    )
    return {
        "log_amplitude": L,
        "riccati_real": x,
        "log_abs_riccati_imag": log_abs_beta,
        "phase": phase,
        "wronskian": fields["wronskian"].copy(),
        "log_wronskian_identity_error": float(np.max(np.abs(identity))),
        "log_wronskian_relative_error": log_wronskian_relative_error,
        "reconstructed_wronskian_relative_error": reconstructed_relative_error,
        "nfev": int(solution.nfev),
    }


def rebase_trajectory(trajectory: ChiTrajectory) -> ChiTrajectory:
    tr = trajectory.validated()
    return ChiTrajectory(
        N=tr.N - tr.N[0], H=tr.H, Hdot=tr.Hdot,
        sigma=tr.sigma, theta=tr.theta,
    ).validated()


def chart_overlap_error(
    trajectory: ChiTrajectory,
    k_comoving: np.ndarray,
    initial: dict[str, np.ndarray],
    rtol: float,
) -> dict[str, float]:
    """Compare direct transport with an exact moving-chart state handoff."""

    tr = trajectory.validated()
    middle = (len(tr.N) - 1) // 2
    if middle < 1 or middle >= len(tr.N) - 1:
        raise ValueError("overlap trajectory requires an interior handoff")
    direct = transport_log_riccati(tr, k_comoving, initial, rtol)
    left = ChiTrajectory(
        N=tr.N[:middle + 1], H=tr.H[:middle + 1],
        Hdot=tr.Hdot[:middle + 1], sigma=tr.sigma[:middle + 1],
        theta=tr.theta[:middle + 1],
    ).validated()
    right_physical = ChiTrajectory(
        N=tr.N[middle:], H=tr.H[middle:], Hdot=tr.Hdot[middle:],
        sigma=tr.sigma[middle:], theta=tr.theta[middle:],
    ).validated()
    left_result = transport_log_riccati(left, k_comoving, initial, rtol)
    transported = {
        key: left_result[key]
        for key in (
            "log_amplitude",
            "riccati_real",
            "log_abs_riccati_imag",
            "phase",
            "wronskian",
        )
    }
    N_handoff = float(right_physical.N[0])
    physical_at_handoff = k_comoving * math.exp(-N_handoff)
    right_local = rebase_trajectory(right_physical)
    rebased = transport_log_riccati(
        right_local, physical_at_handoff, transported, rtol
    )
    phase_delta = np.angle(
        np.exp(1j * (rebased["phase"] - direct["phase"]))
    )
    return {
        "max_abs_log_amplitude": float(
            np.max(np.abs(rebased["log_amplitude"] - direct["log_amplitude"]))
        ),
        "max_abs_riccati_real_scaled": float(
            np.max(np.abs(rebased["riccati_real"] - direct["riccati_real"]))
            / max(float(np.max(np.abs(direct["riccati_real"]))), 1.0e-30)
        ),
        "max_abs_log_riccati_imag": float(
            np.max(
                np.abs(
                    rebased["log_abs_riccati_imag"]
                    - direct["log_abs_riccati_imag"]
                )
            )
        ),
        "max_abs_phase_wrapped": float(np.max(np.abs(phase_delta))),
        "state_reinitialized_at_handoff": False,
    }


def _max_relative(left: np.ndarray, right: np.ndarray) -> float:
    left = np.asarray(left, dtype=float)
    right = np.asarray(right, dtype=float)
    return float(
        np.max(np.abs(left - right) / np.maximum(np.abs(right), 1.0e-300))
    )


def _window_metrics(
    fine: ChiTrajectory,
    coarse: ChiTrajectory,
    arrays: dict[str, np.ndarray],
    indices: tuple[int, ...],
    config: ReferenceConfig,
) -> dict[str, Any]:
    q_seed, k_comoving = seed_comoving_momenta(arrays, indices)
    initial = seed_mode_state(arrays, indices)
    seed_omega2 = np.asarray(arrays["omega2_Mpl2"], dtype=float)[:, list(indices)]
    # The canonical seed state uses the instantaneous positive vacuum
    # frequency k_phys^2+m^2+M_PV^2.  The -9H^2/4-3Hdot/2 terms enter the
    # subsequent u=a^(3/2)chi propagation but are not part of that stored
    # initialization frequency.
    p = ChiParameters()
    mass2_seed = float(portal_terms(fine.sigma[0], fine.theta[0], p)[0])
    muR2 = (p.muR_over_Lambda * p.Lambda) ** 2
    independent_seed_omega2 = (
        k_comoving[None, :] ** 2 * math.exp(-2.0 * float(fine.N[0]))
        + mass2_seed
        + PV_J[:, None] * muR2
    )
    seed_frequency_relative_error = _max_relative(
        independent_seed_omega2, seed_omega2
    )
    decimal_fine_high = propagate_endpoint_decimal(
        fine, k_comoving, digits=config.decimal_high_digits
    )
    decimal_fine_low = propagate_endpoint_decimal(
        fine, k_comoving, digits=config.decimal_low_digits
    )
    decimal_coarse_high = propagate_endpoint_decimal(
        coarse, k_comoving, digits=config.decimal_high_digits
    )
    riccati_fine = transport_log_riccati(
        fine, k_comoving, initial, config.riccati_fine_rtol
    )
    riccati_loose = transport_log_riccati(
        fine, k_comoving, initial, config.riccati_loose_rtol
    )
    decimal_amp = np.asarray(
        decimal_fine_high["all_sector_endpoint_abs_u"], dtype=float
    )
    decimal_low_amp = np.asarray(
        decimal_fine_low["all_sector_endpoint_abs_u"], dtype=float
    )
    decimal_coarse_amp = np.asarray(
        decimal_coarse_high["all_sector_endpoint_abs_u"], dtype=float
    )
    log_decimal = np.log(decimal_amp)
    omega2_samples = np.stack(
        [_omega2_independent(float(N), fine, k_comoving) for N in fine.N]
    )
    chart = chart_overlap_error(
        fine, k_comoving, initial, config.riccati_fine_rtol
    )
    physical_end = k_comoving * math.exp(-float(fine.N[-1]))
    direct_end = q_seed * math.exp(-float(fine.N[-1] - fine.N[0]))
    return {
        "endpoint_relative_N": float(fine.N[-1] - fine.N[0]),
        "mode_indices": list(indices),
        "mode_count_all_PV_sectors": int(4 * len(indices)),
        "seed_frequency_relative_error": seed_frequency_relative_error,
        "physical_momentum_transport_relative_error": _max_relative(
            physical_end, direct_end
        ),
        "min_physical_sector_omega2_Mpl2": float(
            np.min(omega2_samples[:, 0, :])
        ),
        "physical_sector_tachyonic_sampled": bool(
            np.any(omega2_samples[:, 0, :] < 0.0)
        ),
        "max_abs_log_amplitude_vs_decimal80": float(
            np.max(np.abs(riccati_fine["log_amplitude"] - log_decimal))
        ),
        "decimal_time_grid_relative_delta": _max_relative(
            decimal_coarse_amp, decimal_amp
        ),
        "decimal_60_to_80_relative_delta": _max_relative(
            decimal_low_amp, decimal_amp
        ),
        "riccati_loose_to_fine_log_amplitude_delta": float(
            np.max(
                np.abs(
                    riccati_loose["log_amplitude"]
                    - riccati_fine["log_amplitude"]
                )
            )
        ),
        "log_wronskian_identity_error": float(
            riccati_fine["log_wronskian_identity_error"]
        ),
        "log_wronskian_relative_error": float(
            riccati_fine["log_wronskian_relative_error"]
        ),
        "reconstructed_wronskian_relative_error": float(
            riccati_fine["reconstructed_wronskian_relative_error"]
        ),
        "decimal_wronskian_relative_error": float(
            decimal_fine_high["wronskian_relative_error"]
        ),
        "chart_overlap": chart,
        "riccati_nfev": int(riccati_fine["nfev"]),
    }


def _aggregate(windows: list[dict[str, Any]]) -> dict[str, float | int]:
    return {
        "windows": len(windows),
        "max_seed_frequency_relative_error": max(
            item["seed_frequency_relative_error"] for item in windows
        ),
        "max_physical_momentum_transport_relative_error": max(
            item["physical_momentum_transport_relative_error"] for item in windows
        ),
        "max_abs_log_amplitude_vs_decimal80": max(
            item["max_abs_log_amplitude_vs_decimal80"] for item in windows
        ),
        "max_decimal_time_grid_relative_delta": max(
            item["decimal_time_grid_relative_delta"] for item in windows
        ),
        "max_decimal_60_to_80_relative_delta": max(
            item["decimal_60_to_80_relative_delta"] for item in windows
        ),
        "max_riccati_refinement_log_amplitude_delta": max(
            item["riccati_loose_to_fine_log_amplitude_delta"] for item in windows
        ),
        "max_chart_overlap_error": max(
            max(item["chart_overlap"][key] for key in (
                "max_abs_log_amplitude",
                "max_abs_riccati_real_scaled",
                "max_abs_log_riccati_imag",
                "max_abs_phase_wrapped",
            ))
            for item in windows
        ),
        "max_reconstructed_wronskian_relative_error": max(
            item["reconstructed_wronskian_relative_error"] for item in windows
        ),
        "max_log_wronskian_relative_error": max(
            item["log_wronskian_relative_error"] for item in windows
        ),
        "max_decimal_wronskian_relative_error": max(
            item["decimal_wronskian_relative_error"] for item in windows
        ),
        "tachyonic_windows": sum(
            bool(item["physical_sector_tachyonic_sampled"]) for item in windows
        ),
    }


def evaluate_windows(
    arrays: dict[str, np.ndarray],
    endpoints: tuple[float, ...],
    indices: tuple[int, ...],
    config: ReferenceConfig,
) -> tuple[list[dict[str, Any]], dict[str, float | int]]:
    maximum = max(endpoints)
    fine_full = short_diagnostic_trajectory(
        arrays, maximum, config.fine_step_N, config.diagnostic_rtol
    )
    coarse_full = short_diagnostic_trajectory(
        arrays, maximum, config.coarse_step_N, config.diagnostic_rtol
    )
    windows = []
    for endpoint in endpoints:
        windows.append(
            _window_metrics(
                prefix_trajectory(fine_full, endpoint),
                prefix_trajectory(coarse_full, endpoint),
                arrays,
                indices,
                config,
            )
        )
    return windows, _aggregate(windows)


def build_pilot(
    apeiron_root: Path,
    code_path: Path,
    config: ReferenceConfig | None = None,
) -> dict[str, Any]:
    config = ReferenceConfig() if config is None else config
    config.validate()
    seed_report, arrays, authorities = load_inputs(apeiron_root)
    windows, aggregate = evaluate_windows(
        arrays,
        config.pilot_endpoints_relative_N,
        config.pilot_mode_indices,
        config,
    )
    gates = {
        "chronology_safe_seed_hash_bound": True,
        "pilot_and_heldout_shells_disjoint": not bool(
            set(config.pilot_mode_indices) & set(config.heldout_mode_indices)
        ),
        "pilot_windows_pre_tachyonic": aggregate["tachyonic_windows"] == 0,
        "seed_frequency_map_below_predeclared_cap": (
            aggregate["max_seed_frequency_relative_error"]
            <= FREEZE_POLICY["initial_seed_match_cap"]
        ),
        "moving_momentum_map_below_predeclared_cap": (
            aggregate["max_physical_momentum_transport_relative_error"]
            <= FREEZE_POLICY["initial_seed_match_cap"]
        ),
        "finite_pilot_metrics": all(
            np.isfinite(float(value)) for key, value in aggregate.items()
            if key != "windows"
        ),
        "no_background_candidate_or_kernel": True,
    }
    complete = bool(all(gates.values()))
    return {
        "schema": "apeiron-ap1-m1-log-amplitude-transport-pilot-v1.0",
        "updated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        "classification": (
            "LOG_AMPLITUDE_TRANSPORT_PILOT_COMPLETE_TOLERANCE_FREEZE_OPEN"
            if complete else "LOG_AMPLITUDE_TRANSPORT_PILOT_INCOMPLETE"
        ),
        "authority_sha256": {
            **authorities,
            "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py": file_sha256(
                code_path
            ),
        },
        "seed_classification": seed_report["classification"],
        "config": asdict(config),
        "predeclared_freeze_policy": FREEZE_POLICY,
        "window_kind": "seed_anchored_nested_overlapping_prefixes",
        "pilot_windows": windows,
        "pilot_aggregate": aggregate,
        "gates": gates,
        "all_pilot_gates_pass": complete,
        "checkpoint_eligible": False,
        "heldout_evaluated": False,
        "new_AP1_M1_background_runs": 0,
        "diagnostic_envelope_integrations": 2,
        "trajectory_rows_persisted": 0,
        "background_candidate_assessed": False,
        "physical_response_kernel_started": False,
        "seed_released": False,
        "claim_boundary": (
            "non-tachyonic numerical tolerance pilot only; not a checkpoint, "
            "background, physical mode history, kernel, curve, fit or significance"
        ),
        "next_required": (
            "freeze pilot-derived tolerances before evaluating disjoint held-out "
            "momentum shells on the registered tachyonic prefixes"
        ),
    }


def build_freeze(
    pilot_path: Path,
    code_path: Path,
) -> dict[str, Any]:
    pilot = json.loads(pilot_path.read_text(encoding="utf-8"))
    if not pilot.get("all_pilot_gates_pass") or pilot.get("heldout_evaluated"):
        raise LogTransportError("complete pre-heldout pilot required")
    if pilot["authority_sha256"].get(
        "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py"
    ) != file_sha256(code_path):
        raise LogTransportError("transport code changed after pilot")
    aggregate = pilot["pilot_aggregate"]
    thresholds = {
        "max_abs_log_amplitude_vs_decimal80": max(
            FREEZE_POLICY["comparison_floor"],
            FREEZE_POLICY["comparison_multiplier"]
            * float(aggregate["max_abs_log_amplitude_vs_decimal80"]),
        ),
        "max_decimal_time_grid_relative_delta": max(
            FREEZE_POLICY["time_grid_floor"],
            FREEZE_POLICY["time_grid_multiplier"]
            * float(aggregate["max_decimal_time_grid_relative_delta"]),
        ),
        "max_decimal_60_to_80_relative_delta": max(
            FREEZE_POLICY["rounding_floor"],
            FREEZE_POLICY["rounding_multiplier"]
            * float(aggregate["max_decimal_60_to_80_relative_delta"]),
        ),
        "max_riccati_refinement_log_amplitude_delta": max(
            FREEZE_POLICY["riccati_refinement_floor"],
            FREEZE_POLICY["riccati_refinement_multiplier"]
            * float(aggregate["max_riccati_refinement_log_amplitude_delta"]),
        ),
        "max_chart_overlap_error": max(
            FREEZE_POLICY["chart_overlap_floor"],
            FREEZE_POLICY["chart_overlap_multiplier"]
            * float(aggregate["max_chart_overlap_error"]),
        ),
        "max_initial_seed_or_momentum_map_relative_error": (
            FREEZE_POLICY["initial_seed_match_cap"]
        ),
        "max_wronskian_relative_error": FREEZE_POLICY[
            "wronskian_relative_cap"
        ],
        "required_tachyonic_windows": len(
            pilot["config"]["heldout_endpoints_relative_N"]
        ),
    }
    if not all(
        np.isfinite(float(value)) and float(value) > 0.0
        for value in thresholds.values()
    ):
        raise LogTransportError("pilot-derived thresholds are invalid")
    return {
        "schema": "apeiron-ap1-m1-log-amplitude-transport-freeze-v1.0",
        "updated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        "classification": "LOG_AMPLITUDE_TRANSPORT_TOLERANCES_FROZEN_HELDOUT_UNSEEN",
        "input_sha256": {
            pilot_path.name: file_sha256(pilot_path),
            "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py": file_sha256(
                code_path
            ),
        },
        "pilot_updated_utc": pilot["updated_utc"],
        "predeclared_freeze_policy": pilot["predeclared_freeze_policy"],
        "heldout_plan": {
            "endpoints_relative_N": pilot["config"][
                "heldout_endpoints_relative_N"
            ],
            "mode_indices": pilot["config"]["heldout_mode_indices"],
            "window_kind": pilot["window_kind"],
            "decimal_high_digits": pilot["config"]["decimal_high_digits"],
            "fine_step_N": pilot["config"]["fine_step_N"],
        },
        "thresholds": thresholds,
        "all_tolerance_freeze_gates_pass": True,
        "tolerances_frozen_before_heldout": True,
        "heldout_evaluated": False,
        "checkpoint_eligible": False,
        "new_AP1_M1_background_runs": 0,
        "background_candidate_assessed": False,
        "physical_response_kernel_started": False,
        "claim_boundary": (
            "pilot-derived numerical reference tolerances only; held-out "
            "tachyonic windows, long background and physical kernel remain unseen"
        ),
        "next_required": (
            "evaluate the registered held-out tachyonic windows once and release "
            "no checkpoint unless every frozen gate passes"
        ),
    }


def build_heldout(
    apeiron_root: Path,
    pilot_path: Path,
    freeze_path: Path,
    code_path: Path,
) -> dict[str, Any]:
    pilot = json.loads(pilot_path.read_text(encoding="utf-8"))
    freeze = json.loads(freeze_path.read_text(encoding="utf-8"))
    code_hash = file_sha256(code_path)
    expected_inputs = {
        pilot_path.name: file_sha256(pilot_path),
        "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py": code_hash,
    }
    if freeze.get("input_sha256") != expected_inputs:
        raise LogTransportError("held-out freeze inputs changed")
    if not freeze.get("tolerances_frozen_before_heldout") or freeze.get(
        "heldout_evaluated"
    ):
        raise LogTransportError("unseen held-out state required")
    if not pilot.get("all_pilot_gates_pass"):
        raise LogTransportError("pilot is not complete")
    config = ReferenceConfig(**{
        key: tuple(value) if isinstance(value, list) else value
        for key, value in pilot["config"].items()
    })
    config.validate()
    seed_report, arrays, authorities = load_inputs(apeiron_root)
    windows, aggregate = evaluate_windows(
        arrays,
        config.heldout_endpoints_relative_N,
        config.heldout_mode_indices,
        config,
    )
    thresholds = freeze["thresholds"]
    max_wronskian = max(
        float(aggregate["max_log_wronskian_relative_error"]),
        float(aggregate["max_decimal_wronskian_relative_error"]),
    )
    gates = {
        "freeze_hash_bound_before_heldout": True,
        "pilot_and_heldout_shells_disjoint": not bool(
            set(config.pilot_mode_indices) & set(config.heldout_mode_indices)
        ),
        "all_registered_windows_tachyonic": (
            int(aggregate["tachyonic_windows"])
            == int(thresholds["required_tachyonic_windows"])
        ),
        "seed_frequency_and_momentum_map_below_frozen_cap": max(
            float(aggregate["max_seed_frequency_relative_error"]),
            float(aggregate["max_physical_momentum_transport_relative_error"]),
        ) <= float(thresholds["max_initial_seed_or_momentum_map_relative_error"]),
        "log_amplitude_matches_decimal_below_frozen_cap": (
            float(aggregate["max_abs_log_amplitude_vs_decimal80"])
            <= float(thresholds["max_abs_log_amplitude_vs_decimal80"])
        ),
        "decimal_time_grid_delta_below_frozen_cap": (
            float(aggregate["max_decimal_time_grid_relative_delta"])
            <= float(thresholds["max_decimal_time_grid_relative_delta"])
        ),
        "decimal_rounding_delta_below_frozen_cap": (
            float(aggregate["max_decimal_60_to_80_relative_delta"])
            <= float(thresholds["max_decimal_60_to_80_relative_delta"])
        ),
        "riccati_refinement_below_frozen_cap": (
            float(aggregate["max_riccati_refinement_log_amplitude_delta"])
            <= float(thresholds["max_riccati_refinement_log_amplitude_delta"])
        ),
        "moving_chart_overlap_below_frozen_cap": (
            float(aggregate["max_chart_overlap_error"])
            <= float(thresholds["max_chart_overlap_error"])
        ),
        "wronskian_below_existing_frozen_cap": max_wronskian
        <= float(thresholds["max_wronskian_relative_error"]),
        "no_handoff_vacuum_reset": all(
            not bool(window["chart_overlap"]["state_reinitialized_at_handoff"])
            for window in windows
        ),
        "no_background_candidate_or_kernel": True,
    }
    passed = bool(all(gates.values()))
    return {
        "schema": "apeiron-ap1-m1-log-amplitude-transport-heldout-v1.0",
        "updated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        "classification": (
            "LOG_AMPLITUDE_MOVING_CHART_REFERENCE_PASS_LONG_BACKGROUND_BLOCKED"
            if passed else "LOG_AMPLITUDE_MOVING_CHART_REFERENCE_NONPASS"
        ),
        "authority_sha256": {
            **authorities,
            pilot_path.name: file_sha256(pilot_path),
            freeze_path.name: file_sha256(freeze_path),
            "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py": code_hash,
        },
        "seed_classification": seed_report["classification"],
        "heldout_was_unseen_when_tolerances_frozen": True,
        "window_kind": pilot["window_kind"],
        "heldout_windows": windows,
        "heldout_aggregate": aggregate,
        "frozen_thresholds": thresholds,
        "gates": gates,
        "all_heldout_gates_pass": passed,
        "checkpoint_eligible": passed,
        "new_AP1_M1_background_runs": 0,
        "diagnostic_envelope_integrations": 2,
        "trajectory_rows_persisted": 0,
        "background_candidate_assessed": False,
        "background_tolerances_frozen": False,
        "physical_response_kernel_started": False,
        "seed_released": False,
        "nonpass_stored_or_used_as_seed": False,
        "equations_changed": False,
        "physics_changed": False,
        "parameters_changed": False,
        "existing_gate_thresholds_changed": False,
        "validated_scope": (
            "already-seeded mode overlap transport through short tachyonic "
            "prefixes in log-amplitude/Riccati variables"
        ),
        "unvalidated_scope": (
            "adiabatic UV-tail-to-resolved entry of new modes, signed PV source "
            "cancellation on the moving split, coupled backreaction and long background"
        ),
        "wronskian_audit": (
            "independently integrated log_abs_Im_R invariant plus Decimal "
            "Wronskian; direct complex reconstruction retained as a cancellation "
            "diagnostic and excluded from the gate"
        ),
        "claim_boundary": (
            "short diagnostic moving-chart mode reference only; no physical "
            "background, UV-tail entry, kernel, observable, curve, fit or significance"
        ),
        "next_required": (
            "implement and validate the adiabatic UV-tail-to-resolved boundary "
            "handoff and signed high-precision PV cancellation before any long "
            "background pilot"
        ),
    }


def _write_json(report: dict[str, Any], output: Path | None) -> None:
    rendered = json.dumps(report, indent=2) + "\n"
    if output is None:
        print(rendered, end="")
    else:
        output.write_text(rendered, encoding="utf-8")


def main() -> None:
    import argparse

    parser = argparse.ArgumentParser()
    subparsers = parser.add_subparsers(dest="command", required=True)
    pilot_parser = subparsers.add_parser("pilot")
    pilot_parser.add_argument("apeiron_root", type=Path)
    pilot_parser.add_argument("code", type=Path)
    pilot_parser.add_argument("--output", type=Path)
    freeze_parser = subparsers.add_parser("freeze")
    freeze_parser.add_argument("pilot", type=Path)
    freeze_parser.add_argument("code", type=Path)
    freeze_parser.add_argument("--output", type=Path)
    heldout_parser = subparsers.add_parser("heldout")
    heldout_parser.add_argument("apeiron_root", type=Path)
    heldout_parser.add_argument("pilot", type=Path)
    heldout_parser.add_argument("freeze", type=Path)
    heldout_parser.add_argument("code", type=Path)
    heldout_parser.add_argument("--output", type=Path)
    args = parser.parse_args()
    if args.command == "pilot":
        report = build_pilot(args.apeiron_root, args.code)
    elif args.command == "freeze":
        report = build_freeze(args.pilot, args.code)
    else:
        report = build_heldout(
            args.apeiron_root, args.pilot, args.freeze, args.code
        )
    _write_json(report, args.output)


if __name__ == "__main__":
    main()
