diff --git a/CHANGELOG.md b/CHANGELOG.md index 32c8aec3..35a3df29 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,10 @@ All significant changes to this project will be documented in this file. ## Unreleased +### New features + +* Add `UddSketch` behind the `uddsketch` feature for bounded quantile summaries of finite signed values, with weighted updates, compatible merges of sketches or serialized states, a queryable final relative-value error bound, and Rust-specific serialization. With identical parameters and mapping versions, states built through updates, merges, and unmodified serialization roundtrips have the same estimates and serialized bytes regardless of input order, partitioning, or merge tree. + ### Improvements * `BloomFilter::insert` is faster for streams of new items when multiple hashes are configured. diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index c02d8d16..b36f6f2d 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -75,6 +75,13 @@ To run only one sketch or workload, pass a Divan filter directly: cargo bench --package benchmarks --bench benchmarks -- cpc::serde ``` +The reproducible quantile comparison includes UDDSketch, the current t-digest, and a pinned t-digest 0.5.0 baseline. It prints tagged CSV with timing, memory estimates, serialized sizes, value and rank errors, and permutation/partition/merge-tree stability. Inputs are synthetic; the experiment excludes database scheduling and I/O. Run it with: + +```shell +cargo run --release -p benchmarks --example quantile_comparison +QUANTILE_ROWS=1000000 cargo run --release -p benchmarks --example quantile_comparison +``` + ## Public API documentation - Describe types with noun phrases and API behavior with third-person present-tense verbs such as `Creates`, `Updates`, and `Returns`. @@ -119,7 +126,7 @@ When adding a case to an existing sketch target, add it to the appropriate modul Cargo automatically discovers `tests-integration/tests/serde_tests.rs`, which aggregates the sketch-specific modules under `tests-integration/tests/serde_tests`. -To add serialization tests for another sketch, add `serde_tests/.rs` and its module declaration in `serde_tests.rs`. Shared path handling belongs in `serde_tests.rs`, and serialization fixtures belong in the appropriate subdirectory under `serde_tests`. +To add serialization tests for another sketch, add `serde_tests/.rs` and its module declaration in `serde_tests.rs`. UDDSketch uses local Rust-format fixtures in `serde_tests/uddsketch.rs`; it has no cross-language TCK snapshots. Shared path handling belongs in `serde_tests.rs`, and serialization fixtures belong in the appropriate subdirectory under `serde_tests`. ## Manual workflow (without xtask) diff --git a/Cargo.lock b/Cargo.lock index 81862f6d..0f496920 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -83,7 +83,8 @@ checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" name = "benchmarks" version = "0.0.0" dependencies = [ - "datasketches", + "datasketches 0.5.0", + "datasketches 0.6.0", "divan", "rand", ] @@ -238,6 +239,12 @@ dependencies = [ "libc", ] +[[package]] +name = "datasketches" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "11c0bd7d22989969a619bae09147b992a6aa7f31dd4f9d9af6345f533744781f" + [[package]] name = "datasketches" version = "0.6.0" @@ -778,7 +785,7 @@ dependencies = [ name = "tests-integration" version = "0.0.0" dependencies = [ - "datasketches", + "datasketches 0.6.0", "googletest", "insta", "quickcheck", diff --git a/README.md b/README.md index 7ec87caf..e606eb8f 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,7 @@ Enable multiple algorithms by listing their features together, such as `features | `tdigest` | `TDigestMut`, `TDigest` | Quantile and rank estimation, with high accuracy near distribution tails. | | `theta` | `ThetaSketch` and set operations | Distinct counts, set expressions, and Jaccard similarity. | | `tuple` | `TupleSketch` and set operations | Theta-style keys with user-defined summaries attached to retained entries. | +| `uddsketch` | `UddSketch` | Reproducible quantiles with bounded bucket count and a queryable relative-value error bound. | See the [API documentation](https://docs.rs/datasketches) for configuration, accuracy guarantees, serialization, and examples for each algorithm. Runnable end-to-end scenarios live in [`datasketches/examples`](datasketches/examples). @@ -82,7 +83,7 @@ See the [API documentation](https://docs.rs/datasketches) for configuration, acc The minimum supported Rust version is 1.86.0. The crate currently supports little-endian targets only. -Supported serialization formats are tested with fixtures produced by Apache DataSketches Java, C++, and Go through the [DataSketches TCK](https://github.com/apache/datasketches-tck). +UDDSketch has a Rust-specific versioned format without cross-language compatibility. Other supported serialization formats are tested with fixtures produced by Apache DataSketches Java, C++, and Go through the [DataSketches TCK](https://github.com/apache/datasketches-tck). Serialization compatibility does not imply that an ordinary Rust `Hash` implementation produces the same update bytes across languages. When sketches must represent the same inputs across implementations, use `hash::value::{raw_bytes, canonical_float, sign_extend, natural_extend}` (and the constructors within those modules) to select matching input encodings. When interoperating with an implementation that skips empty strings, skip them before updating here as well. diff --git a/benchmarks/Cargo.toml b/benchmarks/Cargo.toml index 25af32a6..8ee0b483 100644 --- a/benchmarks/Cargo.toml +++ b/benchmarks/Cargo.toml @@ -34,6 +34,11 @@ datasketches = { workspace = true, features = [ "tdigest", "theta", "tuple", + "uddsketch", +] } +# Pin the comparison baseline so results remain reproducible. +datasketches-v05 = { package = "datasketches", version = "=0.5.0", features = [ + "tdigest", ] } divan = { workspace = true } rand = { workspace = true } @@ -43,5 +48,9 @@ harness = false name = "benchmarks" path = "main.rs" +[[example]] +name = "quantile_comparison" +path = "uddsketch/comparison.rs" + [lints] workspace = true diff --git a/benchmarks/main.rs b/benchmarks/main.rs index ee288034..bbf6ff1f 100644 --- a/benchmarks/main.rs +++ b/benchmarks/main.rs @@ -31,6 +31,7 @@ mod req; mod tdigest; mod theta; mod tuple; +mod uddsketch; fn main() { divan::main(); diff --git a/benchmarks/uddsketch/comparison.rs b/benchmarks/uddsketch/comparison.rs new file mode 100644 index 00000000..1fb20816 --- /dev/null +++ b/benchmarks/uddsketch/comparison.rs @@ -0,0 +1,571 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +//! Reproducible, single-threaded quantile accuracy and cost experiment. +//! +//! Runs with `cargo run --release -p benchmarks --example quantile_comparison`. +//! `QUANTILE_ROWS` overrides the default 100,000 rows. Output is tagged CSV on stdout. +//! Timing is the median/min/max of seven samples after one warmup, with inputs prepared outside +//! the timer. Build includes finalization; merge includes decoding serialized partials and a +//! final query. Groups match 512 groups of 8/64 rows; shard counts are 1/8/64/256. +//! This measures sketch work, not database scheduling, networking, or storage I/O. + +use std::hint::black_box; +use std::time::Instant; + +use datasketches::tdigest::TDigestMut; +use datasketches::uddsketch::UddSketch; +use datasketches_v05::tdigest::TDigestMut as OldDigest; + +#[derive(Clone, Copy, Debug)] +enum Config { + Udd(u32, f64), + Current(u16), + Old(u16), +} + +impl Config { + fn name(self) -> String { + match self { + Self::Udd(m, a) => format!("udd_m{m}_a{a}"), + Self::Current(k) => format!("td06_k{k}"), + Self::Old(k) => format!("td05_k{k}"), + } + } + + fn empty(self) -> Sketch { + match self { + Self::Udd(m, a) => Sketch::Udd(UddSketch::new(m, a).unwrap()), + Self::Current(k) => Sketch::Current(TDigestMut::new(k).unwrap()), + Self::Old(k) => Sketch::Old(OldDigest::new(k).unwrap()), + } + } + + fn decode(self, bytes: &[u8]) -> Sketch { + match self { + Self::Udd(..) => Sketch::Udd(UddSketch::deserialize(bytes).unwrap()), + Self::Current(..) => Sketch::Current(TDigestMut::deserialize(bytes).unwrap()), + Self::Old(..) => Sketch::Old(OldDigest::deserialize(bytes).unwrap()), + } + } + + fn build(self, values: &[f64]) -> Sketch { + let mut sketch = self.empty(); + for &x in values { + sketch.update(x); + } + black_box(sketch.query(0.5)); + sketch + } +} + +#[derive(Clone)] +enum Sketch { + Udd(UddSketch), + Current(TDigestMut), + Old(OldDigest), +} + +impl Sketch { + fn update(&mut self, value: f64) { + match self { + Self::Udd(s) => s.update(value).unwrap(), + Self::Current(s) => s.update(value), + Self::Old(s) => s.update(value), + } + } + + fn merge(&mut self, other: &Self) { + match (self, other) { + (Self::Udd(s), Self::Udd(t)) => s.merge(t).unwrap(), + (Self::Current(s), Self::Current(t)) => s.merge(t), + (Self::Old(s), Self::Old(t)) => s.merge(t), + _ => unreachable!(), + } + } + + fn query(&mut self, q: f64) -> f64 { + match self { + Self::Udd(s) => s.quantile(q).unwrap().unwrap(), + Self::Current(s) => s.quantile(q).unwrap().unwrap(), + Self::Old(s) => s.quantile(q).unwrap(), + } + } + + fn bytes(&mut self) -> Vec { + match self { + Self::Udd(s) => s.serialize(), + Self::Current(s) => s.serialize(), + Self::Old(s) => s.serialize(), + } + } + + fn memory(&self) -> usize { + match self { + Self::Udd(s) => s.estimated_size(), + // The current API estimates the heap; include the inline state as in the other two. + Self::Current(s) => s.estimated_size() + size_of::(), + Self::Old(s) => s.estimated_size(), + } + } + + fn bound(&self) -> f64 { + match self { + Self::Udd(s) => s.relative_error(), + _ => f64::NAN, + } + } + + fn level(&self) -> u32 { + match self { + Self::Udd(s) => s.compaction_level(), + _ => 0, + } + } +} + +fn next(seed: &mut u64) -> u64 { + *seed ^= *seed << 13; + *seed ^= *seed >> 7; + *seed ^= *seed << 17; + *seed +} + +fn uniform(seed: &mut u64) -> f64 { + ((next(seed) >> 11) as f64 + 0.5) / ((1u64 << 53) as f64) +} + +fn shuffle(values: &mut [T], seed: &mut u64) { + for i in (1..values.len()).rev() { + values.swap(i, next(seed) as usize % (i + 1)); + } +} + +fn data(name: &str, n: usize) -> Vec { + let mut seed = 42; + (0..n) + .map(|i| { + let u = uniform(&mut seed); + match name { + "uniform" => u * 1000.0, + "lognormal" => (2.0 + * (-2.0 * u.ln()).sqrt() + * (std::f64::consts::TAU * uniform(&mut seed)).cos()) + .exp(), + "pareto" => (1.0 - u).powf(-1.0 / 1.2), + "outliers" => { + if i % 1000 == 0 { + 1e12 + } else { + 1.0 + u + } + } + "repeated" => (next(&mut seed) % 8) as f64, + "signed_zero" => { + if i % 4 == 0 { + 0.0 + } else { + (u - 0.5) * 2000.0 + } + } + "wide" => 10.0f64.powf(u * 600.0 - 300.0), + _ => unreachable!(), + } + }) + .collect() +} + +fn measure(mut run: impl FnMut()) -> (f64, f64, f64) { + run(); + let mut ns = Vec::new(); + for _ in 0..7 { + let start = Instant::now(); + run(); + ns.push(start.elapsed().as_secs_f64() * 1e9); + } + ns.sort_by(f64::total_cmp); + (ns[3], ns[0], ns[6]) +} + +fn timing(work: &str, cfg: Config, rows: usize, units: usize, result: (f64, f64, f64)) { + println!( + "time,{work},{},{rows},{units},{:.3},{:.3},{:.3}", + cfg.name(), + result.0 / units as f64, + result.1 / units as f64, + result.2 / units as f64 + ); +} + +const RANKS: [f64; 6] = [0.0, 0.5, 0.9, 0.99, 0.999, 1.0]; + +// A rounded floating product can cross an integer even for small sample sizes (e.g. 0.99 +// and 100,001 rows). Evaluate the binary rational rank exactly for the discrete reference. +fn lower_discrete_index(rank: f64, count: usize) -> usize { + let bits = rank.to_bits(); + let exponent = ((bits >> 52) & 0x7ff) as u32; + if exponent == 0 { + return 0; + } + let significand = (bits & ((1 << 52) - 1)) | (1 << 52); + let product = u128::from(significand) * (count - 1) as u128; + let shift = 1075 - exponent; + if shift >= 128 { + 0 + } else { + (product >> shift) as usize + } +} + +fn metrics(name: &str, cfg: Config, values: &[f64]) { + let mut sorted = values.to_vec(); + sorted.sort_by(f64::total_cmp); + let mut sketch = cfg.build(values); + let memory = sketch.memory(); + let bytes = sketch.bytes().len(); + println!( + "state,{name},{},{},{memory},{bytes},{},{:.10e}", + cfg.name(), + values.len(), + sketch.level(), + sketch.bound() + ); + for q in RANKS { + let discrete = sorted[lower_discrete_index(q, values.len())]; + let position = q * (values.len() - 1) as f64; + let lower = sorted[position.floor() as usize]; + let upper = sorted[position.ceil() as usize]; + let fraction = position.fract(); + let continuous = lower * (1.0 - fraction) + upper * fraction; + let estimate = sketch.query(q); + let relative = if discrete == 0.0 { + if estimate == 0.0 { 0.0 } else { f64::NAN } + } else { + ((estimate - discrete) / discrete).abs() + }; + let continuous_relative = if continuous == 0.0 { + if estimate == 0.0 { 0.0 } else { f64::NAN } + } else { + ((estimate - continuous) / continuous).abs() + }; + // Ties occupy a rank interval. Distance to that interval avoids treating an exact repeated + // value as erroneous solely because F(x) jumps at it. + let below = sorted.partition_point(|x| *x < estimate) as f64 / values.len() as f64; + let at_or_below = sorted.partition_point(|x| *x <= estimate) as f64 / values.len() as f64; + let rank_error = (below - q).max(q - at_or_below).max(0.0); + println!( + "accuracy,{name},{},{q},{estimate:.10e},{discrete:.10e},{continuous:.10e},{relative:.10e},{continuous_relative:.10e},{rank_error:.10e}", + cfg.name() + ); + } +} + +fn partials(cfg: Config, values: &[f64], shards: usize, layout: usize) -> Vec> { + let mut states = (0..shards).map(|_| cfg.empty()).collect::>(); + for (i, &x) in values.iter().enumerate() { + let shard = if layout == 0 { + i * shards / values.len() + } else { + i % shards + }; + states[shard].update(x); + } + states.iter_mut().map(Sketch::bytes).collect() +} + +fn reduce(cfg: Config, partials: &[Vec], tree: usize, seed: &mut u64) -> Sketch { + let mut states = partials + .iter() + .map(|bytes| cfg.decode(bytes)) + .collect::>(); + if tree == 1 { + states.reverse(); + } + if tree >= 2 { + shuffle(&mut states, seed); + } + if tree == 3 { + while states.len() > 1 { + let mut next_level = Vec::new(); + for pair in states.chunks(2) { + let mut a = pair[0].clone(); + if pair.len() == 2 { + a.merge(&pair[1]); + } + next_level.push(cfg.decode(&a.bytes())); + } + states = next_level; + } + states.pop().unwrap() + } else { + let mut result = cfg.empty(); + for s in states { + result.merge(&s); + } + result + } +} + +fn stability(name: &str, cfg: Config, values: &[f64]) { + let mut base = cfg.build(values); + let reference = RANKS.map(|q| base.query(q)); + let base_bytes = base.bytes(); + let mut min = reference; + let mut max = reference; + let mut different = [0; RANKS.len()]; + let mut states_different = 0; + let mut cases = 0; + let mut seed = 2026; + let mut shuffled = values.to_vec(); + shuffle(&mut shuffled, &mut seed); + let mut sorted = values.to_vec(); + sorted.sort_by(f64::total_cmp); + let reversed = sorted.iter().rev().copied().collect::>(); + let orders = [values, &shuffled, &sorted, &reversed]; + for order in orders { + for shards in [1, 8, 64, 256] { + for layout in 0..2 { + let parts = partials(cfg, order, shards, layout); + for tree in 0..4 { + let mut result = reduce(cfg, &parts, tree, &mut seed); + let answers = RANKS.map(|q| result.query(q)); + states_different += usize::from(result.bytes() != base_bytes); + for i in 0..RANKS.len() { + different[i] += usize::from(answers[i].to_bits() != reference[i].to_bits()); + min[i] = min[i].min(answers[i]); + max[i] = max[i].max(answers[i]); + } + cases += 1; + } + } + } + } + for i in 0..RANKS.len() { + println!( + "stability,{name},{},{},{cases},{},{states_different},{:.10e},{:.10e}", + cfg.name(), + RANKS[i], + different[i], + min[i], + max[i] + ); + } + // Isolate merge order/tree variation with exactly the same partial bytes. + let parts = partials(cfg, values, 64, 0); + let mut first = reduce(cfg, &parts, 0, &mut seed); + let expected = RANKS.map(|q| first.query(q)); + let mut changed = 0; + for tree in 1..4 { + let mut result = reduce(cfg, &parts, tree, &mut seed); + changed += usize::from(RANKS.map(|q| result.query(q)) != expected); + } + println!("fixed_partials,{name},{},3,{changed}", cfg.name()); +} + +fn costs(cfg: Config, values: &[f64]) { + timing( + "insert", + cfg, + values.len(), + values.len(), + measure(|| { + let mut sketch = cfg.empty(); + for &value in black_box(values) { + sketch.update(value); + } + black_box(sketch); + }), + ); + + timing( + "build_finalize", + cfg, + values.len(), + values.len(), + measure(|| { + black_box(cfg.build(black_box(values))); + }), + ); + let mut sketch = cfg.build(values); + timing( + "query_p99", + cfg, + values.len(), + 10_000, + measure(|| { + for _ in 0..10_000 { + black_box(sketch.query(black_box(0.99))); + } + }), + ); + timing( + "serialize", + cfg, + values.len(), + 100, + measure(|| { + for _ in 0..100 { + black_box(sketch.bytes()); + } + }), + ); + let bytes = sketch.bytes(); + timing( + "deserialize", + cfg, + values.len(), + 100, + measure(|| { + for _ in 0..100 { + black_box(cfg.decode(black_box(&bytes))); + } + }), + ); + for shards in [1, 8, 64, 256] { + let parts = partials(cfg, values, shards, 0); + timing( + &format!("merge_serialized_{shards}"), + cfg, + values.len(), + 1, + measure(|| { + let mut result = reduce(cfg, black_box(&parts), 0, &mut 42); + black_box(result.query(0.99)); + }), + ); + let decoded = parts + .iter() + .map(|bytes| cfg.decode(bytes)) + .collect::>(); + timing( + &format!("merge_memory_{shards}"), + cfg, + values.len(), + 1, + measure(|| { + let mut result = cfg.empty(); + for part in black_box(&decoded) { + result.merge(part); + } + black_box(result.query(0.99)); + }), + ); + println!( + "wire,{},{},{},{}", + cfg.name(), + values.len(), + shards, + parts.iter().map(Vec::len).sum::() + ); + } + for rows in [8, 64] { + let inputs = (0..512) + .map(|group| { + (0..rows) + .map(|row| row as f64 + group as f64 / 512.0) + .collect::>() + }) + .collect::>(); + timing( + "groups_512", + cfg, + rows, + 512 * rows, + measure(|| { + for input in &inputs { + black_box(cfg.build(black_box(input)).bytes()); + } + }), + ); + timing( + "partial_final_64", + cfg, + rows, + 64 * rows, + measure(|| { + let mut merged = cfg.empty(); + for group in 0..64 { + let mut partial = cfg.empty(); + for row in 0..rows { + partial.update(row as f64 + group as f64 / 64.0); + } + merged.merge(&cfg.decode(&partial.bytes())); + } + black_box(merged.query(0.9)); + }), + ); + let mut states = inputs.iter().map(|v| cfg.build(v)).collect::>(); + let memory = states.iter().map(Sketch::memory).sum::(); + let bytes = states.iter_mut().map(|s| s.bytes().len()).sum::(); + println!("groups_state,{},{rows},512,{memory},{bytes}", cfg.name()); + } +} + +fn main() { + let n = std::env::var("QUANTILE_ROWS") + .ok() + .map(|x| x.parse().unwrap()) + .unwrap_or(100_000); + assert!(n >= 1024); + println!( + "# seed=42; rows={n}; timing=median/min/max of 7 after warmup; nan means undefined/not provided" + ); + println!("# time,work,config,rows,units,median_ns_per_unit,min_ns_per_unit,max_ns_per_unit"); + println!( + "# state,dataset,config,rows,estimated_memory_bytes,serialized_bytes,level,relative_error_bound" + ); + println!( + "# accuracy,dataset,config,q,estimate,lower_discrete,type7_continuous,relative_value_error,continuous_relative_error,rank_interval_distance" + ); + println!( + "# stability,dataset,config,q,cases,changed_answers,changed_states,min_estimate,max_estimate" + ); + let configs = [ + Config::Udd(128, 0.01), + Config::Udd(512, 0.01), + Config::Udd(2048, 0.001), + Config::Current(100), + Config::Current(200), + Config::Current(500), + Config::Old(200), + ]; + for name in [ + "uniform", + "lognormal", + "pareto", + "outliers", + "repeated", + "signed_zero", + "wide", + ] { + let values = data(name, n); + for cfg in configs { + metrics(name, cfg, &values); + } + for cfg in [ + Config::Udd(512, 0.01), + Config::Current(200), + Config::Old(200), + ] { + stability(name, cfg, &values); + } + } + let values = data("lognormal", n); + for cfg in configs { + costs(cfg, &values); + } +} diff --git a/benchmarks/uddsketch/merge.rs b/benchmarks/uddsketch/merge.rs new file mode 100644 index 00000000..06c87098 --- /dev/null +++ b/benchmarks/uddsketch/merge.rs @@ -0,0 +1,65 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::uddsketch::UddSketch; +use divan::Bencher; +use divan::black_box; + +use super::support::build; +use super::support::values; + +#[divan::bench(args = [1, 8, 64, 256])] +fn partials(bencher: Bencher, shards: usize) { + let values = values(131_072); + let partials = values + .chunks(values.len() / shards) + .map(build) + .collect::>(); + bencher.bench_local(|| { + let mut merged = UddSketch::default(); + for partial in black_box(&partials) { + merged.merge(partial).unwrap(); + } + black_box(merged) + }); +} + +#[divan::bench(consts = [false, true], args = [8, 64, 16_384])] +fn serialized_partials(bencher: Bencher, rows: usize) { + let partials = (0..64) + .map(|group| { + let mut sketch = UddSketch::default(); + for row in 0..rows { + sketch.update(row as f64 + group as f64 / 64.0).unwrap(); + } + sketch.serialize() + }) + .collect::>(); + bencher.bench_local(|| { + let mut merged = UddSketch::default(); + for bytes in black_box(&partials) { + if DIRECT { + merged.merge_serialized(bytes).unwrap(); + } else { + merged + .merge(&UddSketch::deserialize(bytes).unwrap()) + .unwrap(); + } + } + black_box(merged) + }); +} diff --git a/benchmarks/uddsketch/mod.rs b/benchmarks/uddsketch/mod.rs new file mode 100644 index 00000000..4ff11ff9 --- /dev/null +++ b/benchmarks/uddsketch/mod.rs @@ -0,0 +1,22 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +mod merge; +mod query; +mod serde; +mod support; +mod update; diff --git a/benchmarks/uddsketch/query.rs b/benchmarks/uddsketch/query.rs new file mode 100644 index 00000000..ce9b23e2 --- /dev/null +++ b/benchmarks/uddsketch/query.rs @@ -0,0 +1,28 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use divan::Bencher; +use divan::black_box; + +use super::support::build; +use super::support::values; + +#[divan::bench] +fn quantile(bencher: Bencher) { + let sketch = build(&values(100_000)); + bencher.bench_local(|| black_box(&sketch).quantile(black_box(0.99)).unwrap()); +} diff --git a/benchmarks/uddsketch/serde.rs b/benchmarks/uddsketch/serde.rs new file mode 100644 index 00000000..1eae038c --- /dev/null +++ b/benchmarks/uddsketch/serde.rs @@ -0,0 +1,45 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::uddsketch::UddSketch; +use divan::Bencher; +use divan::black_box; + +use super::support::build; +use super::support::values; + +#[divan::bench] +fn serialize(bencher: Bencher) { + let sketch = build(&values(100_000)); + bencher.bench_local(|| black_box(&sketch).serialize()); +} + +#[divan::bench] +fn serialize_reusing_buffer(bencher: Bencher) { + let sketch = build(&values(100_000)); + let mut bytes = Vec::new(); + bencher.bench_local(|| { + black_box(&sketch).serialize_into(black_box(&mut bytes)); + black_box(bytes.len()) + }); +} + +#[divan::bench] +fn deserialize(bencher: Bencher) { + let bytes = build(&values(100_000)).serialize(); + bencher.bench_local(|| UddSketch::deserialize(black_box(&bytes)).unwrap()); +} diff --git a/benchmarks/uddsketch/support.rs b/benchmarks/uddsketch/support.rs new file mode 100644 index 00000000..c45c222a --- /dev/null +++ b/benchmarks/uddsketch/support.rs @@ -0,0 +1,32 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::uddsketch::UddSketch; + +pub fn values(n: usize) -> Vec { + (0..n) + .map(|i| (((i * 7919) % n) as f64 / n as f64 * 12.0).exp()) + .collect() +} + +pub fn build(values: &[f64]) -> UddSketch { + let mut sketch = UddSketch::default(); + for &value in values { + sketch.update(value).unwrap(); + } + sketch +} diff --git a/benchmarks/uddsketch/update.rs b/benchmarks/uddsketch/update.rs new file mode 100644 index 00000000..e44a2eed --- /dev/null +++ b/benchmarks/uddsketch/update.rs @@ -0,0 +1,50 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use divan::Bencher; +use divan::black_box; +use divan::counter::ItemsCount; + +use super::support::build; +use super::support::values; + +#[divan::bench(args = [1_000, 100_000])] +fn update(bencher: Bencher, n: usize) { + let values = values(n); + bencher + .counter(ItemsCount::new(n)) + .bench_local(|| build(black_box(&values))); +} + +#[divan::bench(args = [8, 64])] +fn partial_groups(bencher: Bencher, rows: usize) { + let groups = (0..512) + .map(|group| { + (0..rows) + .map(|row| row as f64 + group as f64 / 512.0) + .collect::>() + }) + .collect::>(); + bencher + .counter(ItemsCount::new(512 * rows)) + .bench_local(|| { + groups + .iter() + .map(|values| build(black_box(values))) + .collect::>() + }); +} diff --git a/datasketches/Cargo.toml b/datasketches/Cargo.toml index 2c931fc7..272f13d7 100644 --- a/datasketches/Cargo.toml +++ b/datasketches/Cargo.toml @@ -48,6 +48,7 @@ req = [] tdigest = [] theta = [] tuple = [] +uddsketch = [] [[example]] name = "bloom_difference" diff --git a/datasketches/src/lib.rs b/datasketches/src/lib.rs index 813a841b..516c2c06 100644 --- a/datasketches/src/lib.rs +++ b/datasketches/src/lib.rs @@ -40,6 +40,7 @@ //! * Use `kll`, `req`, or `tdigest` for ranks and quantiles. KLL provides strong general-purpose //! rank accuracy, REQ targets configurable high- or low-rank accuracy, and T-Digest emphasizes //! distribution tails. +//! * Use `uddsketch` for reproducible quantiles with a queryable relative-value error bound. //! * Use `tuple` when retained Theta keys need application-defined summaries. //! //! See each module's documentation for accuracy, memory, serialization, and update examples. @@ -92,6 +93,9 @@ pub use self::thetafamily::theta; #[cfg(feature = "tuple")] pub use self::thetafamily::tuple; +#[cfg(feature = "uddsketch")] +pub mod uddsketch; + // common modules pub mod codec; pub mod common; diff --git a/datasketches/src/uddsketch/mapping.rs b/datasketches/src/uddsketch/mapping.rs new file mode 100644 index 00000000..0eec94f2 --- /dev/null +++ b/datasketches/src/uddsketch/mapping.rs @@ -0,0 +1,361 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +// The UDD1 mapping is defined by these operations and round-to-nearest f64 arithmetic, +// not by a platform libm. All nonzero floating intermediates are normal. +// Keep its constants and rounding rules stable: they determine serialized bucket keys. +const LOG_UNIT: i64 = 1 << 52; +const MIN_LOG: i64 = -1074 * LOG_UNIT; +// floor(log2((16 + i) / 16) * 2^52), i = 0..=16. +const ANCHORS: [u64; 17] = [ + 0, + 393_897_619_263_793, + 765_274_173_176_528, + 1_116_566_257_159_425, + 1_449_835_248_174_818, + 1_766_840_599_037_378, + 2_069_096_066_497_155, + 2_357_913_430_203_731, + 2_634_436_900_273_512, + 2_899_670_496_349_637, + 3_154_500_053_615_717, + 3_399_711_073_450_041, + 3_636_003_326_134_361, + 3_864_002_889_947_501, + 4_084_272_148_448_331, + 4_297_318_147_896_633, + 4_503_599_627_370_496, +]; + +#[derive(Clone, Copy, Debug)] +pub struct Mapping { + step: u64, + reciprocal: u64, + offset: u64, +} + +impl Mapping { + pub fn new(accuracy: f64) -> Self { + let step = log2_fixed(1.0 + accuracy) - log2_fixed(1.0 - accuracy); + // ceil(L / step) - ceil(MIN_LOG / step) is floor((L + offset) / step). + // MIN_LOG is negative; offset makes the numerator nonnegative over the f64 domain. + // For supported accuracies, 0 < step <= 54 * LOG_UNIT, so the numerator is below + // (1074 + 1024 + 54) * LOG_UNIT < 2^64 even though it can exceed i64::MAX. + let step = step as u64; + Self { + step, + reciprocal: u64::MAX / step, + offset: MIN_LOG.unsigned_abs() + step - 1 - MIN_LOG.unsigned_abs() % step, + } + } + + pub fn index(self, magnitude: f64) -> u64 { + // The mathematical sum is in [0, 2^64); wrapping addition represents it even when L < 0. + self.divide((log2_fixed(magnitude) as u64).wrapping_add(self.offset)) + } + + fn divide(self, numerator: u64) -> u64 { + // reciprocal / 2^64 underestimates 1 / step by at most 2^-64. + // For a u64 numerator, the resulting quotient is exact or one too small. + // Its product with step cannot exceed numerator; one remainder check gives exact division. + let quotient = ((u128::from(numerator) * u128::from(self.reciprocal)) >> 64) as u64; + quotient + u64::from(numerator - quotient * self.step >= self.step) + } + + pub fn max_level(self) -> u32 { + 64 - self.index(f64::MAX).leading_zeros() + } +} + +// Approximates log2(x) in signed Q52 for every positive finite f64, including subnormals. +// Normalize exactly to x = 2^exponent * significand / 2^52. Within each interval, +// z = (m - anchor) / (m + anchor) is in [0, 1/33]. The identity +// log2(m / anchor) = 2/ln(2) * (z + z^3/3 + z^5/5 + ...) gives a positive polynomial. +// Five terms leave a truncation error below 5.3e-18 before arithmetic rounding. +// +// Monotonicity also holds after rounding: the numerator conversion is exact; increasing +// significand by one increases the rounded denominator by at most two, and the denominator +// exceeds twice the numerator. Thus even the rounded-operand ratio is nondecreasing. +// Correctly rounded division, positive polynomial operations, and integer truncation preserve +// that order. Capping at the next anchor joins intervals monotonically, including powers of two. +// Integer exponent assembly avoids cancellation near one. No transcendental calls are used. +fn log2_fixed(value: f64) -> i64 { + let bits = value.to_bits(); + let raw_exponent = (bits >> 52) as i64; + let (exponent, significand) = if raw_exponent == 0 { + let top = 63 - bits.leading_zeros(); + (i64::from(top) - 1074, bits << (52 - top)) + } else { + (raw_exponent - 1023, (bits & ((1 << 52) - 1)) | (1 << 52)) + }; + // Normalization puts significand in [2^52, 2^53); masking also makes the table bounds + // explicit to the compiler without changing the segment selected by UDD1. + let segment = ((significand >> 48) & 15) as usize; + // Scaling the original integer numerator and denominator by 2^-52 commutes with + // round-to-nearest: all nonzero operands/results are normal, and subtraction is exact. + // This produces the identical z bits while avoiding integer-to-float conversions. + let normalized = f64::from_bits((1023u64 << 52) | (significand & ((1 << 52) - 1))); + let anchor = f64::from_bits((1023u64 << 52) | ((segment as u64) << 48)); + let z = (normalized - anchor) / (normalized + anchor); + let square = z * z; + let mut polynomial = 1.0 / 9.0; + for coefficient in [1.0 / 7.0, 1.0 / 5.0, 1.0 / 3.0, 1.0] { + polynomial = coefficient + square * polynomial; + } + let offset = (z * polynomial * (2.0 * std::f64::consts::LOG2_E * LOG_UNIT as f64)) as u64; + let fraction = (ANCHORS[segment] + offset).min(ANCHORS[segment + 1]); + exponent * LOG_UNIT + fraction as i64 +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn monotone_at_segment_and_exponent_boundaries() { + let mut boundaries = vec![1, (1 << 52) - 1, 1 << 52, f64::MAX.to_bits()]; + for exponent in 1..2047u64 { + for segment in 0..16 { + boundaries.push((exponent << 52) | (segment << 48)); + } + } + for top in 0..52 { + boundaries.push(1 << top); + } + for bits in boundaries { + let lo = bits.saturating_sub(8).max(1); + let hi = (bits + 8).min(f64::MAX.to_bits()); + for bits in lo..hi { + assert!(log2_fixed(f64::from_bits(bits)) <= log2_fixed(f64::from_bits(bits + 1))); + } + } + } + + #[test] + fn monotone_interior_and_mapping_range() { + let mut seed = 42u64; + let mappings = + [1e-9, 0.001, 0.01, 0.5, f64::from_bits(1.0f64.to_bits() - 1)].map(Mapping::new); + for mapping in mappings { + assert_eq!(mapping.index(f64::from_bits(1)), 0); + assert!(mapping.max_level() < 64); + } + for _ in 0..100_000 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + let bits = 1 + seed % (f64::MAX.to_bits() - 1); + let x = f64::from_bits(bits); + let next = f64::from_bits(bits + 1); + assert!(log2_fixed(x) <= log2_fixed(next)); + for mapping in mappings { + assert!(mapping.index(x) <= mapping.index(next)); + assert!(mapping.index(x) <= mapping.index(f64::MAX)); + } + } + } + + #[test] + fn reciprocal_division_is_exact_at_integer_boundaries() { + let mut seed = 42u64; + let configured_steps = [1e-9, 0.001, 0.01, 0.5, f64::from_bits(1.0f64.to_bits() - 1)] + .map(|accuracy| Mapping::new(accuracy).step); + let mut steps = (1..=255).collect::>(); + for exponent in 1..64 { + let power = 1u64 << exponent; + steps.extend([power - 1, power, power + 1]); + } + steps.extend(configured_steps); + steps.push(u64::MAX); + steps.sort_unstable(); + steps.dedup(); + for step in steps { + let mapping = Mapping { + step, + reciprocal: u64::MAX / step, + offset: 0, + }; + // Small numerators exercise every quotient and remainder in a bounded grid. + if step <= 255 { + for numerator in 0..4096 { + assert_eq!(mapping.divide(numerator), numerator / step); + } + } + let last_multiple = (u64::MAX / step) * step; + for numerator in [ + 0, + 1, + last_multiple.saturating_sub(1), + last_multiple, + last_multiple.saturating_add(1), + u64::MAX - 1, + mapping.step - 1, + mapping.step, + mapping.step.saturating_add(1), + u64::MAX, + ] { + assert_eq!(mapping.divide(numerator), numerator / mapping.step); + } + for _ in 0..256 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + assert_eq!(mapping.divide(seed), seed / mapping.step); + let multiple = (seed / mapping.step) * mapping.step; + for numerator in [ + multiple.saturating_sub(1), + multiple, + multiple.saturating_add(1), + ] { + assert_eq!(mapping.divide(numerator), numerator / mapping.step); + } + } + } + } + + #[test] + fn translated_indices_match_signed_ceiling() { + let mut seed = 42u64; + let mut values = vec![ + f64::from_bits(1), + f64::MIN_POSITIVE, + 0.5, + 1.0, + 2.0, + f64::MAX, + ]; + for _ in 0..4096 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + values.push(f64::from_bits(1 + seed % f64::MAX.to_bits())); + } + let ceil = |value: i128, step: i128| -(-value).div_euclid(step); + for accuracy in [1e-9, 0.001, 0.01, 0.5, f64::from_bits(1.0f64.to_bits() - 1)] { + let mapping = Mapping::new(accuracy); + for &value in &values { + // Signed wide arithmetic is independent of the optimized unsigned translation. + let coordinate = i128::from(log2_fixed(value)); + let numerator = coordinate + i128::from(mapping.offset); + assert!((0..=i128::from(u64::MAX)).contains(&numerator)); + let step = i128::from(mapping.step); + let expected = ceil(coordinate, step) - ceil(i128::from(MIN_LOG), step); + assert_eq!(i128::from(mapping.index(value)), expected); + } + } + } + + #[test] + fn scaled_logarithm_matches_udd1_integer_operands() { + // This reference retains the unscaled UDD1 operands and rounding sequence. + // Comparing Q52 coordinates also protects bucket boundaries for every accuracy. + let reference = |significand: u64| { + let segment = ((significand - (1 << 52)) >> 48) as usize; + let anchor = (16 + segment as u64) << 48; + let z = (significand - anchor) as f64 / (significand + anchor) as f64; + let square = z * z; + let polynomial = 1.0 + + square + * (1.0 / 3.0 + + square * (1.0 / 5.0 + square * (1.0 / 7.0 + square * (1.0 / 9.0)))); + let offset = + (z * polynomial * (2.0 * std::f64::consts::LOG2_E * LOG_UNIT as f64)) as u64; + (ANCHORS[segment] + offset).min(ANCHORS[segment + 1]) as i64 + }; + let check = |significand: u64| { + let value = f64::from_bits((1023 << 52) | (significand - (1 << 52))); + assert_eq!( + log2_fixed(value), + reference(significand), + "{significand:#x}" + ); + }; + for segment in 0..16 { + let first = (1 << 52) + (segment << 48); + let last = first + (1 << 48) - 1; + // Consecutive significands cover both directions of ties-to-even denominator rounding. + for delta in 0..4096 { + check(first + delta); + check(last - delta); + } + } + let mut seed = 42u64; + for _ in 0..100_000 { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + check((1 << 52) | (seed & ((1 << 52) - 1))); + } + for top in 0..=52 { + for bits in [ + (1u64 << top).saturating_sub(1).max(1), + 1 << top, + (1 << top) + 1, + ] { + let top = 63 - bits.leading_zeros(); + let expected = (i64::from(top) - 1074) * LOG_UNIT + reference(bits << (52 - top)); + assert_eq!(log2_fixed(f64::from_bits(bits)), expected); + } + } + } + + #[test] + fn udd1_bucket_boundary_fixtures() { + // Frozen from the original UDD1 mapping before reciprocal division or operand scaling. + // Check adjacent f64s at bucket transitions, where even a one-bit change matters. + let fixtures = [ + (1e-9, 0x000fffffff9d4163u64, 18_021_826_274u64), + (1e-9, 0x3fb99999994f48bdu64, 371_068_734_765u64), + (1e-9, 0x3feffffffeed1f43u64, 372_220_027_284u64), + (1e-9, 0x3ff3fffffff649e8u64, 372_331_599_058u64), + (1e-9, 0x7fefffffff43a5e6u64, 727_111_375_459u64), + (0.001, 0x000ff88fd6b7d652u64, 18_021u64), + (0.001, 0x3fefefa1e2d85c14u64, 372_219u64), + (0.001, 0x54b247a6d904b1c7u64, 487_349u64), + (0.001, 0x7feffc195daab813u64, 727_111u64), + (0.01, 0x000fe2a1ff8747cdu64, 1_802u64), + (0.01, 0x3fef5dc83cd4e932u64, 37_220u64), + (0.01, 0x3ff3eff7bc7a494du64, 37_232u64), + (0.01, 0x7fef65643b949679u64, 72_708u64), + (0.5, 0x000cfd1c8ffd9c2cu64, 33u64), + (0.5, 0x3fd5555555555557u64, 677u64), + (0.5, 0x3ff0000000000001u64, 678u64), + (0.5, 0x7fed906a378b57ffu64, 1_324u64), + (0.9999999999999999, 0x0001000000000001u64, 1u64), + (0.9999999999999999, 0x3ff0000000000001u64, 20u64), + (0.9999999999999999, 0x7cb0000000000001u64, 38u64), + ]; + for (accuracy, bits, index) in fixtures { + let mapping = Mapping::new(accuracy); + assert_eq!(mapping.index(f64::from_bits(bits - 1)), index - 1); + assert_eq!(mapping.index(f64::from_bits(bits)), index); + assert_eq!(mapping.index(f64::from_bits(bits + 1)), index); + } + } + + #[test] + fn powers_of_two_are_exact() { + for exponent in -1074..=1023 { + let bits = if exponent < -1022 { + 1 << (exponent + 1074) + } else { + ((exponent + 1023) as u64) << 52 + }; + assert_eq!(log2_fixed(f64::from_bits(bits)), exponent * LOG_UNIT); + } + } +} diff --git a/datasketches/src/uddsketch/mod.rs b/datasketches/src/uddsketch/mod.rs new file mode 100644 index 00000000..037f8127 --- /dev/null +++ b/datasketches/src/uddsketch/mod.rs @@ -0,0 +1,95 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +//! Uniformly collapsing logarithmic histograms for reproducible quantiles. +//! +//! [`UddSketch`] implements the insertion-only uniform-collapse strategy from +//! [UDDSketch](https://arxiv.org/abs/2004.08604), with a fixed integer bucket hierarchy. +//! Compatible sketches produce identical serialized states and estimates for the same multiset, +//! regardless of input order, partition boundaries, or merge tree. This guarantee assumes the +//! same mapping version and states built through updates, merges, and unmodified serialization +//! roundtrips. The initial mapping is deterministic and monotone. +//! +//! # Accuracy and quantile semantics +//! +//! Queries approximate the lower discrete quantile: sorted item `floor(q * (n - 1))`, using +//! zero-based indexing. There is no interpolation between ranks. Endpoints return exact extrema. +//! [`UddSketch::relative_error`] bounds relative **value** error against that item, not rank error +//! and not error against an interpolated percentile. A true zero quantile returns exactly zero. +//! +//! The initial accuracy is a target, not an unconditional bound. Exceeding the bucket budget +//! coarsens every bucket, including those near the median and tails. The final error bound is +//! computed from retained bucket extrema, including floating-point rounding and subnormal values. +//! Use [`UddSketch::relative_error`] for a guarantee even when no collapse has occurred. +//! +//! # Mapping and compatibility +//! +//! The initial index is `ceil(L(abs(x)) / (L(1 + alpha) - L(1 - alpha)))`, where `L` is +//! a fixed, monotone Q52 approximation to `log2` and `1 +/- alpha` uses `f64` arithmetic. +//! `L` normalizes the exponent and significand exactly, then evaluates a polynomial +//! with positive coefficients within sixteen significand intervals using correctly rounded basic +//! `f64` arithmetic. The resulting Q52 coordinates are combined and divided as integers. Interval +//! endpoints are capped to join monotonically. It does not depend on platform logarithms. Indices +//! are translated so that the smallest positive `f64` has index zero. Each collapse divides these +//! nonnegative indices by two using integer arithmetic. New values always traverse this same +//! hierarchy; they never use a new logarithm base after compaction. Negative values have a mirrored +//! hierarchy and zero has its own bucket. Translation lets all finite magnitudes eventually share +//! one bucket per sign, so three buckets suffice for any finite input. It changes the pairing +//! alignment from the paper's grid at one. +//! +//! Each bucket retains its count and observed extrema. Its representative is their harmonic mean +//! in magnitude, rounded and clamped to those extrema. Unlike representatives derived from powers +//! of gamma, this remains finite across the entire `f64` domain and permits a data-dependent error +//! certificate. Counts and extrema merge exactly; no order-dependent floating-point sum is stored. +//! +//! For insertion-only streams, the occupied bucket count at any fixed level never decreases. +//! A subset cannot require a coarser level than the complete multiset. Consequently every build +//! reaches the same smallest level that fits the budget. Integer parent mapping composes exactly, +//! and combining counts and extrema is associative and commutative, so merging partitions reaches +//! the same state. These arguments require identical initial parameters and exclude deletion. +//! +//! Serialization uses a versioned, Rust-specific format. The mapping constants and rounding +//! rules are part of this version's bucket semantics. Deserialization validates +//! histogram structure, not provenance or minimal compaction level. Externally constructed +//! images can represent coarser histograms and need not reproduce a direct build from original +//! inputs; see [`UddSketch::deserialize`]. +//! Independent groups can reuse the sketch's bucket allocation with [`UddSketch::reset`] and +//! the output buffer with [`UddSketch::serialize_into`]. +//! [`UddSketch::merge_serialized`] validates and merges an image in one pass when the levels +//! match, using already validated intervals to establish extrema membership. It preserves the +//! validation and result of deserializing and merging separately, and leaves the receiver +//! unchanged on every error. +//! +//! # Examples +//! +//! ``` +//! use datasketches::uddsketch::UddSketch; +//! +//! let mut sketch = UddSketch::new(512, 0.01)?; +//! for value in [1.0, 2.0, 3.0] { +//! sketch.update(value)?; +//! } +//! assert_eq!(sketch.quantile(0.5)?, Some(2.0)); +//! let restored = UddSketch::deserialize(&sketch.serialize())?; +//! sketch.merge(&restored)?; +//! assert_eq!(sketch.total_weight(), 6); +//! # Ok::<(), datasketches::error::Error>(()) +//! ``` + +mod mapping; +mod sketch; +pub use self::sketch::UddSketch; diff --git a/datasketches/src/uddsketch/sketch.rs b/datasketches/src/uddsketch/sketch.rs new file mode 100644 index 00000000..ff17da54 --- /dev/null +++ b/datasketches/src/uddsketch/sketch.rs @@ -0,0 +1,639 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use super::mapping::Mapping; +use crate::codec::SketchSlice; +use crate::error::Error; + +/// A bounded, mergeable histogram with reproducible relative-value quantile estimates. +/// +/// See the [module documentation](crate::uddsketch) for mapping and accuracy contracts. +#[derive(Clone, Debug)] +pub struct UddSketch { + max_buckets: u32, + initial_accuracy: f64, + mapping: Mapping, + level: u32, + count: u64, + buckets: Vec, +} + +#[derive(Clone, Copy, Debug)] +struct Bucket { + key: i64, + count: u64, + min: f64, + max: f64, +} + +impl Default for UddSketch { + /// Creates a sketch with at most `512` buckets and initial accuracy `0.01`. + fn default() -> Self { + Self::new(512, 0.01).unwrap() + } +} + +impl UddSketch { + /// Creates a sketch with a shared budget for negative, zero, and positive buckets. + /// + /// # Errors + /// + /// Returns `InvalidArgument` if `max_buckets < 3` or `initial_accuracy` is outside + /// `[1e-9, 1)`, including NaN. The lower accuracy limit bounds the resolution + /// required of the fixed-point logarithmic mapping. + pub fn new(max_buckets: u32, initial_accuracy: f64) -> Result { + if max_buckets < 3 || !(1e-9..1.0).contains(&initial_accuracy) { + return Err(Error::invalid_argument( + "max_buckets must be at least 3 and initial_accuracy must be in [1e-9, 1)", + )); + } + Ok(Self { + max_buckets, + initial_accuracy, + mapping: Mapping::new(initial_accuracy), + level: 0, + count: 0, + buckets: Vec::new(), + }) + } + + /// Removes all values and restores the initial compaction level, retaining the configuration + /// and allocated bucket capacity. + pub fn reset(&mut self) { + self.level = 0; + self.count = 0; + self.buckets.clear(); + } + + /// Inserts one finite value, treating both signed zeros as positive zero. + /// + /// # Errors + /// + /// Returns `InvalidArgument` for non-finite values or total weight overflow, leaving the + /// sketch unchanged. + pub fn update(&mut self, value: f64) -> Result<(), Error> { + self.update_weighted(value, 1) + } + + /// Inserts `weight` copies of a finite value. + /// + /// A zero weight has no effect, but the value is still validated. + /// + /// # Errors + /// + /// Returns `InvalidArgument` for non-finite values or total weight overflow, leaving the + /// sketch unchanged. + pub fn update_weighted(&mut self, value: f64, weight: u64) -> Result<(), Error> { + if !value.is_finite() { + return Err(Error::invalid_argument("value must be finite")); + } + let count = self + .count + .checked_add(weight) + .ok_or_else(|| Error::invalid_argument("total weight overflow"))?; + if weight == 0 { + return Ok(()); + } + let value = if value == 0.0 { 0.0 } else { value }; + // Monotonicity lets observed extrema locate the only possible bucket. A value inside + // an existing interval needs no mapping; a gap can match only its adjacent buckets. + // Ascending input stays at the tail, avoiding a binary search for each observation. + let position = match self.buckets.last() { + Some(last) if value >= last.min => self.buckets.len() - usize::from(value <= last.max), + _ => self.buckets.partition_point(|b| b.max < value), + }; + if let Some(bucket) = self.buckets.get_mut(position).filter(|b| b.min <= value) { + bucket.count += weight; + } else { + let key = self.key(value); + let bucket = Bucket { + key, + count: weight, + min: value, + max: value, + }; + if position < self.buckets.len() && self.buckets[position].key == key { + // The interval lookup already established value < this bucket's minimum. + self.buckets[position].count += weight; + self.buckets[position].min = value; + } else if position > 0 && self.buckets[position - 1].key == key { + // Its predecessor is the only other candidate, with maximum < value. + self.buckets[position - 1].count += weight; + self.buckets[position - 1].max = value; + } else { + self.buckets.insert(position, bucket); + self.enforce_budget(); + } + } + self.count = count; + Ok(()) + } + + /// Merges a compatible sketch without modifying the input. + /// + /// Both sketches must have identical bucket limits and initial accuracy, even when empty. + /// Different compaction levels are aligned automatically. + /// + /// # Errors + /// + /// Returns `InvalidArgument` for incompatible parameters or total weight overflow, leaving + /// the receiver unchanged. + pub fn merge(&mut self, other: &Self) -> Result<(), Error> { + if self.max_buckets != other.max_buckets + || self.initial_accuracy.to_bits() != other.initial_accuracy.to_bits() + { + return Err(Error::invalid_argument("incompatible UDDSketch parameters")); + } + let count = self + .count + .checked_add(other.count) + .ok_or_else(|| Error::invalid_argument("total weight overflow"))?; + if other.is_empty() { + return Ok(()); + } + if self.is_empty() { + self.buckets.clone_from(&other.buckets); + self.level = other.level; + self.count = count; + return Ok(()); + } + let level = self.level.max(other.level); + self.collapse_to(level); + let mut merged: Vec = Vec::with_capacity(self.buckets.len() + other.buckets.len()); + let mut left = self.buckets.as_slice(); + for mut bucket in other.buckets.iter().copied() { + bucket.key = parent(bucket.key, level - other.level); + while let Some((&first, rest)) = left.split_first() { + if first.key > bucket.key { + break; + } + left = rest; + if first.key == bucket.key { + bucket.merge(first); + break; + } + merged.push(first); + } + // Alignment can map several consecutive input buckets to the same parent. + if let Some(last) = merged.last_mut().filter(|last| last.key == bucket.key) { + last.merge(bucket); + } else { + merged.push(bucket); + } + } + merged.extend_from_slice(left); + self.buckets = merged; + self.count = count; + self.enforce_budget(); + Ok(()) + } + + /// Validates and merges a Rust-specific version 1 image. + /// + /// Produces the same state as [`Self::deserialize`] followed by [`Self::merge`]. + /// Already validated bucket intervals can establish extrema membership in the input image. + /// + /// # Errors + /// + /// Returns `InvalidData` for any image rejected by [`Self::deserialize`]. After validating + /// the image, returns `InvalidArgument` for incompatible parameters or total weight overflow. + /// Every error leaves the receiver unchanged. + pub fn merge_serialized(&mut self, bytes: &[u8]) -> Result<(), Error> { + let (other, payload) = Self::deserialize_header(bytes)?; + let count = self.count.checked_add(other.count); + if self.max_buckets != other.max_buckets + || self.initial_accuracy.to_bits() != other.initial_accuracy.to_bits() + || count.is_none() + { + return self.merge(&other.deserialize_body(payload)?); + } + if self.is_empty() { + *self = other.deserialize_body(payload)?; + return Ok(()); + } + if payload.is_empty() { + return Ok(()); + } + if self.level != other.level { + return self.merge(&other.deserialize_body(payload)?); + } + + let mut merged = Vec::with_capacity(self.buckets.len() + payload.len() / 32); + let mut left = self.buckets.as_slice(); + let mut previous = None; + let mut input_count = 0u64; + for chunk in payload.chunks_exact(32) { + let mut bucket = Bucket::read(chunk); + bucket.validate_after(previous)?; + // Bound every intermediate sum before combining it with receiver counts. The + // declared total may be corrupt even when it fits alongside the receiver's total. + input_count = input_count + .checked_add(bucket.count) + .filter(|&count| count <= other.count) + .ok_or_else(invalid_image)?; + previous = Some(bucket); + while let Some((&first, rest)) = left.split_first() { + if first.key >= bucket.key { + break; + } + merged.push(first); + left = rest; + } + let known = left.first().filter(|first| first.key == bucket.key); + self.validate_membership(&bucket, known)?; + if let Some(&first) = known { + bucket.merge(first); + left = &left[1..]; + } + merged.push(bucket); + } + if input_count != other.count { + return Err(invalid_image()); + } + merged.extend_from_slice(left); + // The receiver, including its compaction level, remains untouched until validation + // has covered every input bucket and reconciled the complete input weight. + self.buckets = merged; + self.count = count.unwrap(); + self.enforce_budget(); + Ok(()) + } + + /// Returns the number of inserted values, including weights. + pub fn total_weight(&self) -> u64 { + self.count + } + + /// Returns whether the sketch contains no values. + pub fn is_empty(&self) -> bool { + self.count == 0 + } + + /// Returns the configured bucket limit. + pub fn max_buckets(&self) -> u32 { + self.max_buckets + } + + /// Returns the initial relative accuracy target. + pub fn initial_relative_accuracy(&self) -> f64 { + self.initial_accuracy + } + + /// Returns the number of retained buckets, including the zero bucket if present. + pub fn num_buckets(&self) -> usize { + self.buckets.len() + } + + /// Returns the number of uniform collapses in the current mapping hierarchy. + pub fn compaction_level(&self) -> u32 { + self.level + } + + /// Returns the exact minimum, or `None` when empty. + pub fn min_value(&self) -> Option { + self.buckets.first().map(|b| b.min) + } + + /// Returns the exact maximum, or `None` when empty. + pub fn max_value(&self) -> Option { + self.buckets.last().map(|b| b.max) + } + + /// Returns a conservative bound on relative value error for every quantile of this state. + /// + /// For a nonzero lower discrete quantile `x`, the estimate `y` satisfies + /// `abs(y - x) / abs(x) <= relative_error()`. Zero is returned exactly. Empty sketches and + /// sketches whose buckets each contain a single distinct value return `0.0`. + /// The bound may exceed the initial target after uniform collapse. It uses observed extrema + /// and outward rounding, rather than assuming exact logarithms or representable gamma powers. + /// This method scans the retained buckets; it does not bound rank error. + pub fn relative_error(&self) -> f64 { + self.buckets + .iter() + .map(|b| { + let estimate = b.estimate().abs(); + if b.min == b.max { + return 0.0; + } + let lo = b.min.abs().min(b.max.abs()); + let hi = b.min.abs().max(b.max.abs()); + // Round each nonnegative operation upward; this also covers subnormal arithmetic. + let low_error = relative_difference(estimate, lo, lo); + // The upper-end error cannot exceed one; clamping also handles upward + // rounding of a difference near f64::MAX to infinity. + let high_error = relative_difference(hi, estimate, hi).min(1.0); + low_error.max(high_error) + }) + .fold(0.0, f64::max) + } + + /// Returns an estimate of the lower discrete quantile at `rank`, or `None` when empty. + /// + /// The target is sorted item `floor(rank * (n - 1))` (zero-based), with the multiplication + /// evaluated exactly for the supplied binary floating-point rank. Endpoints are exact. + /// + /// # Errors + /// + /// Returns `InvalidArgument` if `rank` is outside `[0, 1]`, including NaN, even when empty. + pub fn quantile(&self, rank: f64) -> Result, Error> { + validate_rank(rank)?; + if self.is_empty() { + return Ok(None); + } + if rank == 0.0 { + return Ok(self.min_value()); + } + if rank == 1.0 { + return Ok(self.max_value()); + } + let target = quantile_index(rank, self.count); + let mut cumulative = 0; + for bucket in &self.buckets { + cumulative += bucket.count; + if cumulative > target { + return Ok(Some(bucket.estimate())); + } + } + unreachable!("validated counts cover the target rank") + } + + /// Returns estimates for all ranks in input order, or `None` when empty. + /// + /// # Errors + /// + /// Returns `InvalidArgument` if any rank is outside `[0, 1]`, including NaN, even when empty. + pub fn quantiles(&self, ranks: &[f64]) -> Result>, Error> { + for &rank in ranks { + validate_rank(rank)?; + } + if self.is_empty() { + return Ok(None); + } + ranks + .iter() + .map(|&q| self.quantile(q).map(Option::unwrap)) + .collect::, _>>() + .map(Some) + } + + /// Returns an implementation-dependent memory estimate including reserved bucket capacity. + pub fn estimated_size(&self) -> usize { + size_of::() + self.buckets.capacity() * size_of::() + } + + /// Serializes the sketch in canonical bucket order using the Rust-specific version 1 format. + /// + /// The little-endian format starts with `UDD1`, followed by `max_buckets: u32`, initial + /// accuracy as `f64`, level as `u32`, total weight as `u64`, and bucket count as `u32`. + /// Each 32-byte bucket contains `key: i64`, `count: u64`, `min: f64`, and `max: f64`. + /// It has no cross-language compatibility contract. + pub fn serialize(&self) -> Vec { + let mut out = Vec::new(); + self.serialize_into(&mut out); + out + } + + /// Replaces `out` with the version 1 image, reusing its allocated capacity. + /// + /// Produces the same bytes as [`Self::serialize`]. + pub fn serialize_into(&self, out: &mut Vec) { + out.clear(); + out.reserve(32 + 32 * self.buckets.len()); + out.extend_from_slice(b"UDD1"); + out.extend_from_slice(&self.max_buckets.to_le_bytes()); + out.extend_from_slice(&self.initial_accuracy.to_le_bytes()); + out.extend_from_slice(&self.level.to_le_bytes()); + out.extend_from_slice(&self.count.to_le_bytes()); + out.extend_from_slice(&(self.buckets.len() as u32).to_le_bytes()); + for b in &self.buckets { + out.extend_from_slice(&b.key.to_le_bytes()); + out.extend_from_slice(&b.count.to_le_bytes()); + out.extend_from_slice(&b.min.to_le_bytes()); + out.extend_from_slice(&b.max.to_le_bytes()); + } + } + + /// Deserializes a Rust-specific version 1 image. + /// + /// Unmodified images produced by this mapping version preserve estimates, error bounds, + /// and subsequent update/merge behavior. Validation checks the represented histogram, + /// including extrema membership, ordering, and counts; it does not authenticate provenance + /// or prove that the level is the smallest one reachable from an original stream. + /// An externally constructed, over-compacted image may therefore be accepted and retain + /// its coarser resolution when merged. Equivalence to replaying original inputs requires + /// states produced by `new`, updates, compatible merges, and unmodified roundtrips. + /// + /// # Errors + /// + /// Returns `InvalidData` for unsupported versions, truncation, trailing bytes, invalid + /// parameters, invalid bucket ordering or extrema, and inconsistent or overflowing counts. + pub fn deserialize(bytes: &[u8]) -> Result { + let (sketch, payload) = Self::deserialize_header(bytes)?; + sketch.deserialize_body(payload) + } + + fn deserialize_header(bytes: &[u8]) -> Result<(Self, &[u8]), Error> { + let bad = invalid_image; + if bytes.len() < 32 || &bytes[..4] != b"UDD1" { + return Err(bad()); + } + let mut input = SketchSlice::new(&bytes[4..]); + let limit = input.read_u32_le().map_err(|_| bad())?; + let accuracy = input.read_f64_le().map_err(|_| bad())?; + let mut sketch = Self::new(limit, accuracy).map_err(|_| bad())?; + sketch.level = input.read_u32_le().map_err(|_| bad())?; + sketch.count = input.read_u64_le().map_err(|_| bad())?; + let len = input.read_u32_le().map_err(|_| bad())? as usize; + let max_level = sketch.mapping.max_level(); + if sketch.level > max_level + || len > limit as usize + || input.remaining().len() / 32 != len + || input.remaining().len() % 32 != 0 + || (len == 0 && (sketch.count != 0 || sketch.level != 0)) + { + return Err(bad()); + } + Ok((sketch, &bytes[32..])) + } + + fn deserialize_body(mut self, payload: &[u8]) -> Result { + // The header validates the complete payload length before allocating. + self.buckets = Vec::with_capacity(payload.len() / 32); + let mut count = 0u64; + for chunk in payload.chunks_exact(32) { + let bucket = Bucket::read(chunk); + bucket.validate_after(self.buckets.last().copied())?; + self.validate_membership(&bucket, None)?; + count = count.checked_add(bucket.count).ok_or_else(invalid_image)?; + self.buckets.push(bucket); + } + if count != self.count { + return Err(invalid_image()); + } + Ok(self) + } + + // Specialize the interval check at both callers: standalone decoding has no known + // interval, while serialized merging can avoid mapping either or both extrema. + #[inline(always)] + fn validate_membership(&self, bucket: &Bucket, known: Option<&Bucket>) -> Result<(), Error> { + // Callers supply a validated bucket with the same parameters, level, and key. + // Monotonicity proves membership throughout its observed interval, including endpoints. + let contains = |value| known.is_some_and(|b| b.min <= value && value <= b.max); + let min_valid = contains(bucket.min) || self.key(bucket.min) == bucket.key; + let max_valid = + bucket.min == bucket.max || contains(bucket.max) || self.key(bucket.max) == bucket.key; + if min_valid && max_valid { + Ok(()) + } else { + Err(invalid_image()) + } + } + + fn key(&self, value: f64) -> i64 { + if value == 0.0 { + return 0; + } + let index = (self.mapping.index(value.abs()) >> self.level) as i64 + 1; + if value < 0.0 { -index } else { index } + } + + fn enforce_budget(&mut self) { + while self.buckets.len() > self.max_buckets as usize { + self.collapse_to(self.level + 1); + } + } + + fn collapse_to(&mut self, level: u32) { + if self.level == level { + return; + } + let mut len = 0; + for i in 0..self.buckets.len() { + let mut bucket = self.buckets[i]; + bucket.key = parent(bucket.key, level - self.level); + if len > 0 && self.buckets[len - 1].key == bucket.key { + self.buckets[len - 1].merge(bucket); + } else { + self.buckets[len] = bucket; + len += 1; + } + } + self.buckets.truncate(len); + self.level = level; + } +} + +impl Bucket { + fn read(chunk: &[u8]) -> Self { + Self { + key: i64::from_le_bytes(chunk[..8].try_into().unwrap()), + count: u64::from_le_bytes(chunk[8..16].try_into().unwrap()), + min: f64::from_le_bytes(chunk[16..24].try_into().unwrap()), + max: f64::from_le_bytes(chunk[24..].try_into().unwrap()), + } + } + + fn validate_after(self, previous: Option) -> Result<(), Error> { + if self.count == 0 + || !self.min.is_finite() + || !self.max.is_finite() + || self.min > self.max + || (self.count == 1 && self.min != self.max) + || self.min.to_bits() == (-0.0f64).to_bits() + || self.max.to_bits() == (-0.0f64).to_bits() + || previous.is_some_and(|last| last.key >= self.key || last.max >= self.min) + { + Err(invalid_image()) + } else { + Ok(()) + } + } + + fn merge(&mut self, other: Self) { + // The public operations check total weight before mutating any bucket. + self.count += other.count; + self.min = self.min.min(other.min); + self.max = self.max.max(other.max); + } + + fn estimate(&self) -> f64 { + if self.min == self.max { + return self.min; + } + let lo = self.min.abs().min(self.max.abs()); + let hi = self.min.abs().max(self.max.abs()); + let magnitude = (lo / (0.5 + 0.5 * (lo / hi))).clamp(lo, hi); + if self.key < 0 { -magnitude } else { magnitude } + } +} + +fn invalid_image() -> Error { + Error::deserial("invalid UDDSketch image") +} + +fn parent(key: i64, levels: u32) -> i64 { + if key == 0 { + return 0; + } + let magnitude = (((key.unsigned_abs() - 1) >> levels) + 1) as i64; + if key < 0 { -magnitude } else { magnitude } +} + +fn validate_rank(rank: f64) -> Result<(), Error> { + if !(0.0..=1.0).contains(&rank) { + return Err(Error::invalid_argument("rank must be in [0, 1]")); + } + Ok(()) +} + +fn quantile_index(rank: f64, count: u64) -> u64 { + let bits = rank.to_bits(); + let exponent = ((bits >> 52) & 0x7ff) as u32; + if exponent == 0 { + return 0; + } + let significand = (bits & ((1u64 << 52) - 1)) | (1u64 << 52); + let product = u128::from(significand) * u128::from(count - 1); + let shift = 1075 - exponent; + if shift >= 128 { + 0 + } else { + (product >> shift) as u64 + } +} + +fn next_up(value: f64) -> f64 { + if value == f64::INFINITY { + value + } else { + f64::from_bits(value.to_bits() + 1) + } +} + +fn relative_difference(larger: f64, smaller: f64, denominator: f64) -> f64 { + if larger == smaller { + return 0.0; + } + // Sterbenz's lemma makes subtraction exact for nearby operands, including subnormals. + let difference = if larger <= 2.0 * smaller { + larger - smaller + } else { + next_up(larger - smaller) + }; + next_up(difference / denominator) +} diff --git a/tests-integration/Cargo.toml b/tests-integration/Cargo.toml index 79021d6c..e90cb3ba 100644 --- a/tests-integration/Cargo.toml +++ b/tests-integration/Cargo.toml @@ -34,6 +34,7 @@ datasketches = { workspace = true, features = [ "tdigest", "theta", "tuple", + "uddsketch", ] } googletest = { workspace = true } insta = { workspace = true } diff --git a/tests-integration/tests/serde_tests.rs b/tests-integration/tests/serde_tests.rs index bd723168..0df3bb70 100644 --- a/tests-integration/tests/serde_tests.rs +++ b/tests-integration/tests/serde_tests.rs @@ -71,3 +71,6 @@ mod theta; #[path = "serde_tests/tuple.rs"] mod tuple; + +#[path = "serde_tests/uddsketch.rs"] +mod uddsketch; diff --git a/tests-integration/tests/serde_tests/uddsketch.rs b/tests-integration/tests/serde_tests/uddsketch.rs new file mode 100644 index 00000000..e82f79d3 --- /dev/null +++ b/tests-integration/tests/serde_tests/uddsketch.rs @@ -0,0 +1,322 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::error::ErrorKind; +use datasketches::uddsketch::UddSketch; + +fn assert_serialized_merge(target: &UddSketch, bytes: &[u8]) { + let before = target.serialize(); + let mut expected = target.clone(); + let reference = UddSketch::deserialize(bytes).and_then(|other| expected.merge(&other)); + let mut actual = target.clone(); + let result = actual.merge_serialized(bytes); + if result.is_err() { + assert_eq!(actual.serialize(), before); + } + assert_eq!( + result.map_err(|e| e.kind()), + reference.map_err(|e| e.kind()) + ); + assert_eq!(actual.serialize(), expected.serialize()); + assert_eq!(actual.relative_error(), expected.relative_error()); + assert_eq!( + actual.quantiles(&[0.0, 0.01, 0.5, 0.99, 1.0]).unwrap(), + expected.quantiles(&[0.0, 0.01, 0.5, 0.99, 1.0]).unwrap() + ); +} + +#[test] +fn serialized_merge_intervals_and_levels() { + let values = [ + -f64::MAX, + -2.0, + -1.009, + -1.001, + -f64::from_bits(5), + -f64::from_bits(1), + -0.0, + 0.0, + f64::from_bits(1), + f64::from_bits(5), + 1.001, + 1.009, + 2.0, + f64::MAX, + ]; + for limit in [3, 16, 512] { + for accuracy in [1e-9, 0.01, 0.9] { + let mut target = UddSketch::new(limit, accuracy).unwrap(); + for value in values { + target.update(value).unwrap(); + } + assert_serialized_merge(&target, &target.serialize()); + let mut sources = vec![UddSketch::new(limit, accuracy).unwrap()]; + for value in values { + for probe in [value, value.next_up(), value.next_down()] { + if probe.is_finite() { + let mut source = UddSketch::new(limit, accuracy).unwrap(); + source.update_weighted(probe, 2).unwrap(); + sources.push(source); + } + } + } + let mut interior = UddSketch::new(limit, accuracy).unwrap(); + for value in [-1.005, -0.0, 1.005] { + interior.update(value).unwrap(); + } + sources.push(interior); + for source in sources { + assert_serialized_merge(&target, &source.serialize()); + assert_serialized_merge(&source, &target.serialize()); + } + } + } +} + +#[test] +fn serialized_merge_validates_entire_image_atomically() { + let mut target = UddSketch::default(); + for value in [-2.0, -1.009, -1.001, 0.0, 1.001, 1.009, 2.0] { + target.update(value).unwrap(); + } + let bytes = target.serialize(); + let mut saturated = UddSketch::default(); + saturated.update_weighted(1.005, u64::MAX - 1).unwrap(); + let incompatible = UddSketch::new(3, 0.1).unwrap(); + let empty = UddSketch::default(); + for receiver in [&target, &saturated, &incompatible, &empty] { + assert_serialized_merge(receiver, &bytes); + for len in 0..bytes.len() { + assert_serialized_merge(receiver, &bytes[..len]); + } + let mut trailing = bytes.clone(); + trailing.push(0); + assert_serialized_merge(receiver, &trailing); + // Corrupt every field, including buckets following a valid, mergeable prefix. + for offset in 0..bytes.len() { + for bit in 0..8 { + let mut changed = bytes.clone(); + changed[offset] ^= 1 << bit; + assert_serialized_merge(receiver, &changed); + } + } + } + // A false declared total must not allow a bucket addition to overflow before validation. + let mut image = saturated.serialize(); + image[20..28].copy_from_slice(&1u64.to_le_bytes()); + assert_serialized_merge(&saturated, &image); + let mut full = UddSketch::default(); + full.update_weighted(1.005, u64::MAX).unwrap(); + assert_serialized_merge(&UddSketch::default(), &full.serialize()); + assert_serialized_merge(&saturated, &full.serialize()); +} + +#[test] +fn roundtrip_then_continue_and_merge() { + let mut direct = UddSketch::new(16, 0.001).unwrap(); + for i in -1000..1000 { + direct.update(f64::from(i)).unwrap(); + } + let mut restored = UddSketch::deserialize(&direct.serialize()).unwrap(); + let mut other = UddSketch::new(16, 0.001).unwrap(); + for value in [f64::from_bits(1), -f64::MAX, 1e100, 0.0, 1e-100] { + restored.update(value).unwrap(); + other.update(value).unwrap(); + } + direct.merge(&other).unwrap(); + assert_eq!(restored.serialize(), direct.serialize()); + assert_eq!(restored.relative_error(), direct.relative_error()); +} + +#[test] +fn serialize_into_replaces_previous_images() { + let mut sketch = UddSketch::default(); + for i in -1000..1000 { + sketch.update(f64::from(i)).unwrap(); + } + let mut bytes = vec![0xff; 32_768]; + sketch.serialize_into(&mut bytes); + assert_eq!(bytes, sketch.serialize()); + assert_eq!(UddSketch::deserialize(&bytes).unwrap().serialize(), bytes); + // Successively smaller groups must not retain any bytes from the previous image. + for values in [&[-0.0, f64::from_bits(1)][..], &[]] { + sketch.reset(); + for &value in values { + sketch.update(value).unwrap(); + } + sketch.serialize_into(&mut bytes); + assert_eq!(bytes, sketch.serialize()); + assert_eq!(UddSketch::deserialize(&bytes).unwrap().serialize(), bytes); + } +} + +#[test] +fn rust_format_fixture() { + // This is a local format fixture, not a cross-language TCK snapshot. + let mut sketch = UddSketch::new(3, 0.01).unwrap(); + sketch.update_weighted(0.0, 7).unwrap(); + let fixture = [ + 85, 68, 68, 49, 3, 0, 0, 0, 123, 20, 174, 71, 225, 122, 132, 63, 0, 0, 0, 0, 7, 0, 0, 0, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, + ]; + assert_eq!(sketch.serialize(), fixture); + assert_eq!(UddSketch::deserialize(&fixture).unwrap().total_weight(), 7); +} + +#[test] +fn nonzero_mapping_format_fixture() { + // These UDD1 keys pin the mapping across normal, subnormal, and extreme values. + let values = [ + (f64::from_bits(1), 1i64), + (f64::MIN_POSITIVE, 1803), + (0.1, 37106), + (1.0, 37221), + (1.25, 37233), + (1.5, 37242), + (1.75, 37249), + (f64::MAX, 72709), + ]; + for (value, key) in values { + for sign in [-1.0, 1.0] { + let value = sign * value; + let key = if sign < 0.0 { -key } else { key }; + let mut fixture = b"UDD1".to_vec(); + fixture.extend(3u32.to_le_bytes()); + fixture.extend(0.01f64.to_le_bytes()); + fixture.extend(0u32.to_le_bytes()); + fixture.extend(7u64.to_le_bytes()); + fixture.extend(1u32.to_le_bytes()); + fixture.extend(key.to_le_bytes()); + fixture.extend(7u64.to_le_bytes()); + fixture.extend(value.to_le_bytes()); + fixture.extend(value.to_le_bytes()); + let mut sketch = UddSketch::new(3, 0.01).unwrap(); + sketch.update_weighted(value, 7).unwrap(); + assert_eq!(sketch.serialize(), fixture); + let restored = UddSketch::deserialize(&fixture).unwrap(); + assert_eq!(restored.serialize(), fixture); + assert_eq!(restored.quantile(0.5).unwrap(), Some(value)); + } + } +} + +#[test] +fn external_coarse_image_preserves_histogram_not_replay_history() { + let mut one = UddSketch::default(); + one.update(1.0).unwrap(); + let mut bytes = one.serialize(); + let level = 16u32; + let key = i64::from_le_bytes(bytes[32..40].try_into().unwrap()); + bytes[16..20].copy_from_slice(&level.to_le_bytes()); + bytes[32..40].copy_from_slice(&(((key - 1) >> level) + 1).to_le_bytes()); + let mut restored = UddSketch::deserialize(&bytes).unwrap(); + assert_eq!(restored.serialize(), bytes); + let mut other = UddSketch::default(); + for value in [2.0, 4.0] { + other.update(value).unwrap(); + one.update(value).unwrap(); + } + restored.merge(&other).unwrap(); + assert_eq!(restored.compaction_level(), level); + assert_ne!(restored.serialize(), one.serialize()); + let estimate = restored.quantile(0.5).unwrap().unwrap(); + assert!((estimate - 2.0).abs() / 2.0 <= restored.relative_error()); + assert_eq!(one.quantile(0.5).unwrap(), Some(2.0)); + + // Imported coarse states also support the interval fast path and checked weights. + let image = restored.serialize(); + assert_serialized_merge(&restored, &image); + assert_serialized_merge(&other, &image); + let mut interior = image.clone(); + interior[48..56].copy_from_slice(&1.5f64.to_le_bytes()); + interior[56..64].copy_from_slice(&3.0f64.to_le_bytes()); + assert_serialized_merge(&restored, &interior); + // A coarse receiver's interval cannot prove membership in a finer input bucket. + let mut wrong_child = other.serialize(); + let key = i64::from_le_bytes(wrong_child[32..40].try_into().unwrap()); + wrong_child[32..40].copy_from_slice(&(key + 1).to_le_bytes()); + assert_serialized_merge(&restored, &wrong_child); + for value in [ + 0.0, + -1.0, + 1.0, + 1.5, + f64::from_bits(4.0f64.to_bits() - 1), + 4.0, + f64::MAX, + ] { + for weight in [0, 1, u64::MAX - restored.total_weight()] { + let mut updated = UddSketch::deserialize(&image).unwrap(); + updated.update_weighted(value, weight).unwrap(); + let mut singleton = UddSketch::default(); + singleton.update_weighted(value, weight).unwrap(); + let mut merged = restored.clone(); + merged.merge(&singleton).unwrap(); + assert_eq!(updated.serialize(), merged.serialize()); + assert_eq!(updated.relative_error(), merged.relative_error()); + assert_eq!( + updated.quantiles(&[0.0, 0.5, 1.0]).unwrap(), + merged.quantiles(&[0.0, 0.5, 1.0]).unwrap() + ); + if updated.total_weight() == u64::MAX { + let bytes = updated.serialize(); + assert!(updated.update(value).is_err()); + assert_eq!(updated.serialize(), bytes); + } + } + } +} + +#[test] +fn malformed_images() { + let mut sketch = UddSketch::new(3, 0.01).unwrap(); + for value in [-1.0, 0.0, 1.0] { + sketch.update(value).unwrap(); + } + let bytes = sketch.serialize(); + for len in 0..bytes.len() { + assert!(UddSketch::deserialize(&bytes[..len]).is_err()); + } + let mut trailing = bytes.clone(); + trailing.push(0); + assert!(UddSketch::deserialize(&trailing).is_err()); + for (offset, replacement) in [ + (3, vec![2]), + (4, 2u32.to_le_bytes().to_vec()), + (8, f64::NAN.to_le_bytes().to_vec()), + (16, 64u32.to_le_bytes().to_vec()), + (20, u64::MAX.to_le_bytes().to_vec()), + (28, u32::MAX.to_le_bytes().to_vec()), + (32, i64::MIN.to_le_bytes().to_vec()), + (40, 0u64.to_le_bytes().to_vec()), + (48, f64::NEG_INFINITY.to_le_bytes().to_vec()), + (56, 1.0f64.to_le_bytes().to_vec()), + (64, 0i64.wrapping_sub(1).to_le_bytes().to_vec()), + ] { + let mut bad = bytes.clone(); + bad[offset..offset + replacement.len()].copy_from_slice(&replacement); + assert_eq!( + UddSketch::deserialize(&bad).unwrap_err().kind(), + ErrorKind::InvalidData, + "offset {offset}" + ); + } + let mut overflow = bytes; + overflow[40..48].copy_from_slice(&u64::MAX.to_le_bytes()); + assert!(UddSketch::deserialize(&overflow).is_err()); +} diff --git a/tests-integration/tests/uddsketch_test/accuracy.rs b/tests-integration/tests/uddsketch_test/accuracy.rs new file mode 100644 index 00000000..1525f3eb --- /dev/null +++ b/tests-integration/tests/uddsketch_test/accuracy.rs @@ -0,0 +1,69 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use super::build; +use super::shuffle; + +#[test] +fn value_error_certificate_across_distributions() { + let mut values = (0..6000) + .map(|i| (f64::from(i) / 4.0 - 750.0).exp()) + .filter(|x| x.is_finite()) + .collect::>(); + values.extend(values.clone().into_iter().map(|x| -x)); + values.extend([0.0, -0.0, f64::MAX, -f64::MAX]); + for bits in 1..256 { + values.extend([f64::from_bits(bits), -f64::from_bits(bits)]); + } + values.extend(std::iter::repeat_n(1.0, 500)); + shuffle(&mut values, 41); + for limit in [3, 32, 512, 4096] { + let sketch = build(&values, limit, 0.01); + assert!(sketch.num_buckets() <= limit as usize); + let bound = sketch.relative_error(); + assert!(bound.is_finite()); + let mut sorted = values.clone(); + sorted.sort_by(f64::total_cmp); + // Power-of-two ranks have exact products here, keeping the reference definition explicit. + for i in 0..=1024 { + let q = f64::from(i) / 1024.0; + let truth = sorted[(q * (sorted.len() - 1) as f64).floor() as usize]; + let estimate = sketch.quantile(q).unwrap().unwrap(); + assert!(estimate.is_finite()); + if truth == 0.0 { + assert_eq!(estimate, 0.0); + } else { + assert!( + ((estimate - truth) / truth).abs() <= bound, + "q={q}, truth={truth}, estimate={estimate}, bound={bound}" + ); + } + } + } +} + +#[test] +fn global_collapse_can_exceed_initial_target() { + let values = (1..=10_000) + .map(|i| f64::from(i).powi(3)) + .collect::>(); + let sketch = build(&values, 8, 0.001); + assert!(sketch.compaction_level() > 0); + assert!(sketch.relative_error() > sketch.initial_relative_accuracy()); + assert_eq!(sketch.quantile(0.0).unwrap(), Some(1.0)); + assert_eq!(sketch.quantile(1.0).unwrap(), Some(1e12)); +} diff --git a/tests-integration/tests/uddsketch_test/main.rs b/tests-integration/tests/uddsketch_test/main.rs new file mode 100644 index 00000000..4328078c --- /dev/null +++ b/tests-integration/tests/uddsketch_test/main.rs @@ -0,0 +1,39 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +mod accuracy; +mod merge; +mod update; + +use datasketches::uddsketch::UddSketch; + +fn build(values: &[f64], limit: u32, alpha: f64) -> UddSketch { + let mut sketch = UddSketch::new(limit, alpha).unwrap(); + for &value in values { + sketch.update(value).unwrap(); + } + sketch +} + +fn shuffle(values: &mut [T], mut seed: u64) { + for i in (1..values.len()).rev() { + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + values.swap(i, seed as usize % (i + 1)); + } +} diff --git a/tests-integration/tests/uddsketch_test/merge.rs b/tests-integration/tests/uddsketch_test/merge.rs new file mode 100644 index 00000000..27e831ba --- /dev/null +++ b/tests-integration/tests/uddsketch_test/merge.rs @@ -0,0 +1,149 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::uddsketch::UddSketch; +use quickcheck::QuickCheck; + +use super::build; +use super::shuffle; + +fn assert_equivalent(values: &[f64], limit: u32, alpha: f64, seed: u64) { + let expected = build(values, limit, alpha); + if !values.is_empty() { + let mut sorted = values.to_vec(); + sorted.sort_by(f64::total_cmp); + for q in [0.25, 0.5, 0.75] { + let exact = sorted[(q * (sorted.len() - 1) as f64) as usize]; + let estimate = expected.quantile(q).unwrap().unwrap(); + if exact == 0.0 { + assert_eq!(estimate, 0.0); + } else { + assert!(((estimate - exact) / exact).abs() <= expected.relative_error()); + } + } + } + let mut reordered = values.to_vec(); + shuffle(&mut reordered, seed); + assert_eq!( + build(&reordered, limit, alpha).serialize(), + expected.serialize() + ); + for shards in [1, 2, 7, 31] { + let mut partials = (0..shards) + .map(|_| UddSketch::new(limit, alpha).unwrap()) + .collect::>(); + for (i, &value) in reordered.iter().enumerate() { + partials[i % shards].update(value).unwrap(); + } + let mut linear = UddSketch::new(limit, alpha).unwrap(); + let mut serialized = UddSketch::new(limit, alpha).unwrap(); + for partial in partials.iter().rev() { + linear.merge(partial).unwrap(); + serialized.merge_serialized(&partial.serialize()).unwrap(); + assert_eq!(serialized.serialize(), linear.serialize()); + } + assert_eq!(linear.serialize(), expected.serialize()); + while partials.len() > 1 { + let tree_seed = seed.wrapping_add(partials.len() as u64); + shuffle(&mut partials, tree_seed); + let right = partials.pop().unwrap(); + let left = partials.last_mut().unwrap(); + left.merge_serialized(&right.serialize()).unwrap(); + *left = UddSketch::deserialize(&left.serialize()).unwrap(); + } + assert_eq!(partials[0].serialize(), expected.serialize()); + assert_eq!( + partials[0] + .quantiles(&[0.0, 0.5, 0.9, 0.99, 0.999, 1.0]) + .unwrap(), + expected + .quantiles(&[0.0, 0.5, 0.9, 0.99, 0.999, 1.0]) + .unwrap() + ); + } +} + +#[test] +fn floating_boundary_regression() { + // Re-evaluating logarithms after collapse can disagree with the fixed integer hierarchy. + let controls = [ + 1.0100505037878156, + 1.0304555644703977, + 1.0512728486011127, + 1.0725106839263876, + ]; + let target = 1.127_501_361_846_421; + let mut values = controls.to_vec(); + values.push(target); + assert_equivalent(&values, 3, 0.01, 42); + let mut first = build(&[target], 3, 0.01); + let control = build(&controls, 3, 0.01); + first.merge(&control).unwrap(); + assert_eq!(first.serialize(), build(&values, 3, 0.01).serialize()); +} + +#[test] +fn hierarchy_boundaries_and_extreme_values() { + for alpha in [1e-9, 0.001, 0.01, 0.5, 0.999] { + let step = f64::ln_1p(alpha) - f64::ln_1p(-alpha); + let mut values = vec![ + 0.0, + -0.0, + f64::MAX, + -f64::MAX, + f64::MIN_POSITIVE, + -f64::MIN_POSITIVE, + ]; + for exponent in [-1000, -100, -4, -1, 0, 1, 2, 3, 6, 100, 1000] { + let boundary = (step * f64::from(exponent)).exp(); + if boundary.is_normal() { + for bits in [ + boundary.to_bits() - 1, + boundary.to_bits(), + boundary.to_bits() + 1, + ] { + let value = f64::from_bits(bits); + if value.is_finite() { + values.extend([value, -value]); + } + } + } + } + for bits in 1..128 { + values.extend([f64::from_bits(bits), -f64::from_bits(bits)]); + } + for limit in [3, 7, 32, 512] { + assert_equivalent(&values, limit, alpha, 103); + } + } +} + +#[test] +fn arbitrary_multisets_and_merge_trees() { + fn property(bits: Vec, seed: u64) -> bool { + let values = bits + .into_iter() + .map(f64::from_bits) + .filter(|x| x.is_finite()) + .collect::>(); + assert_equivalent(&values, 3 + (seed % 40) as u32, 0.01, seed.max(1)); + true + } + QuickCheck::new() + .tests(200) + .quickcheck(property as fn(Vec, u64) -> bool); +} diff --git a/tests-integration/tests/uddsketch_test/update.rs b/tests-integration/tests/uddsketch_test/update.rs new file mode 100644 index 00000000..21ad8b02 --- /dev/null +++ b/tests-integration/tests/uddsketch_test/update.rs @@ -0,0 +1,234 @@ +// Licensed to the Apache Software Foundation (ASF) under one +// or more contributor license agreements. See the NOTICE file +// distributed with this work for additional information +// regarding copyright ownership. The ASF licenses this file +// to you under the Apache License, Version 2.0 (the +// "License"); you may not use this file except in compliance +// with the License. You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, +// software distributed under the License is distributed on an +// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY +// KIND, either express or implied. See the License for the +// specific language governing permissions and limitations +// under the License. + +use datasketches::error::ErrorKind; +use datasketches::uddsketch::UddSketch; + +use super::build; + +#[test] +fn invalid_updates_and_merges_are_atomic() { + let mut sketch = build(&[-2.0, 0.0, 3.0], 3, 0.01); + let before = sketch.serialize(); + for value in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert_eq!( + sketch.update(value).unwrap_err().kind(), + ErrorKind::InvalidArgument + ); + assert!(sketch.update_weighted(value, 0).is_err()); + } + for other in [ + UddSketch::new(4, 0.01).unwrap(), + UddSketch::new(3, 0.02).unwrap(), + ] { + assert!(sketch.merge(&other).is_err()); + } + assert_eq!(sketch.serialize(), before); + sketch.update_weighted(8.0, 0).unwrap(); + assert_eq!(sketch.serialize(), before); + let mut full = UddSketch::new(3, 0.01).unwrap(); + full.update_weighted(1.0, u64::MAX - 1).unwrap(); + full.update(2.0).unwrap(); + let before = full.serialize(); + assert!(full.update(3.0).is_err()); + assert!(full.merge(&sketch).is_err()); + assert_eq!(full.serialize(), before); + assert_eq!(UddSketch::deserialize(&before).unwrap().serialize(), before); + assert_eq!(full.quantile(1.0).unwrap(), Some(2.0)); + assert_eq!( + full.quantile(f64::from_bits(1.0f64.to_bits() - 1)).unwrap(), + Some(1.0) + ); +} + +#[test] +fn empty_and_discrete_semantics() { + let empty = UddSketch::default(); + assert_eq!(empty.quantile(0.5).unwrap(), None); + assert_eq!(empty.quantiles(&[]).unwrap(), None); + for rank in [f64::NAN, f64::NEG_INFINITY, -0.01, 1.01] { + assert!(empty.quantile(rank).is_err()); + assert!(empty.quantiles(&[0.5, rank]).is_err()); + } + let sketch = build(&[-100.0, 0.0, 10.0, 20.0], 512, 0.01); + assert_eq!(sketch.quantile(0.5).unwrap(), Some(0.0)); + assert_eq!(sketch.quantile(0.0).unwrap(), Some(-100.0)); + assert_eq!(sketch.quantile(1.0).unwrap(), Some(20.0)); + assert_eq!(sketch.quantiles(&[]).unwrap(), Some(vec![])); + assert_eq!( + build(&[-0.0, 0.0], 3, 0.01).serialize(), + build(&[0.0, 0.0], 3, 0.01).serialize() + ); + for alpha in [0.0, 1e-10, 1.0, f64::NAN] { + assert!(UddSketch::new(512, alpha).is_err()); + } + assert!(UddSketch::new(2, 0.01).is_err()); +} + +#[test] +fn reset_starts_an_independent_stream() { + for limit in [3, 16, 512] { + let mut sketch = build( + &(-2000..2000).map(f64::from).collect::>(), + limit, + 0.001, + ); + assert!(sketch.compaction_level() > 0); + for values in [ + &[][..], + &[-f64::MAX, -0.0, f64::from_bits(1), f64::MAX], + &[1.0, 1.001, 1.005, 2.0], + ] { + sketch.reset(); + let fresh = UddSketch::new(limit, 0.001).unwrap(); + assert_eq!(sketch.serialize(), fresh.serialize()); + assert_eq!(sketch.quantile(0.5).unwrap(), None); + for &value in values { + sketch.update(value).unwrap(); + } + let expected = build(values, limit, 0.001); + assert_eq!(sketch.serialize(), expected.serialize()); + assert_eq!(sketch.relative_error(), expected.relative_error()); + } + sketch.reset(); + sketch.update_weighted(1.0, u64::MAX).unwrap(); + sketch.reset(); + sketch.update(2.0).unwrap(); + assert_eq!(sketch.total_weight(), 1); + assert_eq!(sketch.quantile(0.5).unwrap(), Some(2.0)); + } +} + +#[test] +fn exact_weighted_rank_above_float_integer_precision() { + let mut sketch = UddSketch::default(); + sketch.update_weighted(1.0, 1u64 << 63).unwrap(); + sketch.update_weighted(2.0, (1u64 << 63) - 1).unwrap(); + assert_eq!(sketch.quantile(0.5).unwrap(), Some(1.0)); + assert_eq!( + sketch + .quantile(f64::from_bits(0.5f64.to_bits() + 1)) + .unwrap(), + Some(2.0) + ); +} + +#[test] +fn weighted_intervals_and_gaps_match_singleton_merges() { + let values = [ + 1.001, + 1.009, + 1.005, + 1.002, + 1.0095, + 1.0005, + 2.0, + 0.5, + -1.009, + -1.001, + -1.005, + -1.0095, + -1.0005, + -2.0, + 0.0, + -0.0, + f64::from_bits(1), + f64::from_bits(5), + f64::from_bits(3), + f64::MAX, + -f64::MAX, + ]; + for limit in [3, 16, 512] { + for accuracy in [1e-9, 0.01, 0.9] { + let mut updated = UddSketch::new(limit, accuracy).unwrap(); + let mut merged = UddSketch::new(limit, accuracy).unwrap(); + for (i, value) in values.into_iter().enumerate() { + let weight = i as u64 + 1; + updated.update_weighted(value, weight).unwrap(); + let mut singleton = UddSketch::new(limit, accuracy).unwrap(); + singleton.update_weighted(value, weight).unwrap(); + merged.merge(&singleton).unwrap(); + assert_eq!(updated.serialize(), merged.serialize()); + assert_eq!(updated.relative_error(), merged.relative_error()); + } + } + } +} + +#[test] +fn retained_boundaries_and_saturated_weights_match_merges() { + let values = [ + -f64::MAX, + -3.0, + -2.0, + -1.0, + -f64::from_bits(5), + -f64::from_bits(1), + 0.0, + f64::from_bits(1), + f64::from_bits(5), + 1.0, + 2.0, + 3.0, + f64::MAX, + ]; + // Every retained extremum comes from this input. Probe either side of each one, + // along with gaps, interiors, and both zero encodings. + let mut probes = vec![-2.5, -1.5, -0.0, 0.0, 1.5, 2.5]; + for value in values { + probes.push(value); + for bits in [ + value.to_bits().saturating_sub(1), + value.to_bits().saturating_add(1), + ] { + let neighbor = f64::from_bits(bits); + if neighbor.is_finite() { + probes.push(neighbor); + } + } + } + for limit in [3, 16, 512] { + for accuracy in [1e-9, 0.01, 0.9, f64::from_bits(1.0f64.to_bits() - 1)] { + for input in [&[][..], &values[..]] { + let base = build(input, limit, accuracy); + for &value in &probes { + for weight in [0, 1, u64::MAX - base.total_weight()] { + let mut updated = base.clone(); + updated.update_weighted(value, weight).unwrap(); + // A singleton merge bypasses the interval-based insertion path. + let mut singleton = UddSketch::new(limit, accuracy).unwrap(); + singleton.update_weighted(value, weight).unwrap(); + let mut merged = base.clone(); + merged.merge(&singleton).unwrap(); + assert_eq!(updated.serialize(), merged.serialize()); + assert_eq!(updated.relative_error(), merged.relative_error()); + assert_eq!( + updated.quantiles(&[0.0, 0.5, 0.99, 1.0]).unwrap(), + merged.quantiles(&[0.0, 0.5, 0.99, 1.0]).unwrap() + ); + if updated.total_weight() == u64::MAX { + let bytes = updated.serialize(); + assert!(updated.update(value).is_err()); + updated.update_weighted(value, 0).unwrap(); + assert_eq!(updated.serialize(), bytes); + } + } + } + } + } + } +} diff --git a/xtask/src/main.rs b/xtask/src/main.rs index 9b3fb1f2..780d820e 100644 --- a/xtask/src/main.rs +++ b/xtask/src/main.rs @@ -128,6 +128,7 @@ fn datasketches_features() -> Vec { let Metadata { packages, .. } = MetadataCommand::new() .manifest_path(datasketches_manifest) + .no_deps() .exec() .expect("failed to get cargo metadata");