8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 1"""Command-line interface: timeseries-entropy [options]."""
3import argparse
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 4import os
5from concurrent.futures import ProcessPoolExecutor
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 6
7import numpy as np
9from . import estimate_conditional_entropy, level_corrections, kernels
10from .model import ConditionalChain
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 11from .theory import predict_entropy_rate
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 12
14def design_kernel(args):
15 if args.filter == 'none':
16 return kernels.identity()
17 if args.filter == 'moving-average':
18 return kernels.moving_average(args.width)
19 if args.filter == 'lowpass':
20 if args.high is None:
21 raise SystemExit('lowpass needs --high')
22 return kernels.windowed_sinc_lowpass(args.high / args.rate, args.taps)
23 if args.filter == 'bandpass':
24 if args.low is None or args.high is None:
25 raise SystemExit('bandpass needs --low and --high')
26 return kernels.windowed_sinc_bandpass(
27 args.low / args.rate, args.high / args.rate, args.taps)
28 if args.filter == 'first-difference':
29 return kernels.first_difference()
30 raise ValueError(args.filter)
33def main():
34 ap = argparse.ArgumentParser(
35 prog='timeseries-entropy',
36 description='Unbiased Monte-Carlo estimate of H(z_next | M past '
37 'samples), in bits, for x iid N(0, sigma^2) -> h * x '
38 '-> round.')
39 ap.add_argument('--sigma', type=float, required=True,
40 help='input std, in quantization steps')
41 ap.add_argument('--filter', required=True,
42 choices=['none', 'moving-average', 'lowpass', 'bandpass',
43 'first-difference'])
44 ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
45 ap.add_argument('--high', type=float,
46 help='lowpass cutoff / bandpass high edge, Hz')
47 ap.add_argument('--taps', type=int, default=101,
48 help='windowed-sinc kernel length')
49 ap.add_argument('--width', type=int, default=8, help='moving-average width')
50 ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
51 ap.add_argument('--past', type=int,
52 help='conditioning window M (default max(512, 4*L))')
53 ap.add_argument('--pasts', type=int, default=24,
54 help='independent pasts to average')
55 ap.add_argument('--reps', type=int, default=8,
56 help='randomized realizations per past')
57 ap.add_argument('--n0', type=int, default=128, help='base block size')
58 ap.add_argument('--r', type=float, default=1.5,
59 help='truncation exponent: P(N >= m) = 2^(-r m)')
60 ap.add_argument('--thin', type=int, default=1,
61 help='Gibbs sweeps per emitted sample')
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 62 ap.add_argument('--workers', type=int,
63 help='parallel processes over pasts (default: all cores)')
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 64 ap.add_argument('--seed', type=int, default=0)
65 ap.add_argument('--pilot', type=int, metavar='LEVELS',
66 help='instead of estimating, print RMS Delta_m over the '
67 'pasts for m = 1..LEVELS, to help choose --r')
68 args = ap.parse_args()
70 kernel = design_kernel(args)
71 L = len(kernel)
72 M = args.past if args.past is not None else max(512, 4 * L)
73 print(f'model: sigma={args.sigma} filter={args.filter} L={L} M={M}')
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 74 pred = predict_entropy_rate(kernel, args.sigma)
75 print(f'predicted rate: {pred["corrected"]:.4f} bits '
76 f'(quantization-corrected; s*={pred["s_star"]:.4g}) '
77 f'high-res Szego: {pred["highres"]:.4f} '
78 f'(sigma_inf={pred["sigma_inf"]:.4g})')
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 79
80 if args.pilot is not None:
81 run_pilot(kernel, args, M)
82 return
84 print(f'{args.pasts} pasts x {args.reps} reps, n0={args.n0} r={args.r} '
85 f'thin={args.thin}')
87 def progress(i, values):
88 mean = float(np.mean(values))
89 se = (float(np.std(values, ddof=1) / np.sqrt(len(values)))
90 if len(values) > 1 else float('nan'))
91 print(f' past {i + 1:3d}/{args.pasts}: H = {values[-1]:.4f} '
92 f'running mean {mean:.4f} +/- {se:.4f}')
94 est = estimate_conditional_entropy(
95 kernel, args.sigma, past=M, pasts=args.pasts, reps=args.reps,
96 n0=args.n0, r=args.r, thin=args.thin, seed=args.seed,
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 97 progress=progress, workers=args.workers)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 98 ratio = f' (ratio vs int16: {16 / est.mean:.3f}x)' if est.mean > 0 else ''
99 print(f'\nH(z_next | {M} past samples) = {est.mean:.4f} +/- {est.se:.4f} '
100 f'bits/sample{ratio}')
101 print('note: an upper bound on the entropy rate that tightens as --past '
102 'grows.')
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 105def _one_pilot(kernel, sigma, M, thin, n0, levels, seed_seq):
106 rng = np.random.default_rng(seed_seq)
107 return level_corrections(
108 ConditionalChain(kernel, sigma, M, rng, thin).draw, n0, levels)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 111def run_pilot(kernel, args, M):
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 112 seeds = np.random.SeedSequence(args.seed).spawn(args.pasts)
113 workers = args.workers or min(args.pasts, os.cpu_count() or 1)
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 114 print(f'pilot: {args.pasts} pasts, levels 1..{args.pilot}, n0={args.n0}')
4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 115 with ProcessPoolExecutor(max_workers=workers) as pool:
116 deltas = np.array(list(pool.map(
117 _one_pilot,
118 *zip(*[(kernel, args.sigma, M, args.thin, args.n0, args.pilot, s)
119 for s in seeds]))))
8c9a374Unbiased Monte-Carlo entropy estimation for quantized filtered Gaussian seriesJeremy Magland 120 rms = np.sqrt((deltas ** 2).mean(axis=0))
121 for m in range(args.pilot):
122 note = ''
123 if m > 0 and rms[m] > 0:
124 note = f' decay exponent {np.log2(rms[m - 1] / rms[m]) * 2:.2f}'
125 print(f' m={m + 1}: rms Delta = {rms[m]:.5f}{note}')
126 print('choose r safely below the E[Delta^2] decay exponent (and > 1); '
127 'r=1.5 suits decay near 2.')
130if __name__ == '__main__':
131 main()