nx_sim_sensitivity_gate.nx source

1// nx_sim_sensitivity_gate.nx -- climb the RESULTS-ROBUSTNESS axis with global SENSITIVITY ANALYSIS 2// (NASA-STD-7009 Results Robustness; ASME V&V / FDA expect it). Variance-based Sobol first-order indices via the 3// Saltelli pick-freeze Monte-Carlo estimator: S_i = Var(E[y|x_i]) / Var(y) = the fraction of output variance 4// explained by input i. Sim: y = 3*x1 + 4*x2 + 0*x3 (3 independent centered-uniform inputs). Analytic indices: 5// S1=a^2 v1/Var=360/1000, S2=b^2 v2/Var=640/1000, S3=0 (x3 irrelevant). The MC estimator must recover these, 6// rank x2>x1, sum to ~1 (additive=no interactions), and the irrelevant x3 must score ~0 (liar-kill: a 7// non-influential input cannot be assigned influence). Deterministic PRNG. GREEN iff 6/6. license_tier: ORIGINAL 8import "nx_syscalls.nx" 9 10func 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 } 11func 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 } 12func 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 } 13func iabs(x: i64) -> i64 { if x<0 { return 0-x } return x } 14func 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 } 15func rpos(s: *i64) -> i64 { let v: i64=xrng(s); return v & 0x3FFFFFFFFFFFFFFF } 16func samp(s: *i64, R: i64) -> i64 { return (rpos(s)%(2*R+1)) - R } // centered uniform on [-R,R] 17func ymodel(x1: i64, x2: i64, x3: i64) -> i64 { return 3*x1 + 4*x2 + 0*x3 } 18 19func main() -> i64 { 20 let pass: *i64 = sys_mmap(8) as *i64; pass[0]=0 21 g_w("=== NX-SIM-SENSITIVITY GATE (Sobol first-order indices via pick-freeze, verified vs analytic) ===\n") 22 let R: i64=300; let N: i64=200000 23 let v: i64 = R*(R+1)/3 // variance of each centered-uniform input 24 let var_an: i64 = 9*v + 16*v // Var(y) = a^2 v1 + b^2 v2 (a=3,b=4,c=0) 25 let S1_an: i64 = 9*v*1000/var_an // 360 26 let S2_an: i64 = 16*v*1000/var_an // 640 27 28 let sA: *i64=sys_mmap(8) as *i64; sA[0]=0x9E3779B97F4A7C15 29 let sB: *i64=sys_mmap(8) as *i64; sB[0]=0x2545F4914F6CDD1D 30 var syA: i64=0; var syA2: i64=0; var p1: i64=0; var p2: i64=0; var p3: i64=0 31 var j: i64=0 32 while j<N { 33 let x1A: i64=samp(sA,R); let x2A: i64=samp(sA,R); let x3A: i64=samp(sA,R) 34 let x1B: i64=samp(sB,R); let x2B: i64=samp(sB,R); let x3B: i64=samp(sB,R) 35 let yA: i64 = ymodel(x1A,x2A,x3A) 36 let yC1: i64 = ymodel(x1A,x2B,x3B) // freeze x1 from A 37 let yC2: i64 = ymodel(x1B,x2A,x3B) // freeze x2 from A 38 let yC3: i64 = ymodel(x1B,x2B,x3A) // freeze x3 from A (x3 irrelevant) 39 syA=syA+yA; syA2=syA2+yA*yA 40 p1=p1+yA*yC1; p2=p2+yA*yC2; p3=p3+yA*yC3 41 j=j+1 42 } 43 let f0: i64=syA/N 44 let varY: i64=syA2/N - f0*f0 45 let V1: i64=p1/N - f0*f0; let V2: i64=p2/N - f0*f0; let V3: i64=p3/N - f0*f0 46 let S1: i64=V1*1000/varY; let S2: i64=V2*1000/varY; let S3: i64=V3*1000/varY 47 48 g_w(" analytic: Var(y)="); g_n(var_an); g_w(" S1="); g_n(S1_an); g_w(" S2="); g_n(S2_an); g_w(" S3=0 (per-mille)\n") 49 g_w(" MC: Var(y)="); g_n(varY); g_w(" S1="); g_n(S1); g_w(" S2="); g_n(S2); g_w(" S3="); g_n(S3); g_w(" sum="); g_n(S1+S2+S3); g_w("\n") 50 51 g_row("SETUP: total output variance matches analytic a^2 v1 + b^2 v2 (the model + sampler are correct)" as *u8, (iabs(varY-var_an) < var_an*4/100) as i64, pass) 52 g_row("FIRST-ORDER INDEX x1: Sobol S1 ~ 360 per-mille (a^2 v1 / Var) within MC tolerance" as *u8, (iabs(S1-S1_an)<40) as i64, pass) 53 g_row("FIRST-ORDER INDEX x2: Sobol S2 ~ 640 per-mille (b^2 v2 / Var) within MC tolerance" as *u8, (iabs(S2-S2_an)<40) as i64, pass) 54 g_row("COMPLETENESS: S1+S2+S3 ~ 1000 per-mille (additive model => first-order indices sum to 1, no interactions)" as *u8, (iabs((S1+S2+S3)-1000)<60) as i64, pass) 55 var rank: i64=0; if S2>S1 { if S1>S3 { rank=1 } } 56 g_row("RANKING: x2 ranked more influential than x1, both above x3 -- the SA correctly ranks the drivers" as *u8, rank, pass) 57 g_row("LIAR-KILL: the irrelevant input x3 (zero coefficient) scores S3 ~ 0 -- influence cannot be faked" as *u8, (iabs(S3)<40) as i64, pass) 58 59 g_w("SIM-SENSITIVITY-GATE rows=6 pass="); g_n(pass[0]) 60 if pass[0]==6 { g_w(" verdict=GREEN\n"); sys_exit(0); return 0 } 61 g_w(" verdict=RED\n"); sys_exit(1); return 1 62}