c4b8c52Check the model against real ephys, and measure FLACJeremy Magland 1"""The fractional-phase loss as a universal function of s alone.
3Given a real-valued prediction mu and prediction-error std s (both in
4quantization steps), write d = mu - round(mu) for the fractional phase.
5The conditional law of the integer residual r = round(y) - round(mu) is
7 P(r = k | d) = Phi((k + 1/2 - d)/s) - Phi((k - 1/2 - d)/s).
9A coder that knows mu pays the phase-averaged conditional entropy
10 G(s) = E_d[ H(P(. | d)) ] (= theory.gauss_uniform_entropy)
11A coder that codes the integer residual with one pooled histogram pays the
12entropy of the phase-mixture
13 M(s) = H( E_d[ P(. | d) ] ).
14The difference L(s) = M(s) - G(s) >= 0 is a mutual information: what the
15integer residual throws away about the phase. It depends on nothing but s.
16"""
17import numpy as np
18from scipy.special import ndtr
21def phase_entropies(s, kmax=None, nphase=2001):
22 s = float(s)
23 if s <= 0:
24 return 0.0, 0.0, 0.0
25 kmax = kmax or max(4, int(np.ceil(8 * s)) + 2)
26 k = np.arange(-kmax, kmax + 1)[:, None]
27 d = np.linspace(-0.5, 0.5, nphase)[None, :]
28 p = ndtr((k + 0.5 - d) / s) - ndtr((k - 0.5 - d) / s)
29 p = np.clip(p, 1e-300, None)
30 p /= p.sum(axis=0, keepdims=True)
32 def H(q):
33 q = np.clip(q, 1e-300, None)
34 return -(q * np.log2(q)).sum(axis=0)
36 w = np.ones(nphase); w[0] = w[-1] = 0.5; w /= w.sum()
37 cond = float(H(p) @ w) # G(s): knows the phase
38 mix = float(H((p * w).sum(axis=1))) # M(s): pools over phases
39 return mix, cond, mix - cond
42def curve(s_values):
43 return np.array([phase_entropies(s) for s in s_values])
46if __name__ == "__main__":
47 print(f"{'s':>7}{'M(s) mixture':>14}{'G(s) cond':>12}{'L(s) loss':>12}")
48 for s in [0.1, 0.2, 0.28, 0.33, 0.4, 0.5, 0.6, 0.74, 0.9, 1.0, 1.1,
49 1.5, 1.66, 2.0, 3.0, 5.0]:
50 m, g, l = phase_entropies(s)
51 print(f"{s:7.2f}{m:14.4f}{g:12.4f}{l:12.4f}")
53 print("\nAgainst the measured real-data gaps (001290 ch0, bandpassed):")
54 print(f" {'s_*':>6}{'predicted L':>13}{'measured gap':>14}")
55 for s_star, measured in [(1.66, 0.072), (1.10, 0.065), (0.74, 0.101),
56 (0.50, 0.193), (0.33, 0.424)]:
57 print(f" {s_star:6.2f}{phase_entropies(s_star)[2]:13.4f}{measured:14.4f}")