#!/usr/bin/env python3
"""Parse live firmware PPS telemetry and replay it through the Python loop.

The firmware status record is an observation of the already-running board.  It
contains the WLS clock-bias path, Doppler datum predictor, quality selector and
RTL PI diagnostics, but it intentionally does not contain a trusted external
PPS edge.  This module therefore reports two separate results:

* direct statistics of the observed firmware/RTL state; and
* a telemetry-conditioned replay of the bit-accurate fixed-point loop model,
  initialized from the captured locked PI/filter state.

When ``GPS_PPS_STATUS_LOG_INTERVAL`` is greater than one, missing solve
epochs are linearly interpolated.  The replay reconstructs a self-consistent
reference trajectory from the logged filtered bias, PI residual and rate word.
That preserves measured low-frequency drift but is a lower bound for
unobserved high-frequency measurement noise.  It is not an independent trusted
PPS measurement.
"""

from __future__ import annotations

import argparse
import csv
from dataclasses import asdict, dataclass, replace
import math
from pathlib import Path
import re
import statistics
from typing import Iterable

from pps_python.adaptive_pps import (
    FixedPointAdaptivePps,
    LoopProfile,
    Observation,
    PROFILES,
    Profile,
    Q32,
    QualityMetrics,
    SPEED_OF_LIGHT_MPS,
    TRACKING_HZ,
)


UPDATE_PERIOD_S = 0.2
NAV_BIT_PERIOD_S = 0.020


PPS_RE = re.compile(
    r"PPS profile=(?P<profile>\w+) mode=(?P<mode>\w+) "
    r"locked=(?P<locked>\d+) tdop=(?P<tdop>[-+\d.]+) "
    r"gdop=(?P<gdop>[-+\d.]+) "
    r"clock_sigma=(?P<clock_sigma_m>[-+\d.]+) "
    r"cn0=(?P<cn0_mean_dbhz>[-+\d.]+)/(?P<cn0_min_dbhz>[-+\d.]+) "
    r"innovation=(?P<innovation_m>[-+\d.]+)m "
    r"raw_bias=(?P<raw_bias_m>[-+\d.]+)m "
    r"filtered_bias=(?P<filtered_bias_m>[-+\d.]+)m "
    r"clock_drift=(?P<clock_drift_mps>[-+\d.]+)mps "
    r"bias_drift=(?P<bias_drift_mps>[-+\d.]+)mps "
    r"predictor=(?P<predictor>\d+) "
    r"eskf_alpha=(?P<eskf_alpha>[-+\d.]+) "
    r"eskf_v=(?P<eskf_speed_mps>[-+\d.]+)/"
    r"(?P<eskf_velocity_rms_mps>[-+\d.]+)/0x(?P<eskf_fde_mask>[0-9a-fA-F]+) "
    r"pps_v=(?P<pps_speed_mps>[-+\d.]+)/"
    r"(?P<pps_velocity_rms_mps>[-+\d.]+)/0x(?P<pps_fde_mask>[0-9a-fA-F]+) "
    r"ref_svid=(?P<reference_svid>\d+) "
    r"capture_corr=(?P<capture_correction_ns>[-+\d.]+)ns "
    r"hw_phase=(?P<hardware_phase_ns>[-+\d.]+)ns "
    r"rate_cmd=(?P<rate_command_ppm>[-+\d.]+)ppm "
    r"integ=(?P<integrator_ppm>[-+\d.]+)ppm "
    r"pulse=(?P<pulse_count>\d+) seq=(?P<sequence>\d+) "
    r"status=0x(?P<status>[0-9a-fA-F]+)"
)

LOCATION_RE = re.compile(
    r"source=gps .*?gdop=(?P<gdop>[-+\d.]+) "
    r"pdop=(?P<pdop>[-+\d.]+) hdop=(?P<hdop>[-+\d.]+) "
    r"vdop=(?P<vdop>[-+\d.]+) tdop=(?P<tdop>[-+\d.]+) "
    r"tracked=(?P<tracked>\d+) used=(?P<used>\d+) "
    r"velocity_valid=(?P<velocity_valid>\d+) "
    r"velocity_fde=(?P<velocity_fde>\d+) "
    r"speed=(?P<speed_mps>[-+\d.]+) "
    r"velocity_rms=(?P<velocity_rms_mps>[-+\d.]+)"
)

TRACKED_RE = re.compile(r"tracked satellites=(?P<tracked>\d+)")


@dataclass(frozen=True)
class GpsContext:
    host_monotonic_s: float
    tracked_satellites: int
    used_satellites: int
    hdop: float
    vdop: float
    velocity_valid: bool
    velocity_fde_applied: bool
    velocity_residual_rms_mps: float


@dataclass(frozen=True)
class FirmwarePpsSample:
    host_monotonic_s: float
    host_unix_s: float
    time_s: float
    profile: str
    mode: str
    locked: bool
    tdop: float
    gdop: float
    hdop: float
    vdop: float
    clock_sigma_m: float
    cn0_mean_dbhz: float
    cn0_min_dbhz: float
    tracked_satellites: int
    used_satellites: int
    innovation_m: float
    raw_bias_m: float
    filtered_bias_m: float
    clock_drift_mps: float
    bias_drift_mps: float
    predictor_enabled: bool
    eskf_alpha: float
    eskf_speed_mps: float
    eskf_velocity_rms_mps: float
    eskf_fde_mask: int
    pps_speed_mps: float
    pps_velocity_rms_mps: float
    pps_fde_mask: int
    reference_svid: int
    capture_correction_ns: float
    hardware_phase_ns: float
    rate_command_ppm: float
    integrator_ppm: float
    pulse_count: int
    sequence: int
    status: int


@dataclass(frozen=True)
class ObservedMetrics:
    samples: int
    duration_s: float
    status_interval_median_s: float
    sequence_step_median: float
    locked_fraction: float
    profile_transitions: int
    reference_changes: int
    hardware_phase_mean_ns: float
    hardware_phase_centered_rms_ns: float
    hardware_phase_centered_p95_ns: float
    hardware_phase_centered_max_ns: float
    hardware_phase_peak_to_peak_ns: float
    hardware_phase_10s_wander_peak_to_peak_ns: float
    sampled_phase_step_rms_ns: float
    sampled_phase_step_p95_ns: float
    sampled_phase_step_max_ns: float
    hardware_phase_slope_ns_per_s: float
    rate_command_median_ppm: float
    rate_command_span_ppm: float
    integrator_median_ppm: float
    inferred_oscillator_ppm: float
    inferred_active_phase_rms_ns: float
    inferred_active_phase_peak_to_peak_ns: float
    phase_rate_correlation: float
    phase_innovation_correlation: float
    pi_rate_reproduction_rms_ppb: float
    pi_rate_reproduction_p95_ppb: float
    pi_rate_reproduction_max_ppb: float
    innovation_rms_m: float
    innovation_p95_m: float
    raw_bias_slope_mps: float
    filtered_bias_slope_mps: float
    tdop_median: float
    gdop_median: float
    clock_sigma_median_m: float
    cn0_mean_median_dbhz: float
    cn0_min_median_dbhz: float
    tracked_median: float
    used_median: float
    pps_velocity_rms_median_mps: float


@dataclass(frozen=True)
class ReplayPoint:
    time_s: float
    source_sequence: int
    interpolated: bool
    raw_bias_m: float
    observed_filtered_bias_m: float
    modeled_filtered_bias_m: float
    reconstructed_active_phase_ns: float
    modeled_active_phase_ns: float
    observed_hardware_phase_ns: float
    modeled_loop_error_ns: float
    observed_rate_command_ppm: float
    modeled_rate_command_ppm: float


@dataclass(frozen=True)
class ReplayMetrics:
    tuning: str
    expanded_updates: int
    interpolation_factor: float
    first_lock_s: float
    locked_fraction_after_warmup: float
    unlocks: int
    inferred_oscillator_ppm: float
    modeled_edge_centered_rms_ns: float
    modeled_edge_centered_p95_ns: float
    modeled_edge_peak_to_peak_ns: float
    modeled_period_rms_ns: float
    modeled_period_p95_ns: float
    modeled_period_max_ns: float
    modeled_loop_residual_centered_rms_ns: float
    modeled_loop_residual_centered_p95_ns: float
    modeled_loop_residual_peak_to_peak_ns: float
    logged_residual_reproduction_rms_ns: float
    logged_residual_reproduction_p95_ns: float
    filtered_bias_reproduction_rms_m: float
    filtered_bias_reproduction_p95_m: float


def _percentile(values: Iterable[float], fraction: float) -> float:
    ordered = sorted(values)
    if not ordered:
        return math.nan
    position = (len(ordered) - 1) * fraction
    lower = math.floor(position)
    upper = math.ceil(position)
    if lower == upper:
        return ordered[lower]
    weight = position - lower
    return ordered[lower] * (1.0 - weight) + ordered[upper] * weight


def _rms(values: Iterable[float]) -> float:
    items = list(values)
    return (
        math.sqrt(statistics.mean(value * value for value in items))
        if items
        else math.nan
    )


def _linear_fit(x: list[float], y: list[float]) -> tuple[float, float]:
    if len(x) != len(y) or not x:
        raise ValueError("linear fit requires equal non-empty inputs")
    if len(x) == 1:
        return y[0], 0.0
    mean_x = statistics.mean(x)
    mean_y = statistics.mean(y)
    denominator = sum((value - mean_x) ** 2 for value in x)
    if denominator == 0.0:
        return mean_y, 0.0
    slope = sum(
        (x_value - mean_x) * (y_value - mean_y)
        for x_value, y_value in zip(x, y)
    ) / denominator
    return mean_y - slope * mean_x, slope


def _correlation(first: list[float], second: list[float]) -> float:
    if len(first) != len(second) or len(first) < 2:
        return math.nan
    first_mean = statistics.mean(first)
    second_mean = statistics.mean(second)
    first_centered = [value - first_mean for value in first]
    second_centered = [value - second_mean for value in second]
    denominator = math.sqrt(
        sum(value * value for value in first_centered)
        * sum(value * value for value in second_centered)
    )
    if denominator == 0.0:
        return math.nan
    return sum(
        first_value * second_value
        for first_value, second_value in zip(first_centered, second_centered)
    ) / denominator


def _profile(sample: FirmwarePpsSample) -> Profile:
    try:
        return Profile[sample.profile.upper()]
    except KeyError:
        return Profile.HOLDOVER


def _pi_rate_prediction_ppm(
    sample: FirmwarePpsSample, update_period_s: float
) -> float:
    """Recreate RTL ``integrator + Kp * phase`` from diagnostic fields."""
    profile = _profile(sample)
    if profile == Profile.HOLDOVER or sample.mode != "tracking":
        return math.nan
    kp_q46, _, max_slew_q32 = PROFILES[profile].fixed_coefficients(
        update_period_s
    )
    phase_q32 = round(
        sample.hardware_phase_ns * 1.0e-9 * TRACKING_HZ * Q32
    )
    integrator_q32 = round(sample.integrator_ppm * 1.0e-6 * Q32)
    proportional_q32 = (phase_q32 * kp_q46) >> 46
    rate_q32 = max(
        -max_slew_q32,
        min(max_slew_q32, integrator_q32 + proportional_q32),
    )
    return rate_q32 / Q32 * 1.0e6


def _reconstructed_active_phase_ns(
    samples: list[FirmwarePpsSample], oscillator_ppm: float
) -> list[float]:
    """Integrate observed rate under a constant-oscillator assumption."""
    phase_ns = [0.0]
    for first, second in zip(samples, samples[1:]):
        elapsed_s = second.time_s - first.time_s
        mean_command_ppm = 0.5 * (
            first.rate_command_ppm + second.rate_command_ppm
        )
        phase_ns.append(
            phase_ns[-1]
            + (oscillator_ppm + mean_command_ppm) * elapsed_s * 1.0e3
        )
    return phase_ns


def _parse_float(group: dict[str, str], key: str) -> float:
    return float(group[key])


def _parse_int(group: dict[str, str], key: str) -> int:
    return int(group[key])


def _infer_update_period(
    pending: list[dict[str, object]],
    override_s: float | None,
) -> float:
    if override_s is not None:
        if not math.isfinite(override_s) or override_s <= 0.0:
            raise ValueError("update period override must be positive")
        return override_s
    ratios = []
    for first, second in zip(pending, pending[1:]):
        sequence_step = (
            int(second["sequence"]) - int(first["sequence"])
        ) & 0xFFFF_FFFF
        host_step = (
            float(second["host_monotonic_s"])
            - float(first["host_monotonic_s"])
        )
        if 0 < sequence_step < 10_000 and host_step > 0.0:
            ratios.append(host_step / sequence_step)
    if not ratios:
        return UPDATE_PERIOD_S
    inferred_s = statistics.median(ratios)
    snapped_s = max(
        NAV_BIT_PERIOD_S,
        round(inferred_s / NAV_BIT_PERIOD_S) * NAV_BIT_PERIOD_S,
    )
    return (
        snapped_s
        if abs(snapped_s - inferred_s) <= max(0.005, 0.1 * snapped_s)
        else inferred_s
    )


def _sample_update_period(samples: list[FirmwarePpsSample]) -> float:
    ratios = []
    for first, second in zip(samples, samples[1:]):
        sequence_step = (
            second.sequence - first.sequence
        ) & 0xFFFF_FFFF
        elapsed_s = second.time_s - first.time_s
        if 0 < sequence_step < 10_000 and elapsed_s > 0.0:
            ratios.append(elapsed_s / sequence_step)
    return statistics.median(ratios) if ratios else UPDATE_PERIOD_S


def load_firmware_telemetry(
    path: Path, *, update_period_s: float | None = None
) -> list[FirmwarePpsSample]:
    """Load host-timestamped UART CSV produced by ``serial_capture.py``."""
    pending: list[dict[str, object]] = []
    latest_location: GpsContext | None = None
    latest_tracked = 0
    with path.open("r", encoding="utf-8", errors="replace", newline="") as source:
        for row in csv.DictReader(source):
            try:
                host_monotonic_s = float(row["host_monotonic_s"])
                host_unix_s = float(row["host_unix_s"])
                line = row["uart_line"]
            except (KeyError, TypeError, ValueError):
                continue

            tracked_match = TRACKED_RE.search(line)
            if tracked_match:
                latest_tracked = int(tracked_match.group("tracked"))

            location_match = LOCATION_RE.search(line)
            if location_match:
                groups = location_match.groupdict()
                latest_location = GpsContext(
                    host_monotonic_s,
                    int(groups["tracked"]),
                    int(groups["used"]),
                    float(groups["hdop"]),
                    float(groups["vdop"]),
                    bool(int(groups["velocity_valid"])),
                    bool(int(groups["velocity_fde"])),
                    float(groups["velocity_rms_mps"]),
                )

            pps_match = PPS_RE.search(line)
            if not pps_match:
                continue
            groups = pps_match.groupdict()
            context = latest_location
            context_fresh = (
                context is not None
                and host_monotonic_s - context.host_monotonic_s <= 1.0
            )
            pending.append(
                {
                    "host_monotonic_s": host_monotonic_s,
                    "host_unix_s": host_unix_s,
                    "profile": groups["profile"],
                    "mode": groups["mode"],
                    "locked": bool(_parse_int(groups, "locked")),
                    "tdop": _parse_float(groups, "tdop"),
                    "gdop": _parse_float(groups, "gdop"),
                    "hdop": context.hdop if context_fresh else math.nan,
                    "vdop": context.vdop if context_fresh else math.nan,
                    "clock_sigma_m": _parse_float(groups, "clock_sigma_m"),
                    "cn0_mean_dbhz": _parse_float(groups, "cn0_mean_dbhz"),
                    "cn0_min_dbhz": _parse_float(groups, "cn0_min_dbhz"),
                    "tracked_satellites": (
                        context.tracked_satellites
                        if context_fresh
                        else latest_tracked
                    ),
                    "used_satellites": (
                        context.used_satellites if context_fresh else 0
                    ),
                    "innovation_m": _parse_float(groups, "innovation_m"),
                    "raw_bias_m": _parse_float(groups, "raw_bias_m"),
                    "filtered_bias_m": _parse_float(groups, "filtered_bias_m"),
                    "clock_drift_mps": _parse_float(groups, "clock_drift_mps"),
                    "bias_drift_mps": _parse_float(groups, "bias_drift_mps"),
                    "predictor_enabled": bool(_parse_int(groups, "predictor")),
                    "eskf_alpha": _parse_float(groups, "eskf_alpha"),
                    "eskf_speed_mps": _parse_float(groups, "eskf_speed_mps"),
                    "eskf_velocity_rms_mps": _parse_float(
                        groups, "eskf_velocity_rms_mps"
                    ),
                    "eskf_fde_mask": int(groups["eskf_fde_mask"], 16),
                    "pps_speed_mps": _parse_float(groups, "pps_speed_mps"),
                    "pps_velocity_rms_mps": _parse_float(
                        groups, "pps_velocity_rms_mps"
                    ),
                    "pps_fde_mask": int(groups["pps_fde_mask"], 16),
                    "reference_svid": _parse_int(groups, "reference_svid"),
                    "capture_correction_ns": _parse_float(
                        groups, "capture_correction_ns"
                    ),
                    "hardware_phase_ns": _parse_float(
                        groups, "hardware_phase_ns"
                    ),
                    "rate_command_ppm": _parse_float(groups, "rate_command_ppm"),
                    "integrator_ppm": _parse_float(groups, "integrator_ppm"),
                    "pulse_count": _parse_int(groups, "pulse_count"),
                    "sequence": _parse_int(groups, "sequence"),
                    "status": int(groups["status"], 16),
                }
            )

    if not pending:
        return []
    update_period_s = _infer_update_period(pending, update_period_s)
    first_sequence = int(pending[0]["sequence"])
    samples: list[FirmwarePpsSample] = []
    for fields in pending:
        sequence = int(fields["sequence"])
        sequence_delta = (sequence - first_sequence) & 0xFFFF_FFFF
        fields["time_s"] = sequence_delta * update_period_s
        samples.append(FirmwarePpsSample(**fields))
    return samples


def observed_metrics(samples: list[FirmwarePpsSample]) -> ObservedMetrics:
    if len(samples) < 2:
        raise ValueError("at least two firmware PPS samples are required")
    phase = [sample.hardware_phase_ns for sample in samples]
    phase_center = statistics.median(phase)
    centered_phase = [value - phase_center for value in phase]
    phase_steps = [
        phase[index] - phase[index - 1] for index in range(1, len(phase))
    ]
    host_intervals = [
        samples[index].host_monotonic_s - samples[index - 1].host_monotonic_s
        for index in range(1, len(samples))
    ]
    sequence_steps = [
        (samples[index].sequence - samples[index - 1].sequence) & 0xFFFF_FFFF
        for index in range(1, len(samples))
    ]
    _, phase_slope = _linear_fit(
        [sample.time_s for sample in samples], phase
    )
    _, raw_slope = _linear_fit(
        [sample.time_s for sample in samples],
        [sample.raw_bias_m for sample in samples],
    )
    _, filtered_slope = _linear_fit(
        [sample.time_s for sample in samples],
        [sample.filtered_bias_m for sample in samples],
    )
    profiles = [sample.profile for sample in samples]
    references = [sample.reference_svid for sample in samples]
    innovations = [sample.innovation_m for sample in samples]
    rate = [sample.rate_command_ppm for sample in samples]
    status_interval = statistics.median(host_intervals)
    wander_window = max(1, round(10.0 / status_interval))
    phase_10s_means = [
        statistics.mean(centered_phase[index - wander_window + 1:index + 1])
        for index in range(wander_window - 1, len(centered_phase))
    ]
    inferred_oscillator_ppm = -statistics.mean(rate)
    inferred_active_phase = _reconstructed_active_phase_ns(
        samples, inferred_oscillator_ppm
    )
    active_center = statistics.median(inferred_active_phase)
    centered_active = [
        value - active_center for value in inferred_active_phase
    ]
    update_period_s = _sample_update_period(samples)
    pi_rate_errors_ppb = [
        (prediction - sample.rate_command_ppm) * 1.0e3
        for sample in samples
        if math.isfinite(
            prediction := _pi_rate_prediction_ppm(sample, update_period_s)
        )
    ]
    return ObservedMetrics(
        len(samples),
        samples[-1].time_s - samples[0].time_s,
        statistics.median(host_intervals),
        statistics.median(sequence_steps),
        statistics.mean(sample.locked for sample in samples),
        sum(
            profiles[index] != profiles[index - 1]
            for index in range(1, len(profiles))
        ),
        sum(
            references[index] != references[index - 1]
            for index in range(1, len(references))
        ),
        statistics.mean(phase),
        _rms(centered_phase),
        _percentile((abs(value) for value in centered_phase), 0.95),
        max(abs(value) for value in centered_phase),
        max(phase) - min(phase),
        (
            max(phase_10s_means) - min(phase_10s_means)
            if phase_10s_means else math.nan
        ),
        _rms(phase_steps),
        _percentile((abs(value) for value in phase_steps), 0.95),
        max(abs(value) for value in phase_steps),
        phase_slope,
        statistics.median(rate),
        max(rate) - min(rate),
        statistics.median(sample.integrator_ppm for sample in samples),
        inferred_oscillator_ppm,
        _rms(centered_active),
        max(inferred_active_phase) - min(inferred_active_phase),
        _correlation(phase, rate),
        _correlation(phase, innovations),
        _rms(pi_rate_errors_ppb),
        _percentile((abs(value) for value in pi_rate_errors_ppb), 0.95),
        max((abs(value) for value in pi_rate_errors_ppb), default=math.nan),
        _rms(innovations),
        _percentile((abs(value) for value in innovations), 0.95),
        raw_slope,
        filtered_slope,
        statistics.median(sample.tdop for sample in samples),
        statistics.median(sample.gdop for sample in samples),
        statistics.median(sample.clock_sigma_m for sample in samples),
        statistics.median(sample.cn0_mean_dbhz for sample in samples),
        statistics.median(sample.cn0_min_dbhz for sample in samples),
        statistics.median(sample.tracked_satellites for sample in samples),
        statistics.median(sample.used_satellites for sample in samples),
        statistics.median(sample.pps_velocity_rms_mps for sample in samples),
    )


def _lerp(first: float, second: float, fraction: float) -> float:
    return first + (second - first) * fraction


def _quality(sample: FirmwarePpsSample) -> QualityMetrics:
    # chi-square is not present in the level-1 status record.  A finite value
    # within the observed profile's accepted band preserves the firmware's
    # already-reported robust/precision classification without pretending it
    # was measured by this parser.
    inferred_chi_square = 0.5
    solution_valid = sample.profile != "holdover" and sample.used_satellites >= 4
    return QualityMetrics(
        sample.tdop,
        sample.gdop,
        sample.clock_sigma_m,
        sample.cn0_mean_dbhz,
        sample.cn0_min_dbhz,
        inferred_chi_square,
        sample.pps_velocity_rms_mps,
        sample.used_satellites,
        solution_valid=solution_valid,
        residual_test_valid=solution_valid,
        velocity_valid=solution_valid,
        fde_applied=sample.pps_fde_mask != 0,
    )


def _interpolated_sample(
    first: FirmwarePpsSample,
    second: FirmwarePpsSample,
    step: int,
    total_steps: int,
) -> FirmwarePpsSample:
    fraction = step / total_steps
    fields = asdict(first)
    linear_fields = (
        "host_monotonic_s", "host_unix_s", "tdop", "gdop", "hdop", "vdop",
        "clock_sigma_m", "cn0_mean_dbhz", "cn0_min_dbhz", "innovation_m",
        "raw_bias_m", "filtered_bias_m", "clock_drift_mps", "bias_drift_mps",
        "eskf_alpha", "eskf_speed_mps", "eskf_velocity_rms_mps",
        "pps_speed_mps", "pps_velocity_rms_mps", "capture_correction_ns",
        "hardware_phase_ns", "rate_command_ppm", "integrator_ppm",
    )
    for name in linear_fields:
        fields[name] = _lerp(
            float(getattr(first, name)), float(getattr(second, name)), fraction
        )
    fields["sequence"] = (first.sequence + step) & 0xFFFF_FFFF
    fields["time_s"] = _lerp(first.time_s, second.time_s, fraction)
    fields["pulse_count"] = round(
        _lerp(first.pulse_count, second.pulse_count, fraction)
    )
    if fraction >= 0.5:
        for name in (
            "profile", "mode", "locked", "tracked_satellites",
            "used_satellites", "predictor_enabled", "eskf_fde_mask",
            "pps_fde_mask", "reference_svid", "status",
        ):
            fields[name] = getattr(second, name)
    return FirmwarePpsSample(**fields)


def expand_samples(samples: list[FirmwarePpsSample]) -> list[tuple[FirmwarePpsSample, bool]]:
    if not samples:
        return []
    expanded: list[tuple[FirmwarePpsSample, bool]] = [(samples[0], False)]
    for first, second in zip(samples, samples[1:]):
        steps = (second.sequence - first.sequence) & 0xFFFF_FFFF
        if steps == 0 or steps > 10_000:
            continue
        for step in range(1, steps + 1):
            expanded.append(
                (
                    second if step == steps else _interpolated_sample(
                        first, second, step, steps
                    ),
                    step != steps,
                )
            )
    return expanded


def replay_firmware_telemetry(
    samples: list[FirmwarePpsSample],
    *,
    warmup_s: float = 20.0,
    tuning: str = "current",
    profiles: dict[Profile, LoopProfile] | None = None,
    software_drift_gain: float = 1.0 / 4.0,
) -> tuple[list[ReplayPoint], ReplayMetrics]:
    if len(samples) < 2:
        raise ValueError("at least two firmware PPS samples are required")
    expanded = expand_samples(samples)
    update_period_s = _sample_update_period(samples)
    first = expanded[0][0]
    first_profile = _profile(first)
    if not first.locked or first_profile == Profile.HOLDOVER:
        raise ValueError(
            "hardware replay requires an initially locked tracking sample"
        )

    # The UART does not contain a direct oscillator measurement.  Choosing
    # the negative mean command is the unique constant-oscillator convention
    # that gives zero net frequency over the capture.  Integrating the logged
    # command then reconstructs the board's relative active-phase trajectory.
    oscillator_ppm = -statistics.mean(
        sample.rate_command_ppm for sample in samples
    )
    expanded_samples = [sample for sample, _ in expanded]
    active_phase_ns = _reconstructed_active_phase_ns(
        expanded_samples, oscillator_ppm
    )
    reference_timing_s = [
        phase_ns * 1.0e-9 + sample.hardware_phase_ns * 1.0e-9
        for phase_ns, sample in zip(active_phase_ns, expanded_samples)
    ]
    # Firmware's filtered bias is the value used to form the PPS reference.
    # Reversing that operation gives a self-consistent datum trajectory while
    # retaining raw WLS bias as the actual software-filter input.
    true_bias_m = [
        sample.filtered_bias_m + timing_s * SPEED_OF_LIGHT_MPS
        for sample, timing_s in zip(expanded_samples, reference_timing_s)
    ]
    model = FixedPointAdaptivePps(
        update_period_s=update_period_s,
        oscillator_ppm=oscillator_ppm,
        receiver_clock_bias_m=true_bias_m[0],
        receiver_clock_drift_mps=0.0,
        profiles=profiles,
        software_drift_gain=software_drift_gain,
    )
    model.seed_tracking_state(
        time_s=first.time_s,
        profile=first_profile,
        phase_error_s=active_phase_ns[0] * 1.0e-9,
        rate_command_ppm=first.rate_command_ppm,
        integrator_ppm=first.integrator_ppm,
        software_error_s=reference_timing_s[0],
        software_bias_m=first.filtered_bias_m,
        software_clock_drift_mps=first.bias_drift_mps,
        software_previous_measurement_m=first.raw_bias_m,
        reference_id=first.reference_svid,
    )

    points: list[ReplayPoint] = [
        ReplayPoint(
            first.time_s,
            first.sequence,
            False,
            first.raw_bias_m,
            first.filtered_bias_m,
            first.filtered_bias_m,
            active_phase_ns[0],
            active_phase_ns[0],
            first.hardware_phase_ns,
            first.hardware_phase_ns,
            first.rate_command_ppm,
            first.rate_command_ppm,
        )
    ]
    original_sequences = {sample.sequence for sample in samples}
    observed_filtered = {sample.sequence: sample.filtered_bias_m for sample in samples}
    for index, (sample, interpolated) in enumerate(expanded[1:], start=1):
        observation = Observation(
            sample.time_s,
            (true_bias_m[index] - sample.raw_bias_m) / SPEED_OF_LIGHT_MPS,
            sample.bias_drift_mps,
            _quality(sample),
            clock_drift_valid=sample.predictor_enabled,
            reference_range_rate_mps=0.0,
            reference_id=sample.reference_svid,
            reference_datum_offset_m=(
                true_bias_m[index] - true_bias_m[0]
            ),
        )
        model.update(observation)
        modeled_loop_error_ns = model.trace.raw_loop_error_s[-1] * 1.0e9
        points.append(
            ReplayPoint(
                sample.time_s,
                sample.sequence,
                interpolated,
                sample.raw_bias_m,
                observed_filtered.get(sample.sequence, math.nan),
                model.software_bias_m,
                active_phase_ns[index],
                model.trace.phase_error_s[-1] * 1.0e9,
                sample.hardware_phase_ns if sample.sequence in original_sequences else math.nan,
                modeled_loop_error_ns,
                sample.rate_command_ppm if sample.sequence in original_sequences else math.nan,
                model.trace.frequency_command_ppm[-1],
            )
        )
    trace = model.finish(expanded[-1][0].time_s + update_period_s)

    edge_errors_ns = [
        error * 1.0e9
        for time_s, error in trace.integer_second_edge_errors(after_s=warmup_s)
    ]
    edge_center = (
        statistics.median(edge_errors_ns) if edge_errors_ns else math.nan
    )
    centered_edges = [value - edge_center for value in edge_errors_ns]
    periods = [
        edge_errors_ns[index] - edge_errors_ns[index - 1]
        for index in range(1, len(edge_errors_ns))
    ]
    locked_after_warmup = [
        locked
        for time_s, locked in zip(trace.time_s, trace.locked)
        if time_s >= warmup_s
    ]
    bias_errors = [
        point.modeled_filtered_bias_m - point.observed_filtered_bias_m
        for point in points
        if not point.interpolated
        and point.time_s >= warmup_s
        and math.isfinite(point.observed_filtered_bias_m)
    ]
    modeled_residuals = [
        point.modeled_loop_error_ns
        for point in points
        if point.time_s >= warmup_s
    ]
    modeled_residual_center = (
        statistics.median(modeled_residuals)
        if modeled_residuals else math.nan
    )
    centered_modeled_residuals = [
        value - modeled_residual_center for value in modeled_residuals
    ]
    residual_reproduction_errors = [
        point.modeled_loop_error_ns - point.observed_hardware_phase_ns
        for point in points
        if not point.interpolated
        and point.time_s >= warmup_s
        and math.isfinite(point.observed_hardware_phase_ns)
    ]
    metrics = ReplayMetrics(
        tuning,
        len(expanded),
        len(expanded) / len(samples),
        0.0,
        statistics.mean(locked_after_warmup) if locked_after_warmup else 0.0,
        trace.unlock_events,
        oscillator_ppm,
        _rms(centered_edges),
        _percentile((abs(value) for value in centered_edges), 0.95),
        max(centered_edges) - min(centered_edges) if centered_edges else math.nan,
        _rms(periods),
        _percentile((abs(value) for value in periods), 0.95),
        max((abs(value) for value in periods), default=math.nan),
        _rms(centered_modeled_residuals),
        _percentile(
            (abs(value) for value in centered_modeled_residuals), 0.95
        ),
        (
            max(centered_modeled_residuals)
            - min(centered_modeled_residuals)
            if centered_modeled_residuals else math.nan
        ),
        _rms(residual_reproduction_errors),
        _percentile(
            (abs(value) for value in residual_reproduction_errors), 0.95
        ),
        _rms(bias_errors),
        _percentile((abs(value) for value in bias_errors), 0.95),
    )
    return points, metrics


def replay_tuning_sensitivity(
    samples: list[FirmwarePpsSample], *, warmup_s: float = 20.0
) -> list[ReplayMetrics]:
    """Run small, explicitly non-production tuning counterfactuals."""
    _, current = replay_firmware_telemetry(
        samples, warmup_s=warmup_s, tuning="current"
    )
    guarded_profiles = dict(PROFILES)
    guarded_profiles[Profile.ROBUST] = replace(
        guarded_profiles[Profile.ROBUST], software_bias_gain=1.0 / 64.0
    )
    _, guarded = replay_firmware_telemetry(
        samples,
        warmup_s=warmup_s,
        tuning="robust_bias_1_64",
        profiles=guarded_profiles,
    )
    _, faster_drift = replay_firmware_telemetry(
        samples,
        warmup_s=warmup_s,
        tuning="drift_gain_1_2",
        software_drift_gain=1.0 / 2.0,
    )
    return [current, guarded, faster_drift]


def write_samples(path: Path, samples: list[FirmwarePpsSample]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as output:
        writer = csv.DictWriter(
            output,
            fieldnames=list(FirmwarePpsSample.__annotations__),
            lineterminator="\n",
        )
        writer.writeheader()
        for sample in samples:
            writer.writerow(asdict(sample))


def write_replay(path: Path, points: list[ReplayPoint]) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("w", encoding="utf-8", newline="") as output:
        writer = csv.DictWriter(
            output,
            fieldnames=list(ReplayPoint.__annotations__),
            lineterminator="\n",
        )
        writer.writeheader()
        for point in points:
            writer.writerow(asdict(point))


def write_report(
    path: Path,
    source: Path,
    samples: list[FirmwarePpsSample],
    observed: ObservedMetrics,
    replay: ReplayMetrics,
    sensitivity: list[ReplayMetrics],
) -> None:
    update_period_s = _sample_update_period(samples)
    profiles = sorted({sample.profile for sample in samples})
    references = sorted({sample.reference_svid for sample in samples})
    lines = [
        "# Genesys2 live PPS telemetry replay",
        "",
        f"Source: `{source}`",
        "",
        "## Direct firmware/RTL observations",
        "",
        f"- Samples/duration: {observed.samples} / {observed.duration_s:.1f} s.",
        f"- Status interval: {observed.status_interval_median_s:.3f} s; "
        f"median sequence step {observed.sequence_step_median:.0f} "
        f"({observed.sequence_step_median * update_period_s:.1f} s of "
        f"{update_period_s * 1000.0:.0f} ms solves).",
        f"- Profiles: {', '.join(profiles)}; locked fraction "
        f"{observed.locked_fraction:.3f}; transitions "
        f"{observed.profile_transitions}; reference SVIDs {references}.",
        f"- TDOP/GDOP medians: {observed.tdop_median:.3f} / "
        f"{observed.gdop_median:.3f}; clock sigma median "
        f"{observed.clock_sigma_median_m:.2f} m.",
        f"- C/N0 mean/min medians: {observed.cn0_mean_median_dbhz:.1f} / "
        f"{observed.cn0_min_median_dbhz:.1f} dB-Hz; tracked/used medians "
        f"{observed.tracked_median:.0f}/{observed.used_median:.0f}.",
        f"- WLS innovation RMS/p95-absolute: {observed.innovation_rms_m:.2f} / "
        f"{observed.innovation_p95_m:.2f} m.",
        f"- Raw/filtered bias slopes: {observed.raw_bias_slope_mps:.3f} / "
        f"{observed.filtered_bias_slope_mps:.3f} m/s; PPS velocity residual "
        f"median {observed.pps_velocity_rms_median_mps:.3f} m/s.",
        f"- Internal hardware phase mean: {observed.hardware_phase_mean_ns:.2f} ns; "
        f"median-centered RMS/p95/max/p-p: "
        f"{observed.hardware_phase_centered_rms_ns:.2f} / "
        f"{observed.hardware_phase_centered_p95_ns:.2f} / "
        f"{observed.hardware_phase_centered_max_ns:.2f} / "
        f"{observed.hardware_phase_peak_to_peak_ns:.2f} ns.",
        f"- Ten-second averaged phase wander p-p: "
        f"{observed.hardware_phase_10s_wander_peak_to_peak_ns:.2f} ns.",
        f"- Sampled phase-step RMS/p95/max: "
        f"{observed.sampled_phase_step_rms_ns:.2f} / "
        f"{observed.sampled_phase_step_p95_ns:.2f} / "
        f"{observed.sampled_phase_step_max_ns:.2f} ns. These are at the status "
        "interval, not consecutive physical 1 PPS edges.",
        f"- Rate command median/span: {observed.rate_command_median_ppm:.4f} / "
        f"{observed.rate_command_span_ppm:.4f} ppm; integrator median "
        f"{observed.integrator_median_ppm:.4f} ppm.",
        f"- Fixed-point PI rate reproduction RMS/p95/max: "
        f"{observed.pi_rate_reproduction_rms_ppb:.3f} / "
        f"{observed.pi_rate_reproduction_p95_ppb:.3f} / "
        f"{observed.pi_rate_reproduction_max_ppb:.3f} ppb. This is the "
        "difference between logged `rate_cmd` and `integrator + Kp*hw_phase`; "
        "the remaining error is dominated by three-decimal UART rounding.",
        f"- Phase correlation with rate command / WLS innovation: "
        f"{observed.phase_rate_correlation:.3f} / "
        f"{observed.phase_innovation_correlation:.3f}.",
        f"- A constant oscillator inferred from zero net frequency is "
        f"{observed.inferred_oscillator_ppm:.6f} ppm. Integrating the logged "
        f"command under that assumption gives centered phase RMS/p-p "
        f"{observed.inferred_active_phase_rms_ns:.2f} / "
        f"{observed.inferred_active_phase_peak_to_peak_ns:.2f} ns.",
        "",
        "`hw_phase` is the RTL PI residual captured over APB. It is not a "
        "trusted-generator time-interval measurement of the output pin, so "
        "these values cannot certify physical PPS jitter or absolute offset.",
        "",
        "## Telemetry-conditioned end-to-end loop replay",
        "",
        "The replay imports the first locked phase/rate/integrator and software "
        "filter snapshot instead of incorrectly resetting the board to an "
        "acquisition state. It uses raw/filtered clock bias, Doppler datum "
        "drift, TDOP/GDOP, clock sigma, C/N0, satellite count, FDE, reference "
        "SVID and the production fixed-point coefficients. The reference "
        "trajectory is reconstructed from filtered bias, the logged PI "
        "residual and the integrated rate command.",
        "",
        f"- Expanded updates: {replay.expanded_updates}; interpolation factor "
        f"{replay.interpolation_factor:.2f}.",
        f"- Snapshot starts locked at {replay.first_lock_s:.2f} s; post-warmup locked "
        f"fraction {replay.locked_fraction_after_warmup:.3f}; unlocks "
        f"{replay.unlocks}.",
        f"- Constant-oscillator convention: "
        f"{replay.inferred_oscillator_ppm:.6f} ppm.",
        f"- Modeled physical-grid edge centered RMS/p95/p-p: "
        f"{replay.modeled_edge_centered_rms_ns:.2f} / "
        f"{replay.modeled_edge_centered_p95_ns:.2f} / "
        f"{replay.modeled_edge_peak_to_peak_ns:.2f} ns.",
        f"- Modeled one-second period RMS/p95/max: "
        f"{replay.modeled_period_rms_ns:.2f} / "
        f"{replay.modeled_period_p95_ns:.2f} / "
        f"{replay.modeled_period_max_ns:.2f} ns.",
        f"- Modeled PI residual centered RMS/p95/p-p: "
        f"{replay.modeled_loop_residual_centered_rms_ns:.2f} / "
        f"{replay.modeled_loop_residual_centered_p95_ns:.2f} / "
        f"{replay.modeled_loop_residual_peak_to_peak_ns:.2f} ns.",
        f"- Logged-epoch PI residual reproduction RMS/p95 error: "
        f"{replay.logged_residual_reproduction_rms_ns:.2f} / "
        f"{replay.logged_residual_reproduction_p95_ns:.2f} ns.",
        f"- Logged-epoch filtered-bias reproduction RMS/p95: "
        f"{replay.filtered_bias_reproduction_rms_m:.2f} / "
        f"{replay.filtered_bias_reproduction_p95_m:.2f} m.",
        "",
        "## Tuning sensitivity on the same trace",
        "",
        "| tuning | edge RMS (ns) | edge p-p (ns) | period max (ns) | "
        "residual reproduction RMS (ns) | unlocks |",
        "|---|---:|---:|---:|---:|---:|",
    ]
    for item in sensitivity:
        lines.append(
            f"| {item.tuning} | {item.modeled_edge_centered_rms_ns:.2f} | "
            f"{item.modeled_edge_peak_to_peak_ns:.2f} | "
            f"{item.modeled_period_max_ns:.2f} | "
            f"{item.logged_residual_reproduction_rms_ns:.2f} | "
            f"{item.unlocks} |"
        )
    lines.extend([
        "",
        "These counterfactuals reuse a reference trajectory reconstructed "
        "from the current filter, and nine of every ten solve epochs are "
        "interpolated. They are useful as a rejection test, not as evidence "
        "to deploy new gains. On this capture neither `robust_bias_1_64` nor "
        "`drift_gain_1_2` improves the current loop's absolute wander.",
        "",
        "## Limits and next measurement",
        "",
        "- Missing status epochs are linearly interpolated when sequence step "
        "is greater than one. This preserves low-frequency wander but "
        f"understates unlogged {update_period_s * 1000.0:.0f} ms noise; it is "
        "not an exact firmware replay.",
        "- The inferred oscillator and active phase assume one constant crystal "
        "frequency chosen to make net phase slope zero over this capture. "
        "Temperature telemetry was unavailable, so slow crystal motion is not "
        "separable from the reconstructed reference.",
        "- Level-1 telemetry does not expose ABS_SAMPLE, REM_CODE_PHASE, full "
        "capture counter, chi-square ratio, pseudorange residual vector, board "
        "temperature, or a trusted PPS edge timestamp. The bit/subframe part "
        "of the end-to-end model remains synthetic for those fields.",
        "- Set `GPS_PPS_STATUS_LOG_INTERVAL=1` for exact solve-cadence input. "
        "For a max-20 ns pin claim, capture DUT and trusted PPS with a "
        "common-clock TIC/scope after cable/output-delay calibration.",
        "",
    ])
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text("\n".join(lines), encoding="utf-8")


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--input", type=Path, required=True)
    parser.add_argument(
        "--output-dir",
        type=Path,
        default=Path(__file__).resolve().parent / "output/hardware_replay",
    )
    parser.add_argument("--warmup", type=float, default=20.0)
    parser.add_argument(
        "--update-period-ms",
        type=float,
        help="override NAV solve period; default infers it from host/sequence",
    )
    parser.add_argument("--no-write", action="store_true")
    return parser.parse_args()


def main() -> int:
    args = parse_args()
    samples = load_firmware_telemetry(
        args.input,
        update_period_s=(
            args.update_period_ms / 1000.0
            if args.update_period_ms is not None
            else None
        ),
    )
    if len(samples) < 2:
        raise SystemExit("input does not contain at least two firmware PPS records")
    observed = observed_metrics(samples)
    points, replay = replay_firmware_telemetry(samples, warmup_s=args.warmup)
    sensitivity = replay_tuning_sensitivity(
        samples, warmup_s=args.warmup
    )
    print(
        f"samples={observed.samples} duration={observed.duration_s:.1f}s "
        f"locked={observed.locked_fraction:.3f} "
        f"hw_phase_centered_rms={observed.hardware_phase_centered_rms_ns:.2f}ns "
        f"hw_phase_p95={observed.hardware_phase_centered_p95_ns:.2f}ns "
        f"hw_phase_pp={observed.hardware_phase_peak_to_peak_ns:.2f}ns"
    )
    print(
        f"replay_updates={replay.expanded_updates} "
        f"model_edge_rms={replay.modeled_edge_centered_rms_ns:.2f}ns "
        f"model_period_rms={replay.modeled_period_rms_ns:.2f}ns "
        f"model_period_max={replay.modeled_period_max_ns:.2f}ns "
        f"unlocks={replay.unlocks}"
    )
    if not args.no_write:
        write_samples(args.output_dir / "firmware_telemetry.csv", samples)
        write_replay(args.output_dir / "model_replay.csv", points)
        write_report(
            args.output_dir / "hardware_replay_report.md",
            args.input,
            samples,
            observed,
            replay,
            sensitivity,
        )
        print(f"wrote {args.output_dir}")
    return 0


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