c4b8c52Check the model against real ephys, and measure FLACJeremy Magland 1# Does the quantized-filtered-Gaussian model predict real ephys compressibility?
3Test: fit `x ~ iid N(0,1) → h * x → round` to a real ephys trace using nothing
4but its power spectrum, synthesize a surrogate of the same length from the fit,
5and run the same 19 codecs on both. If the model is a good stand-in, every
6codec should land at the same bits/sample on the surrogate as on the real
7trace.
9Data: the three single-channel traces benchcompress benchmarks — DANDI 000876
10ch45, 000409 ch101, 001290 ch0 — 200k samples each, raw and in the `-filtered`
11variant (bandpass 300–6000 Hz, normalized to MAD noise = 1, requantized at step
12v = 0.25). Scripts: `fetch.py`, `fitmodel.py`, `codec_suite.py`, `compare.py`,
13`sweep.py`. Logs: `full_compare.log`, `sweep.log`.
15## Verdict
17**The model is a good surrogate for compressibility, and a bad surrogate for
18the signal.** The surrogate has kurtosis ~3.0 against 4.5–19.2 for the real
19traces — it has no spikes at all — yet it reproduces the measured rate of every
20prediction-based codec to within **0.05 bits/sample**, and of the generic byte
21compressors to within ~0.05 bits/sample on the filtered variants. Compression
22of these recordings is set by the noise spectrum, not by the spikes.
24**But the interesting regime is not where ephys currently sits.** At the
25deployed quantization (noise = 4 steps, `s_* ≈ 1`) the fractional-phase loss is
26only 0.06–0.08 bits/sample (~3%), and LPC+ANS is within ~6% of the analytic
27entropy rate of the fitted model. The loss becomes large only when the
28quantizer is coarsened to ~0.5 steps of noise — and the model predicts exactly
29where that happens.
31## 1. Model vs real, at the benchcompress operating points
33`-filtered` variants (bandpass, noise = 4 quantization steps). bits/sample.
35| method | 000409 real / model | 000876 real / model | 001290 real / model |
36|---|---|---|---|
37| order-0 H(z) | 3.871 / 4.070 | 4.236 / 4.298 | 4.129 / 4.150 |
38| zstd+delta | 2.701 / 2.825 | 3.065 / 3.105 | 3.348 / 3.303 |
39| lzma+delta | 2.507 / 2.601 | 2.823 / 2.854 | 3.032 / 3.056 |
40| bz2 | 2.519 / 2.580 | 2.752 / 2.739 | 2.919 / 2.935 |
41| brotli+delta | 2.692 / 2.806 | 3.019 / 3.043 | 3.223 / 3.238 |
42| **LPC(32)+ANS** | **2.083 / 2.128** | **2.261 / 2.270** | **2.369 / 2.391** |
43| LPC(32) resid H0 | 2.053 / 2.094 | 2.215 / 2.232 | 2.334 / 2.347 |
44| **cond-Gauss(32)** | **2.029 / 2.039** | **2.173 / 2.190** | **2.301 / 2.311** |
45| *G(s\*) analytic rate* | *1.977* | *2.134* | *2.237* |
46| *s\** | *0.908* | *1.022* | *1.104* |
47| *kurtosis real / model* | *17.7 / 3.0* | *6.0 / 3.0* | *5.9 / 3.0* |
49Raw (unfiltered) traces agree nearly as well for the prediction-based methods
50(|Δ| ≤ 0.05 bits) and worse for the weak general-purpose coders (lz4 off by
51up to 0.5 bits) — those are the ones most sensitive to the heavy tails the
52Gaussian model does not have.
54One visible fit artifact: on 000876 raw the surrogate std is 26.1 against 47.9
55real. Welch detrends each segment, so the fit discards drift below ~7 Hz. That
56content is almost perfectly predictable, so it moves `order-0 H(z)` by 0.77
57bits and `LPC(32)+ANS` by 0.01 bits.
59## 2. Walking the s\* axis by coarsening the quantizer
61Same recordings, bandpassed and normalized to noise = 1.0, requantized at step
62`v`. `gap` = LPC+ANS − cond-Gauss, i.e. what integer-residual coding loses by
63throwing away the fractional part of the prediction. 200k samples.
64(001290 ch0; the other two traces agree to ~0.05 bits — see `sweep.log`.)
66| v | noise/step | s\* | G(s\*) | real LPC+ANS | real condG | real gap | model gap | model−real |
67|---|---|---|---|---|---|---|---|---|
68| 0.125 | 8.0 | 1.658 | 2.798 | 2.930 | 2.859 | 0.072 | 0.054 | +0.004 |
69| 0.25 | 4.0 | 1.099 | 2.232 | 2.360 | 2.295 | 0.065 | 0.082 | +0.029 |
70| 0.5 | 2.0 | 0.741 | 1.717 | 1.872 | 1.771 | 0.101 | 0.101 | +0.012 |
71| 1.0 | 1.0 | 0.500 | 1.255 | 1.489 | 1.296 | **0.193** | 0.185 | +0.004 |
72| 2.0 | 0.5 | 0.325 | 0.842 | 1.260 | 0.836 | **0.424** | 0.405 | −0.013 |
73| 4.0 | 0.25 | 0.103 | 0.267 | 0.951 | 1.938 | — | — | −0.943 |
74| 8.0 | 0.12 | 0.001 | 0.001 | 0.102 | 0.069 | — | — | −0.099 |
76Three things to read off this:
781. **The predicted mechanism is real and it is large.** The fractional-phase
79 gap grows monotonically as `s_*` falls — 0.07 → 0.10 → 0.19 → 0.42
80 bits/sample — exactly the behaviour the synthetic study predicted, now
81 measured on real recordings. At `v = 2` the conditional-Gaussian coder is
82 **34% smaller** than LPC+ANS on real data.
832. **The model predicts the gap, not just the rate.** Real and surrogate gaps
84 agree to ~0.02 bits at every step down to `v = 2`. So the model can be used
85 to answer "what would I gain?" without running the codec.
863. **The model fails when the quantizer step exceeds the noise** (`v ≥ 4`,
87 noise ≤ 0.25 steps). There the dither is gone: the real trace becomes a
88 sparse spike train on a bed of zeros while the Gaussian surrogate collapses
89 to all zeros (0.008 bits/sample against 0.95 real). The conditional-Gaussian
90 coder also inverts there and becomes *worse* than LPC+ANS on real data — a
91 single-scale Gaussian is the wrong conditional law for a sparse spiky
92 signal. **`s_* ≳ 0.3`, equivalently noise ≳ 0.5 quantization steps, is the
93 model's domain of validity.**
95## 1b. LPC+ANS beats FLAC by ~21% at the deployed settings
97FLAC is the mature instance of the same architecture — LPC prediction, integer
98residual, memoryless coding — differing only in that it Rice-codes the residual
99where ANS codes its empirical histogram. That one substitution is worth a fifth
100of the file:
102| | 000409 | 000876 | 001290 |
103|---|---|---|---|
104| FLAC (LPC + Rice) | 2.680 | 2.865 | 3.043 |
105| **LPC(32) + ANS** | **2.083** | **2.261** | **2.369** |
106| LPC(32) residual H0 — the memoryless ceiling | 2.053 | 2.215 | 2.333 |
107| conditional-Gaussian | 2.029 | 2.173 | 2.301 |
108| G(s\*) analytic entropy rate | 1.977 | 2.134 | 2.237 |
110Read as distance above the residual-entropy ceiling: **FLAC sits 29–31% above
111it, LPC+ANS 1.5–2.1%.** The loss is entirely in the residual coder, not the
112predictor. On raw traces (`s_*` = 3.8–6.8) the advantage shrinks to 1–3%: Rice
113coding is near-optimal for wide residuals and mismatched for narrow ones, which
114is the same $s_*$ story from a different direction.
116Caveat: libsndfile does not expose the FLAC compression level, so this is its
117default (level 5). `flac -8` raises the max LPC order and improves residual
118partitioning — worth a few percent, not twenty. **The paper must re-measure
119with the reference `flac -8` binary and with WavPack** (whose prebuilt wheels
120do not match this machine's glibc).
122## 2b. The model does not need an explicit filtering step
124`h` in the fit is the minimum-phase spectral factor of whatever spectrum the
125trace has — not a filter anyone applied. Any stationary Gaussian process with
126spectral density `S` can be written `h * x`. Two checks that this is not just a
127formal remark:
129- **The raw traces were never preprocessed by us.** Only the acquisition chain
130 shaped them, and the fit reproduces them as well as the bandpassed versions:
131 LPC+ANS real/model 4.855/4.902, 4.522/4.512, 4.074/4.117; cond-Gauss agrees
132 to ≤ 0.015 bits on all three. Kurtosis up to 19.2.
133- **The whole `s_*` sweep repeated with no bandpass at all** (`sweep_nofilter.log`,
134 `--nofilter`, same normalization and same steps, only the bandpass removed).
135 Model−real stays ≤ 0.06 bits from `s_* = 6.8` down to `s_* ≈ 0.32`, and the
136 fractional-phase gap grows the same way. 000409 ch101:
138 | noise/step | s\* | real gap | model gap | model−real |
139 |---|---|---|---|---|
140 | 8.0 | 4.350 | 0.048 | 0.091 | +0.056 |
141 | 4.0 | 2.190 | 0.014 | 0.061 | +0.058 |
142 | 2.0 | 1.112 | 0.058 | 0.069 | +0.017 |
143 | 1.0 | 0.593 | 0.201 | 0.198 | −0.003 |
144 | 0.5 | 0.319 | 0.431 | 0.477 | −0.024 |
146So the preconditions are **stationarity over the block**, **second-order
147statistics sufficing for the rate** (validated at kurtosis 19), and **noise
148≳ 0.5 quantization steps**. Explicit filtering is not among them.
150## 2c. The gap decomposes into named terms, and the structural one has a
151closed form
153Write the distance from a measured codec rate to the analytic entropy rate as
155```
156R_ANS - G(s*) = [R_ANS - H0] + [H0 - R_condG] + [R_condG - G(s*)]
157 coder residual-model prediction
158```
160- **coder** — rANS overhead against a perfect memoryless coder on its own
161 residual stream (symbol table + arithmetic loss).
162- **residual-model** — the pooled integer-residual histogram against the
163 phase-conditioned law. *This is the structural term, and theory predicts it
164 from `s` alone*: with `M(s)` the entropy of the phase-mixture and `G(s)` the
165 phase-averaged conditional entropy, the loss is `L(s) = M(s) - G(s)`, a
166 universal curve (`phase_loss.py`). It peaks near `s ≈ 0.2` at ~0.27
167 bits/sample and decays in both directions.
168- **prediction** — linear prediction is not optimal from the quantized past,
169 plus the single-scale Gaussian's parametric mismatch.
171Measured against predicted, three recordings, 200k samples (`decompose.log`):
173| s\* | coder | residual-model measured | **L(s) predicted** | prediction |
174|---|---|---|---|---|
175| 1.36–1.66 | 0.052–0.059 | −0.009 – 0.016 | 0.019–0.029 | 0.032–0.061 |
176| 0.92–1.10 | 0.029–0.046 | 0.025–0.042 | 0.042–0.058 | 0.039–0.063 |
177| 0.63–0.74 | 0.023–0.028 | 0.078–0.085 | 0.084–0.109 | 0.033–0.054 |
178| 0.43–0.50 | 0.022–0.034 | 0.159–0.272 | 0.153–0.181 | 0.015–0.042 |
179| 0.28–0.33 | 0.016–0.028 | 0.396–0.492 | 0.235–0.271 | −0.071 – −0.006 |
181Read honestly: **the closed form gets the structural term right to within about
182a factor of 1.5 over `s_* ∈ [0.4, 1.7]`, and underestimates it by ~1.7× at
183`s_* ≈ 0.3`** — where the real residual is heavier-tailed than the single
184Gaussian the theory assumes. The coder term is a stable 1–2% of the rate, as
185claimed. The prediction term is small and positive until `s_*` drops below
186~0.35, where it goes negative because `G(s_*)` itself starts to break down.
188## 3. Consequences for the paper
190- §8's working hypothesis was that deployed ephys sits in the low-`s_*` region
191 where the fractional-phase loss is largest. **It does not.** At the standard
192 filtered-and-quantized settings `s_* ≈ 1`, and LPC+ANS is within a few
193 percent of the limit. The honest headline for real data is: *for the way
194 ephys is stored today, standard prediction + entropy coding is close to
195 optimal, and the model says so without any experiment.*
196- The gain lives one step coarser. Since coarsening the quantizer is exactly
197 the near-lossless knob practitioners already reach for, the useful statement
198 is a joint one: at each step size, here is the entropy rate, here is what
199 standard coding gets, and here is what conditional-Gaussian coding gets.
200- Spikes cost almost nothing in bits. Worth stating plainly — it is the reason
201 a Gaussian surrogate works at all, and it is counterintuitive.
202- The Welch-detrending artifact means the fit should either keep the very low
203 frequencies or the comparison should be stated on high-passed data. Minor,
204 but it affects `order-0 H(z)` a lot.
206## Open
208- Only linear-prediction methods and generic byte compressors so far. FLAC and
209 WavPack are not installed here; they should land near LPC+ANS.
210- `s_*` here is computed from the *observed* spectrum via Szegő + 1/12. It has
211 not been checked against the Monte-Carlo estimator on the fitted kernel.
212- Single channel. Cross-channel redundancy untested.