Compressor/backend/indicator_features/event_bins.py

#!/usr/bin/env python3
"""Same-point before/after comparison of the per-angular-bin and whole-cycle features.

Reads the long CSV produced by ``angle_bins.py`` and, for each valve/piston repair
event, compares the **same measurement point** between the running period before
the event and the running period after it.  Only same-point comparisons are made
(never across cylinders).

Two feature families (see ``bins_common.py``):

* per-bin features  : level/shape split into whole-point shift (all bins move =
  load / sensor) and local deviation (only some bins move = valve / local fault);
* cycle features    : suction/discharge plateaus and valve-event angles/slopes.

Running periods come from ``run_segments.py`` (``wave_file.rpm > 0``); the pre
window is carried back over earlier running segments until it has enough days
(short post windows are used as-is).  Output is a ranked list per event plus a
detail CSV.  Read-only.
"""
from __future__ import annotations

import argparse
import csv
from collections import defaultdict
from datetime import datetime
from pathlib import Path

import numpy as np

from bins_common import (
    ANGLE_METRICS, BIN_METRICS, CYCLE_METRICS, CYCLE_VOLATILE, VOLATILE, apply_normalization,
)
from run_segments import build_segments, unit_segments


def circ_median(values):
    radians = np.radians(np.asarray(values, dtype=float))
    return float(np.degrees(np.arctan2(np.mean(np.sin(radians)), np.mean(np.cos(radians)))) % 360.0)


def circ_diff(a, b):
    return ((b - a + 180.0) % 360.0) - 180.0

Z = 3.0
FLOOR_FRAC = 0.2  # per-bin scale floor = 20% of that metric's median bin scale
MIN_SPAN = 0.3  # a dead/stuck channel has a near-zero per-cycle pressure span
ANGLE_SCALE_FLOOR = 1.0  # degrees; a constant valve-angle must not inflate z

# 故障汇总.xlsx 中的气阀/活塞事件（unit, date, kind）
EVENTS = [
    ("7", "2025-01-10", "活塞"),
    ("7", "2025-08-08", "气阀"),
    ("7", "2025-10-07", "气阀"),
    ("7", "2026-05-21", "气阀"),
    ("8", "2025-02-17", "气阀"),
    ("8", "2025-12-06", "活塞"),
    ("8", "2026-03-29", "气阀"),
    ("8", "2026-04-01", "气阀"),
    ("9", "2026-01-18", "活塞"),
]


def to_float(value):
    try:
        result = float(value)
    except (TypeError, ValueError):
        return None
    return result if np.isfinite(result) else None


def load(path: str) -> list[dict]:
    rows = []
    for row in csv.DictReader(open(path, encoding="utf-8-sig")):
        try:
            row["time"] = datetime.fromisoformat(row["sample_time"])
        except ValueError:
            continue
        row["bin"] = int(row["bin"])
        row["bin_start"] = float(row["bin_start"])
        row["p_span"] = to_float(row.get("p_span"))
        for name in BIN_METRICS + CYCLE_METRICS:
            row[name] = to_float(row.get(name))
        rows.append(row)
    return rows


def robust_scale(values: np.ndarray) -> float:
    if values.size == 0:
        return 0.0
    center = float(np.median(values))
    scale = float(np.median(np.abs(values - center))) * 1.4826
    if scale <= 1e-12:
        scale = float(np.std(values))
    return scale


def cycle_count(rows: list[dict]) -> int:
    return len({r["wave_file_id"] for r in rows})


def median_span(rows: list[dict]) -> float:
    spans = [r["p_span"] for r in rows if r.get("p_span") is not None]
    return float(np.median(spans)) if spans else 0.0


def fallback_segments(rows: list[dict], merge_gap: int) -> list[tuple]:
    days = sorted({r["time"].date() for r in rows})
    return build_segments(days, merge_gap)


def collect_before(rows, t0, segments, min_days):
    segs = sorted([s for s in segments if s[0] <= t0.date()], key=lambda s: s[0], reverse=True)
    chosen, days = [], set()
    for start, end in segs:
        part = [r for r in rows if start <= r["time"].date() <= end and r["time"] < t0]
        chosen.extend(part)
        days |= {r["time"].date() for r in part}
        if len(days) >= min_days:
            break
    return sorted(chosen, key=lambda r: r["time"]), days


def collect_after(rows, t0, segments, min_days, max_gap):
    """First running segment after the event; extend only over small gaps."""
    segs = sorted([s for s in segments if s[1] >= t0.date()], key=lambda s: s[0])
    chosen, days = [], set()
    prev_end = None
    for start, end in segs:
        if prev_end is not None and (start - prev_end).days > max_gap:
            break
        part = [r for r in rows if start <= r["time"].date() <= end and r["time"] >= t0]
        chosen.extend(part)
        days |= {r["time"].date() for r in part}
        prev_end = end
        if len(days) >= min_days:
            break
    return sorted(chosen, key=lambda r: r["time"]), days


def window_rows(rows, t0, segments, pre_days, post_days, max_gap):
    pre, pre_set = collect_before(rows, t0, segments, pre_days)
    post, post_set = collect_after(rows, t0, segments, post_days, max_gap)
    return pre, post, pre_set, post_set


def compare_bins(pre, post, min_cycles):
    """Per-bin level/shape comparison (BIN_METRICS)."""
    pre_bins = defaultdict(lambda: defaultdict(list))
    post_bins = defaultdict(lambda: defaultdict(list))
    bin_start_of = {}
    for row in pre:
        bin_start_of[row["bin"]] = row["bin_start"]
        for name in BIN_METRICS:
            if row[name] is not None:
                pre_bins[row["bin"]][name].append(row[name])
    for row in post:
        for name in BIN_METRICS:
            if row[name] is not None:
                post_bins[row["bin"]][name].append(row[name])

    all_bins = sorted(set(pre_bins) & set(post_bins))
    changes = []
    for name in BIN_METRICS:
        items = []
        for b in all_bins:
            a = np.asarray(pre_bins[b][name])
            c = np.asarray(post_bins[b][name])
            if a.size < min_cycles or c.size < min_cycles:
                continue
            pre_med = float(np.median(a))
            post_med = float(np.median(c))
            items.append({
                "bin": b, "bin_start": bin_start_of.get(b, b * 10.0), "metric": name,
                "pre_med": pre_med, "post_med": post_med, "delta": post_med - pre_med,
                "scale": robust_scale(a), "pre_n": int(a.size), "post_n": int(c.size),
            })
        if len(items) < 4:
            continue
        scales = [it["scale"] for it in items if it["scale"] > 1e-12]
        point_scale = float(np.median(scales)) if scales else 0.0
        floor = FLOOR_FRAC * point_scale
        pre_center = float(np.median([it["pre_med"] for it in items]))
        global_delta = float(np.median([it["delta"] for it in items]))
        global_z = global_delta / point_scale if point_scale > 1e-12 else None
        for it in items:
            denom = max(it["scale"], floor)
            it["local_delta"] = it["delta"] - global_delta
            it["local_z"] = it["local_delta"] / denom if denom > 1e-12 else None
            it["rel"] = it["local_delta"] / abs(pre_center) if abs(pre_center) > 1e-9 else None
            it["global_delta"] = global_delta
            it["global_z"] = global_z
            it["pre_center"] = pre_center
        changes.extend(items)
    return changes


def compare_cycle(pre, post, min_cycles):
    """Whole-cycle (valve timing / plateau) comparison, one record per metric."""
    records = []
    for name in CYCLE_METRICS:
        a = np.asarray([r[name] for r in pre if r[name] is not None])
        b = np.asarray([r[name] for r in post if r[name] is not None])
        if a.size < min_cycles or b.size < min_cycles:
            continue
        if name in ANGLE_METRICS:  # circular
            pre_med = circ_median(a)
            post_med = circ_median(b)
            delta = circ_diff(pre_med, post_med)
            residual = np.asarray([circ_diff(pre_med, v) for v in a])
            scale = max(robust_scale(residual), ANGLE_SCALE_FLOOR)
            rel = None
        else:
            pre_med = float(np.median(a))
            post_med = float(np.median(b))
            delta = post_med - pre_med
            scale = robust_scale(a)
            rel = delta / abs(pre_med) if abs(pre_med) > 1e-9 else None
        records.append({
            "metric": name, "pre_med": pre_med, "post_med": post_med, "delta": delta,
            "scale": scale, "z": delta / scale if scale > 1e-12 else None,
            "rel": rel, "pre_n": int(a.size), "post_n": int(b.size),
        })
    return records


def trend_bins(pre, min_cycles):
    series = defaultdict(lambda: defaultdict(list))
    bin_start_of = {}
    for row in pre:
        bin_start_of[row["bin"]] = row["bin_start"]
        for name in BIN_METRICS:
            if row[name] is not None:
                series[row["bin"]][name].append((row["time"], row[name]))
    results = []
    for name in BIN_METRICS:
        items = []
        for b, per_name in series.items():
            points = sorted(per_name.get(name, []))
            if len(points) < max(min_cycles, 8):
                continue
            x = np.asarray([(p[0] - points[0][0]).total_seconds() / 86400.0 for p in points])
            y = np.asarray([p[1] for p in points])
            if np.ptp(x) <= 0:
                continue
            slope = float(np.polyfit(x, y, 1)[0])
            items.append({
                "bin": b, "bin_start": bin_start_of.get(b, b * 10.0), "metric": name,
                "total_drift": slope * float(np.ptp(x)), "scale": robust_scale(y), "n": len(points),
            })
        if len(items) < 4:
            continue
        scales = [it["scale"] for it in items if it["scale"] > 1e-12]
        floor = FLOOR_FRAC * (float(np.median(scales)) if scales else 0.0)
        global_drift = float(np.median([it["total_drift"] for it in items]))
        for it in items:
            denom = max(it["scale"], floor)
            it["local_drift"] = it["total_drift"] - global_drift
            it["local_mad"] = it["local_drift"] / denom if denom > 1e-12 else None
        results.extend(items)
    return results


def trend_cycle(pre, min_cycles):
    results = []
    for name in CYCLE_METRICS:
        if name in ANGLE_METRICS:  # circular; skip linear trend
            continue
        points = sorted((r["time"], r[name]) for r in pre if r[name] is not None)
        if len(points) < max(min_cycles, 8):
            continue
        x = np.asarray([(p[0] - points[0][0]).total_seconds() / 86400.0 for p in points])
        y = np.asarray([p[1] for p in points])
        if np.ptp(x) <= 0:
            continue
        slope = float(np.polyfit(x, y, 1)[0])
        scale = robust_scale(y)
        drift = slope * float(np.ptp(x))
        results.append({"metric": name, "total_drift": drift, "scale": scale,
                        "mad": drift / scale if scale > 1e-12 else None, "n": len(points)})
    return results


def point_name(unit, cylinder, side):
    return f"{unit}号机组{cylinder}缸压力{side}"


def bin_label(bin_start, bin_deg):
    return f"{bin_start:g}-{bin_start + bin_deg:g}°"


def main() -> None:
    parser = argparse.ArgumentParser(description="同点 分箱+整周期 运行段前后对比 + 趋势")
    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
    parser.add_argument("--csv", nargs="+", required=True, help="angle_bins.py 长表 CSV（可多个）")
    parser.add_argument("--pre-days", type=int, default=15, help="前段至少累计的运行日（不足也可）")
    parser.add_argument("--post-days", type=int, default=15, help="后段目标运行日（不足也可）")
    parser.add_argument("--trend-days", type=int, default=30, help="故障前趋势至少累计的运行日")
    parser.add_argument("--max-gap", type=int, default=30, help="后段累加允许跨越的最大停机天数")
    parser.add_argument("--merge-gap", type=int, default=2, help="运行段桥接允许的停机天数")
    parser.add_argument("--min-cycles", type=int, default=3, help="每段最少周期数")
    parser.add_argument("--z", type=float, default=3.0, help="变化阈值（MAD 倍数）")
    parser.add_argument("--normalize", choices=("none", "span"), default="none",
                        help="span: 用本周期 p_span 归一再比（压工况）")
    parser.add_argument("--top", type=int, default=15, help="每张表打印条数")
    parser.add_argument("--out", default="", help="明细 CSV 输出路径")
    args = parser.parse_args()
    global Z
    Z = args.z

    rows = []
    for path in args.csv:
        part = load(path)
        print(f"  {path}: {len(part):,} 行")
        rows.extend(part)
    apply_normalization(rows, args.normalize)
    if args.normalize != "none":
        print(f"归一化: {args.normalize}")
    starts = sorted({r["bin_start"] for r in rows})
    bin_deg = round(starts[1] - starts[0], 6) if len(starts) >= 2 else 10.0
    print(f"载入 {len(rows):,} 行，测点 {len({r['point_name'] for r in rows})}，角度箱宽 {bin_deg:g}°")

    try:
        segments = unit_segments(args.unit, args.merge_gap)
        print(f"运行段(rpm>0): {len(segments)} 个")
    except Exception as error:
        print(f"运行段取 DB 失败({error})，改用数据日期间隔")
        segments = fallback_segments(rows, args.merge_gap)

    detail_rows: list[dict] = []
    for unit, date, kind in EVENTS:
        if unit != args.unit:
            continue
        t0 = datetime.fromisoformat(date)
        points = sorted({(r["cylinder"], r["side"]) for r in rows})
        bin_changes, cycle_changes, bin_trends, cycle_trends = [], [], [], []
        coverage = None
        for cylinder, side in points:
            point_rows = [r for r in rows if r["cylinder"] == cylinder and r["side"] == side]
            pre, post, pre_set, post_set = window_rows(point_rows, t0, segments, args.pre_days, args.post_days, args.max_gap)
            if cycle_count(pre) < args.min_cycles or cycle_count(post) < args.min_cycles:
                continue
            if median_span(pre) < MIN_SPAN or median_span(post) < MIN_SPAN:
                continue
            if coverage is None:
                coverage = (min(r["time"] for r in pre), max(r["time"] for r in pre), pre_set,
                            min(r["time"] for r in post), max(r["time"] for r in post), post_set)
            point = f"{cylinder}缸{side}"
            label = point_name(unit, cylinder, side)
            pre_dates = sorted({r["time"] for r in pre})
            post_dates = sorted({r["time"] for r in post})
            wmeta = {
                "pre_start": str(pre_dates[0].date()), "pre_end": str(pre_dates[-1].date()),
                "post_start": str(post_dates[0].date()), "post_end": str(post_dates[-1].date()),
                "gap_days": (post_dates[0].date() - pre_dates[-1].date()).days,
                "pre_span_days": (pre_dates[-1].date() - pre_dates[0].date()).days,
            }

            changes = compare_bins(pre, post, args.min_cycles)
            strong_bins = defaultdict(int)
            for c in changes:
                if c["local_z"] is not None and abs(c["local_z"]) >= Z:
                    strong_bins[c["bin"]] += 1
            for c in changes:
                if c["local_z"] is None or abs(c["local_z"]) < Z:
                    continue
                bin_changes.append({
                    "point": point, "label": label, "bin": c["bin"], "bin_start": c["bin_start"],
                    "type": "局部波动" if c["metric"] in VOLATILE else "局部均值",
                    "metric": c["metric"], "pre": c["pre_med"], "post": c["post_med"],
                    "delta": c["local_delta"], "mad": c["local_z"], "rel": c["rel"],
                    "candidate": strong_bins[c["bin"]] >= 2,
                })
            for name in BIN_METRICS:
                items = [c for c in changes if c["metric"] == name]
                if not items or items[0]["global_z"] is None or abs(items[0]["global_z"]) < Z:
                    continue
                pre_c = items[0]["pre_center"]
                bin_changes.append({
                    "point": point, "label": label, "bin": -1, "bin_start": -1.0,
                    "type": "整体波动" if name in VOLATILE else "整体水平", "metric": name,
                    "pre": pre_c, "post": pre_c + items[0]["global_delta"], "delta": items[0]["global_delta"],
                    "mad": items[0]["global_z"], "rel": None, "candidate": False,
                })

            for c in compare_cycle(pre, post, args.min_cycles):
                if c["z"] is None or abs(c["z"]) < Z:
                    continue
                cycle_changes.append({
                    "point": point, "label": label,
                    "type": "整周期波动" if c["metric"] in CYCLE_VOLATILE else "整周期水平",
                    "metric": c["metric"], "pre": c["pre_med"], "post": c["post_med"],
                    "delta": c["delta"], "mad": c["z"], "rel": c["rel"],
                })

            for t in trend_bins(pre, args.min_cycles):
                if t.get("local_mad") is None or abs(t["local_mad"]) < Z:
                    continue
                bin_trends.append({
                    "point": point, "label": label, "bin": t["bin"], "bin_start": t["bin_start"],
                    "type": "局部趋势" + ("波动" if t["metric"] in VOLATILE else "均值"),
                    "metric": t["metric"], "drift": t["local_drift"], "mad": t["local_mad"],
                    "n": t["n"],
                })
            for t in trend_cycle(pre, args.min_cycles):
                if t["mad"] is None or abs(t["mad"]) < Z:
                    continue
                cycle_trends.append({
                    "point": point, "label": label,
                    "type": "整周期趋势" + ("波动" if t["metric"] in CYCLE_VOLATILE else "水平"),
                    "metric": t["metric"], "drift": t["total_drift"], "mad": t["mad"], "n": t["n"],
                })

        print(f"\n{'=' * 100}")
        print(f"== {unit}号机 {date} {kind}")
        if coverage:
            print(f"   前段 {coverage[0]:%Y-%m-%d}~{coverage[1]:%Y-%m-%d}（{len(coverage[2])}运行日）"
                  f"    后段 {coverage[3]:%Y-%m-%d}~{coverage[4]:%Y-%m-%d}（{len(coverage[5])}运行日）")
        if not (bin_changes or cycle_changes or bin_trends or cycle_trends):
            print("   事件前后无足够运行数据")
            continue

        bin_changes.sort(key=lambda r: (r["type"].startswith("整体"), -abs(r["mad"])))
        if bin_changes:
            print(f"   【分箱前后变化】{len(bin_changes)} 条（★=同角度箱≥2指标一致）")
            print(f"   {'#':>3}  {'测点':<9} {'角度箱':<10} {'类型':<8} {'指标':<11} {'前':>9} {'后':>9} {'变化':>9} {'相对':>8}")
            for i, r in enumerate(bin_changes[:args.top], 1):
                binst = "全箱" if r["bin"] < 0 else bin_label(r["bin_start"], bin_deg)
                rel = f"{r['rel'] * 100:+.0f}%" if r["rel"] is not None else "-"
                print(f"   {i:>3}  {r['point']:<9} {binst:<10} {r['type']:<8} {r['metric']:<11} "
                      f"{r['pre']:>9.4g} {r['post']:>9.4g} {r['mad']:>+8.1f}M {rel:>8} {'★' if r['candidate'] else ''}")

        if cycle_changes:
            cycle_changes.sort(key=lambda r: -abs(r["mad"]))
            print(f"   【整周期/阀时序 前后变化】{len(cycle_changes)} 条")
            print(f"   {'#':>3}  {'测点':<9} {'类型':<10} {'指标':<15} {'前':>10} {'后':>10} {'变化':>9} {'相对':>8}")
            for i, r in enumerate(cycle_changes[:args.top], 1):
                rel = f"{r['rel'] * 100:+.0f}%" if r["rel"] is not None else "-"
                print(f"   {i:>3}  {r['point']:<9} {r['type']:<10} {r['metric']:<15} "
                      f"{r['pre']:>10.4g} {r['post']:>10.4g} {r['mad']:>+8.1f}M {rel:>8}")

        if bin_trends:
            bin_trends.sort(key=lambda r: -abs(r["mad"]))
            print(f"   【故障前分箱趋势】{len(bin_trends)} 条")
            print(f"   {'#':>3}  {'测点':<9} {'角度箱':<10} {'类型':<12} {'指标':<11} {'总变化':>11} {'变化':>9}")
            for i, r in enumerate(bin_trends[:args.top], 1):
                print(f"   {i:>3}  {r['point']:<9} {bin_label(r['bin_start'], bin_deg):<10} {r['type']:<12} "
                      f"{r['metric']:<11} {r['drift']:>+11.4g} {r['mad']:>+8.1f}M")

        if cycle_trends:
            cycle_trends.sort(key=lambda r: -abs(r["mad"]))
            print(f"   【故障前整周期趋势】{len(cycle_trends)} 条")
            print(f"   {'#':>3}  {'测点':<9} {'类型':<14} {'指标':<15} {'总变化':>12} {'变化':>9}")
            for i, r in enumerate(cycle_trends[:args.top], 1):
                print(f"   {i:>3}  {r['point']:<9} {r['type']:<14} {r['metric']:<15} "
                      f"{r['drift']:>+12.4g} {r['mad']:>+8.1f}M")

        for r in bin_changes:
            detail_rows.append({
                **wmeta,
                "unit": unit, "date": date, "kind": kind, "family": "bin", "point": r["point"],
                "label": r["label"], "bin": r["bin"], "bin_start": r["bin_start"], "type": r["type"],
                "metric": r["metric"], "pre": r["pre"], "post": r["post"], "delta": r["delta"],
                "mad": r["mad"], "rel": r["rel"] if r["rel"] is not None else "",
                "candidate": int(r["candidate"]),
            })
        for r in cycle_changes:
            detail_rows.append({
                **wmeta,
                "unit": unit, "date": date, "kind": kind, "family": "cycle", "point": r["point"],
                "label": r["label"], "bin": "", "bin_start": "", "type": r["type"],
                "metric": r["metric"], "pre": r["pre"], "post": r["post"], "delta": r["delta"],
                "mad": r["mad"], "rel": r["rel"] if r["rel"] is not None else "", "candidate": "",
            })
        for r in bin_trends:
            detail_rows.append({
                **wmeta,
                "unit": unit, "date": date, "kind": kind, "family": "bin_trend", "point": r["point"],
                "label": r["label"], "bin": r["bin"], "bin_start": r["bin_start"], "type": r["type"],
                "metric": r["metric"], "pre": "", "post": "", "delta": r["drift"], "mad": r["mad"],
                "rel": "", "candidate": "",
            })
        for r in cycle_trends:
            detail_rows.append({
                **wmeta,
                "unit": unit, "date": date, "kind": kind, "family": "cycle_trend", "point": r["point"],
                "label": r["label"], "bin": "", "bin_start": "", "type": r["type"],
                "metric": r["metric"], "pre": "", "post": "", "delta": r["drift"], "mad": r["mad"],
                "rel": "", "candidate": "",
            })

    if args.out and detail_rows:
        out = Path(args.out)
        out.parent.mkdir(parents=True, exist_ok=True)
        with out.open("w", encoding="utf-8-sig", newline="") as handle:
            writer = csv.DictWriter(handle, fieldnames=list(detail_rows[0].keys()))
            writer.writeheader()
            writer.writerows(detail_rows)
        print(f"\n明细 -> {out}（{len(detail_rows)} 行）")
    print(f"时间：{datetime.now():%Y-%m-%d %H:%M:%S}")


if __name__ == "__main__":
    main()