sketch_reservoir.nx source
↩ module page · 182 lines · 6375 B
1// sketch_reservoir.nx -- Vitter reservoir sampling + quantile.
2//
3// Vitter 1985 "Random Sampling with a Reservoir." Algorithm R:
4// maintain a cap-element reservoir. For the i-th input value
5// (1-indexed):
6// - if i <= cap: items[i-1] = value
7// - else: draw j uniformly from [1, i]; if j <= cap,
8// items[j-1] = value (else discard)
9//
10// At end of stream, items[] holds a UNIFORMLY RANDOM sample of
11// the input. Per-quantile estimation: sort items, look up the
12// floor(p * n_items)-th value. Hoeffding bound on quantile from
13// k samples: error ~ sqrt(ln(1/delta) / (2k)) ~ 1/sqrt(k) at fixed
14// delta. We declare 1/sqrt(cap) as the rank_error envelope.
15//
16// OPENS TWO ROADMAP AXES SIMULTANEOUSLY:
17// - Sampling family (Vitter / VarOpt -- DataSketches ships both;
18// this is the first axis primitive)
19// - Quantile family (rank_error envelope; T-Digest / KLL ship
20// this same family with tighter bounds; reservoir is the
21// SIMPLEST sample-based quantile)
22//
23// LCG-based PRNG for deterministic reproducibility. Per the
24// lossless-language discipline (doc 20): randomized sketches MUST
25// declare their RNG state in the typed envelope so two runs with
26// the same seed produce bit-identical output.
27
28import "syscalls.nx"
29import "sketch_types.nx"
30
31const NX_RES_CAP_MIN: i64 = 4
32const NX_RES_CAP_MAX: i64 = 1000000
33
34// LCG constants -- glibc / Numerical Recipes lineage.
35const NX_RES_LCG_A: i64 = 1103515245
36const NX_RES_LCG_C: i64 = 12345
37const NX_RES_LCG_MOD: i64 = 0x7FFFFFFF
38
39struct Reservoir {
40 items: *i64,
41 cap: i64,
42 n_items: i64, // = min(total_seen, cap)
43 total_seen: i64,
44 rng_state: i64,
45}
46
47func nx_reservoir_alloc(cap: i64, seed: i64) -> *Reservoir {
48 if cap < NX_RES_CAP_MIN { return 0 as *Reservoir }
49 if cap > NX_RES_CAP_MAX { return 0 as *Reservoir }
50 let raw: *u8 = sys_mmap(40)
51 let r: *Reservoir = raw as *Reservoir
52 let items_raw: *u8 = sys_mmap(cap * 8)
53 r.items = items_raw as *i64
54 r.cap = cap
55 r.n_items = 0
56 r.total_seen = 0
57 r.rng_state = seed | 1 // never zero (LCG degenerate at 0)
58 return r
59}
60
61// LCG next: returns a 31-bit pseudorandom value. Mutates state.
62func nx_reservoir_rng_next(r: *Reservoir) -> i64 {
63 let next: i64 = ((r.rng_state * NX_RES_LCG_A) + NX_RES_LCG_C) & NX_RES_LCG_MOD
64 r.rng_state = next
65 return next
66}
67
68// Uniform integer in [0, n). n must be > 0.
69func nx_reservoir_rng_below(r: *Reservoir, n: i64) -> i64 {
70 let raw: i64 = nx_reservoir_rng_next(r)
71 return raw % n
72}
73
74func nx_reservoir_add(r: *Reservoir, value: i64) -> i64 {
75 r.total_seen = r.total_seen + 1
76 if r.n_items < r.cap {
77 r.items[r.n_items] = value
78 r.n_items = r.n_items + 1
79 return 0
80 }
81 // r.total_seen > cap: replace with probability cap / total_seen.
82 let j: i64 = nx_reservoir_rng_below(r, r.total_seen)
83 if j < r.cap {
84 r.items[j] = value
85 }
86 return 0
87}
88
89// === sort the sample (insertion sort; cap is small) ==============
90
91func nx_reservoir_sort(r: *Reservoir) -> i64 {
92 var i: i64 = 1
93 while i < r.n_items {
94 let cur: i64 = r.items[i]
95 var j: i64 = i - 1
96 var done: i64 = 0
97 while done == 0 {
98 if j < 0 { done = 1 }
99 if done == 0 {
100 let prev: i64 = r.items[j]
101 if prev <= cur { done = 1 }
102 if done == 0 {
103 r.items[j + 1] = prev
104 j = j - 1
105 }
106 }
107 }
108 r.items[j + 1] = cur
109 i = i + 1
110 }
111 return 0
112}
113
114// === quantile ====================================================
115//
116// p_milli: quantile in parts-per-thousand (500 = median, 990 = p99).
117// Returns the value at index floor(p_milli * n_items / 1000).
118// Caller responsible for calling nx_reservoir_sort first if the
119// reservoir hasn't been finalized -- we do it lazily here.
120
121func nx_reservoir_quantile(r: *Reservoir, p_milli: i64) -> i64 {
122 if r.n_items == 0 { return 0 }
123 nx_reservoir_sort(r)
124 var idx: i64 = (p_milli * r.n_items) / 1000
125 if idx < 0 { idx = 0 }
126 if idx >= r.n_items { idx = r.n_items - 1 }
127 return r.items[idx]
128}
129
130// === rank-of-value (inverse: how many items <= value, as fraction) ==
131
132func nx_reservoir_rank(r: *Reservoir, value: i64) -> i64 {
133 // Returns rank in parts-per-thousand. Sorts then binary-searches.
134 if r.n_items == 0 { return 0 }
135 nx_reservoir_sort(r)
136 // Count items <= value.
137 var lo: i64 = 0
138 var hi: i64 = r.n_items
139 while lo < hi {
140 let mid: i64 = (lo + hi) / 2
141 if r.items[mid] <= value {
142 lo = mid + 1
143 }
144 if r.items[mid] > value {
145 hi = mid
146 }
147 }
148 // lo now equals the number of items <= value.
149 return (lo * 1000) / r.n_items
150}
151
152// === error bound: rank_error ~ 1/sqrt(cap) ======================
153//
154// Hoeffding-style: P(|F_n - F| > eps) <= 2 * exp(-2 * n * eps^2)
155// For 95% confidence (delta=0.05): eps ~ sqrt(ln(40) / (2 * cap))
156// ~ 1.36 / sqrt(cap). We hardcode for typical cap values; for
157// other values, return a conservative bound based on cap class.
158
159func nx_reservoir_rank_error_ppb(cap: i64) -> i64 {
160 if cap <= 16 { return 340000000 } // 0.34
161 if cap <= 64 { return 170000000 } // 0.17
162 if cap <= 256 { return 85000000 } // 0.085
163 if cap <= 1024 { return 42500000 } // 0.0425
164 if cap <= 4096 { return 21250000 } // 0.02125
165 if cap <= 16384 { return 10625000 } // 0.0106
166 return 5000000 // 0.005 fallback
167}
168
169func nx_reservoir_query_quantile(r: *Reservoir, p_milli: i64) -> *ApproxI64 {
170 let v: i64 = nx_reservoir_quantile(r, p_milli)
171 return nx_approx_new(v, NX_ENV_RANK_ERROR,
172 nx_reservoir_rank_error_ppb(r.cap),
173 950000000, // 95% confidence
174 NX_MATURITY_REFERENCE_IMPL,
175 NX_ADV_HONEST)
176}
177
178// === memory introspection =======================================
179
180func nx_reservoir_memory_bytes(r: *Reservoir) -> i64 {
181 return 40 + r.cap * 8
182}