1import json
  2import numpy as np
  3from scipy.optimize import linprog
  4
  5SEED=11
  6rng=np.random.default_rng(SEED)
  7mu,sigma=0.055,0.20
  8U_EDGE,T_EDGE=-0.80,0.80
  9x=np.linspace(-1.,1.,201); N=len(x)
 10U=x<=U_EDGE; T=x>=T_EDGE; C=~(U|T); X0=(x>=-.2)&(x<=.2)
 11z,qw=np.polynomial.hermite.hermgauss(61); z*=np.sqrt(2.); qw/=np.sqrt(np.pi)
 12uix,tix=np.flatnonzero(U)[0],np.flatnonzero(T)[-1]
 13P=np.zeros((N,N))
 14for i,xi in enumerate(x):
 15    if U[i]: P[i,uix]=1.; continue
 16    if T[i]: P[i,tix]=1.; continue
 17    for y,wt in zip(xi+mu+sigma*z,qw):
 18        if y<=U_EDGE: P[i,uix]+=wt
 19        elif y>=T_EDGE: P[i,tix]+=wt
 20        else:
 21            j=np.clip(np.searchsorted(x,y)-1,0,N-2); a=(y-x[j])/(x[j+1]-x[j])
 22            P[i,j]+=wt*(1-a); P[i,j+1]+=wt*a
 23
 24def oracle(delta):
 25    tv=N; c=np.zeros(N+1); c[tv]=1.; A=[]; b=[]
 26    for i in np.flatnonzero(C):
 27        r=np.zeros(N+1); r[:N]=P[i]; r[i]-=1.; A.append(r); b.append(-delta)
 28    for i in np.flatnonzero(X0):
 29        r=np.zeros(N+1); r[i]=1.; r[tv]=-1.; A.append(r); b.append(0.)
 30    bounds=[(1,1) if U[i] else ((0,0) if T[i] else (0,1)) for i in range(N)]+[(0,1)]
 31    r=linprog(c,A_ub=np.asarray(A),b_ub=np.asarray(b),bounds=bounds,method='highs')
 32    return (None,None) if not r.success else (float(r.x[tv]),r.x[:N])
 33
 34def exact_failure():
 35    ci=np.flatnonzero(C); h=np.zeros(N); h[U]=1
 36    h[ci]=np.linalg.solve(np.eye(len(ci))-P[np.ix_(ci,ci)],P[np.ix_(ci,[uix])].ravel())
 37    return h
 38
 39def margin_sweep():
 40    out=[]
 41    for d in [0,.002,.005,.01,.02,.03,.04,.05]:
 42        t,B=oracle(d); out.append({'delta':d,'min_initial_B':t,'p_cert':None if t is None else 1-t,'feasible':t is not None})
 43    return out
 44
 45def mc_scaling():
 46    def bf(y): return np.clip((y+1)/2,0,1)
 47    xx=.1; truth=float(np.sum(qw*bf(xx+mu+sigma*z))); out=[]
 48    for M in [8,16,32,64,128,256,512,1024]:
 49        ee=[]
 50        for _ in range(400): ee.append(abs(bf(xx+mu+sigma*rng.normal(size=M)).mean()-truth))
 51        out.append({'M':M,'mae':float(np.mean(ee))})
 52    slope=float(np.polyfit(np.log([a['M'] for a in out]),np.log([a['mae'] for a in out]),1)[0])
 53    return truth,out,slope
 54
 55def neural_comparison(h,delta=.02,steps=2500):
 56    try:
 57        import torch
 58        torch.manual_seed(SEED); torch.set_num_threads(4)
 59        dev='cuda' if torch.cuda.is_available() else 'cpu'
 60        # Small scalar MLP; P is a deterministic quadrature estimate of the stochastic expectation.
 61        xx=torch.tensor(x[:,None],dtype=torch.float32,device=dev); PP=torch.tensor(P,dtype=torch.float32,device=dev)
 62        UU=torch.tensor(U,device=dev); TT=torch.tensor(T,device=dev); CC=torch.tensor(C,device=dev); XX=torch.tensor(X0,device=dev)
 63        def run(kind):
 64            net=torch.nn.Sequential(torch.nn.Linear(1,32),torch.nn.Tanh(),torch.nn.Linear(32,32),torch.nn.Tanh(),torch.nn.Linear(32,1),torch.nn.Sigmoid()).to(dev)
 65            opt=torch.optim.Adam(net.parameters(),lr=3e-3)
 66            target=torch.tensor(h[:,None],dtype=torch.float32,device=dev)
 67            for k in range(steps):
 68                B=net(xx).squeeze(1); PB=PP@B
 69                boundary=((B[UU]-1)**2).mean()+(B[TT]**2).mean()
 70                if kind=='baseline': loss=boundary+2*((B-target)**2).mean()
 71                else:
 72                    sp=torch.nn.functional.softplus
 73                    loss=8*boundary+8*sp(B[XX]-(1-delta)).mean()+8*sp(PB[CC]-B[CC]+delta).mean()
 74                opt.zero_grad(); loss.backward(); opt.step()
 75            with torch.no_grad():
 76                B=net(xx).squeeze(1); PB=PP@B
 77                return {'loss':float(loss),'max_init_B':float(B[XX].max()),'max_init_limit_violation':float((B[XX]-(1-delta)).max()),'max_drift_residual':float((PB[CC]-B[CC]+delta).max()),'boundary_rmse':float(torch.sqrt(((B[UU]-1)**2).mean()+(B[TT]**2).mean()))}
 78        return {'device':dev,'baseline':run('baseline'),'idea':run('idea')}
 79    except Exception as e: return {'error':repr(e),'fallback':'neural comparison unavailable'}
 80
 81def rollout(n=30000):
 82    ps=[]
 83    for s in np.linspace(-.2,.2,9):
 84        ok=0
 85        for _ in range(n//9):
 86            y=s
 87            for _ in range(300):
 88                if y<=U_EDGE: break
 89                if y>=T_EDGE: ok+=1; break
 90                y+=mu+sigma*rng.normal()
 91        ps.append(ok/(n//9))
 92    return ps
 93
 94def main():
 95    h=exact_failure(); sweep=margin_sweep(); truth,mc,slope=mc_scaling(); nn=neural_comparison(h)
 96    d1,d2=sweep[1],sweep[2]
 97    slope_delta=(d2['min_initial_B']-d1['min_initial_B'])/(d2['delta']-d1['delta'])
 98    feas=[r for r in sweep if r['feasible']]; _,B=oracle(feas[-1]['delta'])
 99    result={'seed':SEED,'system':{'mu':mu,'sigma':sigma,'U':U_EDGE,'T':T_EDGE},'predictions':{'margin_slope_observed':slope_delta,'margin_prediction':'positive approximately linear until infeasible','mc_slope_observed':slope,'mc_prediction':-.5,'feasibility_cliff_between_delta': [.03,.04]},'margin_sweep':sweep,'mc_true_expectation':truth,'mc_scaling':mc,'exact_failure_max_X0':float(h[X0].max()),'rollout_success_by_start':rollout(),'certificate_check':{'delta':feas[-1]['delta'],'max_B_X0':float(B[X0].max()),'max_exact_failure_X0':float(h[X0].max()),'conservative_gap':float(B[X0].max()-h[X0].max())},'neural_comparison':nn}
100    with open('results.json','w') as f: json.dump(result,f,indent=2)
101    print(json.dumps(result,indent=2))
102if __name__=='__main__': main()