guanlan-src/app_ontology/app_ontology_guanlan/sop/stat_clean.py

# -*- coding: utf-8 -*-
"""stat_clean.py — 统计清洗层 (opt-in; CL-1 契约闸之上的第二道, 2026-07-11 并回).

来源: gym CLEAN 赛道 C1-C4 四轮训练的池首组合 full3 (hampel + frozen + scale2) + 备胎 scale(v1)
的生产化移植 — 预注册 EV: EXP-C1-01/C2-01/C3-01/C4-01; 冻结套件 clean_frozen_v0/v1;
三轮独立 lite 审 (审计留痕 outputs/_genebank/审计_2026-07-11_lite.md + EXP 追加区)。
vendored 自 gym/clean_algos/statpack.py (src 不依赖 gym; 数值语义逐位对齐, 见 DP 测试)。

纪律 (硬编码非文档):
- **mark-not-delete**: 只置 NaN, 不删行, 不改好值; 带内持续漂移 (W091 +2.1℃ 型故障前兆) 结构上
  不吃 — 属诊断域 (温度线/NBM), 清洗层保留 (EXP-C2/C3 preserve gate ≥0.98 实证, 幅度无关)。
- **效力域守卫** (每臂的"在哪灵/在哪反"来自 gym 机器记录):
  · hampel: 带内单点尖刺; σ floor 防量化死区。
  · frozen: exact-repeat 卡滞 (含 exact-0 零膨胀); 量化粗通道会假冻结 (CNV 照妖镜已验 ydk=0)。
  · scale(v1, 自身尾随): 块长≲尾随窗/2; 永久型漏 (0.046); **共模免疫** (参照自身)。
  · scale2 (fleet 中位): 长块/永久全捕 (0.992/0.996); 需兄弟机 ≥min_fleet_n;
    **>半数同向共模 = 反转 (坏台全逃+好台全杀, EXP-C4-01)** → 由 kill 熔断兜底。
- **kill 熔断 (C4 教训生产化)**: 任一通道 × 滚动 30d 标记率 > fuse_cap (默认 0.10) →
  该通道本次**拒清洗** (零置 NaN) + 响亮进 report — 宁漏不屠 (生产无 truth mask,
  熔断是反转域的唯一防线)。
- **opt-in**: 不接入 apply_contract_gate 默认路径; 上游应已过 CL-1 (契约 sentinel/valid_range)。
- **channels 选择责任 (EXP-C5-01 实证)**: frozen/hampel 语义只适用**连续物理传感通道** (温度/风速类);
  功率/状态码/设定值类有物理合法恒值与阶跃 (停机恒功率被 frozen 吃 = mima 实测 conc 0.29) — 勿入
  channels。恒值合法性本质是契约属性, 通道级语义域收紧在 C5-02 排期; 当前由调用方按此纪律选通道。
- **混合病理日 (EXP-C6-01 逮 → C6-02 已修)**: scale2 日判据 中位→**占比** (q=0.30) + 孤立混合日
  min_days=1 + 极端格连段 ≥1h 直捕 (带日合格 ≥24 守卫) — 混族块残余 100%→带隙级 (gym 实测)。
  **残余 known-limit (带隙)**: 3-5 格的 5× 连段落在 hampel (1-2 格) 与极端连段 (≥6 格) 之间 —
  E2E 实测可致下游 z 偏 ~0.6σ (EXP-C6-02 s1); 多病并发通道建议 scale2+frozen+hampel 全开并读 report。

接口:
    from src.sop.stat_clean import stat_clean
    df_clean, report = stat_clean(df, channels=[...], time_col="时间", turbine_col="风机名称")
report: per 通道 {flags_by_arm, coverage(scale2), fuse_tripped, n_marked}; 顶层 totals + params。
"""
from __future__ import annotations

import numpy as np
import pandas as pd

DEFAULTS = dict(
    hampel_window=13, hampel_k=6.0, hampel_sigma_floor=0.5,
    frozen_min_run=12,
    scale_log_thresh=0.7, scale_min_days=2, scale_trailing_days=30,
    scale_min_trailing_level=5.0, scale2_min_fleet_n=5,
    scale2_day_frac_q=0.30, scale2_min_cells_day=24, scale2_min_days=1,
    use_hampel=True, use_frozen=True, use_scale2=True, use_scale_v1=False,
    fuse_cap=0.10, fuse_window_days=30, fuse_min_days=7,
)


def _hampel_flags(s: pd.Series, window: int, k: float, floor: float) -> np.ndarray:
    med = s.rolling(window, center=True, min_periods=5).median()
    mad = (s - med).abs().rolling(window, center=True, min_periods=5).median()
    sigma = (1.4826 * mad).clip(lower=floor)
    return ((s - med).abs() > k * sigma).to_numpy() & s.notna().to_numpy()


def _frozen_flags(s: pd.Series, min_run: int) -> np.ndarray:
    grp = (s != s.shift()).cumsum()
    run_len = s.groupby(grp).transform("size")
    return (run_len >= min_run).to_numpy() & s.notna().to_numpy()


def _cell_ratio_gate(s: pd.Series, ref: np.ndarray, log_thresh: float) -> np.ndarray:
    """日级标记内的格级比值闸 (EXP-C5-03 修2, 与 gym statpack 镜像): 边界日正常格释放。"""
    with np.errstate(divide="ignore", invalid="ignore"):
        r = np.abs(np.log10(np.where((ref > 0) & (s.to_numpy() > 0), s.to_numpy() / ref, np.nan)))
    return np.nan_to_num(r, nan=0.0) >= log_thresh


def _trail_two_pass(dm: pd.Series, trailing_days: int, min_trailing_level: float,
                    log_thresh: float) -> pd.Series:
    """两遍尾随 (EXP-C5-03 修2b, gym 镜像): pass1 粗筛异常日, pass2 参照排除之。"""
    mp = max(5, trailing_days // 3)
    trail1 = dm.rolling(trailing_days, min_periods=mp).median().shift(1)
    r1 = np.abs(np.log10((dm / trail1).where((trail1 >= min_trailing_level) & (dm > 0))))
    dm_clean = dm.where(~(r1 >= log_thresh))
    return dm_clean.rolling(trailing_days, min_periods=mp).median().shift(1)


def _scale_v1_flags(s: pd.Series, days: pd.Series, log_thresh: float, min_days: int,
                    trailing_days: int, min_trailing_level: float) -> np.ndarray:
    dm = s.groupby(days.values).median()
    trail = _trail_two_pass(dm, trailing_days, min_trailing_level, log_thresh)
    ratio = dm / trail
    bad_day = (np.abs(np.log10(ratio.where((trail >= min_trailing_level) & (dm > 0)))) >= log_thresh)
    bad_day = bad_day.fillna(False)
    grp = (bad_day != bad_day.shift()).cumsum()
    run = bad_day.groupby(grp).transform("sum")
    bad_day = bad_day & (run >= min_days)
    day_flag = days.map(bad_day).fillna(False).to_numpy() & s.notna().to_numpy()
    cell_ref = days.map(trail).to_numpy(dtype="float64")
    return day_flag & _cell_ratio_gate(s, cell_ref, log_thresh)


def _scale2_flags(d: pd.DataFrame, ch: str, time_col: str, turbine_col: str,
                  log_thresh: float, min_days: int, min_level: float,
                  min_fleet_n: int, day_frac_q: float = 0.30,
                  min_cells_day: int = 24) -> tuple[np.ndarray, float]:
    s = pd.to_numeric(d[ch], errors="coerce")
    g = s.groupby(d[time_col].values)
    fm = g.transform("median")
    fn = g.transform("count")
    ok = (fm >= min_level) & (fn >= min_fleet_n) & (s > 0)
    logr = pd.Series(np.where(ok, np.log10(s / fm), np.nan), index=d.index)
    day = d[time_col].dt.floor("D")
    key = pd.MultiIndex.from_arrays([d[turbine_col].values, day.values])
    # EXP-C6-02 根修 (gym 镜像): 日入围 中位→占比, 去"当日同质"假设 (混合病理日稀释洞)
    grp_keys = [d[turbine_col].values, day.values]
    bad_cell = (logr.abs() >= log_thresh)
    elig = logr.notna()
    frac = bad_cell.groupby(grp_keys).sum() / elig.groupby(grp_keys).sum().clip(lower=1)
    n_elig = elig.groupby(grp_keys).sum()
    bad_day = ((frac >= day_frac_q) & (n_elig >= min_cells_day)).fillna(False)
    bd = bad_day.reset_index()
    bd.columns = ["tid", "day", "bad"]
    bd = bd.sort_values(["tid", "day"])
    grp = (bd["bad"] != bd.groupby("tid")["bad"].shift()).cumsum()
    run = bd.groupby(["tid", grp])["bad"].transform("sum")
    bd["bad2"] = bd["bad"] & (run >= min_days)
    lookup = bd.set_index(["tid", "day"])["bad2"]
    flags = lookup.reindex(key).fillna(False).to_numpy()
    finite = s.notna().to_numpy()
    coverage = float(ok.sum() / max(int(finite.sum()), 1))
    flags = flags & _cell_ratio_gate(s, fm.to_numpy(dtype="float64"), log_thresh)
    # 修E (EXP-C6-02, gym 镜像): 极端格连段 ≥6 (1h) 直接标 — 捕块边缘碎片; 守卫=当日合格格 ≥min_cells_day
    ext = bad_cell.fillna(False)
    run_id = pd.Series(ext.to_numpy(), index=d.index).groupby(
        d[turbine_col].values).transform(lambda x: (x != x.shift()).cumsum())
    run_len = ext.groupby([d[turbine_col].values, run_id.to_numpy()]).transform("size")
    day_ok = n_elig.reindex(key).fillna(0).to_numpy() >= min_cells_day
    flags |= (ext & (run_len >= 6)).to_numpy() & day_ok
    return flags & finite, coverage


def _fuse_check(flags: np.ndarray, finite: np.ndarray, day: np.ndarray,
                cap: float, window_days: int, min_days: int) -> tuple[bool, float, bool]:
    """滚动窗标记率 max > cap → 熔断 (C4: 集中屠杀防线; 全域均值会稀释, 必须用集中度)。
    返回 (tripped, conc_max, evaluable)。终审 F-F1 修: 数据 < min_days 天时滚动窗不可评估 —
    "无法评估 ≠ 无风险": 此时只要有任何标记即视同熔断 (宁漏不屠), evaluable=False 响亮出 report。"""
    s_k = pd.Series(flags.astype("int64")).groupby(day).sum()
    s_n = pd.Series(finite.astype("int64")).groupby(day).sum()
    roll = (s_k.rolling(window_days, min_periods=min_days).sum()
            / s_n.rolling(window_days, min_periods=min_days).sum().clip(lower=1))
    if not roll.notna().any():
        return bool(flags.any()), float("nan"), False
    conc = float(roll.max())
    return conc > cap, conc, True


_POS = "_stat_clean_pos__"
# 修1b (EXP-C5-03): +转速类 — 双稳态 regime-switching 通道, hampel/frozen 语义不适用 (镜像 gym statpack)
_NON_CONTINUOUS_RX = r"功率|power|状态|state|mode|setpoint|设定|指令|码|counter|电量|发电量|转速|rpm|rotor.?speed|gen(erator)?.?speed"


def stat_clean(df: pd.DataFrame, channels: list[str], time_col: str, turbine_col: str,
               params: dict | None = None) -> tuple[pd.DataFrame, dict]:
    """统计清洗 (mark-not-delete)。返回 (df_clean, report); **行数/行序/index 逐位不变**
    (终审 F-F2/F-F3 修: 位置法还原, 对重复/乱序 index 均成立)。"""
    import re as _re
    p = {**DEFAULTS, **(params or {})}
    # F-F4 临时闸 (C5-02 契约语义域落地前): 非连续物理通道 (恒值/阶跃合法) 须显式豁免
    if not p.get("allow_non_continuous"):
        bad = [c for c in channels if _re.search(_NON_CONTINUOUS_RX, str(c), _re.I)]
        if bad:
            raise ValueError(
                f"channels 含非连续物理通道 {bad}: 功率/状态/设定类恒值合法, frozen/hampel 语义不适用 "
                "(EXP-C5-01 缺陷①: 停机恒功率被吃)。确需处理传 params={'allow_non_continuous': True}")
    assert _POS not in df.columns, f"保留列名冲突: {_POS}"
    d = df.copy()
    if not np.issubdtype(d[time_col].dtype, np.datetime64):
        d[time_col] = pd.to_datetime(d[time_col])
    d[_POS] = np.arange(len(d))
    d = d.sort_values([turbine_col, time_col], kind="mergesort")
    inv = np.argsort(d[_POS].to_numpy(), kind="stable")   # 位置逆排 → 输入行序精确还原
    day_arr = d[time_col].dt.floor("D").values
    days_ser = d[time_col].dt.floor("D")
    report = {"params": {k: v for k, v in p.items()}, "channels": {}, "discipline":
              "mark-not-delete; 熔断=宁漏不屠 (>半数共模反转防线, EXP-C4-01); opt-in 于 CL-1 之上"}
    for ch in channels:
        if ch not in d.columns:
            report["channels"][ch] = {"error": "column missing"}
            continue
        col = pd.to_numeric(d[ch], errors="coerce")
        finite = col.notna().to_numpy()
        flags = np.zeros(len(d), dtype=bool)
        by_arm, coverage = {}, None
        for _, idx in d.groupby(turbine_col, sort=False).indices.items():
            s = col.iloc[idx]
            if p["use_hampel"]:
                f = _hampel_flags(s, p["hampel_window"], p["hampel_k"], p["hampel_sigma_floor"])
                flags[idx] |= f
                by_arm["hampel"] = by_arm.get("hampel", 0) + int(f.sum())
            if p["use_frozen"]:
                f = _frozen_flags(s, p["frozen_min_run"])
                flags[idx] |= f
                by_arm["frozen"] = by_arm.get("frozen", 0) + int(f.sum())
            if p["use_scale_v1"]:
                f = _scale_v1_flags(s, days_ser.iloc[idx], p["scale_log_thresh"],
                                    p["scale_min_days"], p["scale_trailing_days"],
                                    p["scale_min_trailing_level"])
                flags[idx] |= f
                by_arm["scale_v1"] = by_arm.get("scale_v1", 0) + int(f.sum())
        if p["use_scale2"]:
            dd = d.copy()
            dd[ch] = col
            f, coverage = _scale2_flags(dd, ch, time_col, turbine_col, p["scale_log_thresh"],
                                        p["scale2_min_days"], p["scale_min_trailing_level"],
                                        p["scale2_min_fleet_n"],
                                        day_frac_q=p["scale2_day_frac_q"],
                                        min_cells_day=p["scale2_min_cells_day"])
            flags |= f
            by_arm["scale2"] = int(f.sum())
        n_coerced = int((col.isna() & pd.notna(d[ch])).sum())   # to_numeric 强转损失, 响亮记账 (F-F5⑤)
        tripped, conc, evaluable = _fuse_check(flags, finite, day_arr, p["fuse_cap"],
                                               p["fuse_window_days"], p["fuse_min_days"])
        if tripped:
            report["channels"][ch] = {
                "fuse_tripped": True, "fuse_evaluable": evaluable,
                "conc_max": (round(conc, 4) if evaluable else None),
                "flags_by_arm": by_arm, "n_marked": 0, "coverage": coverage,
                "n_coerced_kept_raw": n_coerced,
                "note": ("熔断: 滚动窗标记率超 fuse_cap → 本通道拒清洗 (疑共模反转/参数失配), 人工复核"
                         if evaluable else
                         "熔断(不可评估): 数据 < fuse_min_days 天, 集中度无法评估且存在标记 → 拒清洗 (宁漏不屠, F-F1)")}
            continue
        d[ch] = col.where(~flags)
        report["channels"][ch] = {"fuse_tripped": False, "fuse_evaluable": evaluable,
                                  "conc_max": round(conc, 4),
                                  "flags_by_arm": by_arm, "n_marked": int(flags.sum()),
                                  "coverage": coverage, "n_coerced_to_nan": n_coerced}
    report["total_marked"] = int(sum(c.get("n_marked", 0) for c in report["channels"].values()))
    d = d.iloc[inv].drop(columns=[_POS])   # F-F2/F-F3: 位置法还原 — 行数/行序/index 逐位=输入
    return d, report