1"""Achievable rate of a model-based Monte-Carlo codec (the nonlinear optimum).
3The Gaussian-conditional codec (codec_gaussian.py) is linear: it conditions on
4the past only through a linear prediction and charges a fixed Gaussian. The
5true conditional law P(z_t | z_{t-M..t-1}) is the posterior predictive of the
6generative model — computable by Gibbs-sampling the latent x under the box
7constraints round(h*x) = z (exactly the machinery of timeseries-entropy) and
8Rao-Blackwellizing:
10 P_hat(k) = mean over sweeps of Phi((k+1/2-c)/(sigma h0)) - Phi((k-1/2-c)/(sigma h0)),
11 c = hr_head . x[M:]
13An encoder/decoder pair sharing the RNG seed can both compute P_hat from
14already-decoded samples, so rate = E[-log2 P_tilde(z_t)] is achievable
15(P_tilde mixes in a wide floor distribution so no symbol gets probability 0).
16The -log2 of a finite-sweep average is on average >= the -log2 of the true
17predictive (Jensen), so this measures an upper bound that tightens with more
18sweeps: whatever number comes out IS achievable by this codec family.
20Chain start: the generating latents of the window, an exact draw from
21p(x | all z) — near the target p(x | window z); BURN sweeps then wash the
22difference. Reuses ConditionalChain's sweep by overwriting its state.
23"""
24import numpy as np
25from scipy.special import ndtr
27import sys
28from pathlib import Path
30# The companion package, cloned alongside this repo.
31sys.path.insert(0, str(Path(__file__).resolve().parents[2] / 'timeseries-entropy' / 'src'))
32from timeseries_entropy.model import ConditionalChain
34SIGMA = 5.0
35RATE = 30000.0
36LOW, HIGH, TAPS = 300.0, 2000.0, 31
39def windowed_sinc_lowpass(fc, taps):
40 n = taps | 1
41 mid = (n - 1) / 2
42 i = np.arange(n)
43 t = i - mid
44 sinc = np.where(t == 0, 2 * fc,
45 np.sin(2 * np.pi * fc * t) / (np.pi * np.where(t == 0, 1, t)))
46 w = 0.54 - 0.46 * np.cos(2 * np.pi * i / (n - 1))
47 h = sinc * w
48 return h / h.sum()
51def make_kernel():
52 return windowed_sinc_lowpass(HIGH / RATE, TAPS) - windowed_sinc_lowpass(LOW / RATE, TAPS)
55def posterior_predictive_rate(M=512, burn=30, sweeps=60, positions=300,
56 spacing=64, seed=7, eps=1e-3):
57 h = make_kernel()
58 L = len(h)
59 rng = np.random.default_rng(seed)
61 n = positions * spacing + M + 4 * L
62 x = SIGMA * rng.standard_normal(n + L - 1)
63 y = np.convolve(x, h, mode='valid')
64 z = np.floor(y + 0.5)
66 h0 = abs(h[0])
67 print(f'sigma*|h0| (innovation scale given full latent past) = {SIGMA * h0:.5f}')
69 # wide fallback so P_tilde is never 0 (what a real codec would also do)
70 std_z = z.std()
71 kmax = 32
72 ks = np.arange(-kmax, kmax + 1)
73 fallback = ndtr((ks + 0.5) / std_z) - ndtr((ks - 0.5) / std_z)
74 fallback /= fallback.sum()
76 proto = ConditionalChain(h, SIGMA, past=M, rng=rng) # template; state overwritten
77 rates = []
78 hit_floor = 0
79 for i in range(positions):
80 t = M + i * spacing # predict z[t] from z[t-M:t]
81 i0 = t - M
82 proto.z = z[i0:t].copy()
83 proto.x = x[i0:t + L - 1].copy() # exact posterior draw (given all z)
84 P = L - 1
85 proto.lo[P:P + M] = proto.z - 0.5
86 proto.hi[P:P + M] = proto.z + 0.5
88 for _ in range(burn):
89 proto._sweep()
90 acc = 0.0
91 zt = int(z[t])
92 for _ in range(sweeps):
93 proto._sweep()
94 c = float(proto.hr_head @ proto.x[M:])
95 acc += (ndtr((zt + 0.5 - c) / (SIGMA * h0))
96 - ndtr((zt - 0.5 - c) / (SIGMA * h0)))
97 p_hat = acc / sweeps
98 p_tilde = (1 - eps) * p_hat + eps * float(fallback[np.clip(zt, -kmax, kmax) + kmax])
99 if p_hat < eps:
100 hit_floor += 1
101 rates.append(-np.log2(p_tilde))
102 if (i + 1) % 50 == 0:
103 r = np.array(rates)
104 print(f' {i + 1:4d}/{positions}: rate = {r.mean():.4f} '
105 f'+/- {r.std(ddof=1) / np.sqrt(len(r)):.4f} bits/sample '
106 f'(floor hits: {hit_floor})')
107 r = np.array(rates)
108 print(f'\nMC-codec achievable rate (M={M}, {sweeps} sweeps): '
109 f'{r.mean():.4f} +/- {r.std(ddof=1) / np.sqrt(len(r)):.4f} bits/sample '
110 f'-> ratio {16 / r.mean():.1f}x')
111 return r
114if __name__ == '__main__':
115 posterior_predictive_rate()