"""Build the folded light curve, epoch table, and figure.

Run from the repository root:

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

The command is offline. It reads the committed synthetic light curve, the
committed trapezoid fit record, and the committed target docket. It writes:

- ``research/tess-timing/data/epoch-fold-run/epoch_measurements.json``
- ``research/tess-timing/data/epoch-fold-run/folded_lightcurve.csv``
- ``docs/figures/epoch-fold.svg``
- ``docs/EPOCH-SHEET.md``

Every number describes the labelled synthetic fixture. No target is selected,
no observation is claimed, and no detection is made.
"""

from __future__ import annotations

import argparse
import math
import sys
from pathlib import Path

import numpy as np

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
from tess_timing.epoch_fold import (
    DEFAULT_CSV,
    DEFAULT_FIGURE,
    DEFAULT_FIT_JSON,
    DEFAULT_JSON,
    DEFAULT_LIGHTCURVE,
    DEFAULT_MARKDOWN,
    DEFAULT_TARGET_LIST,
    SECONDS_PER_DAY,
    build_epoch_sheet,
    read_lightcurve,
    relative_path,
    write_epoch_sheet,
)
from tess_timing.toy_transit import trapezoid_flux
# isort: on

SVG_NAMESPACE = (
    "http://www.w3.org/2000/svg"  # sec-ok: http-url - XML namespace identifier
)
INK = "#1b1b1b"
MUTED = "#67645c"
GRID = "#dedbd2"
AXIS = "#8c887e"
DATA = "#3d5a80"
FIT = "#a24b2a"
ZERO = "#245c4a"
PANEL = "#fffdf8"
BACKGROUND = "#faf7ef"
PLOT_CADENCES = 2500


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,
    opacity: float | None = None,
) -> str:
    extra = f' stroke-dasharray="{dash}"' if dash else ""
    if opacity is not None:
        extra += f' stroke-opacity="{opacity:.2f}"'
    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,
) -> str:
    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}" />'
    )


def _circle(
    x: float,
    y: float,
    radius: float,
    fill: str,
    opacity: float | None = None,
) -> str:
    extra = f' fill-opacity="{opacity:.2f}"' if opacity is not None else ""
    return f'<circle cx="{x:.1f}" cy="{y:.1f}" r="{radius:.1f}" fill="{fill}"{extra} />'


def _polyline(
    points: list[tuple[float, float]],
    stroke: str,
    width: float,
    dash: str | None = None,
) -> str:
    encoded = " ".join(f"{x:.1f},{y:.1f}" for x, y in points)
    extra = f' stroke-dasharray="{dash}"' if dash else ""
    return (
        f'<polyline points="{encoded}" fill="none" stroke="{stroke}" '
        f'stroke-width="{width:.2f}" stroke-linejoin="round" '
        f'stroke-linecap="round"{extra} />'
    )


def _x_scale(
    value: float, domain: tuple[float, float], left: float, right: float
) -> float:
    low, high = domain
    return left + (value - low) / (high - low) * (right - left)


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 _nice_ticks(low: float, high: float, count: int = 5) -> list[float]:
    span = high - low
    if span <= 0.0:
        return [low]
    raw_step = span / max(count - 1, 1)
    magnitude = 10.0 ** math.floor(math.log10(raw_step))
    residual = raw_step / magnitude
    if residual <= 1.0:
        step = 1.0
    elif residual <= 2.0:
        step = 2.0
    elif residual <= 5.0:
        step = 5.0
    else:
        step = 10.0
    step *= magnitude
    start = math.ceil(low / step) * step
    ticks: list[float] = []
    value = start
    while value <= high + 1.0e-9:
        ticks.append(value)
        value += step
    return ticks


def _panel(
    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,
) -> None:
    parts.append(_rect(x0, y0, width, height, PANEL, AXIS, 0.8))
    for tick in y_ticks:
        y = _y_scale(tick, y_domain, y0, y0 + height)
        parts.append(_line(x0, y, x0 + width, y, GRID, 0.8))
        parts.append(_text(x0 - 10, y + 4, y_format.format(tick), 11.0, MUTED, "end"))
    for tick in x_ticks:
        x = _x_scale(tick, x_domain, x0, x0 + width)
        parts.append(_line(x, y0, x, y0 + height, GRID, 0.8))
        parts.append(
            _text(x, y0 + height + 20, x_format.format(tick), 11.0, MUTED, "middle")
        )
    parts.append(
        _text(
            x0 + width / 2,
            y0 + height + 45,
            x_label,
            12.0,
            INK,
            "middle",
        )
    )
    parts.append(
        f'<text x="{x0 - 62:.1f}" y="{y0 + height / 2:.1f}" '
        f'font-family="Georgia, serif" font-size="12.0" fill="{INK}" '
        f'text-anchor="middle" transform="rotate(-90 {x0 - 62:.1f} '
        f'{y0 + height / 2:.1f})">{_escape(y_label)}</text>'
    )


def _render_svg(sheet: dict, lightcurve) -> str:
    reference = sheet["reference"]
    period = reference["period_days"]
    bins = [
        row
        for row in sheet["fold"]["bins"]
        if row["count"] > 0 and row["flux"] is not None
    ]
    epochs = [row for row in sheet["epochs"] if row["status"] == "measured"]

    width = 1080.0
    height = 880.0
    parts: list[str] = []
    parts.append(
        f'<svg xmlns="{SVG_NAMESPACE}" viewBox="0 0 {width:.0f} {height:.0f}" '
        f'width="{width:.0f}" height="{height:.0f}" role="img" '
        f'aria-labelledby="title desc">'
    )
    parts.append('<title id="title">Phase-folded synthetic light curve</title>')
    parts.append(
        '<desc id="desc">Top panel: relative flux folded on the reference '
        "ephemeris with binned means and a trapezoid model. Bottom panel: "
        "per-transit timing residuals from a linear consistency fit. The data "
        "are a synthetic fixture, not a TESS observation, and no detection is "
        "claimed.</desc>"
    )
    parts.append(_rect(0, 0, width, height, BACKGROUND))
    parts.append(
        _text(60, 55, "Phase-folded synthetic light curve", 22.0, INK, "start", "bold")
    )
    parts.append(
        _text(
            60,
            80,
            "Labelled synthetic fixture, 2-minute cadence, 27 d. "
            "Not a TESS observation. No detection is claimed.",
            12.5,
            MUTED,
        )
    )

    top_x0, top_y0, top_w, top_h = 105.0, 110.0, 900.0, 420.0
    phase_domain = (-0.5 * period, 0.5 * period)
    flux_values = np.asarray(lightcurve.flux, dtype=float)
    flux_min = float(np.min(flux_values))
    flux_max = float(np.max(flux_values))
    margin = max(0.0004, 0.08 * (flux_max - flux_min))
    flux_domain = (flux_min - margin, flux_max + margin)
    phase_ticks = _nice_ticks(phase_domain[0], phase_domain[1], 8)
    flux_ticks = _nice_ticks(flux_domain[0], flux_domain[1], 5)
    _panel(
        parts,
        top_x0,
        top_y0,
        top_w,
        top_h,
        phase_domain,
        flux_domain,
        phase_ticks,
        flux_ticks,
        "Phase (days from reference mid-transit)",
        "Relative flux",
        "{:.1f}",
        "{:.3f}",
    )

    def px(value: float) -> float:
        return _x_scale(value, phase_domain, top_x0, top_x0 + top_w)

    def py(value: float) -> float:
        return _y_scale(value, flux_domain, top_y0, top_y0 + top_h)

    stride = max(1, lightcurve.time.size // PLOT_CADENCES)
    indices = range(0, lightcurve.time.size, stride)
    for index in indices:
        phase = (lightcurve.time[index] - reference["t0_btjd"]) % period
        if phase > 0.5 * period:
            phase -= period
        parts.append(
            _circle(
                px(float(phase)), py(float(lightcurve.flux[index])), 1.5, DATA, 0.30
            )
        )

    for row in bins:
        x = px(row["phase_days"])
        y = py(row["flux"])
        low = py(row["flux"] - row["uncertainty"])
        high = py(row["flux"] + row["uncertainty"])
        parts.append(_line(x, low, x, high, DATA, 1.0, opacity=0.75))
        parts.append(_circle(x, y, 2.2, DATA))

    model_phase = np.linspace(-0.5 * period, 0.5 * period, 900)
    model_flux = trapezoid_flux(
        reference["t0_btjd"] + model_phase,
        reference["t0_btjd"],
        period,
        reference["depth"],
        reference["duration_days"],
        reference["baseline_offset"],
        reference["ingress_fraction"],
    )
    parts.append(
        _polyline(
            [
                (px(float(phase)), py(float(flux)))
                for phase, flux in zip(model_phase, model_flux)
            ],
            FIT,
            2.2,
        )
    )

    legend_y = top_y0 + 22
    parts.append(_circle(top_x0 + 20, legend_y, 2.5, DATA, 0.55))
    parts.append(
        _text(
            top_x0 + 32,
            legend_y + 4,
            f"individual cadences (1 in {stride})",
            11.0,
            MUTED,
        )
    )
    parts.append(_circle(top_x0 + 240, legend_y, 3.0, DATA))
    parts.append(_text(top_x0 + 252, legend_y + 4, "binned weighted mean", 11.0, MUTED))
    parts.append(_line(top_x0 + 400, legend_y, top_x0 + 430, legend_y, FIT, 2.2))
    parts.append(_text(top_x0 + 438, legend_y + 4, "reference trapezoid", 11.0, MUTED))

    bottom_x0, bottom_y0, bottom_w, bottom_h = 105.0, 610.0, 900.0, 200.0
    parts.append(
        _text(
            60,
            600,
            "Per-transit timing residuals from the linear consistency fit",
            14.0,
            INK,
            "start",
            "bold",
        )
    )
    residuals_seconds = [
        float(row["residual_vs_line_days"]) * SECONDS_PER_DAY for row in epochs
    ]
    sigma_seconds = [float(row["uncertainty_days"]) * SECONDS_PER_DAY for row in epochs]
    if residuals_seconds:
        largest = max(
            max(
                abs(value) + sigma
                for value, sigma in zip(residuals_seconds, sigma_seconds)
            ),
            3.0,
        )
        limit = math.ceil(largest / 2.0) * 2.0
        residual_domain = (-limit, limit)
    else:
        residual_domain = (-5.0, 5.0)
    cycles = [row["cycle"] for row in epochs]
    if cycles:
        cycle_domain = (min(cycles) - 0.6, max(cycles) + 0.6)
    else:
        cycle_domain = (-1.0, 1.0)
    cycle_ticks = _nice_ticks(cycle_domain[0], cycle_domain[1], 7)
    residual_ticks = _nice_ticks(residual_domain[0], residual_domain[1], 5)
    _panel(
        parts,
        bottom_x0,
        bottom_y0,
        bottom_w,
        bottom_h,
        cycle_domain,
        residual_domain,
        cycle_ticks,
        residual_ticks,
        "Cycle number relative to the reference t0",
        "Residual (s)",
        "{:.0f}",
        "{:+.0f}",
    )

    def bx(value: float) -> float:
        return _x_scale(value, cycle_domain, bottom_x0, bottom_x0 + bottom_w)

    def by(value: float) -> float:
        return _y_scale(value, residual_domain, bottom_y0, bottom_y0 + bottom_h)

    zero_y = by(0.0)
    parts.append(
        _line(bottom_x0, zero_y, bottom_x0 + bottom_w, zero_y, ZERO, 1.2, dash="5 4")
    )
    for cycle, residual, sigma in zip(cycles, residuals_seconds, sigma_seconds):
        x = bx(float(cycle))
        y = by(residual)
        parts.append(
            _line(
                x,
                by(residual - sigma),
                x,
                by(residual + sigma),
                DATA,
                1.1,
                opacity=0.85,
            )
        )
        parts.append(_circle(x, y, 3.2, DATA))

    parts.append(
        _text(
            60,
            845,
            "Top: every cadence folded on the committed reference ephemeris. "
            "Bottom: measured epoch minus the linear consistency fit; bars are 1 sigma. "
            "Synthetic fixture; not a detection.",
            11.0,
            MUTED,
        )
    )
    parts.append("</svg>")
    return "\n".join(parts) + "\n"


def _parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description=(
            "Fold the committed synthetic light curve, measure one mid-transit "
            "time per event, and write a JSON record, a folded CSV, an SVG "
            "figure, and docs/EPOCH-SHEET.md. The command is offline and makes "
            "no detection claim."
        )
    )
    parser.add_argument("--lightcurve", type=Path, default=DEFAULT_LIGHTCURVE)
    parser.add_argument("--fit-json", type=Path, default=DEFAULT_FIT_JSON)
    parser.add_argument("--target-list", type=Path, default=DEFAULT_TARGET_LIST)
    parser.add_argument("--json-out", type=Path, default=DEFAULT_JSON)
    parser.add_argument("--csv-out", type=Path, default=DEFAULT_CSV)
    parser.add_argument("--figure-out", type=Path, default=DEFAULT_FIGURE)
    parser.add_argument("--markdown-out", type=Path, default=DEFAULT_MARKDOWN)
    parser.add_argument("--bin-count", type=int, default=140)
    parser.add_argument("--window-days", type=float, default=None)
    return parser


def main(argv: list[str] | None = None) -> int:
    args = _parser().parse_args(argv)
    try:
        sheet = build_epoch_sheet(
            args.lightcurve,
            args.fit_json,
            args.target_list,
            bin_count=args.bin_count,
            window_days=args.window_days,
        )
        json_path, csv_path, markdown_path = write_epoch_sheet(
            sheet,
            json_path=args.json_out,
            csv_path=args.csv_out,
            markdown_path=args.markdown_out,
            lightcurve_path=args.lightcurve,
        )
        figure_path = args.figure_out
        figure_path.parent.mkdir(parents=True, exist_ok=True)
        figure_path.write_text(
            _render_svg(sheet, read_lightcurve(args.lightcurve)),
            encoding="utf-8",
            newline="\n",
        )
    except (OSError, ValueError, KeyError, TypeError) as error:
        print(f"status=FAILED error={error}", file=sys.stderr)
        return 1

    linear = sheet["linear_fit"]
    print("status=OK synthetic=yes not_a_detection=yes")
    print(f"lightcurve={relative_path(args.lightcurve)}")
    print(f"cadences={sheet['inputs']['lightcurve']['cadences']}")
    print(f"events_measured={sheet['counts']['events_measured']}")
    print(f"events_refused={sheet['counts']['events_refused']}")
    print(
        f"reference_t0={sheet['reference']['t0_btjd']:.6f} "
        f"reference_period={sheet['reference']['period_days']:.9f}"
    )
    if linear.get("status") == "OK":
        print(
            f"line_t0={linear['t0_btjd']:.6f} +/- "
            f"{linear['t0_uncertainty_days']:.3e} "
            f"line_period={linear['period_days']:.9f} +/- "
            f"{linear['period_uncertainty_days']:.3e} "
            f"chi2_dof={linear['chi_square_per_dof']:.3f}"
        )
    print(f"wrote={relative_path(json_path)}")
    print(f"wrote={relative_path(csv_path)}")
    print(f"wrote={relative_path(figure_path)}")
    print(f"wrote={relative_path(markdown_path)}")
    return 0


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