#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""事后分析 rec_ecat_slots.py 落盘的 ECAT 原始流：相邻数据大跳变检测。

思路（按用户指定）：不看协议判据，只看数值形态。同一台电机的 STATUS 是连续量，
相邻两帧（按 seq 排序）之间出现"跳出去又立刻跳回来"的孤立尖峰，就是这一帧的字节
被拼坏了 —— 真突变会保持住（下一帧仍在新值附近），不会一帧就回头。

指标：每个 (node,字段) 先取全体相邻差的稳健尺度 scale = median(|Δ|)，
  大跳变   : |v[i]-v[i-1]| > K*scale
  孤立尖峰 : 大跳变 且 |v[i+1]-v[i-1]| <= 2*K*scale   ← 撕裂签名
  跳变保持 : 大跳变 但下一帧不回头的（真实突变/新参考值），只计数不算异常

--selftest 会把"前 j 字节取新帧、其余取上一帧"的混合帧注入实测数据，
逐 j 报告检出率：证明这套判据对撕裂本身有灵敏度（否则 0 报警没有意义）。

用法:
  python scripts/anal_ecat_jumps.py ecat_slots_raw.csv [--k 8]
  python scripts/anal_ecat_jumps.py ecat_slots_raw.csv --selftest
"""
import csv
import statistics
import struct
import sys

# (名称, 偏移, 解码器)
FIELDS = [
    ("mcuTemp", 0, lambda b: struct.unpack_from("<h", b, 0)[0] * 0.01),
    ("motTemp", 2, lambda b: struct.unpack_from("<h", b, 2)[0] * 0.01),
    ("uCur", 4, lambda b: struct.unpack_from("<h", b, 4)[0]),
    ("vCur", 6, lambda b: struct.unpack_from("<h", b, 6)[0]),
    ("wCur", 8, lambda b: struct.unpack_from("<h", b, 8)[0]),
    ("flags", 10, lambda b: b[10]),
    ("pwr_V", 11, lambda b: struct.unpack_from("<H", b, 11)[0] * 0.001),
    ("rtAngle", 13, lambda b: struct.unpack_from("<f", b, 13)[0]),
    ("rtVel", 17, lambda b: struct.unpack_from("<f", b, 17)[0]),
    ("rtCur", 21, lambda b: struct.unpack_from("<f", b, 21)[0]),
    ("testF", 25, lambda b: struct.unpack_from("<f", b, 25)[0]),
    ("testU16", 29, lambda b: struct.unpack_from("<H", b, 29)[0]),
    ("cpu%", 31, lambda b: b[31]),
    ("initSt", 32, lambda b: struct.unpack_from("<I", b, 32)[0]),
]


def load(path):
    """返回 {node: [(seq, 36B bytes), ...]} 按 seq 升序。"""
    per = {}
    with open(path, "r", encoding="ascii", errors="ignore") as f:
        for row in csv.reader(f):
            if not row or row[0] != "M":
                continue
            try:
                node, ch, seq = int(row[2]), int(row[3]), int(row[4])
                payload = bytes.fromhex(row[7])
            except (ValueError, IndexError):
                continue
            if node <= 0 or len(payload) < 36:
                continue
            per.setdefault(node, []).append((seq, payload))
    for v in per.values():
        v.sort(key=lambda x: x[0])
    return per


def decode(payload):
    out = {}
    for name, off, fn in FIELDS:
        try:
            out[name] = float(fn(payload))
        except Exception:
            out[name] = None
    return out


def analyze(per, k=8.0, verbose=True):
    """逐 (node,字段) 找孤立尖峰。返回总尖峰数与明细列表。"""
    spikes, holds, details = 0, 0, []
    tails = []
    for node, rows in sorted(per.items()):
        if len(rows) < 40:
            continue
        cols = [decode(p) for _, p in rows]
        for name, _off, _fn in FIELDS:
            v = [c[name] for c in cols]
            finite = [x for x in v if x is not None and x == x and abs(x) != float("inf")]
            if len(finite) < 40:
                continue
            deltas = [abs(v[i] - v[i - 1]) for i in range(1, len(v))
                      if v[i] is not None and v[i - 1] is not None
                      and v[i] == v[i] and v[i - 1] == v[i - 1]]
            if not deltas:
                continue
            # 尺度取 |Δ| 的 p99 分位数（不是 median）：量化字段（mcuTemp 这类）
            # 的 median(|Δ|)=0，用 median 会让阈值塌到 eps，正常的 1LSB 抖动
            # 全被判成尖峰。p99 代表"全体相邻差的尾部"，阈值 = K×p99 才是
            # "远超正常抖动尾巴"的意思。
            ordered = sorted(deltas)
            scale = max(ordered[min(len(ordered) - 1, int(len(ordered) * 0.99))], 1e-9)
            mx, mxctx = 0.0, None
            for i in range(1, len(v) - 1):
                a, b, c = v[i - 1], v[i], v[i + 1]
                if None in (a, b, c) or any(x != x for x in (a, b, c)):
                    continue
                d = abs(b - a)
                if d > mx:
                    mx, mxctx = d, (a, b, c)
                if d <= k * scale:
                    continue
                back = abs(c - a) <= 2 * k * scale
                if back:
                    spikes += 1
                    if len(details) < 25:
                        details.append((node, name, rows[i][0], scale, a, b, c,
                                        rows[i][1].hex()))
                else:
                    holds += 1
            if mxctx and mx > 3 * scale and len(tails) < 60:
                tails.append((node, name, scale, mx, mxctx))
    if verbose:
        print(f"孤立尖峰(撕裂签名) = {spikes}    跳变保持(真实突变) = {holds}")
        for node, name, seq, scale, a, b, c, hexs in details[:12]:
            print(f"  node{node} {name}: seq={seq} 尺度={scale:.4g} "
                  f"{a:.6g} -> {b:.6g} -> {c:.6g}")
            print(f"      payload={hexs}")
        if tails:
            tails.sort(key=lambda x: -(x[3] / max(x[2], 1e-12)))
            print("各字段最大相邻差（倍数 = 最大|Δ| / p99 尺度）：")
            for node, name, scale, mx, (a, b, c) in tails[:14]:
                print(f"  node{node} {name:9s} p99尺度={scale:.4g}  最大|Δ|={mx:.4g} "
                      f"({mx/max(scale,1e-12):.0f}×)  {a:.6g} -> {b:.6g} -> {c:.6g}")
    return spikes, holds, details


def selftest(per, k=8.0, n_inject=40):
    """注入混合帧，测各边界 j 的检出率（证明判据有灵敏度）。"""
    import random
    random.seed(7)
    base_sp, _bh, _bd = analyze(per, k=k, verbose=False)
    nodes = [n for n, r in per.items() if len(r) >= 60]
    if not nodes:
        print("自检跳过：数据太短（每节点需 >=60 帧）")
        return
    print(f"自检：基线孤立尖峰 = {base_sp}；向实测流注入『新帧前缀 j 字节 + 上一帧后缀』"
          f"混合帧，每节点各 {n_inject} 次")
    print("  j(字节)   净增尖峰/注入总数   检出率")
    for j in (1, 4, 8, 12, 13, 16, 17, 20, 21, 24, 25, 28, 31, 32, 35):
        rnd = random.Random(j)
        merged = {}
        for node in nodes:
            rows = per[node]
            idxs = rnd.sample(range(2, len(rows) - 2), min(n_inject, len(rows) - 6))
            new = list(rows)
            for i in idxs:
                cur, prev = rows[i][1], rows[i - 1][1]
                new[i] = (rows[i][0], cur[:j] + prev[j:36])
            merged[node] = new
        sp, _h, det = analyze(merged, k=k, verbose=False)
        inj_n = n_inject * len(nodes)
        net = sp - base_sp
        worst = max((abs(d[5] - d[4]) for d in det), default=0.0)
        print(f"  j={j:<4d}  {net:5d}/{inj_n:<5d}          {100.0*net/inj_n:5.1f}%   "
              f"最大偏离≈{worst:.4g}")
    print("（j 越大越靠尾部字节，混合后数值偏离越小、越难检出 —— 这是形态判据的物理上限）")


def main():
    args = [a for a in sys.argv[1:] if not a.startswith("--")]
    if not args:
        print(__doc__)
        return
    path = args[0]
    k = 8.0
    if "--k" in sys.argv:
        k = float(sys.argv[sys.argv.index("--k") + 1])
    per = load(path)
    tot = sum(len(v) for v in per.values())
    print(f"读入 {path}: {tot} 帧，节点 {sorted(per)}")
    for node, rows in sorted(per.items()):
        print(f"  node{node}: {len(rows)} 个不同 seq，"
              f"seq 范围 {rows[0][0]}..{rows[-1][0]}")
    if "--selftest" in sys.argv:
        selftest(per, k=k)
        return
    analyze(per, k=k)


if __name__ == "__main__":
    main()