#!/usr/bin/env python3
"""Score AI load cycles against the furnace log by PHASE, not by event matching.

The furnace is close to a metronome. In the pilot hour its charge interval was
212-229 s, and load cycles arrive at a similar rate, so a +-90 s matching window
gives ~0.77 recall to *any* detector firing every ~210 s -- including one firing
at random. Measured on the pilot: recall 0.938 at the true offset vs 0.773 mean
at |offset| >= 180 s, and several wrong offsets scored 1.000. Recall/precision
therefore say almost nothing here, and the old 0.84 headline should not be
quoted.

What does carry information is PHASE. If the detector really sees loading, each
AI cycle should land at a repeatable point inside the furnace's charge-to-charge
interval. If it is noise, phases spread uniformly. Because the interval is not
constant across the day (breaks stretch it), phase is normalised per interval:

    phase = (t_ai - charge_prev) / (charge_next - charge_prev)   in [0, 1)

then tested for concentration with the Rayleigh test. `--control` re-runs the
same test on random timestamps to confirm the test is calibrated -- it should
reject at about the nominal 5%.

Pilot result for reference: R=0.754, p<0.001, mean phase +41 s after charge,
circular SD 26 s; the random control passed 5.5% of the time.
"""

from __future__ import annotations

import argparse
import json
from pathlib import Path

import numpy as np


def phase_of(events, charges):
    """Normalised position of each event inside its enclosing charge interval."""
    ch = np.asarray(sorted(charges), float)
    out, spans = [], []
    for t in events:
        i = np.searchsorted(ch, t) - 1
        if i < 0 or i + 1 >= len(ch):
            continue
        a, b = ch[i], ch[i + 1]
        if b <= a:
            continue
        out.append((t - a) / (b - a))
        spans.append(b - a)
    return np.asarray(out), np.asarray(spans)


def rayleigh(phases):
    """Return (R, p, mean_phase). R=0 uniform, R=1 perfectly locked."""
    if len(phases) < 3:
        return 0.0, 1.0, float("nan")
    ang = 2 * np.pi * phases
    C, S = np.cos(ang).mean(), np.sin(ang).mean()
    R = float(np.hypot(C, S))
    n = len(phases)
    Z = n * R * R
    p = float(np.exp(-Z) * (1 + (2 * Z - Z * Z) / (4 * n)))
    p = min(max(p, 0.0), 1.0)
    mu = float((np.degrees(np.arctan2(S, C)) % 360) / 360.0)
    return R, p, mu


def circ_sd_frac(R):
    return float(np.sqrt(-2 * np.log(max(R, 1e-9))) / (2 * np.pi))


def main() -> int:
    ap = argparse.ArgumentParser()
    ap.add_argument("--cycles", required=True, help="load_cycles_v2.json")
    ap.add_argument("--furnace", required=True, help="furnace_20260617.json")
    ap.add_argument("--clip", required=True)
    ap.add_argument("--video-t0", type=int, default=0, help="clip-local offset of the analysed window")
    ap.add_argument("--key", default="t_start", choices=["t_start", "t_peak", "t_end"])
    ap.add_argument("--win-start", type=int, default=None, help="clip-local start of the analysed window")
    ap.add_argument("--win-end", type=int, default=None)
    ap.add_argument("--control", type=int, default=200)
    ap.add_argument("--out", default="")
    args = ap.parse_args()

    cyc = json.loads(Path(args.cycles).read_text(encoding="utf-8"))["cycles"]
    fur = json.loads(Path(args.furnace).read_text(encoding="utf-8"))
    clip = fur["per_clip"][args.clip]
    charges = clip["charge_video_t"]

    ev = [c[args.key] + args.video_t0 for c in cyc]
    if not ev:
        raise SystemExit("no AI cycles")
    # Count comparison must use the analysed window, not the span of the AI
    # events -- padding it lets extra charges in and fakes a count error.
    # Both sides must be clipped: filtering only the charges compares a whole
    # clip's cycles against one hour's batches.
    lo = args.win_start if args.win_start is not None else min(ev)
    hi = args.win_end if args.win_end is not None else max(ev)
    ev = [t for t in ev if lo <= t < hi]
    if not ev:
        raise SystemExit("no AI cycles inside the requested window")
    in_win = [g for g in charges if lo <= g < hi]

    ph, spans = phase_of(ev, charges)
    R, p, mu = rayleigh(ph)
    period = float(np.median(spans)) if len(spans) else float("nan")

    rng = np.random.default_rng(0)
    npass = 0
    for _ in range(args.control):
        fake = rng.uniform(lo, hi, len(ev))
        fp, _ = phase_of(fake, charges)
        _, pc, _ = rayleigh(fp)
        npass += pc < 0.05
    ctrl = npass / max(args.control, 1)

    rep = {
        "clip": args.clip,
        "anchor": args.key,
        "n_ai_cycles": len(ev),
        "n_furnace_charges_in_window": len(in_win),
        "count_ratio": round(len(ev) / len(in_win), 3) if in_win else None,
        "count_error_pct": round(100 * (len(ev) - len(in_win)) / len(in_win), 1) if in_win else None,
        "n_phased": int(len(ph)),
        "median_interval_sec": round(period, 1),
        "phase_R": round(R, 4),
        "rayleigh_p": round(p, 6),
        "mean_offset_sec_after_charge": round(mu * period, 1) if period == period else None,
        "circular_sd_sec": round(circ_sd_frac(R) * period, 1) if period == period else None,
        "control_false_positive_rate": round(ctrl, 3),
        "locked": bool(p < 0.05),
        # raw phases so the viewer can draw the concentration directly; a locked
        # detector shows one clear lobe, noise shows a flat ring
        "phases": [round(float(v), 4) for v in ph],
        "offsets_sec": [round(float(v * s), 1) for v, s in zip(ph, spans)],
        "interpretation": (
            "AI cycles are phase-locked to the machine-logged furnace charges; "
            "the random-timestamp control rejects at the nominal rate, so the "
            "test is calibrated."
            if p < 0.05
            else "Phase is not distinguishable from uniform -- no evidence the "
            "detector tracks the furnace rhythm."
        ),
    }
    if args.out:
        Path(args.out).write_text(json.dumps(rep, ensure_ascii=False, indent=2) + "\n", encoding="utf-8")
    print(json.dumps(rep, ensure_ascii=False, indent=2))
    return 0


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