nx_sim_uq_gate.nx source

1// nx_sim_uq_gate.nx -- climb the UNCERTAINTY QUANTIFICATION axis (NASA-STD-7009 "Results Uncertainty"; ASME V&V 2// 20 / FDA all require UQ). A sovereign Monte-Carlo uncertainty-propagation engine, VERIFIED against the exact 3// analytic propagation law. Sim: projectile range R = (2*vy/g)*vx = k*vx (k=5 with vy=10,g=4). The horizontal 4// launch speed vx is UNCERTAIN: uniform over integers [mu-W, mu+W] (analytic variance W(W+1)/3). MC samples vx, 5// runs the sim, and the OUTPUT distribution must satisfy the analytic law: E[R]=k*mu, Var[R]=k^2*Var[vx]. The 6// sigma interval uses our sovereign isqrt. Liar-kill: with an uncertain input the output variance is provably 7// > 0 -- a "certain/deterministic output" claim is rejected. Deterministic PRNG (reproducible). GREEN iff 6/6. 8// license_tier: ORIGINAL 9import "nx_syscalls.nx" 10 11func g_w(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} sys_write(1,s,n); return 0 } 12func g_n(v: i64) -> i64 { var m: i64=v; if m<0{g_w("-");m=0-m} let t:*u8=sys_mmap(24); var k:i64=0; if m==0{t[0]=48 as u8;k=1}; while m>0{t[k]=(48+(m%10)) as u8;m=m/10;k=k+1}; var i:i64=0; let o:*u8=sys_mmap(24); while i<k{o[i]=t[k-1-i];i=i+1}; sys_write(1,o,k); return 0 } 13func g_row(id: *u8, ok: i64, pass: *i64) -> i64 { g_w(" "); g_w(id); g_w(": "); if ok==1 { g_w("OK\n"); pass[0]=pass[0]+1 } else { g_w("FAIL\n") } return 0 } 14func iabs(x: i64) -> i64 { if x<0 { return 0-x } return x } 15func xrng(s: *i64) -> i64 { var x: i64=s[0]; x = x ^ (x << 13); x = x ^ (x >> 7); x = x ^ (x << 17); s[0]=x; return x } 16func rpos(s: *i64) -> i64 { let v: i64=xrng(s); return v & 0x3FFFFFFFFFFFFFFF } 17func isqrt(N: i64) -> i64 { if N<2 { return N } var x: i64=N; var y: i64=(x+1)/2; while y<x { x=y; y=(x + N/x)/2 } return x } 18 19// the sim: projectile range as a function of horizontal launch speed (linear, gain k). 20func sim_range(vx: i64, k: i64) -> i64 { return k*vx } 21 22// MC propagation: out[0]=mean_R, out[1]=var_R, out[2]=count within mu_R +/- sigma_R, out[3]=var_vx(empirical). 23func mc_run(seed: i64, mu: i64, W: i64, k: i64, N: i64, out: *i64) -> i64 { 24 let st: *i64=sys_mmap(8) as *i64; st[0]=seed 25 let span: i64=2*W+1 26 var sx: i64=0; var sx2: i64=0; var sy: i64=0; var sy2: i64=0 27 let ys: *i64=sys_mmap(8*N) as *i64 28 var i: i64=0 29 while i<N { 30 let vx: i64 = (mu-W) + (rpos(st)%span) 31 let R: i64 = sim_range(vx,k) 32 sx=sx+vx; sx2=sx2+vx*vx; sy=sy+R; sy2=sy2+R*R; ys[i]=R 33 i=i+1 34 } 35 let mean_y: i64=sy/N; let var_y: i64=sy2/N - mean_y*mean_y 36 let mean_x: i64=sx/N; let var_x: i64=sx2/N - mean_x*mean_x 37 let sig: i64=isqrt(var_y) 38 var cnt: i64=0; i=0 39 while i<N { if iabs(ys[i]-mean_y)<=sig { cnt=cnt+1 } i=i+1 } 40 out[0]=mean_y; out[1]=var_y; out[2]=cnt; out[3]=var_x 41 return 0 42} 43 44func main() -> i64 { 45 let pass: *i64 = sys_mmap(8) as *i64; pass[0]=0 46 g_w("=== NX-SIM-UQ GATE (Monte-Carlo uncertainty propagation, verified vs the analytic law) ===\n") 47 let mu: i64=1000; let W: i64=300; let k: i64=5; let N: i64=200000 48 let var_x_analytic: i64 = W*(W+1)/3 // variance of discrete uniform on [mu-W, mu+W] 49 let mean_y_analytic: i64 = k*mu // E[R] = k*E[vx] 50 let var_y_analytic: i64 = k*k*var_x_analytic // Var[R] = k^2 Var[vx] 51 52 let a: *i64=sys_mmap(8*8) as *i64; mc_run(0x9E3779B97F4A7C15, mu, W, k, N, a) 53 let b: *i64=sys_mmap(8*8) as *i64; mc_run(0x2545F4914F6CDD1D, mu, W, k, N, b) 54 let sig_y: i64=isqrt(a[1]) 55 56 g_w(" analytic: Var[vx]="); g_n(var_x_analytic); g_w(" E[R]="); g_n(mean_y_analytic); g_w(" Var[R]="); g_n(var_y_analytic); g_w("\n") 57 g_w(" MC seed1: mean_R="); g_n(a[0]); g_w(" var_R="); g_n(a[1]); g_w(" sigma_R="); g_n(sig_y); g_w(" within1sig%="); g_n(a[2]*100/N); g_w(" var_vx="); g_n(a[3]); g_w("\n") 58 g_w(" MC seed2: mean_R="); g_n(b[0]); g_w(" var_R="); g_n(b[1]); g_w("\n") 59 60 // tolerances (MC error for N=2e5 ~0.3%): variance +/-3%, mean +/-1% 61 let vtol: i64=var_y_analytic*3/100; let vxtol: i64=var_x_analytic*3/100; let mtol: i64=mean_y_analytic/100 62 g_row("SAMPLER: empirical input variance matches analytic W(W+1)/3 (the sampler is correct)" as *u8, (iabs(a[3]-var_x_analytic)<vxtol) as i64, pass) 63 g_row("MEAN PROPAGATION: E[R] == k*E[vx] (the mean propagates through the sim)" as *u8, (iabs(a[0]-mean_y_analytic)<mtol) as i64, pass) 64 g_row("VARIANCE PROPAGATION: Var[R] == k^2 * Var[vx] (UQ law, within MC tolerance)" as *u8, (iabs(a[1]-var_y_analytic)<vtol) as i64, pass) 65 var ci: i64=0; if a[2]*100/N > 52 { if a[2]*100/N < 64 { ci=1 } } // uniform: ~57.7% within +/-1 sigma 66 g_row("CONFIDENCE INTERVAL: ~58% of outputs within mean +/- 1 sigma (sigma via sovereign isqrt) -- usable CI" as *u8, ci, pass) 67 g_row("ROBUSTNESS: an independent PRNG seed reproduces the propagated variance (within MC tolerance)" as *u8, (iabs(a[1]-b[1]) < vtol) as i64, pass) 68 var liar: i64=0; if a[1] > 0 { if iabs(a[1]-var_y_analytic)<vtol { liar=1 } } 69 g_row("LIAR-KILL: with uncertain input the output variance is provably > 0 -- a 'certain output' claim is rejected" as *u8, liar, pass) 70 71 g_w("SIM-UQ-GATE rows=6 pass="); g_n(pass[0]) 72 if pass[0]==6 { g_w(" verdict=GREEN\n"); sys_exit(0); return 0 } 73 g_w(" verdict=RED\n"); sys_exit(1); return 1 74}