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

# -*- coding: utf-8 -*-
"""yaw_validation.py — 偏航静态对风误差: 强估计器 + 逐台置换显著性门 (检出轴).

裁决对象 = per 机组 **静态偏航对风误差** 是否显著 (筛查轴) — 出 flag=候选, 非确诊;
场级弥散/共模判读留调用方 (fleet_common_mode as-is 叠加), 绝对零位须现场独立量 (yaw-static-error skill)。

来源: gym Y 偏航赛道 yaw_estgate (EXP-Y2-01 PI 首轮换届 → 走设计§7四闸并回 src/ 2026-07-12)。
gym 分数≠交付; 本件 = 走完 §7 四闸 (①冻结suite+holdout均显著优现役臂 ②3原则审 ③RV-1 lite独立审
④labeled DP test) 后的生产实现; DP 验收 = tests/sop/test_yaw_estgate_dp.py。

★核心主张 (为何换届, EXP-Y2-01 ledger 复现):
  裸估计器 (domefit/cosp) 检出强但当筛查器 **泄漏** — perm 置换零假设下 FP 1.0/0.833 (健康机群乱报);
  panel 守卫安全 (perm-FP 0) 但检出弱 (INJ 0.185)。本臂 = est_cosp 估计 (as-is 复用铁律) + 穹顶深度
  逐台置换显著性门, 填 "安全∧检出" 的空缺: perm-FP 0.056 (近 panel) ∧ INJ 检出 0.611 (3.3× panel)。
  数字 (ledger_yaw EXP-Y2-01, V-019 RV-2 复现): perm 0.0555 / recall_d9(non-holdout) 0.722 /
  recall_d13(non-holdout) 0.833 / holdout(yuxian) d9 0.889·d13 0.667 (泛化成立) / tys-INV fp 0.061。
  ★口径 (RV-1 逮, 防误读): non-holdout recall 含 pingyin 近全场 flag 抬升 (pyn INV fp 0.80 = 场级共模
  非本臂误报) → 0.833/0.722 非"干净检出力"; 真检出力参考 = tys 干净场 (d13 recall 0.889 ∧ fp 0.033)。

★物理/统计 trace (第一性, no black box):
- 物理: 真偏航失配 → 功率随对风角呈 cos³-型穹顶 (系统性结构); 健康台无穹顶 (浅/无结构)。
- 统计: 固定对风角分箱, shuffle 归一功率 r 得穹顶深度 null 分布 → p=(1+#{D_perm≥D_obs})/(B+1);
  perm 评分场景 (功率已被打乱) → D_obs 本身=null 抽样 → p≈0.5 不 flag (by construction 堵泄漏)。
  flag = |θ̂|≥flag_deg ∧ p<alpha (幅度门 + 显著性门 双条件)。
- 复现: 置换种子 = sha256(f'estgate|{tid}') (非 PYTHONHASHSEED 随机化的 hash()) → 逐台确定性。

★纪律 (gym 教训 + RV-1):
- 本臂 = **检出轴筛查器**, 出 flag=候选级, 非确诊、非绝对零位。绝对对风零位须现场独立量。
- 共模不进本臂 (scalar=检出力轴); pyn/yux 高 perm-FP = 场级风标共模真弥散 (fleet_common_mode 判读,
  归 INSUFFICIENT 场级候选), 非本臂误报 — 调用方须叠共模判读层。
- 采样/数据面: 需 native 对风角通道 (对风角度) + 归一功率 r; 无 native 对风角 → 不适用。

纯核心 (yaw_misalignment_verdict) 无 I/O 可测; per-farm canonical_clean 在调用方 (skill/report 层),
与 sensor_validation 同分工。文献锚 (素材非机制来源, EV 裁决唯一权威): OpenOA StaticYawMisalignment
(Monte-Carlo 重采样 CI) / 置换检验 = 标准 change-point 显著性范式。
"""
from __future__ import annotations

import hashlib

import numpy as np
import pandas as pd

# CFG 默认 (源 scripts/yaw_reliability_panel.CFG 的相关子集; est_cosp/穹顶几何锚定)
_DOME_WIN = 20      # 对风角窗 ±20°
_DOME_BIN = 1.0     # 1° 分箱
_BIN_MIN_N = 30     # 有效仓最小样本
_EDGE_LO, _EDGE_HI = 12.0, 20.0   # 边缘基线取 |center|∈[12,20]


def _binned_mean(wd: np.ndarray, r: np.ndarray, dome_win: int, dome_bin: float, bin_min_n: int):
    """1° 分箱加权均值剖面 (est_cosp 用)。返回 centers, cnt, mean, valid。"""
    edges = np.arange(-dome_win - 0.5, dome_win + 0.5 + 1e-9, dome_bin)
    centers = (edges[:-1] + edges[1:]) / 2
    idx = np.digitize(wd, edges) - 1
    ok = (idx >= 0) & (idx < len(centers))
    cnt = np.bincount(idx[ok], minlength=len(centers))
    ssum = np.bincount(idx[ok], weights=r[ok], minlength=len(centers))
    valid = cnt >= bin_min_n
    mean = np.where(valid, ssum / np.maximum(cnt, 1), np.nan)
    return centers, cnt, mean, valid


def _est_cosp(wd: np.ndarray, r: np.ndarray, dome_win: int, dome_bin: float, bin_min_n: int) -> float:
    """模型法估计对风峰: 网格搜 θ 使 r≈a·cos²(wd−θ)+b 加权最优 (源 est_cosp, as-is 复用)。"""
    centers, cnt, mean, valid = _binned_mean(wd, r, dome_win, dome_bin, bin_min_n)
    if valid.sum() < 8:
        return np.nan
    x = centers[valid]; y = mean[valid]; w = cnt[valid].astype(float)
    best = (np.inf, np.nan)
    for th in np.arange(-dome_win, dome_win + 0.01, 0.5):
        c2 = np.cos(np.radians(x - th)) ** 2
        A = np.column_stack([c2, np.ones_like(c2)])
        W = np.sqrt(w)
        beta, *_ = np.linalg.lstsq(A * W[:, None], y * W, rcond=None)
        if beta[0] <= 0:
            continue
        sse = float(np.sum(w * (y - A @ beta) ** 2))
        if sse < best[0]:
            best = (sse, th)
    return float(best[1])


def _dome_depth(idx_ok: np.ndarray, r_ok: np.ndarray, nbin: int,
                valid_edge: np.ndarray, bin_min_n: int) -> float:
    """向量化穹顶深度 = 平滑 mean 峰 − 边缘基线 mean (纯 numpy, 供置换重复调用)。"""
    cnt = np.bincount(idx_ok, minlength=nbin)
    ssum = np.bincount(idx_ok, weights=r_ok, minlength=nbin)
    with np.errstate(invalid="ignore", divide="ignore"):
        mean = np.where(cnt >= bin_min_n, ssum / np.maximum(cnt, 1), np.nan)
    sm = np.full(nbin, np.nan)
    for i in range(nbin):
        w = mean[max(0, i - 2):i + 3]
        if np.sum(~np.isnan(w)) >= 3:
            sm[i] = np.nanmean(w)
    if np.all(np.isnan(sm)):
        return np.nan
    peak = np.nanmax(sm)
    base_vals = mean[valid_edge & ~np.isnan(mean)]
    base = np.nanmedian(base_vals) if base_vals.size else np.nanmin(sm)
    return float(peak - base)


def yaw_misalignment_verdict(
    per_turbine: dict, *, vane_col: str = "对风角度", r_col: str = "r",
    flag_deg: float = 3.0, alpha: float = 0.05, n_perm: int = 200,
    dome_win: int = _DOME_WIN, dome_bin: float = _DOME_BIN, bin_min_n: int = _BIN_MIN_N,
    edge_lo: float = _EDGE_LO, edge_hi: float = _EDGE_HI, min_n: int = 1500,
) -> dict:
    """per-机组 静态偏航对风误差裁决 (检出轴筛查器)。

    参数
    ----
    per_turbine : {tid: DataFrame} — 每台已 canonical_clean 的数据, 含 vane_col (对风角, 已居中 ±180)
        + r_col (归一功率 = 功率 / per-风速箱中位)。调用方负责清洗 (与 sensor_validation 同分工)。

    返回 {tid: {"flag": bool, "est_deg": 对风峰角(°)|None, "p_perm": 置换p|None,
                "dome_depth": 穹顶深度|None, "tier": None|"INSUFFICIENT"|"EST_FAIL"|"DEPTH_NAN"}}。
    flag = |est_deg|≥flag_deg ∧ p_perm<alpha (幅度门 + 逐台置换显著性门)。tier≠None → 不参评 (数据/估计不足)。
    """
    out: dict = {}
    for tid, d in per_turbine.items():
        if d is None or len(d) < min_n:
            out[tid] = {"flag": False, "est_deg": None, "p_perm": None,
                        "dome_depth": None, "tier": "INSUFFICIENT"}
            continue
        wd = np.asarray(d[vane_col], dtype=float)
        r = np.asarray(d[r_col], dtype=float)
        est = _est_cosp(wd, r, dome_win, dome_bin, bin_min_n)
        if est is None or np.isnan(est):
            out[tid] = {"flag": False, "est_deg": None, "p_perm": None,
                        "dome_depth": None, "tier": "EST_FAIL"}
            continue
        edges = np.arange(-dome_win - 0.5, dome_win + 0.5 + 1e-9, dome_bin)
        centers = (edges[:-1] + edges[1:]) / 2
        nbin = len(centers)
        idx = np.digitize(wd, edges) - 1
        ok = (idx >= 0) & (idx < nbin)
        idx_ok = idx[ok]
        r_ok = r[ok]
        valid_edge = (np.abs(centers) >= edge_lo) & (np.abs(centers) <= edge_hi)
        d_obs = _dome_depth(idx_ok, r_ok, nbin, valid_edge, bin_min_n)
        if np.isnan(d_obs):
            out[tid] = {"flag": False, "est_deg": round(float(est), 2), "p_perm": None,
                        "dome_depth": None, "tier": "DEPTH_NAN"}
            continue
        # 逐台置换零假设 (可复现种子 = sha256(tid), 非 PYTHONHASHSEED 随机化的 hash())
        seed = int(hashlib.sha256(f"estgate|{tid}".encode()).hexdigest()[:8], 16)
        rng = np.random.default_rng(seed)
        ge = 0
        for _ in range(n_perm):
            d_perm = _dome_depth(idx_ok, rng.permutation(r_ok), nbin, valid_edge, bin_min_n)
            if not np.isnan(d_perm) and d_perm >= d_obs:
                ge += 1
        p = (1 + ge) / (n_perm + 1)
        out[tid] = {"flag": bool(abs(est) >= flag_deg and p < alpha),
                    "est_deg": round(float(est), 2), "p_perm": round(p, 4),
                    "dome_depth": round(d_obs, 4), "tier": None}
    return out


# ==================== 动态偏航校核 (yaw-dynamic, 与上静态正交) ====================

def _wrap180(v: np.ndarray) -> np.ndarray:
    """折到 ±180°。"""
    return (np.asarray(v, float) + 180.0) % 360.0 - 180.0


def _active_blocks(active: np.ndarray):
    """命令位 active(0/1) 的极大连续块 [(start,end)] 闭区间 (end 含); 相邻块间必有 0 分隔。
    抗采样: 不依赖 nacelle 精度, 命令位一置位即计机动 (逐秒差分坑规避)。"""
    a = np.asarray(active).astype(int)
    if a.size == 0:
        return []
    d = np.diff(np.concatenate([[0], a, [0]]))
    starts = np.flatnonzero(d == 1)
    ends = np.flatnonzero(d == -1) - 1
    return list(zip(starts.tolist(), ends.tolist()))


def _robust_cv_outliers(vals, cv_floor: float = 0.05, k_mad: float = 3.5):
    """fleet 高离群守卫: hi_thresh = median + k·**地板 scale** (scale=max(1.4826·MAD, cv_floor·|median|))。
    ★RV-1 逮 (Finding2): 原 MAD≈0/紧 bulk → robust_cv<cv_floor 短路成"uniform 无离群" → **漏掉 clean 2× 离群**
      (你最想抓的台)。修: scale 加地板 (cv_floor·|median|, 防 MAD=0 归零) → 紧 bulk 仍抓明显离群;
      **批次盲仍保** (全场一致高 → 无值超 hi_thresh → 无 flag)。uniform 字段降为信息位 (紧 bulk 标记), 不再抑制 flag。"""
    v = np.asarray([x for x in vals if x is not None and np.isfinite(x)], float)
    if v.size < 5:
        return {"uniform": None, "robust_cv": None, "hi_thresh": None}
    med = float(np.median(v)); mad = float(np.median(np.abs(v - med)))
    rcv = (1.4826 * mad) / abs(med) if med != 0 else np.inf
    scale = (max(1.4826 * mad, cv_floor * abs(med)) if med != 0
             else (1.4826 * mad if mad > 0 else float(np.std(v))))   # 地板 scale: MAD=0 不归零
    return {"uniform": bool(rcv < cv_floor), "robust_cv": round(rcv, 4), "median": round(med, 3),
            "hi_thresh": round(med + k_mad * scale, 3)}


def yaw_dynamic_verdict(
    df: pd.DataFrame, *, turbine_col: str, time_col: str,
    cw_col: str | None = None, ccw_col: str | None = None, nacelle_col: str | None = None,
    yaw_to_wind_col: str | None = None, wind_dir_col: str | None = None,
    power_col: str | None = None, gen_flag_col: str | None = None, gen_power_min: float = 50.0,
    coarse_dt_min: float = 0.25, rate_mag_min: float = 5.0, min_gen_n: int = 100,
    uniform_cv_floor: float = 0.05, k_mad: float = 3.5, static_bias_route_deg: float = 3.0,
) -> dict:
    """per-机组 **动态偏航控制行为** 校核 (时变: 偏航次数/活动量/动态σ/死区/速率) — 与
    yaw_misalignment_verdict (静态 vane 零点偏, 标定轴) **正交** (本轴=控制行为学, 时序驱动)。

    返回 {
      "verdict": "CANDIDATE_RANKING" | "INSUFFICIENT",
      "sampling": {"dt_min": 中位采样间隔, "tier": "full(≤~1min)" | "coarse(死区/速率降级)"},
      "per_turbine": {tid: {maneuvers_per_day, active_duty, yaw_travel_deg_per_day, dyn_sigma_deg,
          static_bias_deg, deadband_engage_deg(+_tier), yaw_rate_deg_per_min(+_tier),
          travel_per_wdvar(环境调整参考), n_maneuvers, n, flag, flag_axis, note?}},
      "fleet": {metric: {median,p10,p90,n}}, "guards": {metric: uniform-guard}, "honesty": [...], "note": ...
    }

    ★机制 (yaw-dynamic skill 骨架, no black box):
      - 偏航次数 = 命令位 active=(cw>0)|(ccw>0) 极大连续块数 (命令位缺 → nacelle |Δ|>阈派生, 弱, note)。
      - 活动量 yaw_travel = Σ 块净位移(|Δnacelle|)/天 (actuation 代理); active_duty = active 占比。
      - 动态 σ = std(对风角|gen) = vane 口径**上限/相对** (非 LiDAR 真对风); 静态偏置 = mean(对风角|gen)。
      - 死区[降级] = 块起点前 1 样本 |对风角| 中位; 速率[降级] = 块净位移/时长 (不逐样本差分, latch 坑)。
      - flag = fleet 稳健 CV 守卫下的**高偏航次数/高σ 离群** (候选级筛查, 非确诊)。
    ★confound (headline 诚实界, 主动前置 — 同 fatigue TI-exposure 环境混杂): 偏航次数/travel/σ **被
      风向变率环境混杂** — veering/复杂地形/场缘位置的台**自然偏航更多/σ更大**, 非偏航系统故障。故高偏航
      flag = **actuation 暴露/跟踪候选, 非磨损故障确诊**; 分离故障须归一风向变率暴露或现场核。
      **travel_per_wdvar (travel÷风向std) 参考字段 (非 flag 轴) — RV-1 逮 (F1): coarse dt(60s) 下 std(Δwind_dir)
      抓高频 jitter 非低频 veer/drift (死区滤掉 jitter), 真盘 corr(maneuvers,tpwv)=0.797 = 不可靠分离 → 仅参考,
      别当分离依据。** 环境无关轴 = 静态偏置(标定, 显著则路由静态 skill) + 死区(控制 setpoint, coarse 采样降级)。
    ★采样率闸 (数据面, 承 §47 第五角): dt>coarse_dt_min → nacelle latch → **死区/速率不可测 (标 _tier=DOWNGRADE,
      报但不下发)**; 只偏航次数(命令位抗采样)+σ 粗筛; ≤~1min 才全档。10min 场 → 死区/速率 INSUFFICIENT。
    ★诚实界:
      - **候选/筛查级非确诊**: σ相对(真值须 LiDAR); 活动量=actuation代理非真磨损/寿命 (真载荷=INSUFFICIENT);
        下发校正/寿命值必现场独立量。
      - **需命令位(cw/ccw) 或 nacelle 位序 + 对风角**: 命令位缺→nacelle派生(弱); 对风角缺→无σ/bias; 皆缺→INSUFFICIENT。
      - **fleet-relative→批次盲**: 全场一致高偏航(整场veer环境)→uniform→无一flag (承 relative-criterion 批次盲)。
      - **静态偏置路由**: |static_bias|>static_bias_route_deg → 疑 vane 标定偏 → 路由 yaw_misalignment_verdict
        (静态轴); 本臂 static_bias≈0 = 动态主导无标定偏, 不路由。

    ★产化精化 (对 gym champion) + 数据面坑 (承 §47 第五角):
      - **σ/bias gen-mask**: 本件 dyn_sigma/static_bias 用 **gen 态** (p_active>gen_power_min 或 gen_flag) →
        hailesi σ 12.6° (排停机野对风角) vs gym 全行 24.9° = 更净的**跟踪态**口径; maneuvers/day 命令位口径
        与 gym **逐位一致** (42.35, gen 无关)。
      - **turbine-id 列须 populated**: hailesi `turbine_hash` 列**全 NaN** (坏列) → 须用 `turbine_idx`;
        调用方须核 turbine_col 非全空 (第五角数据面: 空/坏 id 列静默塌成单组)。
      - **coarse_dt_min=0.25 (15s)**: 死区/速率须 sub-15s 分辨 (yaw 机动 deadband 跨越 ~5-15s); >15s → DOWNGRADE。

    来源: gym yawdyn 赛道 (EXP-YD1-01 bootstrap TRUE, 5闸PASS: H1物理3.93×/H2 corr0.77/H3 recovery1.0/H4反臆造0;
    唯 hailesi 60s 命令位数据面 GO) → §7 四闸并回 src/ 2026-07-15。DP = tests/sop/test_yaw_dynamic_dp.py。
    """
    HONESTY = [
        "(a) confound headline: 偏航次数/travel/σ 被风向变率环境混杂 (veer/复杂地形位置的台自然偏航多) → "
        "高偏航=actuation暴露/跟踪候选非磨损故障确诊; 分离须归一风向变率(travel_per_wdvar参考)或现场核。",
        "(b) 环境无关轴 = 静态偏置(标定, 显著则路由静态skill) + 死区(控制setpoint, coarse采样降级); 其余=环境暴露轴。",
        "(c) 候选/筛查级非确诊: σ=vane相对上限(真对风须LiDAR); 活动量=actuation代理非真磨损/寿命(真载荷INSUFFICIENT)。",
        "(d) 采样率闸: dt>coarse→nacelle latch→死区/速率不可测(_tier DOWNGRADE不下发); 只次数+σ粗筛; ≤~1min全档。",
        "(e) 需命令位或nacelle位序+对风角; 命令位缺→nacelle派生(弱); 对风角缺→无σ/bias; 皆缺→INSUFFICIENT。",
        "(f) fleet-relative→批次盲: 全场一致高偏航(整场veer)→uniform→无一flag; 批次嫌疑场须绝对基准/场型识别。",
    ]
    # NaN-safe: Arrow str backend 下 astype(str) 留 float NaN → sorted 混型崩; dropna + 显式 str()
    tids = sorted({str(x) for x in df[turbine_col].dropna().unique()})

    def _insuf(reason):
        return {"verdict": "INSUFFICIENT", "sampling": None,
                "per_turbine": {t: {"flag": False, "note": reason} for t in tids},
                "fleet": None, "guards": None, "honesty": HONESTY, "note": reason}

    has_cmd = cw_col in df.columns and ccw_col in df.columns if (cw_col and ccw_col) else False
    if not has_cmd and not (nacelle_col and nacelle_col in df.columns):
        return _insuf("无命令位(cw/ccw) 且无 nacelle 位序 → 偏航机动不可辨; INSUFFICIENT")
    # bad-ID 守卫 (RV-1 F3 逮, 第五角数据面): turbine_col 全空/无有效值 → astype(str) 塌成单组"nan" 静默 → 弃权
    if df[turbine_col].notna().sum() == 0 or (len(tids) == 1 and tids[0] in ("nan", "None", "")):
        return _insuf("turbine id 列全空/坏列 (如 hailesi turbine_hash 全 NaN) → 无法分组, 塌单组; INSUFFICIENT (换有效 id 列如 turbine_idx)")

    tv = df[turbine_col].astype(str).to_numpy()
    t_all = pd.to_datetime(df[time_col], errors="coerce")
    # 采样率闸 (全局中位)
    dt_min = float(t_all.diff().dt.total_seconds().median()) / 60.0 if len(df) > 1 else np.nan
    coarse = (not np.isfinite(dt_min)) or dt_min > coarse_dt_min
    dtier = "coarse(死区/速率降级)" if coarse else "full(≤~1min)"

    per = {}
    for t in tids:
        m = tv == t
        d = df[m].copy()
        ti = pd.to_datetime(d[time_col], errors="coerce")
        order = np.argsort(ti.to_numpy())
        d = d.iloc[order]; ti = ti.iloc[order]
        n = len(d)
        dtm = float(ti.diff().dt.total_seconds().median()) / 60.0 if n > 1 else np.nan
        days = n * dtm / 1440.0 if (np.isfinite(dtm) and dtm > 0) else np.nan
        # active 命令位 (优先) 或 nacelle 派生
        if has_cmd:
            cw = pd.to_numeric(d[cw_col], errors="coerce").fillna(0).to_numpy()
            ccw = pd.to_numeric(d[ccw_col], errors="coerce").fillna(0).to_numpy()
            active = ((cw > 0) | (ccw > 0)).astype(int); man_src = "命令位"
        else:
            nacd = pd.to_numeric(d[nacelle_col], errors="coerce").to_numpy()
            dn = np.abs(_wrap180(np.diff(nacd, prepend=nacd[:1])))
            active = (dn > rate_mag_min).astype(int); man_src = "nacelle派生(弱)"
        nac = pd.to_numeric(d[nacelle_col], errors="coerce").to_numpy() if (nacelle_col and nacelle_col in d.columns) else None
        ywa = _wrap180(pd.to_numeric(d[yaw_to_wind_col], errors="coerce").to_numpy()) if (yaw_to_wind_col and yaw_to_wind_col in d.columns) else None
        if gen_flag_col and gen_flag_col in d.columns:
            gen = pd.to_numeric(d[gen_flag_col], errors="coerce").to_numpy() > 0
        elif power_col and power_col in d.columns:
            gen = pd.to_numeric(d[power_col], errors="coerce").to_numpy() > gen_power_min
        else:
            gen = np.ones(n, bool)

        blocks = _active_blocks(active)
        n_man = len(blocks)
        travel = 0.0; rates = []; pre = []
        for s, e in blocks:
            if nac is not None and np.isfinite(nac[e]) and np.isfinite(nac[s]):
                disp = abs(_wrap180(np.array([nac[e] - nac[s]]))[0]); travel += disp
                dur = (e - s + 1) * (dtm if np.isfinite(dtm) else 1.0)
                if disp >= rate_mag_min and dur > 0:
                    rates.append(disp / dur)
            if ywa is not None and s >= 1 and np.isfinite(ywa[s - 1]):
                pre.append(abs(ywa[s - 1]))
        yg = ywa[gen & np.isfinite(ywa)] if ywa is not None else np.array([])
        has_align = yg.size >= min_gen_n
        wd = pd.to_numeric(d[wind_dir_col], errors="coerce").to_numpy() if (wind_dir_col and wind_dir_col in d.columns) else None
        wd_var = float(np.nanstd(_wrap180(np.diff(wd)))) if (wd is not None and n > 2) else None
        rec = {
            "n_maneuvers": int(n_man),
            "maneuvers_per_day": round(n_man / days, 2) if (days and days > 0) else None,
            "active_duty": round(float(active.mean()), 4) if n else None,
            "yaw_travel_deg_per_day": round(travel / days, 2) if (days and days > 0 and nac is not None) else None,
            "dyn_sigma_deg": round(float(np.std(yg)), 2) if has_align else None,
            "static_bias_deg": round(float(np.mean(yg)), 2) if has_align else None,
            "deadband_engage_deg": round(float(np.median(pre)), 2) if pre else None,
            "deadband_tier": "DOWNGRADE" if coarse else None,
            "yaw_rate_deg_per_min": round(float(np.median(rates)), 2) if rates else None,
            "yaw_rate_tier": "DOWNGRADE" if coarse else None,
            "travel_per_wdvar": round((travel / days) / wd_var, 3) if (wd_var and wd_var > 0 and days and days > 0 and nac is not None) else None,
            "n": int(n), "man_src": man_src, "flag": False, "flag_axis": None,
        }
        notes = []
        if man_src.startswith("nacelle"):
            notes.append("命令位缺→nacelle派生机动(弱): 次数受 latch/量化影响, 命令位场更准")
        if has_align and abs(rec["static_bias_deg"]) > static_bias_route_deg:
            notes.append(f"静态偏置 {rec['static_bias_deg']}°>{static_bias_route_deg}°: 疑vane标定偏→路由 yaw_misalignment_verdict(静态轴)")
        if notes:
            rec["note"] = "; ".join(notes)
        per[t] = rec

    # fleet 稳健 CV 守卫 + 高离群 flag (次数/σ 双轴, 候选级)
    def _col(k):
        return [per[t].get(k) for t in tids]
    fleet, guards = {}, {}
    for k in ("maneuvers_per_day", "active_duty", "dyn_sigma_deg", "yaw_travel_deg_per_day", "static_bias_deg"):
        v = np.asarray([x for x in _col(k) if x is not None and np.isfinite(x)], float)
        fleet[k] = ({"median": round(float(np.median(v)), 3), "p10": round(float(np.quantile(v, .1)), 3),
                     "p90": round(float(np.quantile(v, .9)), 3), "n": int(v.size)} if v.size else None)
    for k in ("maneuvers_per_day", "dyn_sigma_deg"):                # 高离群 = 候选 (环境暴露 caveat 见 honesty a)
        g = _robust_cv_outliers(_col(k), uniform_cv_floor, k_mad)
        guards[k] = g
        if g.get("hi_thresh") is not None:                          # RV-1 F2 修: hi_thresh 判离群 (不看 uniform, 紧 bulk 也抓)
            for t in tids:
                val = per[t].get(k)
                if val is not None and np.isfinite(val) and val > g["hi_thresh"]:
                    per[t]["flag"] = True
                    per[t]["flag_axis"] = (per[t]["flag_axis"] + "+" + k) if per[t]["flag_axis"] else k

    n_flag = sum(1 for t in tids if per[t]["flag"])
    note = (f"动态偏航控制行为 (候选级筛查, 非确诊): 采样 {round(dt_min,2) if np.isfinite(dt_min) else '?'}min "
            f"[{dtier}]; {len(tids)} 台; 高偏航暴露/σ 离群 flag {n_flag} 台 (环境混杂 caveat 见 honesty a)。")
    if coarse:
        note += " ⚠ coarse 采样: 死区/速率降级 (报不下发)。"
    return {"verdict": "CANDIDATE_RANKING", "sampling": {"dt_min": round(dt_min, 3) if np.isfinite(dt_min) else None, "tier": dtier},
            "per_turbine": per, "fleet": fleet, "guards": guards, "honesty": HONESTY, "note": note}


# ============================================================================
# 电缆扭角 / 解缆健康 (cable-twist / unwind) — 低侧·绝对单侧台判据 (round-2)
# ============================================================================
# gym EXP-CT1-02 (round-2, TRUE, H6∧H7∧H8) → 本件。round-1 (EXP-CT1-01) 检出器成立但
# **通道独立 premise 被 H2 证伪**: twist_ang ≈ nacelle_pos ≈ yaw_ang (逐台 corr 1.0) = 别名同
# 累积机舱转角 (承 LESSONS #50 aliasing 型 / YD1-01 nacelle_pos 同源)。**诚实保留**: CT = 该共享
# 通道上的 **distinct 分析轴** (解缆健康: 单侧漂移/回中节律), **非 distinct 传感通道**; twist_ang
# 不作独立佐证源。round-2 不翻 round-1 此结论, 只补检出的低侧盲 (从不回中的单侧台)。

_CT_MIN_SAMPLES = 500     # 60s 采样, 10 天 ≈14400/台; < = 数据不足 (源 ct_control.MIN_SAMPLES)
_CT_P_MIN = 50.0          # kW; > = 发电 (活跃度门常量)
_CT_RS_MIN = 0.05         # rpm; > = 旋转


def _ct_metrics(v: np.ndarray, pwr, rpm, min_samples: int, p_min: float, rs_min: float):
    """per 台从 twist 时序 (+功率/转速) 一阶统计 (源 ct_control._metrics, 逐位等价)。
    drift=|mean| 单侧不回中 / unwind_returns=过零次数 回中节律 / absmax=峰值累积 /
    twist_travel=Σ|Δ| 偏航驱动行程 (停机冻结→≈0) / active_frac=发电∨旋转占比 / rotor_spin_frac。
    无功率/转速通道 → active_frac/rotor_spin_frac=None (停机门退化仅 travel; 向后兼容 round-1)。"""
    v = v[np.isfinite(v)]
    if v.size < min_samples:
        return None
    zc = int(np.sum(np.diff(np.sign(v)) != 0))
    m = {"drift": round(float(abs(np.mean(v))), 3), "absmax": round(float(np.max(np.abs(v))), 3),
         "span": round(float(np.max(v) - np.min(v)), 3), "net": round(float(abs(v[-1] - v[0])), 3),
         "unwind_returns": zc, "twist_travel": round(float(np.sum(np.abs(np.diff(v)))), 3),
         "n_min": int(v.size), "days": round(v.size / 1440.0, 2)}
    if pwr is not None or rpm is not None:
        p = pwr[np.isfinite(pwr)] if pwr is not None else np.array([])
        r = rpm[np.isfinite(rpm)] if rpm is not None else np.array([])
        gen = (p > p_min) if p.size else np.array([], bool)
        spin = (r > rs_min) if r.size else np.array([], bool)
        if p.size and r.size and p.size == r.size:
            active = float(np.mean(gen | spin))
        else:
            active = max(float(np.mean(gen)) if gen.size else 0.0,
                         float(np.mean(spin)) if spin.size else 0.0)
        m["active_frac"] = round(active, 4)
        m["rotor_spin_frac"] = round(float(np.mean(spin)) if spin.size else 0.0, 4)
    else:
        m["active_frac"] = None
        m["rotor_spin_frac"] = None
    return m


def _ct_parked(m: dict, travel_min: float, act_min: float, rs_min_frac: float) -> bool:
    """停机/冻结门 (排假阳): twist_travel < TRAVEL_MIN ∨ (active_frac < ACT_MIN ∧ rotor_spin_frac < RS_MIN)。
    源 ct_runner._is_parked 逐位等价。低 unwind 的**冻结高偏置台**必被此门排除 (本轮头号证伪风险)。"""
    tr = m.get("twist_travel")
    frozen = tr is not None and np.isfinite(tr) and tr < travel_min
    af, rf = m.get("active_frac"), m.get("rotor_spin_frac")
    offline = (af is not None and rf is not None and af < act_min and rf < rs_min_frac)
    return bool(frozen or offline)


def cable_twist_verdict(
    df: pd.DataFrame, *, twist_col: str, turbine_col: str,
    power_col: str | None = None, rpm_col: str | None = None, nacelle_col: str | None = None,
    drift_abs: float = 180.0, travel_min: float = 2000.0,
    act_min: float = 0.05, rs_min_frac: float = 0.05, unwind_floor: float = 1.0,
    near_miss_margin: float = 2.0,
    p_min: float = _CT_P_MIN, rs_min: float = _CT_RS_MIN, min_samples: int = _CT_MIN_SAMPLES,
) -> dict:
    """per-机组 **电缆扭角/解缆健康** 低侧·绝对单侧台校核 (候选级筛查, 非确诊)。

    捞『从不回中的单侧台』(低 unwind ∧ 高绝对 drift ∧ 排停机) —— 补 round-1 drift 高侧 MAD 相对守卫
    对『全场都偏但某几台从不回中』的结构盲。物理: 机舱主动偏航但扭角全程一侧不过零 = 电缆持续单侧
    应力/解缆失败风险。**低 unwind ≠ 单侧应力** (停机不偏航→twist 冻结→也 0 过零但无累积) → 停机门必备。

    返回 {tid: {"drift": |mean(twist)|°, "unwind_returns": 过零次数, "absmax_turns": 峰值圈数,
                "twist_travel": Σ|Δtwist|°, "active_frac": 发电∨旋转占比, "rotor_spin_frac": ...,
                "candidate": bool, "parked": bool|None, "note"?}}
    + {"_fleet": {"unwind_low_thr": max(fleet_p10(unwind), floor), "candidates": [...], "parked_excluded": [...],
                  "near_miss": [...] (高 drift 但 unwind 刚过刀刃门, 交付必并看; RV-1 #53 精确集断言盲修),
                  "n_candidates": ..., "verdict": CANDIDATE_RANKING|INSUFFICIENT, "alias_corr_nac": ..., "honesty": [...]}}。

    判据 (源 gym ct_runner.classify_lowside, 逐位等价):
      unwind_low = max(fleet_p10(unwind_returns), unwind_floor)         # hailesi p10=1 → unwind<1 即从不回中
      parked     = twist_travel < travel_min ∨ (active_frac<act_min ∧ rotor_spin_frac<rs_min_frac)
      candidate  = (unwind_returns < unwind_low) ∧ (drift ≥ drift_abs) ∧ ¬parked

    ★物理/统计 trace (第一性, no black box): 全是 twist(+功率/转速) 时序一阶统计, 无谱/无振动 (CMS 出圈)。
      flag 只挂 drift(单侧性)+停机门, 不挂 span/absmax → NULL 去均值(drift→0)时不臆造 flag (gym H8 反臆造锁)。

    ★别名前提 (承 LESSONS #50/#51, 数据面第0步硬门 — **本件最重要诚实界**):
      twist_ang ≈ nacelle_pos ≈ yaw_ang (hailesi 逐台 corr **1.0** = 别名同累积机舱转角, gym round-1 H2 证伪
      『distinct 传感通道』)。**CT = 该共享通道上 distinct 分析轴 (单侧/回中), 非独立传感通道**; twist_ang
      不作独立佐证源。传 nacelle_col → _fleet.alias_corr_nac 报逐台中位 corr (≥0.99 → note 别名警示)。

    ★诚实界 (兑现 gym EXP-CT1-02 §诚实界):
      1. **候选级非确诊**: drift/圈数 = 几何代理, 非电缆真应力/疲劳寿命 (需 OEM 解缆 spec + 现场核)。
      2. **触限余量硬限 = INSUFFICIENT**: 无 OEM 解缆触发圈数 → 只报 absmax 圈数, 不给『距限剩 X°』。
      3. **twist=0=neutral 为前提级假设** (数据支持: fleet 绕 0 震荡; 非 OEM 参考零点确认) → 若零点非
         cable-neutral 而是任意偏置, drift=|mean|/过零的『单侧/回中』语义偏, 前提破则候选清单须重算。
      4. **窗口界**: 『N 天 0 次回中』= 窗内观测, 更长窗可能出现迟发解缆 → 候选 window-bounded 非永久单侧断言。
      5. **停机门 n=1 天然样本** (hailesi 仅 tid12 天然全停机) + PARK 合成补冻结高偏置强混淆; 真·停机-高偏置
         天然样本尚缺 (跨场泛化补)。
      6. turbine id 列全空 → INSUFFICIENT (承 yaw turbine_hash 坏列, §48 F3)。

    n=1 场 (hailesi 远景 EN182-6250, 64台×60s×10天; 候选 tid 33/46/30/55) [暂行]。
    DP 验收 = tests/sop/test_cable_twist_dp.py。gym EXP-CT1-02 TRUE (H6∧H7∧H8), RV-1 独立审挂 §7。
    """
    honesty = [
        "候选级非确诊: drift/圈数=几何代理 (非电缆真应力/寿命); 触限硬限=INSUFFICIENT (需 OEM 解缆 spec)。",
        "别名前提 (LESSONS #50/#51): twist_ang≈nacelle_pos≈yaw_ang (corr 1.0) = 共享通道; CT=distinct 分析轴 (单侧/回中) 非独立传感通道; twist_ang 不作独立佐证源。",
        "twist=0=neutral 为前提级假设 (需 OEM 参考零点确认); 前提破则候选清单须重算。",
        "『N 天 0 次回中』= window-bounded 窗内观测, 非永久单侧断言; 更长窗可能出现迟发解缆。",
        "停机门天然样本 n=1 (仅全停机台) + PARK 合成补冻结混淆; 真停机-高偏置天然样本尚缺 (跨场泛化补)。",
        "刀刃阈敏感 · 候选集非穷尽 (RV-1 #53): unwind<max(p10,1) 严门把『穿零 1 次 vs 0 次』(物理近同的单侧台) "
        "隐去 → 全场最高 drift 台可能落 near_miss (hailesi tid52 drift291°>候选首台) 而非 candidates; "
        "交付必并看 _fleet.near_miss (预注册严门不动, near_miss=诊断轴)。",
    ]
    if turbine_col not in df.columns or df[turbine_col].notna().sum() == 0:
        return {"_fleet": {"verdict": "INSUFFICIENT", "reason": "turbine id 列缺失/全空 (承 §48 F3 坏 id 列)",
                           "honesty": honesty}}

    # per-台 metrics
    metrics: dict[str, dict] = {}
    for tid, g in df.groupby(df[turbine_col].astype(str)):
        v = pd.to_numeric(g[twist_col], errors="coerce").to_numpy(float)
        pwr = (pd.to_numeric(g[power_col], errors="coerce").to_numpy(float)
               if power_col and power_col in g else None)
        rpm = (pd.to_numeric(g[rpm_col], errors="coerce").to_numpy(float)
               if rpm_col and rpm_col in g else None)
        m = _ct_metrics(v, pwr, rpm, min_samples, p_min, rs_min)
        metrics[str(tid)] = m if m is not None else {"drift": None, "unwind_returns": None,
                                                      "n_min": int(np.isfinite(v).sum()), "insufficient": True}

    # unwind_low 阈 (bottom decile, 地板 unwind_floor)
    u = np.asarray([m.get("unwind_returns") for m in metrics.values()
                    if m.get("unwind_returns") is not None], float)
    unwind_low = max(float(np.quantile(u, 0.1)), unwind_floor) if u.size >= 5 else None

    # 别名 corr 诚实探针 (可选 nacelle_col)
    alias_corr = None
    if nacelle_col and nacelle_col in df.columns:
        cs = []
        for _tid, g in df.groupby(df[turbine_col].astype(str)):
            a = pd.to_numeric(g[twist_col], errors="coerce").to_numpy(float)
            b = pd.to_numeric(g[nacelle_col], errors="coerce").to_numpy(float)
            mk = np.isfinite(a) & np.isfinite(b)
            if mk.sum() > 200 and np.std(a[mk]) > 0 and np.std(b[mk]) > 0:
                cs.append(float(np.corrcoef(a[mk], b[mk])[0, 1]))
        alias_corr = round(float(np.median(cs)), 5) if cs else None

    per: dict[str, dict] = {}
    cands, parked_excl, near_miss = [], [], []
    for tid, m in metrics.items():
        if m.get("insufficient") or m.get("drift") is None or m.get("unwind_returns") is None:
            per[tid] = {"candidate": False, "parked": None, "insufficient": True,
                        "note": f"数据不足 (<{min_samples} 有效样本)"}
            continue
        d, un = m["drift"], m["unwind_returns"]
        is_parked = _ct_parked(m, travel_min, act_min, rs_min_frac)
        low_unwind = unwind_low is not None and un < unwind_low
        high_drift = d >= drift_abs
        is_cand = bool(low_unwind and high_drift and not is_parked)
        # near-miss: 高 drift ∧ ¬停机 但 unwind 刚过刀刃 (∈[thr, thr+margin)) — 强制 unwind<thr 门吞掉的
        # 高单侧台 (物理: 10天穿零 1 次 vs 0 次 对单侧应力几无差别)。诊断轴, **不进 candidates** (预注册严门不动)。
        is_near = bool(high_drift and not is_parked and not is_cand and unwind_low is not None
                       and un < unwind_low + near_miss_margin)
        rec = {"_tid": tid, "drift": round(d, 2), "unwind_returns": un,
               "absmax_turns": round(m.get("absmax", 0) / 360.0, 2) if m.get("absmax") is not None else None,
               "twist_travel": m.get("twist_travel"), "active_frac": m.get("active_frac"),
               "rotor_spin_frac": m.get("rotor_spin_frac")}
        if is_parked and low_unwind and high_drift:
            parked_excl.append(rec)   # 低unwind+高drift 但停机门排除 (关键反混淆记录)
        if is_cand:
            cands.append(rec)
        if is_near:
            near_miss.append(rec)     # 刀刃邻域高 drift 台 (交付摘要必带, 防精确集断言吞极端台)
        per[tid] = {"drift": round(d, 2), "unwind_returns": un,
                    "absmax_turns": rec["absmax_turns"], "twist_travel": m.get("twist_travel"),
                    "active_frac": m.get("active_frac"), "rotor_spin_frac": m.get("rotor_spin_frac"),
                    "candidate": is_cand, "parked": is_parked, "near_miss": is_near,
                    "low_unwind": bool(low_unwind), "high_drift": bool(high_drift)}

    verdict = "CANDIDATE_RANKING" if unwind_low is not None else "INSUFFICIENT"
    cands = sorted(cands, key=lambda x: -x["drift"])
    near_miss = sorted(near_miss, key=lambda x: -x["drift"])
    note = (f"电缆扭角/解缆健康 (候选级筛查, 非确诊): {len(per)} 台; 低侧单侧候选 {len(cands)} 台 "
            f"(unwind<{unwind_low} ∧ drift≥{drift_abs}° ∧ ¬停机); 停机门排除 {len(parked_excl)} 台; "
            f"刀刃 near-miss {len(near_miss)} 台 (高 drift 但 unwind 刚过门)。")
    if near_miss:
        note += (f" ⚠ 候选集非穷尽: near-miss 首台 {near_miss[0]['_tid']} "
                 f"(drift {near_miss[0]['drift']}°>候选首台, unwind {near_miss[0]['unwind_returns']}) — "
                 f"刀刃阈敏感, 边缘台须并看。")
    if alias_corr is not None and alias_corr >= 0.99:
        note += f" ⚠ twist≈nacelle 别名 (corr {alias_corr}): distinct 分析轴非独立通道, 不作独立佐证。"
    return {"_fleet": {"verdict": verdict, "unwind_low_thr": unwind_low, "n_candidates": len(cands),
                       "candidates": cands, "parked_excluded": sorted(parked_excl, key=lambda x: -x["drift"]),
                       "n_parked_excluded": len(parked_excl), "near_miss": near_miss,
                       "n_near_miss": len(near_miss), "near_miss_margin": near_miss_margin,
                       "alias_corr_nac": alias_corr,
                       "cfg": {"drift_abs": drift_abs, "travel_min": travel_min, "act_min": act_min,
                               "rs_min_frac": rs_min_frac, "unwind_floor": unwind_floor},
                       "honesty": honesty, "note": note},
            **per}


# ---------------------------------------------------------------- §4.4⑪ 归一器独立性三角
def normalizer_independence_verdict(df, *, vane_col, p_col, ws_cols, band=(3.0, 9.5),
                                    bin_w=0.25, mad_k=5.0, vane_abs_max=25.0,
                                    dome_win=20, dome_bin=1.0, bin_min_n=30,
                                    agree_deg=2.0, null_deg=2.0, min_n=2500):
    """§4.4⑪ dome 归一器独立性三角 (fushan/jtxq 2026-07-17 实证抽象, 同行 swap 检验).

    动机: 机舱风速计本身 vane 相关 (转子尾流传函) → nacelle-ws 分箱归一可在真值≈0 时制造
    fleet 级假峰 (fushan -9°, 传感器级坐实); 粒度/分半/扇区/形状四检查全在归一器家族内,
    对归一器共模零判别力 → 独立性 swap 是必做腿。
    设计: **同一批行** (所有归一器带内且非空 = 同行), 逐归一器 (0.25 ws-bin 中位归一 + 5×MAD)
    拟对风峰 (_est_cosp), 比峰位:
      - 各腿峰齐 (max 两两差 ≤ agree_deg)            → INDEPENDENT_CONSISTENT (峰位过独立性门)
      - 主腿 |峰| 大而独立腿 |峰| ≤ null_deg          → NORMALIZER_CONFOUND_SUSPECT (fushan 型, 主腿伪影)
      - 其余 (互不相干/独立腿自身可疑)               → INCONCLUSIVE_DIVERGENT (jtxq 型, 须第三锚/现场)
      - 有效腿 <2 或同行 n<min_n                      → INSUFFICIENT
    ws_cols: dict {name: col}, **第一个 = 主归一器 (通常 nacelle ws), 其余 = 独立归一器**
    (mast 自由来流 / 其他 vane 无关风代理)。调用方自行保证独立腿 QC (mast 结冰周剔除等)。
    诚实界: 峰位精度未经现场独立量验证 (铁律不变); INDEPENDENT_CONSISTENT ≠ 确诊, 只是
    归一器共模被排除的候选门; 邻机功率比作独立腿须先过方向稳定性守卫 (skill ❌条)。
    """
    import pandas as pd
    names = list(ws_cols)
    if len(names) < 2:
        return {"verdict": "INSUFFICIENT", "reason": "须 ≥2 归一器 (主 + ≥1 独立)", "peaks": {}}
    cols = [vane_col, p_col] + [ws_cols[n] for n in names]
    d = df[cols].dropna()
    m = d[vane_col].abs() <= vane_abs_max
    for n in names:
        c = d[ws_cols[n]]
        m &= (c >= band[0]) & (c <= band[1])
    d = d[m]
    if len(d) < min_n:
        return {"verdict": "INSUFFICIENT", "reason": f"同行 n={len(d)}<{min_n}", "peaks": {},
                "n_samerows": int(len(d))}
    peaks = {}
    for n in names:
        g = d.copy()
        wb = (g[ws_cols[n]] / bin_w).astype(int)
        med = g.groupby(wb)[p_col].transform("median")
        mad = (g[p_col] - med).abs().groupby(wb).transform("median")
        g = g[(g[p_col] - med).abs() <= mad_k * np.maximum(mad, 1.0)]
        wb = (g[ws_cols[n]] / bin_w).astype(int)
        r = g[p_col] / np.maximum(g.groupby(wb)[p_col].transform("median"), 1.0)
        th = _est_cosp(g[vane_col].to_numpy(float), r.to_numpy(float), dome_win, dome_bin, bin_min_n)
        peaks[n] = None if not np.isfinite(th) else round(float(th), 2)
    valid = {n: v for n, v in peaks.items() if v is not None}
    if len(valid) < 2:
        return {"verdict": "INSUFFICIENT", "reason": "有效腿<2 (估计器未出值)", "peaks": peaks,
                "n_samerows": int(len(d))}
    vals = list(valid.values())
    spread = round(float(max(vals) - min(vals)), 2)
    primary = names[0]
    indep_ok = [n for n in names[1:] if n in valid]
    if spread <= agree_deg:
        verdict = "INDEPENDENT_CONSISTENT"
    elif (primary in valid and abs(valid[primary]) > agree_deg
          and indep_ok and all(abs(valid[n]) <= null_deg for n in indep_ok)):
        verdict = "NORMALIZER_CONFOUND_SUSPECT"
    else:
        verdict = "INCONCLUSIVE_DIVERGENT"
    return {"verdict": verdict, "peaks": peaks, "spread_deg": spread,
            "primary": primary, "n_samerows": int(len(d)),
            "note": "峰位过/不过归一器独立性门; 绝对零位仍须现场独立量 (§4.4⑥ 铁律)"}