diff --git a/CHANGELOG.md b/CHANGELOG.md index c5c6553b..995aa648 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 5a7d1c65..ef39ffdb 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, 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::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: +/// . +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 1b871ef8..b90b1ca6 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 6e5fbc14..61577efd 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) -> 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::() - 1.).abs() < 1e-12); + TestResult::passed() + } + + QuickCheck::new() + .tests(128) + .min_tests_passed(128) + .rng(Gen::new(512)) + .quickcheck(property as fn(Vec) -> TestResult); +} + #[test] fn prop_quantile_is_non_decreasing_and_within_the_observed_range() { fn property(values: Vec) -> TestResult { diff --git a/tests-integration/tests/tdigest_test/sketch.rs b/tests-integration/tests/tdigest_test/sketch.rs index cae816bb..428d96c7 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::(), 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][..]] {