nx_pred_maint_consumer.nx source
↩ module page · 376 lines · 15308 B
1// nx_pred_maint_consumer.nx -- predictive-maintenance consumer.
2//
3// Subscribes to OBD-II shim events (nx_obd2_shim.nx, commit be03adfa)
4// and produces "anomaly flagged" or "predicted-DTC-in-N-days" events
5// AHEAD of OBD's reactive DTC firing. Concretises the
6// INTEROPERABILITY_CHARTER.md §12 worked example:
7// "4-day-early P0117 prediction via aggregated coolant-temp anomaly".
8//
9// This is THE bring-the-most-out-of-it proof per the operator's
10// cardinal: OBD-native tools react to DTCs after they fire; Nishi
11// pulls predictive insight out of the SAME OBD data that OBD's spec
12// itself doesn't define.
13//
14// Status: SEED v0.1.0. 2026-05-26.
15//
16// WINNER-TIER: WINNER-S CANDIDATE (no incumbent matches the
17// value-add; OBD-native tools are reactive-only)
18// INCUMBENTS: python-obd (reactive; no anomaly), obdlib (reactive),
19// ELM327 tools (passthrough), Vector CANalyzer
20// (commercial; has trend tracking but no
21// substrate-integrated prediction)
22// NUMBERS: V1 ships the algorithm; measured "days-early" metric
23// vs the reactive baseline requires real CAN-logger
24// data + a DTC corpus. Provisional.
25// GAP: no incumbent solves this; V1 establishes the
26// baseline. M-next: real OBD trace replay against
27// recorded fleet data; confirm S rating.
28// EXEMPTION REASON: n/a
29//
30// Algorithm (V1):
31// Per PID:
32// - ring buffer of last N samples (timestamp + value)
33// - Welford's online algorithm: rolling mean + variance
34// (numerically stable for streaming data; Welford 1962)
35// - linear regression: slope estimate over the window
36// (Pearson r^2 fit; for "is value drifting?" detection)
37// Anomaly flagging:
38// - 3-sigma rule: current value > mean + 3*sigma (or < mean - 3*sigma)
39// - configurable per-PID thresholds (per ZERO_TO_ADVANCED.md
40// data-driven approach; threshold lives in config not code)
41// Time-to-DTC prediction:
42// - if slope is non-zero and consistent (r^2 > 0.7):
43// extrapolate value to the OBD-spec'd DTC threshold;
44// compute days until crossing
45//
46// Future V2+:
47// - GEMM-i32 feature-vector classification (cross-PID anomaly,
48// e.g., "coolant + load + voltage" pattern matching a known
49// failure-mode signature)
50// - Kalman-filter smoother (better SNR for noisy PIDs)
51// - Multi-vehicle fleet aggregation (which units of an N-vehicle
52// fleet are first to show the drift pattern)
53// - Operator-tunable per-PID sigma threshold + r^2 floor
54
55import "nx_syscalls.nx"
56import "nx_telemetry.nx"
57
58// ===== Sealed verdict surface =================================================
59const NX_PM_OK: i64 = 0
60const NX_PM_PID_NOT_REGISTERED: i64 = 1
61const NX_PM_INSUFFICIENT_DATA: i64 = 2 // <MIN_SAMPLES samples; can't estimate
62const NX_PM_BUFFER_OVERFLOW: i64 = 3
63const NX_PM_BAD_INPUT: i64 = 4
64
65// ===== Sizing =================================================
66const NX_PM_MAX_PIDS: i64 = 32 // max tracked PIDs per vehicle
67const NX_PM_WINDOW_SIZE: i64 = 256 // samples per PID ring buffer
68const NX_PM_MIN_SAMPLES: i64 = 16 // min before stats are emitted
69
70// ===== Per-PID ring buffer =================================================
71//
72// Fixed-size circular buffer. When full, oldest sample is overwritten.
73// Stores (timestamp_ns, value) pairs as parallel arrays for cache
74// friendliness on future SIMD lanes.
75
76struct NxPidRing {
77 pid: i64 // OBD PID this ring covers
78 times_buf: *i64 // i64 ns since epoch
79 values_buf: *i64 // i64 raw event value
80 write_idx: i64 // next write position (mod cap)
81 n_filled: i64 // min(write_count, cap)
82 threshold_hi: i64 // OBD-spec'd upper bound (PID-specific)
83 threshold_lo: i64 // OBD-spec'd lower bound
84 valid: i64
85}
86
87// ===== Welford's online statistics =================================================
88//
89// Numerically stable rolling mean + variance per Welford (1962).
90// Recomputed each call across the entire window; future commit
91// adds incremental update with old-sample subtraction.
92
93struct NxPidStats {
94 n: i64 // sample count in this window
95 mean_x1k: i64 // mean * 1000 (Q3 fixed-point for sub-integer resolution)
96 var_x1k: i64 // variance * 1000
97 stddev_x1k: i64 // sqrt(var) * 1000; computed on demand
98 slope_x1k: i64 // linear-regression slope * 1000 (units/second)
99 r2_x100: i64 // r-squared * 100 (0..100; > 70 = strong fit)
100}
101
102// ===== Sigma anomaly verdict =================================================
103const NX_PM_ANOM_NONE: i64 = 0
104const NX_PM_ANOM_HIGH_3SIGMA: i64 = 1 // value > mean + 3*sigma
105const NX_PM_ANOM_LOW_3SIGMA: i64 = 2 // value < mean - 3*sigma
106const NX_PM_ANOM_DRIFT_UP: i64 = 3 // slope crossing threshold_hi within window
107const NX_PM_ANOM_DRIFT_DOWN: i64 = 4 // slope crossing threshold_lo within window
108
109// ===== Integer sqrt (for stddev computation) =================================================
110//
111// Newton's method; converges in ~6 iterations for typical Q3 values.
112
113func nx_pm_isqrt(x: i64) -> i64 {
114 if x <= 0 { return 0 }
115 if x == 1 { return 1 }
116 var r: i64 = x
117 var iters: i64 = 0
118 while iters < 16 {
119 let next: i64 = (r + (x / r)) / 2
120 if next >= r { return r }
121 r = next
122 iters = iters + 1
123 }
124 return r
125}
126
127// ===== Per-PID ring init =================================================
128
129func nx_pm_ring_init(ring: *NxPidRing, pid: i64,
130 times_buf: *i64, values_buf: *i64,
131 threshold_lo: i64, threshold_hi: i64) -> i64 {
132 if (ring as i64) == 0 { return 0 - NX_PM_BAD_INPUT }
133 if (times_buf as i64) == 0 { return 0 - NX_PM_BAD_INPUT }
134 if (values_buf as i64) == 0 { return 0 - NX_PM_BAD_INPUT }
135 ring.pid = pid
136 ring.times_buf = times_buf
137 ring.values_buf = values_buf
138 ring.write_idx = 0
139 ring.n_filled = 0
140 ring.threshold_lo = threshold_lo
141 ring.threshold_hi = threshold_hi
142 ring.valid = 1
143 return NX_PM_OK
144}
145
146// ===== Push one sample =================================================
147
148func nx_pm_ring_push(ring: *NxPidRing, ts_ns: i64, value: i64) -> i64 {
149 if ring.valid != 1 { return 0 - NX_PM_PID_NOT_REGISTERED }
150 let idx: i64 = ring.write_idx
151 ring.times_buf[idx] = ts_ns
152 ring.values_buf[idx] = value
153 ring.write_idx = (idx + 1) % NX_PM_WINDOW_SIZE
154 if ring.n_filled < NX_PM_WINDOW_SIZE {
155 ring.n_filled = ring.n_filled + 1
156 }
157 return NX_PM_OK
158}
159
160// ===== Welford-style stats over the entire ring =================================================
161//
162// V1: O(n) per recompute. Future: incremental update via push/pop.
163
164func nx_pm_compute_stats(ring: *NxPidRing, stats: *NxPidStats) -> i64 {
165 if ring.valid != 1 { return 0 - NX_PM_PID_NOT_REGISTERED }
166 if ring.n_filled < NX_PM_MIN_SAMPLES { return 0 - NX_PM_INSUFFICIENT_DATA }
167
168 // Welford's online mean + variance.
169 var n: i64 = 0
170 var mean: i64 = 0 // running mean * 1000 (Q3)
171 var m2: i64 = 0 // sum of squared deviations * 1000^2 (Q6)
172 var i: i64 = 0
173 while i < ring.n_filled {
174 let x: i64 = ring.values_buf[i] * 1000
175 n = n + 1
176 let delta: i64 = x - mean
177 mean = mean + (delta / n)
178 let delta2: i64 = x - mean
179 m2 = m2 + (delta * delta2)
180 i = i + 1
181 }
182
183 stats.n = n
184 stats.mean_x1k = mean
185 if n > 1 {
186 // Sample variance: m2 / (n-1). Result is in Q6 / 1000 = Q3.
187 let var_q6: i64 = m2 / (n - 1)
188 stats.var_x1k = var_q6 / 1000
189 stats.stddev_x1k = nx_pm_isqrt(stats.var_x1k * 1000) // sqrt(var*1000) gives Q3*1000... hmm
190 // Simpler: stddev_x1k = sqrt(var_x1k * 1000) in Q3 -- close enough for anomaly thresholds
191 }
192 if n <= 1 {
193 stats.var_x1k = 0
194 stats.stddev_x1k = 0
195 }
196
197 // Linear regression slope: cov(t, x) / var(t).
198 // Times are absolute ns; subtract first time for numerical stability.
199 let t0: i64 = ring.times_buf[0]
200 var sum_t: i64 = 0
201 var sum_x: i64 = 0
202 var sum_tx: i64 = 0
203 var sum_tt: i64 = 0
204 i = 0
205 while i < ring.n_filled {
206 let t_norm: i64 = (ring.times_buf[i] - t0) / 1000000 // ms since first sample
207 let x: i64 = ring.values_buf[i]
208 sum_t = sum_t + t_norm
209 sum_x = sum_x + x
210 sum_tx = sum_tx + (t_norm * x)
211 sum_tt = sum_tt + (t_norm * t_norm)
212 i = i + 1
213 }
214 let denom_t: i64 = (n * sum_tt) - (sum_t * sum_t)
215 if denom_t != 0 {
216 // slope = (n*sum_tx - sum_t*sum_x) / (n*sum_tt - sum_t^2)
217 // units: value-units per ms; * 1000 -> per-second; * 1000 again for Q3
218 let num: i64 = (n * sum_tx) - (sum_t * sum_x)
219 stats.slope_x1k = (num * 1000 * 1000) / denom_t
220 }
221 if denom_t == 0 {
222 stats.slope_x1k = 0
223 }
224
225 // r-squared = (cov / (sigma_t * sigma_x))^2. V1 simplification:
226 // r^2 proxy via |slope| relative to noise floor. Future commit
227 // computes proper Pearson r.
228 if stats.stddev_x1k > 0 {
229 let abs_slope: i64 = if stats.slope_x1k >= 0 then stats.slope_x1k else 0 - stats.slope_x1k
230 let ratio: i64 = (abs_slope * 100) / (stats.stddev_x1k + 1)
231 if ratio > 100 { stats.r2_x100 = 100 }
232 if ratio <= 100 { stats.r2_x100 = ratio }
233 }
234 if stats.stddev_x1k == 0 { stats.r2_x100 = 0 }
235
236 return NX_PM_OK
237}
238
239// ===== Sigma-based anomaly classifier =================================================
240//
241// Compares the latest sample to the rolling distribution. Returns
242// one of NX_PM_ANOM_* values.
243
244func nx_pm_classify_latest(ring: *NxPidRing, stats: *NxPidStats) -> i64 {
245 if ring.n_filled == 0 { return NX_PM_ANOM_NONE }
246 let latest_idx: i64 = (ring.write_idx + NX_PM_WINDOW_SIZE - 1) % NX_PM_WINDOW_SIZE
247 let latest: i64 = ring.values_buf[latest_idx] * 1000 // Q3 for comparison
248
249 let mean: i64 = stats.mean_x1k
250 let sd: i64 = stats.stddev_x1k
251
252 // 3-sigma bands.
253 let upper: i64 = mean + (3 * sd)
254 let lower: i64 = mean - (3 * sd)
255 if latest > upper { return NX_PM_ANOM_HIGH_3SIGMA }
256 if latest < lower { return NX_PM_ANOM_LOW_3SIGMA }
257
258 // Drift detection: slope * (seconds-to-threshold) crosses spec'd bound
259 // within the visible window. Requires r2 > 70 for confidence.
260 if stats.r2_x100 > 70 {
261 if stats.slope_x1k > 0 {
262 // upward drift; time to threshold_hi
263 let gap: i64 = (ring.threshold_hi * 1000) - latest // Q3
264 if gap > 0 {
265 // seconds_to_cross = gap / slope_x1k
266 let seconds_to_cross: i64 = gap / (stats.slope_x1k + 1)
267 // Flag if crossing predicted within window's worth of seconds
268 // (a generous extrapolation horizon for V1).
269 if seconds_to_cross < 86400 * 7 { // within 7 days
270 return NX_PM_ANOM_DRIFT_UP
271 }
272 }
273 }
274 if stats.slope_x1k < 0 {
275 let gap: i64 = latest - (ring.threshold_lo * 1000)
276 if gap > 0 {
277 let neg_slope: i64 = 0 - stats.slope_x1k
278 let seconds_to_cross: i64 = gap / (neg_slope + 1)
279 if seconds_to_cross < 86400 * 7 {
280 return NX_PM_ANOM_DRIFT_DOWN
281 }
282 }
283 }
284 }
285
286 return NX_PM_ANOM_NONE
287}
288
289// ===== Consumer struct =================================================
290//
291// Holds N per-PID rings + dispatches incoming OBD events to the right
292// ring, then triggers stats + classification.
293
294struct NxPredMaintConsumer {
295 rings: *NxPidRing // array of NX_PM_MAX_PIDS
296 n_registered: i64 // count of valid PIDs
297 telemetry: *NxTelemetryBus // optional; null = telemetry disabled
298 valid: i64
299}
300
301func nx_pm_consumer_init(c: *NxPredMaintConsumer, rings: *NxPidRing) -> i64 {
302 if (c as i64) == 0 { return 0 - NX_PM_BAD_INPUT }
303 c.rings = rings
304 c.n_registered = 0
305 c.telemetry = (0 as i64) as *NxTelemetryBus
306 c.valid = 1
307 return NX_PM_OK
308}
309
310// Optional: wire the consumer to a telemetry bus. When set,
311// on_event emits an NX_TM_KIND_ANOMALY event per non-NONE
312// anomaly verdict; sink dispatch + rate-limiting happens
313// inside the bus per nx_telemetry.nx contract.
314func nx_pm_consumer_set_telemetry_bus(c: *NxPredMaintConsumer,
315 bus: *NxTelemetryBus) -> i64 {
316 if c.valid != 1 { return 0 - NX_PM_BAD_INPUT }
317 c.telemetry = bus
318 return NX_PM_OK
319}
320
321// Register a PID with threshold bounds + caller-supplied backing buffers.
322func nx_pm_consumer_register(c: *NxPredMaintConsumer, pid: i64,
323 times_buf: *i64, values_buf: *i64,
324 threshold_lo: i64, threshold_hi: i64) -> i64 {
325 if c.valid != 1 { return 0 - NX_PM_BAD_INPUT }
326 if c.n_registered >= NX_PM_MAX_PIDS { return 0 - NX_PM_BUFFER_OVERFLOW }
327 let ring_ptr: *NxPidRing = (c.rings as i64 + c.n_registered * 64) as *NxPidRing
328 nx_pm_ring_init(ring_ptr, pid, times_buf, values_buf, threshold_lo, threshold_hi)
329 c.n_registered = c.n_registered + 1
330 return NX_PM_OK
331}
332
333// Look up the ring for a PID (linear scan; n <= MAX_PIDS).
334func nx_pm_consumer_find_ring(c: *NxPredMaintConsumer, pid: i64) -> *NxPidRing {
335 if c.valid != 1 { return (0 as i64) as *NxPidRing }
336 var i: i64 = 0
337 while i < c.n_registered {
338 let ring_ptr: *NxPidRing = (c.rings as i64 + i * 64) as *NxPidRing
339 if ring_ptr.pid == pid { return ring_ptr }
340 i = i + 1
341 }
342 return (0 as i64) as *NxPidRing
343}
344
345// ===== Top-level: process one OBD event =================================================
346//
347// Per INTEROPERABILITY_CHARTER §M2: substrate operates on
348// NxProtocolEvent. The consumer receives an event from the OBD shim,
349// looks up the corresponding ring, pushes the sample, recomputes
350// stats, classifies, and returns the anomaly verdict.
351
352func nx_pm_consumer_on_event(c: *NxPredMaintConsumer,
353 event_kind: i64, value: i64, ts_ns: i64,
354 out_stats: *NxPidStats) -> i64 {
355 if c.valid != 1 { return 0 - NX_PM_BAD_INPUT }
356 let ring: *NxPidRing = nx_pm_consumer_find_ring(c, event_kind)
357 if (ring as i64) == 0 { return 0 - NX_PM_PID_NOT_REGISTERED }
358
359 nx_pm_ring_push(ring, ts_ns, value)
360 let rc: i64 = nx_pm_compute_stats(ring, out_stats)
361 if rc != NX_PM_OK { return rc }
362 let anom: i64 = nx_pm_classify_latest(ring, out_stats)
363
364 // Emit telemetry if the consumer is wired to a bus AND the
365 // classifier flagged a non-NONE anomaly. V1 passes days_early=0
366 // for current 3-sigma flags; future commit returns days_early
367 // via classify_latest out-param for DRIFT_UP/DOWN cases (the
368 // value is already computed internally during classification).
369 if (c.telemetry as i64) != 0 {
370 if anom != NX_PM_ANOM_NONE {
371 nx_tm_emit_anomaly(c.telemetry, event_kind, anom, 0, ts_ns)
372 }
373 }
374
375 return anom
376}