1import math, random, time
 2from collections import deque
 3import numpy as np
 4import torch
 5import torch.nn as nn
 6
 7SEED = 7
 8random.seed(SEED); np.random.seed(SEED); torch.manual_seed(SEED)
 9
10class HankelResidualObserver:
11    def __init__(self, p=1, t_ini=8, horizon=8, max_cols=40, ridge=1e-3, tau=1.0):
12        self.p, self.t_ini, self.horizon = p, t_ini, horizon
13        self.max_cols, self.ridge, self.tau = p*0 + max_cols, ridge, tau
14        self.residuals = deque(maxlen=t_ini + horizon + max_cols + 4)
15    def update(self, residual):
16        self.residuals.append(np.asarray(residual, dtype=np.float64).reshape(self.p))
17    def predict(self):
18        if len(self.residuals) < self.t_ini + self.horizon:
19            return np.zeros((self.horizon, self.p)), 0.0, 0
20        a = np.asarray(self.residuals)
21        n = len(a) - self.t_ini - self.horizon + 1
22        starts = np.arange(max(0, n - self.max_cols), n)
23        dp = np.stack([a[j:j+self.t_ini].reshape(-1) for j in starts], axis=1)
24        df = np.stack([a[j+self.t_ini:j+self.t_ini+self.horizon].reshape(-1) for j in starts], axis=1)
25        d_ini = a[-self.t_ini:].reshape(-1)
26        g = np.linalg.solve(dp.T @ dp + self.ridge*np.eye(dp.shape[1]), dp.T @ d_ini)
27        recon = dp @ g
28        q = np.linalg.norm(recon-d_ini)/(np.linalg.norm(d_ini)+1e-8)
29        gamma = float(np.clip(1.0-q/self.tau, 0.0, 1.0))
30        return (df @ g).reshape(self.horizon, self.p), gamma, dp.shape[1]
31
32class SmallGRU(nn.Module):
33    def __init__(self, hidden=24, horizon=8):
34        super().__init__(); self.gru=nn.GRU(1,hidden,batch_first=True); self.head=nn.Linear(hidden,horizon)
35    def forward(self,x):
36        z,_=self.gru(x); return self.head(z[:,-1]).unsqueeze(-1)
37
38def make_train(n=5000,w=16,h=8):
39    y=np.zeros(n+w+h+1); noise=np.random.default_rng(SEED).normal(0,.08,len(y))
40    for t in range(1,len(y)): y[t]=.82*y[t-1]+noise[t]
41    return (np.asarray([y[t:t+w] for t in range(n)],np.float32)[...,None],
42            np.asarray([y[t+w:t+w+h] for t in range(n)],np.float32)[...,None])
43
44def fit_model(device):
45    X,Y=make_train(); model=SmallGRU().to(device); opt=torch.optim.Adam(model.parameters(),lr=3e-3); g=torch.Generator().manual_seed(SEED)
46    for _ in range(18):
47        perm=torch.randperm(len(X),generator=g)
48        for i in range(0,len(X),128):
49            ix=perm[i:i+128]; xb=torch.from_numpy(X[ix]).to(device); yb=torch.from_numpy(Y[ix]).to(device)
50            loss=(model(xb)-yb).square().mean(); opt.zero_grad(); loss.backward(); opt.step()
51    return model
52
53def run_stream(model,device,disturbed=True,n=700,w=16,h=8):
54    rng=np.random.default_rng(101 if disturbed else 102); y=np.zeros(n+w+h+2)
55    for t in range(1,len(y)):
56        d=.55*math.sin(.42*t+.4) if disturbed else 0.; y[t]=.82*y[t-1]+d+rng.normal(0,.035)
57    obs=HankelResidualObserver(t_ini=8,horizon=h,max_cols=40,ridge=1e-2,tau=1.0)
58    base_err=[]; corr_err=[]; gammas=[]; lat=[]; model.eval()
59    with torch.no_grad():
60        for t in range(w,w+n-h):
61            x=torch.from_numpy(y[t-w:t].astype(np.float32)[None,:,None]).to(device); pred=model(x).cpu().numpy()[0,:,0]
62            if t>w: obs.update([y[t]-prior_pred[0]])
63            tic=time.perf_counter(); dh,gamma,_=obs.predict(); lat.append(time.perf_counter()-tic)
64            target=y[t:t+h]; base_err.extend(pred-target); corr_err.extend(pred+gamma*dh[:,0]-target); gammas.append(gamma); prior_pred=pred
65    be=np.asarray(base_err); ce=np.asarray(corr_err)
66    return {'baseline_rmse':float(np.sqrt(np.mean(be**2))),'observer_rmse':float(np.sqrt(np.mean(ce**2))),
67      'baseline_long_rmse':float(np.sqrt(np.mean(be.reshape(-1,h)[:,-3:]**2))), 'observer_long_rmse':float(np.sqrt(np.mean(ce.reshape(-1,h)[:,-3:]**2))),
68      'mean_gamma':float(np.mean(gammas)),'p95_latency_ms':float(np.percentile(lat,95)*1000),'n_predictions':len(gammas)}
69
70def math_check():
71    # Exactly period-4 residuals: every length-8 past and future window is a dictionary shift.
72    pattern=np.array([1.,.25,-1.,-.25]); r=np.resize(pattern,96)[:,None]
73    o=HankelResidualObserver(t_ini=8,horizon=8,max_cols=40,ridge=1e-10,tau=1.)
74    for x in r: o.update(x)
75    pred,gamma,cols=o.predict(); truth=np.resize(pattern,8)[:,None]
76    # Random unrelated history should fail reconstruction and suppress correction.
77    z=np.random.default_rng(99).normal(size=(96,1)); qobs=HankelResidualObserver(t_ini=8,horizon=8,max_cols=40,ridge=1e-2,tau=1.)
78    for x in z:qobs.update(x)
79    _,bad_gamma,_=qobs.predict()
80    return {'periodic_prediction_rmse':float(np.sqrt(np.mean((pred-truth)**2))), 'periodic_gamma':gamma,
81            'random_mismatch_gamma':bad_gamma,'dictionary_columns':cols}
82
83def main():
84    device='cuda' if torch.cuda.is_available() else 'cpu'
85    try: model=fit_model(device)
86    except Exception: device='cpu'; model=fit_model(device)
87    if device=='cuda': torch.cuda.empty_cache()
88    print(__import__('json').dumps({'device':device,'math_check':math_check(),'disturbed':run_stream(model,device,True),'clean':run_stream(model,device,False)},indent=2))
89if __name__=='__main__': main()