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