From ed11605d36aed8d09a9fdbdd6cd6fc821cc93c05 Mon Sep 17 00:00:00 2001 From: jaideeppyne Date: Sun, 30 Aug 2026 23:52:18 +0530 Subject: [PATCH 1/6] fix(tdigest): interpolate quantiles toward the nearer centroid and normalize the left tail --- CHANGELOG.md | 2 + datasketches/src/tdigest/sketch.rs | 17 +- tests-integration/tests/tdigest_test/main.rs | 1 + .../tests/tdigest_test/property.rs | 154 ++++++++++++++++++ .../tests/tdigest_test/sketch.rs | 69 ++++++++ 5 files changed, 236 insertions(+), 7 deletions(-) create mode 100644 tests-integration/tests/tdigest_test/property.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 04d6c0c4..4e135c18 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,6 +47,8 @@ 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. +* T-Digest `quantile` now interpolates toward the nearer of the two bracketing centroids, so results rise with the requested rank. Every digest built from `update` or `merge` previously returned some non-monotonic quantiles, and mean rank error over uniform streams improves by roughly 20x. Results now differ from the current C++, Java, and Go ports, which share the defect. +* T-Digest `quantile` no longer returns values above `max_value()`, and `rank`, `cdf`, and `pmf` now scale the left tail by the total weight, so ranks stay within `[0, 1]` and PMF masses stay non-negative. These branches are reachable only for deserialized digests whose first or last centroid carries more than unit weight, including images from the reference implementation. ## v0.4.0 (2026-08-18) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 613c2c7d..66d163f9 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 @@ -1378,7 +1379,7 @@ impl TDigestView<'_> { if last_weight > 1. && (centroids_weight - weight <= last_weight / 2.) { 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 +1404,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 825a6281..25ecc796 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 00000000..7fddee35 --- /dev/null +++ b/tests-integration/tests/tdigest_test/property.rs @@ -0,0 +1,154 @@ +// 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 +} + +#[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(); + } + + let mut tdigest = digest_of(&values); + let min = tdigest.min_value().unwrap(); + let max = tdigest.max_value().unwrap(); + + let mut previous = f64::NEG_INFINITY; + for step in 0..=RANK_STEPS { + let rank = step as f64 / RANK_STEPS as f64; + let quantile = tdigest.quantile(rank).unwrap(); + assert!( + (min..=max).contains(&quantile), + "quantile {quantile} at rank {rank} escapes [{min}, {max}]" + ); + assert!( + quantile >= previous, + "quantile {quantile} at rank {rank} is below {previous} at the preceding rank" + ); + previous = quantile; + } + + 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 merged = digest_of(&left); + merged.merge(&digest_of(&right)); + let max = merged.max_value().unwrap(); + + let mut previous = merged.min_value().unwrap(); + for step in 0..=RANK_STEPS { + let rank = step as f64 / RANK_STEPS as f64; + let quantile = merged.quantile(rank).unwrap(); + assert!( + (previous..=max).contains(&quantile), + "merged quantile {quantile} at rank {rank} escapes [{previous}, {max}]" + ); + previous = quantile; + } + + TestResult::passed() + } + + QuickCheck::new() + .tests(64) + .min_tests_passed(64) + .rng(Gen::new(900)) + .quickcheck(property as fn(Vec, Vec) -> TestResult); +} + +#[test] +fn prop_cdf_is_a_distribution() { + fn property(values: Vec) -> TestResult { + if !(500..1500).contains(&values.len()) { + return TestResult::discard(); + } + + let mut tdigest = digest_of(&values); + let min = tdigest.min_value().unwrap(); + let max = tdigest.max_value().unwrap(); + if min == max { + return TestResult::discard(); + } + + let mut split_points = Vec::with_capacity(RANK_STEPS); + for step in 1..RANK_STEPS { + let point = min + (max - min) * (step as f64 / RANK_STEPS as f64); + if split_points.last().is_none_or(|last| point > *last) { + split_points.push(point); + } + } + + let mut previous = 0.0; + for point in &split_points { + let rank = tdigest.rank(*point).unwrap(); + assert!( + (0.0..=1.0).contains(&rank), + "rank {rank} at value {point} escapes [0, 1]" + ); + assert!( + rank >= previous, + "rank {rank} at value {point} is below {previous} at the preceding value" + ); + previous = rank; + } + + let pmf = tdigest.pmf(&split_points).unwrap(); + for mass in &pmf { + assert!(*mass >= 0.0, "negative mass {mass} in {pmf:?}"); + } + let total: f64 = pmf.iter().sum(); + assert!((total - 1.0).abs() < 1e-9, "masses sum to {total}"); + + TestResult::passed() + } + + QuickCheck::new() + .tests(128) + .min_tests_passed(128) + .rng(Gen::new(1200)) + .quickcheck(property as fn(Vec) -> TestResult); +} diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index c0cced92..01e7f59c 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -332,3 +332,72 @@ 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_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)); +} From 1c1f1ece6fc03ca9a327dfb4d3647b89bb6fdaaf Mon Sep 17 00:00:00 2001 From: tison Date: Mon, 31 Aug 2026 17:13:31 +0800 Subject: [PATCH 2/6] fix(tdigest): handle two-sample terminal centroids --- datasketches/src/tdigest/sketch.rs | 101 +++++++++++++++--- .../tests/tdigest_test/sketch.rs | 10 ++ 2 files changed, 94 insertions(+), 17 deletions(-) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index 66d163f9..e3091d7d 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -15,6 +15,7 @@ // specific language governing permissions and limitations // under the License. +use std::borrow::Cow; use std::cmp::Ordering; use std::convert::identity; use std::num::NonZeroU64; @@ -1225,6 +1226,69 @@ struct TDigestView<'a> { } impl TDigestView<'_> { + fn quantile_centroids(&self) -> Cow<'_, [Centroid]> { + if self.centroids.is_empty() { + return Cow::Borrowed(self.centroids); + } + + if self.centroids.len() == 1 { + return if self.centroids[0].weight.get() == 2 { + Cow::Owned(vec![ + Centroid { + mean: self.min, + weight: DEFAULT_WEIGHT, + }, + Centroid { + mean: self.max, + weight: DEFAULT_WEIGHT, + }, + ]) + } else { + Cow::Borrowed(self.centroids) + }; + } + + let split_first = self.centroids[0].weight.get() == 2; + let split_last = self.centroids[self.centroids.len() - 1].weight.get() == 2; + if !split_first && !split_last { + return Cow::Borrowed(self.centroids); + } + + // A two-sample terminal centroid contains the observed extreme and the value reflected + // across its mean. Treating both as singletons preserves their exact quantile steps. + let mut centroids = Vec::with_capacity( + self.centroids.len() + usize::from(split_first) + usize::from(split_last), + ); + if split_first { + centroids.push(Centroid { + mean: self.min, + weight: DEFAULT_WEIGHT, + }); + centroids.push(Centroid { + mean: self.centroids[0].mean.mul_add(2.0, -self.min), + weight: DEFAULT_WEIGHT, + }); + } else { + centroids.push(self.centroids[0]); + } + centroids.extend_from_slice(&self.centroids[1..self.centroids.len() - 1]); + if split_last { + centroids.push(Centroid { + mean: self.centroids[self.centroids.len() - 1] + .mean + .mul_add(2.0, -self.max), + weight: DEFAULT_WEIGHT, + }); + centroids.push(Centroid { + mean: self.max, + weight: DEFAULT_WEIGHT, + }); + } else { + centroids.push(self.centroids[self.centroids.len() - 1]); + } + Cow::Owned(centroids) + } + fn pmf(&self, split_points: &[f64]) -> Option> { let mut buckets = self.cdf(split_points)?; for i in (1..buckets.len()).rev() { @@ -1353,13 +1417,16 @@ impl TDigestView<'_> { return None; } - if self.centroids.len() == 1 { - return Some(self.centroids[0].mean); + let query_centroids = self.quantile_centroids(); + let centroids = query_centroids.as_ref(); + + if centroids.len() == 1 { + return Some(centroids[0].mean); } // at least 2 centroids let centroids_weight = self.centroids_weight as f64; - let num_centroids = self.centroids.len(); + let num_centroids = centroids.len(); let weight = rank * centroids_weight; if weight < 1. { return Some(self.min); @@ -1367,40 +1434,40 @@ impl TDigestView<'_> { if weight > centroids_weight - 1. { return Some(self.max); } - let first_weight = self.centroids[0].weight(); + let first_weight = centroids[0].weight(); if first_weight > 1. && weight < first_weight / 2. { return Some( self.min + (((weight - 1.) / ((first_weight / 2.) - 1.)) - * (self.centroids[0].mean - self.min)), + * (centroids[0].mean - self.min)), ); } - let last_weight = self.centroids[num_centroids - 1].weight(); + let last_weight = centroids[num_centroids - 1].weight(); if last_weight > 1. && (centroids_weight - weight <= last_weight / 2.) { return Some( self.max - (((centroids_weight - weight - 1.) / ((last_weight / 2.) - 1.)) - * (self.max - self.centroids[num_centroids - 1].mean)), + * (self.max - centroids[num_centroids - 1].mean)), ); } // interpolate between extremes let mut weight_so_far = first_weight / 2.; for i in 0..(num_centroids - 1) { - let dw = (self.centroids[i].weight() + self.centroids[i + 1].weight()) / 2.; + let dw = (centroids[i].weight() + centroids[i + 1].weight()) / 2.; if weight_so_far + dw > weight { // the target weight is between centroids i and i+1 let mut left_weight = 0.; - if self.centroids[i].weight.get() == 1 { + if centroids[i].weight.get() == 1 { if weight - weight_so_far < 0.5 { - return Some(self.centroids[i].mean); + return Some(centroids[i].mean); } left_weight = 0.5; } let mut right_weight = 0.; - if self.centroids[i + 1].weight.get() == 1 { + if centroids[i + 1].weight.get() == 1 { if weight_so_far + dw - weight <= 0.5 { - return Some(self.centroids[i + 1].mean); + return Some(centroids[i + 1].mean); } right_weight = 0.5; } @@ -1409,19 +1476,19 @@ impl TDigestView<'_> { 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, + centroids[i].mean, distance_to_right, - self.centroids[i + 1].mean, + centroids[i + 1].mean, distance_from_left, )); } weight_so_far += dw; } - let w1 = weight - (centroids_weight) - ((self.centroids[num_centroids - 1].weight()) / 2.); - let w2 = (self.centroids[num_centroids - 1].weight() / 2.) - w1; + let w1 = weight - (centroids_weight) - ((centroids[num_centroids - 1].weight()) / 2.); + let w2 = (centroids[num_centroids - 1].weight() / 2.) - w1; Some(weighted_average( - self.centroids[num_centroids - 1].mean, + centroids[num_centroids - 1].mean, w1, self.max, w2, diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 01e7f59c..61bf44d9 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -385,6 +385,16 @@ fn test_quantile_right_tail_stays_within_max() { assert_that!(tdigest.quantile(5.0 / 30.0).unwrap(), near(10.0, 1e-12)); } +#[test] +fn test_two_sample_terminal_centroids_are_singletons() { + let mut tdigest = + deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 2), (50.0, 1), (90.0, 2)]); + + for (rank, expected) in [(0.2, 20.0), (0.3, 20.0), (0.7, 80.0), (0.8, 100.0)] { + assert_that!(tdigest.quantile(rank).unwrap(), near(expected, 1e-12)); + } +} + #[test] fn test_rank_left_tail_is_a_fraction_of_the_total_weight() { let mut tdigest = From b5abf98c20d9e006216642b404ad51faf3a32f7e Mon Sep 17 00:00:00 2001 From: tison Date: Mon, 31 Aug 2026 17:13:48 +0800 Subject: [PATCH 3/6] docs(changelog): summarize tdigest query corrections --- CHANGELOG.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4e135c18..663c0cd5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,8 +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. -* T-Digest `quantile` now interpolates toward the nearer of the two bracketing centroids, so results rise with the requested rank. Every digest built from `update` or `merge` previously returned some non-monotonic quantiles, and mean rank error over uniform streams improves by roughly 20x. Results now differ from the current C++, Java, and Go ports, which share the defect. -* T-Digest `quantile` no longer returns values above `max_value()`, and `rank`, `cdf`, and `pmf` now scale the left tail by the total weight, so ranks stay within `[0, 1]` and PMF masses stay non-negative. These branches are reachable only for deserialized digests whose first or last centroid carries more than unit weight, including images from the reference implementation. +* T-Digest `quantile`, `rank`, `cdf`, and `pmf` queries now preserve monotonicity and their documented ranges, including for deserialized digests with weighted terminal centroids. Quantile results may differ from existing C++, Java, and Go releases that use the previous interpolation behavior. ## v0.4.0 (2026-08-18) From fe2c1aeb4f6c1b01489fd48b597daec40121e5ba Mon Sep 17 00:00:00 2001 From: tison Date: Mon, 31 Aug 2026 17:14:54 +0800 Subject: [PATCH 4/6] test(tdigest): focus quantile property coverage --- .../tests/tdigest_test/property.rs | 105 ++++-------------- 1 file changed, 20 insertions(+), 85 deletions(-) diff --git a/tests-integration/tests/tdigest_test/property.rs b/tests-integration/tests/tdigest_test/property.rs index 7fddee35..28ac0b4a 100644 --- a/tests-integration/tests/tdigest_test/property.rs +++ b/tests-integration/tests/tdigest_test/property.rs @@ -32,6 +32,22 @@ fn digest_of(values: &[u32]) -> TDigestMut { 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 { @@ -39,24 +55,7 @@ fn prop_quantile_is_non_decreasing_and_within_the_observed_range() { return TestResult::discard(); } - let mut tdigest = digest_of(&values); - let min = tdigest.min_value().unwrap(); - let max = tdigest.max_value().unwrap(); - - let mut previous = f64::NEG_INFINITY; - for step in 0..=RANK_STEPS { - let rank = step as f64 / RANK_STEPS as f64; - let quantile = tdigest.quantile(rank).unwrap(); - assert!( - (min..=max).contains(&quantile), - "quantile {quantile} at rank {rank} escapes [{min}, {max}]" - ); - assert!( - quantile >= previous, - "quantile {quantile} at rank {rank} is below {previous} at the preceding rank" - ); - previous = quantile; - } + assert_quantiles_are_monotonic(&mut digest_of(&values)); TestResult::passed() } @@ -75,20 +74,9 @@ fn prop_merged_quantile_is_non_decreasing() { return TestResult::discard(); } - let mut merged = digest_of(&left); - merged.merge(&digest_of(&right)); - let max = merged.max_value().unwrap(); - - let mut previous = merged.min_value().unwrap(); - for step in 0..=RANK_STEPS { - let rank = step as f64 / RANK_STEPS as f64; - let quantile = merged.quantile(rank).unwrap(); - assert!( - (previous..=max).contains(&quantile), - "merged quantile {quantile} at rank {rank} escapes [{previous}, {max}]" - ); - previous = quantile; - } + let mut tdigest = digest_of(&left); + tdigest.merge(&digest_of(&right)); + assert_quantiles_are_monotonic(&mut tdigest); TestResult::passed() } @@ -99,56 +87,3 @@ fn prop_merged_quantile_is_non_decreasing() { .rng(Gen::new(900)) .quickcheck(property as fn(Vec, Vec) -> TestResult); } - -#[test] -fn prop_cdf_is_a_distribution() { - fn property(values: Vec) -> TestResult { - if !(500..1500).contains(&values.len()) { - return TestResult::discard(); - } - - let mut tdigest = digest_of(&values); - let min = tdigest.min_value().unwrap(); - let max = tdigest.max_value().unwrap(); - if min == max { - return TestResult::discard(); - } - - let mut split_points = Vec::with_capacity(RANK_STEPS); - for step in 1..RANK_STEPS { - let point = min + (max - min) * (step as f64 / RANK_STEPS as f64); - if split_points.last().is_none_or(|last| point > *last) { - split_points.push(point); - } - } - - let mut previous = 0.0; - for point in &split_points { - let rank = tdigest.rank(*point).unwrap(); - assert!( - (0.0..=1.0).contains(&rank), - "rank {rank} at value {point} escapes [0, 1]" - ); - assert!( - rank >= previous, - "rank {rank} at value {point} is below {previous} at the preceding value" - ); - previous = rank; - } - - let pmf = tdigest.pmf(&split_points).unwrap(); - for mass in &pmf { - assert!(*mass >= 0.0, "negative mass {mass} in {pmf:?}"); - } - let total: f64 = pmf.iter().sum(); - assert!((total - 1.0).abs() < 1e-9, "masses sum to {total}"); - - TestResult::passed() - } - - QuickCheck::new() - .tests(128) - .min_tests_passed(128) - .rng(Gen::new(1200)) - .quickcheck(property as fn(Vec) -> TestResult); -} From e32eb1bd2fd852605d60f17fa787bf30c5beaafd Mon Sep 17 00:00:00 2001 From: tison Date: Mon, 31 Aug 2026 17:38:55 +0800 Subject: [PATCH 5/6] fix(tdigest): handle two-sample tail locally --- datasketches/src/tdigest/sketch.rs | 104 ++++-------------- .../tests/tdigest_test/sketch.rs | 8 +- 2 files changed, 23 insertions(+), 89 deletions(-) diff --git a/datasketches/src/tdigest/sketch.rs b/datasketches/src/tdigest/sketch.rs index e3091d7d..b70ce518 100644 --- a/datasketches/src/tdigest/sketch.rs +++ b/datasketches/src/tdigest/sketch.rs @@ -15,7 +15,6 @@ // specific language governing permissions and limitations // under the License. -use std::borrow::Cow; use std::cmp::Ordering; use std::convert::identity; use std::num::NonZeroU64; @@ -1226,69 +1225,6 @@ struct TDigestView<'a> { } impl TDigestView<'_> { - fn quantile_centroids(&self) -> Cow<'_, [Centroid]> { - if self.centroids.is_empty() { - return Cow::Borrowed(self.centroids); - } - - if self.centroids.len() == 1 { - return if self.centroids[0].weight.get() == 2 { - Cow::Owned(vec![ - Centroid { - mean: self.min, - weight: DEFAULT_WEIGHT, - }, - Centroid { - mean: self.max, - weight: DEFAULT_WEIGHT, - }, - ]) - } else { - Cow::Borrowed(self.centroids) - }; - } - - let split_first = self.centroids[0].weight.get() == 2; - let split_last = self.centroids[self.centroids.len() - 1].weight.get() == 2; - if !split_first && !split_last { - return Cow::Borrowed(self.centroids); - } - - // A two-sample terminal centroid contains the observed extreme and the value reflected - // across its mean. Treating both as singletons preserves their exact quantile steps. - let mut centroids = Vec::with_capacity( - self.centroids.len() + usize::from(split_first) + usize::from(split_last), - ); - if split_first { - centroids.push(Centroid { - mean: self.min, - weight: DEFAULT_WEIGHT, - }); - centroids.push(Centroid { - mean: self.centroids[0].mean.mul_add(2.0, -self.min), - weight: DEFAULT_WEIGHT, - }); - } else { - centroids.push(self.centroids[0]); - } - centroids.extend_from_slice(&self.centroids[1..self.centroids.len() - 1]); - if split_last { - centroids.push(Centroid { - mean: self.centroids[self.centroids.len() - 1] - .mean - .mul_add(2.0, -self.max), - weight: DEFAULT_WEIGHT, - }); - centroids.push(Centroid { - mean: self.max, - weight: DEFAULT_WEIGHT, - }); - } else { - centroids.push(self.centroids[self.centroids.len() - 1]); - } - Cow::Owned(centroids) - } - fn pmf(&self, split_points: &[f64]) -> Option> { let mut buckets = self.cdf(split_points)?; for i in (1..buckets.len()).rev() { @@ -1417,16 +1353,13 @@ impl TDigestView<'_> { return None; } - let query_centroids = self.quantile_centroids(); - let centroids = query_centroids.as_ref(); - - if centroids.len() == 1 { - return Some(centroids[0].mean); + if self.centroids.len() == 1 { + return Some(self.centroids[0].mean); } // at least 2 centroids let centroids_weight = self.centroids_weight as f64; - let num_centroids = centroids.len(); + let num_centroids = self.centroids.len(); let weight = rank * centroids_weight; if weight < 1. { return Some(self.min); @@ -1434,40 +1367,43 @@ impl TDigestView<'_> { if weight > centroids_weight - 1. { return Some(self.max); } - let first_weight = centroids[0].weight(); + let first_weight = self.centroids[0].weight(); if first_weight > 1. && weight < first_weight / 2. { return Some( self.min + (((weight - 1.) / ((first_weight / 2.) - 1.)) - * (centroids[0].mean - self.min)), + * (self.centroids[0].mean - self.min)), ); } - let last_weight = centroids[num_centroids - 1].weight(); + 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.)) - * (self.max - centroids[num_centroids - 1].mean)), + * (self.max - self.centroids[num_centroids - 1].mean)), ); } // interpolate between extremes let mut weight_so_far = first_weight / 2.; for i in 0..(num_centroids - 1) { - let dw = (centroids[i].weight() + centroids[i + 1].weight()) / 2.; + let dw = (self.centroids[i].weight() + self.centroids[i + 1].weight()) / 2.; if weight_so_far + dw > weight { // the target weight is between centroids i and i+1 let mut left_weight = 0.; - if centroids[i].weight.get() == 1 { + if self.centroids[i].weight.get() == 1 { if weight - weight_so_far < 0.5 { - return Some(centroids[i].mean); + return Some(self.centroids[i].mean); } left_weight = 0.5; } let mut right_weight = 0.; - if centroids[i + 1].weight.get() == 1 { + if self.centroids[i + 1].weight.get() == 1 { if weight_so_far + dw - weight <= 0.5 { - return Some(centroids[i + 1].mean); + return Some(self.centroids[i + 1].mean); } right_weight = 0.5; } @@ -1476,19 +1412,19 @@ impl TDigestView<'_> { 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( - centroids[i].mean, + self.centroids[i].mean, distance_to_right, - centroids[i + 1].mean, + self.centroids[i + 1].mean, distance_from_left, )); } weight_so_far += dw; } - let w1 = weight - (centroids_weight) - ((centroids[num_centroids - 1].weight()) / 2.); - let w2 = (centroids[num_centroids - 1].weight() / 2.) - w1; + let w1 = weight - (centroids_weight) - ((self.centroids[num_centroids - 1].weight()) / 2.); + let w2 = (self.centroids[num_centroids - 1].weight() / 2.) - w1; Some(weighted_average( - centroids[num_centroids - 1].mean, + self.centroids[num_centroids - 1].mean, w1, self.max, w2, diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index 61bf44d9..302cc4c7 100644 --- a/tests-integration/tests/tdigest_test/sketch.rs +++ b/tests-integration/tests/tdigest_test/sketch.rs @@ -386,13 +386,11 @@ fn test_quantile_right_tail_stays_within_max() { } #[test] -fn test_two_sample_terminal_centroids_are_singletons() { +fn test_quantile_handles_two_sample_last_centroid() { let mut tdigest = - deserialize_with_centroids(100, 0.0, 100.0, &[(10.0, 2), (50.0, 1), (90.0, 2)]); + deserialize_with_centroids(100, 0.0, 100.0, &[(0.0, 1), (50.0, 1), (90.0, 2)]); - for (rank, expected) in [(0.2, 20.0), (0.3, 20.0), (0.7, 80.0), (0.8, 100.0)] { - assert_that!(tdigest.quantile(rank).unwrap(), near(expected, 1e-12)); - } + assert_eq!(tdigest.quantile(0.75), Some(100.0)); } #[test] From c652bbce97ee1c031b7fad05ae3c35f2c205ed79 Mon Sep 17 00:00:00 2001 From: tison Date: Mon, 31 Aug 2026 17:39:00 +0800 Subject: [PATCH 6/6] docs(changelog): describe tdigest bug fix --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 663c0cd5..3f1f1bcf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -47,7 +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. -* T-Digest `quantile`, `rank`, `cdf`, and `pmf` queries now preserve monotonicity and their documented ranges, including for deserialized digests with weighted terminal centroids. Quantile results may differ from existing C++, Java, and Go releases that use the previous interpolation behavior. +* 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)