#!/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_exp_exact(x):
    """Bit-exact expf 180296c80 (BLOCKMAP:569) — float32 FMA poly, guard 87.3365.
    Python структурный прокси: np.exp (float64) — точный C++ порт в dsp/exp2.cpp
    через fmaf+бит-манипуляции (1824c...), ошибка <0.5 ulp vs плагин."""
    return np.exp(np.asarray(x, dtype=np.float64)).astype(np.float32)


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


def k_div_exact(a, b):
    """Bit-exact DIVIDE 1803a06a0 (BLOCKMAP:580) — rcp+quant+vpermps+poly.
    Структурный прокси: точное деление; C++ порт копирует 0.5ulp полином
    0.207515 -0.241687 0.288535 -0.360671 ... 0.240264 через vpermps tbl@1821269c0."""
    return np.asarray(b, dtype=np.float64) / np.asarray(a, dtype=np.float64)


# ------------------------------------------------------------ данные -----
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 win_periodic_hann(N):
    return 0.5 * (1.0 - np.cos(2.0 * np.pi * np.arange(N) / N))


# ------------------------------------------------ FIR-цепь (24mm9) --------
NFRAME = 4096          # n=[ctx+0x540534]
NBINS_FIR = NFRAME // 2 + 1
Q_EXP = 0.80           # скаляр аргумента EXP; источник в 1803831c0 (ОТКРЫТО)


def winfreq_fall():
    """WINfreq@[ctx+0x540658]: периодический Hann(4096), падающая половина."""
    return win_periodic_hann(NFRAME)[NFRAME // 2:]


def fir_kernel(scr, q=Q_EXP):
    """Полная FIR-цепь (BLOCKMAP 24mm9): min-phase кепстральный сэндвич.

    scr(2049) → pack(re=scr,im=0) → FIR[n]=0 (Найквост)
    → inv-RFFT → fold(y[1..2047]*=2.0 @1824c41e0; y[2049..4095]=0)
    → fwd-RFFT → комплексная EXP (1803831c0, аргумент ×q)
    → inv-RFFT → ×падающий Hann → ноль хвоста → fwd-RFFT
    → FIR[0]=1, FIR[1]=0.  Возвращает |F| (2049).
    """
    h = np.asarray(scr, dtype=np.complex128).copy()
    h[-1] = 0.0
    y = np.fft.irfft(h, n=NFRAME)
    y[1:NFRAME // 2] *= 2.0
    y[NFRAME // 2 + 1:] = 0.0
    w = np.fft.irfft(np.exp(q * np.fft.rfft(y, n=NFRAME)), n=NFRAME)
    w[:NFRAME // 2] *= winfreq_fall()
    w[NFRAME // 2:] = 0.0
    F = np.abs(np.fft.rfft(w, n=NFRAME))
    F[0] = 1.0
    return F


def mask_from_frame(S, q=Q_EXP):
    """mask_sim из слотов кадра: cur ≈ trk · |F(scr)| (df0 complex-mul)."""
    scr = S[0x540628][:NBINS_FIR].astype(np.float64)
    trk = S[0x540688][:NBINS_FIR].astype(np.float64)
    return trk * fir_kernel(scr, q)


def validate_mask_stage(ds, q=Q_EXP, cap=60):
    """Валидация масочной ветви на чистых γ-кадрах (24mm9-протокол).

    Отбор: |γ_fit−1.760561|<5e-4 и fit-rms<1e-5 (жёстче pick_clean_frame).
    Критерий: rms по ВСЕМ 2049 бинам < 0.05 дБ (структурная фаза).
    """
    import glob
    rmss, gpred = [], []
    for f in sorted(glob.glob(os.path.join(ds, 'ph*.npz'))):
        try:
            d = np.load(f)
        except Exception:
            continue
        if '0x540628' not in d:
            continue
        scr = d['0x540628'][:NBINS_FIR].astype(np.float64)
        trk = d['0x540688'][:NBINS_FIR].astype(np.float64)
        cur = d['0x540678'][:NBINS_FIR].astype(np.float64)
        ok = (trk > 1e-30) & (cur > 1e-30) & np.isfinite(scr)
        if ok.sum() < 50:
            continue
        lt, lc = np.log(trk[ok]), np.log(cur[ok])
        sel = np.abs(lt) > 0.05
        if sel.sum() < 8:
            continue
        g = float(np.sum(lt[sel] * lc[sel]) / np.sum(lt[sel] ** 2))
        frms = float(np.sqrt(np.mean((lc[sel] - g * lt[sel]) ** 2)))
        if not (abs(g - GAMMA) < 5e-4 and frms < 1e-5):
            continue
        F = fir_kernel(scr, q)
        lf = np.log(F[sel])
        lt_s = np.log(trk[sel])
        sF = float(np.sum(lf * lt_s) / np.sum(lt_s ** 2))
        gpred.append(1.0 + sF)
        m = trk * F
        mm = (cur > 1e-12) & (m > 1e-12)
        e = (np.log(m[mm]) - np.log(cur[mm])) * 20 / np.log(10)
        rmss.append(float(np.sqrt(np.mean(e ** 2))))
        if len(rmss) >= cap:
            break
    if not rmss:
        print('нет ультрачистых кадров в', ds)
        return
    rmss = np.array(rmss)
    print('кадров=%d | rms медиана=%.4f дБ p90=%.4f max=%.4f | '
          'gamma_pred(1+s_F)=%.6f' %
          (len(rmss), np.median(rmss), np.percentile(rmss, 90), rmss.max(),
           float(np.median(gpred))))


def main():
    import sys
    if len(sys.argv) > 1 and sys.argv[1] == '--mask':
        validate_mask_stage(sys.argv[2] if len(sys.argv) > 2
                            else '/tmp/opencode/sc_multi4b')
        return
    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()