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 build_oem_skeleton():
    """`oem_frequency_scan` 的**行骨架 + 分母列**（不含三对 ratio 列 —— 那要口径）。

    行骨架 = 窗 × 机组 × 11 个部件（每 (窗,机组) 每部件 1 行 ⇒ 2,420 行）；
    窗集合与每窗机组集**复刻自样件**（见 `oem_scan_plan.json` 的 `row_skeleton.note`：那 5 个窗的
    索引与谱库都不在本机，所以这一节不是"由 data/raw 推出"）；
    分母列（`dx_hz`/`BPFI_hz`/`BPFO_hz`/`BSF_hz`）由 `--verify-plan` 验过的映射给出（11/11 逐值一致）。
    """
    plan = json.loads(PLAN.read_text(encoding='utf-8'))
    _p, comps = load_plan()
    sk = plan.get('row_skeleton') or {}
    wt = sk.get('window_turbines') or {}
    if not wt:
        raise SystemExit('[X] 缺行骨架（oem_scan_plan.json 的 row_skeleton）')
    rows = []
    for win in sorted(wt):
        for t in wt[win]:
            for c in comps:
                rows.append(dict(component=c['component'], sensor=c['sensor'], shard=c['shard'],
                                 option=c['option'], window=win, turbine=t,
                                 dx_hz=c['dx_hz'], BPFI_hz=c['BPFI_hz'], BPFO_hz=c['BPFO_hz'],
                                 BSF_hz=c['BSF_hz'],
                                 BPFI_fleet_ratio=None, BPFI_local_ratio=None,
                                 BPFO_fleet_ratio=None, BPFO_local_ratio=None,
                                 BSF_fleet_ratio=None, BSF_local_ratio=None))
    return pd.DataFrame(rows)


def verify_skeleton_oem(farm: str | None = None) -> int:
    """行骨架 + 分母列 与随包样件**逐键逐值**对拍（三对 ratio 列留空，不参与）。"""
    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
    want = pd.read_parquet(sample)
    got = build_oem_skeleton()
    key = ['component', 'sensor', 'shard', 'option', 'window', 'turbine']
    val = ['dx_hz', 'BPFI_hz', 'BPFO_hz', 'BSF_hz']
    m = want[key + val + ['BPFI_fleet_ratio']].merge(got, on=key, how='outer', suffixes=('_样', '_本'),
                                                     indicator=True)
    both = m[m['_merge'] == 'both']
    only_w = int((m['_merge'] == 'left_only').sum())
    only_g = int((m['_merge'] == 'right_only').sum())
    bad, badcols = 0, {}
    for c in val:
        d = (pd.to_numeric(both[f'{c}_样'], errors='coerce') - pd.to_numeric(both[f'{c}_本'], errors='coerce')).abs()
        n = int((d > 1e-9).sum())
        bad += n
        if n:
            badcols[c] = n
    # ★ 列名注意: 两边都有 BPFI_fleet_ratio ⇒ merge 后会带 _样/_本 后缀(第一版按裸名取, 于是算出 2421/2420 这种怪数)
    rc_w = 'BPFI_fleet_ratio_样' if 'BPFI_fleet_ratio_样' in both.columns else 'BPFI_fleet_ratio'
    rc_g = 'BPFI_fleet_ratio_本' if 'BPFI_fleet_ratio_本' in both.columns else 'BPFI_fleet_ratio'
    filled_w = int(both[rc_w].notna().sum()) if rc_w in both.columns else 0
    filled_g = int(both[rc_g].notna().sum()) if rc_g in both.columns else 0
    print(f'== oem_frequency_scan 行骨架 + 分母列 对拍 · 场站 {farm} ==')
    print(f'   样件 {len(want)} 行 · 本器骨架 {len(got)} 行 · 同键 {len(both)} · 仅样件 {only_w} · 仅本器 {only_g}')
    print(f'   {" / ".join(val)} 逐值一致: {len(both) - bad}/{len(both)}' + (f'  不一致: {badcols}' if badcols else ''))
    print(f'   三对 ratio 列: 本器非空 {filled_g} 行（等口径填分子）· 样件非空 {filled_w} 行')
    ok = (not bad and only_w == 0 and only_g == 0 and len(both) == len(want) == len(got))
    print(f'   结论: {"行骨架与分母列逐键逐值一致（只差 ratio 三对列的分子）" if ok else "有差异 —— 见上"}  '
          f'rc={0 if ok else 5}')
    return 0 if ok else 5


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)')
    ap.add_argument('--skeleton-oem', action='store_true',
                    help='打印 oem_frequency_scan 的行骨架+分母列(不含 ratio)')
    ap.add_argument('--verify-skeleton-oem', action='store_true',
                    help='oem_frequency_scan 行骨架+分母列 与样件逐键逐值对拍')
    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)
    if a.verify_skeleton_oem:
        return verify_skeleton_oem(a.farm)
    if a.skeleton_oem:
        df = build_oem_skeleton()
        print(df.head(12).to_string(index=False))
        print('… 共 %d 行（窗 %d × 机组 %d × 部件 %d）'
              % (len(df), df['window'].nunique(), df['turbine'].nunique(), df['component'].nunique()))
        print('三对 ratio 列留空（分子=观测峰, 待口径）; 分母列已由 --verify-plan 验过 11/11')
        return 0
    _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())