The sketch is a list of centroids, each holding a mean and a weight, kept sorted by mean. + * Together they approximate the distribution of everything ever added. What makes a t-digest more + * than a histogram is the rule that decides how big a centroid is allowed to be: instead of a fixed + * bucket width, the size limit is expressed through a scale function + * + *
+ * k(q) = compression / (2 * pi) * asin(2q - 1) + *+ * + *
and a centroid may absorb points only while it spans at most one unit of {@code k}. Because + * {@code asin} is steep near {@code q = 0} and {@code q = 1}, centroids near the tails are forced to + * stay tiny - often a single point - while centroids near the median are allowed to grow large. + * That is exactly the trade one wants in practice: the rank error is of the order of + * {@code 1 / compression} around the median and shrinks as {@code q} approaches 0 or 1, so + * {@code p99.9} comes out nearly exact. On smooth data the interpolation between centroids brings + * the error down by another order of magnitude, and the whole sketch still fits in a few kilobytes. + * + *
Insertions are buffered and folded into the centroid list in batches, which keeps the amortised + * cost per sample low. The minimum and the maximum are tracked separately and reported exactly. + * + *
| Operation | Complexity |
|---|---|
| {@link #add(double)} | O(log d) amortised |
| {@link #quantile(double)}, {@link #cdf(double)} | O(d) |
| {@link #merge(TDigest)} | O(d log d) |
| memory | O(d), independent of the stream length |
{@code
+ * TDigest digest = new TDigest(100);
+ * for (double latency : latencies) {
+ * digest.add(latency);
+ * }
+ * digest.quantile(0.5); // median
+ * digest.quantile(0.999); // deep tail, the case t-digest is built for
+ * digest.cdf(250.0); // share of samples below 250
+ *
+ * // Sketches computed on different shards combine the way partial sums do.
+ * shardA.merge(shardB);
+ * }
+ *
+ * Query methods fold the pending insertion buffer in before answering, so they mutate internal + * state; they are not safe to call concurrently. This class is not thread-safe. + * + * @see P2QuantileEstimator for a constant-memory estimator of one fixed quantile + * @see T. Dunning, O. Ertl, Computing extremely accurate quantiles using t-digests + */ +public final class TDigest { + + /** Below roughly this value the scale function stops leaving room for a useful number of centroids. */ + private static final double MIN_COMPRESSION = 10.0; + + private static final double TWO_PI = 2.0 * Math.PI; + + /** Buffered points per centroid slot; larger buffers mean fewer, bigger merge passes. */ + private static final int BUFFER_FACTOR = 5; + + /** + * How many centroid slots to reserve per unit of compression. Repeated merges make centroids + * atomic, so they pack less tightly than the scale function alone would suggest: in practice the + * sketch settles around {@code 0.6 * compression} centroids rather than the ideal + * {@code compression / 2}. Twice the compression leaves ample headroom. + */ + private static final int CENTROID_SLOTS_PER_UNIT = 2; + + /** Below this length the parallel sort switches to insertion sort. */ + private static final int INSERTION_SORT_THRESHOLD = 16; + + private final double compression; + + private final double[] centroidMeans; + private final double[] centroidWeights; + private int centroidCount; + + private final double[] bufferMeans; + private final double[] bufferWeights; + private int bufferCount; + + private final double[] scratchMeans; + private final double[] scratchWeights; + + private double totalWeight; + private double min = Double.POSITIVE_INFINITY; + private double max = Double.NEGATIVE_INFINITY; + + /** + * Creates an empty digest with a compression of 100, a good default for most streams. + */ + public TDigest() { + this(100.0); + } + + /** + * Creates an empty digest. + * + * @param compression the size/accuracy trade-off; around {@code 0.6 * compression} centroids end + * up being kept, so larger values cost more memory and give smaller quantile errors + * @throws IllegalArgumentException if {@code compression} is smaller than 10 or not finite + */ + public TDigest(double compression) { + if (!(compression >= MIN_COMPRESSION) || !Double.isFinite(compression)) { + throw new IllegalArgumentException("The compression must be finite and at least " + MIN_COMPRESSION + ", but was " + compression); + } + this.compression = compression; + + int scale = (int) Math.ceil(compression) + 2; + int centroidCapacity = CENTROID_SLOTS_PER_UNIT * scale; + this.centroidMeans = new double[centroidCapacity]; + this.centroidWeights = new double[centroidCapacity]; + + int bufferCapacity = BUFFER_FACTOR * scale; + this.bufferMeans = new double[bufferCapacity]; + this.bufferWeights = new double[bufferCapacity]; + + this.scratchMeans = new double[centroidCapacity + bufferCapacity]; + this.scratchWeights = new double[centroidCapacity + bufferCapacity]; + } + + /** + * Adds one sample of unit weight. + * + * @param value the sample to add + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public void add(double value) { + add(value, 1.0); + } + + /** + * Adds a sample that stands for several observations of the same value. + * + * @param value the sample to add + * @param weight how many observations the sample represents, strictly positive + * @throws IllegalArgumentException if {@code value} is not finite or {@code weight} is not strictly positive + */ + public void add(double value, double weight) { + if (!Double.isFinite(value)) { + throw new IllegalArgumentException("Samples must be finite, but was " + value); + } + if (!(weight > 0.0) || !Double.isFinite(weight)) { + throw new IllegalArgumentException("The weight must be finite and positive, but was " + weight); + } + if (bufferCount == bufferMeans.length) { + flushBuffer(); + } + bufferMeans[bufferCount] = value; + bufferWeights[bufferCount] = weight; + bufferCount++; + min = Math.min(min, value); + max = Math.max(max, value); + } + + /** + * Adds every given sample with unit weight. + * + * @param values the samples to add + * @throws IllegalArgumentException if any value is NaN or infinite + * @throws NullPointerException if {@code values} is {@code null} + */ + public void addAll(double... values) { + for (double value : values) { + add(value); + } + } + + /** + * Folds another digest into this one. The result approximates the union of the two streams, which + * is what makes t-digests usable as a map-reduce style aggregate. + * + * @param other the digest to absorb; it is left unchanged apart from its pending buffer being folded in + * @throws NullPointerException if {@code other} is {@code null} + */ + public void merge(TDigest other) { + other.flushBuffer(); + for (int i = 0; i < other.centroidCount; i++) { + add(other.centroidMeans[i], other.centroidWeights[i]); + } + if (!other.isEmpty()) { + min = Math.min(min, other.min); + max = Math.max(max, other.max); + } + } + + /** + * Estimates the value below which the given fraction of the stream lies. + * + * @param q the requested quantile, between 0 and 1 inclusive + * @return the estimated quantile; exactly {@link #min()} for {@code q == 0} and {@link #max()} for {@code q == 1} + * @throws IllegalArgumentException if {@code q} is outside {@code [0, 1]} or is NaN + * @throws IllegalStateException if no sample has been added yet + */ + public double quantile(double q) { + if (!(q >= 0.0) || !(q <= 1.0)) { + throw new IllegalArgumentException("The quantile probability must lie in [0, 1], but was " + q); + } + flushBuffer(); + requireNonEmpty(); + if (centroidCount == 1) { + return centroidMeans[0]; + } + + double index = q * totalWeight; + if (index <= 0.0) { + return min; + } + if (index >= totalWeight) { + return max; + } + + // The centroids are treated as knots of a piecewise linear quantile function, each sitting at + // the centre of the weight it carries, with min and max closing the two ends. + double previousValue = min; + double previousIndex = 0.0; + double currentIndex = 0.0; + for (int i = 0; i < centroidCount; i++) { + currentIndex += i == 0 ? centroidWeights[0] / 2.0 : (centroidWeights[i - 1] + centroidWeights[i]) / 2.0; + if (index < currentIndex) { + return interpolate(previousIndex, previousValue, currentIndex, centroidMeans[i], index); + } + previousValue = centroidMeans[i]; + previousIndex = currentIndex; + } + return interpolate(previousIndex, previousValue, totalWeight, max, index); + } + + /** + * Estimates the fraction of the stream that is smaller than the given value, i.e. the empirical + * cumulative distribution function. + * + * @param value the value to look up + * @return a number in {@code [0, 1]} + * @throws IllegalArgumentException if {@code value} is NaN + * @throws IllegalStateException if no sample has been added yet + */ + public double cdf(double value) { + if (Double.isNaN(value)) { + throw new IllegalArgumentException("The value must not be NaN"); + } + flushBuffer(); + requireNonEmpty(); + if (value < min) { + return 0.0; + } + if (value > max) { + return 1.0; + } + if (max <= min) { + return 0.5; + } + + double previousValue = min; + double previousIndex = 0.0; + double currentIndex = 0.0; + for (int i = 0; i < centroidCount; i++) { + currentIndex += i == 0 ? centroidWeights[0] / 2.0 : (centroidWeights[i - 1] + centroidWeights[i]) / 2.0; + if (value < centroidMeans[i]) { + return interpolate(previousValue, previousIndex, centroidMeans[i], currentIndex, value) / totalWeight; + } + previousValue = centroidMeans[i]; + previousIndex = currentIndex; + } + return interpolate(previousValue, previousIndex, max, totalWeight, value) / totalWeight; + } + + /** + * Returns the smallest sample ever added, tracked exactly. + * + * @return the minimum + * @throws IllegalStateException if no sample has been added yet + */ + public double min() { + requireAnySample(); + return min; + } + + /** + * Returns the largest sample ever added, tracked exactly. + * + * @return the maximum + * @throws IllegalStateException if no sample has been added yet + */ + public double max() { + requireAnySample(); + return max; + } + + /** + * Returns the total weight of everything added, which for unweighted input is the number of + * samples. + * + * @return the accumulated weight + */ + public double totalWeight() { + double pending = 0.0; + for (int i = 0; i < bufferCount; i++) { + pending += bufferWeights[i]; + } + return totalWeight + pending; + } + + /** + * Returns how many centroids the sketch currently holds, after folding in pending insertions. + * + * @return the number of centroids, bounded by roughly {@code compression / 2} + */ + public int centroidCount() { + flushBuffer(); + return centroidCount; + } + + /** + * Returns the compression parameter. + * + * @return the value given at construction time + */ + public double compression() { + return compression; + } + + /** + * Tells whether the digest holds no samples. + * + * @return {@code true} if nothing has been added + */ + public boolean isEmpty() { + return centroidCount == 0 && bufferCount == 0; + } + + /** + * Forgets every sample. + */ + public void reset() { + centroidCount = 0; + bufferCount = 0; + totalWeight = 0.0; + min = Double.POSITIVE_INFINITY; + max = Double.NEGATIVE_INFINITY; + } + + @Override + public String toString() { + return "TDigest{compression=" + compression + ", weight=" + totalWeight() + ", centroids=" + (isEmpty() ? 0 : centroidCount()) + '}'; + } + + /** + * Folds every buffered point into the centroid list. + * + *
Buffer and centroids are concatenated, sorted by mean and then swept once from the smallest + * mean upwards. The sweep keeps filling the current output centroid while doing so would not push + * it past the weight limit given by the scale function at the quantile the centroid starts at; + * otherwise the centroid is closed and a new one is opened. + */ + private void flushBuffer() { + if (bufferCount == 0) { + return; + } + + int size = centroidCount + bufferCount; + System.arraycopy(centroidMeans, 0, scratchMeans, 0, centroidCount); + System.arraycopy(centroidWeights, 0, scratchWeights, 0, centroidCount); + System.arraycopy(bufferMeans, 0, scratchMeans, centroidCount, bufferCount); + System.arraycopy(bufferWeights, 0, scratchWeights, centroidCount, bufferCount); + bufferCount = 0; + sortByMean(scratchMeans, scratchWeights, 0, size - 1); + + double weight = 0.0; + for (int i = 0; i < size; i++) { + weight += scratchWeights[i]; + } + totalWeight = weight; + + int last = 0; + centroidMeans[0] = scratchMeans[0]; + centroidWeights[0] = scratchWeights[0]; + double weightOfClosedCentroids = 0.0; + double weightLimit = totalWeight * quantileLimit(0.0); + + for (int i = 1; i < size; i++) { + double projected = centroidWeights[last] + scratchWeights[i]; + boolean roomLeft = weightOfClosedCentroids + projected <= weightLimit; + if (roomLeft || last == centroidMeans.length - 1) { + centroidWeights[last] = projected; + centroidMeans[last] += (scratchMeans[i] - centroidMeans[last]) * scratchWeights[i] / projected; + } else { + weightOfClosedCentroids += centroidWeights[last]; + weightLimit = totalWeight * quantileLimit(weightOfClosedCentroids / totalWeight); + last++; + centroidMeans[last] = scratchMeans[i]; + centroidWeights[last] = scratchWeights[i]; + } + } + centroidCount = last + 1; + } + + /** + * Returns the quantile at which the centroid starting at quantile {@code q} must be closed, that + * is {@code kInverse(k(q) + 1)} for the scale function {@code k}. + */ + private double quantileLimit(double q) { + // Rounding in the running weight sum could push the ratio a hair outside [0, 1], which would + // turn asin into NaN and silently disable the size limit. + double bounded = Math.min(1.0, Math.max(0.0, q)); + double k = compression / TWO_PI * Math.asin(2.0 * bounded - 1.0) + 1.0; + double angle = k * TWO_PI / compression; + if (angle >= Math.PI / 2.0) { + return 1.0; + } + if (angle <= -Math.PI / 2.0) { + return 0.0; + } + return (Math.sin(angle) + 1.0) / 2.0; + } + + private void requireNonEmpty() { + if (centroidCount == 0) { + throw new IllegalStateException("The digest has not seen any sample yet"); + } + } + + private void requireAnySample() { + if (isEmpty()) { + throw new IllegalStateException("The digest has not seen any sample yet"); + } + } + + /** + * Linear interpolation between two knots, degenerating to the right knot when they coincide. + */ + private static double interpolate(double x0, double y0, double x1, double y1, double x) { + if (x1 <= x0) { + return y1; + } + return y0 + (x - x0) / (x1 - x0) * (y1 - y0); + } + + /** + * Sorts two parallel arrays by the values of the first one, using quicksort with a + * median-of-three pivot and an insertion sort for short ranges. + */ + private static void sortByMean(double[] means, double[] weights, int from, int to) { + int low = from; + int high = to; + while (low < high) { + if (high - low < INSERTION_SORT_THRESHOLD) { + insertionSort(means, weights, low, high); + return; + } + int pivotIndex = partition(means, weights, low, high); + // Recurse into the smaller half and loop on the larger one to keep the stack shallow. + if (pivotIndex - low < high - pivotIndex) { + sortByMean(means, weights, low, pivotIndex - 1); + low = pivotIndex + 1; + } else { + sortByMean(means, weights, pivotIndex + 1, high); + high = pivotIndex - 1; + } + } + } + + private static int partition(double[] means, double[] weights, int low, int high) { + int middle = low + (high - low) / 2; + if (means[middle] < means[low]) { + swap(means, weights, low, middle); + } + if (means[high] < means[low]) { + swap(means, weights, low, high); + } + if (means[high] < means[middle]) { + swap(means, weights, middle, high); + } + swap(means, weights, middle, high - 1); + double pivot = means[high - 1]; + + int store = low; + for (int i = low; i < high - 1; i++) { + if (means[i] < pivot) { + swap(means, weights, store, i); + store++; + } + } + swap(means, weights, store, high - 1); + return store; + } + + private static void insertionSort(double[] means, double[] weights, int low, int high) { + for (int i = low + 1; i <= high; i++) { + double mean = means[i]; + double weight = weights[i]; + int j = i - 1; + while (j >= low && means[j] > mean) { + means[j + 1] = means[j]; + weights[j + 1] = weights[j]; + j--; + } + means[j + 1] = mean; + weights[j + 1] = weight; + } + } + + private static void swap(double[] means, double[] weights, int i, int j) { + double mean = means[i]; + means[i] = means[j]; + means[j] = mean; + double weight = weights[i]; + weights[i] = weights[j]; + weights[j] = weight; + } + + /** + * Returns the centroid means, smallest first. Exposed for tests and for inspecting the shape of + * the sketch; the returned array is a copy. + * + * @return a copy of the centroid means + */ + double[] centroidMeansSnapshot() { + flushBuffer(); + return Arrays.copyOf(centroidMeans, centroidCount); + } + + /** + * Returns the centroid weights, aligned with {@link #centroidMeansSnapshot()}. + * + * @return a copy of the centroid weights + */ + double[] centroidWeightsSnapshot() { + flushBuffer(); + return Arrays.copyOf(centroidWeights, centroidCount); + } +} diff --git a/src/test/java/com/thealgorithms/streaming/TDigestTest.java b/src/test/java/com/thealgorithms/streaming/TDigestTest.java new file mode 100644 index 000000000000..6968513824c9 --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/TDigestTest.java @@ -0,0 +1,391 @@ +package com.thealgorithms.streaming; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +import java.util.Arrays; +import java.util.Random; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.ValueSource; + +class TDigestTest { + + /** + * Fraction of the sample that is strictly smaller than the given value, i.e. the rank the digest + * actually hit. Comparing ranks rather than values is the meaningful way to score a quantile + * sketch: it is insensitive to how steep the distribution happens to be. + */ + private static double trueRank(double[] sortedValues, double value) { + int index = Arrays.binarySearch(sortedValues, value); + if (index < 0) { + index = -(index + 1); + } + return (double) index / sortedValues.length; + } + + private static double[] gaussianSample(int count, long seed) { + Random random = new Random(seed); + double[] values = new double[count]; + for (int i = 0; i < count; i++) { + values[i] = random.nextGaussian(); + } + return values; + } + + private static double[] uniformSample(int count, long seed) { + Random random = new Random(seed); + double[] values = new double[count]; + for (int i = 0; i < count; i++) { + values[i] = random.nextDouble(); + } + return values; + } + + @ParameterizedTest + @ValueSource(doubles = {9.0, 0.0, -5.0, Double.NaN, Double.POSITIVE_INFINITY}) + void rejectsInvalidCompression(double compression) { + assertThrows(IllegalArgumentException.class, () -> new TDigest(compression)); + } + + @Test + void queriesBeforeTheFirstSampleFail() { + TDigest digest = new TDigest(); + assertTrue(digest.isEmpty()); + assertEquals(0.0, digest.totalWeight()); + assertEquals(100.0, digest.compression()); + assertThrows(IllegalStateException.class, () -> digest.quantile(0.5)); + assertThrows(IllegalStateException.class, () -> digest.cdf(0.5)); + assertThrows(IllegalStateException.class, digest::min); + assertThrows(IllegalStateException.class, digest::max); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + TDigest digest = new TDigest(); + assertThrows(IllegalArgumentException.class, () -> digest.add(value)); + } + + @ParameterizedTest + @ValueSource(doubles = {0.0, -1.0, Double.NaN, Double.POSITIVE_INFINITY}) + void rejectsInvalidWeights(double weight) { + TDigest digest = new TDigest(); + assertThrows(IllegalArgumentException.class, () -> digest.add(1.0, weight)); + } + + @ParameterizedTest + @ValueSource(doubles = {-0.1, 1.1, Double.NaN}) + void rejectsInvalidQuantiles(double q) { + TDigest digest = new TDigest(); + digest.add(1.0); + assertThrows(IllegalArgumentException.class, () -> digest.quantile(q)); + } + + @Test + void rejectsNaNInCdf() { + TDigest digest = new TDigest(); + digest.add(1.0); + assertThrows(IllegalArgumentException.class, () -> digest.cdf(Double.NaN)); + } + + @Test + void aSingleSampleAnswersEveryQuantile() { + TDigest digest = new TDigest(); + digest.add(4.0); + assertEquals(4.0, digest.quantile(0.0)); + assertEquals(4.0, digest.quantile(0.5)); + assertEquals(4.0, digest.quantile(1.0)); + assertEquals(4.0, digest.min()); + assertEquals(4.0, digest.max()); + assertEquals(1.0, digest.totalWeight()); + assertEquals(1, digest.centroidCount()); + } + + @Test + void aConstantStreamCollapsesIntoOneValue() { + TDigest digest = new TDigest(); + for (int i = 0; i < 10_000; i++) { + digest.add(3.5); + } + assertEquals(3.5, digest.quantile(0.1), 1e-12); + assertEquals(3.5, digest.quantile(0.9), 1e-12); + assertEquals(0.5, digest.cdf(3.5), 1e-12); + } + + @Test + void extremesAreExact() { + double[] values = gaussianSample(20_000, 11L); + TDigest digest = new TDigest(); + digest.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + assertEquals(sorted[0], digest.min(), 0.0); + assertEquals(sorted[sorted.length - 1], digest.max(), 0.0); + assertEquals(sorted[0], digest.quantile(0.0), 0.0); + assertEquals(sorted[sorted.length - 1], digest.quantile(1.0), 0.0); + assertEquals(0.0, digest.cdf(sorted[0] - 1.0)); + assertEquals(1.0, digest.cdf(sorted[sorted.length - 1] + 1.0)); + } + + @ParameterizedTest + @ValueSource(doubles = {0.001, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999}) + @DisplayName("hits the requested rank of a uniform stream") + void hitsTheRequestedRankOnUniformData(double q) { + double[] values = uniformSample(100_000, 20240517L); + TDigest digest = new TDigest(100.0); + digest.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + assertEquals(q, trueRank(sorted, digest.quantile(q)), 0.01, "requested q=" + q); + } + + @ParameterizedTest + @ValueSource(doubles = {0.001, 0.01, 0.1, 0.5, 0.9, 0.99, 0.999}) + @DisplayName("is far more accurate in the tails than in the middle, by design") + void isMostAccurateInTheTails(double q) { + double[] values = gaussianSample(200_000, 4242L); + TDigest digest = new TDigest(200.0); + digest.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + double rankError = Math.abs(trueRank(sorted, digest.quantile(q)) - q); + + // The tolerance follows the shape of the scale function: it shrinks as q approaches 0 or 1. + double allowedRankError = 0.01 * Math.max(0.05, 4.0 * q * (1.0 - q)); + assertTrue(rankError <= allowedRankError, "q=" + q + " rank error " + rankError + " exceeded " + allowedRankError); + } + + @Test + void quantilesAreMonotonic() { + TDigest digest = new TDigest(); + digest.addAll(gaussianSample(50_000, 777L)); + + double previous = digest.quantile(0.0); + for (int i = 1; i <= 1_000; i++) { + double current = digest.quantile(i / 1_000.0); + assertTrue(current >= previous, "quantile decreased at q=" + i / 1_000.0); + previous = current; + } + } + + @Test + @DisplayName("the cdf inverts the quantile function") + void cdfAgreesWithQuantile() { + TDigest digest = new TDigest(200.0); + digest.addAll(uniformSample(100_000, 314159L)); + + for (int step = 1; step <= 19; step++) { + double q = step * 0.05; + assertEquals(q, digest.cdf(digest.quantile(q)), 0.02, "at q=" + q); + } + } + + @Test + void cdfMatchesTheEmpiricalDistribution() { + double[] values = uniformSample(50_000, 2718L); + TDigest digest = new TDigest(); + digest.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + for (int step = 1; step <= 19; step++) { + double x = step * 0.05; + assertEquals(trueRank(sorted, x), digest.cdf(x), 0.02, "at x=" + x); + } + } + + @Test + @DisplayName("a weighted sample stands for that many observations") + void supportsWeightedSamples() { + TDigest weighted = new TDigest(); + TDigest repeated = new TDigest(); + for (int i = 0; i < 100; i++) { + weighted.add(i, 10.0); + for (int repeat = 0; repeat < 10; repeat++) { + repeated.add(i); + } + } + + assertEquals(1_000.0, weighted.totalWeight(), 1e-9); + assertEquals(repeated.totalWeight(), weighted.totalWeight(), 1e-9); + assertEquals(repeated.quantile(0.5), weighted.quantile(0.5), 2.0); + } + + @Test + @DisplayName("merging shard digests approximates a digest of the whole stream") + void mergesLikeAPartialAggregate() { + double[] values = gaussianSample(120_000, 5150L); + TDigest whole = new TDigest(200.0); + TDigest shardA = new TDigest(200.0); + TDigest shardB = new TDigest(200.0); + TDigest shardC = new TDigest(200.0); + for (int i = 0; i < values.length; i++) { + whole.add(values[i]); + if (i % 3 == 0) { + shardA.add(values[i]); + } else if (i % 3 == 1) { + shardB.add(values[i]); + } else { + shardC.add(values[i]); + } + } + + shardA.merge(shardB); + shardA.merge(shardC); + + assertEquals(whole.totalWeight(), shardA.totalWeight(), 1e-9); + assertEquals(whole.min(), shardA.min(), 0.0); + assertEquals(whole.max(), shardA.max(), 0.0); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + for (double q : new double[] {0.01, 0.1, 0.5, 0.9, 0.99}) { + assertEquals(q, trueRank(sorted, shardA.quantile(q)), 0.01, "merged digest at q=" + q); + } + } + + @Test + void mergingAnEmptyDigestChangesNothing() { + TDigest digest = new TDigest(); + digest.addAll(1.0, 2.0, 3.0, 4.0, 5.0); + double median = digest.quantile(0.5); + + digest.merge(new TDigest()); + assertEquals(median, digest.quantile(0.5), 0.0); + assertEquals(5.0, digest.totalWeight()); + } + + @ParameterizedTest + @ValueSource(doubles = {20.0, 100.0, 500.0}) + @DisplayName("the sketch stays small no matter how long the stream is") + void keepsTheCentroidCountBounded(double compression) { + TDigest digest = new TDigest(compression); + digest.addAll(gaussianSample(200_000, 8L)); + + assertTrue(digest.centroidCount() <= compression, "kept " + digest.centroidCount() + " centroids for compression " + compression); + assertEquals(digest.centroidCount(), digest.centroidMeansSnapshot().length); + assertEquals(digest.centroidCount(), digest.centroidWeightsSnapshot().length); + + double[] means = digest.centroidMeansSnapshot(); + for (int i = 1; i < means.length; i++) { + assertTrue(means[i - 1] <= means[i], "centroids are not sorted at index " + i); + } + assertEquals(digest.totalWeight(), Arrays.stream(digest.centroidWeightsSnapshot()).sum(), 1e-6); + } + + @Test + @DisplayName("a higher compression buys accuracy") + void higherCompressionIsMoreAccurate() { + double[] values = gaussianSample(200_000, 6060L); + double[] sorted = values.clone(); + Arrays.sort(sorted); + + TDigest coarse = new TDigest(20.0); + TDigest fine = new TDigest(500.0); + coarse.addAll(values); + fine.addAll(values); + + double coarseError = Math.abs(trueRank(sorted, coarse.quantile(0.5)) - 0.5); + double fineError = Math.abs(trueRank(sorted, fine.quantile(0.5)) - 0.5); + assertTrue(fineError <= coarseError, "coarse=" + coarseError + " fine=" + fineError); + } + + @Test + void handlesSortedInput() { + TDigest digest = new TDigest(); + for (int i = 1; i <= 100_000; i++) { + digest.add(i); + } + assertEquals(50_000.0, digest.quantile(0.5), 1_000.0); + assertEquals(1.0, digest.min()); + assertEquals(100_000.0, digest.max()); + } + + @Test + void resetForgetsEverything() { + TDigest digest = new TDigest(); + digest.addAll(gaussianSample(1_000, 3L)); + digest.reset(); + + assertTrue(digest.isEmpty()); + assertEquals(0.0, digest.totalWeight()); + assertThrows(IllegalStateException.class, () -> digest.quantile(0.5)); + + digest.addAll(1.0, 2.0, 3.0); + assertEquals(3.0, digest.totalWeight()); + assertEquals(2.0, digest.quantile(0.5), 0.5); + } + + @Test + void toStringMentionsTheState() { + TDigest digest = new TDigest(); + assertTrue(digest.toString().contains("centroids=0"), digest.toString()); + digest.addAll(1.0, 2.0, 3.0); + assertFalse(digest.isEmpty()); + assertTrue(digest.toString().contains("weight=3.0"), digest.toString()); + } + + @Test + @DisplayName("the very edges of the quantile function interpolate towards the tracked extremes") + void interpolatesTowardsTheExtremes() { + double[] values = uniformSample(100_000, 161803L); + TDigest digest = new TDigest(); + digest.addAll(values); + + double[] sorted = values.clone(); + Arrays.sort(sorted); + double deepTail = digest.quantile(0.999999); + assertTrue(deepTail <= digest.max(), "the estimate must not exceed the maximum"); + assertTrue(deepTail >= sorted[sorted.length - 10], "the deep tail should sit among the largest samples"); + + assertEquals(1.0, digest.cdf(digest.max()), 1e-9); + assertEquals(0.0, digest.cdf(digest.min()), 0.01); + } + + @Test + @DisplayName("a stream made of two repeated values keeps its quantiles sane") + void handlesHeavilyDuplicatedValues() { + TDigest digest = new TDigest(); + for (int i = 0; i < 20_000; i++) { + digest.add(i % 2 == 0 ? 1.0 : 2.0); + } + + assertEquals(1.0, digest.min()); + assertEquals(2.0, digest.max()); + assertTrue(digest.quantile(0.25) <= 1.5, "the lower quarter is made of ones"); + assertTrue(digest.quantile(0.75) >= 1.5, "the upper quarter is made of twos"); + assertTrue(digest.cdf(1.5) >= 0.0 && digest.cdf(1.5) <= 1.0); + } + + @Test + @DisplayName("weights spanning eighteen orders of magnitude still give sane, monotone quantiles") + void survivesWildlyDisparateWeights() { + TDigest digest = new TDigest(50.0); + digest.add(0.0, 1e12); + for (int i = 1; i <= 10_000; i++) { + digest.add(i, 1e-6); + } + + // A point mass this dominant is the regime where the packing rule degenerates into one centroid + // per sample, so it is also the regime that decides whether the centroid array is large enough. + assertEquals(1e12 + 10_000 * 1e-6, digest.totalWeight(), 1.0); + assertEquals(0.0, digest.min()); + assertEquals(10_000.0, digest.max()); + assertEquals(0.0, digest.quantile(0.5), 1e-6, "essentially all the weight sits at zero"); + assertTrue(digest.quantile(0.999999) > 0.0, "the tail must not collapse onto the bulk"); + + double previous = digest.quantile(0.0); + for (int i = 1; i <= 1_000; i++) { + double current = digest.quantile(i / 1_000.0); + assertTrue(current >= previous, "quantiles stopped being monotonic at q=" + i / 1_000.0); + previous = current; + } + } +}