1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4
  5SEED = 927
  6np.random.seed(SEED); random.seed(SEED)
  7
  8
  9def math_verification():
 10    # Scalar nested ellipsoids: r_{t+1}=q r_t. This directly tests diam(S_n)<=q^n diam(S_0).
 11    q_values = [0.2, 0.5, 0.8, 0.95]
 12    decay = []
 13    for q in q_values:
 14        r = 1.0
 15        vals = []
 16        for _ in range(20):
 17            vals.append(r); r *= q
 18        # fit log slope, excluding the first point
 19        slope = np.polyfit(np.arange(1, 20), np.log(np.maximum(vals[1:], 1e-30)), 1)[0]
 20        decay.append({"q": q, "observed_ratio": float(np.mean(np.array(vals[1:]) / np.array(vals[:-1]))),
 21                      "predicted_ratio": q, "log_slope": float(slope),
 22                      "predicted_log_slope": math.log(q)})
 23
 24    # Perturbation recurrence e_{t+1}=rho e_t + eta. The claimed bound is eta/(1-rho).
 25    perturb = []
 26    for rho in [0.2, 0.5, 0.8, 0.95]:
 27        eta = 0.01
 28        e = 0.0; trace = []
 29        for _ in range(200):
 30            e = rho * e + eta; trace.append(e)
 31        predicted = eta / (1-rho)
 32        perturb.append({"rho": rho, "observed_limit": float(trace[-1]),
 33                        "predicted_bound": predicted,
 34                        "relative_error": float(abs(trace[-1]-predicted)/predicted)})
 35
 36    # Stability boundary with gain a=gamma*Lambda and additive perturbation.
 37    boundary = []
 38    for a in [0.8, 0.95, 0.99, 1.0, 1.01, 1.1, 1.25]:
 39        r = 1e-3; vals=[]
 40        for _ in range(80):
 41            r = a*r + 1e-4; vals.append(r)
 42        boundary.append({"gain_gamma_times_Lambda": a, "radius_at_80": float(vals[-1]),
 43                         "growth_ratio_last10": float(np.mean(np.array(vals[-10:]) / np.array(vals[-11:-1])))})
 44    return {"decay": decay, "perturbation_bound": perturb, "stability_boundary": boundary}
 45
 46
 47def make_data(n=160, T=24):
 48    # Mildly nonlinear, stable 2-D state-space system.
 49    A=np.array([[.82,.12],[-.08,.76]], dtype=np.float32)
 50    X=[]; U=[]
 51    for _ in range(n):
 52        x=np.random.randn(2).astype(np.float32)*.7
 53        xs=[]; us=[]
 54        for t in range(T):
 55            u=np.random.randn(1).astype(np.float32)*.25
 56            xs.append(x.copy()); us.append(u.copy())
 57            x=(A@x + np.array([.08*np.tanh(x[1]), -.06*np.tanh(x[0])],np.float32)
 58               + np.array([u[0], .4*u[0]],np.float32) + np.random.randn(2).astype(np.float32)*.008)
 59        X.append(xs); U.append(us)
 60    return np.array(X), np.array(U)
 61
 62
 63def mini_experiment():
 64    try:
 65        import torch
 66        import torch.nn as nn
 67        torch.manual_seed(SEED)
 68        device = "cuda" if torch.cuda.is_available() else "cpu"
 69        try:
 70            torch.cuda.empty_cache()
 71        except Exception:
 72            device = "cpu"
 73    except Exception:
 74        return {"error":"torch unavailable"}
 75    X,U=make_data(); split=120
 76    xt=torch.tensor(X[:split],device=device); ut=torch.tensor(U[:split],device=device)
 77    xv=torch.tensor(X[split:],device=device); uv=torch.tensor(U[split:],device=device)
 78
 79    class Transition(nn.Module):
 80        def __init__(self):
 81            super().__init__(); self.net=nn.Sequential(nn.Linear(3,32),nn.Tanh(),nn.Linear(32,32),nn.Tanh(),nn.Linear(32,2))
 82        def forward(self,x,u): return self.net(torch.cat([x,u],-1))
 83
 84    def train_baseline():
 85        m=Transition().to(device); opt=torch.optim.Adam(m.parameters(),lr=.008)
 86        for _ in range(260):
 87            pred=m(xt[:,:-1].reshape(-1,2),ut[:,:-1].reshape(-1,1))
 88            loss=((pred-xt[:,1:].reshape(-1,2))**2).mean(); opt.zero_grad(); loss.backward(); opt.step()
 89        return m
 90
 91    class RegionModel(nn.Module):
 92        def __init__(self):
 93            super().__init__(); self.f=Transition(); self.logit_q=nn.Parameter(torch.tensor(-1.5))
 94        def forward(self,x,u): return self.f(x,u)
 95        def q(self): return torch.sigmoid(self.logit_q)
 96
 97    def train_region():
 98        m=RegionModel().to(device); opt=torch.optim.Adam(m.parameters(),lr=.008)
 99        # Fixed-radius observed-compatible regions; q is learned under inclusion and contraction penalties.
100        R=.22; margin=.025
101        vs=torch.tensor([[1.,0.],[-1.,0.],[0.,1.],[0.,-1.],[.707,.707],[-.707,.707],[.707,-.707],[-.707,-.707]],device=device)
102        for _ in range(300):
103            c=xt[:,:-1].reshape(-1,2); cn=xt[:,1:].reshape(-1,2); uu=ut[:,:-1].reshape(-1,1)
104            # supervised center transition, plus all sampled successor boundary points mapped into predecessor ball
105            center=m(c,uu)
106            q=m.q(); succ=cn[:,None,:] + (R*q)*vs[None,:,:]
107            mapped=m.f(succ.reshape(-1,2),uu[:,None,:].expand(-1,vs.shape[0],-1).reshape(-1,1)).reshape(-1,vs.shape[0],2)
108            d=torch.linalg.vector_norm(mapped-c[:,None,:],dim=-1)
109            inclusion=torch.relu(d-(R-margin)).mean()
110            contraction=torch.relu(q-.93)
111            loss=((center-cn)**2).mean()+2.0*inclusion+0.3*contraction**2
112            opt.zero_grad(); loss.backward(); opt.step()
113        return m
114
115    base=train_baseline(); region=train_region()
116    def rollout(m, x0, u, perturb=0.0, is_region=False):
117        c=x0.clone(); preds=[]; radii=[]
118        R=.22; q=float(m.q().detach().cpu()) if is_region else None
119        for t in range(u.shape[1]):
120            if perturb:
121                c=m(c,u[:,t])+torch.randn_like(c)*perturb
122            else: c=m(c,u[:,t])
123            preds.append(c); radii.append(R*(q**(t+1)) if is_region else 0.)
124        return torch.stack(preds,1), np.array(radii)
125    with torch.no_grad():
126        pb,_=rollout(base,xv[:,0],uv[:,:-1],0.012)
127        pr,rad=rollout(region,xv[:,0],uv[:,:-1],0.012,True)
128        errb=torch.sqrt(((pb-xv[:,1:])**2).mean()).item()
129        errr=torch.sqrt(((pr-xv[:,1:])**2).mean()).item()
130        one_b=((base(xv[:,0],uv[:,0])-xv[:,1])**2).mean().sqrt().item()
131        one_r=((region(xv[:,0],uv[:,0])-xv[:,1])**2).mean().sqrt().item()
132        # empirical inclusion margin on validation boundary points
133        v=torch.tensor([[1.,0.],[-1.,0.],[0.,1.],[0.,-1.]],device=device)
134        c=xv[:,:-1].reshape(-1,2); u=uv[:,:-1].reshape(-1,1); q=region.q(); cn=xv[:,1:].reshape(-1,2)
135        z=cn[:,None,:]+(.22*q)*v[None,:,:]
136        y=region.f(z.reshape(-1,2),u[:,None,:].expand(-1,4,-1).reshape(-1,1)).reshape(-1,4,2)
137        margins=.22-torch.linalg.vector_norm(y-c[:,None,:],dim=-1)
138        return {"device":device,"baseline_one_step_rmse":one_b,"region_one_step_rmse":one_r,
139                "baseline_noisy_rollout_rmse":errb,"region_noisy_rollout_rmse":errr,
140                "learned_q":float(region.q().detach().cpu()),"validation_min_inclusion_margin":float(margins.min().cpu()),
141                "validation_mean_inclusion_margin":float(margins.mean().cpu()),"region_radius_step_20":float(rad[-1])}
142
143if __name__ == '__main__':
144    out={"math":math_verification(),"mini_experiment":mini_experiment()}
145    Path('results.json').write_text(json.dumps(out,indent=2))
146    print(json.dumps(out,indent=2))