import json, math, random from pathlib import Path import numpy as np SEED = 927 np.random.seed(SEED); random.seed(SEED) def math_verification(): # Scalar nested ellipsoids: r_{t+1}=q r_t. This directly tests diam(S_n)<=q^n diam(S_0). q_values = [0.2, 0.5, 0.8, 0.95] decay = [] for q in q_values: r = 1.0 vals = [] for _ in range(20): vals.append(r); r *= q # fit log slope, excluding the first point slope = np.polyfit(np.arange(1, 20), np.log(np.maximum(vals[1:], 1e-30)), 1)[0] decay.append({"q": q, "observed_ratio": float(np.mean(np.array(vals[1:]) / np.array(vals[:-1]))), "predicted_ratio": q, "log_slope": float(slope), "predicted_log_slope": math.log(q)}) # Perturbation recurrence e_{t+1}=rho e_t + eta. The claimed bound is eta/(1-rho). perturb = [] for rho in [0.2, 0.5, 0.8, 0.95]: eta = 0.01 e = 0.0; trace = [] for _ in range(200): e = rho * e + eta; trace.append(e) predicted = eta / (1-rho) perturb.append({"rho": rho, "observed_limit": float(trace[-1]), "predicted_bound": predicted, "relative_error": float(abs(trace[-1]-predicted)/predicted)}) # Stability boundary with gain a=gamma*Lambda and additive perturbation. boundary = [] for a in [0.8, 0.95, 0.99, 1.0, 1.01, 1.1, 1.25]: r = 1e-3; vals=[] for _ in range(80): r = a*r + 1e-4; vals.append(r) boundary.append({"gain_gamma_times_Lambda": a, "radius_at_80": float(vals[-1]), "growth_ratio_last10": float(np.mean(np.array(vals[-10:]) / np.array(vals[-11:-1])))}) return {"decay": decay, "perturbation_bound": perturb, "stability_boundary": boundary} def make_data(n=160, T=24): # Mildly nonlinear, stable 2-D state-space system. A=np.array([[.82,.12],[-.08,.76]], dtype=np.float32) X=[]; U=[] for _ in range(n): x=np.random.randn(2).astype(np.float32)*.7 xs=[]; us=[] for t in range(T): u=np.random.randn(1).astype(np.float32)*.25 xs.append(x.copy()); us.append(u.copy()) x=(A@x + np.array([.08*np.tanh(x[1]), -.06*np.tanh(x[0])],np.float32) + np.array([u[0], .4*u[0]],np.float32) + np.random.randn(2).astype(np.float32)*.008) X.append(xs); U.append(us) return np.array(X), np.array(U) def mini_experiment(): try: import torch import torch.nn as nn torch.manual_seed(SEED) device = "cuda" if torch.cuda.is_available() else "cpu" try: torch.cuda.empty_cache() except Exception: device = "cpu" except Exception: return {"error":"torch unavailable"} X,U=make_data(); split=120 xt=torch.tensor(X[:split],device=device); ut=torch.tensor(U[:split],device=device) xv=torch.tensor(X[split:],device=device); uv=torch.tensor(U[split:],device=device) class Transition(nn.Module): def __init__(self): super().__init__(); self.net=nn.Sequential(nn.Linear(3,32),nn.Tanh(),nn.Linear(32,32),nn.Tanh(),nn.Linear(32,2)) def forward(self,x,u): return self.net(torch.cat([x,u],-1)) def train_baseline(): m=Transition().to(device); opt=torch.optim.Adam(m.parameters(),lr=.008) for _ in range(260): pred=m(xt[:,:-1].reshape(-1,2),ut[:,:-1].reshape(-1,1)) loss=((pred-xt[:,1:].reshape(-1,2))**2).mean(); opt.zero_grad(); loss.backward(); opt.step() return m class RegionModel(nn.Module): def __init__(self): super().__init__(); self.f=Transition(); self.logit_q=nn.Parameter(torch.tensor(-1.5)) def forward(self,x,u): return self.f(x,u) def q(self): return torch.sigmoid(self.logit_q) def train_region(): m=RegionModel().to(device); opt=torch.optim.Adam(m.parameters(),lr=.008) # Fixed-radius observed-compatible regions; q is learned under inclusion and contraction penalties. R=.22; margin=.025 vs=torch.tensor([[1.,0.],[-1.,0.],[0.,1.],[0.,-1.],[.707,.707],[-.707,.707],[.707,-.707],[-.707,-.707]],device=device) for _ in range(300): c=xt[:,:-1].reshape(-1,2); cn=xt[:,1:].reshape(-1,2); uu=ut[:,:-1].reshape(-1,1) # supervised center transition, plus all sampled successor boundary points mapped into predecessor ball center=m(c,uu) q=m.q(); succ=cn[:,None,:] + (R*q)*vs[None,:,:] mapped=m.f(succ.reshape(-1,2),uu[:,None,:].expand(-1,vs.shape[0],-1).reshape(-1,1)).reshape(-1,vs.shape[0],2) d=torch.linalg.vector_norm(mapped-c[:,None,:],dim=-1) inclusion=torch.relu(d-(R-margin)).mean() contraction=torch.relu(q-.93) loss=((center-cn)**2).mean()+2.0*inclusion+0.3*contraction**2 opt.zero_grad(); loss.backward(); opt.step() return m base=train_baseline(); region=train_region() def rollout(m, x0, u, perturb=0.0, is_region=False): c=x0.clone(); preds=[]; radii=[] R=.22; q=float(m.q().detach().cpu()) if is_region else None for t in range(u.shape[1]): if perturb: c=m(c,u[:,t])+torch.randn_like(c)*perturb else: c=m(c,u[:,t]) preds.append(c); radii.append(R*(q**(t+1)) if is_region else 0.) return torch.stack(preds,1), np.array(radii) with torch.no_grad(): pb,_=rollout(base,xv[:,0],uv[:,:-1],0.012) pr,rad=rollout(region,xv[:,0],uv[:,:-1],0.012,True) errb=torch.sqrt(((pb-xv[:,1:])**2).mean()).item() errr=torch.sqrt(((pr-xv[:,1:])**2).mean()).item() one_b=((base(xv[:,0],uv[:,0])-xv[:,1])**2).mean().sqrt().item() one_r=((region(xv[:,0],uv[:,0])-xv[:,1])**2).mean().sqrt().item() # empirical inclusion margin on validation boundary points v=torch.tensor([[1.,0.],[-1.,0.],[0.,1.],[0.,-1.]],device=device) c=xv[:,:-1].reshape(-1,2); u=uv[:,:-1].reshape(-1,1); q=region.q(); cn=xv[:,1:].reshape(-1,2) z=cn[:,None,:]+(.22*q)*v[None,:,:] y=region.f(z.reshape(-1,2),u[:,None,:].expand(-1,4,-1).reshape(-1,1)).reshape(-1,4,2) margins=.22-torch.linalg.vector_norm(y-c[:,None,:],dim=-1) return {"device":device,"baseline_one_step_rmse":one_b,"region_one_step_rmse":one_r, "baseline_noisy_rollout_rmse":errb,"region_noisy_rollout_rmse":errr, "learned_q":float(region.q().detach().cpu()),"validation_min_inclusion_margin":float(margins.min().cpu()), "validation_mean_inclusion_margin":float(margins.mean().cpu()),"region_radius_step_20":float(rad[-1])} if __name__ == '__main__': out={"math":math_verification(),"mini_experiment":mini_experiment()} Path('results.json').write_text(json.dumps(out,indent=2)) print(json.dumps(out,indent=2))