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