Fix and clean percentile computation

Author: JonasKunz

libs/exponential-histogram/src/main/java/org/elasticsearch/exponentialhistogram/ExpHistoPercentiles.java

@@ -0,0 +1,81 @@
+/*
+ * Copyright Elasticsearch B.V. and/or licensed to Elasticsearch B.V. under one
+ * or more contributor license agreements. Licensed under the "Elastic License
+ * 2.0", the "GNU Affero General Public License v3.0 only", and the "Server Side
+ * Public License v 1"; you may not use this file except in compliance with, at
+ * your election, the "Elastic License 2.0", the "GNU Affero General Public
+ * License v3.0 only", or the "Server Side Public License, v 1".
+ */
+
+package org.elasticsearch.exponentialhistogram;
+
+public class ExponentialHistogramQuantile {
+
+    /**
+     * Provides a quantile for the distribution represented by the given histogram.
+     *
+     * It returns the value of the element at rank {@code max(0, min( n-1, (quantile * (n+1))-1))}, where rank starts at 0.
+     * If that value is fractional, we linearly interpolate based on the fraction the values of the two neighboring ranks.
+     *
+     * @param histo the histogram representing the distribution
+     * @param quantile the quantile to query, in the range [0,1]
+     * @return NaN if the histogram is empty, otherwise the quantile
+     */
+    public static double getQuantile(ExponentialHistogram histo, double quantile) {
+        if (quantile < 0 || quantile > 1) {
+            throw new IllegalArgumentException("quantile must be in range [0, 1]");
+        }
+
+        long zeroCount = histo.zeroBucket().count();
+        long negCount = getTotalCount(histo.negativeBuckets());
+        long posCount = getTotalCount(histo.positiveBuckets());
+
+        long totalCount = zeroCount + negCount + posCount;
+        if (totalCount == 0) {
+            // Can't compute quantile on empty histogram
+            return Double.NaN;
+        }
+
+        double exactRank = Math.max(0, Math.min(totalCount - 1, (totalCount + 1) * quantile - 1));
+        long lowerRank = (long) Math.floor(exactRank);
+        long upperRank = (long) Math.ceil(exactRank);
+        double upperFactor = exactRank - lowerRank;
+
+        return getElementAtRank(histo, lowerRank, negCount, zeroCount) * ( 1 - upperFactor)
+            +getElementAtRank(histo, upperRank, negCount, zeroCount) * upperFactor;
+    }
+
+    private static double getElementAtRank(ExponentialHistogram histo, long rank, long negCount, long zeroCount) {
+        if (rank < negCount) {
+            return -getBucketMidpointForRank(histo.negativeBuckets(), (negCount - 1) - rank);
+        } else if (rank < (negCount + zeroCount)) {
+            return 0.0;
+        } else {
+            return getBucketMidpointForRank(histo.positiveBuckets(), rank - (negCount + zeroCount));
+        }
+    }
+
+    private static double getBucketMidpointForRank(ExponentialHistogram.BucketIterator buckets, long rank) {
+        long seenCount = 0;
+        while (buckets.hasNext()) {
+            seenCount += buckets.peekCount();
+            if (rank < seenCount) {
+                double prev = ExponentialHistogramUtils.getLowerBucketBoundary(buckets.peekIndex(), buckets.scale());
+                double next = ExponentialHistogramUtils.getLowerBucketBoundary(buckets.peekIndex()+1, buckets.scale());
+                double result = ExponentialHistogramUtils.getPointOfLeastRelativeError(buckets.peekIndex(), buckets.scale());
+                return result;
+            }
+            buckets.advance();
+        }
+        throw new IllegalStateException("buckets contain in total less elements than the desired rank");
+    }
+
+    private static long getTotalCount(ExponentialHistogram.BucketIterator buckets) {
+        long count = 0;
+        while (buckets.hasNext()) {
+            count += buckets.peekCount();
+            buckets.advance();
+        }
+        return count;
+    }
+}

libs/exponential-histogram/src/main/java/org/elasticsearch/exponentialhistogram/ExponentialHistogramQuantile.java

@@ -38,9 +38,9 @@ public void testPrintBucketsLinearScale() {
             IntStream.range(-1_000_000, 2_000_000).mapToDouble(Double::valueOf)
         );
 
-        double smallPerc = ExpHistoPercentiles.getPercentile(result, 0.00001);
-        double highPerc = ExpHistoPercentiles.getPercentile(result, 0.9999);
-        double median = ExpHistoPercentiles.getPercentile(result, 0.5);
+        double smallPerc = ExponentialHistogramQuantile.getQuantile(result, 0.00001);
+        double highPerc = ExponentialHistogramQuantile.getQuantile(result, 0.9999);
+        double median = ExponentialHistogramQuantile.getQuantile(result, 0.5);
 
         printMidpoints(result);
     }

libs/exponential-histogram/src/test/java/org/elasticsearch/exponentialhistogram/FixedSizeExponentialHistogramTests.java

@@ -23,66 +23,58 @@
 
 import java.util.Arrays;
 import java.util.Random;
+import java.util.stream.DoubleStream;
+import java.util.stream.IntStream;
 
 import static org.hamcrest.Matchers.closeTo;
 import static org.hamcrest.Matchers.lessThan;
 
-public class PercentileAccuracyTests extends ESTestCase {
+public class QuantileAccuracyTests extends ESTestCase {
 
-    public static final double[] PERCENTILES_TO_TEST = { 0, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99, 1.0 };
+    public static final double[] QUANTILES_TO_TEST = { 0, 0.01, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99, 1.0 };
+
+    public void testBasicSmall() {
+        DoubleStream values = IntStream.range(1,10).mapToDouble(Double::valueOf);
+        testQuantileAccuracy(values.toArray(), 100);
+    }
 
-    /**
-     * Test percentile accuracy with uniform distribution
-     */
     public void testUniformDistribution() {
-        testDistributionPercentileAccuracy(new UniformRealDistribution(new Well19937c(42), 0, 100), 50000, 500);
+        testDistributionQuantileAccuracy(new UniformRealDistribution(new Well19937c(42), 0, 100), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with normal distribution
-     */
     public void testNormalDistribution() {
-        testDistributionPercentileAccuracy(new NormalDistribution(new Well19937c(42), 100, 15), 50000, 500);
+        testDistributionQuantileAccuracy(new NormalDistribution(new Well19937c(42), 100, 15), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with exponential distribution
-     */
     public void testExponentialDistribution() {
-        testDistributionPercentileAccuracy(new ExponentialDistribution(new Well19937c(42), 10), 50000, 500);
+        testDistributionQuantileAccuracy(new ExponentialDistribution(new Well19937c(42), 10), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with log-normal distribution
-     */
     public void testLogNormalDistribution() {
-        testDistributionPercentileAccuracy(new LogNormalDistribution(new Well19937c(42), 0, 1), 50000, 500);
+        testDistributionQuantileAccuracy(new LogNormalDistribution(new Well19937c(42), 0, 1), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with gamma distribution
-     */
     public void testGammaDistribution() {
-        testDistributionPercentileAccuracy(new GammaDistribution(new Well19937c(42), 2, 5), 50000, 500);
+        testDistributionQuantileAccuracy(new GammaDistribution(new Well19937c(42), 2, 5), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with beta distribution
-     */
     public void testBetaDistribution() {
-        testDistributionPercentileAccuracy(new BetaDistribution(new Well19937c(42), 2, 5), 50000, 500);
+        testDistributionQuantileAccuracy(new BetaDistribution(new Well19937c(42), 2, 5), 50000, 500);
     }
 
-    /**
-     * Test percentile accuracy with Weibull distribution
-     */
     public void testWeibullDistribution() {
-        testDistributionPercentileAccuracy(new WeibullDistribution(new Well19937c(42), 2, 5), 50000, 500);
+        testDistributionQuantileAccuracy(new WeibullDistribution(new Well19937c(42), 2, 5), 50000, 500);
+    }
+
+    public void testBigJump() {
+        double[] values = DoubleStream.concat(
+            IntStream.range(0,18).mapToDouble(Double::valueOf),
+            DoubleStream.of(1_000_000.0)
+        ).toArray();
+
+        testQuantileAccuracy(values, 500);
     }
 
-    /**
-     * Test how bucket count affects percentile accuracy
-     */
     public void testBucketCountImpact() {
         RealDistribution distribution = new LogNormalDistribution(new Well19937c(42), 0, 1);
         int sampleSize = 50000;
@@ -91,32 +83,26 @@ public void testBucketCountImpact() {
         // Test with different bucket counts
         int[] bucketCounts = { 10, 50, 100, 200, 500 };
         for (int bucketCount : bucketCounts) {
-            double maxError = testPercentileAccuracy(values, bucketCount);
+            double maxError = testQuantileAccuracy(values, bucketCount);
             logger.info("Bucket count: " + bucketCount + ", Max relative error: " + maxError);
         }
 
         // Verify that more buckets generally means better accuracy
-        double errorWithFewBuckets = testPercentileAccuracy(values, 20);
-        double errorWithManyBuckets = testPercentileAccuracy(values, 200);
+        double errorWithFewBuckets = testQuantileAccuracy(values, 20);
+        double errorWithManyBuckets = testQuantileAccuracy(values, 200);
         assertThat("More buckets should improve accuracy", errorWithManyBuckets, lessThan(errorWithFewBuckets));
     }
 
-    /**
-     * Test percentile accuracy with mixed positive and negative values
-     */
     public void testMixedSignValues() {
         Random random = new Random(42);
         double[] values = new double[10000];
         for (int i = 0; i < values.length; i++) {
             values[i] = (random.nextDouble() * 200) - 100; // Range from -100 to 100
         }
 
-        testPercentileAccuracy(values, 100);
+        testQuantileAccuracy(values, 100);
     }
 
-    /**
-     * Test percentile accuracy with skewed data
-     */
     public void testSkewedData() {
         // Create a highly skewed dataset
         Random random = new Random(42);
@@ -131,12 +117,9 @@ public void testSkewedData() {
             }
         }
 
-        testPercentileAccuracy(values, 100);
+        testQuantileAccuracy(values, 100);
     }
 
-    /**
-     * Test percentile accuracy with data containing zeros
-     */
     public void testDataWithZeros() {
         Random random = new Random(42);
         double[] values = new double[10000];
@@ -149,32 +132,23 @@ public void testDataWithZeros() {
             }
         }
 
-        testPercentileAccuracy(values, 100);
+        testQuantileAccuracy(values, 100);
     }
 
-    /**
-     * Helper method to test percentile accuracy for a given distribution
-     */
-    private void testDistributionPercentileAccuracy(RealDistribution distribution, int sampleSize, int bucketCount) {
+    private void testDistributionQuantileAccuracy(RealDistribution distribution, int sampleSize, int bucketCount) {
         double[] values = generateSamples(distribution, sampleSize);
-        testPercentileAccuracy(values, bucketCount);
+        testQuantileAccuracy(values, bucketCount);
     }
 
-    /**
-     * Helper method to generate samples from a distribution
-     */
-    private double[] generateSamples(RealDistribution distribution, int sampleSize) {
+    private static double[] generateSamples(RealDistribution distribution, int sampleSize) {
         double[] values = new double[sampleSize];
         for (int i = 0; i < sampleSize; i++) {
             values[i] = distribution.sample();
         }
         return values;
     }
 
-    /**
-     * Helper method to test percentile accuracy for a given dataset
-     */
-    private double testPercentileAccuracy(double[] values, int bucketCount) {
+    private double testQuantileAccuracy(double[] values, int bucketCount) {
         // Create histogram
         ExponentialHistogram histogram = ExponentialHistogramGenerator.createFor(bucketCount, Arrays.stream(values));
 
@@ -185,17 +159,17 @@ private double testPercentileAccuracy(double[] values, int bucketCount) {
         double allowedError = getMaximumRelativeError(values, bucketCount);
         double maxError = 0;
 
-        // Compare histogram percentiles with exact percentiles
-        for (double p : PERCENTILES_TO_TEST) {
+        // Compare histogram quantiles with exact quantiles
+        for (double q : QUANTILES_TO_TEST) {
             double exactValue;
-            if (p == 0) {
+            if (q == 0) {
                 exactValue = Arrays.stream(values).min().getAsDouble();
-            } else if (p == 1) {
+            } else if (q == 1) {
                 exactValue = Arrays.stream(values).max().getAsDouble();
             } else {
-                exactValue = exactPercentile.evaluate(p * 100);
+                exactValue = exactPercentile.evaluate(q * 100);
             }
-            double histoValue = ExpHistoPercentiles.getPercentile(histogram, p);
+            double histoValue = ExponentialHistogramQuantile.getQuantile(histogram, q);
 
             // Skip comparison if exact value is zero to avoid division by zero
             if (Math.abs(exactValue) < 1e-10) {
@@ -207,8 +181,8 @@ private double testPercentileAccuracy(double[] values, int bucketCount) {
 
             logger.info(
                 String.format(
-                    "Percentile %.2f: Exact=%.6f, Histogram=%.6f, Relative Error=%.8f, Allowed Relative Error=%.8f",
-                    p,
+                    "Quantile %.2f: Exact=%.6f, Histogram=%.6f, Relative Error=%.8f, Allowed Relative Error=%.8f",
+                    q,
                     exactValue,
                     histoValue,
                     relativeError,
@@ -217,7 +191,7 @@ private double testPercentileAccuracy(double[] values, int bucketCount) {
             );
 
             assertThat(
-                String.format("Percentile %.2f should be accurate within %.6f%% relative error", p, allowedError * 100),
+                String.format("Quantile %.2f should be accurate within %.6f%% relative error", q, allowedError * 100),
                 histoValue,
                 closeTo(exactValue, Math.abs(exactValue * allowedError))
             );
@@ -231,7 +205,7 @@ private double testPercentileAccuracy(double[] values, int bucketCount) {
      * The error depends on the raw values put into the histogram and the number of buckets allowed.
      * This is an implementation of the error bound computation proven by Theorem 3 in the <a href="https://arxiv.org/pdf/2004.08604">UDDSketch paper</a>
      */
-    private double getMaximumRelativeError(double[] values, int bucketCount) {
+    private static double getMaximumRelativeError(double[] values, int bucketCount) {
         double smallestAbsNegative = Double.MAX_VALUE;
         double largestAbsNegative = 0;
         double smallestPositive = Double.MAX_VALUE;