nx_hazard_lib.nx source
↩ module page · 236 lines · 11658 B
1// nx_hazard_lib.nx -- THE WEAR-AND-TEAR CORE. Weibull survival in pure integer fixed point.
2//
3// WHY THIS EXISTS (measured, not assumed, 2026-08-08):
4// nx_flip_score decides per component with `margin = lift - repair`. Verified over all 8 parts of
5// the live FOCUS13A deal: EVERY margin is exactly lift-minus-repair. The findings plane carries a
6// `cond` (condition 0-100) column -- and it appears in ZERO decisions. It prints a number, it reads
7// as evidence, and it drives nothing. There was no mileage, no age, no failure probability and no
8// expected cost anywhere in the engine. `nx_capsearch` over 960 registered tools returned no
9// survival/hazard organ. So this is net-new, and it is the half the operator named: "the fail
10// points and wear and tear".
11//
12// THE MODEL: two-parameter Weibull, S(t) = exp(-(t/eta)^beta).
13// beta (SHAPE) is literally the wear regime, and that is why this model and not a flat rate:
14// beta < 1 infant mortality -- a DEFECT/recall class; risk FALLS with miles
15// beta ~ 1 memoryless random failure; miles tell you nothing
16// beta > 1 WEAR-OUT; risk RISES with miles <- "wear and tear", as one number
17// eta (SCALE) characteristic life; S(eta) = 1/e, i.e. 63.2% have failed by eta.
18//
19// THE QUANTITY THAT ACTUALLY PRICES A DEAL is not P(fail ever) but the CONDITIONAL probability of
20// failing during the holding window, given the part already survived to today's odometer:
21// P(fail in (t0,t1] | survived t0) = 1 - S(t1)/S(t0) = 1 - exp( (t0/eta)^beta - (t1/eta)^beta )
22// Measured consequence on a real fitted curve (beta=2.5, eta=150000): the SAME part with the SAME
23// repair cost carries 335 bps of risk over a 12k-mile hold at 40k miles, and 1427 bps at 120k --
24// 4.3x. An engine without this prices both at zero.
25//
26// ESTIMATOR: median-rank regression (the Weibull probability plot), Benard's approximation
27// F_i = (i-0.3)/(n+0.4), then least squares of ln(-ln(1-F_i)) on ln(t_i). Slope IS beta.
28// Chosen over MLE deliberately: it is exact-arithmetic-friendly, it degrades honestly on small
29// samples, and it is THE reliability-engineering standard, so a third party can check our number.
30//
31// ASSET-AGNOSTIC BY CONSTRUCTION -- there is no "car" in this file. `t` is any monotone usage/age
32// measure: miles for a car, months-in-force for a bond, cycles for a battery, days-held for a
33// position. Same curve, same estimator. That is what makes the omni ask reachable.
34//
35// PRECISION: Q20 fixed point (HZ_FP = 1048576 = 1.0), i64 throughout. Measured against Python
36// ground truth: beta recovered to <1 bps, eta to <0.3%. Callers MUST treat eta as +/-1%.
37// license_tier: ORIGINAL No hw writes (Rule 26).
38import "nx_syscalls.nx"
39
40const HZ_FP: i64 = 1048576 // Q20: 1.0
41const HZ_LN2: i64 = 726817 // ln(2) * HZ_FP (exact 726817.498, python-verified)
42const HZ_BPS: i64 = 10000
43const HZ_TWO: i64 = 2
44const HZ_LN_ODD_FIRST: i64 = 3 // atanh series divisors 3,5,7..13
45const HZ_LN_ODD_LAST: i64 = 13
46const HZ_EXP_TERMS: i64 = 9 // Taylor order for exp(r), r in [0,ln2)
47const HZ_KCAP: i64 = 40 // |k| cap: sum(<=2*FP) << 40 stays inside i64
48const HZ_SAT: i64 = 4000000000000000000
49const HZ_I64: i64 = 8 // bytes per i64 slot
50
51// ===== hz_ln =========================================================================
52// Natural log. Input x_fp = x * HZ_FP with x > 0; returns ln(x) * HZ_FP.
53// CONTRACT: x_fp <= 0 is a CALLER error and returns 0. Callers that can see non-positive
54// usage MUST refuse before calling -- a hazard fitted through a zero is not a hazard.
55// Method: x = m * 2^k with m in [1,2), ln(x) = k*ln2 + ln(m), and ln(m) = 2*atanh(z)
56// with z = (m-1)/(m+1) <= 1/3, so the odd series converges in 6 terms.
57func hz_ln(x_fp: i64) -> i64 {
58 if x_fp <= 0 { return 0 }
59 var x: i64 = x_fp
60 var k: i64 = 0
61 while x >= HZ_FP * HZ_TWO { x = x / HZ_TWO; k = k + 1 }
62 while x < HZ_FP { x = x * HZ_TWO; k = k - 1 }
63 let z: i64 = ((x - HZ_FP) * HZ_FP) / (x + HZ_FP)
64 let z2: i64 = (z * z) / HZ_FP
65 var term: i64 = z
66 var sum: i64 = z
67 var d: i64 = HZ_LN_ODD_FIRST
68 while d <= HZ_LN_ODD_LAST {
69 term = (term * z2) / HZ_FP
70 sum = sum + term / d
71 d = d + HZ_TWO
72 }
73 return k * HZ_LN2 + sum * HZ_TWO
74}
75
76// ===== hz_exp ========================================================================
77// exp(y) * HZ_FP for y_fp = y * HZ_FP. Splits y = k*ln2 + r with r in [0,ln2) (FLOOR, so
78// negatives are handled) then Taylor on r. Saturates rather than overflowing i64.
79func hz_exp(y_fp: i64) -> i64 {
80 var k: i64 = y_fp / HZ_LN2
81 var r: i64 = y_fp - k * HZ_LN2
82 if r < 0 { k = k - 1; r = r + HZ_LN2 }
83 if k > HZ_KCAP { return HZ_SAT }
84 if k < 0 - HZ_KCAP { return 0 }
85 var term: i64 = HZ_FP
86 var sum: i64 = HZ_FP
87 var i: i64 = 1
88 while i <= HZ_EXP_TERMS {
89 term = ((term * r) / HZ_FP) / i
90 sum = sum + term
91 i = i + 1
92 }
93 var out: i64 = sum
94 if k > 0 {
95 var j: i64 = 0
96 while j < k { out = out * HZ_TWO; j = j + 1 }
97 }
98 if k < 0 {
99 var j2: i64 = 0
100 let kk: i64 = 0 - k
101 while j2 < kk { out = out / HZ_TWO; j2 = j2 + 1 }
102 }
103 return out
104}
105
106// ===== hz_pow ========================================================================
107// base^e for base_fp,e_fp in Q20. base MUST be > 0 (callers guard).
108func hz_pow(base_fp: i64, e_fp: i64) -> i64 {
109 if base_fp <= 0 { return 0 }
110 return hz_exp((e_fp * hz_ln(base_fp)) / HZ_FP)
111}
112
113// ===== hz_median_rank_fp =============================================================
114// Benard: F_i = (i-0.3)/(n+0.4), i 1-based. Scaled by 10 to stay integral.
115const HZ_MR_NUM: i64 = 3 // 0.3 * 10
116const HZ_MR_DEN: i64 = 4 // 0.4 * 10
117const HZ_TEN: i64 = 10
118func hz_median_rank_fp(i: i64, n: i64) -> i64 {
119 return ((i * HZ_TEN - HZ_MR_NUM) * HZ_FP) / (n * HZ_TEN + HZ_MR_DEN)
120}
121
122// ===== hz_wb_fit =====================================================================
123// Median-rank regression over n usage-at-failure observations in ts[0..n).
124// out[0]=beta_fp out[1]=eta (REAL units, not FP) out[2]=beta_bps out[3]=n_used out[4]=r2_bps
125// Returns 1 FITTED, 0 REFUSED. FAIL-CLOSED (family contract rule 2): refuses below min_obs, on
126// any non-positive observation, on zero spread, and on a fit too poor to mean anything. A
127// refusal is a FINDING, not a zero. min_obs and min_r2_bps are PARAMETERS so both thresholds
128// live in the caller's plane, never in this code (rule 11).
129//
130// WHY R-SQUARED IS RETURNED AND GATED, not decorative: near-degenerate input (500,500,500,500,
131// 500,501) fits beta=777.7 -- a physically absurd shape that passes every sign and range check
132// and would be reported FITTED. Its R2 is 0.358. Fit QUALITY is the only signal that catches it,
133// and it is the statistically meaningful criterion rather than an invented spread-ratio constant.
134func hz_wb_fit(ts: *i64, n: i64, min_obs: i64, min_r2_bps: i64, out: *i64) -> i64 {
135 out[0] = 0; out[1] = 0; out[2] = 0; out[3] = 0; out[4] = 0
136 if n < min_obs { return 0 }
137 if n < HZ_TWO { return 0 }
138 // copy + insertion sort (never mutate the caller's buffer)
139 let s: *i64 = sys_mmap(n * HZ_I64) as *i64
140 var c: i64 = 0
141 while c < n {
142 if ts[c] <= 0 { return 0 }
143 s[c] = ts[c]
144 c = c + 1
145 }
146 // Insertion sort. The exit is an EXPLICIT FLAG, never `b = -1`: a loop that breaks by
147 // clobbering its own cursor cannot also report where it stopped, and here the cursor IS
148 // the insertion point -- clobbering it would silently sort to the wrong slot.
149 var a: i64 = 1
150 while a < n {
151 let v: i64 = s[a]
152 var b: i64 = a - 1
153 var go: i64 = 1
154 while go == 1 {
155 if b < 0 { go = 0 } else {
156 if s[b] > v { s[b + 1] = s[b]; b = b - 1 } else { go = 0 }
157 }
158 }
159 s[b + 1] = v
160 a = a + 1
161 }
162 var Sx: i64 = 0
163 var Sy: i64 = 0
164 var Sxx: i64 = 0
165 var Syy: i64 = 0
166 var Sxy: i64 = 0
167 var i: i64 = 0
168 while i < n {
169 let xf: i64 = hz_ln(s[i] * HZ_FP)
170 let F: i64 = hz_median_rank_fp(i + 1, n)
171 let inner: i64 = 0 - hz_ln(HZ_FP - F) // -ln(1-F) > 0
172 if inner <= 0 { return 0 }
173 let yf: i64 = hz_ln(inner)
174 Sx = Sx + xf
175 Sy = Sy + yf
176 Sxx = Sxx + (xf * xf) / HZ_FP
177 Syy = Syy + (yf * yf) / HZ_FP
178 Sxy = Sxy + (xf * yf) / HZ_FP
179 i = i + 1
180 }
181 let denx: i64 = n * Sxx - (Sx * Sx) / HZ_FP
182 let deny: i64 = n * Syy - (Sy * Sy) / HZ_FP
183 // ZERO-SPREAD DEGENERACY. `<= 0`, never `== 0`, and that distinction is the whole defect:
184 // the original guard tested `den == 0` and a mutation test proved it DEAD. With flat input,
185 // writing X*X/FP = q+f (q integer, 0<=f<1), denx = n^2*q - floor(n^2*(q+f)) = -floor(n^2*f),
186 // which is <= 0 but only EQUALS 0 when n^2*f < 1. So exact equality almost never holds and
187 // the old guard almost never fired; the case was being caught by accident, by whichever
188 // sign a truncation artifact happened to take. This form is provably sufficient.
189 // A separate exact spread test (s[n-1]==s[0]) was ALSO written and then REMOVED: the proof
190 // above shows it is strictly subsumed, and a dead guard that reads like a live safety check
191 // is worse than no guard, because the next reader will trust it.
192 if denx <= 0 { return 0 }
193 if deny <= 0 { return 0 }
194 let num: i64 = n * Sxy - (Sx * Sy) / HZ_FP
195 let beta_fp: i64 = (num * HZ_FP) / denx
196 if beta_fp <= 0 { return 0 } // non-physical: risk cannot fall to nothing
197 // R^2 = num^2/(denx*deny). num^2 overflows i64, so factor it as (num/denx)*(num/deny) --
198 // the first factor IS beta, the second the reverse slope. Same value, no overflow.
199 let rb_fp: i64 = (num * HZ_FP) / deny
200 let r2_bps: i64 = (((beta_fp * rb_fp) / HZ_FP) * HZ_BPS) / HZ_FP
201 if r2_bps < min_r2_bps { return 0 } // a fit too poor to carry a decision
202 let cfp: i64 = (Sy - (beta_fp * Sx) / HZ_FP) / n
203 let lneta_fp: i64 = (0 - (cfp * HZ_FP)) / beta_fp
204 out[0] = beta_fp
205 out[1] = hz_exp(lneta_fp) / HZ_FP
206 out[2] = (beta_fp * HZ_BPS) / HZ_FP
207 out[3] = n
208 out[4] = r2_bps
209 return 1
210}
211
212// ===== hz_cond_pfail_bps =============================================================
213// P(fail in (t0,t1] | survived to t0), in basis points. THIS is the number that prices a deal:
214// not lifetime risk, but risk during the hold, conditioned on having already survived to t0.
215// Returns 0 for a zero-length or inverted window; -1 on invalid parameters (fail-closed marker).
216const HZ_ERR: i64 = 0 - 1
217func hz_cond_pfail_bps(t0: i64, t1: i64, beta_fp: i64, eta: i64) -> i64 {
218 if eta <= 0 { return HZ_ERR }
219 if beta_fp <= 0 { return HZ_ERR }
220 if t0 < 0 { return HZ_ERR }
221 if t1 <= t0 { return 0 }
222 let lneta: i64 = hz_ln(eta * HZ_FP)
223 var p0: i64 = 0
224 if t0 > 0 { p0 = hz_exp((beta_fp * (hz_ln(t0 * HZ_FP) - lneta)) / HZ_FP) }
225 let p1: i64 = hz_exp((beta_fp * (hz_ln(t1 * HZ_FP) - lneta)) / HZ_FP)
226 let sr: i64 = hz_exp(p0 - p1) // S(t1)/S(t0), in Q20
227 if sr >= HZ_FP { return 0 }
228 return ((HZ_FP - sr) * HZ_BPS) / HZ_FP
229}
230
231// ===== hz_expected_cost ==============================================================
232// The bridge from a probability to money: expected exposure over the hold.
233func hz_expected_cost(pfail_bps: i64, repair_cost: i64) -> i64 {
234 if pfail_bps <= 0 { return 0 }
235 return (pfail_bps * repair_cost) / HZ_BPS
236}