sketch_theta.nx source
↩ module page · 371 lines · 12861 B
1// sketch_theta.nx -- Theta sketch (Beyer/Haas 2014).
2//
3// Extends KMV with EXPLICIT theta and SET INTERSECTION + DIFFERENCE.
4// Where KMV's theta is always implied by kth_smallest, Theta carries
5// theta as a first-class field so intersected sketches (which have
6// smaller-than-natural theta) remain valid for further composition.
7//
8// Theta semantics:
9// theta ∈ (0, HASH_MAX]: the bound below which hashes are sampled.
10// Initially HASH_MAX (entire domain). After K hashes seen,
11// theta = kth_smallest. After intersection, theta = min of inputs'.
12//
13// Cardinality estimator (theta-corrected, unbiased per Beyer 2007):
14// est = |hashes < theta| * HASH_MAX / theta
15// For a freshly-filled KMV this reduces to (K-1)*HASH_MAX/kth
16// which matches the KMV formula.
17//
18// COMPLETES THE SET-OP TRIO:
19// - union: runtime/sketch_kmv.nx::nx_kmv_union
20// - intersection: nx_theta_intersect (THIS FILE) -- new capability
21// - difference: nx_theta_difference -- new capability
22//
23// This is the family DataSketches ships as Theta + Tuple sketches.
24// We ship Theta first; Tuple (attaching auxiliary data to each hash)
25// is queued for v2 once the core sketches stabilize.
26//
27// Reference: DataSketches Theta paper -- arXiv:1306.5664
28// "Beyond HyperLogLog: A Survey of Sketches for Cardinality Estimation"
29
30import "syscalls.nx"
31import "murmur3.nx"
32import "sketch_types.nx"
33import "sketch_kmv.nx"
34
35// 2^32 hash domain.
36const NX_THETA_HASH_MAX: i64 = 4294967296
37
38struct ThetaSketch {
39 values: *i64, // sorted ascending; max at n_items-1
40 k: i64, // capacity
41 n_items: i64, // = |hashes stored|, <= k
42 theta: i64, // sample-rate bound, in [1, HASH_MAX]
43 seed: i64,
44}
45
46// === construction =================================================
47
48func nx_theta_alloc(k: i64, seed: i64) -> *ThetaSketch {
49 if k < 16 { return 0 as *ThetaSketch }
50 if k > 65536 { return 0 as *ThetaSketch }
51 let raw: *u8 = sys_mmap(48)
52 let t: *ThetaSketch = raw as *ThetaSketch
53 t.k = k
54 t.n_items = 0
55 t.theta = NX_THETA_HASH_MAX
56 let vals_raw: *u8 = sys_mmap(k * 8)
57 t.values = vals_raw as *i64
58 t.seed = seed
59 return t
60}
61
62// === wrap an existing KMV as a Theta sketch =====================
63//
64// Used so users can promote a KMV (which only does union) into a
65// Theta (which does intersection) without re-streaming.
66
67func nx_theta_from_kmv(kmv: *Kmv) -> *ThetaSketch {
68 let t: *ThetaSketch = nx_theta_alloc(kmv.k, kmv.seed)
69 var i: i64 = 0
70 while i < kmv.n_items {
71 t.values[i] = kmv.values[i]
72 i = i + 1
73 }
74 t.n_items = kmv.n_items
75 // theta is kth_smallest when full, else HASH_MAX.
76 if kmv.n_items == kmv.k {
77 t.theta = kmv.values[kmv.k - 1]
78 }
79 return t
80}
81
82// === add ========================================================
83
84func nx_theta_add(t: *ThetaSketch, key: *u8, len: i64) -> i64 {
85 let raw_h: i64 = murmur3_32(t.seed, key, len) & 0xFFFFFFFF
86 var h: i64 = raw_h
87 if h == 0 { h = 1 }
88 // If h >= theta, item is outside the sample.
89 if h >= t.theta { return 0 }
90
91 // Binary search.
92 var lo: i64 = 0
93 var hi: i64 = t.n_items
94 while lo < hi {
95 let mid: i64 = (lo + hi) / 2
96 let v: i64 = t.values[mid]
97 if v == h { return 0 }
98 if v < h { lo = mid + 1 }
99 if v > h { hi = mid }
100 }
101 if t.n_items < t.k {
102 var j: i64 = t.n_items
103 while j > lo {
104 t.values[j] = t.values[j - 1]
105 j = j - 1
106 }
107 t.values[lo] = h
108 t.n_items = t.n_items + 1
109 // First time we hit cap: theta becomes kth_smallest.
110 if t.n_items == t.k {
111 t.theta = t.values[t.k - 1]
112 }
113 return 0
114 }
115 // Full: replace max if h < max.
116 let cur_max: i64 = t.values[t.k - 1]
117 if h >= cur_max { return 0 }
118 var j2: i64 = t.k - 1
119 while j2 > lo {
120 t.values[j2] = t.values[j2 - 1]
121 j2 = j2 - 1
122 }
123 t.values[lo] = h
124 // theta tightens to new kth_smallest.
125 t.theta = t.values[t.k - 1]
126 return 0
127}
128
129// === cardinality (theta-corrected) ==============================
130//
131// est = |hashes < theta| * HASH_MAX / theta
132// |hashes < theta| = n_items minus the count of "border" hashes
133// equal to theta itself. Since insertions reject h >= theta,
134// all stored hashes are strictly < theta IF we apply the strict
135// rule: but in KMV (and standard Theta) the kth_smallest IS
136// stored and theta == kth_smallest, so n_items includes one
137// element AT theta. Use the canonical unbiased estimator:
138// est = (n_items - 1) * HASH_MAX / theta when n_items > 0
139// and theta < HASH_MAX. Otherwise est = n_items.
140
141func nx_theta_estimate(t: *ThetaSketch) -> i64 {
142 if t.n_items == 0 { return 0 }
143 if t.theta == NX_THETA_HASH_MAX {
144 // No down-sampling yet; exact count.
145 return t.n_items
146 }
147 // Count hashes STRICTLY below theta (exclude border).
148 var below: i64 = 0
149 var i: i64 = 0
150 while i < t.n_items {
151 if t.values[i] < t.theta { below = below + 1 }
152 i = i + 1
153 }
154 if below == 0 { return 0 }
155 return (below * NX_THETA_HASH_MAX) / t.theta
156}
157
158// === typed query ================================================
159
160func nx_theta_stddev_rel_ppb(k: i64) -> i64 {
161 if k <= 16 { return 258000000 }
162 if k <= 64 { return 126000000 }
163 if k <= 256 { return 62700000 }
164 if k <= 1024 { return 31300000 }
165 if k <= 4096 { return 15600000 }
166 if k <= 16384 { return 7810000 }
167 return 3900000
168}
169
170func nx_theta_query(t: *ThetaSketch) -> *ApproxI64 {
171 let est: i64 = nx_theta_estimate(t)
172 return nx_approx_new(est, NX_ENV_REL_STDDEV,
173 nx_theta_stddev_rel_ppb(t.k),
174 682700000,
175 NX_MATURITY_REFERENCE_IMPL,
176 NX_ADV_HONEST)
177}
178
179// === union (theta-aware) ========================================
180//
181// theta_out = min(theta_a, theta_b)
182// hashes_out = K smallest of (A.values ∪ B.values), filtered to < theta_out
183// (filtering keeps the unbiased estimator valid).
184
185func nx_theta_min2(a: i64, b: i64) -> i64 {
186 if a < b { return a }
187 return b
188}
189
190func nx_theta_union(a: *ThetaSketch, b: *ThetaSketch) -> *ThetaSketch {
191 if a.k != b.k { return 0 as *ThetaSketch }
192 if a.seed != b.seed { return 0 as *ThetaSketch }
193 let out: *ThetaSketch = nx_theta_alloc(a.k, a.seed)
194 out.theta = nx_theta_min2(a.theta, b.theta)
195 // Linear merge. Skip any hashes >= out.theta.
196 var i: i64 = 0
197 var j: i64 = 0
198 var w: i64 = 0
199 var done: i64 = 0
200 while done == 0 {
201 if w >= out.k { done = 1 }
202 if done == 0 {
203 let a_done: i64 = i >= a.n_items
204 let b_done: i64 = j >= b.n_items
205 if a_done == 1 {
206 if b_done == 1 { done = 1 }
207 if b_done == 0 {
208 let vb: i64 = b.values[j]
209 j = j + 1
210 if vb < out.theta {
211 out.values[w] = vb
212 w = w + 1
213 }
214 }
215 }
216 if done == 0 {
217 if a_done == 0 {
218 if b_done == 1 {
219 let va: i64 = a.values[i]
220 i = i + 1
221 if va < out.theta {
222 out.values[w] = va
223 w = w + 1
224 }
225 }
226 if b_done == 0 {
227 let va: i64 = a.values[i]
228 let vb: i64 = b.values[j]
229 if va == vb {
230 i = i + 1
231 j = j + 1
232 if va < out.theta {
233 out.values[w] = va
234 w = w + 1
235 }
236 }
237 if va < vb {
238 i = i + 1
239 if va < out.theta {
240 out.values[w] = va
241 w = w + 1
242 }
243 }
244 if va > vb {
245 j = j + 1
246 if vb < out.theta {
247 out.values[w] = vb
248 w = w + 1
249 }
250 }
251 }
252 }
253 }
254 }
255 }
256 out.n_items = w
257 return out
258}
259
260// === intersection ===============================================
261//
262// Beyer 2014 unbiased intersection:
263// theta_out = min(theta_a, theta_b)
264// matches = hashes present in BOTH A and B with value < theta_out
265// cardinality = |matches| * HASH_MAX / theta_out
266//
267// Storage: matched hashes sorted ascending.
268
269func nx_theta_intersect(a: *ThetaSketch, b: *ThetaSketch) -> *ThetaSketch {
270 if a.k != b.k { return 0 as *ThetaSketch }
271 if a.seed != b.seed { return 0 as *ThetaSketch }
272 let out: *ThetaSketch = nx_theta_alloc(a.k, a.seed)
273 out.theta = nx_theta_min2(a.theta, b.theta)
274 // Two-pointer scan over sorted A and B; keep equal hashes < theta_out.
275 var i: i64 = 0
276 var j: i64 = 0
277 var w: i64 = 0
278 while i < a.n_items {
279 if j >= b.n_items { i = a.n_items }
280 if i < a.n_items {
281 if j < b.n_items {
282 let va: i64 = a.values[i]
283 let vb: i64 = b.values[j]
284 if va == vb {
285 if va < out.theta {
286 if w < out.k {
287 out.values[w] = va
288 w = w + 1
289 }
290 }
291 i = i + 1
292 j = j + 1
293 }
294 if va < vb { i = i + 1 }
295 if va > vb { j = j + 1 }
296 }
297 }
298 }
299 out.n_items = w
300 return out
301}
302
303// === difference (A \ B) =========================================
304//
305// theta_out = theta_a
306// matches = hashes in A with value < theta_a that are NOT in B
307// (restricting to < theta_b; if h < theta_b and h ∉ B,
308// we can confidently say h ∉ B. If h >= theta_b,
309// we can't tell whether h ∈ B or not, so we exclude
310// it from the unbiased estimator).
311//
312// Cardinality of A \ B = |matches| * HASH_MAX / theta_out, valid
313// only when theta_a <= theta_b (else we'd be guessing about
314// hashes outside B's sampling region). Returns a sketch with
315// theta_out = min(theta_a, theta_b) to preserve the invariant.
316
317func nx_theta_difference(a: *ThetaSketch, b: *ThetaSketch) -> *ThetaSketch {
318 if a.k != b.k { return 0 as *ThetaSketch }
319 if a.seed != b.seed { return 0 as *ThetaSketch }
320 let out: *ThetaSketch = nx_theta_alloc(a.k, a.seed)
321 out.theta = nx_theta_min2(a.theta, b.theta)
322 var i: i64 = 0
323 var j: i64 = 0
324 var w: i64 = 0
325 var done: i64 = 0
326 while done == 0 {
327 if i >= a.n_items { done = 1 }
328 if done == 0 {
329 let va: i64 = a.values[i]
330 if va >= out.theta { done = 1 }
331 if done == 0 {
332 // Advance j while b.values[j] < va; halt at end of
333 // b or first b[j] >= va. Check bounds BEFORE reading.
334 var jdone: i64 = 0
335 while jdone == 0 {
336 if j >= b.n_items {
337 jdone = 1
338 }
339 if jdone == 0 {
340 if b.values[j] >= va {
341 jdone = 1
342 }
343 if jdone == 0 {
344 j = j + 1
345 }
346 }
347 }
348 // b.values[j] >= va, or j == b.n_items.
349 var in_b: i64 = 0
350 if j < b.n_items {
351 if b.values[j] == va { in_b = 1 }
352 }
353 if in_b == 0 {
354 if w < out.k {
355 out.values[w] = va
356 w = w + 1
357 }
358 }
359 i = i + 1
360 }
361 }
362 }
363 out.n_items = w
364 return out
365}
366
367// === introspection ==============================================
368
369func nx_theta_memory_bytes(t: *ThetaSketch) -> i64 {
370 return 48 + t.k * 8
371}