1import json, math, random
  2from pathlib import Path
  3import numpy as np
  4from sklearn.metrics import roc_auc_score
  5
  6SEED=2467
  7np.random.seed(SEED); random.seed(SEED)
  8
  9
 10def hankel_pair(y,p):
 11    y=np.asarray(y,float).reshape(-1)
 12    n=len(y)-2*p+1
 13    ym=np.stack([y[t:t+p] for t in range(n)],axis=1)
 14    yp=np.stack([y[t+p:t+2*p] for t in range(n)],axis=1)
 15    return ym,yp
 16
 17
 18def ridge_map(X,Y,lam=1e-8):
 19    return Y@[email protected](X@X.T+lam*np.eye(X.shape[0]))
 20
 21
 22def math_checks():
 23    # AR(2): delay coordinates are an exact finite Koopman invariant subspace.
 24    a,b=0.82,-0.21
 25    y=np.zeros(500); y[:2]=[0.7,-0.1]
 26    for t in range(2,len(y)): y[t]=a*y[t-1]+b*y[t-2]
 27    X,Y=hankel_pair(y,2); K=ridge_map(X,Y,1e-10)
 28    exact_rel=np.linalg.norm(Y-K@X)/np.linalg.norm(Y)
 29    # Additive observation noise: irreducible one-step residual should scale linearly in sigma.
 30    sigmas=np.array([0.002,0.005,0.01,0.02,0.05,0.1])
 31    rms=[]
 32    for s in sigmas:
 33        vals=[]
 34        for rep in range(12):
 35            yn=y+s*np.random.default_rng(1000+rep).normal(size=len(y))
 36            xx,yy=hankel_pair(yn,2); kk=ridge_map(xx,yy,1e-3)
 37            vals.append(np.sqrt(np.mean((yy-kk@xx)**2)))
 38        rms.append(float(np.mean(vals)))
 39    slope,inter=np.polyfit(sigmas,rms,1)
 40    pred_noise={"prediction":"Koopman residual is linear in observation noise sigma","slope":float(slope),"intercept":float(inter),"r2":float(1-np.sum((rms-(slope*sigmas+inter))**2)/np.sum((rms-np.mean(rms))**2)),"observed_pairs":[[float(s),float(r)] for s,r in zip(sigmas,rms)]}
 41    # Scalar latent rollout: log norm slope is log|a|, boundary is |a|=1.
 42    factors=np.array([0.7,0.9,0.98,1.0,1.02,1.1,1.3])
 43    slopes=[]
 44    H=60
 45    for q in factors:
 46        vals=np.abs(q**np.arange(H)+1e-30)
 47        slopes.append(float(np.polyfit(np.arange(H),np.log(vals),1)[0]))
 48    boundary_factor=float(factors[np.argmin(np.abs(np.array(slopes)))])
 49    stability={"prediction":"rollout growth-rate log(norm_t)/t equals log|a| and changes sign at |a|=1","factors":[float(x) for x in factors],"observed_log_slopes":slopes,"predicted_log_slopes":[float(np.log(x)) for x in factors],"estimated_boundary":boundary_factor,"predicted_boundary":1.0,"max_boundary_error":abs(boundary_factor-1.0)}
 50    # Hankel delay depth: exact AR(2) residual is zero at p>=2, but p=1 cannot represent it.
 51    ps=[1,2,3,4,6]
 52    delay_res=[]
 53    for p in ps:
 54        xx,yy=hankel_pair(y,p); kk=ridge_map(xx,yy,1e-9)
 55        delay_res.append(float(np.linalg.norm(yy-kk@xx)/np.linalg.norm(yy)))
 56    delay={"prediction":"delay depth p>=2 closes the AR(2) Koopman residual; p=1 does not","p":ps,"observed_relative_residual":delay_res,"p2_to_p4_ratio":delay_res[1]/max(delay_res[3],1e-15)}
 57    return {"noise_scaling":pred_noise,"stability_boundary":stability,"delay_transition":delay}
 58
 59
 60def make_data(ntraj=100,length=34,fault=False):
 61    rng=np.random.default_rng(88 if not fault else 99)
 62    out=[]
 63    for _ in range(ntraj):
 64        y=np.zeros(length); y[:2]=rng.normal(0,0.7,2)
 65        aa,bb=(0.82,-0.21) if not fault else (1.03,-0.21)
 66        for t in range(2,length): y[t]=aa*y[t-1]+bb*y[t-2]+rng.normal(0,0.025)
 67        for t in range(length-7):
 68            out.append((y[t:t+4].astype(np.float32),y[t+4:t+8].astype(np.float32)))
 69    return np.array([x for x,_ in out]),np.array([z for _,z in out])
 70
 71
 72def train_compare():
 73    import torch
 74    import torch.nn as nn
 75    device='cuda' if torch.cuda.is_available() else 'cpu'
 76    try:
 77        torch.manual_seed(SEED)
 78        if device=='cuda': torch.cuda.manual_seed_all(SEED)
 79        xp,xf=make_data(90)
 80        vp,vf=make_data(25)
 81        ap,af=make_data(35, fault=True)
 82        class Model(nn.Module):
 83            def __init__(self,dual):
 84                super().__init__(); self.dual=dual
 85                self.e=nn.Sequential(nn.Linear(4,24),nn.Tanh(),nn.Linear(24,3))
 86                self.dp=nn.Sequential(nn.Linear(3,24),nn.Tanh(),nn.Linear(24,4))
 87                self.df=nn.Sequential(nn.Linear(3,24),nn.Tanh(),nn.Linear(24,4))
 88            def forward(self,x):
 89                z=self.e(x); return self.dp(z),self.df(z)
 90        def fit(dual):
 91            m=Model(dual).to(device); opt=torch.optim.Adam(m.parameters(),lr=3e-3)
 92            X=torch.tensor(xp,device=device); F=torch.tensor(xf,device=device)
 93            g=torch.Generator(device=device); g.manual_seed(SEED)
 94            for step in range(700):
 95                ix=torch.randint(0,len(X),(128,),generator=g,device=device)
 96                hp,hf=m(X[ix]); loss=((hp-X[ix])**2).mean()+(1.0 if dual else 0.0)*((hf-F[ix])**2).mean()
 97                opt.zero_grad(); loss.backward(); opt.step()
 98            return m
 99        def score(m,p,f):
100            with torch.no_grad():
101                hp,hf=m(torch.tensor(p,device=device))
102                # Both models are assessed on the same future prediction task; past AE has no trained future head.
103                if m.dual: return ((hf-torch.tensor(f,device=device))**2).mean(1).cpu().numpy()
104                return ((hp-torch.tensor(p,device=device))**2).mean(1).cpu().numpy()
105        base=fit(False); dual=fit(True)
106        # Calibrate each score on nominal validation, then compare fault-vs-nominal AUC.
107        bnom=score(base,vp,vf); bfault=score(base,ap,af)
108        dnom=score(dual,vp,vf); dfault=score(dual,ap,af)
109        result={"device":device,"baseline_past_AUC":float(roc_auc_score(np.r_[np.zeros(len(bnom)),np.ones(len(bfault))],np.r_[bnom,bfault])),"dual_future_AUC":float(roc_auc_score(np.r_[np.zeros(len(dnom)),np.ones(len(dfault))],np.r_[dnom,dfault])),"baseline_nominal_mean":float(bnom.mean()),"baseline_fault_mean":float(bfault.mean()),"dual_nominal_future_mean":float(dnom.mean()),"dual_fault_future_mean":float(dfault.mean())}
110        return result
111    except Exception as e:
112        return {"error":repr(e),"fallback":"training failed; math checks remain valid"}
113
114if __name__=='__main__':
115    result={"seed":SEED,"math":math_checks(),"mini_experiment":train_compare()}
116    Path('results.json').write_text(json.dumps(result,indent=2))
117    print(json.dumps(result,indent=2))