1import json, math
 2import numpy as np
 3from pk_regularizer import poisson_density, poisson_nll, estimate_attractor, sample_poisson_phases
 4
 5
 6def hist_kl(phases, x, bins=48):
 7    h, _ = np.histogram(phases, bins=bins, range=(-np.pi, np.pi), density=False)
 8    p = (h + 1e-6) / (h.sum() + bins * 1e-6)
 9    centers = -np.pi + (np.arange(bins) + .5) * 2*np.pi/bins
10    q = poisson_density(centers, x)
11    q = q / q.sum()
12    return float(np.sum(p * np.log(p/q)))
13
14
15def main():
16    rng=np.random.default_rng(965)
17    out={'normalization':[],'contraction':[],'likelihood':[],'kl':[]}
18    ph=np.linspace(0,2*np.pi,200000,endpoint=False)
19    for r in [0,.2,.5,.8,.95]:
20        x=r*np.exp(.73j); integ=float(np.mean(poisson_density(ph,x)))
21        out['normalization'].append({'r':r,'integral':integ,'error':abs(integ-1)})
22    x=.37+.19j; probes=np.array([-.7+.1j,.2-.5j,.6+.2j,-.1+.7j])
23    for q in [.2,.5,.8]:
24        ks=np.arange(1,13); es=[]
25        for k in ks:
26            es.append(abs(estimate_attractor([(q,(1-q)*x)]*int(k),probes)-x))
27        slope=float(np.polyfit(ks,np.log(np.maximum(es,1e-16)),1)[0])
28        out['contraction'].append({'q':q,'predicted_log_slope':math.log(q),'observed_log_slope':slope,'relative_error':abs(slope-math.log(q))/abs(math.log(q)),'error_K1':es[0],'error_K12':es[-1]})
29    # Under normalized Lebesgue, uniform density is 1 and NLL is 0.
30    for r in [0,.2,.5,.8]:
31        x=r*np.exp(.41j); ps=sample_poisson_phases(x,200000,rng)
32        pk=poisson_nll(ps,x); out['likelihood'].append({'r':r,'pk_nll':pk,'uniform_nll':0.0,'pk_minus_uniform':pk})
33        out['kl'].append({'r':r,'empirical_to_pk':hist_kl(ps,x),'empirical_to_uniform':hist_kl(ps,0j)})
34    out['summary']={'max_normalization_error':max(v['error'] for v in out['normalization']),'contraction_relative_errors':[v['relative_error'] for v in out['contraction']],'pk_nll_becomes_more_negative_with_radius':all(out['likelihood'][i]['pk_minus_uniform']<=out['likelihood'][i-1]['pk_minus_uniform'] for i in range(1,4))}
35    with open('toy_results.json','w') as f: json.dump(out,f,indent=2)
36    print(json.dumps(out,indent=2))
37
38if __name__=='__main__': main()