#!/usr/bin/env python3
"""cascade_sim.py — структурный симулятор тракта маски soothe2 (float-путь).

Цель (24mm3): воспроизвести scr(f)=дизайн-сигнал детектора; применённая
маска = exp(γ·scr), γ=1.760561 (точно), trk@688=exp(scr@628) бит-в-бит.

Структура канонической цепи (BLOCKMAP 24hh/24ii + 24mm2):
  шаг 9a: vec698 *= (1 − [54087c])            ; zero при дефолтах
  шаг 9b: vec6f8 += [54087c]·0.8              ; xmm10=0.8 @1824c3e28
  шаг 9c: bands_curve_i /= ... divide-ядро    ; dst=678i, A/B уточняются
  шаг 10: vec6f8 = bands_curve_i − ACC_i      ; dc40, ACC @таблицы 0x5407c8
  шаг 11: fma att/rel (тройки re/im/coef)     ; коэф @6c8/6e8
  шаг 12: COPY                                ; memcpy
  шаг 13: зеркало 9
  шаг 14: expf(bands_curve); bands_curve += (−1)   ; ПОРЯДОК исправлен 24mm2
  шаг 15: bands_curve *= track_i              ; th2000 array-mul
  шаг 16: bands_curve *= kWarp@[5406a8]
  шаг 17: expf ещё раз                        ; call-site 52b32c
  пост-17: exp-вариант(140a40) + pow?(140b00)
  FIR-секция: кривая-float(140b30→1803831c0) + sincos-twiddle(140aa0)

ЯДРА (структурная фаза — математически точные numpy-эквиваленты;
канонический C++ порт = инструкци-точная транскрипция, см. BLOCKMAP 24mm2):
"""
import numpy as np
import glob
import os

GAMMA = 1.760561          # 24mm3: показатель степени, rms фита 0 на чистых кадрах
N = 2049                  # число бинов полной сетки


# ---------------------------------------------------------------- ядра ----
def k_exp(x):
    """expf-ядро 180296c80. Структурная фаза: np.exp.
    Каноническая формула (для C++ порта, FMA-точно):
        n = fma(log2e_hi=1.4427f, x, 12582912.0f); k = n - MAGIC
        r = (x - 0.693146f*k) - 1.42861e-06f*k
        p = (((0.00829172f*r+0.0418735f)*r+0.166674f)*r+0.499994f)*r+1)*r+1
        out = bits((k<<23) + bits(p)); guard |x|>87.3365 -> slow path
    """
    return np.exp(x)


def k_div(a, b):
    """divide-ядро 1803a06a0: dst = B/A (~0.5 ulp, rcp+таблицы+полином).
    Структурная фаза: точное деление."""
    return b / a


# ------------------------------------------------------------ данные -----
def load_tract(path):
    """tract_*.txt: k am res lvl_raw band_level prewarp w"""
    t = np.loadtxt(path)
    return {'am': t[:, 1], 'res': t[:, 2], 'lvl': t[:, 3]}


def load_frame(npz):
    """Слоты кадра rendersnap2 → dict[int, np.ndarray]."""
    d = np.load(npz)
    out = {}
    for k in d.keys():
        if k.startswith('0x'):
            out[int(k[2:], 16)] = d[k]
    return out, d['t_snap']


def pick_clean_frame(ds_dir, min_bins=8):
    """Отбор чистых стационарных кадров по фазам (24mm3):
       возвращает лучший на фазе γ* (~1.7606, маска применена)
       и лучший на фазе γ=1 (степень ещё не применена)."""
    classes = {'gamma': None, 'identity': None}
    for f in sorted(glob.glob(os.path.join(ds_dir, 'ph*.npz'))):
        try:
            S, ts = load_frame(f)
        except Exception:
            continue
        if not all(x in S for x in (0x540628, 0x540688, 0x540678)):
            continue
        s = S[0x540628][:1025].astype(np.float64)
        t = S[0x540688][:1025].astype(np.float64)
        c = S[0x540678][:1025].astype(np.float64)
        ok = (t > 1e-30) & (c > 1e-30) & np.isfinite(s)
        if ok.sum() < 50:
            continue
        lt = np.log(t[ok])
        lc = np.log(c[ok])
        sel = np.abs(lt) > 0.05
        if sel.sum() < min_bins:
            continue
        g = float(np.sum(lt[sel] * lc[sel]) / np.sum(lt[sel] ** 2))
        rms = float(np.sqrt(np.mean((lc[sel] - g * lt[sel]) ** 2)))
        depth = float(-lc.min())
        key = 'gamma' if abs(g - GAMMA) < 0.01 else \
              ('identity' if abs(g - 1.0) < 1e-4 else None)
        if key is None or rms > 1e-4:
            continue
        cand = (depth, f, s, t, c, g, rms)
        if classes[key] is None or depth > classes[key][0]:
            classes[key] = cand
    return classes


# ------------------------------------------------------- валидация -------
def validate_scr(sim_scr, cap_scr, tol_db=0.05):
    """rms в дБ между симулированным и захваченным scr."""
    m = np.abs(cap_scr) > 0.02
    err = (sim_scr[m] - cap_scr[m]) * (20 / np.log(10))
    return float(np.sqrt(np.mean(err ** 2))), int(m.sum())


def main():
    import sys
    ds = sys.argv[1] if len(sys.argv) > 1 else '/tmp/opencode/sc_multi6'
    tract = sys.argv[2] if len(sys.argv) > 2 else '/tmp/opencode/tract_multi6.txt'

    classes = pick_clean_frame(ds)
    ph_g = classes['gamma']
    ph_i = classes['identity']
    if not ph_g and not ph_i:
        print('нет чистых кадров в', ds)
        return
    for lbl, best in (('γ-фаза', ph_g), ('identity', ph_i)):
        if not best:
            continue
        _, f, scr, trk, cur, gamma_fit, grms = best
        n = len(scr)
        cut_meas = -20 / np.log(10) * np.log(np.maximum(cur, 1e-30))
        g_use = gamma_fit
        cut_sim = g_use * (-scr) * 20 / np.log(10)
        e = cut_sim - cut_meas
        sel = np.abs(cut_meas) > 0.1
        rms_db = float(np.sqrt(np.mean(e[sel] ** 2))) if sel.any() else 0.0
        print(f'{lbl}: {os.path.basename(f)} γ={gamma_fit:.6f} (rms {grms:.1e}) '
              f'закон: rms={rms_db:.4f} дБ / {int(sel.sum())} бинов')


if __name__ == '__main__':
    main()