Compressor/backend/indicator_features/local_evolution_analysis.py

#!/usr/bin/env python3
"""Find directional 10-degree local feature evolution around repair events."""
from __future__ import annotations

import argparse
import csv
import math
import sys
from collections import defaultdict
from datetime import date, datetime
from pathlib import Path

import numpy as np

BACKEND = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(BACKEND))

from indicator_features.extract_from_db import (  # noqa: E402
    connect, fetch_cycles, fetch_file_rows, point_name, sample_per_day,
)
from indicator_features.indicator_features import _plateau_levels, smooth  # noqa: E402

CYLINDERS = ("一", "二", "三", "四", "五", "六")
SIDES = ("盖侧", "轴侧")
METRICS = ("p50", "slope", "std", "iqr", "d_rms", "d2_abs_max")
METRIC_CN = {
    "p50": "局部中位压力", "slope": "局部斜率", "std": "标准差",
    "iqr": "四分位距", "d_rms": "导数RMS", "d2_abs_max": "最大二阶导",
}
FAMILY = {"p50": "位置", "slope": "斜率", "std": "波动", "iqr": "波动",
          "d_rms": "导数", "d2_abs_max": "导数"}
BAD = {
    "4": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
    "6": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
    "7": {("四", "盖侧"), ("六", "盖侧")},
}
EVENTS = [
    {
        "id": "7号_2025维修过程", "unit": "7", "label": "7号2025停机维修过程",
        "start": "2025-04-29", "end": "2026-02-17", "event": "2025-08-08",
        "pre": [("2025-04-29", "2025-05-04"), ("2025-06-11", "2025-06-14"), ("2025-07-25", "2025-07-28")],
        "pre_segments": [("2025-04-29", "2025-05-04"), ("2025-06-11", "2025-06-14"), ("2025-07-25", "2025-07-28")],
        "early": [("2025-04-29", "2025-05-04")], "late": [("2025-07-25", "2025-07-28")],
        "post_base": [("2025-11-19", "2025-11-21")],
        "post_check": [("2026-02-10", "2026-02-17")],
        "mode": "跨运行段", "markers": "2025-08-08;2025-10-07",
    },
    {
        "id": "7号_2026-05-21", "unit": "7", "label": "7号2026-05-21气阀事件",
        "start": "2026-02-10", "end": "2026-06-07", "event": "2026-05-21",
        "pre": [("2026-05-06", "2026-05-21")],
        "pre_segments": [("2026-02-10", "2026-02-17"), ("2026-03-12", "2026-04-05"), ("2026-05-06", "2026-05-21")],
        "early": [("2026-05-06", "2026-05-09")], "late": [("2026-05-18", "2026-05-21")],
        "post_base": [("2026-05-31", "2026-06-04")],
        "post_check": [("2026-06-04", "2026-06-07")],
        "mode": "连续运行段", "markers": "2026-05-21",
    },
    {
        "id": "4号_2026-06-10", "unit": "4", "label": "4号2026-06-10气阀事件",
        "start": "2026-04-21", "end": "2026-07-31", "event": "2026-06-10",
        "pre": [("2026-05-29", "2026-05-31"), ("2026-06-06", "2026-06-10")],
        "pre_segments": [("2026-04-21", "2026-05-02"), ("2026-05-29", "2026-05-31"), ("2026-06-06", "2026-06-10")],
        "early": [("2026-06-06", "2026-06-08")], "late": [("2026-06-08", "2026-06-10")],
        "post_base": [("2026-06-11", "2026-06-15")],
        "post_check": [("2026-06-21", "2026-06-22"), ("2026-06-30", "2026-07-08"),
                       ("2026-07-11", "2026-07-14"), ("2026-07-22", "2026-07-31")],
        "mode": "短运行段", "markers": "2026-06-10",
    },
    {
        "id": "6号_2026-04-22_26", "unit": "6", "label": "6号2026-04-22~26气阀维修过程",
        "start": "2026-04-12", "end": "2026-08-01", "event": "2026-04-22",
        "pre": [("2026-04-12", "2026-04-22")],
        "pre_segments": [("2026-04-12", "2026-04-22")],
        "early": [("2026-04-12", "2026-04-15")], "late": [("2026-04-19", "2026-04-22")],
        "post_base": [("2026-05-10", "2026-05-17")],
        "post_check": [("2026-05-17", "2026-05-31"), ("2026-06-21", "2026-07-02"),
                       ("2026-07-24", "2026-08-01")],
        "mode": "连续运行段", "markers": "2026-04-22;2026-04-26",
    },
]


def in_ranges(day: str, ranges) -> bool:
    return any(start <= day < end for start, end in ranges)


def robust_scale(values) -> float:
    values = np.asarray(values, dtype=float)
    if values.size == 0:
        return 0.0
    scale = float(np.median(np.abs(values - np.median(values))) * 1.4826)
    return scale if scale > 1e-12 else float(np.std(values))


def rankdata(values):
    values = np.asarray(values)
    order = np.argsort(values, kind="mergesort")
    ranks = np.empty(len(values), dtype=float)
    start = 0
    while start < len(values):
        end = start + 1
        while end < len(values) and values[order[end]] == values[order[start]]:
            end += 1
        ranks[order[start:end]] = (start + end - 1) / 2
        start = end
    return ranks


def spearman(values) -> float:
    if len(values) < 3:
        return float("nan")
    x = np.arange(len(values), dtype=float)
    y = rankdata(values)
    if np.std(y) <= 1e-12:
        return 0.0
    return float(np.corrcoef(x, y)[0, 1])


def normalized_bins(signal):
    signal = np.asarray(signal, dtype=float)
    signal = signal[np.isfinite(signal)]
    if len(signal) < 60:
        return None
    pressure = smooth(signal, 21)
    suction, discharge, _ = _plateau_levels(pressure)
    span = discharge - suction
    if not np.isfinite(span) or span < .3:
        return None
    normalized = (pressure - suction) / span
    angles = np.arange(len(normalized), dtype=float) * 360 / len(normalized)
    result = []
    for index in range(36):
        mask = (angles >= index * 10) & (angles < (index + 1) * 10)
        x, y = angles[mask], normalized[mask]
        if len(y) < 5:
            return None
        slope = float(np.polyfit(x, y, 1)[0])
        dy = np.diff(y)
        d2 = np.diff(dy)
        result.append({
            "bin": index, "bin_start": index * 10,
            "p50": float(np.median(y)), "slope": slope, "std": float(np.std(y)),
            "iqr": float(np.percentile(y, 75) - np.percentile(y, 25)),
            "d_rms": float(np.sqrt(np.mean(dy ** 2))),
            "d2_abs_max": float(np.max(np.abs(d2))) if len(d2) else 0.0,
        })
    return result


def extract(output: Path, per_day: int):
    grouped = defaultdict(list)
    connection = connect()
    fields = ["事件", "机组", "测点", "缸", "腔", "文件ID", "采样时间", "日期", "角度箱"] + list(METRICS)
    output.parent.mkdir(parents=True, exist_ok=True)
    with output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=fields)
        writer.writeheader()
        try:
            for event in EVENTS:
                periods = 0
                for cylinder in CYLINDERS:
                    for side in SIDES:
                        if (cylinder, side) in BAD[event["unit"]]:
                            continue
                        files = fetch_file_rows(connection, point_name(event["unit"], cylinder, side), event["start"], event["end"])
                        picked = sample_per_day(files, per_day)
                        cycles = fetch_cycles(connection, [int(row["id"]) for row in picked], 100) if picked else {}
                        for row in picked:
                            features = normalized_bins(cycles.get(int(row["id"]), np.array([])))
                            if features is None:
                                continue
                            periods += 1
                            day = str(row["sample_time"])[:10]
                            point = f"{cylinder}缸{side}"
                            for item in features:
                                record = {
                                    "事件": event["id"], "机组": event["unit"], "测点": point,
                                    "缸": cylinder, "腔": side, "文件ID": int(row["id"]),
                                    "采样时间": str(row["sample_time"]), "日期": day,
                                    "角度箱": item["bin_start"], **{metric: item[metric] for metric in METRICS},
                                }
                                writer.writerow(record)
                                for metric in METRICS:
                                    grouped[(event["id"], point, day, item["bin_start"], metric)].append(item[metric])
                print(f"{event['label']}: 有效周期 {periods:,}", flush=True)
        finally:
            connection.close()
    return grouped


def load_extracted(path: Path):
    grouped = defaultdict(list)
    with path.open(encoding="utf-8-sig", newline="") as handle:
        for row in csv.DictReader(handle):
            for metric in METRICS:
                grouped[(row["事件"], row["测点"], row["日期"], int(row["角度箱"]), metric)].append(float(row[metric]))
    return grouped


def baseline_stats(series, ranges):
    daily_values = [daily for day, daily, _ in series if in_ranges(day, ranges)]
    cycle_values = [value for day, _, values in series if in_ranges(day, ranges) for value in values]
    center = float(np.median(daily_values)) if daily_values else float("nan")
    # Keep dates equal-weighted for the center, but estimate natural variation
    # from cycles so a 2-3 day baseline cannot produce near-zero scales.
    scale = robust_scale(cycle_values)
    return center, scale, len(daily_values), len(cycle_values)


def merge_bands(candidates):
    grouped = defaultdict(list)
    for row in candidates:
        if row["可信等级"] == "低":
            continue
        start = int(row["角度范围"].split("-")[0])
        grouped[(row["事件"], row["测点"], row["故障主方向"], row["可信等级"])].append((start, row))
    bands = []
    for (event, point, direction, confidence), values in grouped.items():
        values.sort()
        runs = []
        for start, row in values:
            if runs and start == runs[-1][-1][0] + 10:
                runs[-1].append((start, row))
            else:
                runs.append([(start, row)])
        # Merge the circular 350-360 and 0-10 runs.
        if len(runs) > 1 and runs[0][0][0] == 0 and runs[-1][-1][0] == 350:
            runs[0] = runs[-1] + runs[0]
            runs.pop()
        for run in runs:
            starts = [item[0] for item in run]
            rows = [item[1] for item in run]
            if starts[0] > starts[-1]:
                band = f"{starts[0]}-360°及0-{starts[-1] + 10}°"
            else:
                band = f"{starts[0]}-{starts[-1] + 10}°"
            bands.append({
                "事件": event, "测点": point, "角度带": band, "故障主方向": direction,
                "包含10度箱数": len(rows), "可信等级": confidence,
                "最高评分": max(row["候选评分"] for row in rows),
                "演化指标数": max(row["演化指标数"] for row in rows),
                "长期演化指标数": max(row["长期演化指标数"] for row in rows),
                "恢复指标数": max(row["恢复指标数"] for row in rows),
                "长期稳定恢复指标数": max(row["长期稳定恢复指标数"] for row in rows),
                "长期恢复比例平均": max((row["长期恢复比例平均"] for row in rows
                                         if row["长期恢复比例平均"] != ""), default=""),
                "局部幅度平均Z": max(row["局部幅度平均Z"] for row in rows),
                "局部幅度综合RMS": max(row["局部幅度综合RMS"] for row in rows),
                "最早持续日期": min((row["最早持续日期"] for row in rows if row["最早持续日期"]), default=""),
                "代表指标": max(rows, key=lambda row: row["候选评分"])["指标摘要"],
            })
    bands.sort(key=lambda row: (row["事件"], -row["最高评分"], -row["包含10度箱数"]))
    return bands


def analyze(grouped, daily_output: Path, candidate_output: Path, repeat_output: Path, band_output: Path,
            segment_output: Path):
    daily_raw = {key: float(np.median(values)) for key, values in grouped.items()}
    counts = {key: len(values) for key, values in grouped.items()}
    event_by_id = {event["id"]: event for event in EVENTS}
    points_by_event = defaultdict(set)
    for event_id, point, _, _, _ in daily_raw:
        points_by_event[event_id].add(point)
    series_index = defaultdict(list)
    for (event_id, point, day, angle, metric), values in grouped.items():
        series_index[(event_id, point, angle, metric)].append((day, daily_raw[(event_id, point, day, angle, metric)], values))

    stats = {}
    for event in EVENTS:
        for point in points_by_event[event["id"]]:
            for metric in METRICS:
                scales = []
                provisional = {}
                for angle in range(0, 360, 10):
                    series = series_index[(event["id"], point, angle, metric)]
                    post = baseline_stats(series, event["post_base"])
                    early = baseline_stats(series, event["early"])
                    provisional[angle] = (post, early)
                    scales.extend([post[1], early[1]])
                positive = [scale for scale in scales if np.isfinite(scale) and scale > 1e-12]
                floor = .2 * float(np.median(positive)) if positive else 1e-9
                for angle, (post, early) in provisional.items():
                    stats[(event["id"], point, angle, metric)] = {
                        "post_center": post[0], "post_scale": max(post[1], floor), "post_days": post[2],
                        "post_cycles": post[3],
                        "early_center": early[0], "early_scale": max(early[1], floor), "early_days": early[2],
                        "early_cycles": early[3],
                    }

    daily_rows = []
    for key, value in daily_raw.items():
        event_id, point, day, angle, metric = key
        stat = stats[(event_id, point, angle, metric)]
        z_post = (value - stat["post_center"]) / stat["post_scale"] if np.isfinite(stat["post_center"]) else float("nan")
        z_early = (value - stat["early_center"]) / stat["early_scale"] if np.isfinite(stat["early_center"]) else float("nan")
        if not (np.isfinite(z_post) and np.isfinite(z_early)):
            continue
        daily_rows.append({
            "事件": event_id, "测点": point, "日期": day, "角度箱": angle,
            "指标": metric, "指标中文": METRIC_CN[metric], "特征族": FAMILY[metric],
            "每日中位": value, "当日周期数": counts[key],
            "维修后基线中心": stat["post_center"], "维修后基线尺度": stat["post_scale"],
            "维修后基线日数": stat["post_days"], "相对维修后Z": z_post,
            "维修后基线周期数": stat["post_cycles"],
            "早期基线中心": stat["early_center"], "早期基线尺度": stat["early_scale"],
            "早期基线日数": stat["early_days"], "早期基线周期数": stat["early_cycles"],
            "相对早期Z": z_early,
        })
    common = defaultdict(list)
    for row in daily_rows:
        if np.isfinite(row["相对维修后Z"]):
            common[(row["事件"], row["日期"], row["角度箱"], row["指标"])].append(row["相对维修后Z"])
    common = {key: float(np.median(values)) for key, values in common.items()}
    for row in daily_rows:
        public = common.get((row["事件"], row["日期"], row["角度箱"], row["指标"]), float("nan"))
        row["全机公共Z"] = public
        row["局部残差Z"] = row["相对维修后Z"] - public

    segment_rows = []
    for event in EVENTS:
        ranges = [("故障前运行段", start, end) for start, end in event["pre_segments"]]
        ranges += [("维修后基线", start, end) for start, end in event["post_base"]]
        ranges += [("维修后长期验证", start, end) for start, end in event["post_check"]]
        event_rows = [row for row in daily_rows if row["事件"] == event["id"]]
        keys = sorted({(row["测点"], row["角度箱"], row["指标"], row["指标中文"], row["特征族"])
                       for row in event_rows})
        for segment_type, start, end in ranges:
            for point, angle, metric, metric_cn, family in keys:
                selected = [row for row in event_rows if row["测点"] == point and row["角度箱"] == angle
                            and row["指标"] == metric and start <= row["日期"] < end]
                if not selected:
                    continue
                segment_rows.append({
                    "事件": event["id"], "阶段类型": segment_type, "阶段开始": start, "阶段结束": end,
                    "测点": point, "角度箱": angle, "指标": metric, "指标中文": metric_cn, "特征族": family,
                    "有效日期数": len({row["日期"] for row in selected}),
                    "周期数": sum(int(row["当日周期数"]) for row in selected),
                    "阶段每日中位的中位": float(np.median([row["每日中位"] for row in selected])),
                    "相对维修后Z中位": float(np.median([row["相对维修后Z"] for row in selected])),
                    "全机公共Z中位": float(np.median([row["全机公共Z"] for row in selected])),
                    "局部残差Z中位": float(np.median([row["局部残差Z"] for row in selected])),
                })

    daily_output.parent.mkdir(parents=True, exist_ok=True)
    with daily_output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(daily_rows[0]))
        writer.writeheader(); writer.writerows(daily_rows)
    with segment_output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(segment_rows[0]))
        writer.writeheader(); writer.writerows(segment_rows)

    rows_index = defaultdict(list)
    for row in daily_rows:
        rows_index[(row["事件"], row["测点"], row["角度箱"], row["指标"])].append(row)
    metric_candidates = []
    for key, rows in rows_index.items():
        event_id, point, angle, metric = key
        event = event_by_id[event_id]
        rows.sort(key=lambda row: row["日期"])
        pre = [row for row in rows if in_ranges(row["日期"], event["pre"])]
        late = [row for row in rows if in_ranges(row["日期"], event["late"])]
        check = [row for row in rows if in_ranges(row["日期"], event["post_check"])]
        if len(pre) < 3 or not late:
            continue
        late_z = float(np.median([row["相对维修后Z"] for row in late]))
        late_residual = float(np.median([row["局部残差Z"] for row in late]))
        direction = 1 if late_z >= 0 else -1
        directed_late_residual = direction * late_residual
        directed = [direction * row["相对维修后Z"] for row in pre]
        directed_local = [direction * row["局部残差Z"] for row in pre]
        rho = spearman(directed)
        last_days = directed[-min(5, len(directed)):]
        persistent = sum(value >= 3 for value in last_days)
        local_persistent = sum(value >= 2 for value in directed_local[-min(5, len(directed_local)):])
        late_early_z = float(np.median([direction * row["相对早期Z"] for row in late]))
        check_z = float(np.median([abs(row["相对维修后Z"]) for row in check])) if check else float("nan")
        check_segments = []
        for start, end in event["post_check"]:
            segment_check = [row for row in rows if start <= row["日期"] < end]
            if segment_check:
                check_segments.append(float(np.median([abs(row["相对维修后Z"]) for row in segment_check])))
        recovery_flags = [value <= 2 or value <= .5 * abs(late_z) for value in check_segments]
        recovered = bool(recovery_flags and recovery_flags[0])
        long_recovery_ratio = float(np.mean(recovery_flags)) if recovery_flags else float("nan")
        long_recovery_stable = bool(recovery_flags and all(recovery_flags))
        segment_values = []
        segment_local_values = []
        segment_labels = []
        for start, end in event["pre_segments"]:
            segment_rows = [row for row in rows if start <= row["日期"] < end]
            if not segment_rows:
                continue
            segment_values.append(float(np.median([direction * row["相对维修后Z"] for row in segment_rows])))
            segment_local_values.append(float(np.median([direction * row["局部残差Z"] for row in segment_rows])))
            segment_labels.append(f"{start}~{end}")
        long_rho = spearman(segment_local_values) if len(segment_local_values) >= 3 else float("nan")
        long_change = segment_local_values[-1] - segment_local_values[0] if len(segment_local_values) >= 2 else float("nan")
        long_evolution = len(segment_local_values) >= 2 and long_change >= 2 and (
            len(segment_local_values) == 2 or long_rho >= .5
        )
        onset = ""
        for index in range(len(pre)):
            window = directed[index:index + 3]
            if len(window) >= 2 and sum(value >= 3 for value in window) >= 2:
                onset = pre[index]["日期"]
                break
        support = abs(late_z) >= 3 and directed_late_residual >= 2
        evolution = (persistent >= 2 and rho >= .35) or late_early_z >= 2
        score = min(abs(late_residual), 12) + min(abs(late_z), 12) * .5
        score += 3 if evolution else 0
        score += 3 if recovered else 0
        score += min(persistent, 3) + min(local_persistent, 3)
        if event["mode"] != "连续运行段":
            score -= 2
        if min(rows[0]["维修后基线日数"], rows[0]["早期基线日数"]) < 3:
            score -= 2
        metric_candidates.append({
            "事件": event_id, "模式": event["mode"], "测点": point, "角度起点": angle,
            "指标": metric, "指标中文": METRIC_CN[metric], "特征族": FAMILY[metric],
            "故障方向": "升高" if direction > 0 else "降低", "维修前末段Z": late_z,
            "维修前末段局部残差Z": late_residual, "定向局部幅度Z": directed_late_residual,
            "末段相对早期定向Z": late_early_z,
            "故障前Spearman": rho, "末5日超3Z天数": persistent,
            "末5日局部超2Z天数": local_persistent, "首次持续日期": onset,
            "独立维修后验证Z": check_z, "维修后恢复": int(recovered),
            "后期验证运行段数": len(check_segments),
            "后期各段绝对Z": ";".join(f"{value:.3f}" for value in check_segments),
            "长期恢复比例": long_recovery_ratio,
            "长期恢复稳定": int(long_recovery_stable),
            "基础支持": int(support), "演化支持": int(evolution), "单指标评分": score,
            "长期运行段数": len(segment_local_values), "长期段间Spearman": long_rho,
            "长期段间局部变化Z": long_change, "长期演化支持": int(long_evolution),
            "长期各段定向Z": ";".join(f"{label}:{value:.3f}" for label, value in zip(segment_labels, segment_local_values)),
            "维修后基线日数": rows[0]["维修后基线日数"], "早期基线日数": rows[0]["早期基线日数"],
        })

    groups = defaultdict(list)
    for row in metric_candidates:
        if row["基础支持"]:
            groups[(row["事件"], row["测点"], row["角度起点"])].append(row)
    candidates = []
    for (event_id, point, angle), items in groups.items():
        families = {item["特征族"] for item in items}
        directions = defaultdict(list)
        for item in items:
            directions[item["故障方向"]].append(item)
        direction, aligned = max(directions.items(), key=lambda pair: len(pair[1]))
        evolution_n = sum(item["演化支持"] for item in items)
        long_evolution_n = sum(item["长期演化支持"] for item in items)
        recovery_n = sum(item["维修后恢复"] for item in items)
        stable_recovery_n = sum(item["长期恢复稳定"] for item in items)
        recovery_ratios = [item["长期恢复比例"] for item in items if np.isfinite(item["长期恢复比例"])]
        directed_amplitudes = [item["定向局部幅度Z"] for item in items]
        amplitude_mean = float(np.mean(directed_amplitudes))
        amplitude_rms = float(np.sqrt(np.mean(np.square(directed_amplitudes))))
        amplitude_median = float(np.median(directed_amplitudes))
        score = sum(max(0, item["单指标评分"]) for item in items) + 2 * len(families) + 2 * evolution_n + 2 * recovery_n
        score += 2 * long_evolution_n + min(amplitude_mean, 8) + min(amplitude_rms, 10)
        candidates.append({
            "事件": event_id, "模式": items[0]["模式"], "测点": point,
            "角度范围": f"{angle}-{angle + 10}°", "故障主方向": direction,
            "支持指标数": len(items), "特征族数": len(families),
            "演化指标数": evolution_n, "恢复指标数": recovery_n,
            "长期稳定恢复指标数": stable_recovery_n,
            "长期恢复比例平均": round(float(np.mean(recovery_ratios)), 3) if recovery_ratios else "",
            "长期演化指标数": long_evolution_n,
            "局部幅度平均Z": round(amplitude_mean, 3),
            "局部幅度中位Z": round(amplitude_median, 3),
            "局部幅度综合RMS": round(amplitude_rms, 3),
            "候选评分": round(score, 3),
            "最早持续日期": min((item["首次持续日期"] for item in items if item["首次持续日期"]), default=""),
            "指标摘要": ";".join(
                f"{item['指标中文']}{'↑' if item['故障方向']=='升高' else '↓'}"
                f"(局部{item['维修前末段局部残差Z']:.1f}Z,ρ={item['故障前Spearman']:.2f},恢复={item['维修后恢复']})"
                for item in sorted(items, key=lambda item: -item["单指标评分"])
            ),
            "可信等级": "高" if items[0]["模式"] == "连续运行段" and len(families) >= 2 and evolution_n >= 2 and recovery_n >= 2 and amplitude_mean >= 2
                         else "中" if items[0]["模式"] == "连续运行段" and len(families) >= 2 and evolution_n >= 2 and recovery_n >= 1
                         else "短期线索" if items[0]["模式"] == "短运行段" and len(families) >= 2 and evolution_n >= 2 and recovery_n >= 1 and amplitude_mean >= 2
                         else "跨段线索" if items[0]["模式"] == "跨运行段" and len(families) >= 2 and long_evolution_n >= 2 and amplitude_mean >= 2
                         else "低",
        })
    candidates.sort(key=lambda row: (row["事件"], -row["候选评分"]))
    with candidate_output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(candidates[0]))
        writer.writeheader(); writer.writerows(candidates)

    bands = merge_bands(candidates)
    with band_output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(bands[0]))
        writer.writeheader(); writer.writerows(bands)

    repeats = []
    by_band = defaultdict(list)
    for row in candidates:
        if row["可信等级"] in ("高", "中", "短期线索", "跨段线索"):
            by_band[row["角度范围"]].append(row)
    for band, items in by_band.items():
        events = sorted({item["事件"] for item in items})
        if len(events) < 2:
            continue
        repeats.append({
            "角度范围": band, "事件数": len(events), "事件": ";".join(events),
            "测点": ";".join(sorted({item["测点"] for item in items})),
            "最高评分": max(item["候选评分"] for item in items),
            "共同特征": ";".join(item["指标摘要"] for item in sorted(items, key=lambda row: -row["候选评分"])[:4]),
        })
    repeats.sort(key=lambda row: (-row["事件数"], -row["最高评分"]))
    with repeat_output.open("w", encoding="utf-8-sig", newline="") as handle:
        fields = ["角度范围", "事件数", "事件", "测点", "最高评分", "共同特征"]
        writer = csv.DictWriter(handle, fieldnames=fields); writer.writeheader(); writer.writerows(repeats)
    return metric_candidates, candidates, repeats, bands


def write_report(path: Path, metric_candidates, candidates, repeats, bands):
    lines = ["# 局部10度特征演化分析", "", "## 方法", "",
             "- 不使用整周期欧氏距离。每测点每天均匀抽取最多12个周期。",
             "- 每周期先平滑并按吸排气平台归一化，再在原始键相角度上划分36个10度箱。",
             "- 核心特征：局部中位压力、斜率、标准差、四分位距、导数RMS、最大二阶导。",
             "- 同时使用运行段早期基线和维修后稳定基线；逐日结果扣除全机同角度、同指标的公共变化。",
             "- 候选强调故障前定向演化、局部残差和独立维修后日期的反向恢复。",
             "- 幅度评估不增加原始特征：对支持指标的定向局部残差计算平均值、中位数和综合RMS。",
             "- 恢复分为事件后首个验证段的短期恢复，以及所有后续运行段都满足的长期稳定恢复。", ""]
    for event in EVENTS:
        rows = [row for row in candidates if row["事件"] == event["id"]]
        high = [row for row in rows if row["可信等级"] == "高"]
        medium = [row for row in rows if row["可信等级"] in ("中", "短期线索", "跨段线索")]
        lines += [f"## {event['label']}", "", f"- 分析模式：{event['mode']}。", f"- 高可信候选：{len(high)}；其他可复核线索：{len(medium)}。"]
        if event["id"] == "7号_2025维修过程":
            lines.append("- 维修后基线仅2个日期且中间长期缺测，只能解释为跨段状态线索，不能证明连续故障前演化。")
        for index, row in enumerate((high + medium)[:12], 1):
            lines.append(
                f"{index}. {row['测点']} {row['角度范围']}，{row['故障主方向']}，评分{row['候选评分']:.1f}，"
                f"{row['可信等级']}可信；幅度平均{row['局部幅度平均Z']:.1f}Z，综合RMS{row['局部幅度综合RMS']:.1f}，"
                f"近期演化{row['演化指标数']}项，长期跨段演化{row['长期演化指标数']}项，"
                f"长期稳定恢复{row['长期稳定恢复指标数']}项；{row['指标摘要']}"
            )
        event_bands = [row for row in bands if row["事件"] == event["id"]]
        if event_bands:
            lines.append("- 相邻角度合并后的主要变化带：" + "；".join(
                f"{row['测点']} {row['角度带']}({row['故障主方向']}，平均{row['局部幅度平均Z']:.1f}Z，RMS{row['局部幅度综合RMS']:.1f})"
                for row in event_bands[:10]
            ) + "。")
        if not high and not medium:
            lines.append("- 未找到同时满足局部残差、演化或维修后恢复条件的中高可信候选。")
        lines.append("")
    lines += ["## 跨事件重复角度", ""]
    if repeats:
        for row in repeats[:20]:
            lines.append(f"- {row['角度范围']}：{row['事件数']}个事件；{row['事件']}；涉及{row['测点']}。")
    else:
        lines.append("- 暂无完全相同10度箱在至少两个事件中形成中高可信候选。")
    lines += ["", "## 限制", "", "- 角度为键相原始坐标，没有自动相位旋转；不同测点的物理过程区间可能不同。",
              "- 基线日期较少时稳健尺度仍可能不稳定，结果保留原始方向、日数和独立恢复验证供后续复核。",
              "- 当前候选是维修相关局部特征，不自动等于通用气阀故障规则。"]
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")


def main():
    parser = argparse.ArgumentParser(description="10度局部特征双基线演化分析")
    parser.add_argument("--out-dir", default="cache/indicator/局部演化分析")
    parser.add_argument("--per-day", type=int, default=12)
    parser.add_argument("--reuse", action="store_true", help="复用已抽取的逐周期CSV")
    args = parser.parse_args()
    out = Path(args.out_dir); out.mkdir(parents=True, exist_ok=True)
    cycle_path = out / "局部特征_逐周期.csv"
    grouped = load_extracted(cycle_path) if args.reuse and cycle_path.exists() else extract(cycle_path, args.per_day)
    metric, candidates, repeats, bands = analyze(
        grouped, out / "局部特征_逐日.csv", out / "局部演化候选.csv", out / "跨事件重复角度.csv",
        out / "局部演化角度带.csv", out / "局部特征_逐运行段.csv",
    )
    write_report(out / "局部演化结论.md", metric, candidates, repeats, bands)
    print(f"逐日组合 {len(grouped):,}，单指标候选 {len(metric):,}，角度候选 {len(candidates):,}，角度带 {len(bands):,}，重复角度 {len(repeats):,}")
    print(f"结果 -> {out}")


if __name__ == "__main__":
    main()