1"""Command-line interface: timeseries-entropy [options]."""
3import argparse
5import numpy as np
7from . import estimate_conditional_entropy, level_corrections, kernels
8from .model import ConditionalChain
11def design_kernel(args):
12 if args.filter == 'none':
13 return kernels.identity()
14 if args.filter == 'moving-average':
15 return kernels.moving_average(args.width)
16 if args.filter == 'lowpass':
17 if args.high is None:
18 raise SystemExit('lowpass needs --high')
19 return kernels.windowed_sinc_lowpass(args.high / args.rate, args.taps)
20 if args.filter == 'bandpass':
21 if args.low is None or args.high is None:
22 raise SystemExit('bandpass needs --low and --high')
23 return kernels.windowed_sinc_bandpass(
24 args.low / args.rate, args.high / args.rate, args.taps)
25 if args.filter == 'first-difference':
26 return kernels.first_difference()
27 raise ValueError(args.filter)
30def main():
31 ap = argparse.ArgumentParser(
32 prog='timeseries-entropy',
33 description='Unbiased Monte-Carlo estimate of H(z_next | M past '
34 'samples), in bits, for x iid N(0, sigma^2) -> h * x '
35 '-> round.')
36 ap.add_argument('--sigma', type=float, required=True,
37 help='input std, in quantization steps')
38 ap.add_argument('--filter', required=True,
39 choices=['none', 'moving-average', 'lowpass', 'bandpass',
40 'first-difference'])
41 ap.add_argument('--low', type=float, help='bandpass low edge, Hz')
42 ap.add_argument('--high', type=float,
43 help='lowpass cutoff / bandpass high edge, Hz')
44 ap.add_argument('--taps', type=int, default=101,
45 help='windowed-sinc kernel length')
46 ap.add_argument('--width', type=int, default=8, help='moving-average width')
47 ap.add_argument('--rate', type=float, default=30000, help='sample rate, Hz')
48 ap.add_argument('--past', type=int,
49 help='conditioning window M (default max(512, 4*L))')
50 ap.add_argument('--pasts', type=int, default=24,
51 help='independent pasts to average')
52 ap.add_argument('--reps', type=int, default=8,
53 help='randomized realizations per past')
54 ap.add_argument('--n0', type=int, default=128, help='base block size')
55 ap.add_argument('--r', type=float, default=1.5,
56 help='truncation exponent: P(N >= m) = 2^(-r m)')
57 ap.add_argument('--thin', type=int, default=1,
58 help='Gibbs sweeps per emitted sample')
59 ap.add_argument('--seed', type=int, default=0)
60 ap.add_argument('--pilot', type=int, metavar='LEVELS',
61 help='instead of estimating, print RMS Delta_m over the '
62 'pasts for m = 1..LEVELS, to help choose --r')
63 args = ap.parse_args()
65 kernel = design_kernel(args)
66 L = len(kernel)
67 M = args.past if args.past is not None else max(512, 4 * L)
68 print(f'model: sigma={args.sigma} filter={args.filter} L={L} M={M}')
70 if args.pilot is not None:
71 run_pilot(kernel, args, M)
72 return
74 print(f'{args.pasts} pasts x {args.reps} reps, n0={args.n0} r={args.r} '
75 f'thin={args.thin}')
77 def progress(i, values):
78 mean = float(np.mean(values))
79 se = (float(np.std(values, ddof=1) / np.sqrt(len(values)))
80 if len(values) > 1 else float('nan'))
81 print(f' past {i + 1:3d}/{args.pasts}: H = {values[-1]:.4f} '
82 f'running mean {mean:.4f} +/- {se:.4f}')
84 est = estimate_conditional_entropy(
85 kernel, args.sigma, past=M, pasts=args.pasts, reps=args.reps,
86 n0=args.n0, r=args.r, thin=args.thin, seed=args.seed,
87 progress=progress)
88 ratio = f' (ratio vs int16: {16 / est.mean:.3f}x)' if est.mean > 0 else ''
89 print(f'\nH(z_next | {M} past samples) = {est.mean:.4f} +/- {est.se:.4f} '
90 f'bits/sample{ratio}')
91 print('note: an upper bound on the entropy rate that tightens as --past '
92 'grows.')
95def run_pilot(kernel, args, M):
96 rng = np.random.default_rng(args.seed)
97 print(f'pilot: {args.pasts} pasts, levels 1..{args.pilot}, n0={args.n0}')
98 deltas = np.array([
99 level_corrections(
100 ConditionalChain(kernel, args.sigma, M, rng, args.thin).draw,
101 args.n0, args.pilot)
102 for _ in range(args.pasts)])
103 rms = np.sqrt((deltas ** 2).mean(axis=0))
104 for m in range(args.pilot):
105 note = ''
106 if m > 0 and rms[m] > 0:
107 note = f' decay exponent {np.log2(rms[m - 1] / rms[m]) * 2:.2f}'
108 print(f' m={m + 1}: rms Delta = {rms[m]:.5f}{note}')
109 print('choose r safely below the E[Delta^2] decay exponent (and > 1); '
110 'r=1.5 suits decay near 2.')
113if __name__ == '__main__':
114 main()