1"""Decompose the measured gap between LPC+ANS and the entropy rate into named
2terms, and check each against what the theory predicts.
4 R_ANS - Hbar
5 = [R_ANS - H0(resid)] coder overhead: rANS vs a perfect memoryless
6 coder on its own residual stream
7 + [H0(resid) - R_condG] residual-model loss: the pooled integer
8 histogram vs the phase-conditioned law.
9 Theory says this is L(s) = M(s) - G(s).
10 + [R_condG - G(s_*)] prediction suboptimality + parametric
11 mismatch of the single-scale Gaussian
12 + [G(s_*) - Hbar] error of the analytic rate itself
14`s_fit` is the residual scale the conditional-Gaussian model actually fits,
15which is the honest argument to L(.) — `s_*` is its prediction from the
16spectrum, and the two differing is itself informative.
17"""
18import os
19import sys
21import numpy as np
22from scipy.linalg import solve_toeplitz
23from scipy.signal import lfilter
24from scipy.special import ndtr
26sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
27import codec_suite as cc # noqa: E402
28from fitmodel import Fit # noqa: E402
29from phase_loss import phase_entropies # noqa: E402
30from sweep import noise_units, STEPS # noqa: E402
32CACHE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
35def cond_gauss_detail(z, order):
36 """Conditional-Gaussian rate plus the residual scale it fits."""
37 zf = np.asarray(z, dtype=np.float64)
38 n = zf.size
39 r = cc.autocorr(zf, order) / n
40 a = solve_toeplitz(r[:order], r[1:order + 1])
41 pred = lfilter(np.concatenate(([0.0], a)), [1.0], zf)
42 zt, mu = zf[order:], pred[order:]
43 s0 = np.sqrt(max((zt - mu).var() - 1.0 / 12.0, 1e-6))
44 best_bits, best_s = np.inf, s0
45 for s in s0 * np.linspace(0.7, 1.3, 25):
46 p = ndtr((zt + 0.5 - mu) / s) - ndtr((zt - 0.5 - mu) / s)
47 bits = float(-np.log2(np.maximum(p, 1e-12)).mean())
48 if bits < best_bits:
49 best_bits, best_s = bits, float(s)
50 return best_bits + 8 * (4 * order + 14) / n, best_s
53def main(n=200_000, order=32, trace="001290"):
54 path = [p for p in sorted(os.listdir(CACHE))
55 if trace in p and p.endswith(".raw.npy")][0]
56 y = noise_units(np.load(os.path.join(CACHE, path)))[:n]
57 print(f"{path} n = {n} order = {order}\n")
58 print(f" {'v':>5}{'s_*':>7}{'s_fit':>7} | "
59 f"{'R_ANS':>7}{'H0':>7}{'R_cG':>7}{'G(s*)':>7} | "
60 f"{'coder':>7}{'resid':>7}{'pred':>7} | {'L(s_fit)':>9}{'L(s_*)':>8}")
61 print(" " + "-" * 92)
62 for v in STEPS:
63 z = np.round(y / v).astype(np.int16)
64 if np.unique(z).size < 3:
65 continue
66 fit = Fit(z, nfft=4096)
67 nb, h0 = cc.lpc_ans_bytes(z, order)
68 r_ans = 8.0 * nb / z.size
69 r_cg, s_fit = cond_gauss_detail(z, order)
70 g = fit.predicted_rate
71 coder, resid, pred = r_ans - h0, h0 - r_cg, r_cg - g
72 print(f" {v:5.2f}{fit.s_star:7.3f}{s_fit:7.3f} | "
73 f"{r_ans:7.3f}{h0:7.3f}{r_cg:7.3f}{g:7.3f} | "
74 f"{coder:7.3f}{resid:7.3f}{pred:7.3f} | "
75 f"{phase_entropies(s_fit)[2]:9.3f}{phase_entropies(fit.s_star)[2]:8.3f}")
78if __name__ == "__main__":
79 main(trace=sys.argv[1] if len(sys.argv) > 1 else "001290")