# -*- coding: utf-8 -*-
"""wind_coupling.py — 风资源→设备 耦合映射 框架 (farm-agnostic, SOP §0.9 应用范式).

把风资源从"去混基线"升级为贯穿透镜: 分解设备信号里哪些风况/载荷驱动 vs 机器内禀,
并建 风况↔故障↔可靠性 关联矩阵。三张映射 = 纯函数, 吃标准接口(cleaned df + fault df + 列名参数),
OEM/场无关。四条护栏焊进代码(均承第一性原理②confounder排):
  ① 带符号 z (非 |z|)         — map_wind_health: 内禀残差方向不能丢(平陆13#低速轴承 |z|误叙热=冷 blood lesson)
  ② 状态混杂检验(gen_frac)    — map_wind_fault: 无则停机态触发的子系统误报"低风驱动"假映射
  ③ 运行时长归一 + 排除活跃度   — map_wind_reliability: 防"暴露=活跃度"logging伪相关
  ④ 前提门在调用方(限电/投产剥离) — 本模块只吃已剥离的发电态数据(gen_mask), 不自剥

结论强度上限 = 关联/趋势级(机舱风 self_ref + 相关≠因果); 【定论】须现场取证 + 独立风况量。
校准 = scripts/wind_coupling_benchmark.py (纯合成注入-回收已知真值卡)。
设计正本 = outputs/sop_review/风资源设备映射_分析域设计草案.md。
"""
from __future__ import annotations
import numpy as np
import pandas as pd

__all__ = ['map_wind_health', 'map_wind_fault', 'map_wind_reliability', 'map_wind_load',
           'map_device_to_wind', 'rainflow_del', 'default_subsystem', 'spearman', 'ks_stat']


# ----------------------------- 数值 helpers (无 scipy) -----------------------------
def _r2(y, yhat):
    ss_res = np.nansum((y - yhat) ** 2)
    ss_tot = np.nansum((y - np.nanmean(y)) ** 2)
    return float(1 - ss_res / ss_tot) if ss_tot > 0 else np.nan


def _ols_pred(X, y):
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    return X @ beta


def spearman(x, y):
    """Spearman r + 正态近似 p + n (无 scipy). 返回 (r, p, n)."""
    from math import erf
    x = np.asarray(x, float); y = np.asarray(y, float)
    m = ~(np.isnan(x) | np.isnan(y))
    x, y = x[m], y[m]
    n = len(x)
    if n < 5:
        return (np.nan, np.nan, n)
    rx = pd.Series(x).rank().to_numpy(); ry = pd.Series(y).rank().to_numpy()
    if np.std(rx) == 0 or np.std(ry) == 0:
        return (0.0, 1.0, n)
    r = float(np.corrcoef(rx, ry)[0, 1])
    t = r * np.sqrt((n - 2) / max(1e-12, 1 - r * r))
    p = float(2 * (1 - 0.5 * (1 + erf(abs(t) / np.sqrt(2)))))
    return (r, p, n)


def ks_stat(a, b):
    """两样本 KS 统计量 (无 scipy)."""
    a = np.sort(np.asarray(a, float)); a = a[~np.isnan(a)]
    b = np.sort(np.asarray(b, float)); b = b[~np.isnan(b)]
    if len(a) < 20 or len(b) < 20:
        return np.nan
    allv = np.concatenate([a, b])
    ca = np.searchsorted(a, allv, 'right') / len(a)
    cb = np.searchsorted(b, allv, 'right') / len(b)
    return float(np.max(np.abs(ca - cb)))


def default_subsystem(desc: str) -> str:
    """状态码描述 → 子系统 (行业约定关键词; 可被调用方自定 subsystem_fn 覆盖)."""
    d = str(desc)
    if any(k in d for k in ['桨叶', '变桨', '桨距', '91°', 'pitch']): return '变桨/桨叶'
    if '主轴承' in d or 'main bearing' in d.lower(): return '主轴承'
    if any(k in d for k in ['变流器', '变频器', 'converter', 'inverter']): return '变流器'
    if any(k in d for k in ['偏航', '解缆', '扭缆', 'yaw']): return '偏航'
    if any(k in d for k in ['齿轮', '齿箱', 'gearbox']): return '齿轮箱'
    if any(k in d for k in ['发电机', '定子', '转子', 'generator', 'stator']): return '发电机'
    if any(k in d for k in ['安全', '紧停', '急停', 'safety', 'estop']): return '安全链'
    if any(k in d for k in ['暴风', '大风', 'storm', 'high wind']): return '暴风'
    if any(k in d for k in ['限功率', 'SCADA', '停止按钮', 'curtail']): return '控制/限功率'
    return '其他'


# ----------------------------- A. 风况 → 健康 -----------------------------
def map_wind_health(df, *, comp_cols, load_cols, wind_col, turbine_col,
                    sector_col=None, month_col=None, gen_mask=None,
                    min_n=5000, min_per_turbine=1000):
    """A映射: 每健康信号 H ~ load + W(风速[+扇区+季节]) 分层回归.

    产出 per 信号: r2_load / r2_full / dr2_wind(风况在载荷之外边际) +
      per台 **带符号z** 内禀残差(护栏①: +热/−冷, 非|z|) + 条件性(台内扇区残差极差z).
    load_cols = 载荷协变量列表(如 [功率,转速,舱外温]); wind_col=风速; sector_col=机舱位置(可选);
    gen_mask = 发电态布尔(护栏④: 前提门在调用方, 传已剥离的发电态). 返回 dict{signals:{...}}.
    """
    d = df if gen_mask is None else df[gen_mask]
    d = d.copy()
    out = {}
    for H in comp_cols:
        need = [H] + list(load_cols) + [wind_col] + ([sector_col] if sector_col else []) \
               + ([month_col] if month_col else [])
        sub = d[[c for c in need if c in d] + [turbine_col]].dropna(
            subset=[c for c in need if c in d])
        if len(sub) < min_n:
            out[H] = {'n': int(len(sub)), 'verdict': 'INSUFFICIENT'}; continue
        y = sub[H].to_numpy(float)
        ones = np.ones(len(sub))
        # load-only 设计
        Xl_parts = [ones] + [sub[c].to_numpy(float) for c in load_cols]
        Xl = np.column_stack(Xl_parts)
        # full = load + 风速 [+ sin/cos扇区] [+ sin/cos季节]
        Xf_parts = list(Xl_parts) + [sub[wind_col].to_numpy(float)]
        if sector_col:
            th = np.deg2rad(sub[sector_col].to_numpy(float))
            Xf_parts += [np.sin(th), np.cos(th)]
        if month_col:
            mth = np.deg2rad((sub[month_col].to_numpy(float) - 1) / 12 * 360)
            Xf_parts += [np.sin(mth), np.cos(mth)]
        Xf = np.column_stack(Xf_parts)
        yl = _ols_pred(Xl, y); yf = _ols_pred(Xf, y)
        r2l, r2f = _r2(y, yl), _r2(y, yf)
        # (b) per台 带符号z 内禀残差 (护栏①)
        sub = sub.assign(_rl=y - yl)
        pm = sub.groupby(turbine_col)['_rl'].agg(['mean', 'size'])
        pm = pm[pm['size'] >= min_per_turbine]
        mu, sd = pm['mean'].mean(), pm['mean'].std()
        pm['z'] = (pm['mean'] - mu) / sd if sd and sd > 0 else np.nan
        pm_sorted = pm.reindex(pm['z'].abs().sort_values(ascending=False).index)
        intrinsic = [{'tid': str(t), 'z': round(float(pm.loc[t, 'z']), 2),          # ★带符号
                      'resid_degC': round(float(pm.loc[t, 'mean']), 2),
                      'dir': ('热' if pm.loc[t, 'z'] > 0 else '冷')}
                     for t in pm_sorted.index[:5]]
        # (c) 条件性缺陷 (台内扇区残差极差 z)
        conditional = []
        if sector_col:
            sub = sub.assign(_rf=y - yf, _sec=np.clip(np.floor(sub[sector_col].to_numpy(float) / 45), 0, 7))
            ms = sub.groupby([turbine_col, '_sec'])['_rf'].mean().unstack()
            span = (ms.max(axis=1) - ms.min(axis=1))
            smu, ssd = span.mean(), span.std()
            if ssd and ssd > 0:
                sz = ((span - smu) / ssd).sort_values(ascending=False)
                conditional = [{'tid': str(t), 'span_z': round(float(sz[t]), 2),
                                'span_degC': round(float(span[t]), 2)} for t in sz.index[:3]]
        out[H] = {
            'n': int(len(sub)),
            'r2_load': round(r2l, 3), 'r2_full': round(r2f, 3),
            'dr2_wind_marginal': round(r2f - r2l, 4),
            'intrinsic_signed_z': intrinsic,           # 护栏①: 带符号
            'conditional_defect': conditional,
        }
    return {'mapping': 'A_wind_health', 'signals': out,
            'guards': ['signed_z(非|z|)', 'gen_mask前提门在调用方']}


# ----------------------------- B. 风况 → 故障 -----------------------------
def map_wind_fault(df, faults, *, time_col, turbine_col, wind_col, gen_col,
                   fault_time_col, fault_turbine_col, fault_desc_col,
                   subsystem_fn=None, tol_min=30, min_events=30,
                   episode_gap_s=3600, wsbins=(0, 3, 5, 7, 9, 11, 13, 25)):
    """B映射: 故障激活时刻回查风况 → P(W|F) vs 基线P(W) + KS + **状态混杂检验(护栏②)**.

    护栏②硬核: 每子系统报发电态占比(gen_frac); gen_frac≪基线 → 负风偏移=停机态触发伪相关,
    非"低风驱动"(不控制会误报). 返回 dict{subsystems:{sys:{wind_shift,ks,gen_frac,class,...}}}.
    class ∈ {wind_driven, wind_independent, state_confounded, weak}.
    """
    subsystem_fn = subsystem_fn or default_subsystem
    sc = df[[turbine_col, time_col, wind_col, gen_col]].copy()
    sc[time_col] = pd.to_datetime(sc[time_col], errors='coerce')
    sc = sc.dropna(subset=[time_col]).sort_values(time_col)
    w_base = sc[wind_col].dropna().to_numpy(float)
    base_mean = float(np.nanmean(w_base))
    base_gen = float(sc[gen_col].mean())

    fa = faults[[fault_turbine_col, fault_time_col, fault_desc_col]].copy()
    fa[fault_time_col] = pd.to_datetime(fa[fault_time_col], errors='coerce')
    fa = fa.dropna(subset=[fault_time_col, fault_turbine_col])
    lo, hi = sc[time_col].min(), sc[time_col].max()
    fa = fa[(fa[fault_time_col] >= lo) & (fa[fault_time_col] <= hi)]
    fa['_sys'] = fa[fault_desc_col].map(subsystem_fn)

    # per台 merge_asof 回查风况+态
    merged = []
    for m, g in fa.groupby(fault_turbine_col):
        scm = sc[sc[turbine_col] == m][[time_col, wind_col, gen_col]].sort_values(time_col)
        if len(scm) == 0:
            continue
        j = pd.merge_asof(g.sort_values(fault_time_col), scm,
                          left_on=fault_time_col, right_on=time_col,
                          direction='nearest', tolerance=pd.Timedelta(f'{tol_min}min'))
        merged.append(j)
    if not merged:
        return {'mapping': 'B_wind_fault', 'subsystems': {}, 'note': '无匹配'}
    fa2 = pd.concat(merged, ignore_index=True)
    bins = np.asarray(wsbins, float)

    subs = {}
    for sysn, g in fa2.groupby('_sys'):
        gg = g.dropna(subset=[wind_col])
        if len(gg) < min_events:
            subs[sysn] = {'n_matched': int(len(gg)), 'class': 'INSUFFICIENT_n'}; continue
        # episode 去重 (同台>gap 算新)
        gg = gg.sort_values([fault_turbine_col, fault_time_col])
        gap = gg.groupby(fault_turbine_col)[fault_time_col].diff().dt.total_seconds()
        ep = gg[(gap.isna()) | (gap > episode_gap_s)]
        w_f = ep[wind_col].to_numpy(float)
        gen_frac = float(ep[gen_col].mean()) if gen_col in ep else np.nan
        shift = float(np.nanmean(w_f) - base_mean)
        ks = ks_stat(w_f, w_base)
        # 分类 (护栏②): 状态混杂优先 — 故障在停机态触发(gen_frac≪基线)本身即混杂,
        # 其任何风况关联(不论正负偏移)都被运行态污染, 不可归"风驱动"。此判据不依赖偏移符号。
        ks_ok = ks is not None and ks == ks
        if gen_frac == gen_frac and gen_frac < base_gen - 0.2:
            cls = 'state_confounded'
        elif ks_ok and abs(ks) < 0.15 and abs(shift) < 1.0:
            cls = 'wind_independent'            # 风况≈基线 = 调度/时间驱动
        elif shift > 1.0 or (ks_ok and ks > 0.3):
            cls = 'wind_driven'                # 风况显著偏离 (且非停机态混杂)
        else:
            cls = 'weak'
        subs[sysn] = {
            'n_episodes': int(len(ep)), 'wind_mean_at_fault': round(float(np.nanmean(w_f)), 2),
            'wind_shift': round(shift, 2), 'ks_vs_baseline': round(ks, 3) if ks == ks else None,
            'gen_frac': round(gen_frac, 2), 'class': cls,
        }
    return {'mapping': 'B_wind_fault', 'baseline_wind_mean': round(base_mean, 2),
            'baseline_gen_frac': round(base_gen, 2), 'subsystems': subs,
            'guards': ['state_confound(gen_frac)', 'episode去重(削re-fire)']}


# ----------------------------- C. 风况 → 可靠性 -----------------------------
def map_wind_reliability(df, faults, *, time_col, turbine_col, wind_col, gen_col, stop_col,
                         fault_time_col, fault_turbine_col, fault_type_col=None,
                         fault_kind_generating=None, cutin=3.5, episode_gap_s=3600,
                         dt_minutes=10):
    """C映射: 台间 可用率/故障率/MTBF vs 风况暴露 回归 + **运行时长归一(护栏③)**.

    护栏③硬核: 故障率按 gen_hours 归一(每千运行时) + 报 暴露vs运行时长 相关
    (若强正 → 暴露=活跃度, 故障相关是logging伪). 返回 dict{tests, per_turbine, confound}.
    """
    sc = df[[turbine_col, time_col, wind_col, gen_col, stop_col]].copy()
    sc[time_col] = pd.to_datetime(sc[time_col], errors='coerce')
    sc = sc.dropna(subset=[time_col])
    rows = []
    for m, g in sc.groupby(turbine_col):
        windy = g[wind_col] >= cutin
        nw = int(windy.sum())
        gen_h = float(g[gen_col].sum()) * dt_minutes / 60.0
        avail = float((g.loc[windy, gen_col] == 1).mean()) if nw > 0 else np.nan
        fdown = float(((g[stop_col] == 1) & windy).sum()) / max(nw, 1)
        rows.append({'tid': str(m), 'ws_exposure': round(float(g[wind_col].mean()), 3),
                     'availability': round(avail, 4), 'forced_down_frac': round(fdown, 4),
                     'gen_hours': round(gen_h, 1)})
    rel = pd.DataFrame(rows)

    # 故障 episode (可选只故障型) per台
    fa = faults[[fault_turbine_col, fault_time_col] + ([fault_type_col] if fault_type_col else [])].copy()
    fa[fault_time_col] = pd.to_datetime(fa[fault_time_col], errors='coerce')
    fa = fa.dropna(subset=[fault_time_col, fault_turbine_col])
    if fault_type_col and fault_kind_generating:
        fa = fa[fa[fault_type_col] == fault_kind_generating]
    fa = fa.sort_values([fault_turbine_col, fault_time_col])
    gap = fa.groupby(fault_turbine_col)[fault_time_col].diff().dt.total_seconds()
    ep = fa[(gap.isna()) | (gap > episode_gap_s)]
    span_days = (fa[fault_time_col].max() - fa[fault_time_col].min()).total_seconds() / 86400 if len(fa) else np.nan
    ec = ep.groupby(fault_turbine_col).size()
    rel['fault_episodes'] = rel['tid'].map(lambda t: int(ec.get(t, 0)))
    rel['fault_per_1kh'] = rel['fault_episodes'] / rel['gen_hours'].replace(0, np.nan) * 1000
    rel['mtbf_days'] = (span_days / rel['fault_episodes'].replace(0, np.nan)).round(1)

    # 暴露 vs 可靠性 台间 Spearman
    tests = {}
    for metric in ['availability', 'forced_down_frac', 'fault_episodes', 'fault_per_1kh', 'mtbf_days']:
        r, p, n = spearman(rel['ws_exposure'], rel[metric])
        tests[metric] = {'spearman_r': round(r, 3) if r == r else None,
                         'p': round(p, 3) if p == p else None, 'n': n,
                         'sig': bool(p < 0.05) if p == p else False}
    # 护栏③: 暴露 vs 运行时长 (排除"暴露=活跃度")
    er, ep_, en = spearman(rel['ws_exposure'], rel['gen_hours'])
    confound = {'exposure_vs_genhours_r': round(er, 3) if er == er else None,
                'p': round(ep_, 3) if ep_ == ep_ else None,
                'is_activity_proxy': bool(ep_ < 0.05 and abs(er) > 0.4) if ep_ == ep_ else False}
    return {'mapping': 'C_wind_reliability', 'n_turbines': len(rel),
            'exposure_range': [round(float(rel['ws_exposure'].min()), 2), round(float(rel['ws_exposure'].max()), 2)],
            'tests': tests, 'confound_activity': confound,
            'per_turbine': rel.sort_values('availability').to_dict('records'),
            'guards': ['fault率按gen_hours归一', '暴露vs运行时长排除活跃度伪']}


# ----------------------------- E. 风况 → 载荷 → 疲劳/DEL -----------------------------
def rainflow_del(x, m=4.0):
    """雨流计数(ASTM 4点) → 相对 DEL = (Σ rangeᵢ^m)^(1/m). x=去趋势载荷代理序列(塔振/应变).
    m=Wöhler S-N 斜率(焊接钢塔≈4). 相对量(无气弹标定→只台内/相对比)."""
    x = np.asarray(x, float); x = x[np.isfinite(x)]
    if len(x) < 8:
        return np.nan
    d = np.diff(x)
    tp = np.concatenate([[x[0]], x[np.where(np.diff(np.sign(d)) != 0)[0] + 1], [x[-1]]])
    # tp 恒 ≥2 (首末); 平/单调信号 → 残余半循环处理(平=0, 斜坡=range), 非 nan
    ranges, stack = [], []
    for p in tp:
        stack.append(p)
        while len(stack) >= 4:
            s1, s2, s3 = stack[-3], stack[-2], stack[-1]
            if abs(s2 - s3) >= abs(s1 - s2):
                ranges.append(abs(s1 - s2)); del stack[-3]; del stack[-2]
            else:
                break
    ranges += [abs(stack[i] - stack[i + 1]) for i in range(len(stack) - 1)]
    return float(np.sum(np.asarray(ranges) ** m) ** (1.0 / m)) if ranges else 0.0


def map_wind_load(windows, *, del_col, ti_col, load_cols, turbine_col,
                  shear_col=None, min_n=2000, min_per_turbine=50):
    """E映射: 载荷剂量(DEL) ~ load + W(TI[,切变]). 台内自比消传感器增益污染(护栏⑤).

    windows = per-窗 DataFrame (每台每10min: del_col=雨流DEL, ti_col=湍流强度, load_cols=风速/转速/功率,
      shear_col=切变α可选). **数据门硬**: 需秒级塔振/应变→DEL(SCADA-10min无, 仅有高频载荷通道的场可跑, 如阜山)。
    护栏⑤(本映射新增): 加速度计/应变无标定 → 跨台绝对DEL不可比 → **台内自比(z-within)是唯一可信轴**;
      台间报但标 gain_confounded=INSUFFICIENT。返回 dict{within_turbine(主), pooled(增益污染), inter_turbine}.
    """
    d = windows.dropna(subset=[del_col, ti_col] + list(load_cols)).copy()
    d = d[d[del_col] > 0]
    if len(d) < min_n:
        return {'mapping': 'E_wind_load_DEL', 'n': int(len(d)), 'verdict': 'INSUFFICIENT'}
    d['_y'] = np.log(d[del_col].to_numpy(float))    # DEL 长尾 → log
    ones = np.ones(len(d))

    def z(col):
        return (pd.to_numeric(d[col], errors='coerce') - pd.to_numeric(d[col], errors='coerce').mean()) \
               / pd.to_numeric(d[col], errors='coerce').std()

    def zwin(col):
        g = d.groupby(turbine_col)[col]
        return ((d[col] - g.transform('mean')) / g.transform('std')).fillna(0).to_numpy()

    def _r2(y, X):
        b, *_ = np.linalg.lstsq(X, y, rcond=None); yh = X @ b
        ss = np.nansum((y - yh) ** 2); st = np.nansum((y - np.nanmean(y)) ** 2)
        return (float(1 - ss / st) if st > 0 else np.nan), b

    y = d['_y'].to_numpy(float)
    load_pool = np.column_stack([ones] + [z(c).fillna(0).to_numpy() for c in load_cols])
    full_pool = np.column_stack([load_pool, z(ti_col).fillna(0).to_numpy()])
    r2lp, _ = _r2(y, load_pool); r2fp, _ = _r2(y, full_pool)
    # ★护栏⑤ 台内自比 (消增益) = 主轴
    yw = zwin('_y')
    load_w = np.column_stack([ones] + [zwin(c) for c in load_cols])
    full_w = np.column_stack([load_w, zwin(ti_col)])
    r2lw, _ = _r2(yw, load_w); r2fw, bw = _r2(yw, full_w)
    within = {'R2_load': round(r2lw, 3), 'R2_plus_TI': round(r2fw, 3),
              'dR2_TI_marginal': round(r2fw - r2lw, 4),
              'TI_coef': round(float(bw[-1]), 3), 'load_coefs': [round(float(bw[i + 1]), 3) for i in range(len(load_cols))]}
    if shear_col and shear_col in d:
        dm = d.dropna(subset=[shear_col])
        if len(dm) > min_n:
            def zwm(col):
                g = dm.groupby(turbine_col)[col]
                return ((dm[col] - g.transform('mean')) / g.transform('std')).fillna(0).to_numpy()
            ym = zwm('_y'); om = np.ones(len(dm))
            lw = np.column_stack([om] + [zwm(c) for c in load_cols] + [zwm(ti_col)])
            fw = np.column_stack([lw, zwm(shear_col)])
            r2l2, _ = _r2(ym, lw); r2f2, bf = _r2(ym, fw)
            within['dR2_shear_marginal'] = round(r2f2 - r2l2, 4)
            within['shear_coef'] = round(float(bf[-1]), 3)

    # 台间 (增益污染 → INSUFFICIENT)
    per = d.groupby(turbine_col).agg(del_mean=(del_col, 'mean'), ti_mean=(ti_col, 'mean'), n=(del_col, 'size'))
    per = per[per['n'] >= min_per_turbine]
    r_ti, p_ti, n_ti = spearman(per['ti_mean'], per['del_mean'])
    return {
        'mapping': 'E_wind_load_DEL', 'n_windows': int(len(d)), 'n_turbines': int(d[turbine_col].nunique()),
        'within_turbine': within,          # 主轴 (护栏⑤消增益)
        'pooled_gain_confounded': {'R2_load': round(r2lp, 3), 'dR2_TI': round(r2fp - r2lp, 4)},
        'inter_turbine': {'TI_exposure_vs_DEL_r': round(r_ti, 3) if r_ti == r_ti else None,
                          'p': round(p_ti, 3) if p_ti == p_ti else None,
                          'verdict': 'INSUFFICIENT_gain_confounded (绝对DEL跨台不可比, 须加速度计标定)'},
        'guards': ['台内自比z-within消传感器增益(护栏⑤)', 'log(DEL)长尾稳', '相对DEL非绝对寿命(无气弹标定)'],
    }


# ----------------------------- D. 双向反馈 (设备 → 风资源) -----------------------------
def _compass_bearing(dx, dy):
    """罗盘方位 j→i (0=N 顺时针), dx=Δ东, dy=Δ北."""
    return float(np.degrees(np.arctan2(dx, dy)) % 360)


def map_device_to_wind(sector_metric, coords, *, turbine_col, sector_col, metric_col,
                       n_sectors=8, independent_axis=False, z_thresh=1.5, max_dist=None,
                       axis='fleet_relative'):
    """D 反馈映射: 设备异常方向/空间 pattern → 反推风资源特征(尾流). **护栏⑥ 防循环自证**.

    sector_metric: per-台per-扇区 设备指标 DataFrame (metric_col = **独立轴**: DEL/健康残差 by 风向扇区;
      sector_col = 风来向扇区 int [0, n_sectors)). coords: DataFrame[turbine_col, x, y] (米/相对, x东 y北).
    **护栏⑥ 硬门**: metric 必是独立于风况测量的设备响应(DEL/健康残差, 非 power/机舱风派生, 否则用设备欠发
      反推风况亏空=循环自证)。**调用方须显式 independent_axis=True 声明**, 否则拒(REJECTED_circular)。
    **axis (消盛行风混杂)**:
      'fleet_relative' (默认, 推荐): **两向可加残差** = metric − 台主效应 − 扇区主效应 + 总均 → 隔离 台×扇区
         交互(=尾流)。台主效应消增益/offset, 扇区主效应消盛行风+风级共模(盛行风是扇区共模被完全吸收 → 假阳~0)。
         正定量(DEL)自动 log 转乘性增益为加性。**代价**: 尾流须相对同侪可辨(同扇区多台齐受尾流→交互被主效应吸收→recall折损)。
      'within_turbine' (legacy): 单步台内跨扇区 z。**对推力驱动量(fore-aft DEL)被盛行风严重骗**(卡8.6假阳/场)。
    方法: 异常(台,扇区)(z>thresh) → 尾流几何检验: 异常扇区 S(风从S来) → 上风邻机 i(bearing j→i ∈ S, 取最近) → i 尾流 j。
    返回: {wake_pairs, n_anomalous_sectors, independence_ok, axis, ...}.
    """
    if not independent_axis:
        return {'mapping': 'D_device_to_wind', 'verdict': 'REJECTED_circular',
                'note': '护栏⑥: metric 未声明独立轴 → 用设备(可能power派生)反推风况=循环自证, 拒。'
                        '须传独立轴(DEL/健康残差) + independent_axis=True。'}
    sw = 360.0 / n_sectors
    cd = {str(r[turbine_col]): (float(r['x']), float(r['y'])) for _, r in coords.iterrows()}
    df = sector_metric.copy()
    df[turbine_col] = df[turbine_col].astype(str)
    if axis == 'fleet_relative':
        # 两向可加残差 = metric − 台主效应 − 扇区主效应 + 总均 → 隔离 台×扇区交互(=尾流)。
        # 台主效应=增益/offset(消); 扇区主效应=盛行风+风级共模(消, 盛行风被完全吸收 → 假阳~0)。
        x = df[metric_col].to_numpy(float)
        df['_v'] = np.log(np.clip(x, 1e-9, None)) if np.all(x > 0) else x   # 乘性增益(DEL)→log转加性
        gm = float(np.nanmean(df['_v']))
        tm = df.groupby(turbine_col)['_v'].transform('mean')
        sm2 = df.groupby(sector_col)['_v'].transform('mean')
        resid = df['_v'] - tm - sm2 + gm
        sd = float(resid.std())
        df['_z'] = resid / sd if sd > 0 else np.nan
    else:   # legacy within_turbine: 单步台内跨扇区 z (被盛行风骗)
        gt = df.groupby(turbine_col)[metric_col]
        sd1 = gt.transform('std')
        df['_z'] = np.where(sd1 > 0, (df[metric_col] - gt.transform('mean')) / sd1, np.nan)
    anom = df[df['_z'] > z_thresh]
    pairs = {}
    for _, row in anom.iterrows():
        j = str(row[turbine_col]); s = int(row[sector_col])
        if j not in cd:
            continue
        xj, yj = cd[j]; lo, hi = s * sw, (s + 1) * sw
        best = None
        for i, (xi, yi) in cd.items():
            if i == j:
                continue
            dx, dy = xi - xj, yi - yj
            dist = (dx * dx + dy * dy) ** 0.5
            if max_dist and dist > max_dist:
                continue
            brg = _compass_bearing(dx, dy)
            if (lo <= brg < hi) and (best is None or dist < best[1]):    # 上风邻机方位∈异常扇区, 取最近
                best = (i, dist, brg)
        if best:
            key = (j, s)
            pairs[key] = {'upwind': best[0], 'downwind': j, 'sector': s,
                          'z': round(float(row['_z']), 2), 'dist': round(best[1], 1),
                          'bearing_j_to_i': round(best[2], 1)}
    wake_pairs = sorted(pairs.values(), key=lambda p: -p['z'])
    return {
        'mapping': 'D_device_to_wind', 'n_turbines': int(df[turbine_col].nunique()),
        'n_anomalous_sectors': int(len(anom)), 'n_wake_pairs': len(wake_pairs),
        'wake_pairs': wake_pairs, 'independence_ok': True, 'axis': axis,
        'guards': ['独立轴防循环自证(护栏⑥, 拒power派生)',
                   ('两步z消盛行风混杂(fleet_relative默认: 台内→扇区内跨台, 卡实测假阳8.6→~0)'
                    if axis == 'fleet_relative' else
                    '⚠legacy within_turbine轴: 推力驱动量被盛行风骗(卡8.6假阳/场), 宜换 fleet_relative'),
                   '几何尾流检验(上风邻机方位∈异常扇区)', '结论=候选, 须performance-deficit交叉+现场核'],
    }