diff --git a/CHANGELOG.md b/CHANGELOG.md index 04d6c0c..3f1f1bc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,6 +47,7 @@ All significant changes to this project will be documented in this file. * HLL, Theta, and Tuple deserializers now return `InvalidData` for malformed payload sizes and entry counts instead of risking oversized allocations or decoding failures. * Malformed CPC images now return `InvalidData` instead of panicking. * Seeded deserializers now return `InvalidData` rather than panicking when the caller supplies a seed whose hash is the reserved zero value. +* Fix T-Digest interpolation and tail calculations that could produce non-monotonic or out-of-range quantiles and invalid rank, CDF, or PMF values. ## v0.4.0 (2026-08-18) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 613c2c7..b70ce51 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -1278,8 +1278,9 @@ impl TDigestView<'_> { return Some(if value == self.min { 0.5 / centroids_weight } else { - 1. + (((value - self.min) / (first_mean - self.min)) - * ((self.centroids[0].weight() / 2.) - 1.)) + (1. + (((value - self.min) / (first_mean - self.min)) + * ((self.centroids[0].weight() / 2.) - 1.))) + / centroids_weight }); } return Some(0.); // should never happen @@ -1376,9 +1377,12 @@ impl TDigestView<'_> { } let last_weight = self.centroids[num_centroids - 1].weight(); if last_weight > 1. && (centroids_weight - weight <= last_weight / 2.) { + if last_weight == 2. { + return Some(self.max); + } return Some( self.max - + (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.)) + - (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.)) * (self.max - self.centroids[num_centroids - 1].mean)), ); } @@ -1403,13 +1407,15 @@ impl TDigestView<'_> { } right_weight = 0.5; } - let w1 = weight - weight_so_far - left_weight; - let w2 = weight_so_far + dw - weight - right_weight; + // Each centroid is weighted by the distance from the target to the *other* + // centroid, so the estimate approaches the nearer one. + let distance_from_left = weight - weight_so_far - left_weight; + let distance_to_right = weight_so_far + dw - weight - right_weight; return Some(weighted_average( self.centroids[i].mean, - w1, + distance_to_right, self.centroids[i + 1].mean, - w2, + distance_from_left, )); } weight_so_far += dw; diff --git a/tests-integration/tests/tdigest_test/main.rs b/tests-integration/tests/tdigest_test/main.rs index 825a628..25ecc79 100644 --- a/tests-integration/tests/tdigest_test/main.rs +++ b/tests-integration/tests/tdigest_test/main.rs @@ -15,4 +15,5 @@ // specific language governing permissions and limitations // under the License. +mod property; mod sketch; diff --git a/tests-integration/tests/tdigest_test/property.rs b/tests-integration/tests/tdigest_test/property.rs new file mode 100644 index 0000000..28ac0b4 --- /dev/null +++ b/tests-integration/tests/tdigest_test/property.rs @@ -0,0 +1,89 @@ +// 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. + +//! Property-based t-digest tests. + +use datasketches::tdigest::TDigestMut; +use quickcheck::Gen; +use quickcheck::QuickCheck; +use quickcheck::TestResult; + +const RANK_STEPS: usize = 500; + +fn digest_of(values: &[u32]) -> TDigestMut { + let mut tdigest = TDigestMut::new(100).unwrap(); + for value in values { + tdigest.update(f64::from(*value) / 4096.0); + } + tdigest +} + +fn assert_quantiles_are_monotonic(tdigest: &mut TDigestMut) { + let min = tdigest.min_value().unwrap(); + let max = tdigest.max_value().unwrap(); + let mut previous = min; + + for step in 0..=RANK_STEPS { + let rank = step as f64 / RANK_STEPS as f64; + let quantile = tdigest.quantile(rank).unwrap(); + assert!( + (previous..=max).contains(&quantile), + "quantile {quantile} at rank {rank} is outside [{previous}, {max}]" + ); + previous = quantile; + } +} + +#[test] +fn prop_quantile_is_non_decreasing_and_within_the_observed_range() { + fn property(values: Vec) -> TestResult { + if !(500..1500).contains(&values.len()) { + return TestResult::discard(); + } + + assert_quantiles_are_monotonic(&mut digest_of(&values)); + + TestResult::passed() + } + + QuickCheck::new() + .tests(128) + .min_tests_passed(128) + .rng(Gen::new(1200)) + .quickcheck(property as fn(Vec) -> TestResult); +} + +#[test] +fn prop_merged_quantile_is_non_decreasing() { + fn property(left: Vec, right: Vec) -> TestResult { + if left.len() < 300 || right.len() < 300 { + return TestResult::discard(); + } + + let mut tdigest = digest_of(&left); + tdigest.merge(&digest_of(&right)); + assert_quantiles_are_monotonic(&mut tdigest); + + TestResult::passed() + } + + QuickCheck::new() + .tests(64) + .min_tests_passed(64) + .rng(Gen::new(900)) + .quickcheck(property as fn(Vec, Vec) -> TestResult); +} diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index c0cced9..302cc4c 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -332,3 +332,80 @@ fn test_estimate_repeat_values() { } assert_eq!(tdigest.quantile(0.9), Some(1.0)); } + +/// Builds a digest whose centroids carry the given weights. +/// +/// Compression never merges the extreme centroids, so digests built through `update` and `merge` +/// always keep unit-weight tails. Heavier tails arrive only through deserialization, including the +/// reference implementation format, and they select the tail interpolation branches. +fn deserialize_with_centroids(k: u16, min: f64, max: f64, centroids: &[(f64, u64)]) -> TDigestMut { + const PREAMBLE_LONGS: u8 = 2; + const SERIAL_VERSION: u8 = 1; + const FAMILY_TDIGEST: u8 = 20; + + let mut bytes = vec![PREAMBLE_LONGS, SERIAL_VERSION, FAMILY_TDIGEST]; + bytes.extend_from_slice(&k.to_le_bytes()); + bytes.push(0); // flags + bytes.extend_from_slice(&0u16.to_le_bytes()); // unused + bytes.extend_from_slice(&(centroids.len() as u32).to_le_bytes()); + bytes.extend_from_slice(&0u32.to_le_bytes()); // buffered values + bytes.extend_from_slice(&min.to_le_bytes()); + bytes.extend_from_slice(&max.to_le_bytes()); + for (mean, weight) in centroids { + bytes.extend_from_slice(&mean.to_le_bytes()); + bytes.extend_from_slice(&weight.to_le_bytes()); + } + TDigestMut::deserialize(&bytes).unwrap() +} + +#[test] +fn test_quantile_moves_toward_the_nearer_bracketing_centroid() { + let mut tdigest = + deserialize_with_centroids(100, -1.0, 21.0, &[(0.0, 4), (10.0, 4), (20.0, 4)]); + + assert_eq!(tdigest.total_weight(), 12); + // Ranks 2/12 and 6/12 sit exactly on the two centroids bracketing the first interval. + assert_that!(tdigest.quantile(2.0 / 12.0).unwrap(), near(0.0, 1e-12)); + assert_that!(tdigest.quantile(3.0 / 12.0).unwrap(), near(2.5, 1e-12)); + assert_that!(tdigest.quantile(4.0 / 12.0).unwrap(), near(5.0, 1e-12)); + assert_that!(tdigest.quantile(5.0 / 12.0).unwrap(), near(7.5, 1e-12)); + assert_that!(tdigest.quantile(6.0 / 12.0).unwrap(), near(10.0, 1e-12)); +} + +#[test] +fn test_quantile_right_tail_stays_within_max() { + let mut tdigest = + deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 10), (50.0, 10), (90.0, 10)]); + + assert_eq!(tdigest.max_value(), Some(100.0)); + assert_that!(tdigest.quantile(0.9).unwrap(), near(95.0, 1e-12)); + assert_that!(tdigest.quantile(29.0 / 30.0).unwrap(), near(100.0, 1e-12)); + // Mirrors the left tail, which interpolates from min up to the first centroid mean. + assert_that!(tdigest.quantile(1.0 / 30.0).unwrap(), near(0.0, 1e-12)); + assert_that!(tdigest.quantile(5.0 / 30.0).unwrap(), near(10.0, 1e-12)); +} + +#[test] +fn test_quantile_handles_two_sample_last_centroid() { + let mut tdigest = + deserialize_with_centroids(100, 0.0, 100.0, &[(0.0, 1), (50.0, 1), (90.0, 2)]); + + assert_eq!(tdigest.quantile(0.75), Some(100.0)); +} + +#[test] +fn test_rank_left_tail_is_a_fraction_of_the_total_weight() { + let mut tdigest = + deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 10), (50.0, 10), (90.0, 10)]); + + assert_that!(tdigest.rank(5.0).unwrap(), near(0.1, 1e-12)); + assert_that!(tdigest.rank(10.0).unwrap(), near(5.0 / 30.0, 1e-12)); + // The right tail is the mirror image and pins the scale the left tail must match. + assert_that!(tdigest.rank(95.0).unwrap(), near(0.9, 1e-12)); + assert_that!(tdigest.rank(90.0).unwrap(), near(25.0 / 30.0, 1e-12)); + + let pmf = tdigest.pmf(&[5.0, 95.0]).unwrap(); + assert_that!(pmf[0], near(0.1, 1e-12)); + assert_that!(pmf[1], near(0.8, 1e-12)); + assert_that!(pmf[2], near(0.1, 1e-12)); +}