Compressor/backend/indicator_features/pseudo_online_detection.py

#!/usr/bin/env python3
"""Pseudo-online local-feature detection with past-only baselines."""
from __future__ import annotations

import csv
import argparse
import math
import sys
from collections import defaultdict, deque
from datetime import date, datetime, timedelta
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
    CYLINDERS, SIDES, connect, fetch_cycles, fetch_file_rows, point_name, sample_per_day,
)
from indicator_features.local_evolution_analysis import (  # noqa: E402
    FAMILY, METRICS, METRIC_CN, normalized_bins, robust_scale,
)

BAD = {
    "4": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
    "6": {("二", "盖侧"), ("四", "盖侧"), ("六", "盖侧")},
    "7": {("四", "盖侧"), ("六", "盖侧")},
}
FAILURE_TEMPLATES = {
    "7": {
        ("一缸盖侧",240,1),("一缸盖侧",250,1),("一缸盖侧",260,1),
        ("二缸盖侧",140,-1),("三缸轴侧",160,1),("三缸轴侧",170,1),
        ("三缸盖侧",20,1),("三缸盖侧",30,1),("三缸盖侧",310,-1),("三缸盖侧",320,-1),
    },
    "4": {("四缸轴侧",60,-1)},
    "6": {
        ("三缸轴侧",160,1),("三缸轴侧",170,1),("三缸轴侧",120,1),("三缸轴侧",130,1),
        ("一缸轴侧",290,1),("一缸盖侧",220,-1),("三缸轴侧",50,1),
    },
}
EVENTS = [
    {"id": "7号_2025维修过程", "unit": "7", "event": "2025-08-08",
     "baseline": ("2025-04-29", "2025-06-14"), "detect": ("2025-07-25", "2025-07-28"),
     "note": "事件附近缺测，仅验证7月底运行段状态"},
    {"id": "7号_2026-05-21", "unit": "7", "event": "2026-05-21",
     "baseline": ("2026-02-10", "2026-04-05"), "detect": ("2026-05-06", "2026-05-22"), "note": ""},
    {"id": "4号_2026-06-10", "unit": "4", "event": "2026-06-10",
     "baseline": ("2026-04-21", "2026-05-31"), "detect": ("2026-06-06", "2026-06-11"), "note": ""},
    {"id": "6号_2026-04-22_26", "unit": "6", "event": "2026-04-22",
     "baseline": ("2026-04-12", "2026-04-15"), "detect": ("2026-04-15", "2026-04-27"),
     "note": "历史不足，仅3个基线日期"},
]
UNITS = ("4", "6", "7", "8", "9")
MERGE_GAP = 2
HEALTH_SEGMENTS = {
    "7": [("2026-03-12", "2026-04-05")],
    "8": [("2026-04-18", "2026-05-07")],
    "9": [("2026-04-01", "2026-04-19")],
}


def calendar_segments(days):
    segments=[]
    for day in sorted(days):
        if segments and (day-segments[-1][1]).days <= MERGE_GAP+1:
            segments[-1][1]=day
        else:
            segments.append([day,day])
    return [(start,end) for start,end in segments]


def date_ranges(event):
    start,end=event["detect"]
    cursor=date.fromisoformat(start);stop=date.fromisoformat(end)
    while cursor<stop:
        yield cursor.isoformat()
        cursor+=timedelta(days=1)


def process_unit(connection, unit, start, end, per_day_baseline=3):
    points={}
    running_days=set()
    for cylinder in CYLINDERS:
        for side in SIDES:
            if (cylinder,side) in BAD.get(unit,set()):continue
            files=fetch_file_rows(connection,point_name(unit,cylinder,side),start,end)
            running_days.update(row["sample_time"].date() for row in files if row.get("rpm") and row["rpm"]>0)
            baseline_files=sample_per_day(files,per_day_baseline)
            cycles=fetch_cycles(connection,[int(row["id"]) for row in baseline_files],100) if baseline_files else {}
            baseline_records=[]
            for row in baseline_files:
                values=stats_for_cycle(cycles.get(int(row["id"]),np.array([])))
                if values is None:continue
                time=row["sample_time"] if isinstance(row["sample_time"],datetime) else datetime.fromisoformat(str(row["sample_time"]))
                baseline_records.append({"time":time,"day":time.date().isoformat(),"values":values})
            if baseline_records:
                points[f"{cylinder}缸{side}"]={"baseline_records":baseline_records,"detection_records":[]}
    return points,calendar_segments(running_days)


def stats_for_cycle(signal):
    bins = normalized_bins(signal)
    if bins is None:
        return None
    return {(item["bin_start"], metric): item[metric] for item in bins for metric in METRICS}


def point_rows(connection, event, cylinder, side, baseline_per_day=3):
    point = point_name(event["unit"], cylinder, side)
    baseline_files = sample_per_day(
        fetch_file_rows(connection, point, event["baseline"][0], event["baseline"][1]), baseline_per_day
    )
    detection_files = fetch_file_rows(connection, point, event["detect"][0], event["detect"][1])
    cycles = fetch_cycles(connection, [int(row["id"]) for row in baseline_files], 100) if baseline_files else {}
    baseline, detection_meta = [], []
    baseline_ids = {int(row["id"]) for row in baseline_files}
    for row in baseline_files + detection_files:
        file_id = int(row["id"])
        time = row["sample_time"] if isinstance(row["sample_time"], datetime) else datetime.fromisoformat(str(row["sample_time"]))
        if file_id in baseline_ids:
            values = stats_for_cycle(cycles.get(file_id, np.array([])))
            if values is None: continue
            baseline.append({"time":time,"day":time.date().isoformat(),"slot":time.hour//6,"values":values})
        else:
            detection_meta.append({"id":file_id,"time":time,"day":time.date().isoformat(),"slot":time.hour//6})
    return baseline, detection_meta


def build_baseline(records):
    per_key_day = defaultdict(list)
    per_key_cycles = defaultdict(list)
    for record in records:
        for key, value in record["values"].items():
            per_key_day[(key, record["day"])].append(value)
            per_key_cycles[key].append(value)
    daily = defaultdict(list)
    for (key, _), values in per_key_day.items():
        daily[key].append(float(np.median(values)))
    provisional = {}
    for key in [(angle, metric) for angle in range(0, 360, 10) for metric in METRICS]:
        days = daily.get(key, [])
        cycles = per_key_cycles.get(key, [])
        provisional[key] = {
            "center": float(np.median(days)) if days else float("nan"),
            "scale": max(robust_scale(days), robust_scale(cycles)),
            "days": len(days), "cycles": len(cycles),
        }
    for metric in METRICS:
        scales = [provisional[(angle, metric)]["scale"] for angle in range(0, 360, 10)
                  if provisional[(angle, metric)]["scale"] > 1e-12]
        floor = .2 * float(np.median(scales)) if scales else 1e-9
        for angle in range(0, 360, 10):
            provisional[(angle, metric)]["scale"] = max(provisional[(angle, metric)]["scale"], floor)
    return provisional


def aggregate(records, group_field):
    grouped = defaultdict(list)
    for record in records:
        group = record[group_field] if group_field != "day_slot" else f"{record['day']}#{record['slot']}"
        for key, value in record["values"].items():
            grouped[(group, key)].append(value)
    return {(group, key): float(np.median(values)) for (group, key), values in grouped.items()}


def merge_bands(alerts):
    groups = defaultdict(list)
    for row in alerts:
        groups[(row["事件"], row["日期"], row["测点"], row["等级"])].append(row)
    output = []
    for key, rows in groups.items():
        rows.sort(key=lambda row: int(row["角度箱"]))
        runs = []
        for row in rows:
            angle = int(row["角度箱"])
            if runs and angle == int(runs[-1][-1]["角度箱"]) + 10:
                runs[-1].append(row)
            else:
                runs.append([row])
        if len(runs) > 1 and int(runs[0][0]["角度箱"]) == 0 and int(runs[-1][-1]["角度箱"]) == 350:
            runs[0] = runs[-1] + runs[0]; runs.pop()
        for run in runs:
            starts = [int(row["角度箱"]) for row in run]
            band = f"{starts[0]}-360°及0-{starts[-1]+10}°" if starts[0] > starts[-1] else f"{starts[0]}-{starts[-1]+10}°"
            best = max(run, key=lambda row: float(row["综合RMS"]))
            output.append({
                "事件": key[0], "日期": key[1], "测点": key[2], "等级": key[3], "角度带": band,
                "角度箱数": len(run), "最大支持特征数": max(int(row["支持特征数"]) for row in run),
                "最大特征族数": max(int(row["支持特征族数"]) for row in run),
                "最大幅度平均Z": max(float(row["幅度平均Z"]) for row in run),
                "最大综合RMS": max(float(row["综合RMS"]) for row in run),
                "最大异常子窗数": max(int(row["异常6小时子窗数"]) for row in run),
                "代表特征": best["特征摘要"],
            })
    return output


def run(out_dir: Path, selected_ids=None, slow_channel=False):
    all_detail, all_bands, validations = [], [], []
    reuse_dir=Path(str(out_dir)+"_reuse")
    reuse_dir.mkdir(parents=True,exist_ok=True)
    connection = connect()
    try:
        for event in EVENTS:
            if selected_ids and event["id"] not in selected_ids:
                continue
            point_data = {}
            for cylinder in CYLINDERS:
                for side in SIDES:
                    if (cylinder, side) in BAD[event["unit"]]:
                        continue
                    baseline, detection_meta = point_rows(connection, event, cylinder, side)
                    if baseline and detection_meta:
                        point_data[f"{cylinder}缸{side}"]={"baseline":build_baseline(baseline),"baseline_records":baseline,
                                                               "detection_meta":detection_meta}
            early_baselines={}
            for point,data in point_data.items():
                early_days=sorted({record["day"] for record in data["baseline_records"]})[:3]
                early_records=[record for record in data["baseline_records"] if record["day"] in early_days]
                early_baselines[point]=build_baseline(early_records) if early_records else {}
            detection_days=sorted({record["day"] for data in point_data.values() for record in data["detection_meta"]})
            history=defaultdict(deque)
            slow_history=defaultdict(deque)
            event_alerts=[]
            for day in detection_days:
                raw_by_point={}
                slot_z=defaultdict(dict)
                for point,data in point_data.items():
                    meta=[record for record in data["detection_meta"] if record["day"]==day]
                    ids=[record["id"] for record in meta]
                    cycles=fetch_cycles(connection,ids,100) if ids else {}
                    records=[]
                    for record in meta:
                        values=stats_for_cycle(cycles.get(record["id"],np.array([])))
                        if values is not None:
                            records.append({**record,"values":values})
                    if not records:continue
                    window_z={}
                    for angle in range(0,360,10):
                        for metric in METRICS:
                            key=(angle,metric);base=data["baseline"][key]
                            if not np.isfinite(base["center"]):continue
                            vals=[rec["values"][key] for rec in records if key in rec["values"]]
                            if not vals:continue
                            window_z[key]=(float(np.median(vals))-base["center"])/base["scale"]
                            for slot in range(4):
                                slot_vals=[rec["values"][key] for rec in records if rec["slot"]==slot and key in rec["values"]]
                                if slot_vals:slot_z[(slot,key)][point]= (float(np.median(slot_vals))-base["center"])/base["scale"]
                    raw_by_point[point]=window_z
                common={}
                for (slot,key),values in slot_z.items():common[(slot,key)]=float(np.median(list(values.values())))
                common24={}
                for key in {(a,m) for point in raw_by_point for a,m in raw_by_point[point]}:
                    vals=[raw_by_point[point][key] for point in raw_by_point if key in raw_by_point[point]]
                    common24[key]=float(np.median(vals)) if vals else 0.0
                by_candidate=defaultdict(list)
                for point,window_z in raw_by_point.items():
                    for (angle,metric),z in window_z.items():
                        row={"point":point,"angle":angle,"metric":metric,"z":z,"local":z-common24[(angle,metric)],"base":point_data[point]["baseline"][(angle,metric)]}
                        by_candidate[(point,angle)].append(row)
                abnormal_bins_by_point=defaultdict(set)
                for (point,angle),rows in by_candidate.items():
                    if any(abs(row["local"])>=6 for row in rows):
                        abnormal_bins_by_point[point].add(angle)
                channel_fault={point for point,bins in abnormal_bins_by_point.items() if len(bins)>=18}
                for (point,angle),rows in sorted(by_candidate.items(),key=lambda item:(item[0][0],item[0][1])):
                    supported = [row for row in rows if abs(row["local"]) >= 2]
                    families = {FAMILY[row["metric"]] for row in supported}
                    amplitudes = [abs(row["local"]) for row in supported]
                    mean_amp = float(np.mean(amplitudes)) if amplitudes else 0.0
                    rms = float(np.sqrt(np.mean(np.square(amplitudes)))) if amplitudes else 0.0
                    abnormal_slots=0
                    for slot in range(4):
                        slot_values=[]
                        for row in rows:
                            slot_key=(slot,(angle,row["metric"]))
                            if point in slot_z.get(slot_key,{}):
                                slot_values.append(abs(slot_z[slot_key][point]-common.get(slot_key,0.0)))
                        if sum(value>=2 for value in slot_values)>=2:abnormal_slots+=1
                    candidate = len(families) >= 2 and mean_amp >= 2 and rms >= 3 and abnormal_slots >= 2
                    hist = history[(point, angle)]; hist.append(candidate)
                    if len(hist) > 5: hist.popleft()
                    persistent = len(hist) >= 3 and sum(hist) >= 3

                    slow_rows=[]
                    for row in rows:
                        metric=row["metric"];key=(angle,metric)
                        early=early_baselines[point].get(key)
                        if not early or not np.isfinite(early["center"]):continue
                        current_value=row["z"]*row["base"]["scale"]+row["base"]["center"]
                        scale=max(early["scale"],row["base"]["scale"]*.2)
                        slow_z=(current_value-early["center"])/scale
                        slow_rows.append({**row,"slow_z":slow_z})
                    slow_public_values=defaultdict(list)
                    for other_point,other_rows in by_candidate.items():
                        other_angle=other_point[1]
                        if other_angle!=angle:continue
                        for other_row in other_rows:
                            other_early=early_baselines[other_point[0]].get((angle,other_row["metric"]))
                            if not other_early or not np.isfinite(other_early["center"]):continue
                            value=other_row["z"]*other_row["base"]["scale"]+other_row["base"]["center"]
                            scale=max(other_early["scale"],other_row["base"]["scale"]*.2)
                            slow_public_values[other_row["metric"]].append((value-other_early["center"])/scale)
                    slow_supported_rows=[]
                    for row in slow_rows:
                        public=float(np.median(slow_public_values[row["metric"]])) if slow_public_values[row["metric"]] else 0.0
                        row["slow_local"]=row["slow_z"]-public
                        if abs(row["slow_local"])>=2:slow_supported_rows.append(row)
                    slow_families={FAMILY[row["metric"]] for row in slow_supported_rows}
                    slow_amplitudes=[abs(row["slow_local"]) for row in slow_supported_rows]
                    slow_mean=float(np.mean(slow_amplitudes)) if slow_amplitudes else 0.0
                    slow_rms=float(np.sqrt(np.mean(np.square(slow_amplitudes)))) if slow_amplitudes else 0.0
                    slow_candidate=len(slow_families)>=2 and slow_mean>=2 and slow_rms>=3
                    slow_hist=slow_history[(point,angle)];slow_hist.append(slow_candidate)
                    if len(slow_hist)>5:slow_hist.popleft()
                    slow_persistent=len(slow_hist)>=3 and sum(slow_hist)>=3
                    raw_supported=supported or slow_supported_rows
                    direction_values=[row["local"] for row in supported] or [row["slow_local"] for row in slow_supported_rows]
                    direction=1 if direction_values and np.median(direction_values)>=0 else -1
                    signature_match=any(tp==point and ta==angle and td==direction
                                        for tp,ta,td in FAILURE_TEMPLATES.get(event["unit"],set()))
                    candidate_level=(candidate or slow_candidate)
                    level="通道质控" if point in channel_fault else "持续预警" if ((persistent and candidate) or (slow_persistent and slow_candidate)) else "候选" if candidate_level else "观察" if raw_supported else "正常"
                    if level in ("候选","持续预警") and not signature_match:level="未分类异常"
                    base_days=min(row["base"]["days"] for row in rows)
                    detail={
                        "事件":event["id"],"机组":event["unit"],"事件日期":event["event"],"日期":day,
                        "测点":point,"角度箱":angle,"等级":level,
                        "支持特征数":len(supported),"支持特征族数":len(families),
                        "幅度平均Z":round(mean_amp,4),"幅度中位Z":round(float(np.median(amplitudes)) if amplitudes else 0,4),
                        "综合RMS":round(rms,4),"异常6小时子窗数":abnormal_slots,
                        "最近5窗候选数":sum(hist),"基线有效日数":base_days,
                        "慢通道支持特征数":len(slow_supported_rows),"慢通道特征族数":len(slow_families),
                        "慢通道幅度平均Z":round(slow_mean,4),"慢通道综合RMS":round(slow_rms,4),
                        "慢通道最近5窗候选数":sum(slow_hist),"慢通道持续":int(slow_persistent),
                        "历史事件签名匹配":int(signature_match),
                        "基线置信":"低" if base_days<7 else "中" if base_days<20 else "高",
                        "特征摘要":";".join(f"{METRIC_CN[row['metric']]}{'↑' if row['local']>0 else '↓'}{abs(row['local']):.1f}Z" for row in sorted(supported,key=lambda x:-abs(x['local']))),
                        "慢通道特征摘要":";".join(f"{METRIC_CN[row['metric']]}{'↑' if row['slow_local']>0 else '↓'}{abs(row['slow_local']):.1f}Z" for row in sorted(slow_supported_rows,key=lambda x:-abs(x['slow_local']))),
                        "备注":event["note"],
                    }
                    all_detail.append(detail)
                    if level not in ("正常","通道质控"):event_alerts.append(detail)
            bands = merge_bands(event_alerts); all_bands.extend(bands)
            event_date = datetime.fromisoformat(event["event"]).date()
            grouped_bands = defaultdict(list)
            for row in bands: grouped_bands[(row["测点"], row["角度带"])].append(row)
            for (point, band), rows in grouped_bands.items():
                rows.sort(key=lambda row: row["日期"])
                candidates = [row for row in rows if row["等级"] in ("候选", "持续预警", "未分类异常")]
                persistent = [row for row in rows if row["等级"] == "持续预警"]
                if not candidates: continue
                first = datetime.fromisoformat(candidates[0]["日期"]).date()
                validations.append({
                    "事件": event["id"], "事件日期": event["event"], "测点": point, "角度带": band,
                    "首次候选日期": str(first), "提前日数": (event_date-first).days,
                    "首次持续预警日期": persistent[0]["日期"] if persistent else "",
                    "候选窗口数": len(candidates), "持续预警窗口数": len(persistent),
                    "未分类异常窗口数": sum(row["等级"] == "未分类异常" for row in candidates),
                    "最大幅度平均Z": max(float(row["最大幅度平均Z"]) for row in candidates),
                    "最大综合RMS": max(float(row["最大综合RMS"]) for row in candidates),
                    "备注": event["note"],
                })
            if all_detail and any(row["事件"]==event["id"] for row in all_detail):
                event_detail=[row for row in all_detail if row["事件"]==event["id"]]
                event_band=[row for row in all_bands if row["事件"]==event["id"]]
                event_validation=[row for row in validations if row["事件"]==event["id"]]
                for filename,records in (("伪在线_窗口检测明细.csv",event_detail),
                                         ("伪在线_角度带告警.csv",event_band),
                                         ("伪在线_事件验证.csv",event_validation)):
                    if records:
                        with (reuse_dir/filename).open("w",encoding="utf-8-sig",newline="") as handle:
                            writer=csv.DictWriter(handle,fieldnames=list(records[0]));writer.writeheader();writer.writerows(records)
            print(event["id"], "测点", len(point_data), "明细", len([r for r in all_detail if r['事件']==event['id']]),
                  "角度带", len(bands), flush=True)
    finally:
        connection.close()

    out_dir.mkdir(parents=True, exist_ok=True)
    def write(name, rows):
        if not rows: return
        with (out_dir/name).open("w", encoding="utf-8-sig", newline="") as handle:
            writer=csv.DictWriter(handle,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows)
    write("伪在线_窗口检测明细.csv", all_detail)
    write("伪在线_角度带告警.csv", all_bands)
    write("伪在线_事件验证.csv", validations)

    false_rows=[]
    for event in EVENTS:
        rows=[row for row in all_bands if row["事件"]==event["id"]]
        total_days=len({row["日期"] for row in all_detail if row["事件"]==event["id"]})
        for point in sorted({row["测点"] for row in rows}):
            selected=[row for row in rows if row["测点"]==point]
            false_rows.append({"事件":event["id"],"测点":point,"检测有效日数":total_days,
                               "观察或以上角度带窗口":len(selected),"候选或以上窗口":sum(row["等级"] in ('候选','持续预警') for row in selected),
                               "持续预警窗口":sum(row["等级"]=='持续预警' for row in selected)})
    write("伪在线_误报统计.csv",false_rows)
    return all_detail, all_bands, validations


def run_health_segments(out_dir: Path, units=("7", "8", "9")):
    """Replay unlabelled controls; triggers are not confirmed false positives."""
    detail_rows=[]
    connection=None
    cache_dir=out_dir/"unlabelled_control_features_v1"
    cache_dir.mkdir(parents=True,exist_ok=True)
    keys=[(angle,metric) for angle in range(0,360,10) for metric in METRICS]
    try:
        for unit in units:
            unit_detail_start=len(detail_rows)
            for segment_start,segment_end in HEALTH_SEGMENTS.get(unit,[]):
                point_data={}
                for cylinder in CYLINDERS:
                    for side in SIDES:
                        if (cylinder,side) in BAD.get(unit,set()) or (unit=="8" and (cylinder,side)==("六","盖侧")):continue
                        cache_path=cache_dir/f"{unit}_{cylinder}_{side}_{segment_start}_{segment_end}.npz"
                        baseline=[];daily=[];slots=[];cycle_count=0
                        if cache_path.exists():
                            with np.load(cache_path,allow_pickle=False) as cached:
                                if list(cached["metrics"])!=list(METRICS):
                                    raise ValueError(f"Feature schema mismatch: {cache_path}")
                                for kind,records in (("baseline",baseline),("daily",daily),("slots",slots)):
                                    for day,slot,values in zip(cached[f"{kind}_days"],cached[f"{kind}_slots"],cached[f"{kind}_values"]):
                                        records.append({"day":str(day),"slot":int(slot),"values":{key:float(value) for key,value in zip(keys,values) if np.isfinite(value)}})
                                cycle_count=int(cached["cycle_count"].item())
                        else:
                            if connection is None:connection=connect()
                            files=fetch_file_rows(connection,point_name(unit,cylinder,side),segment_start,segment_end)
                            by_day=defaultdict(list)
                            for row in files:
                                time=row["sample_time"] if isinstance(row["sample_time"],datetime) else datetime.fromisoformat(str(row["sample_time"]))
                                by_day[time.date().isoformat()].append({**row,"sample_time":time})
                            valid_days=0
                            for day,day_files in sorted(by_day.items()):
                                records=[]
                                day_files.sort(key=lambda row:(row["sample_time"],int(row["id"])))
                                for offset in range(0,len(day_files),100):
                                    batch=day_files[offset:offset+100]
                                    cycles=fetch_cycles(connection,[int(row["id"]) for row in batch],100)
                                    for row in batch:
                                        values=stats_for_cycle(cycles.get(int(row["id"]),np.array([])))
                                        if values is None:continue
                                        values={key:value for key,value in values.items() if key in keys and np.isfinite(value)}
                                        if not values:continue
                                        time=row["sample_time"]
                                        records.append({"sample_time":time,"day":day,"slot":time.hour//6,"values":values})
                                    del cycles
                                if not records:continue
                                valid_days+=1
                                if valid_days<=3:
                                    baseline.extend(sample_per_day(records,3))
                                    continue
                                cycle_count+=len(records)
                                daily_values=aggregate(records,"day")
                                daily.append({"day":day,"slot":-1,"values":{key:value for (group,key),value in daily_values.items()}})
                                slot_values=aggregate(records,"day_slot")
                                for slot in sorted({record["slot"] for record in records}):
                                    slots.append({"day":day,"slot":slot,"values":{key:value for (group,key),value in slot_values.items() if group==f"{day}#{slot}"}})
                            payload={"metrics":np.asarray(METRICS),"cycle_count":np.asarray(cycle_count)}
                            for kind,records in (("baseline",baseline),("daily",daily),("slots",slots)):
                                payload[f"{kind}_days"]=np.asarray([record["day"] for record in records],dtype="U10")
                                payload[f"{kind}_slots"]=np.asarray([record["slot"] for record in records],dtype=int)
                                payload[f"{kind}_values"]=np.asarray([[record["values"].get(key,np.nan) for key in keys] for record in records],dtype=float).reshape(-1,len(keys))
                            temporary=cache_path.with_suffix(".tmp.npz")
                            np.savez_compressed(temporary,**payload)
                            temporary.replace(cache_path)
                        if baseline and daily:
                            point_data[f"{cylinder}缸{side}"]={"baseline":build_baseline(baseline),"daily":daily,"slots":slots,"cycle_count":cycle_count}
                if not point_data:continue
                z_rows=[];slot_z=defaultdict(dict)
                for point,data in point_data.items():
                    for record in data["daily"]:
                        for key,value in record["values"].items():
                            base=data["baseline"][key]
                            if not np.isfinite(base["center"]):continue
                            z=(value-base["center"])/base["scale"]
                            z_rows.append({"point":point,"day":record["day"],"angle":key[0],"metric":key[1],"z":z,"baseline":base})
                    for record in data["slots"]:
                        for key,value in record["values"].items():
                            base=data["baseline"][key]
                            if not np.isfinite(base["center"]):continue
                            slot_z[(record["day"],record["slot"],key)][point]=(value-base["center"])/base["scale"]
                public=defaultdict(list)
                for row in z_rows:public[(row["day"],row["angle"],row["metric"])].append(row["z"])
                public={key:float(np.median(values)) for key,values in public.items()}
                slot_public={key:float(np.median(list(values.values()))) for key,values in slot_z.items()}
                grouped=defaultdict(list)
                for row in z_rows:
                    row["local"]=row["z"]-public[(row["day"],row["angle"],row["metric"])]
                    grouped[(row["point"],row["day"],row["angle"])].append(row)
                history=defaultdict(deque)
                for (point,day,angle),rows in sorted(grouped.items(),key=lambda item:(item[0][1],item[0][0],item[0][2])):
                    supported=[row for row in rows if abs(row["local"])>=2]
                    families={FAMILY[row["metric"]] for row in supported}
                    amps=[abs(row["local"]) for row in supported]
                    avg=float(np.mean(amps)) if amps else 0.0
                    rms=float(np.sqrt(np.mean(np.square(amps)))) if amps else 0.0
                    abnormal_slots=0
                    for slot in range(4):
                        slot_vals=[]
                        for row in rows:
                            key=(day,slot,(angle,row["metric"]))
                            if point in slot_z.get(key,{}):slot_vals.append(abs(slot_z[key][point]-slot_public[key]))
                        if sum(v>=2 for v in slot_vals)>=2:abnormal_slots+=1
                    candidate=len(families)>=2 and avg>=2 and rms>=3 and abnormal_slots>=2
                    hist=history[(point,angle)];hist.append(candidate)
                    if len(hist)>5:hist.popleft()
                    persistent=len(hist)>=3 and sum(hist)>=3
                    level="持续预警" if persistent and candidate else "候选" if candidate else "观察" if supported else "正常"
                    base_days=min(row["baseline"]["days"] for row in rows)
                    detail_rows.append({
                        "机组":unit,"运行段":f"{segment_start}~{segment_end}","日期":day,"测点":point,
                        "样本标签":"unlabelled control",
                        "备注":"2026-04-18距2026-04-01维修17天，不保证健康" if unit=="8" else "未确认健康，触发不能直接认定为误报",
                        "角度箱":angle,"等级":level,"支持特征数":len(supported),"支持特征族数":len(families),
                        "幅度平均Z":round(avg,4),"幅度中位Z":round(float(np.median(amps)) if amps else 0,4),
                        "综合RMS":round(rms,4),"异常6小时子窗数":abnormal_slots,"最近5窗候选数":sum(hist),
                        "基线有效日数":base_days,"基线置信":"低" if base_days<7 else "中" if base_days<20 else "高",
                        "特征摘要":";".join(f"{METRIC_CN[r['metric']]}{'↑' if r['local']>0 else '↓'}{abs(r['local']):.1f}Z" for r in sorted(supported,key=lambda x:-abs(x['local']))),
                    })
                print(f"unlabelled control {unit}号 {segment_start}~{segment_end}: 测点{len(point_data)}，全量检测周期{sum(x['cycle_count'] for x in point_data.values()):,}",flush=True)
            unit_rows=detail_rows[unit_detail_start:]
            if unit_rows:
                out_dir.mkdir(parents=True,exist_ok=True)
                with (out_dir/f"伪在线_未标注对照触发明细_{unit}号.csv").open("w",encoding="utf-8-sig",newline="") as handle:
                    writer=csv.DictWriter(handle,fieldnames=list(unit_rows[0]));writer.writeheader();writer.writerows(unit_rows)
    finally:
        if connection is not None:connection.close()
    out_dir.mkdir(parents=True,exist_ok=True)
    with (out_dir/"伪在线_未标注对照触发明细.csv").open("w",encoding="utf-8-sig",newline="") as handle:
        writer=csv.DictWriter(handle,fieldnames=list(detail_rows[0]) if detail_rows else ["机组","样本标签"]);writer.writeheader();writer.writerows(detail_rows)
    summary=[]
    by_unit=defaultdict(list)
    for row in detail_rows:by_unit[row["机组"]].append(row)
    for unit,rows in by_unit.items():
        days=len({row["日期"] for row in rows});points=len({row["测点"] for row in rows})
        summary.append({"机组":unit,"样本标签":"unlabelled control","对照检测有效日数":days,"测点数":points,"窗口点箱记录":len(rows),
                        "观察记录":sum(row["等级"]=="观察" for row in rows),
                        "候选记录":sum(row["等级"]=="候选" for row in rows),
                        "持续预警记录":sum(row["等级"]=="持续预警" for row in rows),
                        "候选及以上日测点箱数":len({(row["日期"],row["测点"],row["角度箱"]) for row in rows if row["等级"] in ("候选","持续预警")}),
                        "持续预警日测点箱数":len({(row["日期"],row["测点"],row["角度箱"]) for row in rows if row["等级"]=="持续预警"})})
    with (out_dir/"伪在线_未标注对照触发统计.csv").open("w",encoding="utf-8-sig",newline="") as handle:
        writer=csv.DictWriter(handle,fieldnames=list(summary[0]) if summary else ["机组","样本标签"]);writer.writeheader();writer.writerows(summary)
    return detail_rows,summary


def report(out_dir, detail, bands, validations):
    lines=["# 伪在线故障检测验证", "", "## 方法", "",
           "- 历史基线每天粗采最多3周期；24小时检测窗口处理全部有效周期。",
           "- 周期平台归一化、测点自身历史稳健Z标准化、单机公共变化扣除。",
           "- 检测不使用维修后数据；维修日期仅用于计算提前量。", ""]
    for event in EVENTS:
        rows=[row for row in validations if row["事件"]==event["id"]]
        rows.sort(key=lambda row:(-int(row["持续预警窗口数"]),-float(row["最大综合RMS"])))
        lines += [f"## {event['id']}", "", f"- {event['note'] or '按事件前数据回放'}。", f"- 形成候选的测点角度带：{len(rows)}。"]
        for i,row in enumerate(rows[:20],1):
            lines.append(f"{i}. {row['测点']} {row['角度带']}：首次候选{row['首次候选日期']}，提前{row['提前日数']}天，持续预警{row['持续预警窗口数']}窗，最大RMS{float(row['最大综合RMS']):.1f}。")
        lines.append("")
    lines += ["## 使用限制", "", "- 当前按有数据的自然日形成0~24小时窗口；下一版可改为每小时滚动24小时。",
              "- 6号基线只有3个日期，所有结果应降级解释。", "- 7号2025事件附近缺测，只能验证运行段状态变化，不能证明事件提前量。"]
    (out_dir/"伪在线检测结论.md").write_text("\n".join(lines)+"\n",encoding="utf-8")


if __name__ == "__main__":
    parser=argparse.ArgumentParser(description="局部特征伪在线回放")
    parser.add_argument("--out-dir",default="cache/indicator/伪在线检测")
    parser.add_argument("--event",action="append",choices=[event["id"] for event in EVENTS])
    parser.add_argument("--healthy",action="store_true",help="回放无已知故障的健康运行段")
    parser.add_argument("--unit",action="append",choices=UNITS,help="限制机组，可重复传入")
    args=parser.parse_args()
    output=Path(args.out_dir)
    if args.healthy:
        run_health_segments(output,tuple(args.unit or ("7","8","9")))
    else:
        detail,bands,validations=run(output,set(args.event or []))
        report(output,detail,bands,validations)
    print("结果 ->",output)