1import json, math, random
 2from pathlib import Path
 3import numpy as np
 4
 5SEED=1723
 6np.random.seed(SEED); random.seed(SEED)
 7
 8def R_step(R,q): return R*R/(2*q)-R+2*q
 9
10def run_map(q,n=80):
11    R=2.0; out=[R]
12    for _ in range(n):
13        R=R_step(R,q); out.append(R)
14        if not np.isfinite(R) or R>1e300: break
15    return np.asarray(out)
16
17def slope(x,y):
18    return float(np.polyfit(x,y,1)[0])
19
20def math_checks():
21    # Three quantitative predictions from the paper's exact recursion.
22    qs=[0.8,0.95,1.0,1.05,1.2]
23    rows=[]
24    for q in qs:
25        a=run_map(q,80); fixed=2*q
26        # q>1: after transient, fixed-R approximately 2q/n (paper Eq.18)
27        if q>1:
28            ns=np.arange(20,min(len(a),70))
29            scaled=np.mean(ns*np.abs(a[ns]-fixed))
30            # q<1: log2 R has slope approximately 2^n, estimated by late ratio
31            metric=scaled
32        elif q<1:
33            ns=np.arange(4,min(len(a)-1,12))
34            vals=np.log(np.maximum(a[ns],1.0))
35            metric=slope(ns, np.log(np.maximum(vals,1e-30)))
36        else: metric=abs(a[-1]-2.0)
37        rows.append({'q':q,'R_last':float(a[-1]),'fixed':fixed,'metric':float(metric)})
38    # Boundary classification over a sweep: below diverges, above converges from R0.
39    boundary=[]
40    for q in np.linspace(.85,1.15,13):
41        a=run_map(float(q),35)
42        boundary.append((float(q), bool(a[-1]>1e6), float(abs(a[-1]-2*q))))
43    # q>1 marginal law: n*(R*-R) tends to 2q, using below-fixed iterates.
44    q=1.2; a=run_map(q,120); ns=np.arange(30,100)
45    marginal=float(np.mean(ns*(2*q-a[ns])))
46    # q<1 double exponential signature: log2 log R increments ~ 1 each generation.
47    q=.8; a=run_map(q,12); v=np.log2(np.log(np.maximum(a[2:],math.e)))
48    dd=float(np.mean(np.diff(v[-5:])))
49    return rows,boundary,marginal,dd
50
51# Minimal vector two-channel block. Its neutral branch is exactly the proposed cubic law;
52# RMS normalization is omitted here so the raw Jacobian prediction is inspectable.
53def block_jacobian(a,b,x,lam):
54    # scalar channels, defect map b' = tanh(.25*a + .25*b), neutral a'=lam^3 a^3+2x^3 b^3
55    J=np.array([[3*lam**3*a*a, 6*x**3*b*b],[.25*(1-math.tanh(.25*a+.25*b)**2),.25*(1-math.tanh(.25*a+.25*b)**2)]])
56    return J
57
58def jacobian_check():
59    # Prediction: cubic neutral local gain is 3 lambda^3 a^2 and scales cubically in lambda.
60    a=.5; b=.2; x=.3
61    vals=[]
62    for lam in [.4,.6,.8,1.0]:
63        J=block_jacobian(a,b,x,lam); rho=max(abs(np.linalg.eigvals(J)))
64        vals.append({'lambda':lam,'rho':float(rho),'neutral_gain':3*lam**3*a*a})
65    # finite difference validates the analytic Jacobian at one point
66    lam=.8; J=block_jacobian(a,b,x,lam)
67    def f(v):
68        aa,bb=v; return np.array([lam**3*aa**3+2*x**3*bb**3, math.tanh(.25*aa+.25*bb)])
69    v=np.array([a,b]); eps=1e-6; Jfd=np.column_stack([(f(v+eps*np.eye(2)[i])-f(v-eps*np.eye(2)[i]))/(2*eps) for i in range(2)])
70    err=float(np.max(abs(J-Jfd)))
71    return vals,err
72
73def toy_regression():
74    # Compare a product-preserving ternary scalar recursion against ordinary mean pooling.
75    # Leaves are positive; target is product of all leaves, so the cubic mechanism is directly relevant.
76    rng=np.random.default_rng(SEED); levels=5; n=3**levels; N=5000
77    X=rng.uniform(.8,1.2,(N,n)); y=np.prod(X,axis=1)
78    # Exact fixed-parameter fractal block (lambda=1,x=0), versus mean pooling.
79    def fractal(x):
80        z=x
81        for _ in range(levels): z=np.prod(z.reshape(-1,3),axis=1)
82        return z
83    pred= np.array([fractal(x) for x in X[:1000]])
84    fractal_mse=float(np.mean((pred-y[:1000])**2))
85    mean=np.mean(X[:1000],axis=1)
86    # best scalar linear readout on mean, fit on first half and evaluate second half
87    A=np.c_[mean[:500],np.ones(500)]; coef=np.linalg.lstsq(A,y[:500],rcond=None)[0]
88    bp=np.c_[mean[500:],np.ones(500)]@coef
89    baseline_mse=float(np.mean((bp-y[500:1000])**2))
90    return {'baseline_mean_mse':baseline_mse,'fractal_product_mse':fractal_mse,'n_leaves':n}
91
92def main():
93    rows,boundary,marginal,dd=math_checks(); jacs,jerr=jacobian_check(); toy=toy_regression()
94    result={'seed':SEED,'paper_map_checks':rows,'boundary_sweep':[{'q':q,'diverged':d,'converged_error':e} for q,d,e in boundary], 'predicted_boundary_q':1.0,'observed_boundary_interval':'[0.975, 1.000] (finite-depth classification)','predicted_marginal_n_error': '2q', 'observed_q1.2_n_error_mean':marginal,'predicted_double_exponential_log2log_increment':1.0,'observed_q0.8_increment':dd,'jacobian_checks':jacs,'jacobian_max_fd_error':jerr,'toy_regression':toy}
95    Path('results.json').write_text(json.dumps(result,indent=2))
96    print(json.dumps(result,indent=2))
97if __name__=='__main__': main()