nx_sketch_correlation.nx source
↩ module page · 211 lines · 6669 B
1// sketch_correlation.nx -- streaming Pearson correlation + linear regression.
2//
3// Bivariate streaming primitive. For paired observations (x, y) consumed
4// one at a time, computes:
5// r -- Pearson correlation coefficient, in PPM in [-1_000_000, 1_000_000]
6// slope -- least-squares regression slope m: y = m*x + b
7// intercept -- least-squares regression intercept b
8//
9// SIX ACCUMULATORS (integer-exact under overflow budget):
10// n = count
11// sum_x = Σ x_i
12// sum_y = Σ y_i
13// sum_xy = Σ x_i * y_i
14// sum_xx = Σ x_i²
15// sum_yy = Σ y_i²
16//
17// FORMULAS:
18// cov_xy = n * sum_xy - sum_x * sum_y
19// var_x = n * sum_xx - sum_x²
20// var_y = n * sum_yy - sum_y²
21// r = cov_xy / sqrt(var_x * var_y)
22// slope = cov_xy / var_x
23// intercept = (sum_y - slope * sum_x) / n
24//
25// OVERFLOW BUDGET:
26// sum_xx grows as n * x_max². For x_max=2^15, n_max=2^32 -> sum_xx
27// max ~ 2^62. And the cross-term computations need extra margin.
28// nx_corr_safe_p flags unsafe inputs.
29//
30// USE CASES:
31// - SRE: correlate metric series (CPU vs latency)
32// - finance: returns correlation
33// - science: regression on streamed measurements
34//
35// LOSSLESS-LANGUAGE DISCIPLINE: r returned in PPM with NX_ENV_ABS,
36// param_a = 0 (exact under budget). Production tier.
37
38// nx_safety_envelope:
39// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
40// sil_target: SIL1
41// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
42// verdict: NOT_YET_EVALUATED
43
44import "nx_syscalls.nx"
45import "nx_sketch_types.nx"
46
47const NX_CORR_VALUE_LIMIT: i64 = 32768 // |x|, |y| < 2^15 for safe sum_xx
48
49struct Correlation {
50 n: i64,
51 sum_x: i64,
52 sum_y: i64,
53 sum_xy: i64,
54 sum_xx: i64,
55 sum_yy: i64,
56 has_data: i64,
57}
58
59// === construction =================================================
60
61func nx_corr_alloc() -> *Correlation {
62 let raw: *u8 = sys_mmap(64)
63 let c: *Correlation = raw as *Correlation
64 c.n = 0
65 c.sum_x = 0
66 c.sum_y = 0
67 c.sum_xy = 0
68 c.sum_xx = 0
69 c.sum_yy = 0
70 c.has_data = 0
71 return c
72}
73
74// === isqrt =======================================================
75
76func nx_corr_isqrt(x: i64) -> i64 {
77 if x < 0 { return 0 }
78 if x == 0 { return 0 }
79 if x < 4 { return 1 }
80 var g: i64 = (x >> 1) + 1
81 var iter: i64 = 0
82 while iter < 64 {
83 let next_g: i64 = (g + x / g) / 2
84 if next_g >= g { iter = 64 }
85 if next_g < g {
86 g = next_g
87 iter = iter + 1
88 }
89 }
90 return g
91}
92
93// === overflow guard ==============================================
94
95func nx_corr_safe_p(x: i64, y: i64) -> i64 {
96 var ax: i64 = x
97 if ax < 0 { ax = -ax }
98 var ay: i64 = y
99 if ay < 0 { ay = -ay }
100 if ax >= NX_CORR_VALUE_LIMIT { return 0 }
101 if ay >= NX_CORR_VALUE_LIMIT { return 0 }
102 return 1
103}
104
105// === add ==========================================================
106
107func nx_corr_add(c: *Correlation, x: i64, y: i64) -> i64 {
108 if nx_corr_safe_p(x, y) == 0 { return -1 }
109 c.n = c.n + 1
110 c.sum_x = c.sum_x + x
111 c.sum_y = c.sum_y + y
112 c.sum_xy = c.sum_xy + x * y
113 c.sum_xx = c.sum_xx + x * x
114 c.sum_yy = c.sum_yy + y * y
115 c.has_data = 1
116 return 0
117}
118
119// === correlation coefficient (in PPM) ============================
120//
121// r = cov_xy / sqrt(var_x * var_y) where cov, var are the *centered*
122// values multiplied by n (so we avoid recomputing means).
123// Returns PPM in [-1_000_000, 1_000_000].
124
125func nx_corr_r_ppm(c: *Correlation) -> i64 {
126 if c.n < 2 { return 0 }
127 let cov_xy: i64 = c.n * c.sum_xy - c.sum_x * c.sum_y
128 let var_x: i64 = c.n * c.sum_xx - c.sum_x * c.sum_x
129 let var_y: i64 = c.n * c.sum_yy - c.sum_y * c.sum_y
130 if var_x <= 0 { return 0 }
131 if var_y <= 0 { return 0 }
132 // r = cov_xy / sqrt(var_x * var_y); compute in PPM.
133 // Avoid overflow in (var_x * var_y) by scaling cov upfront.
134 // r_ppm = (cov_xy * 1_000_000) / sqrt(var_x * var_y)
135 // Use staged isqrt: sqrt(var_x) * sqrt(var_y) (approximate but safe).
136 let sx: i64 = nx_corr_isqrt(var_x)
137 let sy: i64 = nx_corr_isqrt(var_y)
138 if sx == 0 { return 0 }
139 if sy == 0 { return 0 }
140 let denom: i64 = sx * sy
141 if denom == 0 { return 0 }
142 let r: i64 = (cov_xy * 1000000) / denom
143 // Clamp to [-1_000_000, 1_000_000].
144 if r > 1000000 { return 1000000 }
145 if r < -1000000 { return -1000000 }
146 return r
147}
148
149// === regression slope and intercept (in PPM) =====================
150//
151// slope_ppm = (n * sum_xy - sum_x * sum_y) / (n * sum_xx - sum_x²) * 1_000_000
152// intercept = mean_y - slope * mean_x (slope here is fractional;
153// scaled-PPM math: intercept = (sum_y * 1_000_000 - slope_ppm * sum_x) / (n * 1_000_000))
154
155func nx_corr_slope_ppm(c: *Correlation) -> i64 {
156 if c.n < 2 { return 0 }
157 let cov_xy: i64 = c.n * c.sum_xy - c.sum_x * c.sum_y
158 let var_x: i64 = c.n * c.sum_xx - c.sum_x * c.sum_x
159 if var_x <= 0 { return 0 }
160 return (cov_xy * 1000000) / var_x
161}
162
163func nx_corr_intercept(c: *Correlation) -> i64 {
164 if c.n < 2 { return 0 }
165 let slope_ppm: i64 = nx_corr_slope_ppm(c)
166 // intercept = mean_y - slope * mean_x
167 // = (sum_y - slope_ppm * sum_x / 1_000_000) / n
168 let slope_times_sum_x: i64 = (slope_ppm * c.sum_x) / 1000000
169 return (c.sum_y - slope_times_sum_x) / c.n
170}
171
172// === typed queries ===============================================
173
174func nx_corr_query_r(c: *Correlation) -> *ApproxI64 {
175 let r: i64 = nx_corr_r_ppm(c)
176 return nx_approx_new(r, NX_ENV_ABS, 1,
177 1000000000,
178 NX_MATURITY_PRODUCTION,
179 NX_ADV_HONEST)
180}
181
182func nx_corr_query_slope(c: *Correlation) -> *ApproxI64 {
183 let s: i64 = nx_corr_slope_ppm(c)
184 return nx_approx_new(s, NX_ENV_ABS, 1,
185 1000000000,
186 NX_MATURITY_PRODUCTION,
187 NX_ADV_HONEST)
188}
189
190// === merge ========================================================
191
192func nx_corr_merge(a: *Correlation, b: *Correlation) -> *Correlation {
193 let out: *Correlation = nx_corr_alloc()
194 out.n = a.n + b.n
195 out.sum_x = a.sum_x + b.sum_x
196 out.sum_y = a.sum_y + b.sum_y
197 out.sum_xy = a.sum_xy + b.sum_xy
198 out.sum_xx = a.sum_xx + b.sum_xx
199 out.sum_yy = a.sum_yy + b.sum_yy
200 if a.has_data == 1 { out.has_data = 1 }
201 if b.has_data == 1 { out.has_data = 1 }
202 return out
203}
204
205func nx_corr_memory_bytes(c: *Correlation) -> i64 {
206 return 64
207}
208
209func nx_corr_count(c: *Correlation) -> i64 {
210 return c.n
211}