guanlan-src/scripts/rudong_model_run.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
r"""六层链 `model_run` 步 —— 窗索引 + 谱库 → **L6 过闸谱线表 + 层小结**（按口径重建）。

## 这一步是什么 / 不是什么（先读这段再读代码）

用户令 2026-09-18「/detail 各页面不许用旧版产出补, 必须基于输入数据重算」之后, CMS 报告里
唯一还空着的是**融合级**那一列, 它的上游就是本步与 `fusion` 步的两件产物:

```
configs:  'model_l6'      = outputs/<场>/m5_cms_tcm/model_run_l6.parquet    ← 从未随包（无标准答案）
          'model_summary' = outputs/<场>/m5_cms_tcm/model_run_summary.json
          'fusion'        = outputs/<场>/m5_cms_tcm/fusion_38.csv           ← 从未随包（无标准答案）
```

**不是**复现: 这两件从未随过包, 盘上没有样件可逐值对拍 ⇒ 无法证明"与振动线那套一致"（见
`docs/振动六层链_接口规格与缺口_v0.1.md` §3）。**是**按口径重建: 判据全部取自**包内既有实现**
（不在本器里另写一套阈值), 每条线都留下"过没过闸、按哪条闸定的级"的来路, 来历在
`model_run_summary.json` 的 `口径` 字段与产物台账里写明（`按口径重建, 无标准答案对拍`）。

## 判据来源（包内单一实现, 本器只做接线）

| 环节 | 包内实现 |
|---|---|
| 候选线（理论频率） | `reference/rudong/oem_scan_plan.json` 的 11 个部件 × {BPFI, BPFO, BSF}（分母已 `--verify-plan` 11/11 逐值验过） |
| 谱前提闸 G1–G8 | `src/sop/discriminators.py::spectral_line_gates`（分辨率/有线/峰稳/非固有/非整阶/选择性/fleet 离群） |
| 六枚举定级 | `src/sop/discriminators.py::vib_verdict_and_writeback`（L0 短路 / 机制未定不命名部件 / 无正样本锚封顶候选） |
| 正样本锚 | `reference/rudong/positive_anchors.json`（人工坐实的实物证据; 本器只读不写） |
| 轴系转频 | `oem_scan_plan.json` 各部件的 `shaft_Hz`（★只传该轴承自己那根轴 —— 传全机会把真线误杀） |

## 输出（消费者的列名是硬契约）

`model_run_l6.parquet` 的列由 `src/windcms/report_std.py::registry()` 与 `report.py` 决定:
`台 / 测点 / 线 / hz / 域 / 绝对量 / 单位 / xfleet / 选择性 / 占比pct / 绝对锚 / 定级`
（`台` 必须是 `f'{n}#'` 形态: registry 的 `_l6_consumed` 闭环断言会核对, 台号格式漂移当场报错）。

用法:
    python scripts/rudong_model_run.py                    # 默认取 M5_WINDOWS / 自动发现的窗
    python scripts/rudong_model_run.py --window w0316
    python scripts/rudong_model_run.py --turbines WTG01,WTG09   # 冒烟
    python scripts/rudong_model_run.py --dry-run          # 只报计划, 不写产物
"""
from __future__ import annotations

import argparse
import json
import pathlib
import sys
import time

import numpy as np
import pandas as pd

ROOT = pathlib.Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))
sys.path.insert(0, str(ROOT / 'scripts'))
from src import paths as P                                              # noqa: E402
from src.sop import discriminators as D                                 # noqa: E402
import rudong_tcm_oem_scan as OS                                        # noqa: E402
import spectra_lookup as SL                                             # noqa: E402

ANCHORS = ROOT / 'reference' / 'rudong' / 'positive_anchors.json'
COMP_CLASS = {                          # 部件 → 正样本锚的部件类 (reference/rudong/positive_anchors.json)
    'generator': 'generator_bearing', 'hs': 'gearbox_hs_bearing', 'ims': 'gearbox_ims_bearing',
    'main': 'main_bearing',
}
# ISO 10816 速度当量 (mm/s, 刚性支撑·>15kW 机组的黄/红带) — 唯一"已校准"的绝对锚
ISO = dict(yellow=4.5, red=11.0)
LINES = (('BPFI', 'BPFI_hz', '内圈'), ('BPFO', 'BPFO_hz', '外圈'), ('BSF', 'BSF_hz', '滚动体'))


def comp_class(component: str) -> str:
    for k, v in COMP_CLASS.items():
        if component.startswith(k):
            return v
    return 'unknown'


def anchors_of() -> list | None:
    if not ANCHORS.is_file():
        return None
    return json.loads(ANCHORS.read_text(encoding='utf-8')).get('anchors')


def windows_of(explicit: str | None) -> list[str]:
    import os
    if explicit:
        return [explicit]
    env = os.environ.get('M5_WINDOWS') or ''
    if env.strip():
        return [w for w in env.split(',') if w.strip()]
    m5 = P.m5()
    got = [p.parent.name for p in sorted((m5 / 'windows').glob('w[0-9][0-9][0-9][0-9]/index.parquet'))]
    if (m5 / 'tcm_index.parquet').is_file():
        got.insert(0, 'w0127')
    return got


_AVAIL: dict[str, set] = {}


def available(win: str) -> set:
    """该窗谱库里**实有**的 (机组, 测点, 测量) 组合 —— 计划里的测量名不一定每台都有。

    ★为什么必须先查: `spectra_lookup.load()` 查不到记录时是 `SystemExit`(它的 CLI 语义), 而
    SystemExit 继承 BaseException, 用 `except Exception` 兜不住 —— 首版跑到第一台就整体退出
    (实测: w0316 里 WTG01 没有 Env_6000_4000_850_Tr 的谱)。缺谱是**常见情形**, 该跳过并计数, 不是崩。
    """
    if win not in _AVAIL:
        mp = SL.meta_path(win)
        s = set()
        if mp.is_file():
            d = pd.read_parquet(mp, columns=['turbine', 'sensor', 'meas_name'])
            s = set(map(tuple, d[['turbine', 'sensor', 'meas_name']].drop_duplicates().values))
        _AVAIL[win] = s
    return _AVAIL[win]


def fleet_median(sensor: str, meas: str, turbines: list[str], win: str) -> np.ndarray | None:
    """同 (测点, 测量) 的**跨台中位谱** —— G4(机型固有) 与 G8(个体离群) 都靠它。"""
    got = []
    for t in turbines:
        if (t, sensor, meas) not in available(win):
            continue
        try:
            r = SL.load(win, t, sensor, meas)
        except BaseException:
            continue
        if not r:
            continue
        got.append(np.asarray(r['v'], dtype=float))
    if len(got) < 5:
        return None
    n = min(len(v) for v in got)
    return np.nanmedian(np.vstack([v[:n] for v in got]), axis=0)


def evaluate(win: str, turbines: list[str], comps: list[dict], anchors, dry=False) -> tuple[list, dict]:
    """逐 (机组 × 部件 × 特征频率) 过闸定级 → L6 行。"""
    rows, stats = [], dict(considered=0, passed=0, no_spectrum=0, by_gate={}, by_level={})
    fmed_cache: dict[tuple, np.ndarray | None] = {}
    avail = available(win)
    for c in comps:
        meas, sensor = c['shard'], c['sensor']
        key = (sensor, meas)
        if key not in fmed_cache:
            fmed_cache[key] = fleet_median(sensor, meas, turbines, win)
        fmed = fmed_cache[key]
        dx = float(c.get('dx_hz') or 0) or None
        shaft = c.get('shaft_Hz') or 0.0
        domain = '包络谱' if str(meas).lower().startswith('env') else '原始谱'
        cls = comp_class(c['component'])
        for t in turbines:
            if (t, sensor, meas) not in avail:
                stats['no_spectrum'] += 1
                continue
            try:
                r = SL.load(win, t, sensor, meas)
            except BaseException:
                r = None
            if r is None:
                stats['no_spectrum'] += 1
                continue
            y = np.asarray(r['v'], dtype=float)
            meta = r['row']
            _dx = float(r.get('x_delta') or dx or 0) or dx
            if not _dx:
                continue
            for label, col, cn in LINES:
                hz = c.get(col)
                if not hz or not np.isfinite(hz) or hz <= 0:
                    continue
                stats['considered'] += 1
                # 观测峰: 目标频率 ±2 bin 内取最大 (与谱前提闸 G7 的搜索窗同宽)
                tol = 2
                k = int(round(hz / _dx))
                if k - tol < 1 or k + tol >= len(y):
                    continue
                seg = y[k - tol:k + tol + 1]
                if not len(seg) or not np.isfinite(seg).max():
                    continue
                obs = float(np.nanmax(seg))
                x_fleet = None
                if fmed is not None and len(fmed) > k + tol:
                    base = float(np.nanmax(fmed[k - tol:k + tol + 1])) or np.nan
                    if np.isfinite(base) and base > 1e-12:
                        x_fleet = obs / base
                sel = None                     # 选择性 = 线 ÷ 同谱宽带本底 (同窗, 与 G6 同语义)
                lo, hi = max(int(0.25 * len(y)), 1), max(int(0.75 * len(y)), 2)
                bg = float(np.nanmedian(np.abs(y[lo:hi]))) if hi > lo else np.nan
                if np.isfinite(bg) and bg > 1e-12:
                    sel = obs / bg
                is_vel = str(meta.get('y_unit') or '') in ('m/s',)
                g = D.spectral_line_gates(y, _dx, float(hz), fleet_med=fmed, shaft_hz=(shaft,) if shaft else (),
                                          x_fleet=x_fleet, selectivity_samewin=sel, domain=domain,
                                          iso_ref_mms=(ISO['yellow'] if is_vel else None))
                gl = g.get('verdict') or '?'
                stats['by_gate'][gl] = stats['by_gate'].get(gl, 0) + 1
                if not g.get('passed'):
                    continue
                # 绝对量: 速度域取 mm/s 当量 (ISO 可直接比), 其余域取谱值原单位
                if is_vel:
                    abs_v, unit = obs * 1000.0, 'mm/s'
                else:
                    abs_v, unit = obs, str(meta.get('y_unit') or '')
                vw = D.vib_verdict_and_writeback(
                    gate_result=g, abs_value=abs_v, abs_unit=unit, x_fleet=x_fleet,
                    iso_yellow=ISO['yellow'] if is_vel else None,
                    iso_red=ISO['red'] if is_vel else None,
                    component_class=cls, anchors=anchors, n_evidence_types=1,
                    mechanism_confirmed=False)         # 机制未定 ⇒ 不命名部件 (硬规则 2)
                lvl = vw.get('level') or 'INSUFFICIENT'
                stats['passed'] += 1
                stats['by_level'][lvl] = stats['by_level'].get(lvl, 0) + 1
                # 解封判据: 锚层给出的最高可达级依据 (无实物锚 ⇒ 写明未解封) —— 消费端 report.py:213 与
                # plugins.verdict_gate 都会显示这一列, 缺列会让逐台页整页崩 (AttributeError, 2026-09-19 实逮)
                try:
                    ci = D.anchor_cap(cls, anchors, n_evidence_types=1, mechanism_confirmed=False)
                    unlock = f"锚 tier={ci.get('tier')}: {str(ci.get('reason'))[:60]}"
                except Exception as e:
                    unlock = f'锚不可判 ({type(e).__name__})'
                rows.append(dict(
                    台=f'{int(t[3:])}#', 测点=sensor, 线=f'{c["component"]}·{label}({cn})',
                    证据族=f'{domain}·{cn}', 解封判据=unlock,
                    hz=round(float(hz), 3), 域=domain,
                    绝对量=round(float(abs_v), 5), 单位=unit,
                    xfleet=(round(float(x_fleet), 2) if x_fleet is not None and np.isfinite(x_fleet) else None),
                    选择性=(round(float(sel), 2) if sel is not None and np.isfinite(sel) else None),
                    占比pct=None,
                    绝对锚=(f'ISO10816 黄{ISO["yellow"]}/红{ISO["red"]} mm/s' if is_vel else '—'),
                    定级=_consumer_level(lvl),
                    _闸=g.get('verdict'), _依据=' · '.join(vw.get('reasons') or [])[:220],
                    _写回=vw.get('writeback'), _窗=win))
    return rows, stats


def _consumer_level(lvl: str) -> str:
    """把 L6 六枚举映射到消费端 `_LINE_STATE` 认得的那一套 (report_std.py 的 line_state 表)。

    `候选·查偶发源` / `候选·记基线` 这类细分在消费端表里**没有条目** ⇒ 会被默认成"优秀"
    (静默降级成假好消息)。这里显式映射, 细分信息留在 `_依据` 里。
    """
    s = str(lvl)
    if s.startswith('确诊') or s == '定论':
        return '定论'
    if s.startswith('准定论'):
        return '准定论·预警'
    if s.startswith('候选·新发'):
        return '候选·新发'
    if s.startswith('候选·记基线'):
        return '候选·记基线'
    if s.startswith('候选'):
        return '候选'
    if s.startswith('参考·上升'):
        return '参考·上升'
    if s.startswith('参考'):
        return '参考'
    if s.startswith('撤回'):
        return '撤回'
    return 'INSUFFICIENT'


def l0_layer(win: str) -> dict:
    """L0 数据质量层 (按测点): 有解析错/过载/无记录 ⇒ 该测点结论不可采信 (report_std 会据它判"不可判")。

    只看**索引自己的字段**(parse_error/overload), 不加任何自造阈值 —— 这层的意义是"别把盲台当健康台"。
    """
    ix = P.m5() / ('tcm_index.parquet' if win == 'w0127' else f'windows/{win}/index.parquet')
    if not ix.is_file():
        return {}
    d = pd.read_parquet(ix, columns=[c for c in ('turbine', 'sensor_name', 'parse_error', 'overload')
                                     if c in pd.read_parquet(ix).columns])
    out = {}
    for (t, s), g in d.groupby(['turbine', 'sensor_name']):
        bad = int(g['parse_error'].notna().sum()) if 'parse_error' in g else 0
        ovl = int((g['overload'] == True).sum()) if 'overload' in g else 0          # noqa: E712
        if bad and bad == len(g):
            out[f'{t}|{s}'] = '不可用: 该测点记录全部解析失败'
        elif ovl > 0.5 * len(g):
            out[f'{t}|{s}'] = '不可用: 该测点过半记录过载'
        elif bad:
            out[f'{t}|{s}'] = f'可用(部分): {bad}/{len(g)} 条解析失败'
    return out


def main() -> int:
    ap = argparse.ArgumentParser(description='六层链 model_run 步: 窗索引+谱库 → L6 过闸谱线 (按口径重建)')
    ap.add_argument('--window', default=None)
    ap.add_argument('--turbines', default=None)
    ap.add_argument('--out', default=None)
    ap.add_argument('--dry-run', action='store_true')
    a = ap.parse_args()
    plan, comps = OS.load_plan()
    anchors = anchors_of()
    wins = windows_of(a.window)
    if not wins:
        print('[X] 没有可分析的窗 (m5_cms_tcm/windows/*/index.parquet 与 tcm_index.parquet 都不在位)')
        return 2
    turbines = ([f'WTG{x.strip().zfill(2)}' for x in a.turbines.split(',')] if a.turbines
                else [f'WTG{i:02d}' for i in range(1, 39)])
    print(f'model_run: 窗 {wins} · 机组 {len(turbines)} 台 · 部件 {len(comps)} 个 · '
          f'候选线 {len(comps) * len(LINES)}/台·窗 · 正样本锚 {len(anchors or [])} 条')
    if a.dry_run:
        for c in comps:
            print(f'   {c["component"]:24} {c["sensor"]:24} {c["shard"]:24} dx={c["dx_hz"]:<9} '
                  f'BPFI={c.get("BPFI_hz")} BPFO={c.get("BPFO_hz")} BSF={c.get("BSF_hz")} '
                  f'shaft={c.get("shaft_Hz")}')
        print('(dry-run, 未写产物)')
        return 0
    t0 = time.time()
    allrows, summary = [], {}
    for win in wins:
        rows, stats = evaluate(win, turbines, comps, anchors)
        allrows += rows
        summary[win] = stats
        print(f'   窗 {win}: 候选 {stats["considered"]} 条 → 过闸 {stats["passed"]} 条 · '
              f'定级 {stats["by_level"]} · 闸分布 {dict(sorted(stats["by_gate"].items(), key=lambda x: -x[1])[:5])}')
    m5 = P.m5()
    l6 = pd.DataFrame(allrows)
    out = pathlib.Path(a.out) if a.out else (m5 / 'model_run_l6.parquet')
    if len(l6):
        l6 = l6.sort_values(['台', '测点', '线']).reset_index(drop=True)
    # 消费端的硬契约 (report_std.registry 的 _l6_consumed 闭环 + report.py 的列单)
    # 列 = 消费端硬契约: registry()/report.py:504 的投影 + report.py:213 与 plugins.verdict_gate 的明细表
    COLS = ['台', '测点', '线', 'hz', '域', '绝对量', '单位', 'xfleet', '选择性', '占比pct', '绝对锚',
            '证据族', '解封判据', '定级', '_闸', '_依据', '_写回', '_窗']
    l6 = l6.reindex(columns=COLS) if len(l6) else pd.DataFrame(columns=COLS)
    bad = sorted(set(l6['台']) - {f'{i}#' for i in range(1, 39)}) if len(l6) else []
    if bad:
        print(f'[X] 台号不在 38 台清单里: {bad} —— 消费端 registry 的闭环断言会当场报错, 本器先拦')
        return 5
    out.parent.mkdir(parents=True, exist_ok=True)
    l6.to_parquet(out, index=False)
    sm_path = m5 / 'model_run_summary.json'
    L0 = l0_layer(wins[0])
    summ = dict(
        status='ok', built=time.strftime('%Y-%m-%d %H:%M:%S'), windows=wins, seconds=round(time.time() - t0, 1),
        n_lines=int(len(l6)), by_gate={w: s['by_gate'] for w, s in summary.items()},
        by_level={w: s['by_level'] for w, s in summary.items()},
        L0=L0,
        口径='按口径重建(**非**逐值对拍的复现 —— model_run_l6.parquet 从未随过包, 无标准答案, 见 '
             'docs/振动六层链_接口规格与缺口_v0.1.md §3)。候选线=reference/rudong/oem_scan_plan.json 的 11 个'
             '部件 × {BPFI,BPFO,BSF}(分母已 11/11 逐值验过); 过闸=src/sop/discriminators.py::spectral_line_gates '
             '(G1..G8); 定级=src/sop/discriminators.py::vib_verdict_and_writeback (L0 短路 / 机制未定不命名部件 / '
             '无正样本锚封顶候选); 轴系转频只传该轴承自己那根轴; ISO 绝对锚只对速度域生效',
        missing_note='峰值拾取仍是**本器自定**: 目标频率 ±2 bin 取最大 (与 G7 搜索窗同宽), 未做阶次跟踪/多记录合并 —— '
                     '振动线那套的"怎么取观测峰"口径仍缺 (docs §7.3), 故本件不可当"复现"用')
    sm_path.write_text(json.dumps(summ, ensure_ascii=False, indent=1), encoding='utf-8')
    print(f'已写 {P.rel(out)}: {len(l6)} 行')
    print(f'已写 {P.rel(sm_path)}: L0 {len(L0)} 条 · 窗 {len(wins)}')
    # 来源自登记 (谁算的谁登记)
    try:
        from src import derived_manifest as DM
        rels = {out.relative_to(P.out_root()).as_posix():
                'scripts/rudong_model_run.py (按口径重建, 无标准答案对拍; 判据=src/sop/discriminators.py)',
                sm_path.relative_to(P.out_root()).as_posix():
                'scripts/rudong_model_run.py (L0 层 + 口径说明)'}
        DM.record(P.out_root(), rels, by='rudong_model_run')
        print('   已自登记 → _derived_manifest.json')
    except Exception as e:
        print(f'   [i] 自登记跳过: {type(e).__name__}: {e}')
    return 0


if __name__ == '__main__':
    for _s in (sys.stdout, sys.stderr):
        try:
            _s.reconfigure(errors='replace')
        except Exception:
            pass
    sys.exit(main())