goodwe-addon/goodwe_controller/tools/replay_recorded_loop.py

#!/usr/bin/env python3
"""SAFETY-48's evidence: the live install's OWN recorded disturbances, replayed
through the shipped control law with and without the card.

    python3 tools/replay_recorded_loop.py                 # the A/B table
    python3 tools/replay_recorded_loop.py --sensitivity   # ...and does it survive
                                                          #    a different plant?
    python3 tools/replay_recorded_loop.py --export --days 16 --token $HASS_TOKEN

WHY A REPLAY AND NOT A STEP RESPONSE
====================================
`test_control.py`'s SAFETY-39 section drives a family of synthetic export steps
through a first-order plant, and it is a good instrument for the question it
asks. It cannot answer this one. The defect this card is about is a MISSING
CASE - the integrator outrunning the command whenever `|error| > slew_w / gain`
- and whether that case matters is a question about how often the house
actually produces such an error, how long it stands, and what else is happening
when it does. A step family answers "what would happen if"; only the recording
answers "how much did this cost".

WHAT IS REAL IN HERE
====================
  * THE DISTURBANCE. `D = grid + batt` at every recorded meter instant: the
    house's own net demand with THIS loop's battery power taken back out. It is
    exogenous - the house does not know what the battery did - which is what
    makes it valid to drive a counterfactual with, for a run of any length.
  * THE TICK. The replay steps on the METER's own recorded sample instants, not
    on a constant. Median 5.0 s over the window. ⚠️ `control.LOOP_TICK_S` says
    2.0 s; that was measured on 2026-09-05 and is NOT what this install has been
    running since. The figure matters to `SLEW_RAIL_FRAC`, and it is wrong in
    the SAFE direction (a slower meter makes validate.py rule 10 conservative).
  * THE CONTROL LAW. `app/control.py`, imported out of this checkout. The
    "before" law is THE SAME FILE with the SAFETY-48 block removed by a named
    substitution and exec'd into a fresh namespace - so it is release/1.0's
    behaviour derived from the shipped bytes, never a transcription. The
    substitution anchor must occur exactly once or this refuses to run.

MODELLED, AND THEREFORE NOT EVIDENCE ABOUT THE THING MODELLED
=============================================================
  * THE INVERTER. A pure delay + first-order lag + affine gain, FITTED on this
    same recording at SoC above the floor. Its residual is printed on every run
    and is large (~375 W RMS): the replay tracks the meter WORSE than the real
    install does, so the absolute export figures below are pessimistic by about
    a third against the recording. ⚠️ THE DELTA IS WHAT THIS TOOL IS FOR - both
    laws run against the same plant, on the same ticks, from the same seed.
  * THE PACK. SoC integrated from the simulated battery power off the recorded
    SoC, discharge refused at or below the floor, charge refused at 100 %.
  * EVERYTHING ELSE THE ARBITER DOES. The recorded setpoint is the ARBITER's
    output - the loop's claim after the peak guard, the SoC clamps, maintenance
    and arbitrage have had their say. This replay runs `compute()` ALONE. That
    is the largest single reason the fidelity section below is not tighter, and
    it is why no absolute figure here should be quoted as "what the install
    did"; the recording itself is that.

⚠️ WHAT THIS CANNOT SHOW, SAID BEFORE THE RESULTS RATHER THAN AFTER THEM
=======================================================================
  * IT CANNOT SHOW THAT THE FIX IS SAFE ON THE REAL PLANT. Nothing here drives
    an inverter. The rig and a live run are what answer that, and SAFETY-47
    records why a shadow run cannot.
  * IT CANNOT SEPARATE THE FIX FROM THE PLANT MODEL ON ONE RUN. `--sensitivity`
    re-runs the whole table on five other plants and tunings; the ranking is the
    claim, not the watt figures.
  * THE BILLED QUANTITY IS A MAXIMUM OVER ~1500 QUARTERS, so it is a tail
    statistic and one quarter decides it. The per-quarter MEAN delta is printed
    beside it for exactly that reason.
"""
import argparse
import bisect
import gzip
import hashlib
import json
import math
import os
import pathlib
import sys
import time
import urllib.error
import urllib.parse
import urllib.request
from datetime import datetime, timedelta

HERE = pathlib.Path(__file__).resolve().parent
ADDON = HERE.parent
sys.path.insert(0, str(ADDON))

SCENARIO = ADDON / "scenarios" / "loop-2026-09-02_2026-09-17"
DEFAULT_URL = "http://192.168.2.158:8123"
TZNAME = "Europe/Brussels"

# The entities, and what each is for. ⚠️ `setpoint` is NOT an input to the
# replay: it is the RECORDING the replay is graded against in §fidelity, and it
# is what `slew_w` and `max_w` are read off rather than assumed.
ENTITIES = {
    "grid": "sensor.p1_meter_active_power",
    "batt": "sensor.goodwe_rs485_controller_goodwe_battery_power",
    "setpoint": "sensor.goodwe_rs485_controller_goodwe_setpoint",
    "soc": "sensor.goodwe_rs485_controller_goodwe_battery_soc",
}

# The live install's shape. `max_w` and `slew_w` are READ OFF THE RECORDING by
# `--export` and re-derived on every run below; `gain` and `target_grid_w` are
# the options and are not observable in the trace.
LIVE = dict(gain=0.6, max_w=5000.0, slew_w=400.0, target_grid_w=-10.0)
SOC_FLOOR = 10.0
CAP_WH = 17379.0        # measured: -dE/dSoC over 1604 segments of this record
GRID_HOLD_S = 30.0      # `meter_max_age_s`: past this a held meter is no reading
TICK_CAP_S = 60.0       # a longer gap is an outage, not a cycle


# ---------------------------------------------------------------------------
# the recording


def _rows(key):
    out = []
    for fn in sorted(os.listdir(SCENARIO)):
        if not fn.endswith(f"-{key}.jsonl.gz"):
            continue
        with gzip.open(SCENARIO / fn, "rb") as fh:
            for line in fh:
                r = json.loads(line)
                out.append((datetime.fromisoformat(r["t"]).timestamp(), r["v"]))
    out.sort(key=lambda r: r[0])
    ded = []
    for t, v in out:
        if ded and abs(ded[-1][0] - t) < 1e-9:
            ded[-1] = (t, v)          # the history endpoint's day-boundary repeat
            continue
        ded.append((t, v))
    return ded


def _steps(key):
    """[(t, float|None)] - None OPENS A HOLE and nothing is held across it.

    ⚠️ A LONG SILENCE IS NOT A HOLE. Home Assistant's recorder stores a row only
    when a state CHANGES, so eleven hours of `setpoint` silence means the
    command held that value for eleven hours - which is exactly what this site
    does overnight at its SoC floor. What IS a hole is a row whose state is
    `unavailable`/`unknown`, and `grid` additionally stops being believed after
    GRID_HOLD_S because that is the age the shipped controller stops believing
    it at.
    """
    out = []
    for t, s in _rows(key):
        try:
            out.append((t, float(s)))
        except (TypeError, ValueError):
            out.append((t, None))
    return out


def stream():
    """[(t, grid, batt, setpoint, soc)] at every recorded METER instant."""
    g = [(t, v) for t, v in _steps("grid") if v is not None]
    others = {}
    for k in ("batt", "setpoint", "soc"):
        s = _steps(k)
        others[k] = ([r[0] for r in s], [r[1] for r in s])
    prev_t = None
    out = []
    for t, v in g:
        if prev_t is not None and t - prev_t > GRID_HOLD_S:
            prev_t = t
            continue
        prev_t = t
        row = [t, v]
        for k in ("batt", "setpoint", "soc"):
            tt, vv = others[k]
            j = bisect.bisect_right(tt, t) - 1
            row.append(None if j < 0 else vv[j])
        if any(x is None for x in row):
            continue
        out.append(tuple(row))
    return out


# ---------------------------------------------------------------------------
# the plant, fitted on this same recording


def _lsq(xs, ys):
    n = len(xs)
    mx, my = sum(xs) / n, sum(ys) / n
    sxy = sum((a - mx) * (b - my) for a, b in zip(xs, ys))
    sxx = sum((a - mx) ** 2 for a in xs)
    syy = sum((b - my) ** 2 for b in ys)
    return sxy / sxx, my - (sxy / sxx) * mx, sxy / math.sqrt(sxx * syy)


def fit_plant(rows, stride=1):
    """(lag_s, tau_s, gain, offset, rms_W), fitted where the pack can MOVE.

    ⚠️ THE SoC FILTER IS THE WHOLE FIT. This install spends 36 % of its ticks at
    or below a 10 % SoC with the command standing near 2 kW and the battery
    delivering nothing, and a fit taken over all ticks reads
    `batt = 0.50 x cmd - 396 W, r = 0.70` - a plant that follows half the
    command, which no inverter does. Above the floor the same data says
    `batt = 1.00 x cmd - 28 W, r = 0.977`. The floor is the PACK and is modelled
    separately.
    """
    t = [r[0] for r in rows]
    ba = [r[2] for r in rows]
    sp = [r[3] for r in rows]
    live = [r[4] > SOC_FLOOR + 2 for r in rows]
    best = None
    for lag in (0.0, 2.5, 5.0, 7.5, 10.0, 12.5):
        j, u = 0, []
        for i, ti in enumerate(t):
            while j + 1 < len(t) and t[j + 1] <= ti - lag:
                j += 1
            u.append(sp[j])
        idx = [i for i in range(0, len(t), stride)
               if live[i] and abs(u[i]) > 50]
        g, o, _ = _lsq([u[i] for i in idx], [ba[i] for i in idx])
        for tau in (0.5, 1.0, 2.0, 4.0, 8.0, 16.0):
            b, err, cnt = 0.0, 0.0, 0
            for i in range(1, len(t)):
                dt = min(t[i] - t[i - 1], TICK_CAP_S)
                b += (1.0 - math.exp(-dt / tau)) * (g * u[i] + o - b)
                if live[i]:
                    err += (b - ba[i]) ** 2
                    cnt += 1
            rms = math.sqrt(err / cnt)
            if best is None or rms < best[-1]:
                best = (lag, tau, g, o, rms)
    return best


class Plant:
    def __init__(self, lag_s, tau_s, g, o, soc0):
        self.lag_s, self.tau_s, self.g, self.o = lag_s, tau_s, g, o
        self.b, self.soc, self.q = 0.0, soc0, []

    def seed(self, t, b, cmd):
        self.b, self.q = b, [(t - self.lag_s - 1e-6, cmd)]

    def step(self, t, dt, cmd):
        self.q.append((t, cmd))
        want = self.q[0][1]
        while len(self.q) > 1 and self.q[1][0] <= t - self.lag_s:
            self.q.pop(0)
            want = self.q[0][1]
        tgt = self.g * want + self.o
        if self.soc <= SOC_FLOOR:
            tgt = min(tgt, 0.0)          # an empty pack cannot discharge
        if self.soc >= 100.0:
            tgt = max(tgt, 0.0)          # a full one cannot charge
        self.b += (1.0 - math.exp(-dt / self.tau_s)) * (tgt - self.b)
        self.soc = max(0.0, min(100.0, self.soc - self.b * dt / 36.0 / CAP_WH))
        return self.b


# ---------------------------------------------------------------------------
# the two laws


SAFETY_48_BLOCK_ANCHOR = "        _lead_cap = SLEW_LEAD_STEPS * slew\n"


def laws():
    """(before, after) - the shipped module, and the same SOURCE without the card.

    ⚠️ NOT A TRANSCRIPTION, AND NOT A FLAG. `test_mqtt.py` part 8's technique:
    the shipped file's bytes with one named substitution, exec'd into a fresh
    namespace. An anchor that does not occur exactly once is a refusal, because
    a substitution that applied to nothing would be reported as "the card
    changes nothing" - MUTATION.md's survivor, wearing an instrument's clothes.
    """
    import app.control as after
    src = (ADDON / "app" / "control.py").read_text()
    if src.count(SAFETY_48_BLOCK_ANCHOR) != 1:
        raise SystemExit(
            "replay_recorded_loop: app/control.py does not contain exactly one\n"
            f"  {SAFETY_48_BLOCK_ANCHOR.strip()!r}\n"
            "NOTHING HAS BEEN MEASURED. This is a REFUSAL, not a result: with no\n"
            "anchor there is no `before` law to compare against, and a run that\n"
            "quietly compared the shipped law with ITSELF would report the card\n"
            "as changing nothing - which is the survivor MUTATION.md is about,\n"
            "wearing an instrument's clothes.")
    start = src.index(SAFETY_48_BLOCK_ANCHOR)
    end = src.index("    target = slewed\n", start)
    block = src[start:end]
    stripped = src.replace(block, "")
    if "SLEW_LEAD_STEPS * slew" in stripped:
        raise SystemExit("replay_recorded_loop: the block was not removed - "
                         "NOTHING HAS BEEN MEASURED")
    before = type(sys)("control_before")
    exec(compile(stripped, "<release/1.0 law>", "exec"), before.__dict__)
    return before, after


def replay(rows, law, fit, tuning):
    """[(t, dt, grid, batt, cmd, i_w, lead, reason, soc)] over the whole record."""
    t0, g0, b0, sp0, soc0 = rows[0]
    pl = Plant(*fit[:4], soc0)
    pl.seed(t0, b0, sp0)
    prev, i_w, sat = sp0, sp0, 0
    out = []
    for i in range(1, len(rows)):
        t, gr, ba, sp, soc = rows[i]
        dt = t - rows[i - 1][0]
        if dt > TICK_CAP_S:
            # An outage. Re-seed from the recording rather than inventing a
            # cycle: nothing this loop did across it is knowable.
            pl = Plant(*fit[:4], soc)
            pl.seed(t, ba, sp)
            prev, i_w, sat = sp, sp, 0
            continue
        D = gr + ba                       # the exogenous half of the meter
        b = pl.step(t, dt, prev)
        grid = D - b
        d = law.compute(prev, grid, b, tuning, sat, i_w)
        out.append((t, dt, grid, b, d.target_w, d.i_w, d.i_w - prev, d.reason,
                    pl.soc))
        prev, i_w, sat = d.target_w, d.i_w, d.sat_count
    return out


# ---------------------------------------------------------------------------
# the five reported quantities - ONE implementation, run on both traces


def export_kwh(tr):
    return sum(max(0.0, -r[2]) * r[1] for r in tr) / 3.6e6


def import_kwh(tr):
    return sum(max(0.0, r[2]) * r[1] for r in tr) / 3.6e6


def reversals(tr, dead_w=50.0):
    n, last = 0, 0
    for r in tr:
        s = 1 if r[3] > dead_w else (-1 if r[3] < -dead_w else 0)
        if s and last and s != last:
            n += 1
        if s:
            last = s
    return n


def wasted_kwh(tr, bin_s=60.0):
    """Throughput that went round rather than anywhere, in FIXED clock bins.

    Per bin, `sum |b|dt - |sum b dt|`. A bin spent charging, or discharging, or
    idle contributes ZERO; only a bin that does both contributes. ⚠️ The bins
    are wall-clock and identical for both laws on purpose: an episode-scoped
    definition would MOVE with the fix and compare two different questions.
    """
    tot, cur, absE, net = 0.0, None, 0.0, 0.0
    for r in tr:
        k = math.floor(r[0] / bin_s)
        if cur is None:
            cur = k
        if k != cur:
            tot += abs(absE) - abs(net)
            cur, absE, net = k, 0.0, 0.0
        absE += abs(r[3]) * r[1]
        net += r[3] * r[1]
    return (tot + abs(absE) - abs(net)) / 3.6e6


def quarters(tr):
    """{quarter -> mean OFFTAKE in W}. The billed quantity.

    Calendar quarter-hours, the average of NET power, floored at 0 - a quarter
    that nets to export is not offtake and is not billed. A quarter with less
    than 80 % coverage is not a quarter and is dropped.
    """
    acc = {}
    for r in tr:
        k = math.floor(r[0] / 900.0)
        e, s = acc.get(k, (0.0, 0.0))
        acc[k] = (e + r[2] * r[1], s + r[1])
    return {k: max(0.0, e / s) for k, (e, s) in acc.items() if s > 720.0}


def recoveries(tr, bias, onset_w, cap_s=900.0):
    """Seconds from the error leaving `onset_w` to it back inside the deadband.

    ⚠️ AN EXCURSION THE PACK CANNOT ANSWER IS NOT A RECOVERY TIME. This install
    spends 36 % of its ticks at the SoC floor with the command standing near
    2 kW and the battery delivering nothing; the error there is not the loop's
    to clear and every law scores the same non-answer. Leaving those in gave a
    p99 of eight hours. An excursion counts only if the pack had headroom in the
    needed direction for its whole length, and one still standing after `cap_s`
    is counted as UNRECOVERED rather than as a large number.
    """
    out, unrec, t_on, ok = [], 0, None, True
    for r in tr:
        err = r[2] - bias
        head = (r[8] > SOC_FLOOR + 1) if err > 0 else (r[8] < 99.0)
        if t_on is None:
            if abs(err) > onset_w:
                t_on, ok = r[0], head
        else:
            ok = ok and head
            if abs(err) < 15.0:
                if ok:
                    out.append(r[0] - t_on)
                t_on = None
            elif r[0] - t_on > cap_s:
                unrec += 1 if ok else 0
                t_on = None
    return sorted(out), unrec


def pct(xs, p):
    return xs[min(len(xs) - 1, int(round(p / 100.0 * (len(xs) - 1))))] if xs else float("nan")


def table(tr, tuning):
    rec, unrec = recoveries(tr, tuning.target_grid_w, tuning.slew_w / tuning.gain)
    leads = sorted(abs(r[6]) for r in tr)
    q = quarters(tr)
    return dict(export_kwh=export_kwh(tr), import_kwh=import_kwh(tr),
                reversals=reversals(tr), wasted_kwh=wasted_kwh(tr),
                worst_quarter_w=max(q.values()), rec_n=len(rec),
                rec_med=pct(rec, 50), rec_p99=pct(rec, 99),
                rec_max=rec[-1] if rec else float("nan"),
                rec_unrecovered=unrec, lead_p99=pct(leads, 99),
                lead_max=leads[-1], slew_ticks=sum(1 for r in tr if r[7] == "slew-limited"),
                ticks=len(tr), _q=q)


ROWS = [("export_kwh", "export kWh", "the cost"),
        ("import_kwh", "import kWh", ""),
        ("reversals", "battery sign reversals", ""),
        ("wasted_kwh", "wasted throughput kWh", ""),
        ("worst_quarter_w", "worst 15-min offtake W", "⚠️ THE BILLED ONE"),
        ("rec_n", "recovery events", ""),
        ("rec_med", "recovery median s", ""),
        ("rec_p99", "recovery p99 s", ""),
        ("rec_max", "recovery max s", ""),
        ("rec_unrecovered", "unrecovered (>900 s)", ""),
        ("lead_p99", "|i_w - prev_w| p99 W", ""),
        ("lead_max", "|i_w - prev_w| max W", ""),
        ("slew_ticks", "slew-limited ticks", "")]


def show(a, b, name_a="release/1.0", name_b="SAFETY-48"):
    print(f"    {'':26} {name_a:>13} {name_b:>13} {'delta':>11}")
    for k, lab, note in ROWS:
        print(f"    {lab:<26} {a[k]:13.2f} {b[k]:13.2f} {b[k] - a[k]:+11.2f}"
              + (f"   {note}" if note else ""))
    d = [b["_q"][k] - a["_q"][k] for k in a["_q"] if k in b["_q"]]
    print(f"    {'mean quarter delta W':<26} {'':13} {'':13} "
          f"{sum(d) / len(d):+11.2f}   over {len(d)} quarters")
    print(f"    {'worst quarter delta W':<26} {'':13} {'':13} {max(d):+11.2f}")


def fidelity(tr, rows):
    """What the BEFORE run gets wrong against the recording it came from.

    ⚠️ PRINTED BEFORE THE A/B TABLE IS READ, NOT AFTER. This is the number that
    says how far any ABSOLUTE figure above may be trusted, and it is large.
    """
    n = len(tr)
    rec_g = [rows[i + 1][1] for i in range(n)]
    rec_sp = [rows[i + 1][3] for i in range(n)]
    sim_g = [r[2] for r in tr]
    sim_c = [r[4] for r in tr]
    dt = [r[1] for r in tr]
    _, _, r_sp = _lsq(sim_c, rec_sp)
    _, _, r_g = _lsq(sim_g, rec_g)
    ex_sim = sum(max(0.0, -g) * d for g, d in zip(sim_g, dt)) / 3.6e6
    ex_rec = sum(max(0.0, -g) * d for g, d in zip(rec_g, dt)) / 3.6e6
    print(f"    setpoint  r={r_sp:.3f}  mean {sum(sim_c) / n:+8.0f} W simulated "
          f"vs {sum(rec_sp) / n:+8.0f} W recorded")
    print(f"    grid      r={r_g:.3f}  mean {sum(sim_g) / n:+8.0f} W simulated "
          f"vs {sum(rec_g) / n:+8.0f} W recorded")
    print(f"    export    {ex_sim:8.2f} kWh simulated vs {ex_rec:8.2f} kWh "
          f"recorded  ->  the model is {100 * (ex_sim / ex_rec - 1):+.0f} % "
          f"PESSIMISTIC about the loop's tracking")


# ---------------------------------------------------------------------------
# --export: the read-only capture


def _get(url, token, timeout=180.0):
    req = urllib.request.Request(url, headers={"Authorization": f"Bearer {token}"})
    last = None
    for attempt in range(1, 5):
        try:
            with urllib.request.urlopen(req, timeout=timeout) as resp:
                if resp.status != 200:
                    raise RuntimeError(f"HTTP {resp.status}")
                return json.loads(resp.read().decode()), attempt
        except Exception as exc:
            last = f"{type(exc).__name__}: {exc}"
            if attempt < 4:
                time.sleep(3 * attempt)
    raise RuntimeError(last)


def export(base, token, days, out_dir):
    """One local day per request, per entity. Read-only: GET only, never a write.

    ⚠️ URL-ENCODE THE TIMESTAMPS. The history endpoint answers `Invalid
    end_time` otherwise, and it answers it with a 400 that reads like a bad
    entity id.
    """
    from zoneinfo import ZoneInfo
    tz = ZoneInfo(TZNAME)
    out_dir = pathlib.Path(out_dir)
    out_dir.mkdir(parents=True, exist_ok=True)
    today = datetime.now(tz).date()
    man = {"generated": datetime.now(tz).isoformat(),
           "generator": "tools/replay_recorded_loop.py --export",
           "source": f"{base}/api/history/period", "timezone": TZNAME,
           "entities": ENTITIES,
           "format": '{"t": <ISO-8601 with offset>, "v": <state string>} per '
                     'line, gzip, mtime 0',
           "days": []}
    rc = 0
    for d in range(days, 0, -1):
        day = today - timedelta(days=d)
        s = datetime(day.year, day.month, day.day, tzinfo=tz)
        e = s + timedelta(days=1)
        rec = {"date": str(day), "window": [s.isoformat(), e.isoformat()],
               "files": {}}
        for key, ent in ENTITIES.items():
            q = urllib.parse.urlencode({"filter_entity_id": ent,
                                        "end_time": e.isoformat(),
                                        "minimal_response": ""})
            url = (f"{base}/api/history/period/"
                   f"{urllib.parse.quote(s.isoformat())}?{q}")
            try:
                doc, attempts = _get(url, token)
            except Exception as exc:
                rec["files"][key] = {"read": "FAILED", "error": str(exc)}
                print(f"  {day} {key:9s} FAILED {exc}", flush=True)
                rc = 1
                continue
            series = doc[0] if doc else []
            body = b""
            n_bad = 0
            for row in series:
                t = row.get("last_changed") or row.get("last_updated")
                v = row.get("state")
                if t is None or v is None:
                    continue
                try:
                    float(v)
                except (TypeError, ValueError):
                    n_bad += 1
                body += (json.dumps({"t": t, "v": v}) + "\n").encode()
            fn = f"{day}-{key}.jsonl.gz"
            with gzip.GzipFile(out_dir / fn, "wb", mtime=0) as fh:
                fh.write(body)
            rec["files"][key] = {
                "read": "OK", "entity": ent, "samples": len(series),
                "non_numeric_states": n_bad, "attempts": attempts, "file": fn,
                "bytes_raw": len(body),
                "sha256_raw": hashlib.sha256(body).hexdigest()}
            print(f"  {day} {key:9s} {len(series):7d} rows, {n_bad} non-numeric",
                  flush=True)
        man["days"].append(rec)
    (out_dir / "manifest.json").write_text(json.dumps(man, indent=1) + "\n")
    return rc


# ---------------------------------------------------------------------------


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--export", action="store_true")
    ap.add_argument("--days", type=int, default=16)
    ap.add_argument("--url", default=DEFAULT_URL)
    ap.add_argument("--token", default=os.environ.get("HASS_TOKEN"))
    ap.add_argument("--out-dir", default=str(SCENARIO))
    ap.add_argument("--sensitivity", action="store_true")
    args = ap.parse_args()

    if args.export:
        if not args.token:
            raise SystemExit("--export needs --token or $HASS_TOKEN")
        return export(args.url, args.token, args.days, args.out_dir)

    from app.control import Tuning
    t0 = time.time()
    rows = stream()
    span = (rows[-1][0] - rows[0][0]) / 86400.0
    dts = sorted(rows[i + 1][0] - rows[i][0] for i in range(len(rows) - 1))
    print(f"\nTHE RECORDING: {len(rows)} meter instants over {span:.2f} days, "
          f"{SCENARIO.name}")
    print(f"    tick   p5 {pct(dts, 5):.2f} s   median {pct(dts, 50):.2f} s   "
          f"p95 {pct(dts, 95):.2f} s   (control.LOOP_TICK_S says 2.0)")

    # ⚠️ READ OFF THE RECORDING, NOT ASSUMED. `max_w` is the setpoint's own
    # extreme and `slew_w` is the step the recorded command takes over and over.
    sp = [r[3] for r in rows]
    steps = {}
    for i in range(1, len(rows)):
        if rows[i][0] - rows[i - 1][0] < 8.0:
            steps[round(abs(sp[i] - sp[i - 1]))] = \
                steps.get(round(abs(sp[i] - sp[i - 1])), 0) + 1
    top = max((c, w) for w, c in steps.items() if w > 100)
    print(f"    the recorded command: extremes {min(sp):+.0f} / {max(sp):+.0f} W"
          f"  ->  max_w {LIVE['max_w']:.0f};  most common step above 100 W is "
          f"{top[1]:.0f} W, {top[0]} times  ->  slew_w {LIVE['slew_w']:.0f}")

    fit = fit_plant(rows, stride=7)
    print(f"    the fitted inverter: lag {fit[0]:.1f} s, tau {fit[1]:.1f} s, "
          f"batt = {fit[2]:.4f} x cmd {fit[3]:+.1f} W, residual {fit[4]:.0f} W "
          f"RMS above the SoC floor")

    before, after = laws()
    T = Tuning(**LIVE)
    tr_b = replay(rows, before, fit, T)
    tr_a = replay(rows, after, fit, T)
    tab_b, tab_a = table(tr_b, T), table(tr_a, T)

    print("\nFIDELITY - the BEFORE law replayed against the recording it came from")
    fidelity(tr_b, rows)

    print("\nTHE A/B, same ticks, same disturbance, same plant, same seed")
    show(tab_b, tab_a)

    if args.sensitivity:
        print("\nSENSITIVITY - the ranking, not the watts")
        variants = [
            ("plant lag 5.0 s", (5.0, fit[1], fit[2], fit[3], 0.0), LIVE),
            ("plant lag 10.0 s", (10.0, fit[1], fit[2], fit[3], 0.0), LIVE),
            ("plant tau 8 s", (fit[0], 8.0, fit[2], fit[3], 0.0), LIVE),
            ("gain 1.0135, offset -26", (fit[0], fit[1], 1.0135, -26.0, 0.0), LIVE),
            ("slew_w 555.6 (the rail)", fit, dict(LIVE, slew_w=555.6)),
            ("gain 0.3", fit, dict(LIVE, gain=0.3)),
        ]
        print(f"      {'variant':<26}{'export kWh':>22}{'worst quarter W':>26}"
              f"{'mean quarter dW':>18}")
        for name, f2, tune in variants:
            T2 = Tuning(**tune)
            b = table(replay(rows, before, f2, T2), T2)
            a = table(replay(rows, after, f2, T2), T2)
            d = [a["_q"][k] - b["_q"][k] for k in b["_q"] if k in a["_q"]]
            print(f"      {name:<26}{b['export_kwh']:9.1f} ->{a['export_kwh']:9.1f}"
                  f"{b['worst_quarter_w']:13.0f} ->{a['worst_quarter_w']:9.0f}"
                  f"{sum(d) / len(d):+18.2f}")

    print(f"\n    ({time.time() - t0:.0f} s)")
    return 0


if __name__ == "__main__":
    sys.exit(main())