sketch_reqsketch.nx source
↩ module page · 380 lines · 12726 B
1// sketch_reqsketch.nx -- relative-error streaming quantile sketch.
2//
3// Cormode-Karnin-Liberty-Thaler-Vesely 2021, "Relative Error
4// Streaming Quantiles" (FOCS / VLDB). ReqSketch tightens KLL on
5// tail queries: where KLL bounds |rank_est - rank_true| <= eps * N
6// (additive in N), ReqSketch bounds |rank_est - rank_true| <=
7// eps * min(rank_true, N - rank_true) (multiplicative in the rank's
8// distance from the nearer extreme). Tail quantiles (p99 / p99.9
9// / p0.1) get strictly tighter guarantees as N grows.
10//
11// COMPLEMENTS KLL + T-DIGEST:
12// - KLL: provable uniform additive bound; loose at tails.
13// - T-Digest: empirically tail-tight; no global rank guarantee.
14// - ReqSketch: provable multiplicative tail bound; strictly
15// stronger than KLL for tail queries at large N.
16//
17// CONSTRUCTION (simplified reference impl, hra = "high-rank-
18// accurate" upper-tail-tight variant -- mirror with hra=0 for
19// lower-tail-tight):
20//
21// Compactor cascade like KLL, but the compact() rule preserves
22// the upper section of each compactor exactly (no down-sampling)
23// while KLL-coin-flipping the lower section. Section size shrinks
24// linearly with level -- at level h, preserve_count_h = max(1, k -
25// 2*h), so at high levels almost all items get evicted, while at
26// low levels half stay in place. Items preserved at level h carry
27// weight 2^h same as KLL; coin-flipped survivors get weight 2^(h+1)
28// upon promotion.
29//
30// SIMPLIFICATION VS PAPER (honest scope):
31// Paper uses numSections in {3, 6, 12, ...} growing geometrically
32// with N and section_size = sqrt(2) * k_base * 2^(-h/2). We use
33// numSections=2 fixed and section_size = max(1, k/2 - h) linear
34// shrinkage. Same asymptotic shape (eps ~ 1/sqrt(k)), slightly
35// looser constant. v2 will swap in the geometric schedule.
36//
37// LOSSLESS-LANGUAGE DISCIPLINE (doc 20):
38// Query returns ApproxI64 with envelope_kind = NX_ENV_REL_RANK_ERROR,
39// param_a = eps_ppb (relative rank error * 1e9), conf = 0.95,
40// maturity = ReferenceImpl, adv = Honest.
41
42import "syscalls.nx"
43import "sketch_types.nx"
44
45const NX_REQ_K_MIN: i64 = 8
46const NX_REQ_K_MAX: i64 = 1024
47const NX_REQ_MAX_LEVELS: i64 = 24
48
49// LCG constants (matches KLL for cross-comparison determinism).
50const NX_REQ_LCG_A: i64 = 1103515245
51const NX_REQ_LCG_C: i64 = 12345
52const NX_REQ_LCG_MOD: i64 = 0x7FFFFFFF
53
54// Compactor capacity = 2 * k (sorted on compact; preserve_count
55// stays, the rest get coin-flipped to next level).
56
57struct Req {
58 k: i64,
59 n_levels: i64,
60 level_items: *i64, // contiguous: max_levels * (2 * k_max) entries
61 level_count: *i64, // current count per level
62 total_items: i64,
63 min_val: i64,
64 max_val: i64,
65 rng_state: i64,
66 hra: i64, // 1 = upper-tail-tight, 0 = lower-tail-tight
67 view_vals: *i64, // bits-up: hoisted from nx_req_quantile
68 view_wts: *i64,
69}
70
71// === construction =================================================
72
73func nx_req_alloc(k: i64, hra: i64, seed: i64) -> *Req {
74 if k < NX_REQ_K_MIN { return 0 as *Req }
75 if k > NX_REQ_K_MAX { return 0 as *Req }
76 let raw: *u8 = sys_mmap(88)
77 let s: *Req = raw as *Req
78 s.k = k
79 s.n_levels = 1
80 let items_bytes: i64 = NX_REQ_MAX_LEVELS * (2 * k) * 8
81 let items_raw: *u8 = sys_mmap(items_bytes)
82 s.level_items = items_raw as *i64
83 let count_raw: *u8 = sys_mmap(NX_REQ_MAX_LEVELS * 8)
84 s.level_count = count_raw as *i64
85 let view_cap_bytes: i64 = NX_REQ_MAX_LEVELS * 2 * k * 8
86 s.view_vals = sys_mmap(view_cap_bytes) as *i64
87 s.view_wts = sys_mmap(view_cap_bytes) as *i64
88 var i: i64 = 0
89 while i < NX_REQ_MAX_LEVELS {
90 s.level_count[i] = 0
91 i = i + 1
92 }
93 s.total_items = 0
94 s.min_val = 0x7FFFFFFFFFFFFFFF
95 s.max_val = -1 - 0x7FFFFFFFFFFFFFFF
96 s.rng_state = seed | 1
97 var hra_norm: i64 = 1
98 if hra == 0 { hra_norm = 0 }
99 s.hra = hra_norm
100 return s
101}
102
103func nx_req_rng_next(s: *Req) -> i64 {
104 let next: i64 = ((s.rng_state * NX_REQ_LCG_A) + NX_REQ_LCG_C) & NX_REQ_LCG_MOD
105 s.rng_state = next
106 return next
107}
108
109// Level h occupies a 2*k-sized window in level_items.
110func nx_req_level_addr(s: *Req, level: i64) -> *i64 {
111 let stride: i64 = 2 * s.k
112 return (s.level_items as i64 + level * stride * 8) as *i64
113}
114
115// preserve_count_h = max(1, k - 2*h). When h >= k/2, preserve = 1
116// (almost everything gets evicted upward; this is what gives the
117// relative-rank-error property at high levels).
118func nx_req_preserve_count(s: *Req, level: i64) -> i64 {
119 let raw: i64 = s.k - 2 * level
120 if raw < 1 { return 1 }
121 return raw
122}
123
124// === sort one level's buffer ======================================
125
126func nx_req_sort_level(s: *Req, level: i64) -> i64 {
127 let count: i64 = s.level_count[level]
128 let buf: *i64 = nx_req_level_addr(s, level)
129 var i: i64 = 1
130 while i < count {
131 let cur: i64 = buf[i]
132 var j: i64 = i - 1
133 var done: i64 = 0
134 while done == 0 {
135 if j < 0 { done = 1 }
136 if done == 0 {
137 let prev: i64 = buf[j]
138 if prev <= cur { done = 1 }
139 if done == 0 {
140 buf[j + 1] = prev
141 j = j - 1
142 }
143 }
144 }
145 buf[j + 1] = cur
146 i = i + 1
147 }
148 return 0
149}
150
151// === compact ======================================================
152//
153// On full (count >= 2*k):
154// 1. Sort.
155// 2. preserve = preserve_count_h items. If hra=1, preserve the
156// TOP preserve items (indices [count-preserve, count)); these
157// stay at this level. Coin-flip-compact the BOTTOM (count -
158// preserve) items (indices [0, count-preserve)); every other
159// one (at random parity) promotes to level h+1.
160// 3. If hra=0, mirror: preserve bottom, compact top.
161
162func nx_req_compact(s: *Req, level: i64) -> i64 {
163 if level >= NX_REQ_MAX_LEVELS - 1 {
164 // No room to promote -- drop level (rare, only at extreme N).
165 s.level_count[level] = 0
166 return 0
167 }
168 nx_req_sort_level(s, level)
169 let count: i64 = s.level_count[level]
170 let buf: *i64 = nx_req_level_addr(s, level)
171 let preserve: i64 = nx_req_preserve_count(s, level)
172 var compact_start: i64 = 0
173 var compact_end: i64 = 0
174 var preserve_start: i64 = 0
175 var preserve_end: i64 = 0
176 if s.hra == 1 {
177 compact_start = 0
178 compact_end = count - preserve
179 preserve_start = count - preserve
180 preserve_end = count
181 }
182 if s.hra == 0 {
183 preserve_start = 0
184 preserve_end = preserve
185 compact_start = preserve
186 compact_end = count
187 }
188 // Coin-flip parity for promotion.
189 let pick_odd: i64 = nx_req_rng_next(s) & 1
190 let next_level: i64 = level + 1
191 let next_buf: *i64 = nx_req_level_addr(s, next_level)
192 let next_count_before: i64 = s.level_count[next_level]
193 let next_cap: i64 = 2 * s.k
194 var i: i64 = compact_start + pick_odd
195 var promoted: i64 = 0
196 while i < compact_end {
197 let dst_idx: i64 = next_count_before + promoted
198 if dst_idx < next_cap {
199 next_buf[dst_idx] = buf[i]
200 promoted = promoted + 1
201 }
202 i = i + 2
203 }
204 s.level_count[next_level] = next_count_before + promoted
205 // Collapse preserved items down to start of buf.
206 if s.hra == 1 {
207 var j: i64 = 0
208 while j < preserve {
209 buf[j] = buf[preserve_start + j]
210 j = j + 1
211 }
212 s.level_count[level] = preserve
213 }
214 if s.hra == 0 {
215 // Preserved already at indices [0, preserve), nothing to move.
216 s.level_count[level] = preserve
217 }
218 if next_level >= s.n_levels {
219 s.n_levels = next_level + 1
220 }
221 // Cascade if next level is now full.
222 if s.level_count[next_level] >= next_cap {
223 nx_req_compact(s, next_level)
224 }
225 return 0
226}
227
228// === add ==========================================================
229
230func nx_req_add(s: *Req, value: i64) -> i64 {
231 if value < s.min_val { s.min_val = value }
232 if value > s.max_val { s.max_val = value }
233 s.total_items = s.total_items + 1
234 let buf: *i64 = nx_req_level_addr(s, 0)
235 let count: i64 = s.level_count[0]
236 let cap: i64 = 2 * s.k
237 buf[count] = value
238 s.level_count[0] = count + 1
239 if s.level_count[0] >= cap {
240 nx_req_compact(s, 0)
241 }
242 return 0
243}
244
245// === total weight ================================================
246
247func nx_req_total_weight(s: *Req) -> i64 {
248 var total: i64 = 0
249 var h: i64 = 0
250 while h < s.n_levels {
251 total = total + s.level_count[h] * (1 << h)
252 h = h + 1
253 }
254 return total
255}
256
257// === quantile =====================================================
258//
259// Materialise (value, weight) view across all levels, sort by value,
260// scan accumulating cumulative weight, return value at target rank.
261
262func nx_req_quantile(s: *Req, p_milli: i64) -> i64 {
263 if s.total_items == 0 { return 0 }
264 let total_weight: i64 = nx_req_total_weight(s)
265 if total_weight == 0 { return 0 }
266 // Bits-up: scratch hoisted to struct (was per-query sys_mmap).
267 let view_vals: *i64 = s.view_vals
268 let view_wts: *i64 = s.view_wts
269 var n: i64 = 0
270 var h: i64 = 0
271 while h < s.n_levels {
272 let buf: *i64 = nx_req_level_addr(s, h)
273 let cnt: i64 = s.level_count[h]
274 let w: i64 = 1 << h
275 var i: i64 = 0
276 while i < cnt {
277 view_vals[n] = buf[i]
278 view_wts[n] = w
279 n = n + 1
280 i = i + 1
281 }
282 h = h + 1
283 }
284 // Insertion sort the view.
285 var i: i64 = 1
286 while i < n {
287 let cur_v: i64 = view_vals[i]
288 let cur_w: i64 = view_wts[i]
289 var j: i64 = i - 1
290 var done: i64 = 0
291 while done == 0 {
292 if j < 0 { done = 1 }
293 if done == 0 {
294 if view_vals[j] <= cur_v { done = 1 }
295 if done == 0 {
296 view_vals[j + 1] = view_vals[j]
297 view_wts[j + 1] = view_wts[j]
298 j = j - 1
299 }
300 }
301 }
302 view_vals[j + 1] = cur_v
303 view_wts[j + 1] = cur_w
304 i = i + 1
305 }
306 let target: i64 = (p_milli * total_weight) / 1000
307 var cum: i64 = 0
308 var k: i64 = 0
309 while k < n {
310 cum = cum + view_wts[k]
311 if cum >= target { return view_vals[k] }
312 k = k + 1
313 }
314 return view_vals[n - 1]
315}
316
317// === rank-of-value ===============================================
318
319func nx_req_rank(s: *Req, value: i64) -> i64 {
320 if s.total_items == 0 { return 0 }
321 if value < s.min_val { return 0 }
322 if value >= s.max_val { return 1000 }
323 let total_weight: i64 = nx_req_total_weight(s)
324 if total_weight == 0 { return 0 }
325 var cum: i64 = 0
326 var h: i64 = 0
327 while h < s.n_levels {
328 let buf: *i64 = nx_req_level_addr(s, h)
329 let cnt: i64 = s.level_count[h]
330 let w: i64 = 1 << h
331 var i: i64 = 0
332 while i < cnt {
333 if buf[i] <= value {
334 cum = cum + w
335 }
336 i = i + 1
337 }
338 h = h + 1
339 }
340 return (cum * 1000) / total_weight
341}
342
343// === relative-rank-error bound ===================================
344//
345// Reference-impl bound: eps ~ 2.0 / sqrt(k) at conf 0.95. This is
346// the relative bound vs the nearer extreme (rank_true or N -
347// rank_true), so for p=0.99 with N=10000 the effective absolute
348// error is eps * 100 = ~2 ranks at k=128, vs KLL's eps_additive *
349// N = ~1200 ranks at k=128 -- 600x tighter on tail.
350//
351// Constants tabulated for representative k; consistent with KLL's
352// 1.36/sqrt(k) shape, looser constant accounts for the simplified
353// section schedule.
354
355func nx_req_rel_rank_error_ppb(k: i64) -> i64 {
356 if k <= 8 { return 707000000 } // 0.707 (2 / sqrt(8))
357 if k <= 32 { return 354000000 } // 0.354
358 if k <= 128 { return 177000000 } // 0.177
359 if k <= 512 { return 88000000 } // 0.088
360 return 44000000 // 0.044 for k > 512
361}
362
363func nx_req_query_quantile(s: *Req, p_milli: i64) -> *ApproxI64 {
364 let v: i64 = nx_req_quantile(s, p_milli)
365 return nx_approx_new(v, NX_ENV_REL_RANK_ERROR,
366 nx_req_rel_rank_error_ppb(s.k),
367 950000000,
368 NX_MATURITY_REFERENCE_IMPL,
369 NX_ADV_HONEST)
370}
371
372// === introspection ================================================
373
374func nx_req_memory_bytes(s: *Req) -> i64 {
375 return 72 + NX_REQ_MAX_LEVELS * (2 * s.k) * 8 + NX_REQ_MAX_LEVELS * 8
376}
377
378func nx_req_levels_used(s: *Req) -> i64 {
379 return s.n_levels
380}