1#!/usr/bin/env python3
2"""Monte-Carlo ground truth for the theoretical rate R shown in the app.
4Model (matching the app): x iid N(0, sigma^2) -> y = h * x -> optional
5additive uniform dither on [-1/2, 1/2) -> z = round(.). The entropy rate
6H = E[-log2 P(z_next | past)] is the true lossless limit in bits/sample.
7(The app applies the kernel zero-phase; a time shift does not change the
8law of the process, so causal convolution is used here.)
10Method
11 1. Draw a past: sample x (and dither) from the prior and push it through
12 the pipeline to get z_1..z_M.
13 2. The generating x is itself an exact draw from p(x | z_1..z_M), so a
14 Gibbs chain started there is already in stationarity - no burn-in
15 bias, only autocorrelation.
16 3. Gibbs-sample x | z: this posterior is a box-truncated multivariate
17 normal, and each conditional x_i | rest is N(0, sigma^2) truncated to
18 an interval read off the <= L constraint boxes x_i appears in. With
19 dither on, the dither values are extra latents with uniform
20 conditionals. Coordinates a multiple of L apart share no constraint,
21 so each of the L "colors" is updated as one vectorized block.
22 4. Rao-Blackwellization: given a chain state, the next sample is
23 z_next = round(c + h_0 * x_free (+ d)) with x_free ~ N(0, sigma^2)
24 still unconstrained, so P(z_next = j | state) has a closed form.
25 Averaging these pmfs over the chain gives P(z_next | past) exactly in
26 the limit; its entropy is the conditional entropy for that past.
27 5. Average over independent pasts; report mean +/- standard error.
29The estimate targets H(z_{M+1} | z_1..z_M), which is an upper bound on the
30rate and decreases toward it as --past grows beyond the memory of the
31process. More --sweeps reduces the (downward) plug-in bias from noise in
32the averaged pmf; narrowband filters and large sigma mix more slowly and
33deserve more sweeps.
35Requires numpy only.
36"""
38import argparse
39import math
41import numpy as np
43SQRT2PI = math.sqrt(2 * math.pi)
46# ---------------------------------------------------------------- kernels
47# Ported from src/model/filters.ts; must stay in step with it.
49def windowed_sinc_lowpass(fc, taps):
50 n = taps | 1
51 mid = (n - 1) / 2
52 i = np.arange(n)
53 t = i - mid
54 sinc = np.where(t == 0, 2 * fc, np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
55 w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
56 h = sinc * w
57 return h / h.sum()
60def design_kernel(args):
61 if args.filter == 'none':
62 return np.array([1.0])
63 if args.filter == 'moving-average':
64 return np.full(args.width, 1.0 / args.width)
65 if args.filter == 'lowpass':
66 return windowed_sinc_lowpass(args.high / args.rate, args.taps)
67 if args.filter == 'bandpass':
68 lo = windowed_sinc_lowpass(args.low / args.rate, args.taps)
69 hi = windowed_sinc_lowpass(args.high / args.rate, args.taps)
70 return hi - lo
71 if args.filter == 'first-difference':
72 return np.array([1.0, -1.0])
73 raise ValueError(args.filter)
76# ------------------------------------------------- the app's formula for R
77# Ported from src/model/theory.ts. Phi via math.erf (machine precision).
79TWO_PI_E = 2 * math.pi * math.e
82def quantized_gaussian_entropy(s):
83 if s <= 0.02:
84 return 0.0
85 # The discrete entropy approaches 1/2 log2(2 pi e s^2) from above like
86 # log2(e)/(24 s^2) (the delta^2/24 Fisher-information correction); at
87 # s >= 6 the corrected asymptote is within 2e-6 bits.
88 if s >= 6:
89 return 0.5 * math.log2(TWO_PI_E * s * s) + math.log2(math.e) / (24 * s * s)
90 zmax = int(math.ceil(8 * s + 4))
91 H = 0.0
92 prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
93 for z in range(-zmax, zmax + 1):
94 cur = 0.5 * (1 + math.erf((z + 0.5) / (s * math.sqrt(2))))
95 p = cur - prev
96 prev = cur
97 if p > 0:
98 H -= p * math.log2(p)
99 return H
102def dithered_quantized_gaussian_entropy(s):
103 """Exact entropy of round(N(0, s^2) + U[-1/2, 1/2)) — the marginal of a
104 stored sample with dither on. The pmf has the closed form
105 p_j = s * (G((j+1)/s) - 2 G(j/s) + G((j-1)/s)) with G(t) = t Phi(t) + phi(t)
106 the antiderivative of Phi; machine-exact at every s via math.erf."""
107 if s <= 0:
108 return 0.0
110 def G(t):
111 return t * 0.5 * (1 + math.erf(t / math.sqrt(2))) + math.exp(-0.5 * t * t) / SQRT2PI
113 jmax = int(math.ceil(8 * s + 2))
114 H = 0.0
115 for j in range(-jmax, jmax + 1):
116 p = s * (G((j + 1) / s) - 2 * G(j / s) + G((j - 1) / s))
117 if p > 0:
118 H -= p * math.log2(p)
119 return H
122def formula_rates(kernel, sigma, dither, points=8192):
123 """The app's R = min(Rspec, Rsamp). Rspec charges the rounding(+dither)
124 noise at its full variance nu per Fourier mode; Rsamp is the exact
125 marginal entropy of one stored sample, a subadditivity upper bound that
126 takes over when the whole process is sub-threshold. Midpoint grid on
127 [0, 1/2]; |H| is symmetric, so the grid mean equals the unit-circle
128 integral. Returns (rspec, rsamp)."""
129 nu = 1 / 6 if dither else 1 / 12
130 f = (0.5 * (np.arange(points) + 0.5)) / points
131 w = -2j * np.pi * np.outer(f, np.arange(len(kernel)))
132 S = sigma * sigma * np.abs(np.exp(w) @ kernel) ** 2
133 rspec = float(np.mean(0.5 * np.log2(TWO_PI_E * (S + nu))))
134 v = sigma * sigma * float(np.sum(np.asarray(kernel) ** 2))
135 rsamp = (dithered_quantized_gaussian_entropy(math.sqrt(v)) if dither
136 else quantized_gaussian_entropy(math.sqrt(v)))
137 return rspec, rsamp
140# ------------------------------------------------ vectorized normal helpers
142def norm_pdf(t):
143 return np.exp(-0.5 * t * t) / SQRT2PI
146def norm_cdf(t):
147 """Abramowitz-Stegun 26.2.17; |error| < 7.5e-8, plenty for sampling and
148 for pmf weights whose entropy is wanted to ~1e-4 bits."""
149 t = np.asarray(t, dtype=float)
150 z = np.abs(t)
151 k = 1.0 / (1.0 + 0.2316419 * z)
152 poly = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + k * 1.330274429))))
153 tail = norm_pdf(z) * poly
154 return np.where(t >= 0, 1.0 - tail, tail)
157def norm_ppf(p):
158 """Acklam's rational approximation to the standard normal quantile."""
159 a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
160 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00)
161 b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
162 6.680131188771972e+01, -1.328068155288572e+01)
163 c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
164 -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00)
165 d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
166 3.754408661907416e+00)
167 p = np.asarray(p, dtype=float)
168 x = np.empty_like(p)
169 plow, phigh = 0.02425, 1 - 0.02425
171 lo = p < plow
172 hi = p > phigh
173 mid = ~(lo | hi)
175 if mid.any():
176 q = p[mid] - 0.5
177 r = q * q
178 x[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \
179 (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
180 if lo.any():
181 q = np.sqrt(-2 * np.log(p[lo]))
182 x[lo] = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
183 ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
184 if hi.any():
185 q = np.sqrt(-2 * np.log(1 - p[hi]))
186 x[hi] = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
187 ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
188 return x
191def trunc_std_normal(lo, hi, rng):
192 """Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into
193 the lower tail so the CDF differences keep precision."""
194 flip = (lo + hi) > 0
195 a = np.where(flip, -hi, lo)
196 b = np.where(flip, -lo, hi)
197 Fa = norm_cdf(a)
198 Fb = norm_cdf(b)
199 u = Fa + (Fb - Fa) * rng.random(a.shape)
200 x = norm_ppf(np.clip(u, 1e-300, 1 - 1e-16))
201 x = np.where(flip, -x, x)
202 return np.clip(x, lo, hi)
205# --------------------------------------------------------- the RB next-pmf
207def big_g(t):
208 """G(t) = t Phi(t) + phi(t), the antiderivative of Phi."""
209 return t * norm_cdf(t) + norm_pdf(t)
212def next_pmf(c, s0, dither):
213 """P(z_next = j | chain state): round(c + N(0, s0^2) (+ U(-1/2,1/2)))."""
214 s = max(s0, 1e-12)
215 half = 8 * s + (1.0 if dither else 0.0) + 1.0
216 js = np.arange(math.floor(c - half), math.ceil(c + half) + 1)
217 if dither:
218 # Integrating the Gaussian bin probability over the dither gives a
219 # second difference of G; as s -> 0 it degrades gracefully to the
220 # uniform-overlap width.
221 p = s * (big_g((js + 1 - c) / s) - 2 * big_g((js - c) / s) + big_g((js - 1 - c) / s))
222 else:
223 edges = norm_cdf((np.append(js, js[-1] + 1) - 0.5 - c) / s)
224 p = np.diff(edges)
225 return js, np.maximum(p, 0.0)
228# ---------------------------------------------------------------- one past
230def conditional_entropy_of_one_past(kernel, sigma, dither, M, sweeps, rng):
231 h = np.asarray(kernel, dtype=float)
232 L = len(h)
233 hr = h[::-1]
234 N = M + L - 1 # latents covering the windows of z_1..z_M
236 # The past, with its true latents as the (stationary) chain start.
237 x = sigma * rng.standard_normal(N)
238 y = np.convolve(x, h, mode='valid')
239 d = (rng.random(M) - 0.5) if dither else None
240 z = np.floor(y + (d if dither else 0.0) + 0.5)
242 # Boxes and y live in padded arrays so that every coordinate x_i sees
243 # exactly L constraint rows (rows outside the data are unconstrained).
244 P = L - 1
245 ypad = np.zeros(M + 2 * P)
246 lo = np.full(M + 2 * P, -np.inf)
247 hi = np.full(M + 2 * P, np.inf)
249 def set_boxes():
250 dd = d if dither else 0.0
251 lo[P:P + M] = z - 0.5 - dd
252 hi[P:P + M] = z + 0.5 - dd
254 set_boxes()
256 # Color classes: coordinates L apart share no constraint row, so a class
257 # updates as one vectorized block. Row i+j (padded) carries coefficient
258 # h[j] for coordinate i.
259 classes = [np.arange(c0, N, L) for c0 in range(L)]
260 rowmats = [idx[:, None] + np.arange(L)[None, :] for idx in classes]
261 nonzero = h != 0
263 s0 = sigma * abs(h[0])
264 pmf = {}
266 for _ in range(sweeps):
267 ypad[P:P + M] = np.convolve(x, h, mode='valid') # kill fp drift
268 if dither:
269 ycur = ypad[P:P + M]
270 dlo = np.maximum(-0.5, z - 0.5 - ycur)
271 dhi = np.minimum(0.5, z + 0.5 - ycur)
272 d = dlo + np.maximum(dhi - dlo, 0.0) * rng.random(M)
273 set_boxes()
274 for idx, rows in zip(classes, rowmats):
275 r = ypad[rows] - np.outer(x[idx], h)
276 with np.errstate(divide='ignore', invalid='ignore'):
277 b1 = (lo[rows] - r) / h[None, :]
278 b2 = (hi[rows] - r) / h[None, :]
279 xlo = np.where(h[None, :] > 0, b1, b2)
280 xhi = np.where(h[None, :] > 0, b2, b1)
281 xlo[:, ~nonzero] = -np.inf
282 xhi[:, ~nonzero] = np.inf
283 xlo = xlo.max(axis=1)
284 xhi = xhi.min(axis=1)
285 xnew = trunc_std_normal(xlo / sigma, xhi / sigma, rng) * sigma
286 delta = xnew - x[idx]
287 x[idx] = xnew
288 ypad[rows] += delta[:, None] * h[None, :]
290 c = float(hr[:-1] @ x[M:]) if L > 1 else 0.0
291 js, p = next_pmf(c, s0, dither)
292 for j, pj in zip(js, p):
293 if pj > 1e-15:
294 pmf[int(j)] = pmf.get(int(j), 0.0) + pj
296 total = sum(pmf.values())
297 return -sum((p / total) * math.log2(p / total) for p in pmf.values() if p > 0)
300# --------------------------------------------------------------------- main
302def main():
303 ap = argparse.ArgumentParser(
304 description='Monte-Carlo estimate of the true entropy rate, for '
305 'checking the analytic rate R shown in the app.')
306 ap.add_argument('--sigma', type=float, required=True, help='input std, in quantization steps')
307 ap.add_argument('--filter', required=True,
308 choices=['none', 'moving-average', 'lowpass', 'bandpass', 'first-difference'])
309 ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
310 ap.add_argument('--high', type=float, help='lowpass cutoff / bandpass high edge, Hz')
311 ap.add_argument('--taps', type=int, default=101, help='windowed-sinc kernel length')
312 ap.add_argument('--width', type=int, default=8, help='moving-average width')
313 ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
314 ap.add_argument('--dither', action='store_true')
315 ap.add_argument('--past', type=int, help='conditioning window M (default max(512, 4·taps))')
316 ap.add_argument('--pasts', type=int, default=24, help='independent pasts to average')
317 ap.add_argument('--sweeps', type=int, default=600, help='Gibbs sweeps per past')
318 ap.add_argument('--seed', type=int, default=0)
319 args = ap.parse_args()
321 if args.filter == 'bandpass' and (args.low is None or args.high is None):
322 ap.error('bandpass needs --low and --high')
323 if args.filter == 'lowpass' and args.high is None:
324 ap.error('lowpass needs --high')
326 kernel = design_kernel(args)
327 L = len(kernel)
328 M = args.past if args.past is not None else max(512, 4 * L)
330 R_spec, R_samp = formula_rates(kernel, args.sigma, args.dither)
331 R = min(R_spec, R_samp)
332 print(f'model: sigma={args.sigma} filter={args.filter} L={L} dither={args.dither}')
333 ratio = f' (ratio vs int16: {16 / R:.3f}x)' if R > 0 else ''
334 print(f'formula R = min(spec {R_spec:.4f}, samp {R_samp:.4f}) = {R:.4f} bits/sample{ratio}')
335 print(f'MC: {args.pasts} pasts x {args.sweeps} sweeps, conditioning on M={M} samples')
337 rng = np.random.default_rng(args.seed)
338 Hs = []
339 for i in range(args.pasts):
340 Hs.append(conditional_entropy_of_one_past(
341 kernel, args.sigma, args.dither, M, args.sweeps, rng))
342 mean = float(np.mean(Hs))
343 se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs))) if len(Hs) > 1 else float('nan')
344 print(f' past {i + 1:3d}/{args.pasts}: H = {Hs[-1]:.4f} running mean {mean:.4f} +/- {se:.4f}')
346 mean = float(np.mean(Hs))
347 se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs)))
348 print(f'\nMC entropy rate: {mean:.4f} +/- {se:.4f} bits/sample'
349 f' (ratio vs int16: {16 / mean:.3f}x)')
350 print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f}; '
351 f'spec {R_spec:.4f}, samp {R_samp:.4f})')
352 print('note: the MC value estimates H(z_next | M past samples), an upper bound '
353 'on the rate that tightens as --past grows.')
356if __name__ == '__main__':
357 main()