This is an automated email from the ASF dual-hosted git repository.
tisonkun pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/datasketches-rust.git
The following commit(s) were added to refs/heads/main by this push:
new e7f1b9f fix(tdigest): stabilize numerical interpolation (#287)
e7f1b9f is described below
commit e7f1b9fe559335dc855c587f6dd233349de85cc9
Author: tison <[email protected]>
AuthorDate: Sun Oct 4 23:22:59 2026 +0800
fix(tdigest): stabilize numerical interpolation (#287)
---
CHANGELOG.md | 6 +-
datasketches/src/tdigest/sketch.rs | 341 +++++++++++-------
tests-integration/tests/serde_tests/tdigest.rs | 84 ++++-
tests-integration/tests/tdigest_test/property.rs | 44 +++
tests-integration/tests/tdigest_test/sketch.rs | 420 ++++++++++++++++++++++-
5 files changed, 751 insertions(+), 144 deletions(-)
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c5c6553..995aa64 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -18,6 +18,7 @@ All significant changes to this project will be documented in
this file.
### Improvements
+* Serializing buffered `TDigestMut` inputs is faster when their total weight
is at most `k / 2`.
* Improve the readability of `Debug` output for HLL and CPC sketches and
unions.
* `KllSketch::quantile` and `ReqSketch::quantile` are faster at ranks `0.0`
and `1.0`, with query work independent of the number of retained samples.
* `BloomFilter::insert` is faster for integer and raw-byte inputs.
`BloomFilter::contains_and_insert` is also faster when checking already-present
integer values.
@@ -45,8 +46,9 @@ All significant changes to this project will be documented in
this file.
* `ReqSketch::cdf` and its sorted view now preserve the rank at each split
point and end at exactly `1.0`, avoiding extra rounding from accumulating PMF
values.
* `TDigestMut` updates and merges now panic without changing the digest if the
total weight would exceed `u64::MAX`, including in release builds.
* `TDigestMut::merge` now preserves the true minimum and maximum from both
inputs, including compressed digests.
-* `TDigest` and `TDigestMut` quantile queries no longer extrapolate above the
stored maximum when large total weights cause rounding near rank `1.0`.
-* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for
invalid flags or extrema, out-of-range values, unsorted centroids, or non-empty
images with no stored values.
+* Fix non-finite, out-of-range, or decreasing `TDigest` and `TDigestMut` query
results for extreme finite values, large weights, and updated or merged
deserialized digests. Single-centroid queries now honor stored extrema, and PMF
buckets remain nonnegative.
+* T-Digest compression reduces rounding loss in centroid means when merging
very unequal weights, including opposite-sign centroids.
+* `TDigest` and `TDigestMut` deserialization now returns `InvalidData` for
malformed images, including unsorted centroids, values outside stored extrema,
and single-sample images with distinct extrema.
## v0.5.0 (2026-09-04)
diff --git a/datasketches/src/tdigest/sketch.rs
b/datasketches/src/tdigest/sketch.rs
index 5a7d1c6..ef39ffd 100644
--- a/datasketches/src/tdigest/sketch.rs
+++ b/datasketches/src/tdigest/sketch.rs
@@ -490,8 +490,7 @@ impl TDigestMut {
if value > self.max {
return Some(1.0);
}
- // one centroid and value == min == max
- if self.buffer.len() == 1 {
+ if self.min == self.max {
return Some(0.5);
}
@@ -765,7 +764,8 @@ impl TDigestMut {
"malformed data: centroid means must be finite, within
extrema, and nondecreasing",
));
}
- checked_weight_sum(compressed_weight, num_buffered as u64)?;
+ let total_weight = checked_weight_sum(compressed_weight, num_buffered
as u64)?;
+ check_single_sample_extrema(min, max, total_weight)?;
let mut buffered_values_valid = true;
for bytes in buffered_payload.chunks_exact(buffered_value_bytes) {
let value = if is_f32 {
@@ -846,6 +846,7 @@ impl TDigestMut {
"malformed data: centroid means in compat double
format must be finite, within extrema, and nondecreasing",
));
}
+ check_single_sample_extrema(min, max, total_weight)?;
Ok(TDigestMut::make(
k,
false,
@@ -900,6 +901,7 @@ impl TDigestMut {
"malformed data: centroid means in compat float format
must be finite, within extrema, and nondecreasing",
));
}
+ check_single_sample_extrema(min, max, total_weight)?;
Ok(TDigestMut::make(
k,
false,
@@ -939,30 +941,38 @@ impl TDigestMut {
let mut num_centroids = 1;
let len = centroids.len();
- let compressed_weight = self.compressed_weight as f64;
- let normalizer = scale_function::normalizer(2.0 * f64::from(self.k),
compressed_weight);
- let mut current = 1;
- let mut weight_so_far = 0.;
- while current < len {
- let c = centroids[current];
- let proposed_weight = centroids[num_centroids - 1].weight() +
c.weight();
- if should_merge_centroid(
- current,
- len,
- weight_so_far,
- proposed_weight,
- compressed_weight,
- normalizer,
- ) {
- // merge into existing centroid
- centroids[num_centroids - 1].add(c);
- } else {
- // copy to a new centroid
- weight_so_far += centroids[num_centroids - 1].weight();
- centroids[num_centroids] = c;
- num_centroids += 1;
+ // No two positive integer weights fit when n <= k / 2. For a positive
K_2
+ // normalizer, the maximum limit is n / (4 * normalizer) = r * (ln(r)
+ 6),
+ // where r = n / (2 * k) <= 1/4: the limit stays below 1.154, hence
below 2.
+ // A nonpositive normalizer also rejects all merges at these small
counts.
+ if total_weight <= u64::from(self.k) / 2 {
+ num_centroids = len;
+ } else {
+ let compressed_weight = self.compressed_weight as f64;
+ let normalizer = scale_function::normalizer(2.0 *
f64::from(self.k), compressed_weight);
+ let mut current = 1;
+ let mut weight_so_far = 0.;
+ while current < len {
+ let c = centroids[current];
+ let proposed_weight = centroids[num_centroids - 1].weight() +
c.weight();
+ if should_merge_centroid(
+ current,
+ len,
+ weight_so_far,
+ proposed_weight,
+ compressed_weight,
+ normalizer,
+ ) {
+ // merge into existing centroid
+ centroids[num_centroids - 1].add(c);
+ } else {
+ // copy to a new centroid
+ weight_so_far += centroids[num_centroids - 1].weight();
+ centroids[num_centroids] = c;
+ num_centroids += 1;
+ }
+ current += 1;
}
- current += 1;
}
centroids.truncate(num_centroids);
@@ -1510,6 +1520,8 @@ impl TDigest {
/// Computes the approximate quantile for the given normalized rank.
///
+ /// Ranks `0.0` and `1.0` return the stored minimum and maximum,
respectively.
+ ///
/// Returns `None` if this t-digest is empty.
///
/// # Panics
@@ -1640,43 +1652,44 @@ impl TDigestView<'_> {
if value > self.max {
return Some(1.0);
}
- // one centroid and value == min == max
- if self.centroids.len() == 1 {
+ if self.min == self.max {
return Some(0.5);
}
let centroids_weight = self.centroids_weight as f64;
let num_centroids = self.centroids.len();
- // left tail
- let first_mean = self.centroids[0].mean;
- if value < first_mean {
- if first_mean - self.min > 0. {
- 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.)))
- / centroids_weight
- });
- }
- return Some(0.); // should never happen
+ // Reserve at most one sample for each stored extremum, without
passing the adjacent
+ // centroid center. Updates can place a unit-weight centroid inside
the stored extrema,
+ // leaving no mass to interpolate in that tail.
+ let first = &self.centroids[0];
+ if value < first.mean {
+ let center = centroid_center(0, first.weight.get());
+ return Some(if value == self.min {
+ 0.5 / centroids_weight
+ } else {
+ interpolate(
+ center.min(1.),
+ center,
+ interpolation_fraction(value, self.min, first.mean),
+ ) / centroids_weight
+ });
}
- // right tail
- let last_mean = self.centroids[num_centroids - 1].mean;
- if value > last_mean {
- if self.max - last_mean > 0. {
- return Some(if value == self.max {
- 1. - (0.5 / centroids_weight)
- } else {
- 1.0 - ((1.0
- + (((self.max - value) / (self.max - last_mean))
- * ((self.centroids[num_centroids - 1].weight() /
2.) - 1.)))
- / centroids_weight)
- });
- }
- return Some(1.); // should never happen
+ let last = &self.centroids[num_centroids - 1];
+ if value > last.mean {
+ let center =
+ centroid_center(self.centroids_weight - last.weight.get(),
last.weight.get());
+ return Some(if value == self.max {
+ // Round the maximum's mass in the same coordinates as the
centroid centers.
+ centroid_center(self.centroids_weight - 1, 1) /
centroids_weight
+ } else {
+ interpolate(
+ center,
+ center.max((self.centroids_weight - 1) as f64),
+ interpolation_fraction(value, last.mean, self.max),
+ ) / centroids_weight
+ });
}
let mut lower = self
@@ -1696,30 +1709,24 @@ impl TDigestView<'_> {
upper -= 1;
}
- let mut weight_below = 0.;
- let mut i = 0;
- while i < lower {
- weight_below += self.centroids[i].weight();
- i += 1;
- }
- weight_below += self.centroids[lower].weight() / 2.;
-
- let mut weight_delta = 0.;
- while i < upper {
- weight_delta += self.centroids[i].weight();
- i += 1;
- }
- weight_delta -= self.centroids[lower].weight() / 2.;
- weight_delta += self.centroids[upper].weight() / 2.;
+ let weight_below = self.centroids[..lower].iter().map(|c|
c.weight.get()).sum();
+ let weight_between: u64 = self.centroids[lower..upper]
+ .iter()
+ .map(|c| c.weight.get())
+ .sum();
+ let left = &self.centroids[lower];
+ let right = &self.centroids[upper];
+ let fraction = if left.mean < right.mean {
+ interpolation_fraction(value, left.mean, right.mean)
+ } else {
+ 0.5
+ };
Some(
- if self.centroids[upper].mean - self.centroids[lower].mean > 0. {
- (weight_below
- + (weight_delta * (value - self.centroids[lower].mean)
- / (self.centroids[upper].mean -
self.centroids[lower].mean)))
- / centroids_weight
- } else {
- (weight_below + weight_delta / 2.) / centroids_weight
- },
+ interpolate(
+ centroid_center(weight_below, left.weight.get()),
+ centroid_center(weight_below + weight_between,
right.weight.get()),
+ fraction,
+ ) / centroids_weight,
)
}
@@ -1771,8 +1778,8 @@ impl TDigestView<'_> {
ranks: impl DoubleEndedIterator<Item = (usize, f64)>,
quantiles: &mut [f64],
) {
- if self.centroids.len() == 1 {
- quantiles.fill(self.centroids[0].mean);
+ if self.min == self.max {
+ quantiles.fill(self.min);
return;
}
@@ -1783,7 +1790,7 @@ impl TDigestView<'_> {
.peekable();
// Consume the right tail from the back, so the centroid scan can omit
all tail checks.
- while let Some((index, _)) = queries.next_if(|&(_, weight)| weight >
centroids_weight - 1.)
+ while let Some((index, _)) = queries.next_if(|&(_, weight)|
centroids_weight - weight < 1.)
{
quantiles[index] = self.max;
}
@@ -1792,17 +1799,22 @@ impl TDigestView<'_> {
}
let last = self.centroids.last().unwrap();
- let last_weight = last.weight();
- if last_weight > 1. {
- while let Some((index, weight)) =
- queries.next_if(|&(_, weight)| centroids_weight - weight <=
last_weight / 2.)
+ let last_weight = last.weight.get();
+ let last_center = centroid_center(self.centroids_weight - last_weight,
last_weight);
+ if last_weight > 1 {
+ let tail_end = (self.centroids_weight - 1) as f64;
+ while let Some((index, weight)) = queries.next_if(|&(_, weight)|
weight >= last_center)
{
- quantiles[index] = if last_weight == 2. {
+ // Integer counts can round a wider tail to zero width as
well. At this
+ // boundary prefer the maximum, as for an exact two-sample
tail.
+ quantiles[index] = if last_center == tail_end {
self.max
} else {
- self.max
- - (((centroids_weight - weight - 1.) / ((last_weight /
2.) - 1.))
- * (self.max - last.mean))
+ interpolate(
+ last.mean,
+ self.max,
+ (weight - last_center) / (tail_end - last_center),
+ )
};
}
}
@@ -1816,13 +1828,14 @@ impl TDigestView<'_> {
}
let first = &self.centroids[0];
- let first_weight = first.weight();
- if first_weight > 1. {
- while let Some((index, weight)) =
- queries.next_if(|&(_, weight)| weight < first_weight / 2.)
+ let first_weight = first.weight.get();
+ let first_center = centroid_center(0, first_weight);
+ if first_weight > 1 {
+ // At first_weight == 2, the weight < 1 pre-pass already consumed
this interval.
+ while let Some((index, weight)) = queries.next_if(|&(_, weight)|
weight < first_center)
{
- quantiles[index] = self.min
- + (((weight - 1.) / ((first_weight / 2.) - 1.)) *
(first.mean - self.min));
+ quantiles[index] =
+ interpolate(self.min, first.mean, (weight - 1.) /
(first_center - 1.));
}
}
if queries.peek().is_none() {
@@ -1830,43 +1843,42 @@ impl TDigestView<'_> {
}
// Answer the remaining queries by centroid interval; both streams
advance only forward.
- let mut weight_so_far = first_weight / 2.;
+ let mut weight_before = 0;
+ let mut left_center = first_center;
for pair in self.centroids.windows(2) {
let left = &pair[0];
let right = &pair[1];
- let dw = (left.weight() + right.weight()) / 2.;
+ let right_before = weight_before + left.weight.get();
+ let right_center = centroid_center(right_before,
right.weight.get());
loop {
let Some(&(index, weight)) = queries.peek() else {
return;
};
- if weight_so_far + dw <= weight {
+ if right_center <= weight {
break;
}
queries.next();
- let mut left_weight = 0.;
+ let mut start = left_center;
if left.weight.get() == 1 {
- if weight - weight_so_far < 0.5 {
+ start = right_before as f64;
+ if weight < start {
quantiles[index] = left.mean;
continue;
}
- left_weight = 0.5;
}
- let mut right_weight = 0.;
+ let mut end = right_center;
if right.weight.get() == 1 {
- if weight_so_far + dw - weight <= 0.5 {
+ end = right_before as f64;
+ if weight >= end {
quantiles[index] = right.mean;
continue;
}
- right_weight = 0.5;
}
- // 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;
quantiles[index] =
- weighted_average(left.mean, distance_to_right, right.mean,
distance_from_left);
+ interpolate(left.mean, right.mean, (weight - start) / (end
- start));
}
- weight_so_far += dw;
+ weight_before = right_before;
+ left_center = right_center;
}
// Rounding at large total weights can exhaust the scan near the
maximum.
@@ -2006,22 +2018,27 @@ impl Centroid {
.checked_add(other.weight.get())
.expect("weight overflow");
- let (self_mean, other_mean) = (self.mean, other.mean);
- let ratio_other = other_weight / total_weight;
- let delta = other_mean - self_mean;
- self.mean = if delta.is_finite() {
- delta.mul_add(ratio_other, self_mean)
+ // Start at the heavier centroid so a small contribution is not
rounded away in 1 - ratio.
+ let (start, end, fraction) = if self_weight >= other_weight {
+ (self.mean, other.mean, other_weight / total_weight)
} else {
- let ratio_self = self_weight / total_weight;
- self_mean.mul_add(ratio_self, other_mean * ratio_other)
+ (other.mean, self.mean, self_weight / total_weight)
+ };
+ debug_assert!(start.is_finite() && end.is_finite() && fraction > 0. &&
fraction <= 0.5);
+
+ // Same-sign means within a factor of two in magnitude subtract
exactly. For wider
+ // gaps, a step of at most 1/2 cannot cross either endpoint even if
subtraction rounds.
+ // Queries can interpolate farther and still need the general bounds.
+ self.mean = if start.is_sign_positive() == end.is_sign_positive() {
+ (end - start).mul_add(fraction, start)
+ } else if 1. - fraction == 1. {
+ // Both corrections can matter even when 1 - fraction rounds to 1.
+ // Halving before subtraction avoids overflow; the FMA combines the
+ // correction with start before rounding the result.
+ (end * 0.5 - start * 0.5).mul_add(2. * fraction, start)
+ } else {
+ interpolate(start, end, fraction)
};
-
- debug_assert!(
- self.mean.is_finite(),
- "Centroid's mean must be finite; self: {}, other: {}",
- self_mean,
- other_mean
- );
}
fn weight(&self) -> f64 {
@@ -2064,6 +2081,16 @@ fn check_extrema(min: f64, max: f64, format: &'static
str) -> Result<(), Error>
Ok(())
}
+fn check_single_sample_extrema(min: f64, max: f64, total_weight: u64) ->
Result<(), Error> {
+ // Tail interpolation assumes that distinct extrema represent at least two
samples.
+ if total_weight == 1 && min != max {
+ return Err(Error::deserial(
+ "malformed data: a single sample must have equal extrema",
+ ));
+ }
+ Ok(())
+}
+
fn check_nonzero(value: u64, tag: &'static str) -> Result<NonZeroU64, Error> {
NonZeroU64::new(value)
.ok_or_else(|| Error::deserial(format!("malformed data: {tag} cannot
be zero")))
@@ -2111,14 +2138,64 @@ mod scale_function {
}
}
-fn weighted_average(x1: f64, w1: f64, x2: f64, w2: f64) -> f64 {
- let total_weight = w1 + w2;
- let ratio = w2 / total_weight;
- if x1.is_sign_positive() != x2.is_sign_positive() {
- // Subtracting opposite-signed finite extremes can overflow.
- x1 * (1. - ratio) + x2 * ratio
+/// Keeps cumulative counts exact until a centroid's midpoint is needed for
interpolation.
+/// All intervals and tails use this same center to avoid inconsistent
boundary rounding.
+fn centroid_center(weight_before: u64, weight: u64) -> f64 {
+ (weight_before + weight / 2) as f64 + (weight % 2) as f64 * 0.5
+}
+
+/// Interpolates finite endpoints at a fraction in `[0, 1]`, preserving their
bounds.
+///
+/// The sign split follows WG21 P0811R3:
+/// <https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2019/p0811r3.html>.
+fn interpolate(start: f64, end: f64, fraction: f64) -> f64 {
+ debug_assert!(start.is_finite() && end.is_finite() &&
(0.0..=1.0).contains(&fraction));
+ // Callers guarantee the fraction's range, but division may round it to an
endpoint.
+ // Even an FMA need not reproduce `end` at 1: its subtracted difference
already rounded.
+ if fraction <= 0. {
+ return start;
+ }
+ if fraction >= 1. {
+ return end;
+ }
+ if start.is_sign_positive() != end.is_sign_positive() {
+ // Avoid the potentially overflowing difference of opposite-signed
endpoints.
+ start * (1. - fraction) + end * fraction
} else {
- // Same-sign subtraction is finite and avoids summing two near-maximum
terms.
- (x2 - x1).mul_add(ratio, x1)
+ // Same-sign subtraction is finite; rounding must not overshoot the
end.
+ let value = (end - start).mul_add(fraction, start);
+ if start < end {
+ value.min(end)
+ } else {
+ value.max(end)
+ }
+ }
+}
+
+/// Locates a value between distinct finite endpoints.
+fn interpolation_fraction(value: f64, start: f64, end: f64) -> f64 {
+ let width = end - start;
+ if width.is_finite() {
+ (value - start) / width
+ } else {
+ // Scaling is needed only when opposite-signed endpoints overflow
their difference.
+ (value * 0.5 - start * 0.5) / (end * 0.5 - start * 0.5)
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::interpolate;
+
+ #[test]
+ fn interpolation_preserves_endpoints_after_subtraction_rounds() {
+ let start = f64::EPSILON / 2.;
+ let end = 1_f64.next_up();
+ // end - start rounds to 1, and even an FMA at fraction 1 then returns
1, not end.
+ // The endpoint branch also retains the sign of a zero endpoint.
+ for (start, end) in [(start, end), (end, start), (-0., 1.), (-1., 0.)]
{
+ assert_eq!(interpolate(start, end, 0.).to_bits(), start.to_bits());
+ assert_eq!(interpolate(start, end, 1.).to_bits(), end.to_bits());
+ }
}
}
diff --git a/tests-integration/tests/serde_tests/tdigest.rs
b/tests-integration/tests/serde_tests/tdigest.rs
index 1b871ef..b90b1ca 100644
--- a/tests-integration/tests/serde_tests/tdigest.rs
+++ b/tests-integration/tests/serde_tests/tdigest.rs
@@ -247,7 +247,7 @@ fn
test_serialized_bytes_stable_for_full_and_merged_digests() {
let mut full_buffer = patterned_digest(200, 1_641, 0);
let bytes = full_buffer.serialize();
assert_eq!(bytes.len(), 2_864);
- assert_eq!(fnv1a(&bytes), 0x5c01_c50d_d1c8_fdbb);
+ assert_eq!(fnv1a(&bytes), 0xaffe_62f2_3716_59b8);
let mut left = patterned_digest(10, 201, 2);
let mut right = patterned_digest(10, 199, 3);
@@ -255,12 +255,12 @@ fn
test_serialized_bytes_stable_for_full_and_merged_digests() {
left.merge(&right);
let bytes = left.serialize();
assert_eq!(bytes.len(), 272);
- assert_eq!(fnv1a(&bytes), 0x7d2e_a927_9b9e_f559);
+ assert_eq!(fnv1a(&bytes), 0x6037_b570_d640_d31f);
for &(left_len, right_len, expected_len, expected_hash) in &[
- (8, 201, 272, 0x8522_1f3f_152f_24e5),
- (201, 8, 256, 0xe60d_1f6f_f4b0_73e0),
- (201, 401, 288, 0x4cb8_4037_5e68_ca4b),
+ (8, 201, 272, 0x4ec5_2661_c556_c10e),
+ (201, 8, 256, 0x2b43_056d_9a57_b2a5),
+ (201, 401, 288, 0x1e43_34b7_ba13_4fe0),
(401, 201, 288, 0x6f6d_e965_77a7_a53f),
] {
let mut left = patterned_digest(10, left_len, 2);
@@ -362,6 +362,80 @@ fn assert_invalid_tdigest(bytes: &[u8]) {
assert_eq!(error.kind(), datasketches::error::ErrorKind::InvalidData);
}
+#[test]
+fn test_deserialize_single_sample_requires_equal_extrema() {
+ // The multiple-value format can store a single centroid or a single
buffered value.
+ for buffered in [false, true] {
+ for is_f32 in [false, true] {
+ for (min, max) in [(0., 100.), (0., 50.), (50., 100.), (50., 50.)]
{
+ let mut bytes = vec![2, 1, 20, 100, 0, 0, 0, 0];
+ bytes.extend(u32::from(!buffered).to_le_bytes());
+ bytes.extend(u32::from(buffered).to_le_bytes());
+ if is_f32 {
+ for value in [min, max, 50.] {
+ bytes.extend((value as f32).to_le_bytes());
+ }
+ if !buffered {
+ bytes.extend(1_u32.to_le_bytes());
+ }
+ } else {
+ for value in [min, max, 50_f64] {
+ bytes.extend(value.to_le_bytes());
+ }
+ if !buffered {
+ bytes.extend(1_u64.to_le_bytes());
+ }
+ }
+ let result = if is_f32 {
+ TDigestMut::deserialize_f32(&bytes)
+ } else {
+ TDigestMut::deserialize(&bytes)
+ };
+ if min == max {
+ let mut digest = result.unwrap();
+ assert_eq!(digest.total_weight(), 1);
+ assert_eq!(digest.quantiles(&[0., 0.5, 1.]),
Some(vec![50.; 3]));
+ } else {
+ assert_eq!(
+ result.unwrap_err().kind(),
+ datasketches::error::ErrorKind::InvalidData
+ );
+ }
+ }
+ }
+ }
+}
+
+#[test]
+fn test_deserialize_compat_single_sample_requires_equal_extrema() {
+ for compact in [false, true] {
+ for (min, max) in [(0., 100.), (0., 50.), (50., 100.), (50., 50.)] {
+ let mut bytes = if compact { 2_u32 } else { 1_u32
}.to_be_bytes().to_vec();
+ bytes.extend(f64::to_be_bytes(min));
+ bytes.extend(f64::to_be_bytes(max));
+ if compact {
+ bytes.extend(100_f32.to_be_bytes());
+ bytes.extend(0_u32.to_be_bytes()); // unused capacities
+ bytes.extend(1_u16.to_be_bytes());
+ bytes.extend(1_f32.to_be_bytes());
+ bytes.extend(50_f32.to_be_bytes());
+ } else {
+ bytes.extend(100_f64.to_be_bytes());
+ bytes.extend(1_u32.to_be_bytes());
+ bytes.extend(1_f64.to_be_bytes());
+ bytes.extend(50_f64.to_be_bytes());
+ }
+ if min == max {
+ let mut digest = TDigestMut::deserialize(&bytes).unwrap();
+ assert_eq!(digest.total_weight(), 1);
+ assert_eq!(digest.quantiles(&[0., 0.5, 1.]), Some(vec![50.;
3]));
+ } else {
+ assert_invalid_tdigest(&bytes);
+ }
+ }
+ }
+}
+
#[test]
fn test_deserialize_rejects_unknown_or_conflicting_flags() {
let mut unknown = serialized_two_value_digest();
diff --git a/tests-integration/tests/tdigest_test/property.rs
b/tests-integration/tests/tdigest_test/property.rs
index 6e5fbc1..61577ef 100644
--- a/tests-integration/tests/tdigest_test/property.rs
+++ b/tests-integration/tests/tdigest_test/property.rs
@@ -51,6 +51,50 @@ fn assert_quantiles_are_monotonic(tdigest: &mut TDigestMut) {
}
}
+#[test]
+fn prop_finite_values_survive_partial_aggregation() {
+ fn property(bits: Vec<u64>) -> TestResult {
+ let mut values: Vec<_> = bits
+ .into_iter()
+ .map(f64::from_bits)
+ .filter(|value| value.is_finite())
+ .collect();
+ if values.is_empty() {
+ return TestResult::discard();
+ }
+ let mut digest = TDigestMut::default();
+ for chunk in values.chunks(64) {
+ let mut partial = TDigestMut::default();
+ for &value in chunk {
+ partial.update(value);
+ }
+
digest.merge(&TDigestMut::deserialize(&partial.serialize()).unwrap());
+ }
+ let mut digest = TDigestMut::deserialize(&digest.serialize()).unwrap();
+ assert_eq!(digest.total_weight(), values.len() as u64);
+ assert_quantiles_are_monotonic(&mut digest);
+
+ // Use the input samples as split points, independently of the
computed quantiles.
+ values.sort_unstable_by(f64::total_cmp);
+ values.dedup();
+ assert_eq!(digest.quantile(0.), values.first().copied());
+ assert_eq!(digest.quantile(1.), values.last().copied());
+ let cdf = digest.cdf(&values).unwrap();
+ assert!(cdf.is_sorted());
+ assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank)));
+ let pmf = digest.pmf(&values).unwrap();
+ assert!(pmf.iter().all(|p| (0.0..=1.0).contains(p)));
+ assert!((pmf.iter().sum::<f64>() - 1.).abs() < 1e-12);
+ TestResult::passed()
+ }
+
+ QuickCheck::new()
+ .tests(128)
+ .min_tests_passed(128)
+ .rng(Gen::new(512))
+ .quickcheck(property as fn(Vec<u64>) -> TestResult);
+}
+
#[test]
fn prop_quantile_is_non_decreasing_and_within_the_observed_range() {
fn property(values: Vec<u32>) -> TestResult {
diff --git a/tests-integration/tests/tdigest_test/sketch.rs
b/tests-integration/tests/tdigest_test/sketch.rs
index cae816b..428d96c 100644
--- a/tests-integration/tests/tdigest_test/sketch.rs
+++ b/tests-integration/tests/tdigest_test/sketch.rs
@@ -576,6 +576,44 @@ fn test_extreme_values_produce_finite_quantiles() {
}
}
+#[test]
+fn test_rank_interpolation_across_extreme_finite_values() {
+ for scale in [1., f64::MAX] {
+ let mut digest = TDigestMut::default();
+ digest.update(-scale);
+ digest.update(scale);
+ let digest = digest.freeze();
+ // The two centroid centers have ranks 1/4 and 3/4, independently of
scale.
+ let points = [-scale, -scale / 2., 0., scale / 2., scale];
+ for (point, expected) in points.into_iter().zip([0.25, 0.375, 0.5,
0.625, 0.75]) {
+ assert_that!(
+ digest.rank(point).unwrap(),
+ near(expected, 4. * f64::EPSILON)
+ );
+ }
+ }
+}
+
+#[test]
+fn test_tail_interpolation_across_extreme_finite_values() {
+ let max = f64::MAX;
+ let left = deserialize_with_centroids(100, -max, max, &[(max / 2., 10),
(max, 1)]);
+ let right = deserialize_with_centroids(100, -max, max, &[(-max, 1), (-max
/ 2., 10)]);
+ // Weight 3 is halfway from the left tail's endpoint (weight 1) to its
center (weight 5).
+ // The expected value is therefore (-MAX + MAX/2)/2 = -MAX/4; the right
tail mirrors it.
+ for (digest, rank, expected) in [(left, 3. / 11., -0.25), (right, 8. /
11., 0.25)] {
+ let digest = digest.freeze();
+ assert_that!(
+ digest.quantile(rank).unwrap() / max,
+ near(expected, 4. * f64::EPSILON)
+ );
+ assert_that!(
+ digest.rank(expected * max).unwrap(),
+ near(rank, 4. * f64::EPSILON)
+ );
+ }
+}
+
#[test]
fn test_estimate_repeat_values() {
let mut tdigest = TDigestMut::default();
@@ -585,11 +623,7 @@ 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.
+/// Builds a weighted image, including tail centroids unavailable through raw
updates.
fn deserialize_with_centroids(k: u16, min: f64, max: f64, centroids: &[(f64,
u64)]) -> TDigestMut {
const PREAMBLE_LONGS: u8 = 2;
const SERIAL_VERSION: u8 = 1;
@@ -663,6 +697,265 @@ fn test_quantile_handles_two_sample_last_centroid() {
assert_eq!(tdigest.quantile(0.75), Some(100.0));
}
+#[test]
+fn test_single_centroid_preserves_stored_tail_information() {
+ let mut digest = deserialize_with_centroids(100, 0., 100., &[(50., 10)]);
+ // Ten samples have a centroid center at weight 5. The tail spans weights
1 through 5,
+ // so weight 2.5 lies 3/8 of the way from the stored minimum to the mean.
+ assert_quantile_queries(
+ &digest,
+ &[0., 0.1, 0.25, 0.5, 0.75, 0.9, 1.],
+ &[0., 0., 18.75, 50., 81.25, 100., 100.],
+ );
+ let frozen = digest.clone().freeze();
+ let points = [-1., 0., 18.75, 50., 81.25, 100., 101.];
+ let expected = [0., 0.05, 0.25, 0.5, 0.75, 0.95, 1.];
+
+ // Repeated queries agree before and after creating a mutable view.
+ for _ in 0..2 {
+ for (&point, &rank) in points.iter().zip(&expected) {
+ assert_eq!(digest.rank(point), Some(rank));
+ assert_eq!(frozen.rank(point), Some(rank));
+ }
+ assert_eq!(&digest.cdf(&points).unwrap()[..points.len()], &expected);
+ }
+}
+
+#[test]
+fn test_rank_tails_stay_monotonic_after_updating_weighted_centroids() {
+ let points = [0., 0_f64.next_up(), 0.5, 1., 99., 99.5, 100.];
+ let expected = [0.5, 0.5, 0.5, 0.5, 11.5, 11.5, 11.5].map(|mass| mass /
12.);
+ for operation in 0..3 {
+ let mut digest = deserialize_with_centroids(100, 0., 100., &[(50.,
10)]);
+ let mut updates = TDigestMut::new(100).unwrap();
+ updates.update(1.);
+ updates.update(99.);
+ match operation {
+ 0 => {
+ digest.update(1.);
+ digest.update(99.);
+ }
+ 1 => digest.merge(&updates),
+ _ => digest = [digest, updates].into_iter().collect(),
+ }
+ let frozen = digest.clone().freeze();
+
+ // The new unit-weight tails lie inside the retained extrema. Their
half-sample
+ // centers leave no mass for interpolation between min/first or
last/max.
+ for (&point, &expected) in points.iter().zip(&expected) {
+ assert_eq!(digest.rank(point), Some(expected));
+ assert_eq!(frozen.rank(point), Some(expected));
+ }
+ let cdf = digest.cdf(&points).unwrap();
+ assert_eq!(cdf, frozen.cdf(&points).unwrap());
+ assert!(cdf.is_sorted(), "CDF at {points:?}: {cdf:?}");
+ let pmf = digest.pmf(&points).unwrap();
+ assert_eq!(pmf, frozen.pmf(&points).unwrap());
+ assert!(pmf.iter().all(|&mass| mass >= 0.), "PMF: {pmf:?}");
+ assert_that!(pmf.iter().sum::<f64>(), near(1., 1e-12));
+ }
+}
+
+#[test]
+fn test_small_compression_preserves_weighted_centroids_and_tie_order() {
+ // Exercise the no-merge bound at two k values and a tiny input whose K_2
normalizer
+ // is negative. Existing weighted centroids must remain intact in every
case.
+ for (k, middle_weight) in [(200, 95), (u16::MAX, 2), (u16::MAX, 32_762)] {
+ for reverse in [false, true] {
+ let mut digest = deserialize_with_centroids(
+ k,
+ -10.,
+ 10.,
+ &[(-10., 1), (-0., 1), (0., middle_weight), (1., 1), (10., 1)],
+ );
+ let mut bytes = digest.serialize();
+ if reverse {
+ bytes[5] |= 1 << 2;
+ }
+ let mut digest = TDigestMut::deserialize(&bytes).unwrap();
+ digest.update(0.);
+
+ // Buffered values precede retained centroids on ties, including
signed zeros.
+ let mut expected = deserialize_with_centroids(
+ k,
+ -10.,
+ 10.,
+ &[
+ (-10., 1),
+ (0., 1),
+ (-0., 1),
+ (0., middle_weight),
+ (1., 1),
+ (10., 1),
+ ],
+ )
+ .serialize();
+ if !reverse {
+ expected[5] |= 1 << 2;
+ }
+ assert_eq!(digest.total_weight(), middle_weight + 5);
+ assert_eq!(digest.serialize(), expected);
+ }
+ }
+}
+
+// Exercise both compression paths and inspect the mean directly; query
interpolation can hide
+// rounding errors. The heavy outer centroids keep the middle pair away from
the protected tails.
+fn merged_centroid_mean(centroids: [(f64, u64); 4], reverse: bool, owned:
bool) -> f64 {
+ let min = centroids[0].0;
+ let max = centroids[3].0;
+ let mut bytes = deserialize_with_centroids(10, min, max,
¢roids).serialize();
+ if reverse {
+ bytes[5] |= 1 << 2;
+ }
+ let mut digest = TDigestMut::deserialize(&bytes).unwrap();
+ if owned {
+ let other = deserialize_with_centroids(10, max, max, &[(max, 1)]);
+ digest = [digest, other].into_iter().collect();
+ } else {
+ digest.update(min);
+ }
+ let bytes = digest.serialize();
+ let merged_weight = centroids[1].1 + centroids[2].1;
+ let centroid = bytes[32..]
+ .chunks_exact(16)
+ .find(|c| u64::from_le_bytes(c[8..].try_into().unwrap()) ==
merged_weight)
+ .unwrap_or_else(|| {
+ panic!(
+ "middle pair {:?} did not merge (reverse={reverse},
owned={owned})",
+ ¢roids[1..3]
+ )
+ });
+ f64::from_le_bytes(centroid[..8].try_into().unwrap())
+}
+
+#[test]
+fn test_compression_preserves_a_small_centroid_weight() {
+ for weight in [1_u64 << 20, 1_u64 << 40, 1_u64 << 54] {
+ for reverse in [false, true] {
+ let mean = merged_centroid_mean(
+ [(0., weight), (1., weight), (1e300, 1), (f64::MAX, weight)],
+ reverse,
+ false,
+ );
+ // The exact mean is (weight * 1 + 1e300) / (weight + 1). The
first term is far
+ // below one ULP of the numerator, so this division is an
independent reference.
+ let expected = 1e300 / (weight + 1) as f64;
+ assert_that!(mean / expected, near(1., 4. * f64::EPSILON));
+ }
+ }
+}
+
+#[test]
+fn test_opposite_sign_merge_preserves_a_representable_correction() {
+ let outer_weight = 1_u64 << 56;
+ for magnitude in [1., f64::MAX] {
+ for exponent in [54, 55] {
+ let total_weight = 1_u64 << exponent;
+ // The exact mean is magnitude * (1 - 2 / total_weight). At 2^54 it
+ // rounds to the next value toward zero; at 2^55 it rounds back to
magnitude.
+ let expected = if exponent == 54 {
+ magnitude.next_down()
+ } else {
+ magnitude
+ };
+ for negative_heavy in [false, true] {
+ let (left_weight, right_weight, expected) = if negative_heavy {
+ (total_weight - 1, 1, -expected)
+ } else {
+ (1, total_weight - 1, expected)
+ };
+ for reverse in [false, true] {
+ for owned in [false, true] {
+ let mean = merged_centroid_mean(
+ [
+ (-f64::MAX, outer_weight),
+ (-magnitude, left_weight),
+ (magnitude, right_weight),
+ (f64::MAX, outer_weight),
+ ],
+ reverse,
+ owned,
+ );
+ assert_eq!(
+ mean, expected,
+ "magnitude={magnitude:?}, exponent={exponent},
negative_heavy={negative_heavy}, reverse={reverse}, owned={owned}"
+ );
+ }
+ }
+ }
+ }
+ }
+}
+
+#[test]
+fn test_compression_keeps_merged_means_within_their_endpoints() {
+ let tiny = f64::from_bits(1);
+ let heavy = 1_u64 << 54;
+ for (left, right) in [
+ (0., tiny),
+ (tiny, 2. * tiny),
+ (f64::MIN_POSITIVE.next_down(), f64::MIN_POSITIVE),
+ (1., 1_f64.next_up()),
+ (f64::MAX.next_down(), f64::MAX),
+ (f64::MAX / 4., f64::MAX),
+ (-f64::MAX, -f64::MAX / 4.),
+ (-2. * tiny, -tiny),
+ (-tiny, tiny),
+ (-f64::MAX, f64::MAX),
+ ] {
+ for (left_weight, right_weight) in [(1, 1), (1, heavy), (heavy, 1)] {
+ for reverse in [false, true] {
+ for owned in [false, true] {
+ let mean = merged_centroid_mean(
+ [
+ (-f64::MAX, heavy),
+ (left, left_weight),
+ (right, right_weight),
+ (f64::MAX, heavy),
+ ],
+ reverse,
+ owned,
+ );
+ assert!(
+ mean.is_finite() && (left..=right).contains(&mean),
+ "{mean:?} is outside [{left:?}, {right:?}]
(weights={left_weight}:{right_weight}, reverse={reverse}, owned={owned})"
+ );
+ }
+ }
+ }
+ }
+}
+
+#[test]
+fn test_signed_zero_and_symmetric_centroids_survive_interpolation() {
+ for magnitude in [0., f64::from_bits(1), 1., f64::MAX] {
+ let heavy = 1_u64 << 54;
+ for reverse in [false, true] {
+ let mean = merged_centroid_mean(
+ [
+ (-f64::MAX, heavy),
+ (-magnitude, 1),
+ (magnitude, 1),
+ (f64::MAX, heavy),
+ ],
+ reverse,
+ false,
+ );
+ // Equal and opposite contributions cancel exactly, including
subnormal values.
+ assert_eq!(mean, 0.);
+ }
+ }
+
+ let mut digest =
+ deserialize_with_centroids(100, -1., 1., &[(-1., 1), (-0., 2), (0.,
2), (1., 1)]);
+ assert_eq!(digest.rank(-0.), Some(0.5));
+ assert_eq!(digest.rank(0.), Some(0.5));
+ assert_quantile_queries(&digest, &[1. / 3., 0.5, 2. / 3.], &[0., 0., 0.]);
+ assert_eq!(digest.cdf(&[-0.]), Some(vec![0.5, 1.]));
+ assert_eq!(digest.pmf(&[0.]), Some(vec![0.5, 0.5]));
+}
+
#[test]
fn test_batch_quantiles_match_scalar_queries_in_input_order() {
let mut tdigest = TDigestMut::new(100).unwrap();
@@ -791,6 +1084,123 @@ fn
test_quantiles_stay_within_extrema_at_large_total_weights() {
}
}
+#[test]
+fn test_queries_stay_ordered_at_large_centroid_centers() {
+ for total_weight in [1_u64 << 53, 1_u64 << 54, u64::MAX] {
+ for centroids in [
+ vec![(10., total_weight - 5), (50., 1), (90., 4)],
+ vec![(0., 1), (90., total_weight - 1)],
+ ] {
+ let digest = deserialize_with_centroids(100, 0., 100.,
¢roids).freeze();
+ // Probe either side of a heavy centroid's center and the rounded
upper endpoint.
+ let ranks = [
+ 0.,
+ 0.5_f64.next_down(),
+ 0.5,
+ 0.5_f64.next_up(),
+ 1.0_f64.next_down(),
+ 1.,
+ ];
+ let quantiles = digest.quantiles(&ranks).unwrap();
+ assert_eq!(quantiles[0], 0.);
+ assert_eq!(quantiles[5], 100.);
+ assert!(
+ quantiles.is_sorted(),
+ "centroids {centroids:?}: {quantiles:?}"
+ );
+ assert!(quantiles.iter().all(|value|
(0.0..=100.0).contains(value)));
+
+ let cdf = digest.cdf(&[0., 5., 10., 50., 90., 95., 100.]).unwrap();
+ assert!(cdf.is_sorted(), "centroids {centroids:?}: {cdf:?}");
+ assert!(cdf.iter().all(|rank| (0.0..=1.0).contains(rank)));
+ }
+ }
+}
+
+#[test]
+fn test_quantile_right_tail_uses_the_same_center_as_rank() {
+ for (total_weight, last_weight, rank) in [
+ (1_u64 << 53, 5, 1. - f64::EPSILON),
+ ((1_u64 << 53) + 1, 4, 1_f64.next_down()),
+ (1_u64 << 54, 9, 1. - f64::EPSILON),
+ (u64::MAX, 4097, 1_f64.next_down()),
+ ] {
+ let mut digest = deserialize_with_centroids(
+ 100,
+ 0.,
+ 100.,
+ &[(0., total_weight - last_weight), (90., last_weight)],
+ );
+ // These ranks land exactly on the rounded mass N - last_weight/2.
Quantile must
+ // return the centroid mean there, before interpolating toward the
stored maximum.
+ assert_eq!(digest.rank(90.), Some(rank));
+ assert_quantile_queries(&digest, &[rank, rank, 1.], &[90., 90., 100.]);
+ }
+
+ // At N = 2^53 + 2 the last center and N - 1 round to the same mass. The
zero-width
+ // tail must stay finite and ordered, just as a two-sample tail does at
smaller counts.
+ let total_weight = (1_u64 << 53) + 2;
+ let digest = deserialize_with_centroids(100, 0., 100., &[(0., total_weight
- 3), (90., 3)]);
+ let rank = 1. - f64::EPSILON;
+ let ranks = [rank.next_down(), rank, rank.next_up(), 1.];
+ let quantiles = digest.clone().quantiles(&ranks).unwrap();
+ assert!(quantiles.is_sorted());
+ assert!(quantiles.iter().all(|value| (0.0..=100.0).contains(value)));
+ assert_eq!(quantiles[3], 100.);
+ assert_quantile_queries(&digest, &ranks, &quantiles);
+}
+
+#[test]
+fn test_quantile_uses_max_when_last_center_rounds_down() {
+ let total_weight = (1_u64 << 53) - 1;
+ let first_weight = (total_weight - 1) / 2;
+ let middle_weight = total_weight - first_weight - 1;
+ let digest = deserialize_with_centroids(
+ 100,
+ 0.,
+ 100.,
+ &[(10., first_weight), (50., middle_weight), (90., 1)],
+ );
+
+ // The last center's +0.5 ties down to the even count 2^53 - 2. The
highest rank
+ // below 1 lands there and exhausts the scan; the adjacent rank stays
inside it.
+ assert_quantile_queries(
+ &digest,
+ &[1. - f64::EPSILON, 1_f64.next_down()],
+ &[90_f64.next_down(), 100.],
+ );
+}
+
+#[test]
+fn test_rank_max_uses_the_same_rounded_mass_as_last_centroid() {
+ // The last half-sample center can round up or down. At N = 2^53 - 2 it
rounds to N,
+ // while 1 - 0.5/N rounds below 1 and would give max a lower rank than the
last mean.
+ for (total_weight, expected) in [
+ ((1_u64 << 53) - 2, 1.),
+ ((1_u64 << 53) - 1, 1_f64.next_down()),
+ (u64::MAX, 1.),
+ ] {
+ let mut digest =
+ deserialize_with_centroids(100, 0., 100., &[(0., total_weight -
1), (90., 1)]);
+ let frozen = digest.clone().freeze();
+ let points = [90., 90_f64.next_up(), 95., 100_f64.next_down(), 100.];
+ let cdf = digest.cdf(&points).unwrap();
+ assert_eq!(&cdf[..points.len()], &[expected; 5]);
+ assert_eq!(cdf, frozen.cdf(&points).unwrap());
+ assert!(frozen.pmf(&points).unwrap().iter().all(|&mass| mass >= 0.));
+ }
+}
+
+#[test]
+fn test_rank_stays_monotonic_at_the_right_tail_boundary() {
+ let digest = deserialize_with_centroids(100, 0., 300., &[(0., 1), (100.,
29)]).freeze();
+ let points = [100., 100_f64.next_up(), 150., 300.];
+ let cdf = digest.cdf(&points).unwrap();
+ // The last center is at weight 1 + 29/2 = 15.5. Entering the tail must
not lower its rank.
+ assert_eq!(cdf[0], 15.5 / 30.);
+ assert!(cdf.is_sorted(), "CDF at {points:?}: {cdf:?}");
+}
+
#[test]
fn test_batch_quantiles_reject_invalid_ranks() {
for values in [&[][..], &[1.0, 2.0, 3.0][..]] {
---------------------------------------------------------------------
To unsubscribe, e-mail: [email protected]
For additional commands, e-mail: [email protected]
