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).
79def quantized_gaussian_entropy(s):
80 if s <= 0.02:
81 return 0.0
82 zmax = int(math.ceil(8 * s + 4))
83 H = 0.0
84 prev = 0.5 * (1 + math.erf((-zmax - 0.5) / (s * math.sqrt(2))))
85 for z in range(-zmax, zmax + 1):
86 cur = 0.5 * (1 + math.erf((z + 0.5) / (s * math.sqrt(2))))
87 p = cur - prev
88 prev = cur
89 if p > 0:
90 H -= p * math.log2(p)
91 return H
94def formula_rate(kernel, sigma, dither, points=8192):
95 noise_var = 1 / 6 if dither else 1 / 12
96 f = (0.5 * (np.arange(points) + 0.5)) / points
97 w = -2j * np.pi * np.outer(f, np.arange(len(kernel)))
98 S = sigma * sigma * np.abs(np.exp(w) @ kernel) ** 2
99 integral = float(np.log2(S + noise_var).mean()) * 0.5
100 return quantized_gaussian_entropy(2.0 ** integral)
103# ------------------------------------------------ vectorized normal helpers
105def norm_pdf(t):
106 return np.exp(-0.5 * t * t) / SQRT2PI
109def norm_cdf(t):
110 """Abramowitz-Stegun 26.2.17; |error| < 7.5e-8, plenty for sampling and
111 for pmf weights whose entropy is wanted to ~1e-4 bits."""
112 t = np.asarray(t, dtype=float)
113 z = np.abs(t)
114 k = 1.0 / (1.0 + 0.2316419 * z)
115 poly = k * (0.319381530 + k * (-0.356563782 + k * (1.781477937 + k * (-1.821255978 + k * 1.330274429))))
116 tail = norm_pdf(z) * poly
117 return np.where(t >= 0, 1.0 - tail, tail)
120def norm_ppf(p):
121 """Acklam's rational approximation to the standard normal quantile."""
122 a = (-3.969683028665376e+01, 2.209460984245205e+02, -2.759285104469687e+02,
123 1.383577518672690e+02, -3.066479806614716e+01, 2.506628277459239e+00)
124 b = (-5.447609879822406e+01, 1.615858368580409e+02, -1.556989798598866e+02,
125 6.680131188771972e+01, -1.328068155288572e+01)
126 c = (-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e+00,
127 -2.549732539343734e+00, 4.374664141464968e+00, 2.938163982698783e+00)
128 d = (7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e+00,
129 3.754408661907416e+00)
130 p = np.asarray(p, dtype=float)
131 x = np.empty_like(p)
132 plow, phigh = 0.02425, 1 - 0.02425
134 lo = p < plow
135 hi = p > phigh
136 mid = ~(lo | hi)
138 if mid.any():
139 q = p[mid] - 0.5
140 r = q * q
141 x[mid] = (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q / \
142 (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
143 if lo.any():
144 q = np.sqrt(-2 * np.log(p[lo]))
145 x[lo] = (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
146 ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
147 if hi.any():
148 q = np.sqrt(-2 * np.log(1 - p[hi]))
149 x[hi] = -(((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5]) / \
150 ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
151 return x
154def trunc_std_normal(lo, hi, rng):
155 """Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into
156 the lower tail so the CDF differences keep precision."""
157 flip = (lo + hi) > 0
158 a = np.where(flip, -hi, lo)
159 b = np.where(flip, -lo, hi)
160 Fa = norm_cdf(a)
161 Fb = norm_cdf(b)
162 u = Fa + (Fb - Fa) * rng.random(a.shape)
163 x = norm_ppf(np.clip(u, 1e-300, 1 - 1e-16))
164 x = np.where(flip, -x, x)
165 return np.clip(x, lo, hi)
168# --------------------------------------------------------- the RB next-pmf
170def big_g(t):
171 """G(t) = t Phi(t) + phi(t), the antiderivative of Phi."""
172 return t * norm_cdf(t) + norm_pdf(t)
175def next_pmf(c, s0, dither):
176 """P(z_next = j | chain state): round(c + N(0, s0^2) (+ U(-1/2,1/2)))."""
177 s = max(s0, 1e-12)
178 half = 8 * s + (1.0 if dither else 0.0) + 1.0
179 js = np.arange(math.floor(c - half), math.ceil(c + half) + 1)
180 if dither:
181 # Integrating the Gaussian bin probability over the dither gives a
182 # second difference of G; as s -> 0 it degrades gracefully to the
183 # uniform-overlap width.
184 p = s * (big_g((js + 1 - c) / s) - 2 * big_g((js - c) / s) + big_g((js - 1 - c) / s))
185 else:
186 edges = norm_cdf((np.append(js, js[-1] + 1) - 0.5 - c) / s)
187 p = np.diff(edges)
188 return js, np.maximum(p, 0.0)
191# ---------------------------------------------------------------- one past
193def conditional_entropy_of_one_past(kernel, sigma, dither, M, sweeps, rng):
194 h = np.asarray(kernel, dtype=float)
195 L = len(h)
196 hr = h[::-1]
197 N = M + L - 1 # latents covering the windows of z_1..z_M
199 # The past, with its true latents as the (stationary) chain start.
200 x = sigma * rng.standard_normal(N)
201 y = np.convolve(x, h, mode='valid')
202 d = (rng.random(M) - 0.5) if dither else None
203 z = np.floor(y + (d if dither else 0.0) + 0.5)
205 # Boxes and y live in padded arrays so that every coordinate x_i sees
206 # exactly L constraint rows (rows outside the data are unconstrained).
207 P = L - 1
208 ypad = np.zeros(M + 2 * P)
209 lo = np.full(M + 2 * P, -np.inf)
210 hi = np.full(M + 2 * P, np.inf)
212 def set_boxes():
213 dd = d if dither else 0.0
214 lo[P:P + M] = z - 0.5 - dd
215 hi[P:P + M] = z + 0.5 - dd
217 set_boxes()
219 # Color classes: coordinates L apart share no constraint row, so a class
220 # updates as one vectorized block. Row i+j (padded) carries coefficient
221 # h[j] for coordinate i.
222 classes = [np.arange(c0, N, L) for c0 in range(L)]
223 rowmats = [idx[:, None] + np.arange(L)[None, :] for idx in classes]
224 nonzero = h != 0
226 s0 = sigma * abs(h[0])
227 pmf = {}
229 for _ in range(sweeps):
230 ypad[P:P + M] = np.convolve(x, h, mode='valid') # kill fp drift
231 if dither:
232 ycur = ypad[P:P + M]
233 dlo = np.maximum(-0.5, z - 0.5 - ycur)
234 dhi = np.minimum(0.5, z + 0.5 - ycur)
235 d = dlo + np.maximum(dhi - dlo, 0.0) * rng.random(M)
236 set_boxes()
237 for idx, rows in zip(classes, rowmats):
238 r = ypad[rows] - np.outer(x[idx], h)
239 with np.errstate(divide='ignore', invalid='ignore'):
240 b1 = (lo[rows] - r) / h[None, :]
241 b2 = (hi[rows] - r) / h[None, :]
242 xlo = np.where(h[None, :] > 0, b1, b2)
243 xhi = np.where(h[None, :] > 0, b2, b1)
244 xlo[:, ~nonzero] = -np.inf
245 xhi[:, ~nonzero] = np.inf
246 xlo = xlo.max(axis=1)
247 xhi = xhi.min(axis=1)
248 xnew = trunc_std_normal(xlo / sigma, xhi / sigma, rng) * sigma
249 delta = xnew - x[idx]
250 x[idx] = xnew
251 ypad[rows] += delta[:, None] * h[None, :]
253 c = float(hr[:-1] @ x[M:]) if L > 1 else 0.0
254 js, p = next_pmf(c, s0, dither)
255 for j, pj in zip(js, p):
256 if pj > 1e-15:
257 pmf[int(j)] = pmf.get(int(j), 0.0) + pj
259 total = sum(pmf.values())
260 return -sum((p / total) * math.log2(p / total) for p in pmf.values() if p > 0)
263# --------------------------------------------------------------------- main
265def main():
266 ap = argparse.ArgumentParser(
267 description='Monte-Carlo estimate of the true entropy rate, for '
268 'checking the analytic rate R shown in the app.')
269 ap.add_argument('--sigma', type=float, required=True, help='input std, in quantization steps')
270 ap.add_argument('--filter', required=True,
271 choices=['none', 'moving-average', 'lowpass', 'bandpass', 'first-difference'])
272 ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
273 ap.add_argument('--high', type=float, help='lowpass cutoff / bandpass high edge, Hz')
274 ap.add_argument('--taps', type=int, default=101, help='windowed-sinc kernel length')
275 ap.add_argument('--width', type=int, default=8, help='moving-average width')
276 ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
277 ap.add_argument('--dither', action='store_true')
278 ap.add_argument('--past', type=int, help='conditioning window M (default max(512, 4·taps))')
279 ap.add_argument('--pasts', type=int, default=24, help='independent pasts to average')
280 ap.add_argument('--sweeps', type=int, default=600, help='Gibbs sweeps per past')
281 ap.add_argument('--seed', type=int, default=0)
282 args = ap.parse_args()
284 if args.filter == 'bandpass' and (args.low is None or args.high is None):
285 ap.error('bandpass needs --low and --high')
286 if args.filter == 'lowpass' and args.high is None:
287 ap.error('lowpass needs --high')
289 kernel = design_kernel(args)
290 L = len(kernel)
291 M = args.past if args.past is not None else max(512, 4 * L)
293 R = formula_rate(kernel, args.sigma, args.dither)
294 print(f'model: sigma={args.sigma} filter={args.filter} L={L} dither={args.dither}')
295 print(f'formula R (as shown in the app): {R:.4f} bits/sample'
296 f' (ratio vs int16: {16 / R:.3f}x)' if R > 0 else f'formula R: {R:.4f} bits/sample')
297 print(f'MC: {args.pasts} pasts x {args.sweeps} sweeps, conditioning on M={M} samples')
299 rng = np.random.default_rng(args.seed)
300 Hs = []
301 for i in range(args.pasts):
302 Hs.append(conditional_entropy_of_one_past(
303 kernel, args.sigma, args.dither, M, args.sweeps, rng))
304 mean = float(np.mean(Hs))
305 se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs))) if len(Hs) > 1 else float('nan')
306 print(f' past {i + 1:3d}/{args.pasts}: H = {Hs[-1]:.4f} running mean {mean:.4f} +/- {se:.4f}')
308 mean = float(np.mean(Hs))
309 se = float(np.std(Hs, ddof=1) / math.sqrt(len(Hs)))
310 print(f'\nMC entropy rate: {mean:.4f} +/- {se:.4f} bits/sample'
311 f' (ratio vs int16: {16 / mean:.3f}x)')
312 print(f'formula R: {R:.4f} bits/sample (formula - MC = {R - mean:+.4f})')
313 print('note: the MC value estimates H(z_next | M past samples), an upper bound '
314 'on the rate that tightens as --past grows.')
317if __name__ == '__main__':
318 main()