1"""Fit the paper's generative model x ~ iid N(0, 1) -> h * x -> round to a
2real integer trace, and synthesize a surrogate from the fit.
4Only the spectrum is identifiable (sigma and the kernel scale are the same
5knob), so sigma is fixed at 1 and h carries everything.
7Fit:
8 1. Welch PSD of z, in the convention mean_{f in [0,1)} S(f) = var(z).
9 2. Dither model: S_z = S_y + 1/12, so S_y = max(S_z - 1/12, floor).
10 3. h = the minimum-phase spectral factor of S_y (real-cepstrum method).
12Prediction (no sampling required):
13 s_*^2 = exp( mean_f ln(S_y + 1/12) ) - 1/12 [Szego on the observed
14 spectrum], and the predicted entropy rate is G(s_*), the entropy of
15 N(0, s_*^2) + U(-1/2, 1/2).
16"""
17import numpy as np
18from scipy.signal import welch, fftconvolve
20FLOOR_FRAC = 1e-4 # S_y is floored at this fraction of its own mean
23def psd(z, nfft=4096):
24 """Two-sided PSD on the rfft grid f = k/nfft, k = 0..nfft/2, normalized so
25 that the mean over f in [0, 1) equals var(z)."""
26 f, pxx = welch(np.asarray(z, dtype=np.float64), fs=1.0, nperseg=nfft,
27 noverlap=nfft // 2, window="hann", detrend="constant",
28 return_onesided=True, scaling="density")
29 return f, pxx / 2.0
32def _grid_mean(v):
33 """Mean over f in [0, 1) of a quantity given on the rfft half-grid
34 (trapezoid on [0, 1/2]; the spectrum is symmetric)."""
35 return float((v.sum() - 0.5 * (v[0] + v[-1])) / (v.size - 1))
38def minimum_phase(s_half, n_taps=None):
39 """Minimum-phase impulse response whose |H(f)|^2 matches s_half on the
40 rfft grid. Real-cepstrum construction."""
41 s_full = np.concatenate([s_half, s_half[-2:0:-1]])
42 n = s_full.size
43 c = np.fft.ifft(0.5 * np.log(s_full)).real # real cepstrum of |H|
44 cm = np.zeros(n)
45 cm[0] = c[0]
46 cm[1:n // 2] = 2.0 * c[1:n // 2]
47 cm[n // 2] = c[n // 2]
48 h = np.fft.ifft(np.exp(np.fft.fft(cm))).real # causal, min phase
49 return h if n_taps is None else h[:n_taps]
52class Fit:
53 def __init__(self, z, nfft=4096, n_taps=None):
54 z = np.asarray(z, dtype=np.float64)
55 self.n = z.size
56 self.nfft = nfft
57 self.freq, self.s_z = psd(z, nfft)
58 floor = FLOOR_FRAC * self.s_z.mean()
59 self.s_y = np.maximum(self.s_z - 1.0 / 12.0, floor)
60 self.kernel = minimum_phase(self.s_y, n_taps)
61 # Szego on the observed spectrum, with roundoff as a 1/12 dither floor
62 gm = np.exp(_grid_mean(np.log(self.s_y + 1.0 / 12.0)))
63 self.s_star = float(np.sqrt(max(gm - 1.0 / 12.0, 0.0)))
64 self.sigma_inf = float(np.sqrt(np.exp(_grid_mean(np.log(self.s_y)))))
67 def predicted_rate(self):
68 """G(s_*) in bits/sample — the analytic entropy-rate prediction."""
69 from timeseries_entropy.theory import gauss_uniform_entropy
70 return float(gauss_uniform_entropy(self.s_star))
72 def kernel_error_db(self):
73 """How well the (possibly truncated) kernel reproduces the fitted
74 spectrum: RMS error in dB over the grid."""
75 h = np.abs(np.fft.rfft(self.kernel, self.nfft)) ** 2
76 return float(np.sqrt(np.mean((10 * np.log10(h / self.s_y)) ** 2)))
78 def synthesize(self, n=None, seed=0):
79 """A surrogate integer trace from the fitted model."""
80 n = self.n if n is None else n
81 rng = np.random.default_rng(seed)
82 k = self.kernel
83 x = rng.standard_normal(n + k.size - 1)
84 y = fftconvolve(x, k, mode="valid")
85 return np.round(y).astype(np.int16)