nx_sketch_varopt.nx source
↩ module page · 245 lines · 8507 B
1// sketch_varopt.nx -- Weighted reservoir sampling (Efraimidis-Spirakis 2006 / Cohen 2011).
2//
3// Given a stream of (item, weight) arrivals, maintain a reservoir
4// of size k such that the inclusion probability is proportional
5// to weight. Heavy items more likely to be retained than light.
6//
7// COMPLETES THE SAMPLING FAMILY:
8// - Reservoir (Vitter 1985): uniform sampling, equal weights.
9// - VarOpt (here): weighted sampling.
10//
11// ALGORITHM (A-ExpJ variant of A-Res):
12// For each (item, w): generate key = u^(1/w) where u ~ Uniform(0,1).
13// Keep top-k items by key (descending). Higher w concentrates the
14// key near 1; lower w near 0. Inclusion probability is provably
15// proportional to w under reasonable conditions.
16//
17// INTEGER-FIXED-POINT APPROXIMATION (because C anchor lacks f64,
18// but we want the sibling to be Wheeler-comparable against C):
19// - u_31 = LCG draw in (0, 2^31)
20// - approx_log2_u = bitlen(u_31) - 31 // in [-30, 0]
21// - approx_log_key = (approx_log2_u * 1_000_000) / w
22// - key = approx_log_key
23// Higher weight pulls key closer to 0; lower weight pulls more
24// negative. Sort descending by key; top-k by key are retained.
25//
26// PROPERTIES UNDER THE APPROXIMATION:
27// - Heavy items (w >> 1) get keys near 0 (large in descending sort).
28// - Light items (w = 1) get keys uniformly distributed in
29// [-30_000_000, 0], breaking ties uniformly.
30// - Quality degrades for very large weight variance. v2 will swap
31// in a finer-grained log approximation.
32//
33// LOSSLESS-LANGUAGE DISCIPLINE: per-item retain probability declared
34// in the typed envelope as conf_ppb proportional to weight_total /
35// weight_sum_sampled.
36
37// nx_safety_envelope:
38// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
39// sil_target: SIL1
40// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
41// verdict: NOT_YET_EVALUATED
42
43import "nx_syscalls.nx"
44import "nx_sketch_types.nx"
45
46const NX_VOPT_K_MIN: i64 = 4
47const NX_VOPT_K_MAX: i64 = 100000
48
49const NX_VOPT_LCG_A: i64 = 1103515245
50const NX_VOPT_LCG_C: i64 = 12345
51const NX_VOPT_LCG_MOD: i64 = 0x7FFFFFFF
52
53struct VoptEntry {
54 item: i64,
55 weight: i64,
56 key: i64, // negative; higher (closer to 0) = more retain-worthy
57}
58
59struct VarOpt {
60 entries: *VoptEntry, // sorted descending by key
61 k: i64,
62 n_items: i64, // = min(seen, k)
63 total_seen: i64,
64 total_weight: i64, // cumulative weight observed
65 rng_state: i64,
66}
67
68// === construction =================================================
69
70func nx_varopt_alloc(k: i64, seed: i64) -> *VarOpt {
71 if k < NX_VOPT_K_MIN { return 0 as *VarOpt }
72 if k > NX_VOPT_K_MAX { return 0 as *VarOpt }
73 let raw: *u8 = sys_mmap(48)
74 let v: *VarOpt = raw as *VarOpt
75 let ent_raw: *u8 = sys_mmap(k * 24)
76 v.entries = ent_raw as *VoptEntry
77 v.k = k
78 v.n_items = 0
79 v.total_seen = 0
80 v.total_weight = 0
81 v.rng_state = seed | 1
82 return v
83}
84
85func nx_varopt_rng_next(v: *VarOpt) -> i64 {
86 let next: i64 = ((v.rng_state * NX_VOPT_LCG_A) + NX_VOPT_LCG_C) & NX_VOPT_LCG_MOD
87 v.rng_state = next
88 return next
89}
90
91// === bit-length of an i64 (highest set bit position + 1) ==========
92//
93// For x in [1, 2^31), returns value in [1, 31].
94
95func nx_varopt_bitlen(x: i64) -> i64 {
96 if x <= 0 { return 0 }
97 var n: i64 = 0
98 var t: i64 = x
99 while t > 0 {
100 t = t >> 1
101 n = n + 1
102 }
103 return n
104}
105
106// === key generation ===============================================
107//
108// approx_log2_u = bitlen(u_31) - 31 ∈ [-30, 0]
109// approx_log_key = (approx_log2_u * 1_000_000) / weight
110
111func nx_varopt_key(v: *VarOpt, weight: i64) -> i64 {
112 if weight <= 0 { return -2147483647 } // refuse zero/negative weight
113 let u: i64 = nx_varopt_rng_next(v)
114 if u <= 0 { return -2147483647 }
115 let lg: i64 = nx_varopt_bitlen(u) // 1..31
116 let approx_log2: i64 = lg - 31 // [-30, 0]
117 return (approx_log2 * 1000000) / weight
118}
119
120// === entry access =================================================
121
122func nx_varopt_entry_at(v: *VarOpt, i: i64) -> *VoptEntry {
123 return (v.entries as i64 + i * 24) as *VoptEntry
124}
125
126// === add ===========================================================
127//
128// New (item, weight) -> generate key -> insert into sorted-by-key
129// descending list (insertion sort). If list at cap, replace last
130// (smallest-key) entry if new key is larger.
131
132func nx_varopt_add(v: *VarOpt, item: i64, weight: i64) -> i64 {
133 if weight <= 0 { return -1 }
134 v.total_seen = v.total_seen + 1
135 v.total_weight = v.total_weight + weight
136 let key: i64 = nx_varopt_key(v, weight)
137
138 if v.n_items < v.k {
139 // Insert in sorted position (descending by key).
140 var pos: i64 = v.n_items
141 var done: i64 = 0
142 while done == 0 {
143 if pos == 0 { done = 1 }
144 if done == 0 {
145 let prev: *VoptEntry = nx_varopt_entry_at(v, pos - 1)
146 if prev.key >= key { done = 1 }
147 if done == 0 {
148 let dst: *VoptEntry = nx_varopt_entry_at(v, pos)
149 dst.item = prev.item
150 dst.weight = prev.weight
151 dst.key = prev.key
152 pos = pos - 1
153 }
154 }
155 }
156 let slot: *VoptEntry = nx_varopt_entry_at(v, pos)
157 slot.item = item
158 slot.weight = weight
159 slot.key = key
160 v.n_items = v.n_items + 1
161 return 0
162 }
163 // At cap: only replace if new key > current minimum (entries[k-1]).
164 let last: *VoptEntry = nx_varopt_entry_at(v, v.k - 1)
165 if key <= last.key { return 0 }
166 // Drop the last; insert new in sorted position.
167 var pos: i64 = v.k - 1
168 var done: i64 = 0
169 while done == 0 {
170 if pos == 0 { done = 1 }
171 if done == 0 {
172 let prev: *VoptEntry = nx_varopt_entry_at(v, pos - 1)
173 if prev.key >= key { done = 1 }
174 if done == 0 {
175 let dst: *VoptEntry = nx_varopt_entry_at(v, pos)
176 dst.item = prev.item
177 dst.weight = prev.weight
178 dst.key = prev.key
179 pos = pos - 1
180 }
181 }
182 }
183 let slot: *VoptEntry = nx_varopt_entry_at(v, pos)
184 slot.item = item
185 slot.weight = weight
186 slot.key = key
187 return 0
188}
189
190// === queries =======================================================
191//
192// Compute the weighted-mean of an attribute the caller wishes to
193// estimate. In this v1 sketch we expose the raw item array; callers
194// using i64 items can compute custom aggregates. Estimator:
195// E[f(stream)] ≈ (1/n_items) * Σ f(item_i) * (total_weight / weight_i)
196// where the (total_weight / weight_i) is the importance-sampling
197// inverse-probability weight.
198
199func nx_varopt_n_items(v: *VarOpt) -> i64 {
200 return v.n_items
201}
202
203func nx_varopt_total_seen(v: *VarOpt) -> i64 {
204 return v.total_seen
205}
206
207func nx_varopt_total_weight(v: *VarOpt) -> i64 {
208 return v.total_weight
209}
210
211// Weighted mean of items in the sample, with inverse-probability
212// reweighting. Returns sum(item_i * total_weight / weight_i) / n_items.
213// (For a constant-attribute estimator; users compute their own when
214// item is structured.)
215func nx_varopt_weighted_mean(v: *VarOpt) -> i64 {
216 if v.n_items == 0 { return 0 }
217 var sum: i64 = 0
218 var i: i64 = 0
219 while i < v.n_items {
220 let e: *VoptEntry = nx_varopt_entry_at(v, i)
221 sum = sum + (e.item * v.total_weight) / e.weight
222 i = i + 1
223 }
224 return sum / v.n_items
225}
226
227// === envelope ======================================================
228//
229// rel_stddev for the importance-weighted estimator. Cohen 2011
230// shows VarOpt achieves variance ~ N · Σ w_i^2 / k. We declare a
231// conservative 1/sqrt(k) bound that's accurate for low-variance
232// weight distributions; high-variance weights degrade quality.
233
234func nx_varopt_query_mean(v: *VarOpt) -> *ApproxI64 {
235 let mean: i64 = nx_varopt_weighted_mean(v)
236 let stddev_ppb: i64 = (1000000000 / nx_varopt_bitlen(v.k)) / 8 // rough 1/sqrt(k)
237 return nx_approx_new(mean, NX_ENV_REL_STDDEV, stddev_ppb,
238 682700000, // 1-sigma
239 NX_MATURITY_REFERENCE_IMPL,
240 NX_ADV_HONEST)
241}
242
243func nx_varopt_memory_bytes(v: *VarOpt) -> i64 {
244 return 48 + v.k * 24
245}