1import json, math
  2import numpy as np
  3
  4SEED = 1675
  5A = np.array([[[1.18, 0.22], [0.04, 0.82]],
  6              [[0.86, 0.16], [-0.10, 1.26]]], dtype=float)
  7
  8
  9def transition_matrix(p):
 10    # symmetric two-state chain: persistence p; primitive for 0<p<1
 11    return np.array([[p, 1-p], [1-p, p]], dtype=float)
 12
 13
 14def simulate(p, T=12000, burn=1000, seed=0, target=None, lr=0.0, epochs=0):
 15    """Finite-time top exponent; optional scalar persistence controller."""
 16    rng = np.random.default_rng(seed)
 17    v = np.array([1.0, 0.37]); v /= np.linalg.norm(v)
 18    z = 0
 19    vals = []
 20    for _ in range(epochs):
 21        v = np.array([1.0, 0.37]); v /= np.linalg.norm(v); z = 0; s = 0.0
 22        for t in range(T):
 23            w = A[z] @ v
 24            n = np.linalg.norm(w)
 25            s += np.log(n)
 26            v = w / n
 27            z = int(rng.random() >= p) if z == 0 else int(rng.random() < p)
 28        lam = s / T
 29        vals.append(lam)
 30        if target is not None:
 31            # finite-difference d lambda/dp, then gradient descent on squared target error
 32            hp=0.01
 33            lp=lambda_estimate(min(.98,p+hp), T=max(3000,T//3), seed=700+_)
 34            lm=lambda_estimate(max(.02,p-hp), T=max(3000,T//3), seed=900+_)
 35            grad=(lp-lm)/(2*hp)
 36            p = float(np.clip(p - lr * 2.0 * (lam-target) * grad, 0.02, 0.98))
 37    # independent estimate, collecting state frequencies and transitions
 38    counts = np.zeros(2, dtype=int); trans = np.zeros((2,2), dtype=int)
 39    s = 0.0; v = np.array([1.0, 0.37]); v /= np.linalg.norm(v); z = 0
 40    for t in range(T + burn):
 41        old = z
 42        w = A[z] @ v; n = np.linalg.norm(w); v = w/n
 43        z = int(rng.random() >= p) if z == 0 else int(rng.random() < p)
 44        if t >= burn:
 45            s += np.log(n); counts[old] += 1; trans[old,z] += 1
 46    return dict(lam=float(s/T), p=float(p), counts=counts.tolist(), trans=trans.tolist(), controller_trace=vals)
 47
 48
 49def lambda_estimate(p, T=8000, seed=1):
 50    return simulate(p, T=T, burn=1000, seed=seed)['lam']
 51
 52
 53def common_random_lambda(p, uniforms, burn=1000):
 54    """Coupled estimate: identical random uniforms make FD noise small."""
 55    z=0; v=np.array([1.0,.37]); v/=np.linalg.norm(v); s=0.0
 56    for t,u in enumerate(uniforms):
 57        old=z
 58        w=A[z]@v; n=np.linalg.norm(w); v=w/n
 59        z = int(u >= p) if z == 0 else int(u < p)
 60        if t >= burn: s += np.log(n)
 61    return s/(len(uniforms)-burn)
 62
 63
 64def finite_difference_scaling(p=0.5):
 65    # At symmetry p=.5, lambda'(p)=0 is predicted by mode-label exchange.
 66    rng=np.random.default_rng(991)
 67    uniforms=rng.random(8000)
 68    rows=[]
 69    for h in [0.08,0.04,0.02,0.01,0.005]:
 70        d=(common_random_lambda(p+h,uniforms)-common_random_lambda(p-h,uniforms))/(2*h)
 71        rows.append({'h':h, 'derivative':float(d), 'abs_error_vs_predicted_zero':float(abs(d))})
 72    # Away from symmetry, check derivative consistency under h halving.
 73    p2=.35; u2=np.random.default_rng(992).random(8000)
 74    drows=[]
 75    for h in [.04,.02,.01,.005]:
 76        d=(common_random_lambda(p2+h,u2)-common_random_lambda(p2-h,u2))/(2*h)
 77        drows.append({'h':h,'derivative':float(d)})
 78    return {'symmetry_point':{'p':p,'predicted_derivative':0.0,'rows':rows},
 79            'interior_point':{'p':p2,'rows':drows}}
 80
 81def stationary_check(p):
 82    x=simulate(p,T=8000,burn=2000,seed=31)
 83    freq=np.array(x['counts'],float)/sum(x['counts'])
 84    pi=np.array([.5,.5])
 85    # q_j mass equation reduces to pi=pi P; compare empirical frequencies.
 86    residual=float(np.max(np.abs(freq - pi @ transition_matrix(p))))
 87    return {'empirical_pi':freq.tolist(), 'theoretical_pi':[.5,.5], 'mass_equation_max_residual':residual}
 88
 89
 90def boundary_sweep():
 91    rows=[]
 92    for p in [0.05,0.10,0.25,0.50,0.75,0.90,0.95]:
 93        vals=[lambda_estimate(p,10000,s) for s in [2,3,4,5]]
 94        rows.append({'p':p,'mean_lambda':float(np.mean(vals)), 'std_over_seeds':float(np.std(vals,ddof=1))})
 95    return rows
 96
 97
 98def controller_demo():
 99    # Target is the exponent at p=.5; controller starts at p=.2.
100    target=lambda_estimate(.5,20000,77)
101    x=simulate(.2,T=5000,burn=500,seed=100,target=target,lr=.08,epochs=12)
102    return {'target_lambda':target,'initial_p':.2,'final_p':x['p'],'trace':x['controller_trace']}
103
104
105def main():
106    result={
107      'matrices':A.tolist(),
108      'stationary_check_p05':stationary_check(.5),
109      'stationary_check_p095':stationary_check(.95),
110      'smoothness_fd':finite_difference_scaling(.5),
111      'boundary_sweep':boundary_sweep(),
112      'controller':controller_demo(),
113      'predictions':[
114        'For primitive p in (0,1), empirical stationary masses should approach pi=(.5,.5).',
115        'In the smooth interior, central finite-difference derivative error should decrease approximately quadratically with h after accounting for Monte Carlo noise.',
116        'Near p=0 or p=1, mixing slows and finite-time Lyapunov estimates should have larger seed variance than at p=.5.'
117      ]
118    }
119    with open('results.json','w') as f: json.dump(result,f,indent=2)
120    print(json.dumps(result,indent=2))
121
122if __name__=='__main__': main()