Check the model against real ephys, and measure FLAC
Fit the generative model to three DANDI traces from their power spectrum
alone (Welch PSD, minus the 1/12 dither floor, minimum-phase spectral
factor) and run the same codecs on the recording and on a surrogate drawn
from the fit.
The surrogate has no spikes at all — kurtosis 3.0 against 4.5–19.2 — yet
reproduces every prediction-based codec to within 0.05 bits/sample, on raw
traces as well as bandpassed ones. Compressibility of these recordings is
set by the noise spectrum, not by the events.
Also measured: FLAC sits 29–31% above the residual-entropy ceiling at the
deployed quantization where LPC+ANS sits 1.5–2.1% above it, and the
fractional-phase loss follows a universal curve L(s) = M(s) - G(s) that
depends on nothing but s.
Numbers and caveats in exploration/ephys/RESULTS.md. The downloaded traces
are cached but not tracked; fetch.py regenerates them.
14 changed files+1248−1
.gitignoremodified+1−1View file
@@ -2,4 +2,4 @@ node_modules
22 dist
33 __pycache__
44 tmp
5-paper
5+exploration/ephys/cache
exploration/ephys/RESULTS.mdadded+212−0View file
@@ -0,0 +1,212 @@
1+# Does the quantized-filtered-Gaussian model predict real ephys compressibility?
2+
3+Test: fit `x ~ iid N(0,1) → h * x → round` to a real ephys trace using nothing
4+but its power spectrum, synthesize a surrogate of the same length from the fit,
5+and run the same 19 codecs on both. If the model is a good stand-in, every
6+codec should land at the same bits/sample on the surrogate as on the real
7+trace.
8+
9+Data: the three single-channel traces benchcompress benchmarks — DANDI 000876
10+ch45, 000409 ch101, 001290 ch0 — 200k samples each, raw and in the `-filtered`
11+variant (bandpass 300–6000 Hz, normalized to MAD noise = 1, requantized at step
12+v = 0.25). Scripts: `fetch.py`, `fitmodel.py`, `codec_suite.py`, `compare.py`,
13+`sweep.py`. Logs: `full_compare.log`, `sweep.log`.
14+
15+## Verdict
16+
17+**The model is a good surrogate for compressibility, and a bad surrogate for
18+the signal.** The surrogate has kurtosis ~3.0 against 4.5–19.2 for the real
19+traces — it has no spikes at all — yet it reproduces the measured rate of every
20+prediction-based codec to within **0.05 bits/sample**, and of the generic byte
21+compressors to within ~0.05 bits/sample on the filtered variants. Compression
22+of these recordings is set by the noise spectrum, not by the spikes.
23+
24+**But the interesting regime is not where ephys currently sits.** At the
25+deployed quantization (noise = 4 steps, `s_* ≈ 1`) the fractional-phase loss is
26+only 0.06–0.08 bits/sample (~3%), and LPC+ANS is within ~6% of the analytic
27+entropy rate of the fitted model. The loss becomes large only when the
28+quantizer is coarsened to ~0.5 steps of noise — and the model predicts exactly
29+where that happens.
30+
31+## 1. Model vs real, at the benchcompress operating points
32+
33+`-filtered` variants (bandpass, noise = 4 quantization steps). bits/sample.
34+
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* |
48+
49+Raw (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
51+up to 0.5 bits) — those are the ones most sensitive to the heavy tails the
52+Gaussian model does not have.
53+
54+One visible fit artifact: on 000876 raw the surrogate std is 26.1 against 47.9
55+real. Welch detrends each segment, so the fit discards drift below ~7 Hz. That
56+content is almost perfectly predictable, so it moves `order-0 H(z)` by 0.77
57+bits and `LPC(32)+ANS` by 0.01 bits.
58+
59+## 2. Walking the s\* axis by coarsening the quantizer
60+
61+Same 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
63+throwing away the fractional part of the prediction. 200k samples.
64+(001290 ch0; the other two traces agree to ~0.05 bits — see `sweep.log`.)
65+
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 |
75+
76+Three things to read off this:
77+
78+1. **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.
83+2. **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.
86+3. **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.**
94+
95+## 1b. LPC+ANS beats FLAC by ~21% at the deployed settings
96+
97+FLAC is the mature instance of the same architecture — LPC prediction, integer
98+residual, memoryless coding — differing only in that it Rice-codes the residual
99+where ANS codes its empirical histogram. That one substitution is worth a fifth
100+of the file:
101+
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 |
109+
110+Read as distance above the residual-entropy ceiling: **FLAC sits 29–31% above
111+it, LPC+ANS 1.5–2.1%.** The loss is entirely in the residual coder, not the
112+predictor. On raw traces (`s_*` = 3.8–6.8) the advantage shrinks to 1–3%: Rice
113+coding is near-optimal for wide residuals and mismatched for narrow ones, which
114+is the same $s_*$ story from a different direction.
115+
116+Caveat: libsndfile does not expose the FLAC compression level, so this is its
117+default (level 5). `flac -8` raises the max LPC order and improves residual
118+partitioning — worth a few percent, not twenty. **The paper must re-measure
119+with the reference `flac -8` binary and with WavPack** (whose prebuilt wheels
120+do not match this machine's glibc).
121+
122+## 2b. The model does not need an explicit filtering step
123+
124+`h` in the fit is the minimum-phase spectral factor of whatever spectrum the
125+trace has — not a filter anyone applied. Any stationary Gaussian process with
126+spectral density `S` can be written `h * x`. Two checks that this is not just a
127+formal remark:
128+
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:
137+
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 |
145+
146+So the preconditions are **stationarity over the block**, **second-order
147+statistics sufficing for the rate** (validated at kurtosis 19), and **noise
148+≳ 0.5 quantization steps**. Explicit filtering is not among them.
149+
150+## 2c. The gap decomposes into named terms, and the structural one has a
151+closed form
152+
153+Write the distance from a measured codec rate to the analytic entropy rate as
154+
155+```
156+R_ANS - G(s*) = [R_ANS - H0] + [H0 - R_condG] + [R_condG - G(s*)]
157+ coder residual-model prediction
158+```
159+
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.
170+
171+Measured against predicted, three recordings, 200k samples (`decompose.log`):
172+
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 |
180+
181+Read honestly: **the closed form gets the structural term right to within about
182+a 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
184+Gaussian the theory assumes. The coder term is a stable 1–2% of the rate, as
185+claimed. 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.
187+
188+## 3. Consequences for the paper
189+
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.
205+
206+## Open
207+
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.
exploration/ephys/codec_suite.pyadded+183−0View file
@@ -0,0 +1,183 @@
1+"""Lossless codecs measured in bits/sample on an int16 trace.
2+
3+Every entry returns the full encoded size including whatever the decoder needs.
4+The generic byte compressors are run on the raw int16 buffer, on the
5+int16-wrapped first difference, and on the byte-split ("shuffled") buffer that
6+Blosc-style pipelines use.
7+"""
8+import bz2
9+import lzma
10+import math
11+import zlib
12+
13+import numpy as np
14+from scipy.linalg import solve_toeplitz
15+from scipy.signal import lfilter
16+from scipy.special import ndtr
17+
18+COEFF_PRECISION = 15
19+
20+
21+# ------------------------------------------------------------- preprocessing
22+def as_int16(z):
23+ return np.asarray(z, dtype=np.int16)
24+
25+
26+def delta(z):
27+ d = np.empty_like(z)
28+ d[0] = z[0]
29+ d[1:] = (z[1:].astype(np.int32) - z[:-1].astype(np.int32)).astype(np.int16)
30+ return d
31+
32+
33+def byteshuffle(z):
34+ b = z.tobytes()
35+ a = np.frombuffer(b, dtype=np.uint8).reshape(-1, 2)
36+ return np.concatenate([a[:, 0], a[:, 1]]).tobytes()
37+
38+
39+# ------------------------------------------------------- generic compressors
40+def _zlib(b):
41+ return len(zlib.compress(b, 9))
42+
43+
44+def _lzma(b):
45+ return len(lzma.compress(b, preset=9 | lzma.PRESET_EXTREME))
46+
47+
48+def _bz2(b):
49+ return len(bz2.compress(b, 9))
50+
51+
52+def _zstd(b):
53+ import zstandard
54+ return len(zstandard.ZstdCompressor(level=19).compress(b))
55+
56+
57+def _brotli(b):
58+ import brotli
59+ return len(brotli.compress(b, quality=11))
60+
61+
62+def _lz4(b):
63+ import lz4.frame
64+ return len(lz4.frame.compress(b, compression_level=12))
65+
66+
67+GENERIC = {
68+ "zlib": _zlib, "zstd": _zstd, "lzma": _lzma,
69+ "bz2": _bz2, "brotli": _brotli, "lz4": _lz4,
70+}
71+
72+
73+# ------------------------------------------------------- lossless audio codecs
74+def flac_bytes(z, rate=30000):
75+ """FLAC via libsndfile: LPC + Rice-coded integer residual, the mature
76+ instance of the architecture LPC+ANS also belongs to. libsndfile does not
77+ expose the compression level; this is its default (level 5)."""
78+ import io
79+ import soundfile as sf
80+ buf = io.BytesIO()
81+ sf.write(buf, np.asarray(z, dtype=np.int16), int(rate),
82+ format="FLAC", subtype="PCM_16")
83+ return buf.getbuffer().nbytes
84+
85+
86+# --------------------------------------------------------------------- rates
87+def entropy0(v):
88+ counts = np.unique(np.asarray(v), return_counts=True)[1]
89+ p = counts / counts.sum()
90+ return float(-(p * np.log2(p)).sum())
91+
92+
93+# ------------------------------------------------ integer LPC (the app's own)
94+def autocorr(z, order, ridge=1e-8):
95+ """Autocorrelation lags 0..order, with a small ridge on lag 0 so Levinson
96+ stays non-singular on degenerate (near-constant, near-empty) blocks."""
97+ zf = np.asarray(z, dtype=np.float64)
98+ r = np.array([zf @ zf if lag == 0 else zf[lag:] @ zf[:-lag]
99+ for lag in range(order + 1)])
100+ r[0] = max(r[0], 1e-12 * zf.size) * (1.0 + ridge)
101+ return r
102+
103+
104+def fit_lpc_quantized(z, order):
105+ r = autocorr(z, order)
106+ a = solve_toeplitz(r[:order], r[1:order + 1])
107+ peak = float(np.abs(a).max())
108+ if peak <= 0:
109+ return np.zeros(order, dtype=np.int64), 0
110+ shift = COEFF_PRECISION - 1 - int(np.floor(np.log2(peak))) - 1
111+ shift = max(0, min(15, shift))
112+ limit = 2 ** (COEFF_PRECISION - 1)
113+ q = np.clip(np.round(a * 2.0 ** shift), -limit, limit - 1).astype(np.int64)
114+ return q, shift
115+
116+
117+def lpc_residual(z, q, shift):
118+ z = np.asarray(z, dtype=np.int64)
119+ order = len(q)
120+ acc = np.convolve(z, q, mode="full")[:len(z)]
121+ pred = np.zeros_like(z)
122+ pred[1:] = acc[:-1] >> shift # floor division
123+ e = z.copy()
124+ e[order:] = z[order:] - pred[order:]
125+ return (((e + (1 << 15)) & 0xFFFF) - (1 << 15)).astype(np.int16)
126+
127+
128+def lpc_ans_bytes(z, order):
129+ """Real encoded size of integer-LPC + rANS, side information included."""
130+ import simple_ans
131+ q, shift = fit_lpc_quantized(z, order)
132+ e = lpc_residual(z, q, shift)
133+ enc = simple_ans.ans_encode(e)
134+ payload = 4 * enc.words.size
135+ table = 4 * enc.symbol_counts.size + 2 * enc.symbol_values.size
136+ header = 2 * order + 4 + 8 # coefficients, shift, n
137+ return payload + table + header, entropy0(e)
138+
139+
140+# ------------------------------- conditional-Gaussian coding (achievable rate)
141+def conditional_gaussian_rate(z, order):
142+ """Cross-entropy of the real-valued-prediction conditional-Gaussian model,
143+ in bits/sample, plus header cost. This is what the arithmetic coder of
144+ exploration/codec_gaussian.py achieves to within ~0.1%."""
145+ zf = np.asarray(z, dtype=np.float64)
146+ n = zf.size
147+ r = autocorr(zf, order) / n
148+ a = solve_toeplitz(r[:order], r[1:order + 1])
149+ pred = lfilter(np.concatenate(([0.0], a)), [1.0], zf)
150+ zt, mu = zf[order:], pred[order:]
151+ s0 = math.sqrt(max((zt - mu).var() - 1.0 / 12.0, 1e-6))
152+ best = np.inf
153+ for s in s0 * np.linspace(0.7, 1.3, 25):
154+ p = ndtr((zt + 0.5 - mu) / s) - ndtr((zt - 0.5 - mu) / s)
155+ best = min(best, float(-np.log2(np.maximum(p, 1e-12)).mean()))
156+ header_bits = 8 * (4 * order + 14) / n
157+ return best + header_bits
158+
159+
160+# ------------------------------------------------------------------ the suite
161+def measure(z, lpc_order=32, generic=True):
162+ """bits/sample for every method, as an ordered dict."""
163+ z = as_int16(z)
164+ n = z.size
165+ out = {}
166+ out["raw int16"] = 16.0
167+ out["order-0 H(z)"] = entropy0(z)
168+ if generic:
169+ buffers = {"": z.tobytes(), "+delta": delta(z).tobytes(),
170+ "+shuffle": byteshuffle(z)}
171+ for name, fn in GENERIC.items():
172+ for suffix, buf in buffers.items():
173+ out[f"{name}{suffix}"] = 8.0 * fn(buf) / n
174+ try:
175+ out["FLAC"] = 8.0 * flac_bytes(z) / n
176+ except Exception as exc: # pragma: no cover
177+ out["FLAC"] = float("nan")
178+ print(f" (FLAC unavailable: {exc})")
179+ nbytes, resid_h0 = lpc_ans_bytes(z, lpc_order)
180+ out[f"LPC({lpc_order})+ANS"] = 8.0 * nbytes / n
181+ out[f"LPC({lpc_order}) resid H0"] = resid_h0
182+ out[f"cond-Gauss({lpc_order})"] = conditional_gaussian_rate(z, lpc_order)
183+ return out
exploration/ephys/compare.pyadded+81−0View file
@@ -0,0 +1,81 @@
1+"""Does the paper's model explain real ephys compressibility?
2+
3+For each cached trace: fit x ~ N(0,1) -> h * x -> round to its spectrum,
4+synthesize a surrogate of the same length from the fit, and run the same codec
5+suite on both. If the model is a good stand-in, every codec should land at the
6+same bits/sample on the surrogate as on the real trace.
7+
8+Usage: python compare.py [--n N] [--nfft NFFT] [--taps T] [--order O] [--fast]
9+"""
10+import argparse
11+import glob
12+import os
13+import sys
14+import time
15+
16+import numpy as np
17+
18+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
19+import codec_suite as cc # noqa: E402
20+from fitmodel import Fit # noqa: E402
21+
22+CACHE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
23+
24+
25+def 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)
28+
29+ fit = Fit(z_real, nfft=args.nfft, n_taps=args.taps)
30+ z_syn = fit.synthesize(seed=0)
31+
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)")
42+
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]")
47+
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}
56+
57+
58+def 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()
68+
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")
74+
75+ for path in paths:
76+ name = os.path.basename(path)[: -len(".npy")]
77+ report(name, np.load(path), args)
78+
79+
80+if __name__ == "__main__":
81+ main()
exploration/ephys/decompose.logadded+39−0View file
@@ -0,0 +1,39 @@
1+ecephys-000409-ch101.raw.npy n = 200000 order = 32
2+
3+ v s_* s_fit | R_ANS H0 R_cG G(s*) | coder resid pred | L(s_fit) L(s_*)
4+ --------------------------------------------------------------------------------------------
5+ 0.12 1.359 1.402 | 2.613 2.561 2.569 2.521 | 0.052 -0.009 0.048 | 0.029 0.030
6+ 0.25 0.922 0.954 | 2.102 2.073 2.048 1.997 | 0.029 0.025 0.051 | 0.058 0.062
7+ 0.50 0.634 0.654 | 1.677 1.649 1.570 1.526 | 0.028 0.078 0.044 | 0.109 0.114
8+ 1.00 0.433 0.453 | 1.417 1.392 1.120 1.105 | 0.024 0.272 0.015 | 0.181 0.190
9+ 2.00 0.283 0.236 | 1.208 1.180 0.688 0.736 | 0.027 0.492 -0.047 | 0.271 0.261
10+ 4.00 0.114 0.001 | 0.631 0.622 0.884 0.297 | 0.009 -0.262 0.587 | 0.010 0.233
11+ 8.00 0.001 0.001 | 0.087 0.082 0.143 0.002 | 0.005 -0.061 0.141 | 0.010 0.007
12+ 16.00 0.000 0.001 | 0.069 0.064 0.054 0.001 | 0.005 0.010 0.053 | 0.006 0.004
13+
14+ecephys-000876-ch45.raw.npy n = 200000 order = 32
15+
16+ v s_* s_fit | R_ANS H0 R_cG G(s*) | coder resid pred | L(s_fit) L(s_*)
17+ --------------------------------------------------------------------------------------------
18+ 0.12 1.514 1.544 | 2.774 2.720 2.704 2.672 | 0.054 0.016 0.032 | 0.024 0.025
19+ 0.25 1.023 1.049 | 2.263 2.217 2.175 2.136 | 0.046 0.042 0.039 | 0.049 0.051
20+ 0.50 0.708 0.724 | 1.806 1.778 1.693 1.660 | 0.028 0.085 0.033 | 0.093 0.096
21+ 1.00 0.487 0.490 | 1.469 1.447 1.240 1.225 | 0.022 0.207 0.015 | 0.165 0.166
22+ 2.00 0.327 0.336 | 1.243 1.227 0.777 0.848 | 0.016 0.451 -0.071 | 0.239 0.243
23+ 4.00 0.156 0.001 | 0.938 0.918 2.100 0.408 | 0.020 -1.182 1.692 | 0.010 0.260
24+ 8.00 0.001 0.001 | 0.288 0.280 0.231 0.002 | 0.008 0.050 0.229 | 0.010 0.007
25+ 16.00 0.000 0.001 | 0.011 0.007 0.009 0.000 | 0.004 -0.002 0.009 | 0.006 0.003
26+
27+ecephys-001290-ch0.raw.npy n = 200000 order = 32
28+
29+ v s_* s_fit | R_ANS H0 R_cG G(s*) | coder resid pred | L(s_fit) L(s_*)
30+ --------------------------------------------------------------------------------------------
31+ 0.12 1.658 1.724 | 2.930 2.872 2.859 2.798 | 0.059 0.013 0.061 | 0.019 0.021
32+ 0.25 1.099 1.147 | 2.360 2.328 2.295 2.232 | 0.032 0.033 0.063 | 0.042 0.045
33+ 0.50 0.741 0.770 | 1.872 1.849 1.770 1.717 | 0.023 0.078 0.054 | 0.084 0.089
34+ 1.00 0.500 0.517 | 1.489 1.455 1.296 1.255 | 0.034 0.159 0.042 | 0.153 0.160
35+ 2.00 0.325 0.345 | 1.260 1.232 0.836 0.842 | 0.028 0.396 -0.006 | 0.235 0.244
36+ 4.00 0.103 0.001 | 0.951 0.930 1.938 0.267 | 0.021 -1.008 1.670 | 0.010 0.223
37+ 8.00 0.001 0.001 | 0.102 0.096 0.069 0.001 | 0.005 0.027 0.068 | 0.006 0.005
38+ 16.00 0.000 0.001 | 0.015 0.011 0.011 0.000 | 0.004 0.000 0.010 | 0.006 0.003
39+
exploration/ephys/decompose.pyadded+79−0View file
@@ -0,0 +1,79 @@
1+"""Decompose the measured gap between LPC+ANS and the entropy rate into named
2+terms, and check each against what the theory predicts.
3+
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
13+
14+`s_fit` is the residual scale the conditional-Gaussian model actually fits,
15+which is the honest argument to L(.) — `s_*` is its prediction from the
16+spectrum, and the two differing is itself informative.
17+"""
18+import os
19+import sys
20+
21+import numpy as np
22+from scipy.linalg import solve_toeplitz
23+from scipy.signal import lfilter
24+from scipy.special import ndtr
25+
26+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
27+import codec_suite as cc # noqa: E402
28+from fitmodel import Fit # noqa: E402
29+from phase_loss import phase_entropies # noqa: E402
30+from sweep import noise_units, STEPS # noqa: E402
31+
32+CACHE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
33+
34+
35+def 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
51+
52+
53+def 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}")
76+
77+
78+if __name__ == "__main__":
79+ main(trace=sys.argv[1] if len(sys.argv) > 1 else "001290")
exploration/ephys/fetch.pyadded+79−0View file
@@ -0,0 +1,79 @@
1+"""Download single-channel ephys traces from the benchcompress datasets and
2+cache them as .npy, together with the bandpass+requantized "filtered" variant
3+that benchcompress benchmarks.
4+
5+Usage: python fetch.py [outdir]
6+"""
7+import os
8+import sys
9+
10+import numpy as np
11+
12+CACHE = sys.argv[1] if len(sys.argv) > 1 else os.path.join(
13+ os.path.dirname(os.path.abspath(__file__)), "cache")
14+
15+NUM_SAMPLES = 500_000
16+RATE = 30000.0
17+
18+# (name, dandi asset url, dataset path in the nwb, channel)
19+SOURCES = [
20+ ("ecephys-000876-ch45",
21+ "https://api.dandiarchive.org/api/assets/7e1de06d-d478-40e2-9b64-9dd04eafaa4c/download/",
22+ "/acquisition/ElectricalSeriesAP/data", 45),
23+ ("ecephys-000409-ch101",
24+ "https://api.dandiarchive.org/api/assets/c04f6b30-82bf-40e1-9210-34f0bcd8be24/download/",
25+ "/acquisition/ElectricalSeriesAp/data", 101),
26+ ("ecephys-001290-ch0",
27+ "https://api.dandiarchive.org/api/assets/78c99d23-da88-4ecd-9086-c488a126eac5/download/",
28+ "/acquisition/ElectricalSeriesAPImec/data", 0),
29+]
30+
31+
32+def bandpass(x, lowcut, highcut, rate):
33+ from scipy.signal import butter, lfilter
34+ nyq = 0.5 * rate
35+ b, a = butter(5, [lowcut / nyq, highcut / nyq], btype="band")
36+ return lfilter(b, a, x)
37+
38+
39+def highpass(x, lowcut, rate):
40+ from scipy.signal import butter, lfilter
41+ nyq = 0.5 * rate
42+ b, a = butter(5, lowcut / nyq, btype="high")
43+ return lfilter(b, a, x)
44+
45+
46+def noise_level(x, rate):
47+ xf = highpass(x, 300.0, rate)
48+ return float(np.median(np.abs(xf - np.median(xf))) / 0.6745)
49+
50+
51+def filtered_variant(x, rate=RATE, v=0.25, lowcut=300.0, highcut=6000.0):
52+ """benchcompress's `-filtered` transform: bandpass, normalize by the MAD
53+ noise level, requantize at step v (so the noise std is ~1/v = 4 steps)."""
54+ xf = bandpass(x - np.median(x), lowcut, highcut, rate)
55+ nl = noise_level(xf, rate)
56+ return np.round(xf / nl / v).astype(np.int16)
57+
58+
59+def main():
60+ os.makedirs(CACHE, exist_ok=True)
61+ for name, url, path, ch in SOURCES:
62+ raw_path = os.path.join(CACHE, f"{name}.raw.npy")
63+ filt_path = os.path.join(CACHE, f"{name}.filtered.npy")
64+ if os.path.exists(raw_path) and os.path.exists(filt_path):
65+ print(f"{name}: cached")
66+ continue
67+ print(f"{name}: downloading {NUM_SAMPLES} samples ...", flush=True)
68+ import lindi
69+ h5f = lindi.LindiH5pyFile.from_hdf5_file(url)
70+ ds = h5f[path]
71+ raw = np.asarray(ds[:NUM_SAMPLES, ch : ch + 1]).flatten().astype(np.int16)
72+ np.save(raw_path, raw)
73+ np.save(filt_path, filtered_variant(raw.astype(np.float64)))
74+ print(f"{name}: raw std={raw.std():.1f} "
75+ f"filtered std={np.load(filt_path).std():.2f}", flush=True)
76+
77+
78+if __name__ == "__main__":
79+ main()
exploration/ephys/fitmodel.pyadded+85−0View file
@@ -0,0 +1,85 @@
1+"""Fit the paper's generative model x ~ iid N(0, 1) -> h * x -> round to a
2+real integer trace, and synthesize a surrogate from the fit.
3+
4+Only the spectrum is identifiable (sigma and the kernel scale are the same
5+knob), so sigma is fixed at 1 and h carries everything.
6+
7+Fit:
8+ 1. Welch PSD of z, in the convention mean_{f in [0,1)} S(f) = var(z).
9+ 2. Dither model: S_z = S_y + 1/12, so S_y = max(S_z - 1/12, floor).
10+ 3. h = the minimum-phase spectral factor of S_y (real-cepstrum method).
11+
12+Prediction (no sampling required):
13+ s_*^2 = exp( mean_f ln(S_y + 1/12) ) - 1/12 [Szego on the observed
14+ spectrum], and the predicted entropy rate is G(s_*), the entropy of
15+ N(0, s_*^2) + U(-1/2, 1/2).
16+"""
17+import numpy as np
18+from scipy.signal import welch, fftconvolve
19+
20+FLOOR_FRAC = 1e-4 # S_y is floored at this fraction of its own mean
21+
22+
23+def psd(z, nfft=4096):
24+ """Two-sided PSD on the rfft grid f = k/nfft, k = 0..nfft/2, normalized so
25+ that the mean over f in [0, 1) equals var(z)."""
26+ f, pxx = welch(np.asarray(z, dtype=np.float64), fs=1.0, nperseg=nfft,
27+ noverlap=nfft // 2, window="hann", detrend="constant",
28+ return_onesided=True, scaling="density")
29+ return f, pxx / 2.0
30+
31+
32+def _grid_mean(v):
33+ """Mean over f in [0, 1) of a quantity given on the rfft half-grid
34+ (trapezoid on [0, 1/2]; the spectrum is symmetric)."""
35+ return float((v.sum() - 0.5 * (v[0] + v[-1])) / (v.size - 1))
36+
37+
38+def minimum_phase(s_half, n_taps=None):
39+ """Minimum-phase impulse response whose |H(f)|^2 matches s_half on the
40+ rfft grid. Real-cepstrum construction."""
41+ s_full = np.concatenate([s_half, s_half[-2:0:-1]])
42+ n = s_full.size
43+ c = np.fft.ifft(0.5 * np.log(s_full)).real # real cepstrum of |H|
44+ cm = np.zeros(n)
45+ cm[0] = c[0]
46+ cm[1:n // 2] = 2.0 * c[1:n // 2]
47+ cm[n // 2] = c[n // 2]
48+ h = np.fft.ifft(np.exp(np.fft.fft(cm))).real # causal, min phase
49+ return h if n_taps is None else h[:n_taps]
50+
51+
52+class Fit:
53+ def __init__(self, z, nfft=4096, n_taps=None):
54+ z = np.asarray(z, dtype=np.float64)
55+ self.n = z.size
56+ self.nfft = nfft
57+ self.freq, self.s_z = psd(z, nfft)
58+ floor = FLOOR_FRAC * self.s_z.mean()
59+ self.s_y = np.maximum(self.s_z - 1.0 / 12.0, floor)
60+ self.kernel = minimum_phase(self.s_y, n_taps)
61+ # Szego on the observed spectrum, with roundoff as a 1/12 dither floor
62+ gm = np.exp(_grid_mean(np.log(self.s_y + 1.0 / 12.0)))
63+ self.s_star = float(np.sqrt(max(gm - 1.0 / 12.0, 0.0)))
64+ self.sigma_inf = float(np.sqrt(np.exp(_grid_mean(np.log(self.s_y)))))
65+
66+ @property
67+ def predicted_rate(self):
68+ """G(s_*) in bits/sample — the analytic entropy-rate prediction."""
69+ from timeseries_entropy.theory import gauss_uniform_entropy
70+ return float(gauss_uniform_entropy(self.s_star))
71+
72+ def kernel_error_db(self):
73+ """How well the (possibly truncated) kernel reproduces the fitted
74+ spectrum: RMS error in dB over the grid."""
75+ h = np.abs(np.fft.rfft(self.kernel, self.nfft)) ** 2
76+ return float(np.sqrt(np.mean((10 * np.log10(h / self.s_y)) ** 2)))
77+
78+ def synthesize(self, n=None, seed=0):
79+ """A surrogate integer trace from the fitted model."""
80+ n = self.n if n is None else n
81+ rng = np.random.default_rng(seed)
82+ k = self.kernel
83+ x = rng.standard_normal(n + k.size - 1)
84+ y = fftconvolve(x, k, mode="valid")
85+ return np.round(y).astype(np.int16)
exploration/ephys/full_compare.logadded+216−0View file
@@ -0,0 +1,216 @@
1+
2+==============================================================================
3+ecephys-000409-ch101.filtered n = 200000
4+==============================================================================
5+ real: std 4.104 kurtosis 17.70
6+ surrogate: std 4.064 kurtosis 3.03
7+ fit: taps 4096 spectrum RMS error 0.00 dB
8+ s_* = 0.9076 quantization steps -> predicted entropy rate G(s_*) = 1.9767 bits/sample (8.1x)
9+ [9.3s]
10+
11+ method real model diff real x model x
12+ ------------------------------------------------------------------
13+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
14+ order-0 H(z) 3.8712 4.0698 +0.1986 4.13 3.93
15+ zlib 3.5021 3.6375 +0.1354 4.57 4.40
16+ zlib+delta 3.0835 3.1808 +0.0973 5.19 5.03
17+ zlib+shuffle 3.8026 3.8942 +0.0916 4.21 4.11
18+ zstd 2.9785 3.1001 +0.1216 5.37 5.16
19+ zstd+delta 2.7010 2.8254 +0.1244 5.92 5.66
20+ zstd+shuffle 3.4699 3.5757 +0.1058 4.61 4.47
21+ lzma 2.7819 2.9048 +0.1229 5.75 5.51
22+ lzma+delta 2.5067 2.6013 +0.0946 6.38 6.15
23+ lzma+shuffle 3.3894 3.5037 +0.1142 4.72 4.57
24+ bz2 2.5189 2.5797 +0.0608 6.35 6.20
25+ bz2+delta 2.6231 2.7260 +0.1029 6.10 5.87
26+ bz2+shuffle 3.1706 3.1916 +0.0210 5.05 5.01
27+ brotli 2.8846 3.0133 +0.1286 5.55 5.31
28+ brotli+delta 2.6920 2.8064 +0.1144 5.94 5.70
29+ brotli+shuffle 3.4607 3.5651 +0.1044 4.62 4.49
30+ lz4 3.8891 4.0757 +0.1866 4.11 3.93
31+ lz4+delta 3.5282 3.6736 +0.1454 4.53 4.36
32+ lz4+shuffle 4.5850 4.7857 +0.2007 3.49 3.34
33+ FLAC 2.6798 2.7252 +0.0454 5.97 5.87
34+ LPC(32)+ANS 2.0827 2.1277 +0.0450 7.68 7.52
35+ LPC(32) resid H0 2.0528 2.0938 +0.0410 7.79 7.64
36+ cond-Gauss(32) 2.0288 2.0390 +0.0102 7.89 7.85
37+
38+==============================================================================
39+ecephys-000409-ch101.raw n = 200000
40+==============================================================================
41+ real: std 14.440 kurtosis 19.18
42+ surrogate: std 14.347 kurtosis 2.98
43+ fit: taps 4096 spectrum RMS error 0.00 dB
44+ s_* = 6.7593 quantization steps -> predicted entropy rate G(s_*) = 4.8053 bits/sample (3.3x)
45+ [6.9s]
46+
47+ method real model diff real x model x
48+ ------------------------------------------------------------------
49+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
50+ order-0 H(z) 5.5107 5.8893 +0.3786 2.90 2.72
51+ zlib 6.6606 7.0779 +0.4173 2.40 2.26
52+ zlib+delta 6.5159 6.6210 +0.1050 2.46 2.42
53+ zlib+shuffle 6.5794 6.8628 +0.2833 2.43 2.33
54+ zstd 6.1070 6.4985 +0.3915 2.62 2.46
55+ zstd+delta 6.0940 6.1588 +0.0649 2.63 2.60
56+ zstd+shuffle 6.3090 6.6619 +0.3529 2.54 2.40
57+ lzma 5.5096 5.8578 +0.3482 2.90 2.73
58+ lzma+delta 5.6008 5.6872 +0.0864 2.86 2.81
59+ lzma+shuffle 5.8715 6.1821 +0.3106 2.73 2.59
60+ bz2 5.1064 5.3492 +0.2428 3.13 2.99
61+ bz2+delta 5.3751 5.4582 +0.0831 2.98 2.93
62+ bz2+shuffle 5.8671 6.0611 +0.1940 2.73 2.64
63+ brotli 5.5399 5.5734 +0.0334 2.89 2.87
64+ brotli+delta 5.2531 5.2589 +0.0058 3.05 3.04
65+ brotli+shuffle 5.7376 5.8984 +0.1608 2.79 2.71
66+ lz4 7.9504 8.4726 +0.5222 2.01 1.89
67+ lz4+delta 7.9536 8.1194 +0.1658 2.01 1.97
68+ lz4+shuffle 8.3784 8.8310 +0.4526 1.91 1.81
69+ FLAC 4.9060 4.9645 +0.0585 3.26 3.22
70+ LPC(32)+ANS 4.8550 4.9015 +0.0465 3.30 3.26
71+ LPC(32) resid H0 4.7790 4.8271 +0.0481 3.35 3.31
72+ cond-Gauss(32) 4.8181 4.8317 +0.0136 3.32 3.31
73+
74+==============================================================================
75+ecephys-000876-ch45.filtered n = 200000
76+==============================================================================
77+ real: std 4.756 kurtosis 5.98
78+ surrogate: std 4.759 kurtosis 3.03
79+ fit: taps 4096 spectrum RMS error 0.00 dB
80+ s_* = 1.0220 quantization steps -> predicted entropy rate G(s_*) = 2.1338 bits/sample (7.5x)
81+ [8.9s]
82+
83+ method real model diff real x model x
84+ ------------------------------------------------------------------
85+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
86+ order-0 H(z) 4.2362 4.2975 +0.0612 3.78 3.72
87+ zlib 4.0053 4.0200 +0.0147 3.99 3.98
88+ zlib+delta 3.4830 3.5000 +0.0170 4.59 4.57
89+ zlib+shuffle 4.1611 4.1991 +0.0380 3.85 3.81
90+ zstd 3.4278 3.4267 -0.0012 4.67 4.67
91+ zstd+delta 3.0654 3.1046 +0.0392 5.22 5.15
92+ zstd+shuffle 3.8712 3.8989 +0.0277 4.13 4.10
93+ lzma 3.1691 3.1901 +0.0210 5.05 5.02
94+ lzma+delta 2.8229 2.8539 +0.0310 5.67 5.61
95+ lzma+shuffle 3.7782 3.8112 +0.0330 4.23 4.20
96+ bz2 2.7516 2.7393 -0.0123 5.81 5.84
97+ bz2+delta 2.8424 2.8723 +0.0299 5.63 5.57
98+ bz2+shuffle 3.4249 3.4278 +0.0029 4.67 4.67
99+ brotli 3.2793 3.3048 +0.0254 4.88 4.84
100+ brotli+delta 3.0191 3.0434 +0.0243 5.30 5.26
101+ brotli+shuffle 3.8470 3.8901 +0.0431 4.16 4.11
102+ lz4 4.4432 4.5126 +0.0694 3.60 3.55
103+ lz4+delta 3.9691 4.0342 +0.0651 4.03 3.97
104+ lz4+shuffle 5.1227 5.2399 +0.1172 3.12 3.05
105+ FLAC 2.8648 2.9087 +0.0439 5.58 5.50
106+ LPC(32)+ANS 2.2614 2.2702 +0.0088 7.08 7.05
107+ LPC(32) resid H0 2.2151 2.2324 +0.0173 7.22 7.17
108+ cond-Gauss(32) 2.1731 2.1898 +0.0167 7.36 7.31
109+
110+==============================================================================
111+ecephys-000876-ch45.raw n = 200000
112+==============================================================================
113+ real: std 47.855 kurtosis 4.49
114+ surrogate: std 26.107 kurtosis 2.90
115+ fit: taps 4096 spectrum RMS error 0.00 dB
116+ s_* = 5.2245 quantization steps -> predicted entropy rate G(s_*) = 4.4346 bits/sample (3.6x)
117+ [5.4s]
118+
119+ method real model diff real x model x
120+ ------------------------------------------------------------------
121+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
122+ order-0 H(z) 7.5227 6.7494 -0.7734 2.13 2.37
123+ zlib 7.4830 7.2800 -0.2029 2.14 2.20
124+ zlib+delta 6.1132 6.1128 -0.0003 2.62 2.62
125+ zlib+shuffle 7.0132 7.0086 -0.0046 2.28 2.28
126+ zstd 6.5590 6.5471 -0.0119 2.44 2.44
127+ zstd+delta 5.6966 5.7258 +0.0292 2.81 2.79
128+ zstd+shuffle 6.3085 6.3790 +0.0705 2.54 2.51
129+ lzma 5.6602 5.7149 +0.0547 2.83 2.80
130+ lzma+delta 5.1886 5.2382 +0.0496 3.08 3.05
131+ lzma+shuffle 5.3280 5.6688 +0.3408 3.00 2.82
132+ bz2 4.9190 4.9421 +0.0231 3.25 3.24
133+ bz2+delta 5.0644 5.0874 +0.0230 3.16 3.15
134+ bz2+shuffle 5.1363 5.3262 +0.1899 3.12 3.00
135+ brotli 5.8924 5.8680 -0.0244 2.72 2.73
136+ brotli+delta 4.9019 4.8860 -0.0160 3.26 3.27
137+ brotli+shuffle 5.5244 5.5517 +0.0272 2.90 2.88
138+ lz4 8.8184 8.5224 -0.2960 1.81 1.88
139+ lz4+delta 7.4188 7.5147 +0.0959 2.16 2.13
140+ lz4+shuffle 8.0256 8.2186 +0.1930 1.99 1.95
141+ FLAC 4.6271 4.6990 +0.0719 3.46 3.41
142+ LPC(32)+ANS 4.5216 4.5119 -0.0097 3.54 3.55
143+ LPC(32) resid H0 4.4219 4.4508 +0.0290 3.62 3.59
144+ cond-Gauss(32) 4.4463 4.4525 +0.0062 3.60 3.59
145+
146+==============================================================================
147+ecephys-001290-ch0.filtered n = 200000
148+==============================================================================
149+ real: std 4.322 kurtosis 5.87
150+ surrogate: std 4.296 kurtosis 3.03
151+ fit: taps 4096 spectrum RMS error 0.00 dB
152+ s_* = 1.1039 quantization steps -> predicted entropy rate G(s_*) = 2.2374 bits/sample (7.2x)
153+ [8.7s]
154+
155+ method real model diff real x model x
156+ ------------------------------------------------------------------
157+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
158+ order-0 H(z) 4.1290 4.1499 +0.0209 3.87 3.86
159+ zlib 4.1723 4.1764 +0.0041 3.83 3.83
160+ zlib+delta 3.7257 3.7364 +0.0107 4.29 4.28
161+ zlib+shuffle 4.4166 4.4142 -0.0025 3.62 3.62
162+ zstd 3.6014 3.5934 -0.0079 4.44 4.45
163+ zstd+delta 3.3478 3.3034 -0.0444 4.78 4.84
164+ zstd+shuffle 4.1235 4.1213 -0.0022 3.88 3.88
165+ lzma 3.3408 3.3544 +0.0136 4.79 4.77
166+ lzma+delta 3.0323 3.0560 +0.0237 5.28 5.24
167+ lzma+shuffle 4.0374 4.0413 +0.0038 3.96 3.96
168+ bz2 2.9192 2.9349 +0.0157 5.48 5.45
169+ bz2+delta 2.9805 2.9936 +0.0130 5.37 5.34
170+ bz2+shuffle 3.7126 3.7067 -0.0059 4.31 4.32
171+ brotli 3.4559 3.4662 +0.0103 4.63 4.62
172+ brotli+delta 3.2234 3.2378 +0.0144 4.96 4.94
173+ brotli+shuffle 4.1373 4.1834 +0.0462 3.87 3.82
174+ lz4 4.7044 4.7448 +0.0404 3.40 3.37
175+ lz4+delta 4.2923 4.3249 +0.0326 3.73 3.70
176+ lz4+shuffle 5.5315 5.5628 +0.0313 2.89 2.88
177+ FLAC 3.0432 3.0751 +0.0319 5.26 5.20
178+ LPC(32)+ANS 2.3687 2.3914 +0.0227 6.75 6.69
179+ LPC(32) resid H0 2.3335 2.3469 +0.0134 6.86 6.82
180+ cond-Gauss(32) 2.3008 2.3109 +0.0100 6.95 6.92
181+
182+==============================================================================
183+ecephys-001290-ch0.raw n = 200000
184+==============================================================================
185+ real: std 5.327 kurtosis 6.33
186+ surrogate: std 5.284 kurtosis 3.00
187+ fit: taps 4096 spectrum RMS error 0.00 dB
188+ s_* = 3.8062 quantization steps -> predicted entropy rate G(s_*) = 3.9796 bits/sample (4.0x)
189+ [7.4s]
190+
191+ method real model diff real x model x
192+ ------------------------------------------------------------------
193+ raw int16 16.0000 16.0000 +0.0000 1.00 1.00
194+ order-0 H(z) 4.4112 4.4485 +0.0373 3.63 3.60
195+ zlib 5.6398 5.6495 +0.0097 2.84 2.83
196+ zlib+delta 5.6191 5.6258 +0.0067 2.85 2.84
197+ zlib+shuffle 5.8701 5.8890 +0.0189 2.73 2.72
198+ zstd 5.2041 5.2108 +0.0067 3.07 3.07
199+ zstd+delta 5.2158 5.2163 +0.0004 3.07 3.07
200+ zstd+shuffle 5.5396 5.5896 +0.0500 2.89 2.86
201+ lzma 4.6978 4.7408 +0.0430 3.41 3.37
202+ lzma+delta 4.7106 4.7165 +0.0059 3.40 3.39
203+ lzma+shuffle 5.2611 5.3035 +0.0424 3.04 3.02
204+ bz2 4.6195 4.6238 +0.0043 3.46 3.46
205+ bz2+delta 4.6766 4.6846 +0.0080 3.42 3.42
206+ bz2+shuffle 5.5546 5.5349 -0.0196 2.88 2.89
207+ brotli 4.4975 4.8334 +0.3359 3.56 3.31
208+ brotli+delta 4.5357 4.5447 +0.0090 3.53 3.52
209+ brotli+shuffle 5.2336 5.3565 +0.1228 3.06 2.99
210+ lz4 6.7834 6.8391 +0.0557 2.36 2.34
211+ lz4+delta 6.8257 6.8540 +0.0282 2.34 2.33
212+ lz4+shuffle 7.6468 7.6958 +0.0491 2.09 2.08
213+ FLAC 4.1928 4.2078 +0.0150 3.82 3.80
214+ LPC(32)+ANS 4.0743 4.1166 +0.0422 3.93 3.89
215+ LPC(32) resid H0 4.0199 4.0253 +0.0054 3.98 3.97
216+ cond-Gauss(32) 4.0247 4.0272 +0.0025 3.98 3.97
exploration/ephys/phase_loss.logadded+25−0View file
@@ -0,0 +1,25 @@
1+ s M(s) mixture G(s) cond L(s) loss
2+ 0.10 0.4812 0.2606 0.2206
3+ 0.20 0.7929 0.5212 0.2717
4+ 0.28 0.9900 0.7285 0.2615
5+ 0.33 1.0971 0.8556 0.2415
6+ 0.40 1.2338 1.0269 0.2069
7+ 0.50 1.4146 1.2544 0.1602
8+ 0.60 1.5844 1.4603 0.1241
9+ 0.74 1.8043 1.7149 0.0895
10+ 0.90 2.0301 1.9657 0.0643
11+ 1.00 2.1583 2.1048 0.0535
12+ 1.10 2.2777 2.2326 0.0450
13+ 1.50 2.6836 2.6583 0.0253
14+ 1.66 2.8206 2.7998 0.0209
15+ 2.00 3.0765 3.0620 0.0146
16+ 3.00 3.6453 3.6387 0.0066
17+ 5.00 4.3738 4.3714 0.0024
18+
19+Against the measured real-data gaps (001290 ch0, bandpassed):
20+ s_* predicted L measured gap
21+ 1.66 0.0209 0.0720
22+ 1.10 0.0450 0.0650
23+ 0.74 0.0895 0.1010
24+ 0.50 0.1602 0.1930
25+ 0.33 0.2415 0.4240
exploration/ephys/phase_loss.pyadded+57−0View file
@@ -0,0 +1,57 @@
1+"""The fractional-phase loss as a universal function of s alone.
2+
3+Given a real-valued prediction mu and prediction-error std s (both in
4+quantization steps), write d = mu - round(mu) for the fractional phase.
5+The conditional law of the integer residual r = round(y) - round(mu) is
6+
7+ P(r = k | d) = Phi((k + 1/2 - d)/s) - Phi((k - 1/2 - d)/s).
8+
9+A coder that knows mu pays the phase-averaged conditional entropy
10+ G(s) = E_d[ H(P(. | d)) ] (= theory.gauss_uniform_entropy)
11+A coder that codes the integer residual with one pooled histogram pays the
12+entropy of the phase-mixture
13+ M(s) = H( E_d[ P(. | d) ] ).
14+The difference L(s) = M(s) - G(s) >= 0 is a mutual information: what the
15+integer residual throws away about the phase. It depends on nothing but s.
16+"""
17+import numpy as np
18+from scipy.special import ndtr
19+
20+
21+def phase_entropies(s, kmax=None, nphase=2001):
22+ s = float(s)
23+ if s <= 0:
24+ return 0.0, 0.0, 0.0
25+ kmax = kmax or max(4, int(np.ceil(8 * s)) + 2)
26+ k = np.arange(-kmax, kmax + 1)[:, None]
27+ d = np.linspace(-0.5, 0.5, nphase)[None, :]
28+ p = ndtr((k + 0.5 - d) / s) - ndtr((k - 0.5 - d) / s)
29+ p = np.clip(p, 1e-300, None)
30+ p /= p.sum(axis=0, keepdims=True)
31+
32+ def H(q):
33+ q = np.clip(q, 1e-300, None)
34+ return -(q * np.log2(q)).sum(axis=0)
35+
36+ w = np.ones(nphase); w[0] = w[-1] = 0.5; w /= w.sum()
37+ cond = float(H(p) @ w) # G(s): knows the phase
38+ mix = float(H((p * w).sum(axis=1))) # M(s): pools over phases
39+ return mix, cond, mix - cond
40+
41+
42+def curve(s_values):
43+ return np.array([phase_entropies(s) for s in s_values])
44+
45+
46+if __name__ == "__main__":
47+ print(f"{'s':>7}{'M(s) mixture':>14}{'G(s) cond':>12}{'L(s) loss':>12}")
48+ for s in [0.1, 0.2, 0.28, 0.33, 0.4, 0.5, 0.6, 0.74, 0.9, 1.0, 1.1,
49+ 1.5, 1.66, 2.0, 3.0, 5.0]:
50+ m, g, l = phase_entropies(s)
51+ print(f"{s:7.2f}{m:14.4f}{g:12.4f}{l:12.4f}")
52+
53+ print("\nAgainst the measured real-data gaps (001290 ch0, bandpassed):")
54+ print(f" {'s_*':>6}{'predicted L':>13}{'measured gap':>14}")
55+ for s_star, measured in [(1.66, 0.072), (1.10, 0.065), (0.74, 0.101),
56+ (0.50, 0.193), (0.33, 0.424)]:
57+ print(f" {s_star:6.2f}{phase_entropies(s_star)[2]:13.4f}{measured:14.4f}")
exploration/ephys/sweep.logadded+42−0View file
@@ -0,0 +1,42 @@
1+
2+================================================================================================
3+ecephys-000409-ch101 n = 200000 (bandpass 300-6000 Hz, noise std = 1.0 before quantization)
4+================================================================================================
5+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
6+ real model surrogate
7+ ----------------------------------------------------------------------------------------------
8+ 0.125 8.00 1.359 2.5215 2.6132 2.5695 0.0437 2.6526 2.5921 0.0605 +0.0394 zstd+delta real 3.654 model 3.792
9+ 0.250 4.00 0.922 1.9968 2.1019 2.0479 0.0540 2.1486 2.0617 0.0870 +0.0467 zstd+delta real 2.740 model 2.852
10+ 0.500 2.00 0.634 1.5263 1.6774 1.5705 0.1069 1.7090 1.5809 0.1282 +0.0317 zstd+delta real 2.022 model 2.108
11+ 1.000 1.00 0.433 1.1047 1.4166 1.1202 0.2964 1.4085 1.1386 0.2699 -0.0081 zstd+delta real 1.435 model 1.536
12+ 2.000 0.50 0.283 0.7356 1.2075 0.6885 0.5190 1.2034 0.6834 0.5200 -0.0042 zstd+delta real 0.872 model 0.986
13+ 4.000 0.25 0.114 0.2969 0.6309 0.8837-0.2528 0.0220 0.0212 0.0008 -0.6089 zstd+delta real 0.242 model 0.013
14+ 8.000 0.12 0.001 0.0023 0.0870 0.1429-0.0558 0.0033 0.0057-0.0024 -0.0838 zstd+delta real 0.043 model 0.001
15+
16+================================================================================================
17+ecephys-000876-ch45 n = 200000 (bandpass 300-6000 Hz, noise std = 1.0 before quantization)
18+================================================================================================
19+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
20+ real model surrogate
21+ ----------------------------------------------------------------------------------------------
22+ 0.125 8.00 1.514 2.6716 2.7744 2.7039 0.0705 2.7943 2.7339 0.0604 +0.0199 zstd+delta real 4.063 model 4.099
23+ 0.250 4.00 1.023 2.1357 2.2628 2.1752 0.0876 2.2711 2.1909 0.0802 +0.0083 zstd+delta real 3.075 model 3.102
24+ 0.500 2.00 0.708 1.6596 1.8056 1.6929 0.1127 1.8218 1.7065 0.1153 +0.0162 zstd+delta real 2.261 model 2.329
25+ 1.000 1.00 0.487 1.2250 1.4693 1.2405 0.2288 1.4654 1.2625 0.2030 -0.0038 zstd+delta real 1.649 model 1.691
26+ 2.000 0.50 0.327 0.8480 1.2434 0.7768 0.4665 1.2458 0.8294 0.4163 +0.0024 zstd+delta real 1.103 model 1.165
27+ 4.000 0.25 0.156 0.4076 0.9376 2.1000-1.1624 0.4970 0.4829 0.0140 -0.4406 zstd+delta real 0.459 model 0.170
28+ 8.000 0.12 0.001 0.0024 0.2885 0.2309 0.0576 0.0033 0.0057-0.0024 -0.2852 zstd+delta real 0.085 model 0.001
29+
30+================================================================================================
31+ecephys-001290-ch0 n = 200000 (bandpass 300-6000 Hz, noise std = 1.0 before quantization)
32+================================================================================================
33+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
34+ real model surrogate
35+ ----------------------------------------------------------------------------------------------
36+ 0.125 8.00 1.658 2.7980 2.9302 2.8586 0.0716 2.9342 2.8806 0.0536 +0.0039 zstd+delta real 4.354 model 4.466
37+ 0.250 4.00 1.099 2.2318 2.3595 2.2947 0.0648 2.3884 2.3066 0.0818 +0.0289 zstd+delta real 3.293 model 3.308
38+ 0.500 2.00 0.741 1.7165 1.8719 1.7705 0.1014 1.8838 1.7826 0.1012 +0.0118 zstd+delta real 2.404 model 2.455
39+ 1.000 1.00 0.500 1.2545 1.4891 1.2963 0.1929 1.4930 1.3081 0.1848 +0.0038 zstd+delta real 1.745 model 1.780
40+ 2.000 0.50 0.325 0.8423 1.2603 0.8361 0.4242 1.2469 0.8424 0.4045 -0.0134 zstd+delta real 1.170 model 1.186
41+ 4.000 0.25 0.103 0.2674 0.9514 1.9378-0.9864 0.0082 0.0095-0.0012 -0.9431 zstd+delta real 0.430 model 0.006
42+ 8.000 0.12 0.001 0.0014 0.1018 0.0689 0.0329 0.0033 0.0057-0.0024 -0.0985 zstd+delta real 0.034 model 0.001
exploration/ephys/sweep.pyadded+105−0View file
@@ -0,0 +1,105 @@
1+"""Walk real ephys along the s_* axis by changing the quantization step.
2+
3+The filtered ephys variant that benchcompress benchmarks fixes the step at
4+v = 0.25 noise units, which puts the noise at 4 quantization steps — the
5+high-resolution corner where integer-residual coding is already near optimal.
6+Coarser steps move the same recording down the s_* axis into the regime where
7+the fractional-phase loss is supposed to bite.
8+
9+For each step size: requantize, fit the model, synthesize a surrogate, and
10+measure the prediction-based methods on both.
11+
12+Usage: python sweep.py [--n N] [--order O] [--zstd]
13+"""
14+import argparse
15+import glob
16+import os
17+import sys
18+
19+import numpy as np
20+from scipy.signal import butter, lfilter
21+
22+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
23+import codec_suite as cc # noqa: E402
24+from fitmodel import Fit # noqa: E402
25+
26+CACHE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "cache")
27+RATE = 30000.0
28+STEPS = [0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0, 32.0]
29+
30+
31+def bandpass(x, lowcut=300.0, highcut=6000.0, rate=RATE):
32+ nyq = 0.5 * rate
33+ b, a = butter(5, [lowcut / nyq, highcut / nyq], btype="band")
34+ return lfilter(b, a, x)
35+
36+
37+def highpass(x, lowcut=300.0, rate=RATE):
38+ nyq = 0.5 * rate
39+ b, a = butter(5, lowcut / nyq, btype="high")
40+ return lfilter(b, a, x)
41+
42+
43+def noise_units(raw, rate=RATE, do_bandpass=True):
44+ """Trace scaled so that the MAD noise level is 1.0, optionally bandpassed
45+ first. With do_bandpass=False nothing shapes the spectrum but the
46+ acquisition hardware — the test of whether the model needs an explicit
47+ filtering step or only a spectrum."""
48+ x = np.asarray(raw, dtype=np.float64) - np.median(raw)
49+ xf = bandpass(x, rate=rate) if do_bandpass else x
50+ ref = xf if do_bandpass else highpass(x, rate=rate)
51+ nl = float(np.median(np.abs(ref - np.median(ref))) / 0.6745)
52+ return xf / nl
53+
54+
55+def main():
56+ p = argparse.ArgumentParser()
57+ p.add_argument("--n", type=int, default=200_000)
58+ p.add_argument("--order", type=int, default=32)
59+ p.add_argument("--zstd", action="store_true",
60+ help="also measure zstd+delta as the practical baseline")
61+ p.add_argument("--nofilter", action="store_true",
62+ help="skip the bandpass: sweep the unfiltered trace")
63+ args = p.parse_args()
64+
65+ for path in sorted(glob.glob(os.path.join(CACHE, "*.raw.npy"))):
66+ name = os.path.basename(path)[: -len(".raw.npy")]
67+ y = noise_units(np.load(path), do_bandpass=not args.nofilter)[: args.n]
68+ shaping = ("unfiltered (acquisition spectrum only)" if args.nofilter
69+ else "bandpass 300-6000 Hz")
70+ print(f"\n{'=' * 96}\n{name} n = {y.size} "
71+ f"({shaping}, noise std = 1.0 before quantization)\n{'=' * 96}")
72+ print(f" {'step':>6}{'noise/step':>11}{'s_*':>8}{'G(s_*)':>9}"
73+ f"{'| LPC+ANS':>12}{'condG':>8}{'gap':>7}"
74+ f"{'| LPC+ANS':>12}{'condG':>8}{'gap':>7}{'| model-real':>13}")
75+ print(f" {' ' * 34}{'real':>27}{'model surrogate':>27}")
76+ print(f" {'-' * 94}")
77+ for v in STEPS:
78+ z = np.round(y / v).astype(np.int16)
79+ if np.unique(z).size < 3:
80+ continue
81+ fit = Fit(z, nfft=4096)
82+ zs = fit.synthesize(seed=0)
83+
84+ def block(zz):
85+ nb, h0 = cc.lpc_ans_bytes(zz, args.order)
86+ return (8.0 * nb / zz.size,
87+ cc.conditional_gaussian_rate(zz, args.order))
88+
89+ ra, rc = block(z)
90+ ma, mc = block(zs)
91+ extra = ""
92+ if args.zstd:
93+ nz = cc.GENERIC["zstd"](cc.delta(z).tobytes())
94+ ns = cc.GENERIC["zstd"](cc.delta(zs).tobytes())
95+ extra = (f" zstd+delta real {8.0 * nz / z.size:6.3f} "
96+ f"model {8.0 * ns / zs.size:6.3f}")
97+ print(f" {v:6.3f}{1.0 / v:11.2f}{fit.s_star:8.3f}"
98+ f"{fit.predicted_rate:9.4f}"
99+ f"{ra:12.4f}{rc:8.4f}{ra - rc:7.4f}"
100+ f"{ma:12.4f}{mc:8.4f}{ma - mc:7.4f}"
101+ f"{ma - ra:+13.4f}{extra}")
102+
103+
104+if __name__ == "__main__":
105+ main()
exploration/ephys/sweep_nofilter.logadded+44−0View file
@@ -0,0 +1,44 @@
1+
2+================================================================================================
3+ecephys-000409-ch101 n = 200000 (unfiltered (acquisition spectrum only), noise std = 1.0 before quantization)
4+================================================================================================
5+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
6+ real model surrogate
7+ ----------------------------------------------------------------------------------------------
8+ 0.125 8.00 4.350 4.1712 4.2325 4.1846 0.0479 4.2885 4.1978 0.0907 +0.0560
9+ 0.250 4.00 2.190 3.1904 3.2185 3.2050 0.0135 3.2768 3.2161 0.0607 +0.0583
10+ 0.500 2.00 1.112 2.2468 2.3228 2.2653 0.0575 2.3398 2.2713 0.0686 +0.0170
11+ 1.000 1.00 0.593 1.4458 1.6728 1.4722 0.2006 1.6703 1.4727 0.1976 -0.0025
12+ 2.000 0.50 0.319 0.8283 1.3068 0.8761 0.4307 1.2828 0.8060 0.4768 -0.0240
13+ 4.000 0.25 0.135 0.3519 0.7882 1.6277-0.8395 0.2355 0.1652 0.0704 -0.5527
14+ 8.000 0.12 0.001 0.0021 0.1821 0.1177 0.0644 0.0033 0.0057-0.0024 -0.1788
15+
16+================================================================================================
17+ecephys-000876-ch45 n = 200000 (unfiltered (acquisition spectrum only), noise std = 1.0 before quantization)
18+================================================================================================
19+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
20+ real model surrogate
21+ ----------------------------------------------------------------------------------------------
22+ 0.125 8.00 5.047 4.3849 4.4697 4.3968 0.0729 4.5026 4.4023 0.1004 +0.0330
23+ 0.250 4.00 2.539 3.4007 3.4474 3.4122 0.0352 3.4866 3.4185 0.0680 +0.0392
24+ 0.500 2.00 1.300 2.4606 2.5330 2.4724 0.0606 2.5622 2.4752 0.0871 +0.0292
25+ 1.000 1.00 0.696 1.6387 1.8077 1.6505 0.1571 1.8154 1.6549 0.1605 +0.0077
26+ 2.000 0.50 0.384 0.9876 1.3380 1.0121 0.3259 1.3541 0.9799 0.3742 +0.0161
27+ 4.000 0.25 0.234 0.6102 1.1036 0.6667 0.4369 1.1032 1.1548-0.0516 -0.0004
28+ 8.000 0.12 0.160 0.4171 1.0171 2.7210-1.7039 0.9432 1.8993-0.9561 -0.0739
29+ 16.000 0.06 0.076 0.1973 0.6390 1.0511-0.4120 0.1870 0.1058 0.0813 -0.4520
30+
31+================================================================================================
32+ecephys-001290-ch0 n = 200000 (unfiltered (acquisition spectrum only), noise std = 1.0 before quantization)
33+================================================================================================
34+ step noise/step s_* G(s_*) | LPC+ANS condG gap | LPC+ANS condG gap | model-real
35+ real model surrogate
36+ ----------------------------------------------------------------------------------------------
37+ 0.125 8.00 6.811 4.8163 4.9351 4.8613 0.0738 4.9411 4.8635 0.0776 +0.0060
38+ 0.250 4.00 3.404 3.8197 3.9443 3.8645 0.0798 3.9458 3.8665 0.0793 +0.0015
39+ 0.500 2.00 1.703 2.8357 2.9405 2.8797 0.0608 2.9326 2.8819 0.0507 -0.0079
40+ 1.000 1.00 0.843 1.8812 2.0226 1.9226 0.0999 2.0128 1.9238 0.0890 -0.0098
41+ 2.000 0.50 0.439 1.1177 1.4340 1.1556 0.2784 1.4080 1.1468 0.2612 -0.0260
42+ 4.000 0.25 0.139 0.3627 1.1339 2.7291-1.5952 0.3558 0.3274 0.0284 -0.7781
43+ 8.000 0.12 0.001 0.0018 0.1950 0.1620 0.0330 0.0033 0.0057-0.0024 -0.1918
44+ 16.000 0.06 0.000 0.0004 0.0190 0.0137 0.0053 0.0033 0.0057-0.0024 -0.0157