Compressor/backend/indicator_features/trend.py

#!/usr/bin/env python3
"""Long-term indicator-feature deviation vs the scanner's anomaly segments.

For each cylinder/side it reduces the feature CSV to daily medians, scores them
as robust z against the point's own baseline (median/MAD over its valid days),
and marks days that fall inside a rule anomaly segment from cache/pks_scanner/
evolution (pressure 阶段段 / 振动异常段).  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

UNIT_BATCH = {"7": 30, "8": 31, "9": 32}
TREND_FEATURES = ("p_discharge", "pressure_ratio", "n_comp", "efficiency")


def parse_dt(value):
    return datetime.strptime(value[:19], "%Y-%m-%d %H:%M:%S")


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


def load_segments(evolution_dir, unit):
    """Rule anomaly segments for a unit from the two stage CSVs."""
    segments = []
    for name, kind in (("润滑油压力小时阶段段.csv", "压力"), ("振动异常段.csv", "振动")):
        path = Path(evolution_dir) / name
        if not path.is_file():
            continue
        for row in csv.DictReader(open(path, encoding="utf-8-sig")):
            if row["机组"] != f"{unit}号机":
                continue
            segments.append({
                "kind": kind, "level": row.get("压力阶段") or row.get("振动阶段"),
                "start": parse_dt(row["开始小时"]), "end": parse_dt(row["结束小时"]),
            })
    return segments


def mark_segments(day: str, segments):
    moment = datetime.strptime(day, "%Y-%m-%d")
    return [s["kind"] + "/" + str(s["level"]) for s in segments if s["start"] <= moment <= s["end"]]


def peer_sudden_changes(rows, feature, segments):
    """Day-to-day jumps of the peer-relative deviation (cancels operating level)."""
    by_key = defaultdict(dict)
    for row in rows:
        if row[feature] is not None:
            by_key[(row["stage"], row["side"], row["day"])][row["cylinder"]] = row[feature]
    series = defaultdict(dict)
    for (_stage, side, day), members in by_key.items():
        values = [v for v in members.values() if v is not None]
        if len(values) < 2:
            continue
        med = float(np.median(values))
        if abs(med) < 1e-9:
            continue
        for cylinder, value in members.items():
            if value is not None:
                series[(cylinder, side)][day] = (value - med) / med
    print(f"\n== 同级缸间偏差 突发变化（{feature}）==")
    for (cylinder, side), daily in sorted(series.items()):
        days = sorted(daily)
        if len(days) < 3:
            continue
        values = np.array([daily[d] for d in days], dtype=float)
        scale = float(np.median(np.abs(values - np.median(values)))) * 1.4826 or float(np.std(values))
        if scale <= 1e-9:
            continue
        for i in range(1, len(days)):
            gap = (datetime.strptime(days[i], "%Y-%m-%d") - datetime.strptime(days[i - 1], "%Y-%m-%d")).days
            if gap > 3:
                continue
            delta = values[i] - values[i - 1]
            if abs(delta) / scale >= 3.0:
                flags = mark_segments(days[i], segments)
                print(f"   {cylinder}缸{side}: 偏差 {values[i-1]*100:+.1f}% → {values[i]*100:+.1f}% 于 {days[i]} ({delta/scale:+.1f}×)  {' '.join(flags)}")


def _fit_drift(days, values):
    x = np.array([(datetime.strptime(d, "%Y-%m-%d") - datetime.strptime(days[0], "%Y-%m-%d")).days for d in days], dtype=float)
    y = np.asarray(values, dtype=float)
    if x[-1] <= 0 or len(y) < 8:
        return None
    slope, intercept = np.polyfit(x, y, 1)
    pred = slope * x + 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
    center = float(np.median(y))
    rel30 = (slope * 30.0 / center) if abs(center) > 1e-9 else 0.0
    return rel30, r2, float(y[0]), float(y[-1])


def slow_trends(rows, features=("p_discharge", "pressure_ratio", "efficiency", "n_comp", "indicated_work")):
    """Slow, weak drifts: linear fit per (cylinder, side) and per peer-relative series."""
    print("\n== 慢趋势：自基线（斜率 %/30天，要求 |%|>=3 且 R2>=0.3） ==")
    by_series = defaultdict(list)
    for row in rows:
        by_series[(row["cylinder"], row["side"])].append(row)
    for (cylinder, side), series in sorted(by_series.items()):
        days = sorted({r["day"] for r in series})
        for feature in features:
            daily = {d: np.median([r[feature] for r in series if r["day"] == d and r[feature] is not None]) for d in days}
            present = [d for d in days if daily.get(d) is not None]
            fit = _fit_drift(present, [daily[d] for d in present])
            if fit and abs(fit[0]) >= 0.03 and fit[1] >= 0.3:
                print(f"   {cylinder}缸{side} {feature:<16} {fit[0]*100:+.1f}%/30天 R2={fit[1]:.2f} ({fit[2]:.3f}->{fit[3]:.3f})")

    print("\n== 慢趋势：同级 peer 相对（%/30天） ==")
    by_ps = defaultdict(dict)
    for row in rows:
        if row["pressure_ratio"] is not None:
            by_ps[(row["stage"], row["side"], row["day"])][row["cylinder"]] = row["pressure_ratio"]
    rel = defaultdict(dict)
    for (stage, side, day), members in by_ps.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():
            rel[(cylinder, side)][day] = (value - med) / med
    for (cylinder, side), daily in sorted(rel.items()):
        days = sorted(daily)
        fit = _fit_drift(days, [daily[d] for d in days])
        if fit and abs(fit[0]) >= 0.03 and fit[1] >= 0.3:
            print(f"   {cylinder}缸{side} pressure_ratio(peer) {fit[0]*100:+.1f}%/30天 R2={fit[1]:.2f} ({fit[2]*100:+.1f}%->{fit[3]*100:+.1f}%)")


def main():
    parser = argparse.ArgumentParser(description="功图特征长期偏离检测（对齐 evolution 异常段）")
    parser.add_argument("csv_path")
    parser.add_argument("--evolution-dir", default=str(Path(__file__).resolve().parents[2] / "cache" / "pks_scanner" / "evolution"))
    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
    args = parser.parse_args()

    rows = load_features(args.csv_path)
    segments = load_segments(args.evolution_dir, args.unit)
    by_series = defaultdict(list)
    for row in rows:
        by_series[(row["cylinder"], row["side"])].append(row)

    print(f"文件 {args.csv_path}  行 {len(rows)}  evolution段 {len(segments)}")
    for (cylinder, side), series in sorted(by_series.items()):
        days = sorted({r["day"] for r in series})
        print(f"\n== {cylinder}缸{side}  ({len(days)} 天)")
        for feature in TREND_FEATURES:
            daily = {}
            for day in days:
                values = [r[feature] for r in series if r["day"] == day and r[feature] is not None]
                if values:
                    daily[day] = float(np.median(values))
            if len(daily) < 5:
                continue
            values = np.array([daily[d] for d in days if d in daily], dtype=float)
            center = float(np.median(values))
            mad = float(np.median(np.abs(values - center)))
            scale = mad * 1.4826 or float(np.std(values))
            if scale <= 1e-9:
                continue
            print(f"   {feature:<16} 基线 {center:.3f}±{scale:.3f}")
            for day in days:
                if day not in daily:
                    continue
                z = (daily[day] - center) / scale
                flags = mark_segments(day, segments)
                if abs(z) >= 2.0 or flags:
                    tag = "  ".join(flags)
                    print(f"      {day}  {daily[day]:>8.3f}  z={z:+.2f}  {tag}")
            # Abrupt day-to-day change (sudden step), independent of any trend.
            present = [d for d in days if d in daily]
            jumps = []
            for i in range(1, len(present)):
                gap = (datetime.strptime(present[i], "%Y-%m-%d") - datetime.strptime(present[i - 1], "%Y-%m-%d")).days
                if gap <= 3:
                    jumps.append((present[i], daily[present[i]] - daily[present[i - 1]]))
            if jumps:
                day_jump, delta = max(jumps, key=lambda item: abs(item[1]))
                if abs(delta) / scale >= 3.0:
                    flags = mark_segments(day_jump, segments)
                    tag = "  ".join(flags)
                    print(f"      [突发] {feature}: Δ={delta:+.3f} ({delta / scale:+.1f}×基线尺度) 于 {day_jump}  {tag}")
    peer_sudden_changes(rows, "p_discharge", segments)
    slow_trends(rows)


if __name__ == "__main__":
    main()