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

# -*- coding: utf-8 -*-
"""validate_pair.py — 配对金标准验证 harness (Tier-1 健康 yes/no vs 故障告警 + CMH 跨场池化).

把一次性脚本 (changyi/yuelinbai/feicheng) 泛化成 farm-agnostic 模块: 新配对场登记进
PAIRED_FARMS → `python -m src.sop.validate_pair` 一命令跑各场 Tier-1 2×2 + CMH 池化。
= 数据获取瓶颈的操作化件 (docs/业主数据需求清单_v1.md §4): 每多一个配对场, Tier-1 跨场立/证伪近一步。

★纪律 (验证程序五次实证):
- **严格主传动链热口径** (剔桨叶/变桨润滑污染, 肥城 raw OR=0.17→严格 3.33)。口径摊开可审 (THERMAL_KW)。
- **红台非忙台** 已证伪共因; 但 **filter-search 陷阱**: 换口径/分界结论会松动 → 报多口径不挑最优。
- **跨场需 n≥3** 才算 systemic (CLAUDE.md §跨场); 单场显著 ≠ 跨场立。CMH 附单场驱动检查。

纯核心 (is_thermal / red_from_findings / tier1_table / cmh) 无 I/O, 可测; per-farm 加载在 driver。
"""

try:
    from app_common.app_common_guanlan.api import install_root as _install_root
except ImportError:                     # 理论不可达；包结构异常时回退到按位置上跳
    from pathlib import Path as _P
    def _install_root(_f): return _P(_f).resolve().parents[4]
import re
import numpy as np
from scipy.stats import fisher_exact

# 严格主传动链热告警口径 (可审; 剔桨叶润滑/电气/环境/控制)
THERMAL_KW = {
    'inc': ['齿轮箱', '齿轮油', '主轴承', '低速轴', '高速轴', '驱动端', '自由端', '非驱动端',
            '定子', '发电机', '变流器', '轴承', '绕组'],
    'temp': ['温', '冷却'],
    'exc': ['桨', '润滑油位', '润滑出错', '润滑压力', '偏航', '电网', '电压', '电流', '通讯',
            '复位', '待机', '环境', '风速', '油位', '液压', '刹车', '制动'],
}


def is_thermal(desc, kw=THERMAL_KW):
    """告警描述 → 是否严格主传动链热 (纯). 剔 exc 污染, 需 inc 部件 ∧ temp 热词。"""
    s = str(desc)
    if any(e in s for e in kw['exc']):
        return False
    return any(k in s for k in kw['inc']) and any(t in s for t in kw['temp'])


def red_from_findings(findings_doc, exclude_dead_verdict='参考'):
    """findings.json → agent 温度标红台并集 (纯). 剔死列/低于候选门的 finding。"""
    red = set()
    for f in (findings_doc.get('findings') if isinstance(findings_doc, dict) else findings_doc) or []:
        t = str(f.get('title', ''))
        if not ('残差' in t or 'NBM' in t or '温度' in t):
            continue
        v = str(f.get('verdict', '')).strip('【】')
        if v in ('候选', '准定论·预警', '定论'):     # 剔参考/INSUFFICIENT (含死列)
            for m in (f.get('affected_machines') or []):
                red.add(str(m).strip())
    return red


def tier1_table(thermal_counts, red_set, cutoff='median'):
    """per台严格热告警数 + 红台集 → Tier-1 2×2 + Fisher (纯).
    thermal_counts: {tid: count} 全台; red_set: set[tid]; cutoff: median|presence|mean.
    返回 {table:[[a,b],[c,d]], or, p, cut, n} (a=红×热高)。"""
    tids = list(thermal_counts)
    vals = np.array([thermal_counts[t] for t in tids], dtype=float)
    if cutoff == 'presence':
        thr = 0.0
    elif cutoff == 'mean':
        thr = float(vals.mean())
    else:
        thr = float(np.median(vals))
    tab = [[0, 0], [0, 0]]
    for t in tids:
        r = 0 if str(t) in {str(x) for x in red_set} else 1
        h = 0 if thermal_counts[t] > thr else 1
        tab[r][h] += 1
    try:
        orr, p = fisher_exact(tab)
    except ValueError:
        orr, p = float('nan'), 1.0
    return {'table': tab, 'or': orr, 'p': p, 'cutoff': f'{cutoff}={thr:.4g}', 'n': len(tids)}


def cmh(tables):
    """Cochran-Mantel-Haenszel 跨场池化 (纯, 手算). tables: list[[[a,b],[c,d]]].
    返回 {mh_or, chi2, p, per_stratum_or, single_farm_driven}。防 statsmodels 依赖。"""
    from scipy.stats import chi2 as chi2dist
    num = den = 0.0
    a_sum = e_sum = v_sum = 0.0
    per = []
    for T in tables:
        (a, b), (c, d) = T
        n = a + b + c + d
        if n == 0:
            continue
        num += a * d / n
        den += b * c / n
        a_sum += a
        e_sum += (a + b) * (a + c) / n
        if n > 1:
            v_sum += (a + b) * (c + d) * (a + c) * (b + d) / (n * n * (n - 1))
        orr = (a * d) / (b * c) if b * c else float('inf')
        per.append(round(orr, 2))
    mh_or = num / den if den else float('inf')
    chi2v = (abs(a_sum - e_sum) - 0.5) ** 2 / v_sum if v_sum else 0.0  # 连续性校正
    p = float(chi2dist.sf(chi2v, 1))
    # 单场驱动检查: 去掉每一场后是否仍 OR>1 方向 (够粗)
    return {'mh_or': round(mh_or, 2), 'chi2': round(chi2v, 3), 'p': round(p, 4),
            'per_stratum_or': per, 'n_strata': len(per),
            'note': f'{len(per)}场; systemic门槛n≥3' + (' (未达)' if len(per) < 3 else '')}


# ══ Tier-2 故障锁定 vs 检修记录 (先兆核; 审-derived 方法烤进纯函数) ══════════
# changyi 五次过度包装教训烤进 lead_signal: 仅运行日残差 vs 对照(非raw分位) + 恒偏vs发展 + 样本闸。

# 缺陷名 → 域: temp(温度可探) / blind(SCADA温度盲, 判 INSUFFICIENT 非漏检) / other
DEFECT_SCOPE_KW = {
    'temp': ['水冷', '冷却', '轴承', '齿轮', '油温', '油泵', '过热', '碳刷', '发电机', '定子', '温度', '绕组'],
    'blind': ['叶片', '变桨', '桨距', '不变桨', '熔断', '690v', '电缆', '绝缘', '接触器'],
}


def defect_scope(name):
    """缺陷名 → temp|blind|other (纯). blind = SCADA温度看不见, agent 判 INSUFFICIENT 是对的非漏检。"""
    s = str(name).lower()
    if any(k in s for k in DEFECT_SCOPE_KW['blind']):
        return 'blind'
    if any(k in str(name) for k in DEFECT_SCOPE_KW['temp']):
        return 'temp'
    return 'other'


def lead_signal(resid_window, control_resids, z_dev=1.0, dev_frac=0.5, min_points=8):
    """故障锁定先兆核 (纯, 审-derived; **只出客观描述, 不下命中**).

    resid_window: 该台该通道 发生前窗的 NBM 残差序列 (仅运行日, 按时序; 季节已在 NBM 扣)。
    control_resids: 对照台同通道 per-台均值残差 list (基线云)。
    kind 判据 (烤进 changyi 五次教训):
      insufficient: n<min_points (A19/A30: 发生前多停机, 运行日点极少, 趋势外推自空气)。
      no_signal: z_vs_control < z_dev (A25: 落对照云内, raw分位高但残差无特异)。
      constant_offset: z≥z_dev ∧ 窗首≈窗尾 (A28/A13: 恒偏一直高, 撞时间窗非预测先兆)。
      developing: z≥z_dev ∧ 窗尾显著高于窗首 (真发生前发展信号)。
    返回 {n, z_vs_control, early_mean, late_mean, kind}。
    """
    w = [x for x in resid_window if x == x]     # 去 NaN
    n = len(w)
    cr = [x for x in control_resids if x == x]
    cmean = float(np.mean(cr)) if cr else 0.0
    cstd = float(np.std(cr)) if cr and np.std(cr) > 0 else 1.0
    if n < min_points:
        return {'n': n, 'z_vs_control': None, 'early_mean': None, 'late_mean': None, 'kind': 'insufficient'}
    z = (float(np.mean(w)) - cmean) / cstd
    half = max(1, n // 2)
    early, late = float(np.mean(w[:half])), float(np.mean(w[half:]))
    if z < z_dev:
        kind = 'no_signal'
    elif (late - early) > dev_frac * cstd:
        kind = 'developing'
    else:
        kind = 'constant_offset'
    return {'n': n, 'z_vs_control': round(z, 2), 'early_mean': round(early, 2),
            'late_mean': round(late, 2), 'kind': kind}


def tier2_summary(results):
    """一场 Tier-2 逐缺陷 lead_signal 结果 → 客观汇总 (纯; 不下准确度定论)。
    results: list[{defect, scope, kind}]。返回各 kind 计数 + 明确"blind=出域非漏检"。"""
    from collections import Counter
    temp = [r for r in results if r.get('scope') == 'temp']
    kc = Counter(r.get('kind') for r in temp)
    return {
        'n_defects': len(results),
        'n_temp_detectable': len(temp),
        'n_blind': sum(1 for r in results if r.get('scope') == 'blind'),
        'developing': kc.get('developing', 0), 'constant_offset': kc.get('constant_offset', 0),
        'no_signal': kc.get('no_signal', 0), 'insufficient': kc.get('insufficient', 0),
        'note': 'developing=干净发生前信号; constant_offset=恒偏撞窗非预测; blind=SCADA盲区(agent判INSUFFICIENT非漏检)。'
                'developing才算"温度锁定命中", 且须独立审(filter/口径)。',
    }


# ── per-farm 加载适配 (I/O; 新配对场登记于此) ──────────────────────────────
def _bx(s):
    m = re.search(r'B(\d+)号', str(s))
    return m.group(1) if m else re.sub(r'\D', '', str(s))[:3]


PAIRED_FARMS = {
    'yuelinbai': {
        'findings': 'outputs/yuelinbai/sop/findings.json',
        'fault_glob': '/Volumes/WINDDATA/DATA2/华电/华电山东/岳林柏风电场-山东-华电/收资数据/故障记录/**/*',
        'enc': 'gbk', 'desc_col': '状态码描述', 'tid_col': '风机名', 'tid_norm': _bx,
    },
    'feicheng': {
        'findings': 'outputs/feicheng/sop/findings.json',
        'fault_path': '/Volumes/WINDDATA/DATA2/华电/华电山东/历史/收资数据-华电山东/肥城风电场/肥城上汽数据/故障记录_20250101_20251124.csv',
        'enc': 'gbk', 'desc_col': '状态码描述', 'tid_col': '风机名', 'tid_norm': lambda s: str(int(float(s))) if str(s).replace('.', '').isdigit() else str(s),
    },
}


def _load_farm(cfg):
    """加载一场: findings→红台 + 故障记录→per台严格热计数。返回 (red, thermal_counts)。"""
    import json
    import glob
    import pandas as pd
    from pathlib import Path
    root = _install_root(__file__)
    fd = json.loads((root / cfg['findings']).read_text(encoding='utf-8'))
    red = {cfg['tid_norm'](m) for m in red_from_findings(fd)}
    # 故障记录 (单文件或 glob 多文件)
    if cfg.get('fault_path'):
        files = [cfg['fault_path']]
    else:
        files = [f for f in glob.glob(cfg['fault_glob'], recursive=True)
                 if f.lower().endswith(('.csv', '.xls', '.xlsx'))]
    parts = []
    for f in files:
        try:
            d = pd.read_excel(f) if f.lower().endswith(('.xls', '.xlsx')) else pd.read_csv(f, encoding=cfg['enc'])
            parts.append(d)
        except Exception:
            pass
    fa = pd.concat(parts, ignore_index=True)
    fa['tid'] = fa[cfg['tid_col']].map(cfg['tid_norm'])
    fa['th'] = fa[cfg['desc_col']].map(is_thermal)
    th = fa[fa['th']].groupby('tid').size()
    counts = {t: int(th.get(t, 0)) for t in fa['tid'].dropna().unique()}
    return red, counts


if __name__ == '__main__':
    # 控制台可能是 GBK(中文 Windows 代码页 936): 正文里的 ✔ ✗ ✅ ⚠ 这类字符编不出来会抛
    # UnicodeEncodeError, 脚本干成了事却以退出码 1 结束(同类坑见 src/console.py)。降级为 '?' 而不是崩;
    # 不用 import 是为了兼顾 python -m 与直接当脚本跑两种启动方式。
    import sys as _sys
    for _s in (_sys.stdout, _sys.stderr):
        try: _s.reconfigure(errors='replace')
        except Exception: pass
    # CLI: python -m src.sop.validate_pair  → 各场 Tier-1 + CMH 池化 (复现验证程序)
    import argparse
    ap = argparse.ArgumentParser()
    ap.add_argument('--cutoff', default='median', choices=['median', 'presence', 'mean'])
    a = ap.parse_args()
    tables = []
    print(f"=== Tier-1 配对验证 harness (严格主传动链热, cutoff={a.cutoff}) ===")
    for farm, cfg in PAIRED_FARMS.items():
        try:
            red, counts = _load_farm(cfg)
            r = tier1_table(counts, red, a.cutoff)
            tables.append(r['table'])
            (ta, tb), (tc, td) = r['table']
            print(f"  {farm:12s}: 红[{ta},{tb}] 非红[{tc},{td}] OR={r['or']:.2f} p={r['p']:.3f} (红{len(red)}台/共{r['n']})")
        except Exception as e:
            print(f"  {farm:12s}: 加载失败 {str(e)[:50]}")
    if len(tables) >= 2:
        pool = cmh(tables)
        print(f"\n  CMH 池化: M-H OR={pool['mh_or']} chi2={pool['chi2']} p={pool['p']} | 各场OR={pool['per_stratum_or']} | {pool['note']}")
        if pool['n_strata'] < 3:
            print("  ⚠ 单场驱动风险 + 未达 n≥3 systemic 门槛 → 跨场未立 (需更多配对场, 见 业主数据需求清单)")