1"""Does the paper's model explain real ephys compressibility?
3For each cached trace: fit x ~ N(0,1) -> h * x -> round to its spectrum,
4synthesize a surrogate of the same length from the fit, and run the same codec
5suite on both. If the model is a good stand-in, every codec should land at the
6same bits/sample on the surrogate as on the real trace.
8Usage: python compare.py [--n N] [--nfft NFFT] [--taps T] [--order O] [--fast]
9"""
10import argparse
11import glob
12import os
13import sys
14import time
16import numpy as np
18sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
19import codec_suite as cc # noqa: E402
20from fitmodel import Fit # noqa: E402
22CACHE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
25def report(name, z_real, args):
26 z_real = np.asarray(z_real, dtype=np.int16)[: args.n]
27 z_real = (z_real - int(np.median(z_real))).astype(np.int16)
29 fit = Fit(z_real, nfft=args.nfft, n_taps=args.taps)
30 z_syn = fit.synthesize(seed=0)
32 print(f"\n{'=' * 78}\n{name} n = {z_real.size}\n{'=' * 78}")
33 print(f" real: std {z_real.std():8.3f} "
34 f"kurtosis {float(((z_real - z_real.mean()) ** 4).mean() / z_real.var() ** 2):6.2f}")
35 print(f" surrogate: std {z_syn.std():8.3f} "
36 f"kurtosis {float(((z_syn - z_syn.mean()) ** 4).mean() / z_syn.var() ** 2):6.2f}")
37 print(f" fit: taps {fit.kernel.size} spectrum RMS error "
38 f"{fit.kernel_error_db():.2f} dB")
39 print(f" s_* = {fit.s_star:.4f} quantization steps -> "
40 f"predicted entropy rate G(s_*) = {fit.predicted_rate:.4f} bits/sample "
41 f"({16 / max(fit.predicted_rate, 1e-9):.1f}x)")
43 t0 = time.time()
44 a = cc.measure(z_real, lpc_order=args.order, generic=not args.fast)
45 b = cc.measure(z_syn, lpc_order=args.order, generic=not args.fast)
46 print(f" [{time.time() - t0:.1f}s]")
48 print(f"\n {'method':<22}{'real':>9}{'model':>9}{'diff':>8}"
49 f"{'real x':>9}{'model x':>9}")
50 print(f" {'-' * 66}")
51 for k in a:
52 d = b[k] - a[k]
53 print(f" {k:<22}{a[k]:9.4f}{b[k]:9.4f}{d:+8.4f}"
54 f"{16 / a[k]:9.2f}{16 / b[k]:9.2f}")
55 return {"name": name, "fit": fit, "real": a, "model": b}
58def main():
59 p = argparse.ArgumentParser()
60 p.add_argument("--n", type=int, default=100_000)
61 p.add_argument("--nfft", type=int, default=4096)
62 p.add_argument("--taps", type=int, default=None)
63 p.add_argument("--order", type=int, default=32)
64 p.add_argument("--fast", action="store_true",
65 help="skip the generic byte compressors")
66 p.add_argument("--only", type=str, default=None)
67 args = p.parse_args()
69 paths = sorted(glob.glob(os.path.join(CACHE, "*.npy")))
70 if args.only:
71 paths = [q for q in paths if args.only in os.path.basename(q)]
72 if not paths:
73 sys.exit("no cached traces; run fetch.py first")
75 for path in paths:
76 name = os.path.basename(path)[: -len(".npy")]
77 report(name, np.load(path), args)
80if __name__ == "__main__":
81 main()