x2 is greater than
+ * x1.
+ */
+ static double weightedAverage(double x1, double w1, double x2, double w2) {
+ if (x1 <= x2) {
+ return weightedAverageSorted(x1, w1, x2, w2);
+ } else {
+ return weightedAverageSorted(x2, w2, x1, w1);
+ }
+ }
+
+ /**
+ * Compute the weighted average between x1 with a weight of
+ * w1 and x2 with a weight of w2.
+ * This expects x1 to be less than or equal to x2
+ * and is guaranteed to return a number between x1 and
+ * x2.
+ */
+ private static double weightedAverageSorted(double x1, double w1, double x2, double w2) {
+ assert x1 <= x2;
+ final double x = (x1 * w1 + x2 * w2) / (w1 + w2);
+ return Math.max(x1, Math.min(x, x2));
+ }
+
+ static double interpolate(double x, double x0, double x1) {
+ return (x - x0) / (x1 - x0);
+ }
+
+ static void encode(ByteBuffer buf, int n) {
+ int k = 0;
+ while (n < 0 || n > 0x7f) {
+ byte b = (byte) (0x80 | (0x7f & n));
+ buf.put(b);
+ n = n >>> 7;
+ k++;
+ if (k >= 6) {
+ throw new IllegalStateException("Size is implausibly large");
+ }
+ }
+ buf.put((byte) n);
+ }
+
+ static int decode(ByteBuffer buf) {
+ int v = buf.get();
+ int z = 0x7f & v;
+ int shift = 7;
+ while ((v & 0x80) != 0) {
+ if (shift > 28) {
+ throw new IllegalStateException("Shift too large in decode");
+ }
+ v = buf.get();
+ z += (v & 0x7f) << shift;
+ shift += 7;
+ }
+ return z;
+ }
+
+ abstract void add(double x, int w, Centroid base);
+
+ /**
+ * Computes an interpolated value of a quantile that is between two centroids.
+ *
+ * Index is the quantile desired multiplied by the total number of samples - 1.
+ *
+ * @param index Denormalized quantile desired
+ * @param previousIndex The denormalized quantile corresponding to the center of the previous centroid.
+ * @param nextIndex The denormalized quantile corresponding to the center of the following centroid.
+ * @param previousMean The mean of the previous centroid.
+ * @param nextMean The mean of the following centroid.
+ * @return The interpolated mean.
+ */
+ static double quantile(double index, double previousIndex, double nextIndex, double previousMean, double nextMean) {
+ final double delta = nextIndex - previousIndex;
+ final double previousWeight = (nextIndex - index) / delta;
+ final double nextWeight = (index - previousIndex) / delta;
+ return previousMean * previousWeight + nextMean * nextWeight;
+ }
+
+ /**
+ * Sets up so that all centroids will record all data assigned to them. For testing only, really.
+ */
+ @Override
+ public TDigest recordAllData() {
+ recordAllData = true;
+ return this;
+ }
+
+ @Override
+ public boolean isRecording() {
+ return recordAllData;
+ }
+
+ /**
+ * Adds a sample to a histogram.
+ *
+ * @param x The value to add.
+ */
+ @Override
+ public void add(double x) {
+ add(x, 1);
+ }
+
+ @Override
+ public void add(TDigest other) {
+ List+ * This can be very fast because the cost of sorting and merging is amortized over several insertion. If + * we keep N centroids total and have the input array is k long, then the amortized cost is something like + *
+ * N/k + log k + *
+ * These costs even out when N/k = log k. Balancing costs is often a good place to start in optimizing an + * algorithm. For different values of compression factor, the following table shows estimated asymptotic + * values of N and suggested values of k: + *
| Compression | N | k |
| 50 | 78 | 25 |
| 100 | 157 | 42 |
| 200 | 314 | 73 |
+ * The virtues of this kind of t-digest implementation include: + *
+ * The current implementation takes the liberty of using ping-pong buffers for implementing the merge resulting
+ * in a substantial memory penalty, but the complexity of an in place merge was not considered as worthwhile
+ * since even with the overhead, the memory cost is less than 40 bytes per centroid which is much less than half
+ * what the AVLTreeDigest uses and no dynamic allocation is required at all.
+ */
+public class MergingDigest extends AbstractTDigest {
+ private int mergeCount = 0;
+
+ private final double publicCompression;
+ private final double compression;
+
+ // points to the first unused centroid
+ private int lastUsedCell;
+
+ // sum_i weight[i] See also unmergedWeight
+ private double totalWeight = 0;
+
+ // number of points that have been added to each merged centroid
+ private final double[] weight;
+ // mean of points added to each merged centroid
+ private final double[] mean;
+
+ // history of all data added to centroids (for testing purposes)
+ private List
+ * The base forms (K_0, K_1, K_2 and K_3) all result in t-digests limited to a number of clusters equal to the
+ * compression factor. The K_2_NO_NORM and K_3_NO_NORM versions result in the cluster count increasing roughly with
+ * log(n).
+ */
+public enum ScaleFunction {
+ /**
+ * Generates uniform cluster sizes. Used for comparison only.
+ */
+ K_0 {
+ @Override
+ public double k(double q, double compression, double n) {
+ return compression * q / 2;
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ return normalizer * q;
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ return 2 * k / compression;
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ return k / normalizer;
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ return 2 / compression;
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ return 1 / normalizer;
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression / 2;
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to sqrt(q*(1-q)). This gives constant relative accuracy if accuracy is
+ * proportional to squared cluster size. It is expected that K_2 and K_3 will give better practical results.
+ */
+ K_1 {
+ @Override
+ public double k(double q, double compression, double n) {
+ return compression * Math.asin(2 * q - 1) / (2 * Math.PI);
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ return normalizer * Math.asin(2 * q - 1);
+ }
+
+
+ @Override
+ public double q(double k, double compression, double n) {
+ return (Math.sin(k * (2 * Math.PI / compression)) + 1) / 2;
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ return (Math.sin(k / normalizer) + 1) / 2;
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ if (q <= 0) {
+ return 0;
+ } else if (q >= 1) {
+ return 0;
+ } else {
+ return 2 * Math.sin(Math.PI / compression) * Math.sqrt(q * (1 - q));
+ }
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ if (q <= 0) {
+ return 0;
+ } else if (q >= 1) {
+ return 0;
+ } else {
+ return 2 * Math.sin(0.5 / normalizer) * Math.sqrt(q * (1 - q));
+ }
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression / (2 * Math.PI);
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to sqrt(q*(1-q)) but avoids computation of asin in the critical path by
+ * using an approximate version.
+ */
+ K_1_FAST {
+ @Override
+ public double k(double q, double compression, double n) {
+ return compression * fastAsin(2 * q - 1) / (2 * Math.PI);
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ return normalizer * fastAsin(2 * q - 1);
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ return (Math.sin(k * (2 * Math.PI / compression)) + 1) / 2;
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ return (Math.sin(k / normalizer) + 1) / 2;
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ if (q <= 0) {
+ return 0;
+ } else if (q >= 1) {
+ return 0;
+ } else {
+ return 2 * Math.sin(Math.PI / compression) * Math.sqrt(q * (1 - q));
+ }
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ if (q <= 0) {
+ return 0;
+ } else if (q >= 1) {
+ return 0;
+ } else {
+ return 2 * Math.sin(0.5 / normalizer) * Math.sqrt(q * (1 - q));
+ }
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression / (2 * Math.PI);
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to q*(1-q). This makes tail error bounds tighter than for K_1. The use of a
+ * normalizing function results in a strictly bounded number of clusters no matter how many samples.
+ */
+ K_2 {
+ @Override
+ public double k(double q, double compression, double n) {
+ if (n <= 1) {
+ if (q <= 0) {
+ return -10;
+ } else if (q >= 1) {
+ return 10;
+ } else {
+ return 0;
+ }
+ }
+ if (q == 0) {
+ return 2 * k(1 / n, compression, n);
+ } else if (q == 1) {
+ return 2 * k((n - 1) / n, compression, n);
+ } else {
+ return compression * Math.log(q / (1 - q)) / Z(compression, n);
+ }
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ if (q < 1e-15) {
+ // this will return something more extreme than q = 1/n
+ return 2 * k(1e-15, normalizer);
+ } else if (q > 1 - 1e-15) {
+ // this will return something more extreme than q = (n-1)/n
+ return 2 * k(1 - 1e-15, normalizer);
+ } else {
+ return Math.log(q / (1 - q)) * normalizer;
+ }
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ double w = Math.exp(k * Z(compression, n) / compression);
+ return w / (1 + w);
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ double w = Math.exp(k / normalizer);
+ return w / (1 + w);
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ return Z(compression, n) * q * (1 - q) / compression;
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ return q * (1 - q) / normalizer;
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression / Z(compression, n);
+ }
+
+ private double Z(double compression, double n) {
+ return 4 * Math.log(n / compression) + 24;
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to min(q, 1-q). This makes tail error bounds tighter than for K_1 or K_2.
+ * The use of a normalizing function results in a strictly bounded number of clusters no matter how many samples.
+ */
+ K_3 {
+ @Override
+ public double k(double q, double compression, double n) {
+ if (q < 0.9 / n) {
+ return 10 * k(1 / n, compression, n);
+ } else if (q > 1 - 0.9 / n) {
+ return 10 * k((n - 1) / n, compression, n);
+ } else {
+ if (q <= 0.5) {
+ return compression * Math.log(2 * q) / Z(compression, n);
+ } else {
+ return -k(1 - q, compression, n);
+ }
+ }
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ if (q < 1e-15) {
+ return 10 * k(1e-15, normalizer);
+ } else if (q > 1 - 1e-15) {
+ return 10 * k(1 - 1e-15, normalizer);
+ } else {
+ if (q <= 0.5) {
+ return Math.log(2 * q) / normalizer;
+ } else {
+ return -k(1 - q, normalizer);
+ }
+ }
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ if (k <= 0) {
+ return Math.exp(k * Z(compression, n) / compression) / 2;
+ } else {
+ return 1 - q(-k, compression, n);
+ }
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ if (k <= 0) {
+ return Math.exp(k / normalizer) / 2;
+ } else {
+ return 1 - q(-k, normalizer);
+ }
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ return Z(compression, n) * Math.min(q, 1 - q) / compression;
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ return Math.min(q, 1 - q) / normalizer;
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression / Z(compression, n);
+ }
+
+ private double Z(double compression, double n) {
+ return 4 * Math.log(n / compression) + 21;
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to q*(1-q). This makes the tail error bounds tighter. This version does not
+ * use a normalizer function and thus the number of clusters increases roughly proportional to log(n). That is good
+ * for accuracy, but bad for size and bad for the statically allocated MergingDigest, but can be useful for
+ * tree-based implementations.
+ */
+ K_2_NO_NORM {
+ @Override
+ public double k(double q, double compression, double n) {
+ if (q == 0) {
+ return 2 * k(1 / n, compression, n);
+ } else if (q == 1) {
+ return 2 * k((n - 1) / n, compression, n);
+ } else {
+ return compression * Math.log(q / (1 - q));
+ }
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ if (q <= 1e-15) {
+ return 2 * k(1e-15, normalizer);
+ } else if (q >= 1 - 1e-15) {
+ return 2 * k(1 - 1e-15, normalizer);
+ } else {
+ return normalizer * Math.log(q / (1 - q));
+ }
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ double w = Math.exp(k / compression);
+ return w / (1 + w);
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ double w = Math.exp(k / normalizer);
+ return w / (1 + w);
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ return q * (1 - q) / compression;
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ return q * (1 - q) / normalizer;
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression;
+ }
+ },
+
+ /**
+ * Generates cluster sizes proportional to min(q, 1-q). This makes the tail error bounds tighter. This version does
+ * not use a normalizer function and thus the number of clusters increases roughly proportional to log(n). That is
+ * good for accuracy, but bad for size and bad for the statically allocated MergingDigest, but can be useful for
+ * tree-based implementations.
+ */
+ K_3_NO_NORM {
+ @Override
+ public double k(double q, double compression, double n) {
+ if (q < 0.9 / n) {
+ return 10 * k(1 / n, compression, n);
+ } else if (q > 1 - 0.9 / n) {
+ return 10 * k((n - 1) / n, compression, n);
+ } else {
+ if (q <= 0.5) {
+ return compression * Math.log(2 * q);
+ } else {
+ return -k(1 - q, compression, n);
+ }
+ }
+ }
+
+ @Override
+ public double k(double q, double normalizer) {
+ if (q <= 1e-15) {
+ return 10 * k(1e-15, normalizer);
+ } else if (q > 1 - 1e-15) {
+ return 10 * k(1 - 1e-15, normalizer);
+ } else {
+ if (q <= 0.5) {
+ return normalizer * Math.log(2 * q);
+ } else {
+ return -k(1 - q, normalizer);
+ }
+ }
+ }
+
+ @Override
+ public double q(double k, double compression, double n) {
+ if (k <= 0) {
+ return Math.exp(k / compression) / 2;
+ } else {
+ return 1 - q(-k, compression, n);
+ }
+ }
+
+ @Override
+ public double q(double k, double normalizer) {
+ if (k <= 0) {
+ return Math.exp(k / normalizer) / 2;
+ } else {
+ return 1 - q(-k, normalizer);
+ }
+ }
+
+ @Override
+ public double max(double q, double compression, double n) {
+ return Math.min(q, 1 - q) / compression;
+ }
+
+ @Override
+ public double max(double q, double normalizer) {
+ return Math.min(q, 1 - q) / normalizer;
+ }
+
+ @Override
+ public double normalizer(double compression, double n) {
+ return compression;
+ }
+ }; // max weight is min(q,1-q), should improve tail accuracy even more
+
+ /**
+ * Converts a quantile to the k-scale. The total number of points is also provided so that a normalizing function
+ * can be computed if necessary.
+ *
+ * @param q The quantile
+ * @param compression Also known as delta in literature on the t-digest
+ * @param n The total number of samples
+ * @return The corresponding value of k
+ */
+ abstract public double k(double q, double compression, double n);
+
+ /**
+ * Converts a quantile to the k-scale. The normalizer value depends on compression and (possibly) number of points
+ * in the digest. #normalizer(double, double)
+ *
+ * @param q The quantile
+ * @param normalizer The normalizer value which depends on compression and (possibly) number of points in the
+ * digest.
+ * @return The corresponding value of k
+ */
+ abstract public double k(double q, double normalizer);
+
+ /**
+ * Computes q as a function of k. This is often faster than finding k as a function of q for some scales.
+ *
+ * @param k The index value to convert into q scale.
+ * @param compression The compression factor (often written as δ)
+ * @param n The number of samples already in the digest.
+ * @return The value of q that corresponds to k
+ */
+ abstract public double q(double k, double compression, double n);
+
+ /**
+ * Computes q as a function of k. This is often faster than finding k as a function of q for some scales.
+ *
+ * @param k The index value to convert into q scale.
+ * @param normalizer The normalizer value which depends on compression and (possibly) number of points in the
+ * digest.
+ * @return The value of q that corresponds to k
+ */
+ abstract public double q(double k, double normalizer);
+
+ /**
+ * Computes the maximum relative size a cluster can have at quantile q. Note that exactly where within the range
+ * spanned by a cluster that q should be isn't clear. That means that this function usually has to be taken at
+ * multiple points and the smallest value used.
+ *
+ * Note that this is the relative size of a cluster. To get the max number of samples in the cluster, multiply this
+ * value times the total number of samples in the digest.
+ *
+ * @param q The quantile
+ * @param compression The compression factor, typically delta in the literature
+ * @param n The number of samples seen so far in the digest
+ * @return The maximum number of samples that can be in the cluster
+ */
+ abstract public double max(double q, double compression, double n);
+
+ /**
+ * Computes the maximum relative size a cluster can have at quantile q. Note that exactly where within the range
+ * spanned by a cluster that q should be isn't clear. That means that this function usually has to be taken at
+ * multiple points and the smallest value used.
+ *
+ * Note that this is the relative size of a cluster. To get the max number of samples in the cluster, multiply this
+ * value times the total number of samples in the digest.
+ *
+ * @param q The quantile
+ * @param normalizer The normalizer value which depends on compression and (possibly) number of points in the
+ * digest.
+ * @return The maximum number of samples that can be in the cluster
+ */
+ abstract public double max(double q, double normalizer);
+
+ /**
+ * Computes the normalizer given compression and number of points.
+ */
+ abstract public double normalizer(double compression, double n);
+
+ /**
+ * Approximates asin to within about 1e-6. This approximation works by breaking the range from 0 to 1 into 5 regions
+ * for all but the region nearest 1, rational polynomial models get us a very good approximation of asin and by
+ * interpolating as we move from region to region, we can guarantee continuity and we happen to get monotonicity as
+ * well. for the values near 1, we just use Math.asin as our region "approximation".
+ *
+ * @param x sin(theta)
+ * @return theta
+ */
+ static double fastAsin(double x) {
+ if (x < 0) {
+ return -fastAsin(-x);
+ } else if (x > 1) {
+ return Double.NaN;
+ } else {
+ // Cutoffs for models. Note that the ranges overlap. In the
+ // overlap we do linear interpolation to guarantee the overall
+ // result is "nice"
+ double c0High = 0.1;
+ double c1High = 0.55;
+ double c2Low = 0.5;
+ double c2High = 0.8;
+ double c3Low = 0.75;
+ double c3High = 0.9;
+ double c4Low = 0.87;
+ if (x > c3High) {
+ return Math.asin(x);
+ } else {
+ // the models
+ double[] m0 = {0.2955302411, 1.2221903614, 0.1488583743, 0.2422015816, -0.3688700895, 0.0733398445};
+ double[] m1 = {-0.0430991920, 0.9594035750, -0.0362312299, 0.1204623351, 0.0457029620, -0.0026025285};
+ double[] m2 = {-0.034873933724, 1.054796752703, -0.194127063385, 0.283963735636, 0.023800124916, -0.000872727381};
+ double[] m3 = {-0.37588391875, 2.61991859025, -2.48835406886, 1.48605387425, 0.00857627492, -0.00015802871};
+
+ // the parameters for all of the models
+ double[] vars = {1, x, x * x, x * x * x, 1 / (1 - x), 1 / (1 - x) / (1 - x)};
+
+ // raw grist for interpolation coefficients
+ double x0 = bound((c0High - x) / c0High);
+ double x1 = bound((c1High - x) / (c1High - c2Low));
+ double x2 = bound((c2High - x) / (c2High - c3Low));
+ double x3 = bound((c3High - x) / (c3High - c4Low));
+
+ // interpolation coefficients
+ //noinspection UnnecessaryLocalVariable
+ double mix0 = x0;
+ double mix1 = (1 - x0) * x1;
+ double mix2 = (1 - x1) * x2;
+ double mix3 = (1 - x2) * x3;
+ double mix4 = 1 - x3;
+
+ // now mix all the results together, avoiding extra evaluations
+ double r = 0;
+ if (mix0 > 0) {
+ r += mix0 * eval(m0, vars);
+ }
+ if (mix1 > 0) {
+ r += mix1 * eval(m1, vars);
+ }
+ if (mix2 > 0) {
+ r += mix2 * eval(m2, vars);
+ }
+ if (mix3 > 0) {
+ r += mix3 * eval(m3, vars);
+ }
+ if (mix4 > 0) {
+ // model 4 is just the real deal
+ r += mix4 * Math.asin(x);
+ }
+ return r;
+ }
+ }
+ }
+
+ private static double eval(double[] model, double[] vars) {
+ double r = 0;
+ for (int i = 0; i < model.length; i++) {
+ r += model[i] * vars[i];
+ }
+ return r;
+ }
+
+ private static double bound(double v) {
+ if (v <= 0) {
+ return 0;
+ } else if (v >= 1) {
+ return 1;
+ } else {
+ return v;
+ }
+ }
+}
\ No newline at end of file
diff --git a/presto-main/src/main/java/com/facebook/presto/operator/scalar/Sort.java b/presto-main/src/main/java/com/facebook/presto/operator/scalar/Sort.java
new file mode 100644
index 0000000000000..602c98430a55c
--- /dev/null
+++ b/presto-main/src/main/java/com/facebook/presto/operator/scalar/Sort.java
@@ -0,0 +1,482 @@
+/*
+ * Licensed to Ted Dunning under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package com.facebook.presto.operator.scalar;
+
+import java.util.Random;
+
+/**
+ * Static sorting methods
+ */
+public class Sort {
+ private static final Random prng = new Random(); // for choosing pivots during quicksort
+ /**
+ * Quick sort using an index array. On return,
+ * values[order[i]] is in order as i goes 0..values.length
+ *
+ * @param order Indexes into values
+ * @param values The values to sort.
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void sort(int[] order, double[] values) {
+ sort(order, values, 0, values.length);
+ }
+
+ /**
+ * Quick sort using an index array. On return,
+ * values[order[i]] is in order as i goes 0..n
+ *
+ * @param order Indexes into values
+ * @param values The values to sort.
+ * @param n The number of values to sort
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void sort(int[] order, double[] values, int n) {
+ sort(order, values, 0, n);
+ }
+
+ /**
+ * Quick sort using an index array. On return,
+ * values[order[i]] is in order as i goes start..n
+ *
+ * @param order Indexes into values
+ * @param values The values to sort.
+ * @param start The first element to sort
+ * @param n The number of values to sort
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void sort(int[] order, double[] values, int start, int n) {
+ for (int i = start; i < start + n; i++) {
+ order[i] = i;
+ }
+ quickSort(order, values, start, start + n, 64);
+ insertionSort(order, values, start, start + n, 64);
+ }
+
+ /**
+ * Standard quick sort except that sorting is done on an index array rather than the values themselves
+ *
+ * @param order The pre-allocated index array
+ * @param values The values to sort
+ * @param start The beginning of the values to sort
+ * @param end The value after the last value to sort
+ * @param limit The minimum size to recurse down to.
+ */
+ private static void quickSort(int[] order, double[] values, int start, int end, int limit) {
+ // the while loop implements tail-recursion to avoid excessive stack calls on nasty cases
+ while (end - start > limit) {
+
+ // pivot by a random element
+ int pivotIndex = start + prng.nextInt(end - start);
+ double pivotValue = values[order[pivotIndex]];
+
+ // move pivot to beginning of array
+ swap(order, start, pivotIndex);
+
+ // we use a three way partition because many duplicate values is an important case
+
+ int low = start + 1; // low points to first value not known to be equal to pivotValue
+ int high = end; // high points to first value > pivotValue
+ int i = low; // i scans the array
+ while (i < high) {
+ // invariant: values[order[k]] == pivotValue for k in [0..low)
+ // invariant: values[order[k]] < pivotValue for k in [low..i)
+ // invariant: values[order[k]] > pivotValue for k in [high..end)
+ // in-loop: i < high
+ // in-loop: low < high
+ // in-loop: i >= low
+ double vi = values[order[i]];
+ if (vi == pivotValue) {
+ if (low != i) {
+ swap(order, low, i);
+ } else {
+ i++;
+ }
+ low++;
+ } else if (vi > pivotValue) {
+ high--;
+ swap(order, i, high);
+ } else {
+ // vi < pivotValue
+ i++;
+ }
+ }
+ // invariant: values[order[k]] == pivotValue for k in [0..low)
+ // invariant: values[order[k]] < pivotValue for k in [low..i)
+ // invariant: values[order[k]] > pivotValue for k in [high..end)
+ // assert i == high || low == high therefore, we are done with partition
+
+ // at this point, i==high, from [start,low) are == pivot, [low,high) are < and [high,end) are >
+ // we have to move the values equal to the pivot into the middle. To do this, we swap pivot
+ // values into the top end of the [low,high) range stopping when we run out of destinations
+ // or when we run out of values to copy
+ int from = start;
+ int to = high - 1;
+ for (i = 0; from < low && to >= low; i++) {
+ swap(order, from++, to--);
+ }
+ if (from == low) {
+ // ran out of things to copy. This means that the the last destination is the boundary
+ low = to + 1;
+ } else {
+ // ran out of places to copy to. This means that there are uncopied pivots and the
+ // boundary is at the beginning of those
+ low = from;
+ }
+
+// checkPartition(order, values, pivotValue, start, low, high, end);
+
+ // now recurse, but arrange it so we handle the longer limit by tail recursion
+ if (low - start < end - high) {
+ quickSort(order, values, start, low, limit);
+
+ // this is really a way to do
+ // quickSort(order, values, high, end, limit);
+ start = high;
+ } else {
+ quickSort(order, values, high, end, limit);
+ // this is really a way to do
+ // quickSort(order, values, start, low, limit);
+ end = low;
+ }
+ }
+ }
+
+ /**
+ * Quick sort in place of several paired arrays. On return,
+ * keys[...] is in order and the values[] arrays will be
+ * reordered as well in the same way.
+ *
+ * @param key Values to sort on
+ * @param values The auxilliary values to sort.
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void sort(double[] key, double[] ... values) {
+ sort(key, 0, key.length, values);
+ }
+
+ /**
+ * Quick sort using an index array. On return,
+ * values[order[i]] is in order as i goes start..n
+ * @param key Values to sort on
+ * @param start The first element to sort
+ * @param n The number of values to sort
+ * @param values The auxilliary values to sort.
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void sort(double[] key, int start, int n, double[]... values) {
+ quickSort(key, values, start, start + n, 8);
+ insertionSort(key, values, start, start + n, 8);
+ }
+
+ /**
+ * Standard quick sort except that sorting rearranges parallel arrays
+ *
+ * @param key Values to sort on
+ * @param values The auxilliary values to sort.
+ * @param start The beginning of the values to sort
+ * @param end The value after the last value to sort
+ * @param limit The minimum size to recurse down to.
+ */
+ private static void quickSort(double[] key, double[][] values, int start, int end, int limit) {
+ // the while loop implements tail-recursion to avoid excessive stack calls on nasty cases
+ while (end - start > limit) {
+
+ // median of three values for the pivot
+ int a = start;
+ int b = (start + end) / 2;
+ int c = end - 1;
+
+ int pivotIndex;
+ double pivotValue;
+ double va = key[a];
+ double vb = key[b];
+ double vc = key[c];
+ //noinspection Duplicates
+ if (va > vb) {
+ if (vc > va) {
+ // vc > va > vb
+ pivotIndex = a;
+ pivotValue = va;
+ } else {
+ // va > vb, va >= vc
+ if (vc < vb) {
+ // va > vb > vc
+ pivotIndex = b;
+ pivotValue = vb;
+ } else {
+ // va >= vc >= vb
+ pivotIndex = c;
+ pivotValue = vc;
+ }
+ }
+ } else {
+ // vb >= va
+ if (vc > vb) {
+ // vc > vb >= va
+ pivotIndex = b;
+ pivotValue = vb;
+ } else {
+ // vb >= va, vb >= vc
+ if (vc < va) {
+ // vb >= va > vc
+ pivotIndex = a;
+ pivotValue = va;
+ } else {
+ // vb >= vc >= va
+ pivotIndex = c;
+ pivotValue = vc;
+ }
+ }
+ }
+
+ // move pivot to beginning of array
+ swap(start, pivotIndex, key, values);
+
+ // we use a three way partition because many duplicate values is an important case
+
+ int low = start + 1; // low points to first value not known to be equal to pivotValue
+ int high = end; // high points to first value > pivotValue
+ int i = low; // i scans the array
+ while (i < high) {
+ // invariant: values[order[k]] == pivotValue for k in [0..low)
+ // invariant: values[order[k]] < pivotValue for k in [low..i)
+ // invariant: values[order[k]] > pivotValue for k in [high..end)
+ // in-loop: i < high
+ // in-loop: low < high
+ // in-loop: i >= low
+ double vi = key[i];
+ if (vi == pivotValue) {
+ if (low != i) {
+ swap(low, i, key, values);
+ } else {
+ i++;
+ }
+ low++;
+ } else if (vi > pivotValue) {
+ high--;
+ swap(i, high, key, values);
+ } else {
+ // vi < pivotValue
+ i++;
+ }
+ }
+ // invariant: values[order[k]] == pivotValue for k in [0..low)
+ // invariant: values[order[k]] < pivotValue for k in [low..i)
+ // invariant: values[order[k]] > pivotValue for k in [high..end)
+ // assert i == high || low == high therefore, we are done with partition
+
+ // at this point, i==high, from [start,low) are == pivot, [low,high) are < and [high,end) are >
+ // we have to move the values equal to the pivot into the middle. To do this, we swap pivot
+ // values into the top end of the [low,high) range stopping when we run out of destinations
+ // or when we run out of values to copy
+ int from = start;
+ int to = high - 1;
+ for (i = 0; from < low && to >= low; i++) {
+ swap(from++, to--, key, values);
+ }
+ if (from == low) {
+ // ran out of things to copy. This means that the the last destination is the boundary
+ low = to + 1;
+ } else {
+ // ran out of places to copy to. This means that there are uncopied pivots and the
+ // boundary is at the beginning of those
+ low = from;
+ }
+
+// checkPartition(order, values, pivotValue, start, low, high, end);
+
+ // now recurse, but arrange it so we handle the longer limit by tail recursion
+ if (low - start < end - high) {
+ quickSort(key, values, start, low, limit);
+
+ // this is really a way to do
+ // quickSort(order, values, high, end, limit);
+ start = high;
+ } else {
+ quickSort(key, values, high, end, limit);
+ // this is really a way to do
+ // quickSort(order, values, start, low, limit);
+ end = low;
+ }
+ }
+ }
+
+
+ /**
+ * Limited range insertion sort. We assume that no element has to move more than limit steps
+ * because quick sort has done its thing. This version works on parallel arrays of keys and values.
+ *
+ * @param key The array of keys
+ * @param values The values we are sorting
+ * @param start The starting point of the sort
+ * @param end The ending point of the sort
+ * @param limit The largest amount of disorder
+ */
+ @SuppressWarnings("SameParameterValue")
+ private static void insertionSort(double[] key, double[][] values, int start, int end, int limit) {
+ // loop invariant: all values start ... i-1 are ordered
+ for (int i = start + 1; i < end; i++) {
+ double v = key[i];
+ int m = Math.max(i - limit, start);
+ for (int j = i; j >= m; j--) {
+ if (j == m || key[j - 1] <= v) {
+ if (j < i) {
+ System.arraycopy(key, j, key, j + 1, i - j);
+ key[j] = v;
+ for (double[] value : values) {
+ double tmp = value[i];
+ System.arraycopy(value, j, value, j + 1, i - j);
+ value[j] = tmp;
+ }
+ }
+ break;
+ }
+ }
+ }
+ }
+
+ private static void swap(int[] order, int i, int j) {
+ int t = order[i];
+ order[i] = order[j];
+ order[j] = t;
+ }
+
+ private static void swap(int i, int j, double[] key, double[]...values) {
+ double t = key[i];
+ key[i] = key[j];
+ key[j] = t;
+
+ for (int k = 0; k < values.length; k++) {
+ t = values[k][i];
+ values[k][i] = values[k][j];
+ values[k][j] = t;
+ }
+ }
+
+ /**
+ * Check that a partition step was done correctly. For debugging and testing.
+ *
+ * @param order The array of indexes representing a permutation of the keys.
+ * @param values The keys to sort.
+ * @param pivotValue The value that splits the data
+ * @param start The beginning of the data of interest.
+ * @param low Values from start (inclusive) to low (exclusive) are < pivotValue.
+ * @param high Values from low to high are equal to the pivot.
+ * @param end Values from high to end are above the pivot.
+ */
+ @SuppressWarnings("UnusedDeclaration")
+ public static void checkPartition(int[] order, double[] values, double pivotValue, int start, int low, int high, int end) {
+ if (order.length != values.length) {
+ throw new IllegalArgumentException("Arguments must be same size");
+ }
+
+ if (!(start >= 0 && low >= start && high >= low && end >= high)) {
+ throw new IllegalArgumentException(String.format("Invalid indices %d, %d, %d, %d", start, low, high, end));
+ }
+
+ for (int i = 0; i < low; i++) {
+ double v = values[order[i]];
+ if (v >= pivotValue) {
+ throw new IllegalArgumentException(String.format("Value greater than pivot at %d", i));
+ }
+ }
+
+ for (int i = low; i < high; i++) {
+ if (values[order[i]] != pivotValue) {
+ throw new IllegalArgumentException(String.format("Non-pivot at %d", i));
+ }
+ }
+
+ for (int i = high; i < end; i++) {
+ double v = values[order[i]];
+ if (v <= pivotValue) {
+ throw new IllegalArgumentException(String.format("Value less than pivot at %d", i));
+ }
+ }
+ }
+
+ /**
+ * Limited range insertion sort. We assume that no element has to move more than limit steps
+ * because quick sort has done its thing.
+ *
+ * @param order The permutation index
+ * @param values The values we are sorting
+ * @param start Where to start the sort
+ * @param n How many elements to sort
+ * @param limit The largest amount of disorder
+ */
+ @SuppressWarnings("SameParameterValue")
+ private static void insertionSort(int[] order, double[] values, int start, int n, int limit) {
+ for (int i = start + 1; i < n; i++) {
+ int t = order[i];
+ double v = values[order[i]];
+ int m = Math.max(i - limit, start);
+ for (int j = i; j >= m; j--) {
+ if (j == 0 || values[order[j - 1]] <= v) {
+ if (j < i) {
+ System.arraycopy(order, j, order, j + 1, i - j);
+ order[j] = t;
+ }
+ break;
+ }
+ }
+ }
+ }
+
+ /**
+ * Reverses an array in-place.
+ *
+ * @param order The array to reverse
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void reverse(int[] order) {
+ reverse(order, 0, order.length);
+ }
+
+ /**
+ * Reverses part of an array. See {@link #reverse(int[])}
+ *
+ * @param order The array containing the data to reverse.
+ * @param offset Where to start reversing.
+ * @param length How many elements to reverse
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static void reverse(int[] order, int offset, int length) {
+ for (int i = 0; i < length / 2; i++) {
+ int t = order[offset + i];
+ order[offset + i] = order[offset + length - i - 1];
+ order[offset + length - i - 1] = t;
+ }
+ }
+
+ /**
+ * Reverses part of an array. See {@link #reverse(int[])}
+ *
+ * @param order The array containing the data to reverse.
+ * @param offset Where to start reversing.
+ * @param length How many elements to reverse
+ */
+ @SuppressWarnings({"WeakerAccess", "SameParameterValue"})
+ public static void reverse(double[] order, int offset, int length) {
+ for (int i = 0; i < length / 2; i++) {
+ double t = order[offset + i];
+ order[offset + i] = order[offset + length - i - 1];
+ order[offset + length - i - 1] = t;
+ }
+ }
+}
\ No newline at end of file
diff --git a/presto-main/src/main/java/com/facebook/presto/operator/scalar/TDigest.java b/presto-main/src/main/java/com/facebook/presto/operator/scalar/TDigest.java
new file mode 100644
index 0000000000000..ec41b52329a5e
--- /dev/null
+++ b/presto-main/src/main/java/com/facebook/presto/operator/scalar/TDigest.java
@@ -0,0 +1,226 @@
+/*
+ * Licensed to Ted Dunning under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package com.facebook.presto.operator.scalar;
+
+import java.io.Serializable;
+import java.nio.ByteBuffer;
+import java.util.Collection;
+import java.util.List;
+
+/**
+ * Adaptive histogram based on something like streaming k-means crossed with Q-digest.
+ *
+ * The special characteristics of this algorithm are:
+ *
+ * - smaller summaries than Q-digest
+ *
+ * - works on doubles as well as integers.
+ *
+ * - provides part per million accuracy for extreme quantiles and typically <1000 ppm accuracy for middle quantiles
+ *
+ * - fast
+ *
+ * - simple
+ *
+ * - test coverage roughly at 90%
+ *
+ * - easy to adapt for use with map-reduce
+ */
+public abstract class TDigest implements Serializable {
+ protected ScaleFunction scale = ScaleFunction.K_2;
+ double min = Double.POSITIVE_INFINITY;
+ double max = Double.NEGATIVE_INFINITY;
+
+ /**
+ * Creates an {@link MergingDigest}. This is generally the best known implementation right now.
+ *
+ * @param compression The compression parameter. 100 is a common value for normal uses. 1000 is extremely large.
+ * The number of centroids retained will be a smallish (usually less than 10) multiple of this number.
+ * @return the MergingDigest
+ */
+ @SuppressWarnings("WeakerAccess")
+ public static TDigest createMergingDigest(double compression) {
+ return new MergingDigest(compression);
+ }
+
+ /**
+ * Creates a TDigest of whichever type is the currently recommended type. MergingDigest is generally the best
+ * known implementation right now.
+ *
+ * @param compression The compression parameter. 100 is a common value for normal uses. 1000 is extremely large.
+ * The number of centroids retained will be a smallish (usually less than 10) multiple of this number.
+ * @return the TDigest
+ */
+ @SuppressWarnings({"unused", "WeakerAccess", "SameParameterValue"})
+ public static TDigest createDigest(double compression) {
+ return createMergingDigest(compression);
+ }
+
+ /**
+ * Adds a sample to a histogram.
+ *
+ * @param x The value to add.
+ * @param w The weight of this point.
+ */
+ public abstract void add(double x, int w);
+
+ final void checkValue(double x) {
+ if (Double.isNaN(x)) {
+ throw new IllegalArgumentException("Cannot add NaN");
+ }
+ }
+
+ public abstract void add(List others);
+
+ /**
+ * Re-examines a t-digest to determine whether some centroids are redundant. If your data are
+ * perversely ordered, this may be a good idea. Even if not, this may save 20% or so in space.
+ *
+ * The cost is roughly the same as adding as many data points as there are centroids. This
+ * is typically < 10 * compression, but could be as high as 100 * compression.
+ *
+ * This is a destructive operation that is not thread-safe.
+ */
+ public abstract void compress();
+
+ /**
+ * Returns the number of points that have been added to this TDigest.
+ *
+ * @return The sum of the weights on all centroids.
+ */
+ public abstract long size();
+
+ /**
+ * Returns the fraction of all points added which are ≤ x.
+ *
+ * @param x The cutoff for the cdf.
+ * @return The fraction of all data which is less or equal to x.
+ */
+ public abstract double cdf(double x);
+
+ /**
+ * Returns an estimate of the cutoff such that a specified fraction of the data
+ * added to this TDigest would be less than or equal to the cutoff.
+ *
+ * @param q The desired fraction
+ * @return The value x such that cdf(x) == q
+ */
+ public abstract double quantile(double q);
+
+ /**
+ * A {@link Collection} that lets you go through the centroids in ascending order by mean. Centroids
+ * returned will not be re-used, but may or may not share storage with this TDigest.
+ *
+ * @return The centroids in the form of a Collection.
+ */
+ public abstract Collection
+ * If the value in the digest needs to be converted before the mean is calculated,
+ * valueFunction can be used. Else, pass in the identity function.
+ */
+ public Double getTruncatedMean(double lowerQuantile, double upperQuantile, TruncMeanValueFunction valueFunction)
+ {
+ if (weightedCount == 0 || lowerQuantile >= upperQuantile) {
+ return null;
+ }
+
+ AtomicDouble meanResult = new AtomicDouble();
+ double lowerRank = lowerQuantile * weightedCount;
+ double upperRank = upperQuantile * weightedCount;
+
+ postOrderTraversal(root, new Callback()
+ {
+ private double sum;
+ private double count;
+ private double mean;
+
+ public boolean process(int node)
+ {
+ double nodeCount = counts[node];
+ if (nodeCount == 0) {
+ return true;
+ }
+
+ sum += nodeCount;
+
+ double amountOverLower = Math.max(sum - lowerRank, 0);
+ double amountOverUpper = Math.max(sum - upperRank, 0);
+ // Constrain node count to the rank of the lower and upper bound quantiles
+ nodeCount = Math.max(0, Math.min(Math.min(nodeCount, amountOverLower), nodeCount - amountOverUpper));
+
+ if (amountOverLower > 0) {
+ double value = valueFunction.value(Math.min(upperBound(node), max));
+ mean = (mean * count + nodeCount * value) / (count + nodeCount);
+ count += nodeCount;
+ }
+
+ if (amountOverUpper > 0 || sum == weightedCount) {
+ meanResult.set(mean);
+ return false;
+ }
+ return true;
+ }
+ });
+
+ return meanResult.get();
+ }
+
+ public long estimatedInMemorySizeInBytes()
+ {
+ return ((long) QUANTILE_DIGEST_SIZE + SizeOf.sizeOf(this.counts) + SizeOf.sizeOf(this.levels) + SizeOf.sizeOf(this.values) + SizeOf.sizeOf(this.lefts) + SizeOf.sizeOf(this.rights));
+ }
+
+ public int estimatedSerializedSizeInBytes()
+ {
+ int nodeSize = 17;
+ return 45 + this.getNodeCount() * nodeSize;
+ }
+
+ public Slice serialize()
+ {
+ this.compress();
+ SliceOutput output = new DynamicSliceOutput(this.estimatedSerializedSizeInBytes());
+ output.writeByte(0);
+ output.writeDouble(this.maxError);
+ output.writeDouble(this.alpha);
+ output.writeLong(this.landmarkInSeconds);
+ output.writeLong(this.min);
+ output.writeLong(this.max);
+ output.writeInt(this.getNodeCount());
+ final int[] nodes = new int[this.getNodeCount()];
+ this.postOrderTraversal(this.root, new QuantileDigest.Callback()
+ {
+ int index;
+
+ public boolean process(int node)
+ {
+ nodes[this.index++] = node;
+ return true;
+ }
+ });
+ int[] var3 = nodes;
+ int var4 = nodes.length;
+
+ for (int var5 = 0; var5 < var4; ++var5) {
+ int node = var3[var5];
+ byte nodeStructure = (byte) (Math.max(this.levels[node] - 1, 0) << 2);
+ if (this.lefts[node] != -1) {
+ nodeStructure = (byte) (nodeStructure | 1);
+ }
+
+ if (this.rights[node] != -1) {
+ nodeStructure = (byte) (nodeStructure | 2);
+ }
+
+ output.writeByte(nodeStructure);
+ output.writeDouble(this.counts[node]);
+ output.writeLong(this.values[node]);
+ }
+
+ return output.slice();
+ }
+
+ @VisibleForTesting
+ int getNodeCount()
+ {
+ return this.nextNode - this.freeCount;
+ }
+
+ @VisibleForTesting
+ void compress()
+ {
+ double bound = Math.floor(this.weightedCount / (double) this.calculateCompressionFactor());
+ this.postOrderTraversal(this.root, (node) -> {
+ int left = this.lefts[node];
+ int right = this.rights[node];
+ if (left == -1 && right == -1) {
+ return true;
+ }
+ else {
+ double leftCount = left == -1 ? 0.0D : this.counts[left];
+ double rightCount = right == -1 ? 0.0D : this.counts[right];
+ boolean shouldCompress = this.counts[node] + leftCount + rightCount < bound;
+ double[] var10000;
+ if (left != -1 && (shouldCompress || leftCount < 1.0E-5D)) {
+ this.lefts[node] = this.tryRemove(left);
+ var10000 = this.counts;
+ var10000[node] += leftCount;
+ }
+
+ if (right != -1 && (shouldCompress || rightCount < 1.0E-5D)) {
+ this.rights[node] = this.tryRemove(right);
+ var10000 = this.counts;
+ var10000[node] += rightCount;
+ }
+
+ return true;
+ }
+ });
+ if (this.root != -1 && this.counts[this.root] < 1.0E-5D) {
+ this.root = this.tryRemove(this.root);
+ }
+ }
+
+ private double weight(long timestamp)
+ {
+ return Math.exp(this.alpha * (double) (timestamp - this.landmarkInSeconds));
+ }
+
+ private void rescale(long newLandmarkInSeconds)
+ {
+ double factor = Math.exp(-this.alpha * (double) (newLandmarkInSeconds - this.landmarkInSeconds));
+ this.weightedCount *= factor;
+
+ for (int i = 0; i < this.nextNode; ++i) {
+ double[] var10000 = this.counts;
+ var10000[i] *= factor;
+ }
+
+ this.landmarkInSeconds = newLandmarkInSeconds;
+ }
+
+ private int calculateCompressionFactor()
+ {
+ return this.root == -1 ? 1 : Math.max((int) ((double) (this.levels[this.root] + 1) / this.maxError), 1);
+ }
+
+ private void insert(long value, double count)
+ {
+ if (count >= 1.0E-5D) {
+ long lastBranch = 0L;
+ int parent = -1;
+ int current = this.root;
+
+ while (current != -1) {
+ long currentValue = this.values[current];
+ byte currentLevel = this.levels[current];
+ if (!inSameSubtree(value, currentValue, currentLevel)) {
+ this.setChild(parent, lastBranch, this.makeSiblings(current, this.createLeaf(value, count)));
+ return;
+ }
+
+ if (currentLevel == 0 && currentValue == value) {
+ double[] var10000 = this.counts;
+ var10000[current] += count;
+ this.weightedCount += count;
+ return;
+ }
+
+ long branch = value & this.getBranchMask(currentLevel);
+ parent = current;
+ lastBranch = branch;
+ if (branch == 0L) {
+ current = this.lefts[current];
+ }
+ else {
+ current = this.rights[current];
+ }
+ }
+
+ this.setChild(parent, lastBranch, this.createLeaf(value, count));
+ }
+ }
+
+ private void setChild(int parent, long branch, int child)
+ {
+ if (parent == -1) {
+ this.root = child;
+ }
+ else if (branch == 0L) {
+ this.lefts[parent] = child;
+ }
+ else {
+ this.rights[parent] = child;
+ }
+ }
+
+ private int makeSiblings(int first, int second)
+ {
+ long firstValue = this.values[first];
+ long secondValue = this.values[second];
+ int parentLevel = 64 - Long.numberOfLeadingZeros(firstValue ^ secondValue);
+ int parent = this.createNode(firstValue, parentLevel, 0.0D);
+ long branch = firstValue & this.getBranchMask(this.levels[parent]);
+ if (branch == 0L) {
+ this.lefts[parent] = first;
+ this.rights[parent] = second;
+ }
+ else {
+ this.lefts[parent] = second;
+ this.rights[parent] = first;
+ }
+
+ return parent;
+ }
+
+ private int createLeaf(long value, double count)
+ {
+ return this.createNode(value, 0, count);
+ }
+
+ private int createNode(long value, int level, double count)
+ {
+ int node = this.popFree();
+ if (node == -1) {
+ if (this.nextNode == this.counts.length) {
+ int newSize = this.counts.length + Math.min(this.counts.length, this.calculateCompressionFactor() / 5 + 1);
+ this.counts = Arrays.copyOf(this.counts, newSize);
+ this.levels = Arrays.copyOf(this.levels, newSize);
+ this.values = Arrays.copyOf(this.values, newSize);
+ this.lefts = Arrays.copyOf(this.lefts, newSize);
+ this.rights = Arrays.copyOf(this.rights, newSize);
+ }
+
+ node = this.nextNode++;
+ }
+
+ this.weightedCount += count;
+ this.values[node] = value;
+ this.levels[node] = (byte) level;
+ this.counts[node] = count;
+ this.lefts[node] = -1;
+ this.rights[node] = -1;
+ return node;
+ }
+
+ private int merge(int node, QuantileDigest other, int otherNode)
+ {
+ if (otherNode == -1) {
+ return node;
+ }
+ else if (node == -1) {
+ return this.copyRecursive(other, otherNode);
+ }
+ else if (!inSameSubtree(this.values[node], other.values[otherNode], Math.max(this.levels[node], other.levels[otherNode]))) {
+ return this.makeSiblings(node, this.copyRecursive(other, otherNode));
+ }
+ else {
+ int left;
+ long branch;
+ if (this.levels[node] > other.levels[otherNode]) {
+ branch = other.values[otherNode] & this.getBranchMask(this.levels[node]);
+ if (branch == 0L) {
+ left = this.merge(this.lefts[node], other, otherNode);
+ this.lefts[node] = left;
+ }
+ else {
+ left = this.merge(this.rights[node], other, otherNode);
+ this.rights[node] = left;
+ }
+
+ return node;
+ }
+ else if (this.levels[node] < other.levels[otherNode]) {
+ branch = this.values[node] & this.getBranchMask(other.levels[otherNode]);
+ int right;
+ if (branch == 0L) {
+ left = this.merge(node, other, other.lefts[otherNode]);
+ right = this.copyRecursive(other, other.rights[otherNode]);
+ }
+ else {
+ left = this.copyRecursive(other, other.lefts[otherNode]);
+ right = this.merge(node, other, other.rights[otherNode]);
+ }
+
+ int result = this.createNode(other.values[otherNode], other.levels[otherNode], other.counts[otherNode]);
+ this.lefts[result] = left;
+ this.rights[result] = right;
+ return result;
+ }
+ else {
+ this.weightedCount += other.counts[otherNode];
+ double[] var10000 = this.counts;
+ var10000[node] += other.counts[otherNode];
+ left = this.merge(this.lefts[node], other, other.lefts[otherNode]);
+ int right = this.merge(this.rights[node], other, other.rights[otherNode]);
+ this.lefts[node] = left;
+ this.rights[node] = right;
+ return node;
+ }
+ }
+ }
+
+ private static boolean inSameSubtree(long bitsA, long bitsB, int level)
+ {
+ return level == 64 || bitsA >>> level == bitsB >>> level;
+ }
+
+ private int copyRecursive(QuantileDigest other, int otherNode)
+ {
+ if (otherNode == -1) {
+ return otherNode;
+ }
+ else {
+ int node = this.createNode(other.values[otherNode], other.levels[otherNode], other.counts[otherNode]);
+ int right;
+ if (other.lefts[otherNode] != -1) {
+ right = this.copyRecursive(other, other.lefts[otherNode]);
+ this.lefts[node] = right;
+ }
+
+ if (other.rights[otherNode] != -1) {
+ right = this.copyRecursive(other, other.rights[otherNode]);
+ this.rights[node] = right;
+ }
+
+ return node;
+ }
+ }
+
+ private int tryRemove(int node)
+ {
+ Preconditions.checkArgument(node != -1, "node is -1");
+ int left = this.lefts[node];
+ int right = this.rights[node];
+ if (left == -1 && right == -1) {
+ this.remove(node);
+ return -1;
+ }
+ else if (left != -1 && right != -1) {
+ this.counts[node] = 0.0D;
+ return node;
+ }
+ else {
+ this.remove(node);
+ return left != -1 ? left : right;
+ }
+ }
+
+ private void remove(int node)
+ {
+ if (node == this.nextNode - 1) {
+ --this.nextNode;
+ }
+ else {
+ this.pushFree(node);
+ }
+
+ if (node == this.root) {
+ this.root = -1;
+ }
+ }
+
+ private void pushFree(int node)
+ {
+ this.lefts[node] = this.firstFree;
+ this.firstFree = node;
+ ++this.freeCount;
+ }
+
+ private int popFree()
+ {
+ int node = this.firstFree;
+ if (node == -1) {
+ return node;
+ }
+ else {
+ this.firstFree = this.lefts[this.firstFree];
+ --this.freeCount;
+ return node;
+ }
+ }
+
+ private void postOrderTraversal(int node, QuantileDigest.Callback callback)
+ {
+ this.postOrderTraversal(node, callback, QuantileDigest.TraversalOrder.FORWARD);
+ }
+
+ private void postOrderTraversal(int node, QuantileDigest.Callback callback, QuantileDigest.TraversalOrder order)
+ {
+ if (order == QuantileDigest.TraversalOrder.FORWARD) {
+ this.postOrderTraversal(node, callback, this.lefts, this.rights);
+ }
+ else {
+ this.postOrderTraversal(node, callback, this.rights, this.lefts);
+ }
+ }
+
+ private boolean postOrderTraversal(int node, QuantileDigest.Callback callback, int[] lefts, int[] rights)
+ {
+ if (node == -1) {
+ return false;
+ }
+ else {
+ int first = lefts[node];
+ int second = rights[node];
+ if (first != -1 && !this.postOrderTraversal(first, callback, lefts, rights)) {
+ return false;
+ }
+ else {
+ return second != -1 && !this.postOrderTraversal(second, callback, lefts, rights) ? false : callback.process(node);
+ }
+ }
+ }
+
+ public double getConfidenceFactor()
+ {
+ return this.computeMaxPathWeight(this.root) * 1.0D / this.weightedCount;
+ }
+
+ @VisibleForTesting
+ boolean equivalent(QuantileDigest other)
+ {
+ return this.getNodeCount() == other.getNodeCount() && this.min == other.min && this.max == other.max && this.weightedCount == other.weightedCount && this.alpha == other.alpha;
+ }
+
+ private void rescaleToCommonLandmark(QuantileDigest one, QuantileDigest two)
+ {
+ long nowInSeconds = TimeUnit.NANOSECONDS.toSeconds(this.ticker.read());
+ long targetLandmark = Math.max(one.landmarkInSeconds, two.landmarkInSeconds);
+ if (nowInSeconds - targetLandmark >= 50L) {
+ targetLandmark = nowInSeconds;
+ }
+
+ if (targetLandmark != one.landmarkInSeconds) {
+ one.rescale(targetLandmark);
+ }
+
+ if (targetLandmark != two.landmarkInSeconds) {
+ two.rescale(targetLandmark);
+ }
+ }
+
+ private double computeMaxPathWeight(int node)
+ {
+ if (node != -1 && this.levels[node] != 0) {
+ double leftMaxWeight = this.computeMaxPathWeight(this.lefts[node]);
+ double rightMaxWeight = this.computeMaxPathWeight(this.rights[node]);
+ return Math.max(leftMaxWeight, rightMaxWeight) + this.counts[node];
+ }
+ else {
+ return 0.0D;
+ }
+ }
+
+ @VisibleForTesting
+ void validate()
+ {
+ AtomicDouble sum = new AtomicDouble();
+ AtomicInteger nodeCount = new AtomicInteger();
+ Set