1import json, math, random
 2import numpy as np
 3
 4# Toy conjugate expanding map F=h^{-1} o A o h, h(x)=x+eps*sin(x), A(z)=a*z.
 5# Exact identity: ell(x)=log(a)+u(F(x))-u(x), u(x)=-log(h'(x)).
 6def h(x, eps): return x + eps*np.sin(x)
 7def hp(x, eps): return 1.0 + eps*np.cos(x)
 8def hinv(z, eps):
 9    x=np.asarray(z,dtype=float).copy()
10    for _ in range(40): x -= (h(x,eps)-z)/hp(x,eps)
11    return x
12def F(x,a,eps): return hinv(a*h(x,eps),eps)
13def ell(x,a,eps):
14    y=F(x,a,eps)
15    return math.log(a)-math.log(hp(y,eps))+math.log(hp(x,eps))
16def u(x,eps): return -np.log(hp(x,eps))
17
18def mechanism_sweep():
19    rng=np.random.default_rng(7); a=1.03; c=math.log(a)
20    xs=rng.uniform(-math.pi,math.pi,1000)
21    # P1: exact potential gives zero finite-horizon residual at every k.
22    telescoping=[]
23    for k in [1,2,4,8,16,32]:
24        vals=[]
25        for x0 in xs[:300]:
26            x=x0; s=0.0
27            for _ in range(k): s+=ell(x,a,.6); x=F(x,a,.6)
28            vals.append(s-(u(x,.6)-u(x0,.6)+k*c))
29        z=np.asarray(vals)
30        telescoping.append({'k':k,'max_abs_R':float(np.max(np.abs(z))),
31                            'std_R_over_k':float(np.std(z/k))})
32    # P2: statewise variation grows quadratically for small eps.
33    epsvals=np.array([.05,.1,.2,.4])
34    variation=[]
35    for e in [0,.05,.1,.2,.4,.6,.8]:
36        z=np.array([ell(x,a,e)-c for x in xs])
37        variation.append({'eps':e,'pointwise_var':float(np.var(z)),
38                          'pointwise_max_abs':float(np.max(np.abs(z)))})
39    log_slope=float(np.polyfit(np.log(epsvals),np.log([
40        np.var([ell(x,a,e)-c for x in xs]) for e in epsvals]),1)[0])
41    # P3: if u is scaled by alpha, R_k/k=(1-alpha)*(u(x0)-u(xk))/k.
42    scaling=[]; e=.6
43    for alpha in [0,.25,.5,.75,1.0]:
44        row={'alpha':alpha}
45        for k in [1,4,16,64]:
46            vals=[]
47            for x0 in xs[:200]:
48                x=x0; s=0.
49                for _ in range(k): s+=ell(x,a,e); x=F(x,a,e)
50                vals.append(s-(alpha*(u(x,e)-u(x0,e))+k*c))
51            row[str(k)]=float(np.std(np.asarray(vals)/k))
52        scaling.append(row)
53    rows=[r for r in scaling if r['alpha']<1]
54    ratios=[]
55    for r in rows:
56        ratios.append({'alpha':r['alpha'],'ratio_to_alpha0_k1':r['1']/scaling[0]['1']})
57    return {'telescoping':telescoping,'variation':variation,
58            'small_eps_loglog_slope_variance_vs_eps':log_slope,
59            'scaled_potential':scaling,
60            'scaled_residual_ratio_prediction':ratios}
61
62def learned_potential_demo():
63    # Baseline: best constant log-Jacobian. Idea: fit a scalar potential and c.
64    rng=np.random.default_rng(11); a=1.03; eps=.6
65    x=rng.uniform(-math.pi,math.pi,1200); y=np.array([F(v,a,eps) for v in x])
66    l=np.array([ell(v,a,eps) for v in x])
67    cols=[]
68    for n in range(1,9):
69        cols += [np.sin(n*y)-np.sin(n*x), np.cos(n*y)-np.cos(n*x)]
70    X=np.column_stack(cols+[np.ones(len(x))])
71    coef=np.linalg.lstsq(X,l,rcond=None)[0]; r=l-X@coef
72    point=np.mean((l-np.mean(l))**2)
73    return {'cohomological_mse':float(np.mean(r*r)),
74            'pointwise_const_mse':float(point),
75            'fitted_c':float(coef[-1]), 'true_c':math.log(a),
76            'mse_ratio_coh_over_point':float(np.mean(r*r)/point)}
77
78def main():
79    random.seed(0); np.random.seed(0)
80    out={'predictions':mechanism_sweep(),'baseline_vs_idea':learned_potential_demo(),
81         'prediction_statements':[
82          'P1 exact conjugacy predicts R_k=0 for every horizon k.',
83          'P2 pointwise variance is predicted to scale as eps^2 near eps=0.',
84          'P3 imperfect-potential residual per step scales as (1-alpha) and has a 1/k boundary decay.'
85         ]}
86    with open('results.json','w') as f: json.dump(out,f,indent=2)
87    print(json.dumps(out,indent=2))
88if __name__=='__main__': main()