Compressor/_tmp_temperature/backtest.py

#!/usr/bin/env python3
"""Backtest the cylinder-temperature scan against labelled 4#/5#/6# faults.

Reads ``_tmp_temperature/out/温度小时特征.csv`` and ``温度异常段.csv`` and
reports, for every temperature-related fault:

* the maximum same-stage spread change / peer deviation in the 7 days before;
* whether a mild/abnormal/severe segment hit the pre-fault window;
* the lead time.

It also prints normal-period percentiles and a small threshold sweep so the
"how many degC is a fault" question can be answered from data instead of a
hand-picked 3 degC rule.
"""
from __future__ import annotations

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

import pandas as pd

ROOT = Path(__file__).resolve().parents[1]
if str(ROOT) not in sys.path:
    sys.path.insert(0, str(ROOT))

from models.common.db import connect

OUT = Path(__file__).resolve().parent / "out"
LOOKBACK_DAYS = 7


def load_events() -> pd.DataFrame:
    connection = connect()
    try:
        with connection.cursor() as cursor:
            cursor.execute(
                "SELECT id, unit_name, fault_date, fault_category, fault_description "
                "FROM compressor_fault "
                "WHERE unit_name IN ('4号机组','5号机组','6号机组') "
                "AND (fault_description LIKE %s OR fault_description LIKE %s "
                "     OR fault_description LIKE %s OR fault_description LIKE %s) "
                "ORDER BY unit_name, fault_date",
                ("%缸温度%", "%排气温度%", "%排温%", "%温度异常%"),
            )
            rows = cursor.fetchall()
    finally:
        connection.close()
    frame = pd.DataFrame(rows)
    frame["event_date"] = pd.to_datetime(frame["fault_date"])
    return frame


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--dir", default=str(OUT))
    args = parser.parse_args()
    directory = Path(args.dir)

    features = pd.read_csv(directory / "温度小时特征.csv", low_memory=False)
    features["dt"] = pd.to_datetime(features["小时"])
    features["unit"] = features["机组"].str.replace("号机", "号机组", regex=False)
    for column in ("级差变化", "最大同级偏离", "最大同级偏低", "最大稳健z", "级差", "最高温度"):
        features[column] = pd.to_numeric(features[column], errors="coerce")

    segments = pd.read_csv(directory / "温度异常段.csv")
    segments = segments[segments["特征类型"] != "高频波动"].copy()
    segments["start"] = pd.to_datetime(segments["开始小时"])
    segments["end"] = pd.to_datetime(segments["结束小时"])
    segments["unit"] = segments["机组"].str.replace("号机", "号机组", regex=False)

    channel_path = directory / "温度通道异常段.csv"
    if channel_path.exists():
        channels = pd.read_csv(channel_path)
        channels["start"] = pd.to_datetime(channels["开始小时"])
        channels["end"] = pd.to_datetime(channels["结束小时"])
        channels["unit"] = channels["机组"].str.replace("号机", "号机组", regex=False)
    else:
        channels = pd.DataFrame(columns=["unit", "start", "end", "异常类型", "点位描述"])

    events = load_events()
    print(f"温度类故障事件：{len(events)} 起\n")

    def _closest_start(frame, moment):
        if frame.empty:
            return ""
        values = frame["start"].tolist()
        before = [value for value in values if value <= moment]
        chosen = max(before) if before else min(values)
        return chosen.strftime("%Y-%m-%d %H:%M")

    rows = []
    for _, event in events.iterrows():
        unit = event["unit_name"]
        date = event["event_date"]
        start = date - timedelta(days=LOOKBACK_DAYS)
        end = date + timedelta(days=1)
        window = features[(features["unit"] == unit) & (features["dt"] >= start) & (features["dt"] <= end)]
        hits = segments[(segments["unit"] == unit) & (segments["start"] <= end) & (segments["end"] >= start)]
        channel_hits = channels[(channels["unit"] == unit) & (channels["start"] <= end) & (channels["end"] >= start)]
        starts = []
        if not hits.empty:
            starts.extend(hits["start"].tolist())
        if not channel_hits.empty:
            starts.extend(channel_hits["start"].tolist())
        if starts:
            before = [value for value in starts if value <= end]
            first = max(before) if before else min(starts)
        else:
            first = pd.NaT
        lead = (date - first).total_seconds() / 3600 if starts else None

        def peak(stage: str, column: str):
            part = window[window["级"] == stage][column].dropna()
            return round(float(part.max()), 2) if not part.empty else None

        def peak_row(stage: str, column: str):
            part = window[(window["级"] == stage) & window[column].notna()]
            if part.empty:
                return None
            row = part.loc[part[column].idxmax()]
            return row

        stage = "一级" if "一级" in str(event["fault_description"]) else "二级" if "二级" in str(event["fault_description"]) else ""
        stages = [stage] if stage else ["一级", "二级"]
        best = None
        for candidate in stages:
            peak_peer = peak(candidate, "最大同级偏离")
            peak_spread = peak(candidate, "级差变化")
            score = max(value for value in (peak_peer, peak_spread) if value is not None) if (peak_peer is not None or peak_spread is not None) else -1
            if best is None or score > best[0]:
                best = (score, candidate, peak_peer, peak_spread)

        hit = not hits.empty or not channel_hits.empty
        status = "数据不足" if not window.empty and window["最高温度"].dropna().empty else ("命中" if hit else "未命中")
        rows.append({
            "机组": unit,
            "故障日期": date.date().isoformat(),
            "描述": event["fault_description"],
            "关注级": best[1] if best else "",
            "最大同级偏离℃": best[2] if best else None,
            "最大级差变化℃": best[3] if best else None,
            "命中": status,
            "首个事件开始": first.strftime("%Y-%m-%d %H:%M") if not hits.empty else "",
            "提前小时": round(lead, 1) if lead is not None else "",
            "命中等级": "、".join(sorted(set(pd.concat([hits["温度阶段"]]).astype(str)))) if not hits.empty else ("通道异常" if not channel_hits.empty else ""),
            "通道异常命中": "是" if not channel_hits.empty else "否",
            "通道异常类型": "、".join(sorted(set(channel_hits["异常类型"].astype(str)))) if not channel_hits.empty else "",
            "通道异常开始": _closest_start(channel_hits, end),
            "同日命中": "是" if starts and first.date() == date.date() else "否",
        })

    report = pd.DataFrame(rows)
    print(report.to_string(index=False))
    report.to_csv(directory / "回测_温度事件.csv", index=False, encoding="utf-8-sig")

    # ---- normal-period distribution ----
    event_mask = pd.Series(False, index=features.index)
    for _, event in events.iterrows():
        start = event["event_date"] - timedelta(days=3)
        end = event["event_date"] + timedelta(days=3)
        event_mask |= (features["unit"] == event["unit_name"]) & (features["dt"] >= start) & (features["dt"] <= end)
    normal = features[~event_mask]
    print("\n正常运行期分位数（排除故障日±3天）：")
    for column in ("级差变化", "最大同级偏离", "最大同级偏低", "最大稳健z", "最高温度"):
        series = normal[column].dropna()
        if series.empty:
            continue
        print(f"  {column:>6}: "
              f"P50={series.quantile(0.5):6.2f}  P90={series.quantile(0.9):6.2f}  "
              f"P99={series.quantile(0.99):6.2f}  P99.5={series.quantile(0.995):6.2f}  "
              f"P99.9={series.quantile(0.999):6.2f}  max={series.max():6.2f}")

    # ---- threshold sweep ----
    print("\n阈值扫描（事件召回 vs 正常期小时触发率）：")
    pre = []
    for _, event in events.iterrows():
        start = event["event_date"] - timedelta(days=LOOKBACK_DAYS)
        end = event["event_date"] + timedelta(days=1)
        part = features[(features["unit"] == event["unit_name"]) & (features["dt"] >= start) & (features["dt"] <= end)]
        pre.append(part)
    print("  peer\\spread | " + " ".join(f"{spread:>6}" for spread in (2, 3, 4, 5, 6, 8)))
    for peer in (2, 3, 4, 5, 6, 8, 10):
        cells = []
        for spread in (2, 3, 4, 5, 6, 8):
            hit_count = 0
            for part in pre:
                mask = (part["最大同级偏离"] >= peer) | (part["级差变化"] >= spread)
                if mask.any():
                    hit_count += 1
            recall = hit_count / len(pre) if pre else 0
            normal_mask = (normal["最大同级偏离"] >= peer) | (normal["级差变化"] >= spread)
            normal_rate = normal_mask.mean() if not normal.empty else 0
            cells.append(f"{recall:.2f}/{normal_rate*100:.2f}%")
        print(f"  {peer:>4}        | " + " ".join(f"{cell:>7}" for cell in cells))
    print("  单元格 = 事件召回 / 正常小时触发率（百分比）")


if __name__ == "__main__":
    main()