From 51eb79f8e30c18298723cee77a745b308030f1db Mon Sep 17 00:00:00 2001 From: zz_y Date: Wed, 14 Oct 2026 21:31:56 +0000 Subject: [PATCH] fix(univmon_optimized): unbiased distinct count and direct L2 for UnivMonPyramid UnivMonPyramid had the same two estimator problems as UnivMon: - calc_card dropped heap keys whose estimate was not positive, or below L2/sqrt(k) on an incomplete layer. - calc_l2 ran the heap recursion with g = x^2 instead of reading the Count Sketch. calc_card now counts every heap key. calc_l2 reads layer 0's Count Sketch in standard mode and sums the strata's squared estimates in terminal mode. Adds top_layer_complete(). The doc comments state the bounds in the pyramid's own terms: - L2 uses layer 0's dimensions in standard mode, and each stratum's own in terminal mode. - The distinct-count bound is independent of the elephant/mouse split. Before the fix, a 5x64 layer 0 missed by more than eps = 0.5 in 39 of 150 Zipf trials, against delta = 0.016. The direct estimate stays within the bound. The Monte Carlo checks have ignored long variants. The e2e univmon_pyramid_weighted_metrics test checked the distinct count within 20% at heap 32, where the unbiased estimator's standard deviation can reach 25% of F0. It now uses the default tiers with heap 128 and asserts top_layer_complete(), as #221 did for UnivMon's documented matrix. --- docs/api/api_univmon_optimized.md | 10 + docs/tests.md | 8 +- src/sketch_framework/univmon_optimized.rs | 305 ++++++++++++++++++++-- tests/e2e/frameworks.rs | 10 +- 4 files changed, 309 insertions(+), 24 deletions(-) diff --git a/docs/api/api_univmon_optimized.md b/docs/api/api_univmon_optimized.md index 2a28521b..5ff7c39f 100644 --- a/docs/api/api_univmon_optimized.md +++ b/docs/api/api_univmon_optimized.md @@ -46,6 +46,7 @@ fn calc_l2(&self) -> f64 fn calc_entropy(&self) -> f64 fn calc_card(&self) -> f64 fn calc_g_sum(&self, g: F, is_card: bool) -> f64 +fn top_layer_complete(&self) -> bool // UnivSketchPool fn available(&self) -> usize @@ -56,6 +57,15 @@ For the supported non-negative update stream, `calc_l1()` returns the exact tracked total weight rather than estimating the linear sum through the hierarchy. +The pyramid's estimators are `UnivMon`'s (see [UnivMon](api_univmon.md)). +`calc_l2()` reads layer 0's Count Sketch directly in standard mode, so its +`(ε, δ)` bound uses layer 0's dimensions (the elephant ones whenever +`elephant_layers ≥ 1`); in terminal mode it sums the strata, each held to its +own layer's dimensions. `calc_card()` counts every heap key and, when +`top_layer_complete()` holds, misses by more than `ε F0` with probability at +most `2/(kε²)` for `k = heap_size ≥ 8`, independent of the elephant/mouse +split. + ## Merge ```rust diff --git a/docs/tests.md b/docs/tests.md index 96bb16eb..537314ce 100644 --- a/docs/tests.md +++ b/docs/tests.md @@ -875,6 +875,12 @@ Wire tests: [`src/sketch_framework/univmon_optimized/wire.rs`](../src/sketch_fra | `pyramid_accuracy_zipf` | Pyramid accuracy Zipf. | On heavy/medium/light Zipf-like workload, requires relative error `<15%` for `L1`, `L2`, cardinality, and entropy. | | `pyramid_fast_insert_accuracy` | Pyramid fast insert accuracy. | Using `fast_insert` only, requires relative error `<15%` for `L1`, `L2`, cardinality, and entropy versus exact frequency map. | | `pyramid_memory_savings_vs_uniform` | Pyramid memory savings vs uniform. | Verifies pyramid column budget is smaller than uniform baseline and computed memory savings exceed `30%`. | +| `pyramid_calc_card_counts_heap_keys_with_nonpositive_estimates` | A heap key counts toward the distinct count whatever its estimate. | In a one-layer `1x1` pyramid, inserts a `+1`-signed key with weight `5` and a `-1`-signed key with weight `1`, verifies the light key's estimate is not positive and the top layer is complete, and that `calc_card()` is `2.0`. | +| `pyramid_calc_l2_reads_layer_zero_directly_in_standard_mode` | Standard-mode L2 is layer 0's own Count Sketch estimate. | After 2000 weighted keys into `UnivMonPyramid::new(16, 2, 5, 256, 3, 64, 8)` through `insert`, verifies `calc_l2()` equals layer 0's `get_l2()` within `1e-9`. | +| `pyramid_calc_l2_sums_stratum_estimates_in_terminal_mode` | Terminal-mode L2 sums the disjoint strata. | After the same keys through `fast_insert`, verifies `calc_l2()` equals `sqrt` of the summed per-layer `get_l2_sqr()` within `1e-9`. | +| `pyramid_top_layer_complete_tracks_the_top_heap` | The distinct-count precondition is readable. | On a one-layer heap-4 pyramid, verifies `top_layer_complete()` holds after 4 keys, fails after a fifth, and holds again after `free()`. | +| `pyramid_calc_l2_meets_its_documented_bound` | `calc_l2` meets the `(ε, δ)` bound on its doc comment. | 150 seeded trials of 400 fresh random Zipf-weighted keys into `UnivMonPyramid::new(8, 1, 5, 64, 3, 16, 3)`: trials with `calc_l2()^2` off by more than `ε = 0.5` relative stay within `n·δ + 4σ + 1` for `δ = P(Bin(5, 0.125) >= 3)`. `_long` (ignored) runs 3000 trials. | +| `pyramid_calc_card_meets_its_documented_bound` | `calc_card` meets its bound when the top layer is complete. | 60 seeded trials of 2000 fresh random keys into a heap-32 pyramid with `3x256` elephant and `2x32` mouse layers, `8 = ceil(log2(2000/32)) + 2` deep: incomplete top layers stay rare, misses beyond `ε = 0.5` stay within the Chebyshev `2/(kε²)` bound, the variance stays under `2 F0²/k`, and the mean error is within four standard errors of zero. `_long` (ignored) runs 2000 trials. | | `pyramid_round_trip_serialization` | The envelope frames a pyramid under kind_id `0x11 0x00`, and the bytes are stable across a round trip. | For `UnivMonPyramid::new(4, 2, 3, 16, 2, 8, 4)` fed four weighted string keys, verifies the bytes open with the ASAPv1 magic, a `kind_id_len` of `2`, and `0x11 0x00`, that the metadata carries that layout and `key_type` `string`, and that the decode matches on `bucket_size`, L1, L2, cardinality, entropy, and the candidate flags, re-serializing identically. | | `pyramid_two_tier_geometry_survives` | The two tiers are derived from the layer's position. | Verifies every decoded layer keeps its tier's dimensions - `3x16` for the two elephant layers and `2x8` for the mouse layers - its own seed index, and its counter array. | | `pyramid_layers_with_different_heap_loads_round_trip` | Each layer's heap contents survive the round trip. | For `UnivMonPyramid::new(8, 2, 3, 32, 2, 16, 4)` fed 40 weighted keys, verifies the layers carry different heap loads, that every layer's length and per-key counts match after a decode, and that the bytes re-serialize identically. | @@ -1858,7 +1864,7 @@ Test file: [`tests/e2e/frameworks.rs`](../tests/e2e/frameworks.rs) | test_name | test_description | what_is_tested | | --- | --- | --- | | `univmon_weighted_metrics_and_fast_insert_parity` | UnivMon's weighted L1 is exact and its L2, entropy and cardinality hold their bands, with `fast_insert` tracking the standard path. | Drives two `UnivMon::init_univmon(32, 5, 2048, 8)` sketches over `zipf_u64(20_000, 1000, 1.2, 4003)` with weight `1 + i % 7`, one through `insert` and one through `fast_insert`, then asserts both `calc_l1()` equal the exact weighted total with `assert_eq!`; `calc_l2()` inside `univmon_l2_band(exact, 5, 2048)` (`b = sqrt(2*3/2048)`, about `0.054`); `calc_entropy()` in `[0.88, 1.12]` times the exact bits; `calc_card()` in `[0.94, 1.06]` times the exact distinct count; and the fast path's `calc_l2()` in `[0.85, 1.15]` times exact. | -| `univmon_pyramid_weighted_metrics` | UnivMonPyramid's weighted L1 is exact and its L2 and cardinality hold their bands. | Drives `UnivMonPyramid::with_defaults()` over `zipf_u64(15_000, 800, 1.3, 4004)` with weight `1 + i % 3` and asserts `calc_l1()` equals the exact weighted total with `assert_eq!`, `calc_l2()` is in `[0.85, 1.15]` times exact, and `calc_card()` in `[0.80, 1.20]` times the exact distinct count; entropy is not checked. | +| `univmon_pyramid_weighted_metrics` | UnivMonPyramid's weighted L1 is exact and its L2 and cardinality hold their bands. | Drives `UnivMonPyramid::new(128, 8, 3, 2048, 3, 512, 16)` (the default tiers with heap 128, so the 20% band is over 1.6 of `calc_card`'s documented standard deviations) over `zipf_u64(15_000, 800, 1.3, 4004)` with weight `1 + i % 3` and asserts `calc_l1()` equals the exact weighted total with `assert_eq!`, `calc_l2()` is in `[0.85, 1.15]` times exact, that `top_layer_complete()` holds, and `calc_card()` in `[0.80, 1.20]` times the exact distinct count; entropy is not checked. | | `univmon_generic_g_sum_reproduces_every_named_estimator` | The generic g-sum reproduces the cardinality and entropy estimators and its own heuristic delegate, and tracks L2. | Over `UnivMon::init_univmon(32, 5, 2_048, 8)` fed `zipf_u64(20_000, 1_000, 1.2, 4_201)` with weight `1 + i % 7`, `assert_eq!`s `calc_l2()` against layer 0's `get_l2()` (the L2 estimator reads layer 0 directly), requires `calc_g_sum(\|x\| x*x, false).sqrt()` within `[0.85, 1.15]` of the exact L2, `assert_eq!`s `calc_g_sum(\|_\| 1.0, true)` against `calc_card()`, `calc_g_sum(x^2, false)` against `calc_g_sum_heuristic(x^2, false)`, and `calc_g_sum(\|_\| 0.0, false)` against `0.0`; requires `calc_entropy()` to equal `log2(L1) - gsum(x*log2 x)/L1` within `1e-9`; and checks `calc_g_sum(\|x\| x, false)` only against `[0.50, 1.50]` times the exact L1, so the identity g-sum's agreement with L1 is not pinned tightly. | | `univmon_g_sum_is_zero_on_an_empty_sketch` | An empty UnivMon reports zero for every estimator. | On a fresh `UnivMon::init_univmon(32, 5, 256, 4)`, `assert_eq!`s `calc_l1()`, `calc_entropy()`, `calc_g_sum(\|x\| x*x, false)` and `calc_card()` each against `0.0`. | | `a_univmon_pool_recycles_a_returned_sketch_and_hands_back_a_cleared_one` | A returned sketch becomes available again and comes back cleared without a new allocation. | On `UnivSketchPool::new(2, 16, 3, 256, 4)`: asserts `available() == 2` and `total_allocated() == 2` when fresh, `available() == 0` after two `take()`s with `total_allocated()` still `2`, and `total_allocated() == 3` after a third `take()` past capacity; then fills the first sketch with `zipf_u64(5_000, 512, 1.1, 4_202)` inserts (checking `calc_l1() > 0.0`), `put()`s it back and asserts `available() == 1`, and requires the next `take()` to report `calc_l1() == 0.0` with `total_allocated()` still `3`. | diff --git a/src/sketch_framework/univmon_optimized.rs b/src/sketch_framework/univmon_optimized.rs index aea4cd57..d07d1d92 100644 --- a/src/sketch_framework/univmon_optimized.rs +++ b/src/sketch_framework/univmon_optimized.rs @@ -270,6 +270,11 @@ impl UnivMonPyramid { // -- Query methods (identical to UnivMon) -------------------------------- /// Computes a g-sum estimate. + /// + /// With `is_card` every heap key contributes `g` of its estimate; that + /// is [`Self::calc_card`]'s estimator. Otherwise a key whose estimate is + /// not positive is skipped, and in terminal mode a layer whose candidate + /// set is incomplete also skips keys below `L2 / √heap_size`. pub fn calc_g_sum(&self, g: F, is_card: bool) -> f64 where F: Fn(f64) -> f64, @@ -278,23 +283,19 @@ impl UnivMonPyramid { return 0.0; } if self.update_mode == UnivMonUpdateMode::Terminal { - return self.calc_terminal_g_sum(g); + return self.calc_terminal_g_sum(g, is_card); } + // A distinct count weighs every heap key 1 whatever its estimate: the + // recursion stays unbiased only if no key is dropped by its value. + let counted = |count: i64| is_card || count > 0; let mut y = vec![0.0; self.layer_size]; - let l2_top = self.l2_sketch_layers[self.layer_size - 1].get_l2(); - let threshold_top = if is_card { - self.heavy_threshold(l2_top, self.candidate_complete[self.layer_size - 1]) - } else { - 0 - }; - let mut tmp = 0.0; for item in self.hh_layers[self.layer_size - 1].heap() { let input = heap_item_to_sketch_input(&item.key); let count = self.l2_sketch_layers[self.layer_size - 1].estimate(&input) as i64; - if count > threshold_top { + if counted(count) { tmp += g(count as f64); } } @@ -302,17 +303,10 @@ impl UnivMonPyramid { for i in (0..(self.layer_size - 1)).rev() { tmp = 0.0; - let l2_val = self.l2_sketch_layers[i].get_l2(); - let threshold = if is_card { - self.heavy_threshold(l2_val, self.candidate_complete[i]) - } else { - 0 - }; - for item in self.hh_layers[i].heap() { let input = heap_item_to_sketch_input(&item.key); let count = self.l2_sketch_layers[i].estimate(&input) as i64; - if count > threshold { + if counted(count) { let hash = (hash_item64_seeded(BOTTOM_LAYER_FINDER, &item.key) >> (i + 1)) & 1; let coe = 1.0 - 2.0 * (hash as f64); tmp += coe * g(count as f64); @@ -323,11 +317,13 @@ impl UnivMonPyramid { y[0] } - fn calc_terminal_g_sum(&self, g: F) -> f64 + fn calc_terminal_g_sum(&self, g: F, is_card: bool) -> f64 where F: Fn(f64) -> f64, { let (candidates, complete) = self.logical_terminal_candidates(); + // As in standard mode, a distinct count weighs every candidate 1. + let counted = |count: i64, threshold: i64| is_card || count > threshold; let mut logical_l2 = vec![0.0; self.layer_size]; let mut suffix_l2_squared = 0.0; for level in (0..self.layer_size).rev() { @@ -341,7 +337,7 @@ impl UnivMonPyramid { let threshold = self.heavy_threshold(logical_l2[last], complete[last]); y[last] = candidates[last] .iter() - .filter(|(_, count)| *count > threshold) + .filter(|(_, count)| counted(*count, threshold)) .map(|(_, count)| g(*count as f64)) .sum(); @@ -349,7 +345,7 @@ impl UnivMonPyramid { let threshold = self.heavy_threshold(logical_l2[level], complete[level]); let correction = candidates[level] .iter() - .filter(|(_, count)| *count > threshold) + .filter(|(_, count)| counted(*count, threshold)) .map(|(key, count)| { let hash = (hash_item64_seeded(BOTTOM_LAYER_FINDER, key) >> (level + 1)) & 1; (1.0 - 2.0 * hash as f64) * g(*count as f64) @@ -408,8 +404,45 @@ impl UnivMonPyramid { } /// Returns the estimated L2 norm. + /// + /// In standard mode layer 0 receives every key, so its Count Sketch + /// estimates `F2 = Σ f_i²` directly (each row's sum of squared counters, + /// median across rows). In terminal mode the layers hold disjoint strata, + /// so their `F2` estimates add. The heap recursion of + /// [`Self::calc_g_sum`] is not used. + /// + /// # Error bound (standard mode) + /// + /// The derivation is [`UnivMon::calc_l2`]'s, applied to layer 0, which is + /// an elephant layer (`elephant_row x elephant_col`) whenever + /// `elephant_layers ≥ 1` and a mouse layer otherwise. With `w` and `r` + /// that layer's columns and rows: one row's estimate has + /// `Var ≤ 2q F2²`, where `q = 1/w` for a power-of-two `w` and + /// `q ≤ 9/(8w)` otherwise, so it misses by more than `ε F2` with + /// probability at most `p = 2q/ε²`. The rows read disjoint slices of one + /// 128-bit hash and are independent, and their median misses only when + /// `⌈r/2⌉` rows miss on the same side: + /// + /// `P(|F2_est − F2| > ε F2) ≤ δ = P(Binomial(r, p) ≥ ⌈r/2⌉)`. + /// + /// The same `ε` and `δ` bound the returned norm's relative error. The + /// mouse layers do not enter the standard-mode estimate. + /// + /// # Terminal mode + /// + /// Each stratum satisfies the bound above with its own layer's `w` and + /// `r` (mouse layers are narrower, so their `δ` is larger). The sum is + /// within `ε` of the total whenever every non-empty stratum is within `ε` + /// of its own `F2`, so the failure probability is at most the sum of the + /// non-empty strata's `δ`s. pub fn calc_l2(&self) -> f64 { - self.calc_g_sum(|x| x * x, false).sqrt() + if self.update_mode == UnivMonUpdateMode::Terminal { + return (0..self.layer_size) + .map(|level| self.l2_sketch_layers[level].get_l2_sqr()) + .sum::() + .sqrt(); + } + self.l2_sketch_layers[0].get_l2() } /// Returns the estimated entropy. @@ -421,11 +454,41 @@ impl UnivMonPyramid { (self.bucket_size as f64).log2() - tmp / (self.bucket_size as f64) } - /// Returns the estimated cardinality. + /// Returns the estimated number of distinct keys, `F0`. + /// + /// This is the recursion of [`Self::calc_g_sum`] with `g = 1`, and every + /// heap key counts whatever its Count Sketch estimate: a key is in a heap + /// only because the stream inserted it. + /// + /// # Error bound (standard mode) + /// + /// The derivation is [`UnivMon::calc_card`]'s and does not depend on the + /// sketch dimensions, so the elephant/mouse split does not enter it. With + /// `k = heap_size` and `S_j`, `H_j`, `M_j = S_j \ H_j` the keys reaching, + /// held at and missed at layer `j`, a complete top layer + /// ([`Self::top_layer_complete`]) gives + /// `Y_0 − F0 = −Σ_{j=0}^{L−2} 2^j Σ_{i∈M_j} (1 − 2 b_{j+1}(i))`. + /// `M_j` is decided by layer `j`'s own sketch over `S_j` and does not + /// depend on bit `j+1`, so the estimate is unbiased with + /// `Var[Y_0] = Σ_j 4^j E(|S_j| − k)^+ ≤ 1.35 F0²/k ≤ 2 F0²/k` for `k ≥ 8`. + /// By Chebyshev, `P(|Y_0 − F0| > ε F0) ≤ 2/(k ε²)`. + /// + /// With `layer_size ≥ ⌈log2(F0_max / k)⌉ + 2` the top layer overflows + /// with probability below `(e/4)^{k/2}`; read + /// [`Self::top_layer_complete`] rather than assume it. Terminal mode + /// reconstructs its candidate sets from the strata and is not covered. pub fn calc_card(&self) -> f64 { self.calc_g_sum(|_| 1.0, true) } + /// Whether the top layer's heap still holds every key that layer received. + /// + /// This is the precondition of [`Self::calc_card`]'s error bound: with an + /// incomplete top layer the distinct count is biased low. + pub fn top_layer_complete(&self) -> bool { + self.candidate_complete[self.layer_size - 1] + } + // -- Lifecycle ----------------------------------------------------------- /// Resets all counters and heaps without deallocating. @@ -507,6 +570,204 @@ impl UnivMonPyramid { mod tests { use super::*; use crate::DataInput; + use rand::{Rng, SeedableRng, rngs::StdRng}; + + /// Sign layer 0 gives `key` in its single row (`+1` or `-1`). + fn layer0_sign(key: &DataInput, cols: usize) -> i64 { + use crate::common::hash128_seeded; + use crate::common::structures::matrix_storage::cols_mask_bits; + use crate::sketches::countsketch_topk::l2hh_cell_for_row; + l2hh_cell_for_row(hash128_seeded(0, key), 0, cols, cols_mask_bits(cols)).1 + } + + /// A key whose Count Sketch estimate is not positive is still a key the + /// stream inserted; a complete layer must count it. + #[test] + fn pyramid_calc_card_counts_heap_keys_with_nonpositive_estimates() { + let plus = (0_u64..) + .map(DataInput::U64) + .find(|key| layer0_sign(key, 1) == 1) + .expect("a key with sign +1"); + let minus = (0_u64..) + .map(DataInput::U64) + .find(|key| layer0_sign(key, 1) == -1) + .expect("a key with sign -1"); + let mut pyramid = UnivMonPyramid::new(4, 1, 1, 1, 1, 1, 1); + pyramid.insert(&plus, 5); + pyramid.insert(&minus, 1); + assert!(pyramid.l2_sketch_layers[0].estimate(&minus) <= 0.0); + assert!(pyramid.top_layer_complete()); + assert_eq!(pyramid.calc_card(), 2.0); + } + + /// Layer 0 sees every key in standard mode, so its own Count Sketch is the + /// L2 estimate. + #[test] + fn pyramid_calc_l2_reads_layer_zero_directly_in_standard_mode() { + let mut pyramid = UnivMonPyramid::new(16, 2, 5, 256, 3, 64, 8); + for key in 0..2000_u64 { + pyramid.insert(&DataInput::U64(key), (key % 7 + 1) as i64); + } + let direct = pyramid.l2_sketch_layers[0].get_l2(); + assert!( + (pyramid.calc_l2() - direct).abs() < 1e-9, + "calc_l2 {} differs from layer 0's direct estimate {direct}", + pyramid.calc_l2() + ); + } + + /// Terminal mode keeps disjoint strata, so the squared L2 estimates add. + #[test] + fn pyramid_calc_l2_sums_stratum_estimates_in_terminal_mode() { + let mut pyramid = UnivMonPyramid::new(16, 2, 5, 256, 3, 64, 8); + for key in 0..2000_u64 { + pyramid.fast_insert(&DataInput::U64(key), (key % 7 + 1) as i64); + } + let direct = (0..pyramid.layer_size) + .map(|level| pyramid.l2_sketch_layers[level].get_l2_sqr()) + .sum::() + .sqrt(); + assert!( + (pyramid.calc_l2() - direct).abs() < 1e-9, + "calc_l2 {} differs from the summed strata {direct}", + pyramid.calc_l2() + ); + } + + #[test] + fn pyramid_top_layer_complete_tracks_the_top_heap() { + let mut pyramid = UnivMonPyramid::new(4, 1, 3, 64, 2, 16, 1); + for key in 0..4_u64 { + pyramid.insert(&DataInput::U64(key), 1); + } + assert!(pyramid.top_layer_complete()); + pyramid.insert(&DataInput::U64(4), 1); + assert!(!pyramid.top_layer_complete()); + pyramid.free(); + assert!(pyramid.top_layer_complete()); + } + + /// Fresh random key ids with Zipf weights: one trial's stream. The seed + /// is fixed, so each trial draws fresh ids instead; a uniformly random + /// hash makes those trials independent. + fn zipf_stream(n: usize, rng: &mut StdRng) -> Vec<(u64, i64)> { + let mut ids = HashSet::with_capacity(n); + while ids.len() < n { + ids.insert(rng.random::()); + } + ids.into_iter() + .enumerate() + .map(|(rank, id)| (id, (10_000.0 / ((rank + 1) as f64).powf(1.1)).ceil() as i64)) + .collect() + } + + /// Failures allowed out of `trials` at a true rate of at most `delta`: + /// the mean plus four standard deviations, plus one. + fn allowed_failures(trials: usize, delta: f64) -> f64 { + let n = trials as f64; + n * delta + 4.0 * (n * delta * (1.0 - delta)).sqrt() + 1.0 + } + + /// `calc_l2`'s documented bound, on a pyramid whose layer 0 is `5 x 64` + /// and whose mouse layers are `3 x 16`: + /// `p = 2/(64 · 0.25) = 0.125`, `δ = P(Bin(5, p) ≥ 3) ≈ 0.016`. + fn check_l2(trials: usize) { + const EPS: f64 = 0.5; + const DELTA: f64 = 0.0161; + let mut rng = StdRng::seed_from_u64(0x5EED_0013); + let mut failures = 0; + for _ in 0..trials { + let keys = zipf_stream(400, &mut rng); + let mut pyramid = UnivMonPyramid::new(8, 1, 5, 64, 3, 16, 3); + for &(id, weight) in &keys { + pyramid.insert(&DataInput::U64(id), weight); + } + let f2: f64 = keys.iter().map(|&(_, f)| (f as f64).powi(2)).sum(); + if (pyramid.calc_l2().powi(2) - f2).abs() > EPS * f2 { + failures += 1; + } + } + assert!( + failures as f64 <= allowed_failures(trials, DELTA), + "{failures}/{trials} misses beyond ε={EPS}, δ={DELTA}" + ); + } + + /// `calc_card`'s documented bound with a complete top layer, on mouse + /// layers much narrower than the elephant ones: unbiased, variance at most + /// `2 F0²/k`, Chebyshev `δ = 2/(k ε²)`. + fn check_card(trials: usize) { + const HEAP: usize = 32; + const KEYS: usize = 2000; + const EPS: f64 = 0.5; + // L ≥ ⌈log2(F0 / k)⌉ + 2. + let layers = (KEYS as f64 / HEAP as f64).log2().ceil() as usize + 2; + let mut rng = StdRng::seed_from_u64(0x5EED_00F1); + let f0 = KEYS as f64; + let (mut failures, mut kept, mut sum, mut sum_sq) = (0, 0, 0.0, 0.0); + for _ in 0..trials { + let mut pyramid = UnivMonPyramid::new(HEAP, 2, 3, 256, 2, 32, layers); + for (id, weight) in zipf_stream(KEYS, &mut rng) { + pyramid.insert(&DataInput::U64(id), weight); + } + if !pyramid.top_layer_complete() { + continue; + } + kept += 1; + let error = pyramid.calc_card() - f0; + sum += error; + sum_sq += error * error; + if error.abs() > EPS * f0 { + failures += 1; + } + } + // (e/4)^{k/2} ≈ 2e-3 for k = 32. + let incomplete = trials - kept; + assert!( + incomplete as f64 <= allowed_failures(trials, 0.01), + "top layer incomplete in {incomplete}/{trials} trials" + ); + let n = kept as f64; + let delta = 2.0 / (HEAP as f64 * EPS * EPS); + assert!( + failures as f64 <= allowed_failures(kept, delta), + "{failures}/{kept} misses beyond ε={EPS}" + ); + let variance = sum_sq / n; + assert!( + variance <= 2.0 * f0 * f0 / HEAP as f64, + "variance {:.3}·F0²/k above 2", + variance / (f0 * f0 / HEAP as f64) + ); + // Unbiased: the mean error is within four standard errors. + assert!( + (sum / n).abs() <= 4.0 * (variance / n).sqrt(), + "mean error {:.1} looks biased", + sum / n + ); + } + + #[test] + fn pyramid_calc_l2_meets_its_documented_bound() { + check_l2(150); + } + + #[test] + #[ignore = "long Monte Carlo run; the fast variant runs by default"] + fn pyramid_calc_l2_meets_its_documented_bound_long() { + check_l2(3000); + } + + #[test] + fn pyramid_calc_card_meets_its_documented_bound() { + check_card(60); + } + + #[test] + #[ignore = "long Monte Carlo run; the fast variant runs by default"] + fn pyramid_calc_card_meets_its_documented_bound_long() { + check_card(2000); + } #[test] fn pool_basic_take_put() { diff --git a/tests/e2e/frameworks.rs b/tests/e2e/frameworks.rs index 66185857..af75299a 100644 --- a/tests/e2e/frameworks.rs +++ b/tests/e2e/frameworks.rs @@ -1529,7 +1529,11 @@ fn univmon_weighted_metrics_and_fast_insert_parity() { #[test] fn univmon_pyramid_weighted_metrics() { - let mut up = UnivMonPyramid::with_defaults(); + // The default tiers with heap 128 instead of 32. With a complete top layer + // the distinct count is unbiased with variance at most `2 F0² / heap`: at + // heap 32 the 20% band is under one of those standard deviations, at heap + // 128 it is over 1.6 (and the exact identity's spread is ~0.6x that bound). + let mut up = UnivMonPyramid::new(128, 8, 3, 2048, 3, 512, 16); let mut truth = common::FreqTruth::default(); let stream = zipf_u64(15_000, 800, 1.3, 4004); for (i, k) in stream.iter().enumerate() { @@ -1545,6 +1549,10 @@ fn univmon_pyramid_weighted_metrics() { let l2_truth = truth.l2_norm(); assert_between(up.calc_l2(), l2_truth * 0.85, l2_truth * 1.15, "pyramid L2"); let card_truth = truth.distinct() as f64; + assert!( + up.top_layer_complete(), + "the top layer must hold every key it received for the bound to apply" + ); assert_between( up.calc_card(), card_truth * 0.80,