guanlan-src/scripts/rudong_tcm_oem_scan.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
r"""六层链 `oem_scan` 步 —— **逆向工程实现（分母已就位/分子待口径）**，2026-09-17。

## 现状：一半已经落地并逐值验过，另一半明确挡住

壳 `src/windcms/pipeline.py::analyze` 会调本脚本（无参数，读 `M5_WINDOWS` / `SCAN_BANDS` 环境变量），
产出特征频率扫描件（`bearing_freq_scan` / `gear_freq_scan` / `cage_slip_scan` / `oem_frequency_scan`）。

**已落地（分母 + 配置，逐值对拍通过）**：

1. **扫描配置**：每个部件用哪条谱、什么分辨率 —— `reference/rudong/oem_scan_plan.json`
   （11 个部件：`shard` 取 `Env_6000_4000_850_Tr` / `FFT_6000_Tr` / `FFT_2000_fast_Tr` / `FFT_62_Tr`，
   `dx_hz` 1.0625 / 0.9375 / 0.625 / 0.03875）。
2. **理论频率（分母）**：`option` 串 → `reference/rudong/oem_bearing_freqs.json` 条目，规则
   = 位置 token（`Pos.150R`→`150r`、`139.1/1060 f-576853.PRL`→`1391/1060`、`6338 M/C3`→`6338m/c3`）
   + 同位置多候选时按型号 token 定夺（NSK 284 vs FAG 273；NSK 137 vs FAG 151）。
   **实测：11/11 个部件的 `BPFI_hz`/`BPFO_hz`/`BSF_hz` 与随包样件逐值一致**（`--verify-plan`）。

**未落地（分子）**：`*_fleet_ratio` / `*_local_ratio` 要的是**观测频率 ÷ 理论频率**，而"在谱里怎么取观测峰"
（带宽 / 取什么量最大 / 是否插值 / 是否按 rpm 阶次跟踪 / 多记录怎么合并）是缺的口径
（见 `docs/向振动线取料单_v0.1.md` 第 1 项、`docs/振动六层链_接口规格与缺口_v0.1.md` §7）。

**因此本脚本默认行为是"响亮失败"**（rc=3），绝不写半成品：
在口径到位前，`windcms.py analyze --steps oem_scan` 会明确告诉你缺什么、去哪儿要，而不是产出一堆
只有一半是真的扫描件让人误以为重算成功了。

## 用法

    python scripts/rudong_tcm_oem_scan.py --verify-plan   # 分母+配置 与随包样件逐值对拍（不写盘）
    python scripts/rudong_tcm_oem_scan.py --plan          # 打印扫描计划（部件/测点/谱/分辨率/理论频率）
    python scripts/rudong_tcm_oem_scan.py --status        # 机器可读状态（供 chain_gap_check 用）
    python scripts/rudong_tcm_oem_scan.py                 # 真跑（口径未到位 ⇒ rc=3, 不写盘）

退出码: 0 成功/对拍通过 · 2 缺参考件 · 3 **峰值拾取口径未到位**（不写盘） · 5 对拍不一致
"""
from __future__ import annotations

import argparse
import json
import pathlib
import re
import sys

import pandas as pd

ROOT = pathlib.Path(__file__).resolve().parents[1]
sys.path.insert(0, str(ROOT))
from src import paths as P                                              # noqa: E402

PLAN = ROOT / 'reference' / 'rudong' / 'oem_scan_plan.json'
THEORY = ROOT / 'reference' / 'rudong' / 'oem_bearing_freqs.json'
# 口径闸: 峰值拾取口径到位后把它翻成 True, 本脚本才会真的写产物
PEAK_PICKING_READY = False
MISSING = ('峰值拾取口径: 理论频率附近**怎么取观测峰**(带宽 / 取什么量最大 / 是否插值 / '
           '是否按 rpm 阶次跟踪 / 多记录怎么合并)',)      # ★ 末尾逗号: 这是"缺失项清单", 不是一整句话


def _norm(s: str) -> str:
    return re.sub(r'[^0-9a-z/]', '', str(s).lower())


def pos_token(option: str) -> str:
    """option 串 → 位置 token（见模块头第 2 条的实测规则）。"""
    s = str(option)
    m = re.search(r'pos\.?\s*([0-9]+[LR]?|[0-9.]+/[0-9]+)', s, re.I)
    if m:
        return _norm(m.group(1))
    m = re.search(r'([0-9.]+/[0-9]+)', s)
    if m:
        return _norm(m.group(1))
    m = re.search(r'([0-9]{3,4}\s*[A-Z]?\s*/\s*[A-Z0-9]+)', s)
    if m:
        return _norm(m.group(1))
    return _norm(s)


def theory_of(option: str, bearings: dict):
    tok = pos_token(option)
    cands = [k for k, v in bearings.items() if pos_token(v.get('pos', '')) == tok and tok]
    if not cands:
        return None, tok, []
    if len(cands) == 1:
        return cands[0], tok, cands
    no = _norm(option)
    for k in cands:
        mo = _norm(bearings[k].get('model', ''))
        if mo and (mo[:10] in no or mo[:8] in no):
            return k, tok, cands
    for k in cands:
        mk = str(bearings[k].get('model', '')).lower()
        for v in ('nsk', 'fag', 'ina'):
            if v in mk and v in no:
                return k, tok, cands
    return cands[0], tok, cands


def load_plan():
    if not PLAN.is_file():
        raise SystemExit(f'[X] 缺扫描配置 {P.rel(PLAN)}')
    if not THEORY.is_file():
        raise SystemExit(f'[X] 缺理论频率表 {P.rel(THEORY)}')
    plan = json.loads(PLAN.read_text(encoding='utf-8'))
    th = json.loads(THEORY.read_text(encoding='utf-8'))['bearings']
    out = []
    for c in plan['components']:
        k, tok, cands = theory_of(c['option'], th)
        e = th.get(k) or {}
        out.append(dict(**c, theory_key=k, pos_token=tok, n_cands=len(cands),
                        BPFI_hz=e.get('BPFI_Hz'), BPFO_hz=e.get('BPFO_Hz'), BSF_hz=e.get('BSF_Hz'),
                        shaft_Hz=e.get('shaft_Hz'), model=e.get('model', '')))
    return plan, out


def verify_plan(farm: str | None = None) -> int:
    """分母 + 配置 与随包样件逐值对拍（样件 = outputs/<场>/m5_cms_tcm/oem_frequency_scan.parquet）。"""
    farm = farm or P.farm()
    sample = P.m5(farm) / 'oem_frequency_scan.parquet'
    if not sample.is_file():
        print(f'[X] 没有样件可对拍: {P.rel(sample)}')
        return 2
    _plan, rows = load_plan()
    s = pd.read_parquet(sample)
    got = pd.DataFrame(rows)
    m = s[['component', 'sensor', 'shard', 'option', 'dx_hz', 'BPFI_hz', 'BPFO_hz', 'BSF_hz']].drop_duplicates()
    j = m.merge(got, on=['component', 'sensor', 'shard', 'option'], how='outer', suffixes=('_样', '_本'),
                indicator=True)
    both = j[j['_merge'] == 'both']
    ok = 0
    bad = []
    for _, r in both.iterrows():
        good = all(abs(float(r[f'{c}_样']) - float(r[f'{c}_本'])) < 1e-9 for c in ('BPFI_hz', 'BPFO_hz', 'BSF_hz')) \
            and abs(float(r['dx_hz_样']) - float(r['dx_hz_本'])) < 1e-9
        ok += good
        if not good:
            bad.append(r['component'])
    print(f'== oem_scan 分母+配置 逐值对拍 · 场站 {farm} ==')
    print(f'   样件组件 {len(m)} 个 · 本器 {len(got)} 个 · 同键 {len(both)} 个'
          f' · 仅样件 {int((j["_merge"] == "left_only").sum())} · 仅本器 {int((j["_merge"] == "right_only").sum())}')
    print(f'   BPFI_hz/BPFO_hz/BSF_hz + dx_hz 逐值一致: {ok}/{len(both)}' + (f'  不一致: {bad}' if bad else ''))
    for _, r in got.iterrows():
        print(f'     {r["component"]:24s} {r["sensor"]:22s} {r["shard"]:22s} dx={r["dx_hz"]:<8} '
              f'BPFI={r["BPFI_hz"]:<6} BPFO={r["BPFO_hz"]:<6} BSF={r["BSF_hz"]:<6} ← {str(r["theory_key"])[:34]}')
    rc = 0 if (ok == len(both) and len(both) == len(m) == len(got)) else 5
    print(f'   结论: {"分母映射成立（理论频率 + 扫描配置逐值一致）" if rc == 0 else "有不一致, 见上"}  rc={rc}')
    return rc


def status(farm: str | None = None) -> int:
    """机器可读状态（JSON）—— 供 scripts/chain_gap_check.py 汇总。"""
    _plan, rows = load_plan()
    st = dict(step='oem_scan', complete=bool(PEAK_PICKING_READY), components=len(rows),
              denominator='ok（理论频率 + 扫描配置已落地, --verify-plan 逐值对拍）',
              missing=list(MISSING) if not PEAK_PICKING_READY else [],
              request_sheet='docs/向振动线取料单_v0.1.md')
    print(json.dumps(st, ensure_ascii=False))
    return 0 if PEAK_PICKING_READY else 3


def skeleton_from_index(win: str = 'w0127'):
    r"""**从包内索引真推行骨架**（不是从样件抄）：窗 × 机组 × 部件/测点 × 分箱。

    骨架分两种成分, 必须分开说清楚（2026-09-17）：
      · **可从数据推**：某窗里"哪些机组有哪个测点的记录""某机组在该窗出现过哪些功率分箱" —— 索引里
        的 `turbine` / `sensor_name` / `condition_key` 直接给出；w0127 的索引在包内（`m5/tcm_index.parquet`）。
      · **只能从样件复刻**：**窗集合**与"哪些窗有哪几台"（w07xx/w08xx 的索引与谱库都不在本机）。
    本函数把前者算出来交给 `--skeleton-verify` 与样件对拍；后者在 `oem_scan_plan.json` 里记为"复刻自样件"。
    """
    ix = P.m5() / ('tcm_index.parquet' if win == 'w0127' else f'windows/{win}/index.parquet')
    if not ix.is_file():
        raise SystemExit(f'[X] 该窗索引不在位: {P.rel(ix)}（只有 w0127 与 w0316 的索引在包内）')
    d = pd.read_parquet(ix, columns=['turbine', 'sensor_name', 'condition_key'])
    out = {}
    for c in plan_components():
        sub = d[(d['sensor_name'] == c['sensor'])]
        turbine_sensor = set(map(tuple, sub[['turbine', 'sensor_name']].drop_duplicates().values))
        bins = (sub.groupby(['turbine', 'sensor_name', 'condition_key']).size().reset_index()[['turbine', 'sensor_name', 'condition_key']])
        out[c['component']] = dict(
            turbine_sensor=turbine_sensor,
            turbine_sensor_bin=set(map(tuple, bins.values)))
    return out


def plan_components():
    if not PLAN.is_file():
        raise SystemExit(f'[X] 缺扫描配置 {P.rel(PLAN)}')
    return json.loads(PLAN.read_text(encoding='utf-8'))['components']


def verify_skeleton(win: str = 'w0127', farm: str | None = None) -> int:
    """把"从索引推出来的骨架"与样件该窗的键对拍（w0127 是样件与索引都在的那一窗）。"""
    farm = farm or P.farm()
    idx_sk = skeleton_from_index(win)
    # 样件侧键
    out = {}
    b = pd.read_parquet(P.m5(farm) / 'bearing_freq_scan.parquet')
    b = b[b['window'] == win]
    bs = {}
    for comp, sens in (('GEN_bearing', ('Generator_DE', 'Generator_NDE')),
                       ('HS_bearing', ('Gear_HS_generator_side', 'Gear_HS_rotor_side'))):
        bs[comp] = set(map(tuple, b[b['sensor'].isin(sens)][['turbine', 'sensor']].drop_duplicates().values))
    out['bearing_freq_scan'] = bs
    c = pd.read_parquet(P.m5(farm) / 'cage_slip_scan.parquet')
    c = c[c['window'] == win]
    cs = {}
    for comp, sens in (('发电机DE', 'Generator_DE'), ('发电机NDE', 'Generator_NDE'), ('高速轴', 'Gear_HS_generator_side')):
        cs[comp] = set(map(tuple, c[c['sensor'] == sens][['turbine', 'sensor', 'bin']].drop_duplicates().values))
    out['cage_slip_scan'] = cs

    print(f'== 行骨架 从索引推 vs 样件 · 窗 {win} ==')
    rc = 0
    rows = []

    for comp, sens_set in (('GEN_bearing', ('Generator_DE', 'Generator_NDE')),
                           ('HS_bearing', ('Gear_HS_generator_side', 'Gear_HS_rotor_side'))):
        want = out['bearing_freq_scan'][comp]
        got = set()
        for c in plan_components():
            if c['sensor'] in sens_set:
                got |= {(t, s) for (t, s) in idx_sk[c['component']]['turbine_sensor']}
        rows.append((f'bearing_freq_scan/{comp}', len(want), len(got), len(want & got)))
    for comp, sens in (('发电机DE', 'Generator_DE'), ('发电机NDE', 'Generator_NDE'), ('高速轴', 'Gear_HS_generator_side')):
        want = out['cage_slip_scan'][comp]
        res = {}
        for c in plan_components():
            if c['sensor'] == sens:
                res[c['component']] = idx_sk[c['component']]['turbine_sensor_bin']
        got = set().union(*res.values()) if res else set()
        rows.append((f'cage_slip_scan/{comp}(机组×分箱)', len(want), len(got), len(want & got)))
    print('   %-38s %6s %6s %6s' % ('键集合', '样件', '索引推', '交集'))
    for name, w, g, i in rows:
        exact = (w == g == i)
        subset = name.startswith('cage_slip') and w <= g and i == w      # 样件是索引的子集
        rc = rc or (0 if (exact or subset) else 5)
        tag = '一致' if exact else ('样件是索引的子集(缺稳态筛选口径)' if subset else '有差异')
        print('   %-38s %6d %6d %6d  %s' % (name, w, g, i, tag))
    print('   结论: bearing/gear 的行骨架**可由包内索引直接推出**（74/74 逐键一致）；'
          'cage_slip 的键是索引的**子集**（样件 83 / 索引 282），少掉的那些由"稳态段筛选"决定 —— '
          '口径未知（n_rec 只有 2~7，与索引里的组记录数 47~1426 完全不是一回事）')
    return rc


def main() -> int:
    ap = argparse.ArgumentParser(description='六层链 oem_scan 步（逆向实现: 分母已就位/分子待口径）')
    ap.add_argument('--farm', default=None)
    ap.add_argument('--verify-plan', action='store_true', help='分母+配置 与随包样件逐值对拍（不写盘）')
    ap.add_argument('--plan', action='store_true', help='打印扫描计划')
    ap.add_argument('--status', action='store_true', help='机器可读状态（JSON）')
    ap.add_argument('--skeleton-verify', action='store_true', help='行骨架: 由索引推 vs 样件逐键对拍(w0127)')
    a = ap.parse_args()
    if a.status:
        return status(a.farm)
    if a.verify_plan:
        return verify_plan(a.farm)
    if a.skeleton_verify:
        return verify_skeleton('w0127', a.farm)
    _plan, rows = load_plan()
    if a.plan:
        print(f'== 扫描计划（{len(rows)} 个部件）==')
        for r in rows:
            print(f'  {r["component"]:24s} {r["sensor"]:22s} {r["shard"]:22s} dx={r["dx_hz"]:<8} '
                  f'BPFI={r["BPFI_hz"]:<6} BPFO={r["BPFO_hz"]:<6} BSF={r["BSF_hz"]:<6} shaft={r["shaft_Hz"]}')
        return 0
    if not PEAK_PICKING_READY:
        print('[X] 峰值拾取口径未到位 —— 本步**拒绝**产出半成品（只写分母会很危险: 那些件看着像扫描结果, 实则没有分子）')
        print(f'    缺的是: {MISSING[0]}')
        print('    要料: docs/向振动线取料单_v0.1.md 第 1 项（源码最佳, 一句话亦可）')
        print('    已就位: 扫描配置 + 理论频率 —— 用 --verify-plan 可逐值对拍（11/11 通过）')
        return 3
    raise SystemExit('[!] PEAK_PICKING_READY=True 但实现未完成 —— 请先补上分子侧实现再翻这个开关')


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