Compressor/backend/indicator_features/common_residual_analysis.py

#!/usr/bin/env python3
"""Separate common operating-condition changes from local diagram changes."""
from __future__ import annotations

import argparse
import csv
import sys
from collections import defaultdict
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))
from event_bins import EVENTS  # noqa: E402

BAD = {("7", "四", "盖侧"), ("7", "六", "盖侧"), ("8", "六", "盖侧")}
SCALARS = {
    "p_suction": "吸气压力", "p_discharge": "排气压力", "p_span": "压力幅度",
    "suction_frac": "吸气占比", "discharge_frac": "排气占比",
    "compression_slope": "压缩斜率", "expansion_slope": "膨胀斜率",
    "corner_sharpness": "拐角变化",
}


def load(path: str) -> list[dict]:
    with open(path, encoding="utf-8-sig", newline="") as handle:
        rows = list(csv.DictReader(handle))
    for row in rows:
        row["time"] = datetime.fromisoformat(row["sample_time"])
        row["span"] = float(row.get("p_span") or 0)
        for key in list(SCALARS) + [k for k in row if k.startswith("prof_")]:
            try:
                row[key] = float(row[key])
            except (TypeError, ValueError):
                row[key] = np.nan
    return rows


def median(rows: list[dict], key: str) -> float | None:
    values = [r[key] for r in rows if np.isfinite(r[key])]
    return float(np.median(values)) if values else None


def profile_median(rows: list[dict], names: list[str]) -> np.ndarray | None:
    if not rows:
        return None
    values = np.asarray([[r[n] for n in names] for r in rows], dtype=float)
    values = values[np.all(np.isfinite(values), axis=1)]
    return np.median(values, axis=0) if len(values) else None


def pct(value: float | None) -> str:
    return "" if value is None else f"{value:.5f}"


def main() -> None:
    parser = argparse.ArgumentParser(description="公共工况变化与局部功图残差")
    parser.add_argument("--csv", nargs="+", required=True)
    parser.add_argument("--unit", nargs="+", required=True)
    parser.add_argument("--pre-days", type=int, default=30)
    parser.add_argument("--post-days", type=int, default=30)
    parser.add_argument("--outdir", required=True)
    parser.add_argument("--event", nargs=3, action="append", metavar=("UNIT", "DATE", "KIND"),
                        help="只分析指定事件；可重复传入")
    args = parser.parse_args()

    by_unit = {unit: load(path) for unit, path in zip(args.unit, args.csv)}
    outdir = Path(args.outdir)
    outdir.mkdir(parents=True, exist_ok=True)
    all_rows: list[dict] = []
    summary: list[dict] = []
    normal: list[dict] = []

    events = [tuple(item) for item in args.event] if args.event else EVENTS
    for unit, date_text, kind in events:
        source = [r for r in by_unit.get(unit, []) if r["span"] >= 0.3]
        t0 = datetime.fromisoformat(date_text)
        pre_start = t0 - timedelta(days=args.pre_days)
        post_end = t0 + timedelta(days=args.post_days)
        points = sorted({(r["cylinder"], r["side"]) for r in source if (unit, r["cylinder"], r["side"]) not in BAD})
        profiles = sorted(k for k in source[0] if k.startswith("prof_")) if source else []
        deltas: dict[tuple[str, str], np.ndarray] = {}
        profile_pairs = {}
        scalar_changes: dict[tuple[str, str], dict[str, float]] = {}
        for cylinder, side in points:
            rows = [r for r in source if r["cylinder"] == cylinder and r["side"] == side]
            pre = [r for r in rows if pre_start <= r["time"] < t0]
            post = [r for r in rows if t0 <= r["time"] <= post_end]
            if not pre or not post:
                continue
            before = profile_median(pre, profiles)
            after = profile_median(post, profiles)
            if before is None or after is None:
                continue
            delta = after - before
            deltas[(cylinder, side)] = delta
            profile_pairs[(cylinder, side)] = (before, after, len(pre), len(post))
            changes = {}
            for key in SCALARS:
                a, b = median(pre, key), median(post, key)
                changes[key] = np.nan if a is None or b is None else b - a
            scalar_changes[(cylinder, side)] = changes

        if not deltas:
            summary.append({"机组": unit, "事件": date_text, "故障类型": kind, "可比较测点": 0,
                            "公共形状变化": "", "局部残差最大": "", "判断": "修前无有效功图"})
            continue

        matrix = np.vstack(list(deltas.values()))
        common = np.median(matrix, axis=0)
        common_norm = float(np.linalg.norm(common))
        local_scores = {}
        for point, delta in deltas.items():
            residual = delta - common
            local = float(np.linalg.norm(residual))
            total = float(np.linalg.norm(delta))
            local_scores[point] = (local, total, residual)
            cylinder, side = point
            before, after, pre_n, post_n = profile_pairs[point]
            for idx, name in enumerate(profiles):
                all_rows.append({
                    "机组": unit, "事件": date_text, "故障类型": kind,
                    "测点": f"{cylinder}缸{side}", "角度": f"{idx * 360 / len(profiles):.1f}°",
                    "前后变化": round(float(after[idx] - before[idx]), 6),
                    "公共变化": round(float(common[idx]), 6),
                    "局部残差": round(float(after[idx] - before[idx] - common[idx]), 6),
                    "修前周期": pre_n, "修后周期": post_n,
                })

        ranked = sorted(local_scores.items(), key=lambda item: item[1][0], reverse=True)
        top = ranked[0]
        ratio = top[1][0] / max(common_norm, 1e-12)
        judgment = "局部差异优先" if ratio >= 1.2 else "公共工况变化优先"
        summary.append({
            "机组": unit, "事件": date_text, "故障类型": kind,
            "可比较测点": len(deltas), "公共形状变化": round(common_norm, 4),
            "局部残差最大": f"{top[0][0]}缸{top[0][1]}",
            "最大局部残差": round(top[1][0], 4), "该点总变化": round(top[1][1], 4),
            "局部/公共": round(ratio, 2), "判断": judgment,
        })

        # Normal candidates: closest to the pre-event median profile for each point.
        for point, (before, after, pre_n, post_n) in profile_pairs.items():
            cylinder, side = point
            rows = [r for r in source if r["cylinder"] == cylinder and r["side"] == side]
            pre_rows = [r for r in rows if pre_start <= r["time"] < t0]
            distances = []
            for row in pre_rows:
                vector = np.asarray([row[n] for n in profiles], dtype=float)
                if np.all(np.isfinite(vector)):
                    distances.append((float(np.linalg.norm(vector - before)), row))
            if distances:
                distance, row = min(distances, key=lambda item: item[0])
                local, total, _ = local_scores[point]
                normal.append({
                    "机组": unit, "事件": date_text, "测点": f"{cylinder}缸{side}",
                    "类别": "正常候选", "日期": row["sample_time"],
                    "距修前模板": round(distance, 4), "事件后该点变化": round(total, 4),
                    "局部残差": round(local, 4),
                })

        # One compact scalar table per event and point.
        for point, changes in scalar_changes.items():
            for key, label in SCALARS.items():
                all_rows.append({
                    "机组": unit, "事件": date_text, "故障类型": kind,
                    "测点": f"{point[0]}缸{point[1]}", "角度": label,
                    "前后变化": "" if np.isnan(changes[key]) else round(changes[key], 6),
                    "公共变化": "", "局部残差": "", "修前周期": "", "修后周期": "",
                })

    def write(path: Path, rows: list[dict]) -> None:
        fields = []
        for row in rows:
            for field in row:
                if field not in fields:
                    fields.append(field)
        with path.open("w", encoding="utf-8-sig", newline="") as handle:
            writer = csv.DictWriter(handle, fieldnames=fields)
            writer.writeheader()
            writer.writerows(rows)

    write(outdir / "三事件_公共与局部变化.csv", all_rows)
    write(outdir / "三事件_公共工况判断.csv", summary)
    write(outdir / "三事件_正常功图候选.csv", normal)
    print(f"公共/局部变化 -> {outdir / '三事件_公共与局部变化.csv'}")
    print(f"事件判断 -> {outdir / '三事件_公共工况判断.csv'}")
    print(f"正常候选 -> {outdir / '三事件_正常功图候选.csv'}")
    for row in summary:
        print(row)


if __name__ == "__main__":
    main()