Compressor/backend/indicator_features/event_scan.py

#!/usr/bin/env python3
"""Scan every repair event for a single-cylinder slow drift that leads it.

For each event it looks back --days before the repair, builds the daily
peer-relative series (pressure_ratio / p_discharge / efficiency within the same
stage and side) and the self series, fits a linear drift, and ranks by strength.
Read-only; uses feature CSVs produced by extract_from_db.py.
"""
from __future__ import annotations

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

import numpy as np

EVENTS = [
    ("7", "2025-08-08", "气阀"),
    ("7", "2025-10-07", "气阀"),
    ("7", "2026-05-21", "气阀"),
    ("8", "2025-02-17", "气阀"),
    ("8", "2025-12-06", "活塞"),
    ("8", "2026-03-10", "十字头"),
    ("8", "2026-03-29", "气阀"),
    ("8", "2026-04-01", "气阀"),
    ("9", "2026-01-18", "活塞"),
]
STAGE_MEMBERS = {"一级": ("二", "四", "六"), "二级": ("三", "五"), "三级": ("一",)}
CYL_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一级", "五": "二级", "六": "一级"}
FEATURES = ("pressure_ratio", "p_discharge", "efficiency")


def load(path):
    rows = []
    for row in csv.DictReader(open(path, encoding="utf-8-sig")):
        row["day"] = row["sample_time"][:10]
        for name in FEATURES:
            try:
                row[name] = float(row[name])
            except (TypeError, ValueError):
                row[name] = None
        rows.append(row)
    return rows


def daily_median(rows, cylinder, side, feature):
    groups = defaultdict(list)
    for row in rows:
        if row["cylinder"] == cylinder and row["side"] == side and row[feature] is not None:
            groups[row["day"]].append(row[feature])
    return {day: float(np.median(v)) for day, v in groups.items()}


def peer_relative(rows, feature):
    by_key = defaultdict(dict)
    for row in rows:
        if row[feature] is not None:
            by_key[(row["day"], row["side"])][row["cylinder"]] = row[feature]
    series = defaultdict(dict)
    for (day, side), members in by_key.items():
        clean = {c: v for c, v in members.items() if v is not None}
        if len(clean) < 2:
            continue
        med = float(np.median(list(clean.values())))
        if abs(med) < 1e-9:
            continue
        for cylinder, value in clean.items():
            series[(cylinder, side)][day] = (value - med) / med
    return series


def fit(daily):
    days = sorted(daily)
    if len(days) < 10:
        return None
    dates = [datetime.strptime(d, "%Y-%m-%d") for d in days]
    span = (dates[-1] - dates[0]).days
    if span <= 0 or span / (len(dates) - 1) > 15:
        return None
    xd = np.array([(d - dates[0]).days for d in dates], dtype=float)
    y = np.array([daily[d] for d in days], dtype=float)
    if xd[-1] <= 0:
        return None
    slope, intercept = np.polyfit(xd, y, 1)
    pred = slope * xd + intercept
    ss_res = float(np.sum((y - pred) ** 2))
    ss_tot = float(np.sum((y - y.mean()) ** 2))
    r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0.0
    scale = float(np.median(np.abs(y - np.median(y)))) * 1.4826 or float(np.std(y))
    if scale <= 1e-9:
        return None
    total_mads = abs(slope) * xd[-1] / scale  # total drift over the window in MADs
    return total_mads, r2, len(days), float(y[0]), float(y[-1])


def main():
    parser = argparse.ArgumentParser(description="扫描所有维修事件的单缸慢漂移")
    parser.add_argument("--data-dir", default=str(Path(__file__).resolve().parents[1] / "cache" / "indicator"))
    parser.add_argument("--days", type=int, default=180, help="维修前回看天数")
    args = parser.parse_args()

    cache = {}
    for unit in ("7", "8", "9"):
        path = Path(args.data_dir) / f"all_{unit}_2025_2026.csv"
        cache[unit] = load(path) if path.is_file() else []

    for unit, date, kind in EVENTS:
        event_day = datetime.strptime(date, "%Y-%m-%d")
        start = (event_day - timedelta(days=args.days)).strftime("%Y-%m-%d")
        rows = [r for r in cache[unit] if start <= r["day"] <= date]
        days = sorted({r["day"] for r in rows})
        print(f"\n== {unit}号机 {date} {kind}  回看{args.days}天  数据日 {len(days)} ({days[0] if days else '-'}~{days[-1] if days else '-'})")
        if len(days) < 8:
            print("   数据不足，跳过")
            continue
        ranked = []
        for feature in FEATURES:
            for (cylinder, side), daily in peer_relative(rows, feature).items():
                result = fit(daily)
                if result and result[0] >= 2.0 and result[1] >= 0.4:
                    ranked.append((result[0], f"peer {cylinder}缸{side} {feature}", result))
            for cylinder in ("一", "二", "三", "四", "五", "六"):
                for side in ("盖侧", "轴侧"):
                    result = fit(daily_median(rows, cylinder, side, feature))
                    if result and result[0] >= 3.0 and result[1] >= 0.5:
                        ranked.append((result[0], f"self {cylinder}缸{side} {feature}", result))
        ranked.sort(reverse=True)
        if not ranked:
            print("   无显著慢漂移")
        for _, label, (mads, r2, n, first, last) in ranked[:5]:
            print(f"   {label:<34} 漂移{mads:4.1f}MAD R2={r2:.2f} (n={n}, {first:.3f}->{last:.3f})")


if __name__ == "__main__":
    main()