sketch_kmv.nx source
↩ module page · 310 lines · 10291 B
1// sketch_kmv.nx -- KMV (K-Minimum-Values) sketch.
2//
3// Bar-Yossef et al. 2002 + Beyer 2007. Maintain the K smallest
4// 32-bit hash values seen. Direct ancestor of DataSketches'
5// Theta sketches; we ship the foundation primitive first then
6// build Theta on top.
7//
8// CAPABILITIES (set-operation cardinality):
9// - cardinality: (K-1) * 2^32 / kth_smallest when |set| > K
10// total_inserted_unique when |set| <= K
11// - union: K smallest across A.values ∪ B.values
12// - Jaccard: |A.kmins ∩ B.kmins| / |A.kmins ∪ B.kmins|
13// unbiased estimator of true Jaccard
14//
15// Error: rel_stddev ~ 1/sqrt(K-1) at 1-sigma.
16// For K=4096: 1.56% stddev. For K=16384: 0.78%.
17//
18// COMPLEMENTS HLL:
19// - HLL: better space per accuracy point for pure cardinality.
20// - KMV: SAME sketch supports union AND intersection (Jaccard).
21// The trade-off is well-studied; both ship in DataSketches
22// (HLL.java and Theta.java) -- we ship the same axis pair.
23//
24// Hash domain: 32-bit (murmur3_32). We store hashes as i64 to
25// keep arithmetic simple; the value space is [0, 2^32). Hash 0
26// is mapped to 1 (matches Cuckoo's empty-sentinel pattern; KMV
27// uses the special "all-slots-filled-with-max" sentinel via
28// the n_items counter instead, but mapping 0->1 avoids edge
29// cases in the cardinality formula).
30
31import "syscalls.nx"
32import "murmur3.nx"
33import "sketch_types.nx"
34
35const NX_KMV_K_MIN: i64 = 16
36const NX_KMV_K_MAX: i64 = 65536
37
38// 2^32 in i64 -- the hash domain size used in cardinality estimation.
39const NX_KMV_HASH_MAX: i64 = 4294967296
40
41struct Kmv {
42 values: *i64, // sorted ascending; max at index n_items-1
43 k: i64,
44 n_items: i64, // = min(unique_inserts, k)
45 seed: i64,
46}
47
48// === construction =================================================
49
50func nx_kmv_alloc(k: i64, seed: i64) -> *Kmv {
51 if k < NX_KMV_K_MIN { return 0 as *Kmv }
52 if k > NX_KMV_K_MAX { return 0 as *Kmv }
53 let raw: *u8 = sys_mmap(40)
54 let kmv: *Kmv = raw as *Kmv
55 kmv.k = k
56 kmv.n_items = 0
57 let values_raw: *u8 = sys_mmap(k * 8)
58 kmv.values = values_raw as *i64
59 kmv.seed = seed
60 return kmv
61}
62
63// === hash helper =================================================
64
65func nx_kmv_hash(kmv: *Kmv, key: *u8, len: i64) -> i64 {
66 let h: i64 = murmur3_32(kmv.seed, key, len) & 0xFFFFFFFF
67 if h == 0 { return 1 } // remap 0 -> 1
68 return h
69}
70
71// === insertion =================================================
72//
73// Binary-search for position; if hash already present, no-op
74// (KMV is a SET of hashes). Else if set not full, insert and
75// shift tail right. Else compare with current max (values[k-1]);
76// if hash < max, replace max + sift down to maintain sort.
77
78func nx_kmv_add(kmv: *Kmv, key: *u8, len: i64) -> i64 {
79 let h: i64 = nx_kmv_hash(kmv, key, len)
80
81 // Binary search for h in values[0..n_items).
82 var lo: i64 = 0
83 var hi: i64 = kmv.n_items
84 while lo < hi {
85 let mid: i64 = (lo + hi) / 2
86 let v: i64 = kmv.values[mid]
87 if v == h { return 0 } // already in set
88 if v < h { lo = mid + 1 }
89 if v > h { hi = mid }
90 }
91 // `lo` is the insertion index.
92
93 if kmv.n_items < kmv.k {
94 // Shift values[lo..n_items) right by one, insert.
95 var j: i64 = kmv.n_items
96 while j > lo {
97 kmv.values[j] = kmv.values[j - 1]
98 j = j - 1
99 }
100 kmv.values[lo] = h
101 kmv.n_items = kmv.n_items + 1
102 return 0
103 }
104 // Full set: only insert if h < current max (values[k-1]).
105 let cur_max: i64 = kmv.values[kmv.k - 1]
106 if h >= cur_max { return 0 }
107 // Shift values[lo..k-1) right by one (overwriting old max), insert h.
108 var j2: i64 = kmv.k - 1
109 while j2 > lo {
110 kmv.values[j2] = kmv.values[j2 - 1]
111 j2 = j2 - 1
112 }
113 kmv.values[lo] = h
114 return 0
115}
116
117// Insert pre-hashed i64 value directly (for testing + union).
118func nx_kmv_add_hash(kmv: *Kmv, raw_h: i64) -> i64 {
119 var h: i64 = raw_h & 0xFFFFFFFF
120 if h == 0 { h = 1 }
121 var lo: i64 = 0
122 var hi: i64 = kmv.n_items
123 while lo < hi {
124 let mid: i64 = (lo + hi) / 2
125 let v: i64 = kmv.values[mid]
126 if v == h { return 0 }
127 if v < h { lo = mid + 1 }
128 if v > h { hi = mid }
129 }
130 if kmv.n_items < kmv.k {
131 var j: i64 = kmv.n_items
132 while j > lo {
133 kmv.values[j] = kmv.values[j - 1]
134 j = j - 1
135 }
136 kmv.values[lo] = h
137 kmv.n_items = kmv.n_items + 1
138 return 0
139 }
140 let cur_max: i64 = kmv.values[kmv.k - 1]
141 if h >= cur_max { return 0 }
142 var j2: i64 = kmv.k - 1
143 while j2 > lo {
144 kmv.values[j2] = kmv.values[j2 - 1]
145 j2 = j2 - 1
146 }
147 kmv.values[lo] = h
148 return 0
149}
150
151// === cardinality estimation ====================================
152//
153// Bar-Yossef 2002 estimator:
154// if n_items < k: cardinality = n_items (exact)
155// else: cardinality = (k - 1) * HASH_MAX / kth_min
156// kth_min is values[k-1] (the largest of the k smallest).
157//
158// For HASH_MAX = 2^32 = 4294967296, the multiplication
159// (k-1) * 2^32 fits in i64 well past k=2^30.
160
161func nx_kmv_estimate(kmv: *Kmv) -> i64 {
162 if kmv.n_items < kmv.k { return kmv.n_items }
163 let kth: i64 = kmv.values[kmv.k - 1]
164 if kth == 0 { return 0 }
165 return ((kmv.k - 1) * NX_KMV_HASH_MAX) / kth
166}
167
168// === typed query ===============================================
169//
170// rel_stddev = 1 / sqrt(k - 1). Tabulate by k class.
171
172func nx_kmv_stddev_rel_ppb(k: i64) -> i64 {
173 if k <= 16 { return 258000000 } // 1/sqrt(15) = 0.258
174 if k <= 64 { return 126000000 } // 1/sqrt(63) = 0.126
175 if k <= 256 { return 62700000 } // 1/sqrt(255)
176 if k <= 1024 { return 31300000 } // 1/sqrt(1023)
177 if k <= 4096 { return 15600000 } // 1/sqrt(4095) = 0.0156
178 if k <= 16384 { return 7810000 } // 1/sqrt(16383)
179 return 3900000 // 1/sqrt(65535)
180}
181
182func nx_kmv_query(kmv: *Kmv) -> *ApproxI64 {
183 let est: i64 = nx_kmv_estimate(kmv)
184 return nx_approx_new(est, NX_ENV_REL_STDDEV,
185 nx_kmv_stddev_rel_ppb(kmv.k),
186 682700000, // 0.6827 = 1-sigma
187 NX_MATURITY_REFERENCE_IMPL,
188 NX_ADV_HONEST)
189}
190
191// === union =====================================================
192//
193// Pre: a.k == b.k and a.seed == b.seed.
194// Result: new KMV holding the K smallest across A ∪ B.
195//
196// Linear merge of two sorted streams, K-bounded.
197
198func nx_kmv_union(a: *Kmv, b: *Kmv) -> *Kmv {
199 if a.k != b.k { return 0 as *Kmv }
200 if a.seed != b.seed { return 0 as *Kmv }
201 let out: *Kmv = nx_kmv_alloc(a.k, a.seed)
202 var i: i64 = 0
203 var j: i64 = 0
204 var w: i64 = 0
205 while w < out.k {
206 let a_done: i64 = i >= a.n_items
207 let b_done: i64 = j >= b.n_items
208 if a_done == 1 {
209 if b_done == 1 { w = out.k; return out }
210 out.values[w] = b.values[j]
211 j = j + 1
212 w = w + 1
213 }
214 if a_done == 0 {
215 if b_done == 1 {
216 out.values[w] = a.values[i]
217 i = i + 1
218 w = w + 1
219 }
220 if b_done == 0 {
221 let va: i64 = a.values[i]
222 let vb: i64 = b.values[j]
223 if va < vb {
224 out.values[w] = va
225 i = i + 1
226 w = w + 1
227 }
228 if va == vb {
229 out.values[w] = va
230 i = i + 1
231 j = j + 1
232 w = w + 1
233 }
234 if va > vb {
235 out.values[w] = vb
236 j = j + 1
237 w = w + 1
238 }
239 }
240 }
241 // Early-exit when both sides exhausted.
242 if i >= a.n_items {
243 if j >= b.n_items { w = out.k }
244 }
245 }
246 // Set n_items based on whether we filled the cap or ran out.
247 // If both sources exhausted before w == k, n_items == w.
248 // Otherwise we wrote exactly k values.
249 out.n_items = w
250 return out
251}
252
253// === Jaccard similarity ========================================
254//
255// J(A, B) = |A ∩ B| / |A ∪ B|
256// ≈ |kmins(A) ∩ kmins(B)| / |kmins(A) ∪ kmins(B)|
257// where intersection/union are computed over the K-min sets.
258//
259// Two-pointer linear scan over sorted hash lists. Returns
260// Jaccard in parts-per-million (0 = disjoint, 1_000_000 = identical).
261
262func nx_kmv_jaccard_ppm(a: *Kmv, b: *Kmv) -> i64 {
263 if a.k != b.k { return 0 }
264 if a.seed != b.seed { return 0 }
265 if a.n_items == 0 {
266 if b.n_items == 0 { return 1000000 } // both empty: defined as identical
267 return 0
268 }
269 if b.n_items == 0 { return 0 }
270 var i: i64 = 0
271 var j: i64 = 0
272 var intersect: i64 = 0
273 var unite: i64 = 0
274 while i < a.n_items {
275 if j >= b.n_items { i = a.n_items }
276 if i < a.n_items {
277 if j < b.n_items {
278 let va: i64 = a.values[i]
279 let vb: i64 = b.values[j]
280 if va == vb {
281 intersect = intersect + 1
282 unite = unite + 1
283 i = i + 1
284 j = j + 1
285 }
286 if va < vb {
287 unite = unite + 1
288 i = i + 1
289 }
290 if va > vb {
291 unite = unite + 1
292 j = j + 1
293 }
294 }
295 }
296 }
297 // Drain remaining b.
298 while j < b.n_items {
299 unite = unite + 1
300 j = j + 1
301 }
302 if unite == 0 { return 0 }
303 return (intersect * 1000000) / unite
304}
305
306// === introspection ==============================================
307
308func nx_kmv_memory_bytes(kmv: *Kmv) -> i64 {
309 return 40 + kmv.k * 8
310}