1import json
  2from pathlib import Path
  3import numpy as np
  4
  5
  6def gamma_formula(k):
  7    k = float(k)
  8    return (4-k*k + k*np.sqrt(k*k+8*k))/(2*(k*k+2*k+2))
  9
 10
 11def paper_gains(lam2, lamN):
 12    # Article theorem: kappa = lambda2/lambdaN.
 13    k = lam2 / lamN
 14    g = gamma_formula(k)
 15    return g, 1-g, 2*(1-g)/lam2
 16
 17
 18def idea_gains(lam2, lamN):
 19    # Literal task statement: kappa = lambdaN/lambda2 with the same formula.
 20    k = lamN / lam2
 21    g = gamma_formula(k)
 22    return g, 1-g, 2*(1-g)/lam2
 23
 24
 25def modal_matrix(lam, alpha, eps):
 26    q = 1.0 - eps*lam
 27    return np.array([[q, -alpha], [-eps*lam, q-alpha]], dtype=float)
 28
 29
 30def pole_radius(lam, alpha, eps):
 31    return float(np.max(np.abs(np.linalg.eigvals(modal_matrix(lam, alpha, eps)))))
 32
 33
 34def worst_radius(lam2, lamN, alpha, eps, n=10001):
 35    ls = np.linspace(lam2, lamN, n)
 36    rs = np.array([pole_radius(x, alpha, eps) for x in ls])
 37    i = int(np.argmax(rs))
 38    return float(rs[i]), float(ls[i])
 39
 40
 41def run():
 42    pole=[]
 43    for ratio in [1, 2, 5, 10, 30, 100]:
 44        l2, lN = 1., float(ratio)
 45        pg, pa, pe = paper_gains(l2,lN)
 46        ig, ia, ie = idea_gains(l2,lN)
 47        pr, pl = worst_radius(l2,lN,pa,pe)
 48        ir, il = worst_radius(l2,lN,ia,ie)
 49        pole.append({'lambdaN_over_lambda2':ratio,
 50            'paper_predicted_gamma':pg, 'paper_measured_radius':pr,
 51            'paper_abs_error':abs(pr-pg), 'paper_worst_lambda':pl,
 52            'idea_predicted_gamma':ig, 'idea_measured_radius':ir,
 53            'idea_abs_error':abs(ir-ig), 'idea_worst_lambda':il})
 54
 55    endpoint=[]
 56    for ratio in [2,5,10,30]:
 57        g,a,e=paper_gains(1.,ratio)
 58        endpoint.append({'ratio':ratio,'predicted_gamma':g,
 59            'radius_lambda2':pole_radius(1.,a,e),
 60            'radius_lambdaN':pole_radius(float(ratio),a,e)})
 61
 62    scaling=[]
 63    for l2,lN in [(0.5,5),(1,10),(2,20),(4,40)]:
 64        g,a,e=paper_gains(l2,lN)
 65        scaling.append({'lambda2':l2,'lambdaN':lN,'epsilon':e,
 66                        'epsilon_times_lambda2':e*l2,'alpha':a,'gamma':g})
 67
 68    boundary=[]
 69    l2,lN=1.,10.; g,a,e=paper_gains(l2,lN)
 70    for mult in [.5,.8,1.,1.2]:
 71        ep=mult*e
 72        r,w=worst_radius(l2,lN,a,ep)
 73        boundary.append({'epsilon_over_lower_bound':mult,'radius':r,
 74                         'predicted_at_or_below_gamma':bool(r<=g+1e-8),'worst_lambda':w})
 75
 76    rng=np.random.default_rng(7); n=8
 77    W=np.zeros((n,n))
 78    for i in range(n-1): W[i,i+1]=W[i+1,i]=1.
 79    L=np.diag(W.sum(1))-W
 80    ev=np.linalg.eigvalsh(L); l2,lN=ev[1],ev[-1]
 81    def simulate(alpha,eps,steps=100):
 82        c=rng.normal(size=n); x=rng.normal(size=n); gg=x-c; y=gg.copy(); vals=[]
 83        for _ in range(steps):
 84            xn=x-eps*(L@x)-alpha*y; gn=xn-c
 85            yn=y-eps*(L@y)+gn-gg; x,y,gg=xn,yn,gn
 86            val=float(np.mean(.5*(x-c)**2)) if np.all(np.isfinite(x)) else float('inf')
 87            vals.append(val)
 88            if not np.isfinite(val) or val>1e100: break
 89        return float(vals[-1]), bool(np.isfinite(vals[-1]))
 90    pg,pa,pe=paper_gains(l2,lN); ig,ia,ie=idea_gains(l2,lN)
 91    rng=np.random.default_rng(7); bl, bok=simulate(.1,.1)
 92    rng=np.random.default_rng(7); pl, pok=simulate(pa,pe)
 93    rng=np.random.default_rng(7); il, iok=simulate(ia,ie)
 94    comparison={'graph_lambda2':float(l2),'graph_lambdaN':float(lN),
 95      'baseline_final_loss':bl,'baseline_finite':bok,
 96      'paper_tuned_final_loss':pl,'paper_tuned_finite':pok,
 97      'idea_literal_final_loss':il,'idea_literal_finite':iok}
 98    out={'pole_convention_sweep':pole,'paper_endpoint_prediction':endpoint,
 99         'epsilon_scaling_prediction':scaling,'lower_bound_prediction':boundary,
100         'quadratic_comparison':comparison,
101         'conclusion':'The authoritative idea reverses the paper theorem kappa definition; literal gains are unstable for condition number > 1.'}
102    text=json.dumps(out,indent=2)
103    Path('results.json').write_text(text)
104    print(text)
105
106if __name__=='__main__': run()