sketch_count_sketch.nx source
↩ module page · 216 lines · 7758 B
1// sketch_count_sketch.nx -- CountSketch (Charikar-Chen-Farach-Colton 2002).
2//
3// Frequency-counting variant that gives UNBIASED estimates via ±1
4// sign hashing. Each (key, count) update adds s_j(key) * count to
5// row j column h_j(key), where s_j ∈ {-1, +1} is a sign hash.
6// Query: median of (s_j(key) * counter[j][h_j(key)]) across rows.
7//
8// COMPLEMENTS CMS (runtime/sketch_cms.nx):
9// - CMS: always OVERESTIMATES (false positives via positive collisions).
10// Use when "upper bound" matters.
11// - CountSketch: UNBIASED. Use when expectation matters (statistical
12// summaries, F_2 estimation, sparse-approximation literature).
13//
14// CAPABILITY STOMP: shipping both means callers pick by error semantics.
15// DataSketches ships CMS as "FrequentLongs"; CountSketch is queued in
16// their docs but not implemented.
17//
18// ERROR BOUND:
19// |estimate(x) - f(x)| <= ε · ||f||_2 / sqrt(d)
20// where ||f||_2 = sqrt(Σ f_i²) is the L2 norm of the frequency vector.
21// For w = O(1/ε²) and d = O(log(1/δ)) rows: confidence 1-δ.
22//
23// Default params: d=5 rows, w=2048 columns -> conf >99%, ε ~ 0.022.
24
25import "syscalls.nx"
26import "sketch_types.nx"
27import "nx_vecmath.nx"
28
29const NX_CS_MIN_D: i64 = 3
30const NX_CS_MAX_D: i64 = 32
31const NX_CS_MIN_W: i64 = 64
32const NX_CS_MAX_W: i64 = 65536
33
34struct CountSketch {
35 counters: *i64, // d * w grid, in i64 (signed)
36 d: i64,
37 w: i64,
38 seed: i64,
39 total: i64,
40 scratch: *i64, // d-sized scratch for median (bits-up: hoisted
41 // from per-call sys_mmap to alloc-time). Saves
42 // one syscall per nx_cs_estimate call.
43}
44
45// === construction =================================================
46
47func nx_cs_alloc(d: i64, w: i64, seed: i64) -> *CountSketch {
48 if d < NX_CS_MIN_D { return 0 as *CountSketch }
49 if d > NX_CS_MAX_D { return 0 as *CountSketch }
50 if w < NX_CS_MIN_W { return 0 as *CountSketch }
51 if w > NX_CS_MAX_W { return 0 as *CountSketch }
52 if (w & (w - 1)) != 0 { return 0 as *CountSketch } // power of 2
53 let raw: *u8 = sys_mmap(56)
54 let c: *CountSketch = raw as *CountSketch
55 let cells: i64 = d * w
56 let cells_raw: *u8 = sys_mmap(cells * 8)
57 c.counters = cells_raw as *i64
58 var i: i64 = 0
59 while i < cells {
60 c.counters[i] = 0
61 i = i + 1
62 }
63 c.d = d
64 c.w = w
65 c.seed = seed
66 c.total = 0
67 let scratch_raw: *u8 = sys_mmap(d * 8)
68 c.scratch = scratch_raw as *i64
69 return c
70}
71
72// === hash + sign helpers =========================================
73//
74// Kirsch-Mitzenmacher double-hashing: h_j(x) = (h1 + j*h2) & mask.
75// Sign: top bit of a separately-derived hash, mapped to ±1.
76
77func nx_cs_h1(c: *CountSketch, key: i64) -> i64 {
78 let mixed: i64 = (key * 0x9E3779B97F4A7C15 + c.seed) & 0xFFFFFFFFFFFFFFFF
79 return mixed & 0xFFFFFFFF
80}
81
82func nx_cs_h2(c: *CountSketch, key: i64) -> i64 {
83 let mixed: i64 = (key * 0xBF58476D1CE4E5B9 + c.seed) & 0xFFFFFFFFFFFFFFFF
84 return mixed & 0xFFFFFFFF
85}
86
87func nx_cs_column(c: *CountSketch, key: i64, j: i64) -> i64 {
88 let combined: i64 = (nx_cs_h1(c, key) + j * nx_cs_h2(c, key)) & 0xFFFFFFFF
89 return combined & (c.w - 1)
90}
91
92// Sign for row j: hash key+j and check the top bit.
93func nx_cs_sign(c: *CountSketch, key: i64, j: i64) -> i64 {
94 let mixed: i64 = ((key + j * 0x9E3779B9) * 0xC2B2AE3D27D4EB4F + c.seed) & 0xFFFFFFFFFFFFFFFF
95 if (mixed & (1 << 63)) == 0 { return 1 }
96 return -1
97}
98
99// === cell access =================================================
100
101func nx_cs_cell_idx(c: *CountSketch, j: i64, col: i64) -> i64 {
102 return j * c.w + col
103}
104
105// === add ==========================================================
106
107func nx_cs_add(c: *CountSketch, key: i64, count: i64) -> i64 {
108 if count == 0 { return 0 }
109 c.total = c.total + count
110 // Bits-up: compute h1 and h2 ONCE per add; were recomputed per row.
111 // Saves (d-1)*2 multiplies per add.
112 let h1: i64 = nx_cs_h1(c, key)
113 let h2: i64 = nx_cs_h2(c, key)
114 let w_mask: i64 = c.w - 1
115 var j: i64 = 0
116 while j < c.d {
117 let col: i64 = ((h1 + j * h2) & 0xFFFFFFFF) & w_mask
118 let sign: i64 = nx_cs_sign(c, key, j)
119 let idx: i64 = j * c.w + col
120 c.counters[idx] = c.counters[idx] + sign * count
121 j = j + 1
122 }
123 return 0
124}
125
126// === query (median of signed reads) ==============================
127//
128// For each row j: signed_read = sign_j(key) * counter[j][h_j(key)]
129// Return median of d signed reads. Median is approximate median via
130// insertion-sort of small d.
131
132func nx_cs_estimate(c: *CountSketch, key: i64) -> i64 {
133 // Bits-up: scratch hoisted to struct; h1+h2 cached.
134 let scratch: *i64 = c.scratch
135 let h1: i64 = nx_cs_h1(c, key)
136 let h2: i64 = nx_cs_h2(c, key)
137 let w_mask: i64 = c.w - 1
138 var j: i64 = 0
139 while j < c.d {
140 let col: i64 = ((h1 + j * h2) & 0xFFFFFFFF) & w_mask
141 let sign: i64 = nx_cs_sign(c, key, j)
142 let idx: i64 = j * c.w + col
143 scratch[j] = sign * c.counters[idx]
144 j = j + 1
145 }
146 // Insertion sort.
147 var i: i64 = 1
148 while i < c.d {
149 let cur: i64 = scratch[i]
150 var k: i64 = i - 1
151 var done: i64 = 0
152 while done == 0 {
153 if k < 0 { done = 1 }
154 if done == 0 {
155 if scratch[k] <= cur { done = 1 }
156 if done == 0 {
157 scratch[k + 1] = scratch[k]
158 k = k - 1
159 }
160 }
161 }
162 scratch[k + 1] = cur
163 i = i + 1
164 }
165 return scratch[c.d / 2]
166}
167
168// === typed envelope ===============================================
169//
170// Error bound is |est - f| <= ||f||_2 * sqrt(2/w) at confidence
171// (3/4)^(d/2). We declare a CONSERVATIVE absolute bound = total/sqrt(w)
172// (this is the worst case when L2 = total, i.e. one massive item).
173// param_a holds the bound; conf = (3/4)^(d/2) tabulated.
174
175func nx_cs_isqrt(x: i64) -> i64 { return vm_isqrt(x) }
176
177// Confidence depends on d: (3/4)^(d/2). Approximate by tabulating.
178func nx_cs_conf_ppb(d: i64) -> i64 {
179 if d <= 3 { return 562500000 } // (3/4)^1.5 = 0.65
180 if d <= 5 { return 421900000 } // (3/4)^2.5
181 if d <= 7 { return 316400000 } // (3/4)^3.5
182 if d <= 11 { return 177900000 } // (3/4)^5.5
183 return 100000000 // higher d converges
184}
185
186func nx_cs_query(c: *CountSketch, key: i64) -> *ApproxI64 {
187 let est: i64 = nx_cs_estimate(c, key)
188 let bound: i64 = c.total / nx_cs_isqrt(c.w)
189 // Returned envelope: 1e9 - conf_ppb confidence the bound holds;
190 // i.e. there's a small chance the estimate is further than `bound`.
191 return nx_approx_new(est, NX_ENV_ABS, bound,
192 1000000000 - nx_cs_conf_ppb(c.d),
193 NX_MATURITY_REFERENCE_IMPL,
194 NX_ADV_HONEST)
195}
196
197// === merge ========================================================
198
199func nx_cs_merge(a: *CountSketch, b: *CountSketch) -> *CountSketch {
200 if a.d != b.d { return 0 as *CountSketch }
201 if a.w != b.w { return 0 as *CountSketch }
202 if a.seed != b.seed { return 0 as *CountSketch }
203 let out: *CountSketch = nx_cs_alloc(a.d, a.w, a.seed)
204 let cells: i64 = a.d * a.w
205 var i: i64 = 0
206 while i < cells {
207 out.counters[i] = a.counters[i] + b.counters[i]
208 i = i + 1
209 }
210 out.total = a.total + b.total
211 return out
212}
213
214func nx_cs_memory_bytes(c: *CountSketch) -> i64 {
215 return 48 + c.d * c.w * 8
216}