Compressor/backend/indicator_features/watch.py

#!/usr/bin/env python3
"""单个角度箱看板：波动"幅度"与"次数"各两条曲线 + 自适应正常带。

- 波动幅度：① 极差 range/span   ② 粗糙度 roughness/span（去趋势残差）
- 波动次数：① 导数过零数 zero_cross   ② 极值数 n_extrema
每个量用过去 `--window` 天的中位数 ± k×MAD 作正常带；越界点标红。
"""
from __future__ import annotations

import argparse
import csv
import sys
from datetime import datetime, timedelta
from pathlib import Path

import numpy as np

BACKEND = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(BACKEND))

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt  # noqa: E402

matplotlib.rcParams["font.sans-serif"] = ["Arial Unicode MS", "Heiti TC", "Hiragino Sans GB", "DejaVu Sans"]
matplotlib.rcParams["axes.unicode_minus"] = False

from event_bins import EVENTS, load  # noqa: E402


def band_series(times, values, window_days, k, min_base=10):
    center = np.full(len(values), np.nan)
    lo = np.full(len(values), np.nan)
    hi = np.full(len(values), np.nan)
    out = np.zeros(len(values), dtype=bool)
    for i in range(len(values)):
        t0 = times[i] - timedelta(days=window_days)
        idx = [j for j in range(i) if times[j] >= t0 and np.isfinite(values[j])]
        if len(idx) < min_base:
            continue
        hist = values[idx]
        med = float(np.median(hist))
        mad = float(np.median(np.abs(hist - med))) * 1.4826
        if mad <= 1e-12:
            mad = float(np.std(hist))
        center[i], lo[i], hi[i] = med, med - k * mad, med + k * mad
        if np.isfinite(values[i]):
            out[i] = values[i] > hi[i] or values[i] < lo[i]
    return center, lo, hi, out


def plot_panel(ax, times, values, center, lo, hi, out, ylabel, events, plot_start, plot_end):
    ax.fill_between(times, lo, hi, color="0.85", label="正常带(中位±3MAD)")
    ax.plot(times, center, color="0.5", lw=1, label="带中位")
    ax.plot(times, values, "-", color="tab:blue", lw=1, alpha=0.6)
    ax.plot(times, values, "o", ms=3, color="tab:blue")
    if out.any():
        ax.plot(np.asarray(times)[out], np.asarray(values)[out], "o", ms=5, color="red", label="越界")
    for ev in events:
        ax.axvline(datetime.fromisoformat(ev), color="k", ls="--", lw=0.7)
    ax.set_ylabel(ylabel)
    ax.grid(alpha=0.25)
    ax.legend(fontsize=7)
    if plot_start or plot_end:
        mask = np.ones(len(times), dtype=bool)
        if plot_start:
            mask &= np.asarray(times) >= plot_start
        if plot_end:
            mask &= np.asarray(times) < plot_end
        val = np.asarray(values)[mask]
        band = np.r_[np.asarray(lo)[mask], np.asarray(hi)[mask]]
        allv = np.r_[val[np.isfinite(val)], band[np.isfinite(band)]]
        if allv.size:
            pad = 0.1 * (np.nanmax(allv) - np.nanmin(allv) + 1e-9)
            ax.set_ylim(np.nanmin(allv) - pad, np.nanmax(allv) + pad)
        if plot_start and plot_end:
            ax.set_xlim(plot_start, plot_end)


def main() -> None:
    parser = argparse.ArgumentParser(description="单箱 波动幅度/次数 看板")
    parser.add_argument("--csv", required=True)
    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
    parser.add_argument("--cylinder", required=True)
    parser.add_argument("--side", required=True)
    parser.add_argument("--bin-start", type=float, required=True)
    parser.add_argument("--window", type=int, default=30)
    parser.add_argument("--k", type=float, default=3.0)
    parser.add_argument("--plot-start", default="", help="只显示该日期之后（含）；基线仍用全部历史")
    parser.add_argument("--plot-end", default="")
    parser.add_argument("--out", required=True)
    args = parser.parse_args()

    rows = load(args.csv)
    sel = [r for r in rows if r["cylinder"] == args.cylinder and r["side"] == args.side
           and abs(r["bin_start"] - args.bin_start) < 1e-6]
    sel.sort(key=lambda r: r["time"])
    if not sel:
        print("没有匹配数据")
        return
    times = [r["time"] for r in sel]
    span = np.asarray([r["p_span"] or np.nan for r in sel], dtype=float)

    def col(name):
        return np.asarray([r[name] if r[name] is not None else np.nan for r in sel], dtype=float)

    series = [
        (col("range") / span, "波动幅度·极差\nrange/span"),
        (col("roughness") / span, "波动幅度·粗糙度\nroughness/span"),
        (col("zero_cross"), "波动次数·导数过零数"),
        (col("n_extrema"), "波动次数·极值数"),
    ]
    events = [e[1] for e in EVENTS if e[0] == args.unit]
    plot_start = datetime.fromisoformat(args.plot_start) if args.plot_start else None
    plot_end = datetime.fromisoformat(args.plot_end) if args.plot_end else None

    fig, axes = plt.subplots(4, 1, figsize=(13, 10), sharex=True)
    bands = []
    for ax, (values, label) in zip(axes, series):
        c, lo, hi, out = band_series(times, values, args.window, args.k)
        bands.append((values, c, lo, hi, out))
        plot_panel(ax, times, values, c, lo, hi, out, label, events, plot_start, plot_end)
    fig.suptitle(f"{args.unit}号机 {args.cylinder}缸{args.side} 箱{args.bin_start:g}°"
                 f"（基线窗口 {args.window} 天，±{args.k:g}MAD）")
    fig.tight_layout()
    out_path = Path(args.out)
    out_path.parent.mkdir(parents=True, exist_ok=True)
    fig.savefig(out_path, dpi=120)
    plt.close(fig)

    header = ["时间", "极差/span", "粗糙度/span", "导数过零数", "极值数",
              "极差越界", "粗糙度越界", "过零越界", "极值越界"]
    with out_path.with_suffix(".csv").open("w", encoding="utf-8-sig", newline="") as fh:
        wr = csv.writer(fh)
        wr.writerow(header)
        for i, t in enumerate(times):
            wr.writerow([t.isoformat()] + [f"{bands[j][0][i]:.4g}" for j in range(4)]
                        + [int(bands[j][4][i]) for j in range(4)])
    print("图:", out_path)
    print("表:", out_path.with_suffix(".csv"))
    for label, (values, c, lo, hi, out) in zip(
            ["极差", "粗糙度", "过零数", "极值数"], bands):
        print(f"  {label}: 越界 {int(out.sum())} 次")


if __name__ == "__main__":
    main()