1import json
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 2698
  6rng = np.random.default_rng(SEED)
  7
  8
  9def fit_dynamics(h, z, ridge=1e-6):
 10    # h[t+1] = A h[t] + B z[t] + e[t]
 11    X = np.concatenate([h[:-1], z], axis=1)
 12    Y = h[1:]
 13    W = Y.T @ X @ np.linalg.inv(X.T @ X + ridge*np.eye(X.shape[1]))
 14    d, p = h.shape[1], z.shape[1]
 15    return W[:, :d], W[:, d:], Y - X @ W.T
 16
 17
 18def q_factor(F, P):
 19    # Similarity-invariant generalized contraction factor.
 20    Pm = np.linalg.inv(np.linalg.cholesky(P))
 21    M = Pm @ F @ P @ F.T @ Pm.T
 22    return float(np.sqrt(max(np.linalg.eigvalsh(M).max(), 0)))
 23
 24
 25def scalar_sweep():
 26    # h+=q h + bounded disturbance. Exact worst-case radius is sigma/(1-q).
 27    qs = np.array([0.50, 0.70, 0.85, 0.95, 0.99, 1.01, 1.10])
 28    sigma = 0.1
 29    rows = []
 30    for q in qs:
 31        # adversarial disturbance maximizes absolute state and exposes the formula.
 32        h = 0.0
 33        vals = []
 34        for _ in range(600):
 35            e = sigma if h >= 0 else -sigma
 36            h = q*h + e
 37            vals.append(abs(h))
 38        if q < 1:
 39            predicted = sigma/(1-q)
 40            observed = float(np.mean(vals[-50:]))
 41            relerr = abs(observed-predicted)/predicted
 42        else:
 43            predicted = float('inf')
 44            observed = float(np.mean(vals[-20:]))
 45            relerr = None
 46        growth = (float(np.log(max(vals[-1], 1e-300)/max(vals[-21], 1e-300))/20) if q > 1 else None)
 47        predicted_growth = float(np.log(q)) if q > 1 else None
 48        rows.append(dict(q=float(q), predicted_radius=predicted,
 49                         observed_tail_radius=observed, relative_error=relerr,
 50                         diverged=bool(max(vals)>1e6), observed_log_growth_per_step=growth,
 51                         predicted_log_growth_per_step=predicted_growth))
 52    return rows
 53
 54
 55def boundary_sweep():
 56    # Random bounded disturbances: classify boundedness across a fine q sweep.
 57    qs = np.linspace(.85, 1.15, 31)
 58    sigma = .02
 59    rows=[]
 60    for q in qs:
 61        h=0.; maxabs=0.
 62        for _ in range(400):
 63            h=q*h+rng.uniform(-sigma,sigma)
 64            maxabs=max(maxabs,abs(h))
 65        rows.append((float(q), float(maxabs)))
 66    # operational divergence threshold is first q whose final state exceeds 10x stable q=.99 scale
 67    threshold = next((q for q,m in rows if m > 2.0), None)
 68    return rows, threshold
 69
 70
 71def data_fit_and_projection():
 72    # Offline trajectory fit, covariance ellipsoid, and same forcing with/without projection.
 73    d=3; p=2
 74    A=np.array([[.72,.08,0],[-.04,.64,.05],[0,.03,.58]])
 75    B=np.array([[.35,.0],[0,.28],[.05,.12]])
 76    T=5000
 77    z=rng.normal(size=(T,p))*0.3
 78    h=np.zeros((T+1,d))
 79    for t in range(T):
 80        h[t+1]=A@h[t]+B@z[t]+rng.normal(size=d)*.006
 81    Ah,Bh,res=fit_dynamics(h,z)
 82    # diagonal inflated residual ellipsoid (99.5% absolute coordinate containment)
 83    rad=np.quantile(np.abs(res), .995, axis=0)*1.05
 84    P=np.diag((rad/(1-.72))**2) # conservative axis box-to-ellipsoid scale
 85    fit_err=float(np.linalg.norm(A-Ah)/np.linalg.norm(A))
 86    # Deliberately unstable learned/recurrent matrix to make projection effect visible.
 87    F=Ah*1.55
 88    qp=q_factor(F,P)
 89    x=np.zeros(d); xp=np.zeros(d); raw=[]; proj=[]
 90    for _ in range(300):
 91        noise=rng.normal(size=d)*rad
 92        x=F@x+noise
 93        xp=F@xp+noise
 94        energy=float(xp@np.linalg.inv(P)@xp)
 95        if energy>1: xp=xp/np.sqrt(energy)
 96        raw.append(float(x@np.linalg.inv(P)@x)); proj.append(float(xp@np.linalg.inv(P)@xp))
 97    return {
 98        'fit_relative_A_error':fit_err, 'residual_axis_radii':rad.tolist(),
 99        'q_of_test_matrix':qp, 'raw_violation_rate':float(np.mean(np.array(raw)>1)),
100        'projected_violation_rate':float(np.mean(np.array(proj)>1+1e-10)),
101        'raw_max_ellipsoid_energy':float(max(raw)), 'projected_max_ellipsoid_energy':float(max(proj)),
102        'projected_violation_tolerance':1e-10}
103
104
105def main():
106    radii=scalar_sweep()
107    b_rows,boundary=boundary_sweep()
108    fit=data_fit_and_projection()
109    stable=[r for r in radii if r['q']<1]
110    mean_err=float(np.mean([r['relative_error'] for r in stable]))
111    report={
112      'seed':SEED,
113      'predictions':{
114        'P1_boundary_predicted_q':1.0,
115        'P1_observed_operational_boundary_q':boundary,
116        'P2_radius_formula_mean_relative_error':mean_err,
117        'P2_radius_sweep':radii,
118        'P3_projection_predicted_zero_violation':True,
119        'P3_projection_result':fit
120      },
121      'boundary_sweep': [{'q':q,'max_abs_state':m} for q,m in b_rows]
122    }
123    Path('results.json').write_text(json.dumps(report, indent=2))
124    print(json.dumps(report, indent=2))
125
126if __name__=='__main__': main()