"""Run a small synthetic transit injection-recovery campaign.

Run from the repository root:

    python research/tess-timing/scripts/run_injection.py

The script builds a deterministic grid of synthetic light curves with the
transit model, adds seeded Gaussian noise, and fits each one with the existing
estimator. It writes ``data/injection-run/injection.json`` and
``injection.svg``. All inputs are synthetic: no archive product, real target,
detection, or discovery is involved.

The grid has three injected mid-transit offsets, two noise levels, and an
antithetic noise pair at each cell. The measured runtime is pinned in
``runtime-ledger.json`` so repeated runs are byte-identical.
"""

from __future__ import annotations

import argparse
import hashlib
import json
import math
import sys
import time
from pathlib import Path

PACKAGE_ROOT = Path(__file__).resolve().parents[1]
REPO_ROOT = PACKAGE_ROOT.parents[1]
SRC = PACKAGE_ROOT / "src"
if str(SRC) not in sys.path:
    sys.path.insert(0, str(SRC))

# isort: off
import numpy as np
from tess_timing.fit import (
    FitError,
    duration_to_a_over_rs,
    fit_transit,
    software_versions,
)
from tess_timing.transit import transit_flux
# isort: on

DEFAULT_OUT_DIR = PACKAGE_ROOT / "data" / "injection-run"
DEFAULT_LEDGER = DEFAULT_OUT_DIR / "runtime-ledger.json"

SEED = 20261022
REFERENCE_T0_DAYS = 2450000.5
PERIOD_DAYS = 5.0
RADIUS_RATIO = 0.10
DURATION_DAYS = 0.04
IMPACT_PARAMETER = 0.30
U1 = 0.30
U2 = 0.25
BASELINE_OFFSET = 3.0e-4
CADENCE_DAYS = 2.0 / 1440.0
CADENCE_COUNT = 50
T0_OFFSETS_DAYS = (-0.005, 0.0, 0.005)
NOISE_SIGMAS = (6.0e-5, 1.0e-4)
SIGNS = (("plus", 1.0), ("minus", -1.0))
RUNTIME_BUDGET_SECONDS = 300.0
GATE_BIAS_OVER_SIGMA = 0.2

SECONDS_PER_DAY = 86400.0

SVG_NAMESPACE = (
    "http://www.w3.org/2000/svg"  # sec-ok: http-url - XML namespace identifier
)

INK = "#1b1b1b"
MUTED = "#67645c"
GRID = "#dedbd2"
AXIS = "#8c887e"
LOW_NOISE = "#245c4a"
HIGH_NOISE = "#a24b2a"
GATE = "#f0ece2"


def _display(path: Path) -> str:
    try:
        return path.resolve().relative_to(REPO_ROOT.resolve()).as_posix()
    except ValueError:
        return str(path)


def _escape(value: object) -> str:
    return str(value).replace("&", "&amp;").replace("<", "&lt;").replace(">", "&gt;")


def _text(
    x: float,
    y: float,
    value: object,
    size: float = 12.0,
    fill: str = INK,
    anchor: str = "start",
    weight: str = "normal",
) -> str:
    return (
        f'<text x="{x:.1f}" y="{y:.1f}" font-family="Georgia, serif" '
        f'font-size="{size:.1f}" fill="{fill}" text-anchor="{anchor}" '
        f'font-weight="{weight}">{_escape(value)}</text>'
    )


def _line(
    x1: float,
    y1: float,
    x2: float,
    y2: float,
    stroke: str,
    width: float = 1.0,
    dash: str | None = None,
) -> str:
    extra = f' stroke-dasharray="{dash}"' if dash else ""
    return (
        f'<line x1="{x1:.1f}" y1="{y1:.1f}" x2="{x2:.1f}" y2="{y2:.1f}" '
        f'stroke="{stroke}" stroke-width="{width:.2f}"{extra} />'
    )


def _rect(
    x: float,
    y: float,
    width: float,
    height: float,
    fill: str,
    stroke: str = "none",
    stroke_width: float = 0.0,
    opacity: float | None = None,
) -> str:
    extra = f' fill-opacity="{opacity:.2f}"' if opacity is not None else ""
    return (
        f'<rect x="{x:.1f}" y="{y:.1f}" width="{width:.1f}" '
        f'height="{height:.1f}" fill="{fill}" stroke="{stroke}" '
        f'stroke-width="{stroke_width:.2f}"{extra} />'
    )


def _circle(
    x: float,
    y: float,
    radius: float,
    fill: str,
    stroke: str = "none",
    stroke_width: float = 0.0,
) -> str:
    return (
        f'<circle cx="{x:.1f}" cy="{y:.1f}" r="{radius:.1f}" fill="{fill}" '
        f'stroke="{stroke}" stroke-width="{stroke_width:.2f}" />'
    )


def _scale(
    value: float, domain: tuple[float, float], span: tuple[float, float]
) -> float:
    low, high = domain
    return span[0] + (value - low) / (high - low) * (span[1] - span[0])


def _y_scale(
    value: float, domain: tuple[float, float], top: float, bottom: float
) -> float:
    return bottom - (value - domain[0]) / (domain[1] - domain[0]) * (bottom - top)


def _panel_axes(
    parts: list[str],
    x0: float,
    y0: float,
    width: float,
    height: float,
    x_domain: tuple[float, float],
    y_domain: tuple[float, float],
    x_ticks: list[float],
    y_ticks: list[float],
    x_label: str,
    y_label: str,
    x_format: str,
    y_format: str,
) -> tuple:
    def xmap(value: float) -> float:
        return _scale(value, x_domain, (x0, x0 + width))

    def ymap(value: float) -> float:
        return _y_scale(value, y_domain, y0, y0 + height)

    parts.append(_rect(x0, y0, width, height, "#ffffff", AXIS, 0.8))
    for value in y_ticks:
        y = ymap(value)
        parts.append(_line(x0, y, x0 + width, y, GRID, 0.6))
        parts.append(_text(x0 - 9, y + 4, y_format.format(value), 10.5, MUTED, "end"))
    for value in x_ticks:
        x = xmap(value)
        parts.append(_line(x, y0, x, y0 + height, GRID, 0.6))
        parts.append(
            _text(x, y0 + height + 17, x_format.format(value), 10.5, MUTED, "middle")
        )
    parts.append(
        _text(x0 + width / 2, y0 + height + 38, x_label, 11.5, MUTED, "middle")
    )
    parts.append(
        f'<text x="{x0 - 52:.1f}" y="{y0 + height / 2:.1f}" font-family="Georgia, serif" '
        f'font-size="11.5" fill="{MUTED}" text-anchor="middle" '
        f'transform="rotate(-90 {x0 - 52:.1f} {y0 + height / 2:.1f})">'
        f"{_escape(y_label)}</text>"
    )
    return xmap, ymap


def _grid() -> dict[str, object]:
    return {
        "injected_t0_offsets_days": list(T0_OFFSETS_DAYS),
        "noise_sigma": list(NOISE_SIGMAS),
        "signs": [name for name, _ in SIGNS],
        "injection_count": len(T0_OFFSETS_DAYS) * len(NOISE_SIGMAS) * len(SIGNS),
        "seed_rule": "base_seed + 100 * offset_index + noise_index",
    }


def _settings() -> dict[str, object]:
    return {
        "estimator": "tess_timing.fit.fit_transit",
        "model": "tess_timing.transit.transit_flux",
        "fixed_parameters": {
            "period_days": PERIOD_DAYS,
            "radius_ratio": RADIUS_RATIO,
            "duration_days": DURATION_DAYS,
            "impact_parameter": IMPACT_PARAMETER,
            "u1": U1,
            "u2": U2,
            "dilution": 0.0,
            "baseline_offset": BASELINE_OFFSET,
        },
        "time": {
            "reference_t0_days": REFERENCE_T0_DAYS,
            "cadence_days": CADENCE_DAYS,
            "cadence_seconds": CADENCE_DAYS * SECONDS_PER_DAY,
            "cadence_count": CADENCE_COUNT,
            "time_span_days": (CADENCE_COUNT - 1) * CADENCE_DAYS,
        },
        "noise": {
            "distribution": "Gaussian",
            "unit": "relative flux",
            "generator": "numpy.random.default_rng",
            "antithetic_pairs": True,
            "description": (
                "Each cell draws one standard-normal vector and adds it with a "
                "positive and a negative sign, so the pair has zero mean noise."
            ),
        },
        "fit": {
            "free_parameters": [
                "t0",
                "radius_ratio",
                "duration",
                "baseline_offset",
                "jitter",
            ],
            "initial_guesses": {
                "t0": None,
                "radius_ratio": RADIUS_RATIO,
                "duration_days": DURATION_DAYS,
            },
            "max_nfev": 4000,
        },
    }


def _campaign_fingerprint(
    seed: int, grid: dict[str, object], settings: dict[str, object]
) -> str:
    canonical = json.dumps(
        {"seed": seed, "grid": grid, "settings": settings},
        sort_keys=True,
        separators=(",", ":"),
    )
    return hashlib.sha256(canonical.encode("utf-8")).hexdigest()


def _run_campaign(
    seed: int,
    grid: dict[str, object],
    settings: dict[str, object],
) -> list[dict[str, object]]:
    time_values = (
        np.arange(CADENCE_COUNT, dtype=float) - CADENCE_COUNT / 2.0
    ) * CADENCE_DAYS + REFERENCE_T0_DAYS
    a_over_rs = duration_to_a_over_rs(
        DURATION_DAYS, PERIOD_DAYS, RADIUS_RATIO, IMPACT_PARAMETER
    )
    records: list[dict[str, object]] = []
    index = 0
    for offset_index, offset in enumerate(T0_OFFSETS_DAYS):
        injected_t0 = REFERENCE_T0_DAYS + offset
        model = (
            transit_flux(
                time_values,
                injected_t0,
                PERIOD_DAYS,
                RADIUS_RATIO,
                a_over_rs,
                IMPACT_PARAMETER,
                U1,
                U2,
            )
            + BASELINE_OFFSET
        )
        for noise_index, sigma in enumerate(NOISE_SIGMAS):
            cell_seed = seed + 100 * offset_index + noise_index
            rng = np.random.default_rng(cell_seed)
            standard_normal = rng.normal(0.0, 1.0, size=time_values.size)
            for sign_name, sign in SIGNS:
                index += 1
                flux = model + sigma * sign * standard_normal
                flux_err = np.full(time_values.size, sigma)
                result = fit_transit(
                    time_values,
                    flux,
                    flux_err,
                    period=PERIOD_DAYS,
                    u1=U1,
                    u2=U2,
                    impact_parameter=IMPACT_PARAMETER,
                    radius_ratio_guess=RADIUS_RATIO,
                    duration_guess=DURATION_DAYS,
                )
                error_days = result.t0 - injected_t0
                uncertainty_days = result.uncertainties["t0"]
                records.append(
                    {
                        "index": index,
                        "synthetic": True,
                        "grid_cell": {
                            "injected_t0_offset_days": offset,
                            "noise_sigma": sigma,
                            "noise_sign": sign_name,
                        },
                        "seed": cell_seed,
                        "injected_t0_days": injected_t0,
                        "recovered_t0_days": result.t0,
                        "reported_sigma_t0_days": uncertainty_days,
                        "error_days": error_days,
                        "error_seconds": error_days * SECONDS_PER_DAY,
                        "pull": error_days / uncertainty_days,
                        "within_0.2_sigma": abs(error_days)
                        <= GATE_BIAS_OVER_SIGMA * uncertainty_days,
                        "fit": {
                            "converged": result.converged,
                            "status": result.status,
                            "in_transit_cadences": result.in_transit_cadences,
                            "cadence_count": result.cadence_count,
                            "rms_ppm": result.rms_ppm,
                            "chi_square_per_dof": result.chi_square_per_dof,
                            "radius_ratio": result.radius_ratio,
                            "duration_days": result.duration,
                            "baseline_offset": result.baseline_offset,
                            "jitter": result.jitter,
                        },
                    }
                )
    return records


def _aggregates(records: list[dict[str, object]]) -> dict[str, object]:
    errors = np.array([float(record["error_days"]) for record in records])
    uncertainties = np.array(
        [float(record["reported_sigma_t0_days"]) for record in records]
    )
    pulls = np.array([float(record["pull"]) for record in records])
    bias_days = float(np.mean(errors))
    median_sigma_days = float(np.median(uncertainties))
    bias_over_sigma = abs(bias_days) / median_sigma_days
    scatter_days = float(np.std(errors, ddof=1))
    within_gate = np.abs(errors) <= GATE_BIAS_OVER_SIGMA * uncertainties
    within_one = np.abs(pulls) <= 1.0
    within_two = np.abs(pulls) <= 2.0
    return {
        "injection_count": len(records),
        "bias_days": bias_days,
        "bias_seconds": bias_days * SECONDS_PER_DAY,
        "median_reported_sigma_t0_days": median_sigma_days,
        "median_reported_sigma_t0_seconds": median_sigma_days * SECONDS_PER_DAY,
        "bias_over_median_sigma": bias_over_sigma,
        "scatter_days": scatter_days,
        "scatter_seconds": scatter_days * SECONDS_PER_DAY,
        "scatter_over_median_sigma": scatter_days / median_sigma_days,
        "pull_mean": float(np.mean(pulls)),
        "pull_rms": float(np.sqrt(np.mean(pulls**2))),
        "max_abs_pull": float(np.max(np.abs(pulls))),
        "fraction_within_0.2_sigma": float(np.mean(within_gate)),
        "fraction_within_1_sigma": float(np.mean(within_one)),
        "fraction_within_2_sigma": float(np.mean(within_two)),
        "count_within_0.2_sigma": int(np.count_nonzero(within_gate)),
        "all_fits_converged": all(
            bool(record["fit"]["converged"]) for record in records
        ),
    }


def _resolve_runtime(
    ledger_path: Path,
    fingerprint: str,
    measured_seconds: float,
    budget_seconds: float,
) -> float:
    if measured_seconds > budget_seconds:
        raise RuntimeError(
            f"campaign took {measured_seconds:.1f} s, above the "
            f"{budget_seconds:.0f} s budget; the grid was not written"
        )
    if ledger_path.exists():
        try:
            ledger = json.loads(ledger_path.read_text(encoding="utf-8"))
            recorded = float(ledger["runtime_seconds"])
            if (
                ledger.get("campaign_fingerprint") == fingerprint
                and math.isfinite(recorded)
                and 0.0 < recorded <= budget_seconds
            ):
                return recorded
        except (OSError, ValueError, KeyError, TypeError, json.JSONDecodeError):
            pass
    recorded = round(measured_seconds, 3)
    ledger_path.parent.mkdir(parents=True, exist_ok=True)
    ledger_path.write_text(
        json.dumps(
            {
                "campaign_fingerprint": fingerprint,
                "runtime_seconds": recorded,
                "budget_seconds": budget_seconds,
                "method": (
                    "Wall-clock time of the first run of this campaign. It is "
                    "reused for byte-identical output on later runs."
                ),
            },
            indent=2,
            sort_keys=True,
        )
        + "\n",
        encoding="utf-8",
    )
    return recorded


def build_svg(records: list[dict[str, object]], aggregates: dict[str, object]) -> str:
    width = 1080
    height = 720
    injected_offsets = [
        (float(record["injected_t0_days"]) - REFERENCE_T0_DAYS) * SECONDS_PER_DAY
        for record in records
    ]
    recovered_offsets = [
        (float(record["recovered_t0_days"]) - REFERENCE_T0_DAYS) * SECONDS_PER_DAY
        for record in records
    ]
    errors = [float(record["error_seconds"]) for record in records]
    uncertainties = [
        float(record["reported_sigma_t0_days"]) * SECONDS_PER_DAY for record in records
    ]
    limit = 500.0
    residual_limit = 3.0
    sigma_limit = 5.0
    median_sigma = float(aggregates["median_reported_sigma_t0_seconds"])
    gate = GATE_BIAS_OVER_SIGMA * median_sigma

    parts = [
        (
            f'<svg xmlns="{SVG_NAMESPACE}" width="{width}" height="{height}" '
            f'viewBox="0 0 {width} {height}" role="img" '
            'aria-labelledby="injection-title injection-desc">'
        ),
        '<title id="injection-title">Synthetic transit injection-recovery curve</title>',
        (
            '<desc id="injection-desc">Synthetic injection recovery, not a '
            "detection. Twelve light curves are generated from the transit model "
            "with seeded Gaussian noise and fitted with the existing estimator. "
            "The left panel compares recovered and injected mid-transit times; "
            "the right panel shows the timing residual against the reported "
            "one-sigma uncertainty. No real target is measured.</desc>"
        ),
        _rect(0, 0, width, height, "#ffffff"),
        _text(
            90,
            34,
            "Synthetic transit injection-recovery curve",
            22,
            INK,
            weight="bold",
        ),
        _text(
            90,
            58,
            "Synthetic inputs. This is an injection-recovery test, not a detection.",
            12.5,
            MUTED,
        ),
        _text(
            90,
            80,
            "Twelve light curves from tess_timing.transit with seeded Gaussian noise, "
            "fitted by tess_timing.fit.fit_transit.",
            11.5,
            MUTED,
        ),
    ]

    xmap, ymap = _panel_axes(
        parts,
        90.0,
        140.0,
        600.0,
        300.0,
        (-limit, limit),
        (-limit, limit),
        [-400.0, -200.0, 0.0, 200.0, 400.0],
        [-400.0, -200.0, 0.0, 200.0, 400.0],
        "injected mid-transit offset (s)",
        "recovered mid-transit offset (s)",
        "{:.0f}",
        "{:.0f}",
    )
    parts.append(
        _line(xmap(-limit), ymap(-limit), xmap(limit), ymap(limit), AXIS, 1.1, "6 4")
    )
    parts.append(_text(xmap(350), ymap(300), "1:1", 11.5, MUTED, "middle"))
    for record, x_value, y_value in zip(records, injected_offsets, recovered_offsets):
        color = (
            LOW_NOISE
            if float(record["grid_cell"]["noise_sigma"]) == NOISE_SIGMAS[0]
            else HIGH_NOISE
        )
        parts.append(_circle(xmap(x_value), ymap(y_value), 3.8, color))

    xmap2, ymap2 = _panel_axes(
        parts,
        770.0,
        140.0,
        280.0,
        300.0,
        (0.0, sigma_limit),
        (-residual_limit, residual_limit),
        [0.0, 1.0, 2.0, 3.0, 4.0, 5.0],
        [-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0],
        "reported 1-sigma t0 uncertainty (s)",
        "recovered minus injected t0 (s)",
        "{:.0f}",
        "{:.0f}",
    )
    parts.append(
        _rect(
            xmap2(0.0),
            ymap2(gate),
            xmap2(sigma_limit) - xmap2(0.0),
            ymap2(-gate) - ymap2(gate),
            GATE,
            "none",
            0.0,
            0.85,
        )
    )
    parts.append(
        _line(xmap2(0.0), ymap2(0.0), xmap2(sigma_limit), ymap2(0.0), AXIS, 1.0)
    )
    parts.append(
        _line(
            xmap2(0.0),
            ymap2(gate),
            xmap2(sigma_limit),
            ymap2(gate),
            HIGH_NOISE,
            0.8,
            "4 3",
        )
    )
    parts.append(
        _line(
            xmap2(0.0),
            ymap2(-gate),
            xmap2(sigma_limit),
            ymap2(-gate),
            HIGH_NOISE,
            0.8,
            "4 3",
        )
    )
    for record, x_value, y_value, uncertainty in zip(
        records, uncertainties, errors, uncertainties
    ):
        color = (
            LOW_NOISE
            if float(record["grid_cell"]["noise_sigma"]) == NOISE_SIGMAS[0]
            else HIGH_NOISE
        )
        parts.append(
            _line(
                xmap2(x_value),
                ymap2(y_value - uncertainty),
                xmap2(x_value),
                ymap2(y_value + uncertainty),
                color,
                0.8,
            )
        )
        parts.append(_circle(xmap2(x_value), ymap2(y_value), 3.4, color))
    parts.append(_text(782, 160, "Timing residual", 13, INK, weight="bold"))
    parts.append(_circle(782, 184, 3.4, LOW_NOISE))
    parts.append(_text(792, 188, "sigma = 6.0e-5", 10.5, MUTED))
    parts.append(_circle(782, 202, 3.4, HIGH_NOISE))
    parts.append(_text(792, 206, "sigma = 1.0e-4", 10.5, MUTED))
    parts.append(_text(782, 226, "band: 0.2 x median sigma", 10.5, MUTED))

    parts.append(_line(90, 488, 1050, 488, AXIS, 0.8))
    parts.append(
        _text(90, 514, "Aggregate over the synthetic grid", 13, INK, weight="bold")
    )
    count = int(aggregates["injection_count"])
    within_gate = int(aggregates["count_within_0.2_sigma"])
    gate_pass = bool(aggregates["bias_over_median_sigma"] < GATE_BIAS_OVER_SIGMA)
    metric_lines = [
        (
            f"{count} injections: 3 injected offsets x 2 noise levels x 2 "
            "antithetic signs"
        ),
        (
            f"timing bias {float(aggregates['bias_seconds']):+.2f} s "
            f"({float(aggregates['bias_days']):+.2e} d)"
        ),
        (
            f"median reported sigma {float(aggregates['median_reported_sigma_t0_seconds']):.2f} s; "
            f"|bias| / sigma = {float(aggregates['bias_over_median_sigma']):.3f} "
            f"(gate < {GATE_BIAS_OVER_SIGMA:.2f}: {'PASS' if gate_pass else 'FAIL'})"
        ),
        (
            f"scatter {float(aggregates['scatter_seconds']):.2f} s; "
            f"pull RMS {float(aggregates['pull_rms']):.2f}; "
            f"all fits converged: {'yes' if aggregates['all_fits_converged'] else 'no'}"
        ),
        (
            f"within 0.2 sigma: {within_gate}/{count} "
            f"({float(aggregates['fraction_within_0.2_sigma']):.3f}); "
            f"within 1 sigma: {float(aggregates['fraction_within_1_sigma']):.3f}; "
            f"within 2 sigma: {float(aggregates['fraction_within_2_sigma']):.3f}"
        ),
    ]
    y = 538.0
    for line in metric_lines:
        parts.append(_text(90, y, line, 11.5, MUTED))
        y += 20.0

    parts.append(_line(90, 660, 1050, 660, AXIS, 0.8))
    parts.append(
        _text(
            90,
            680,
            "Inputs are synthetic light curves generated by this script. "
            "No archive data, real target, detection, or discovery is involved.",
            10.2,
            MUTED,
        )
    )
    parts.append(
        _text(
            90,
            696,
            "Generated by research/tess-timing/scripts/run_injection.py. "
            "Values are the same records written to injection.json.",
            10.2,
            MUTED,
        )
    )
    parts.append("</svg>")
    return "\n".join(parts) + "\n"


def _parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description=(
            "Run a small synthetic transit injection-recovery campaign with "
            "seeded noise and write a JSON record and an SVG recovery curve. "
            "No network request is made."
        )
    )
    parser.add_argument(
        "--out-dir",
        type=Path,
        default=DEFAULT_OUT_DIR,
        help=f"output directory (default: {_display(DEFAULT_OUT_DIR)})",
    )
    parser.add_argument(
        "--runtime-ledger",
        type=Path,
        default=DEFAULT_LEDGER,
        help=f"runtime ledger path (default: {_display(DEFAULT_LEDGER)})",
    )
    parser.add_argument("--seed", type=int, default=SEED)
    parser.add_argument(
        "--budget-seconds",
        type=float,
        default=RUNTIME_BUDGET_SECONDS,
        help=f"wall-clock budget in seconds (default: {RUNTIME_BUDGET_SECONDS:.0f})",
    )
    return parser


def main(argv: list[str] | None = None) -> int:
    arguments = _parser().parse_args(argv)
    grid = _grid()
    settings = _settings()
    fingerprint = _campaign_fingerprint(arguments.seed, grid, settings)

    start = time.perf_counter()
    try:
        records = _run_campaign(arguments.seed, grid, settings)
        aggregates = _aggregates(records)
        runtime_seconds = _resolve_runtime(
            arguments.runtime_ledger,
            fingerprint,
            time.perf_counter() - start,
            arguments.budget_seconds,
        )
        svg = build_svg(records, aggregates)
    except (FitError, RuntimeError, ValueError) as error:
        print(f"error: {error}", file=sys.stderr)
        return 1

    payload = {
        "task": "jdp-tess-injection",
        "schema_version": 1,
        "synthetic": True,
        "label": "SYNTHETIC INJECTION RECOVERY",
        "label_text": (
            "Synthetic injection-recovery campaign. This is not a detection. "
            "No observation, real target, or discovery is involved."
        ),
        "seed": arguments.seed,
        "campaign_fingerprint": fingerprint,
        "grid": grid,
        "settings": settings,
        "runtime_seconds": runtime_seconds,
        "runtime": {
            "seconds": runtime_seconds,
            "budget_seconds": arguments.budget_seconds,
            "within_budget": runtime_seconds <= arguments.budget_seconds,
            "method": (
                "Wall-clock time from the start of the run through the injection "
                "fits and figure generation. The first measured value for this "
                "campaign is pinned in runtime-ledger.json so repeated runs are "
                "byte-identical."
            ),
        },
        "software": software_versions(),
        "gate": {
            "source": "docs/SIMULATION-ROADMAP.md section 5.4",
            "quantity": "absolute timing bias divided by median reported t0 uncertainty",
            "threshold": GATE_BIAS_OVER_SIGMA,
            "value": aggregates["bias_over_median_sigma"],
            "passes": aggregates["bias_over_median_sigma"] < GATE_BIAS_OVER_SIGMA,
        },
        "aggregates": aggregates,
        "injections": records,
        "limits": [
            "Every light curve is synthetic and generated in this script.",
            "No archive product, real target, detection, or discovery is used.",
            "Period, limb darkening, impact parameter, and dilution are fixed.",
            "The fit is a demonstration on seeded noise, not a real timing measurement.",
        ],
    }

    try:
        arguments.out_dir.mkdir(parents=True, exist_ok=True)
        json_path = arguments.out_dir / "injection.json"
        svg_path = arguments.out_dir / "injection.svg"
        json_path.write_text(
            json.dumps(payload, indent=2, ensure_ascii=False, allow_nan=False) + "\n",
            encoding="utf-8",
        )
        svg_path.write_text(svg, encoding="utf-8")
    except OSError as error:
        print(f"error: {error}", file=sys.stderr)
        return 1

    bias_over_sigma = float(aggregates["bias_over_median_sigma"])
    print(
        f"status=OK synthetic=yes injections={len(records)} "
        f"bias={float(aggregates['bias_seconds']):+.3f}s "
        f"median_sigma={float(aggregates['median_reported_sigma_t0_seconds']):.3f}s "
        f"bias_over_sigma={bias_over_sigma:.3f} "
        f"gate={'PASS' if bias_over_sigma < GATE_BIAS_OVER_SIGMA else 'FAIL'} "
        f"runtime={runtime_seconds:.3f}s "
        f"wrote={_display(json_path)},{_display(svg_path)}"
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())