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

# -*- coding: utf-8 -*-
"""thermreg_validation.py — 温控执行器卡死族分型: 反环温签名判别器 (加热器卡ON vs 泛过热).

裁决对象 = per 机组 **过热机制分型** — 上游 (temp_fleet_nbm / temp-nbm skill) 报"这台过热"但分不清:
  · 加热器卡ON = **反环温** (舱外温↓ 时残差↑, 残差 ∝ −T_amb) — memory thermal-regulation-stuck-family 模态①
  · 冷却卡死 / 泛过热 = 环温无关持续过热 (残差 ≈ 常偏置)
两者残差都升 → 需二阶判别量。**判别量 = per台 [fleet-NBM残差 vs 舱外温] 日均回归斜率**。

物理 trace (第一性, no black box): 温控加热器卡满档 → 恒温器追 T_set 注热 ∝ (T_set − T_amb)
  → 舱外越冷注热越多 → 残差 ∂/∂T_amb < 0 (反环温)。泛过热 (冷却退化/常偏置) → ∂/∂T_amb ≈ 0。
  NBM 已含 fleet-共模环温项 (ambient ∈ OLS) → pooled 残差 ⊥ 环温 (整场), 但 **per台** 残差仍可与
  环温相关 ⟺ 该台偏离 fleet-共模环温响应 (EXP-M3-16 §8.2 教训)。

机制 (双门 AND, 同 sensor_drift_verdict 主-fp-cap+显著性模式):
  1. fleet-池化 OLS 残差 (内联 gym temp_fleet_fullz.pooled_resid 语义, 见下) → 日均聚合 (降 10min 自相关)
  2. per台 OLS 斜率 β_i = slope(resid_day ~ ambient_day)
  3. **置换显著性门**: 逐台 shuffle 环温标签 × n_perm → null β 分布 → 双侧 p (+1 平滑) < p_gate
  4. **健康 fp cap**: 斜率跨台稳健 z (median + 1.4826·MAD) → |z| ≥ z_cap
  flag = 两门 AND; **两侧**: β<0 = anti_ambient(加热器卡ON, 主目标) / β>0 = pro_ambient(异常正环温敏感)。
  泛过热 (β 量级正常于 fleet) 被 cap 正确拒 → 不 flag = 本判别器价值。

★关键发现 (EXP-M3-22 §6, manipulation check 实证, DP test 锁回归): **双门皆 load-bearing, 缺一不可**。
  常偏置后残留的 per台环温斜率仍可与 0 区分 → 泛过热台置换 p 亦显著 → **置换单门会误报泛过热 (2/2)
  + 全健康台 (8/8)**。分离主力 = fp cap ("斜率异常于 fleet 否", 特异); 置换门 = "斜率非零否" (灵敏不特异)。
  `cap_off=True` 仅为 manipulation-check/debug 保留 (gym 版 env GYM_ANTIAMB_CAPOFF), 生产禁用。

★残差为何内联而非复用 sensor_validation._fleet_residual (比对后决定, 语义不一致):
  _fleet_residual = OLS[1, P, rpm, T_amb] 无风速项、无月哑元、缺列零占位、fit≥50 行即走;
  gym pooled_resid = OLS[1, P, rpm, T_amb, wind, 月哑元 m2..m12] + dropna + 池化行数 ≥ MIN_POOL_ROWS(5000)。
  月哑元 + 风速 = fleet-共模季节/风况吸收项, 对"per台环温斜率 = 偏离共模响应"这一判别量 load-bearing
  (少吸收 → 健康台残差留季节结构 → 斜率底噪抬 → cap 定标漂)。忠实并回 (阈值/口径不重调) → 内联
  `_pooled_fleet_resid` 逐语句镜像 gym `temp_fleet_fullz.pooled_resid`, 不复用 _fleet_residual。

★对 gym 版 (gym/algos/thermreg_antiamb.py, EXP-M3-22 PROMOTED) 的产化改动 (其余逐语句一致):
  (a) **种子策略**: gym = 单 rng(seed·7919+2203) 顺序流 (p 值依赖台遍历顺序/机群构成) → 产化改
      **sha256 逐台确定性** (seed = sha256(f"antiamb|{tid}")[:8], 同 yaw_validation estgate 既有产线模式)
      → 同输入恒同输出、增删台互不扰。**后果 (诚实, RV-1 独立审订正措辞)**: slope/z_slope/direction
      恒逐字段相等 (斜率链路无 rng); 但 **perm_sig/flag 在 p≈p_gate 边界台可因种子策略差异翻转**
      (MC 噪声, 审核员对抗构造实测门口 |Δp|≈0.02-0.04、翻转 7/180 — 非语义变化, 两种 rng 都是合法
      MC 抽样; 高信噪台 p 打 1/201 地板不受影响)。置换 p 计算体与 gym 逐语句一致 (同 rng 状态 bitwise 等)。
  (b) 单窗口 (无 gym `_phase` 场景机器): 调用方切好一致性窗再喂; 多相位需求 = 调用方分窗多次调用。
  (c) 弃权显式化 (gym 静默 skip): 无环温列 / 池化行数不足 / 有效日数<min_days / fleet 斜率台数<min_turb
      → 全部带 note 的 abstain (INSUFFICIENT), 非静默全零"干净"。
  (d) cap_off 由 env 改显式关键字参数。

★纪律 / 诚实边界 (EXP-M3-22 §6 ledger, 未调参掩盖):
  - **候选级筛查**, 非确诊: flag = "反环温签名成立" → 现场核温控执行器/加热器才定谳。
  - **无变频器 gold**: 验证 = 合成 probe + INJ-only + feicheng 净场 fp 上界 (0/25); 无真实卡ON 标注召回。
  - **fleet-相对假设少数台异常**: majority-fleet 卡ON 会污染 median → fp cap 失效 (场景外声明;
    memory relative-criterion-batch-blindness 同族)。批次嫌疑场先做场型识别再用本判别器。
  - **全窗静态斜率 → 无 onset 维度** (何时卡的答不了); 需提前量走时序化 (backlog)。
  - **合成 probe 信噪比高于真注入** (z≈−44 vs §8.2 真注入 z≈−4.9): 真盘幅度更温和, 方向/分离一致。
  - 置换为日级 shuffle (10min 行级自相关已由日均聚合大幅降, 非严格块 bootstrap)。

来源: gym G4 家族③ thermreg_antiamb (EXP-M3-22, PROMOTED 2026-07-12) → §7 并回 2026-07-15。
并回闸门: gate-1 gym suite 复现 (8 passed) + gate-2 忠实度对拍/DP test (tests/sop/test_thermreg_antiamb_dp.py)
本会话完成; gate-3 RV-1 独立审由主会话另行 spawn (本文件不自证已审)。
纯核心无 I/O 可测; per-farm 加载在调用方 (skill/report 层)。
"""
from __future__ import annotations

import hashlib

import numpy as np
import pandas as pd

N_PERM = 200          # 逐台置换重数 (null β 分布)
P_GATE = 0.05         # 置换双侧 p 门
Z_CAP = 3.0           # 斜率跨台稳健 |z| 门 = 健康 fp cap (3.0 控净场极值台假阳)
MIN_DAYS = 30         # per台最少有效日数 (斜率回归 n)
MIN_TURB = 5          # fleet-相对 z 最少台数
MIN_POOL_ROWS = 5000  # 池化 OLS 最少行数 (镜像 temp_fleet_fullz.MIN_POOL_ROWS)


def _pooled_fleet_resid(sub: pd.DataFrame, channel: str, load_cols: list,
                        wind_col: str | None, turbine_col: str, time_col: str,
                        min_pool_rows: int):
    """池化 OLS (载荷 + 风速 + 月哑元) → 行级残差列 _r。None = 样本不足。

    逐语句镜像 gym temp_fleet_fullz.pooled_resid (列序: 1, load_cols..., wind, m2..m12);
    唯一放宽: wind_col 可 None (无风速列时降级掉该项, 调用方须知残差模型变弱)。"""
    need = [turbine_col, time_col, channel] + ([wind_col] if wind_col else []) + load_cols
    d = sub[need].dropna()
    if len(d) < min_pool_rows:
        return None
    y = d[channel].to_numpy(float)
    cols = [np.ones(len(d))] + [d[c].to_numpy(float) for c in load_cols]
    if wind_col:
        cols.append(d[wind_col].to_numpy(float))
    mth = d[time_col].dt.month
    for m in range(2, 13):
        cols.append((mth == m).to_numpy(float))
    X = np.column_stack(cols)
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    return d.assign(_r=y - X @ beta)


def _slope(x: np.ndarray, y: np.ndarray) -> float:
    """OLS 斜率 slope(y ~ x) = cov(x,y)/var(x)。x 已保证非常量 (调用前 check)。"""
    xc = x - x.mean()
    vx = float(xc @ xc)
    if vx <= 0:
        return float("nan")
    return float((xc @ (y - y.mean())) / vx)


def _perm_p(x: np.ndarray, y: np.ndarray, beta_obs: float, rng, n_perm: int) -> float:
    """逐台置换零: shuffle 环温标签 → null β 分布 → 双侧 p = P(|β_null| ≥ |β_obs|), +1 平滑。

    计算体逐语句镜像 gym thermreg_antiamb._perm_p (同 rng 状态下 bitwise 等); rng 由调用方
    以 sha256 逐台种子构造 (产化确定性, 见模块 docstring 改动 a)。"""
    if not np.isfinite(beta_obs):
        return 1.0
    xc = x - x.mean()
    vx = float(xc @ xc)
    if vx <= 0:
        return 1.0
    yc = y - y.mean()
    ge = 1                                   # +1 平滑 (含观测本身, 免 p=0)
    for _ in range(n_perm):
        yp = rng.permutation(yc)
        b = float(xc @ yp) / vx
        if abs(b) >= abs(beta_obs):
            ge += 1
    return ge / (n_perm + 1)


def thermreg_antiamb_verdict(
    df: pd.DataFrame, *, channel: str, ambient_col: str, time_col: str, turbine_col: str,
    power_col: str, rpm_col: str | None = None, gen_flag_col: str | None = None,
    wind_col: str | None = None, gen_power_min: float = 50.0,
    n_perm: int = N_PERM, p_gate: float = P_GATE, z_cap: float = Z_CAP,
    min_days: int = MIN_DAYS, min_turb: int = MIN_TURB, min_pool_rows: int = MIN_POOL_ROWS,
    cap_off: bool = False,
) -> dict:
    """per-机组 反环温签名裁决 (温控卡死族分型: 加热器卡ON vs 泛过热)。

    返回 {tid: {"flag": bool, "score": |z_slope|, "slope": 残差-环温斜率(℃/℃),
                "z_slope": 斜率跨台稳健z, "p_perm": 置换双侧p, "perm_sig": bool, "cap_hit": bool,
                "direction": "anti_ambient"|"pro_ambient"|None, "n_days": 有效日数, "note"?: 弃权原因}}。
    flag = (p_perm < p_gate) ∧ (|z_slope| ≥ z_cap) [cap_off=True 关后门, 仅 manipulation-check]。
    direction: slope<0 → anti_ambient (加热器卡ON 主目标) / slope>0 → pro_ambient / 弃权 → None。
    df 要求: time_col 为 datetime64; 调用方保证清洗后单一一致性窗 (阈值/口径见模块 docstring, 不重调)。
    """
    # NaN-safe: Arrow str backend 下 astype(str) 留 float NaN → sorted 混型崩; dropna + 显式 str()
    tids = sorted({str(x) for x in df[turbine_col].dropna().unique()})
    NULL = {"flag": False, "score": 0.0, "slope": 0.0, "z_slope": 0.0, "p_perm": None,
            "perm_sig": False, "cap_hit": False, "direction": None, "n_days": 0}
    if ambient_col is None or ambient_col not in df.columns:
        return {t: {**NULL, "note": "无 ambient_col → 反环温判别不适用; abstain (INSUFFICIENT)"}
                for t in tids}
    if channel not in df.columns or power_col not in df.columns:
        return {t: {**NULL, "note": "channel/power 列缺; abstain"} for t in tids}
    if not pd.api.types.is_datetime64_any_dtype(df[time_col]):
        return {t: {**NULL, "note": "time_col 非 datetime64; abstain (调用方先转)"} for t in tids}
    # gen 掩码 (gym 语义: flag≥0.5, 否则 P>gen_power_min)
    if gen_flag_col and gen_flag_col in df.columns:
        gen = df[gen_flag_col].astype(float) >= 0.5
    else:
        gen = df[power_col] > float(gen_power_min)
    sub = df[gen]
    out = {t: dict(NULL) for t in tids}
    load_cols = ([power_col] + ([rpm_col] if rpm_col and rpm_col in df.columns else [])
                 + [ambient_col])
    rd = _pooled_fleet_resid(sub, channel, load_cols,
                             wind_col if (wind_col and wind_col in df.columns) else None,
                             turbine_col, time_col, min_pool_rows)
    if rd is None:
        for t in tids:
            out[t]["note"] = f"池化行数<{min_pool_rows}; abstain (INSUFFICIENT)"
        return out
    # 日均聚合 (降 10min 自相关) → per台 斜率 + 逐台置换 p
    daily = (rd.assign(_day=rd[time_col].dt.floor("D"))
               .groupby([turbine_col, "_day"])
               .agg(_rd=("_r", "mean"), _ad=(ambient_col, "mean"))
               .reset_index())
    slopes: dict = {}
    pvals: dict = {}
    for tid, g in daily.groupby(turbine_col):
        key = str(tid)
        x = g["_ad"].to_numpy(float)
        y = g["_rd"].to_numpy(float)
        m = np.isfinite(x) & np.isfinite(y)
        x, y = x[m], y[m]
        out[key]["n_days"] = int(len(x))
        if len(x) < min_days:
            out[key]["note"] = f"有效日数 {len(x)}<{min_days}; abstain"
            continue
        if x.std() <= 0:
            out[key]["note"] = "环温日均零方差; abstain"
            continue
        b = _slope(x, y)
        seed = int(hashlib.sha256(f"antiamb|{key}".encode()).hexdigest()[:8], 16)
        slopes[key] = b
        pvals[key] = _perm_p(x, y, b, np.random.default_rng(seed), n_perm)
    for t in tids:
        # RV-1 must-fix: 整台被 dropna/gen 掩码淘汰的机组从不进 daily 循环 → 原静默无 note
        # (违 docstring (c) "全部带 note"; 独立审逮)。补显式弃权。
        if out[t]["n_days"] == 0 and "note" not in out[t]:
            out[t]["note"] = f"有效日数 0<{min_days} (通道/环温全NaN或非gen态, 整台无有效日); abstain"
    keys = [k for k in slopes]
    if len(keys) < min_turb:
        for t in tids:
            if t in slopes:                    # 斜率/p 仍报 (信息), 但无跨台 z → 不可裁
                out[t].update({"slope": round(slopes[t], 4), "p_perm": round(pvals[t], 4),
                               "perm_sig": bool(pvals[t] < p_gate)})
            out[t]["note"] = (f"fleet 有效斜率台数 {len(keys)}<{min_turb} → 跨台稳健 z 不可估; "
                              "abstain (INSUFFICIENT)")
        return out
    # 健康 fp cap 轴: 斜率跨台稳健 z (median + 1.4826·MAD; MAD=0 退化 std)
    arr = np.array([slopes[k] for k in keys], float)
    med = float(np.median(arr))
    mad = float(np.median(np.abs(arr - med))) * 1.4826
    sd = mad if mad > 0 else float(arr.std(ddof=1))
    for k in keys:
        z = (slopes[k] - med) / sd if sd and sd > 0 else 0.0
        perm_sig = pvals[k] < p_gate
        cap_hit = abs(z) >= z_cap
        out[k].update({
            "flag": bool(perm_sig and (cap_hit or cap_off)),
            "score": round(abs(float(z)), 4),
            "slope": round(slopes[k], 4),
            "z_slope": round(float(z), 2),
            "p_perm": round(pvals[k], 4),
            "perm_sig": bool(perm_sig),
            "cap_hit": bool(cap_hit),
            "direction": "anti_ambient" if slopes[k] < 0 else "pro_ambient",
        })
    return out