RandomPercentile.java

/*
 * Licensed to the Hipparchus project under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The Hipparchus project 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
 *
 *      https://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 org.hipparchus.stat.descriptive.rank;

import java.io.Serializable;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
import java.util.NoSuchElementException;
import java.util.UUID;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.exception.MathIllegalArgumentException;
import org.hipparchus.exception.MathIllegalStateException;
import org.hipparchus.exception.NullArgumentException;
import org.hipparchus.random.RandomGenerator;
import org.hipparchus.random.Well19937c;
import org.hipparchus.stat.StatUtils;
import org.hipparchus.stat.descriptive.AbstractStorelessUnivariateStatistic;
import org.hipparchus.stat.descriptive.AggregatableStatistic;
import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathArrays;

/**
 * A {@link StorelessUnivariateStatistic} estimating percentiles using the
 * <a href=http://dimacs.rutgers.edu/~graham/pubs/papers/nquantiles.pdf>RANDOM</a>
 * Algorithm.
 * <p>
 * Storage requirements for the RANDOM algorithm depend on the desired
 * accuracy of quantile estimates. Quantile estimate accuracy is defined as follows.
 * <p>
 * Let \(X\) be the set of all data values consumed from the stream and let \(q\)
 * be a quantile (measured between 0 and 1) to be estimated. If <ul>
 * <li>\(\epsilon\) is the configured accuracy</li>
 * <li> \(\hat{q}\) is a RandomPercentile estimate for \(q\) (what is returned
 *      by {@link #getResult()} or {@link #getResult(double)}) with \(100q\) as
 *      actual parameter)</li>
 * <li> \(rank(\hat{q}) = |\{x \in X : x &lt; \hat{q}\}|\) is the actual rank of
 *      \(\hat{q}\) in the full data stream</li>
 * <li>\(n = |X|\) is the number of observations</li></ul>
 * then we can expect \((q - \epsilon)n &lt; rank(\hat{q}) &lt; (q + \epsilon)n\).
 * <p>
 * The algorithm maintains \(\left\lceil {log_{2}(1/\epsilon)}\right\rceil + 1\) buffers
 * of size \(\left\lceil {1/\epsilon \sqrt{log_2(1/\epsilon)}}\right\rceil\).  When
 * {@code epsilon} is set to the default value of \(10^{-4}\), this makes 15 buffers
 * of size 36,453.
 * <p>
 * The algorithm uses the buffers to maintain samples of data from the stream.  Until
 * all buffers are full, the entire sample is stored in the buffers.
 * If one of the {@code getResult} methods is called when all data are available in memory
 * and there is room to make a copy of the data (meaning the combined set of buffers is
 * less than half full), the {@code getResult} method delegates to a {@link Percentile}
 * instance to compute and return the exact value for the desired quantile.
 * For default {@code epsilon}, this means exact values will be returned whenever fewer than
 * \(\left\lceil {15 \times 36453 / 2} \right\rceil = 273,398\) values have been consumed
 * from the data stream.
 * <p>
 * When buffers become full, the algorithm merges buffers so that they effectively represent
 * a larger set of values than they can hold. Subsequently, data values are sampled from the
 * stream to fill buffers freed by merge operations.  Both the merging and the sampling
 * require random selection, which is done using a {@code RandomGenerator}.  To get
 * repeatable results for large data streams, users should provide {@code RandomGenerator}
 * instances with fixed seeds. {@code RandomPercentile} itself does not reseed or otherwise
 * initialize the {@code RandomGenerator} provided to it.  By default, it uses a
 * {@link Well19937c} generator with the default seed.
 * <p>
 * Note: This implementation is not thread-safe.
 */
public class RandomPercentile
    extends AbstractStorelessUnivariateStatistic implements StorelessUnivariateStatistic,
    AggregatableStatistic<RandomPercentile>, Serializable {

    /** Default quantile estimation error setting */
    public static final double DEFAULT_EPSILON = 1e-4;
    /** Serialization version id */
    private static final long serialVersionUID = 1L;
    /** Storage size of each buffer */
    private final int s;
    /** Maximum number of buffers minus 1 */
    private final int h;
    /** Data structure used to manage buffers */
    private final BufferMap bufferMap;
    /** Bound on the quantile estimation error */
    private final double epsilon;
    /** Source of random data */
    private final RandomGenerator randomGenerator;
    /** Number of elements consumed from the input data stream */
    private long n;
    /** Buffer currently being filled */
    private Buffer currentBuffer;

    /**
     * Constructs a {@code RandomPercentile} with quantile estimation error
     * {@code epsilon} using {@code randomGenerator} as its source of random data.
     *
     * @param epsilon bound on quantile estimation error (see class javadoc)
     * @param randomGenerator PRNG used in sampling and merge operations
     * @throws MathIllegalArgumentException if percentile is not in the range [0, 100]
     */
    public RandomPercentile(double epsilon, RandomGenerator randomGenerator) {
        if (epsilon <= 0) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
                                                   epsilon, 0);
        }
        this.h = (int) FastMath.ceil(log2(1/epsilon));
        this.s = (int) FastMath.ceil(FastMath.sqrt(log2(1/epsilon)) / epsilon);
        this.randomGenerator = randomGenerator;
        bufferMap = new BufferMap(h + 1, s, randomGenerator);
        currentBuffer = bufferMap.create(0);
        this.epsilon = epsilon;
    }

    /**
     * Constructs a {@code RandomPercentile} with default estimation error
     * using {@code randomGenerator} as its source of random data.
     *
     * @param randomGenerator PRNG used in sampling and merge operations
     * @throws MathIllegalArgumentException if percentile is not in the range [0, 100]
     */
    public RandomPercentile(RandomGenerator randomGenerator) {
        this(DEFAULT_EPSILON, randomGenerator);
    }

    /**
     * Constructs a {@code RandomPercentile} with quantile estimation error
     * {@code epsilon} using the default PRNG as source of random data.
     *
     * @param epsilon bound on quantile estimation error (see class javadoc)
     * @throws MathIllegalArgumentException if percentile is not in the range [0, 100]
     */
    public RandomPercentile(double epsilon) {
        this(epsilon, new Well19937c());
    }

    /**
     * Constructs a {@code RandomPercentile} with quantile estimation error
     * set to the default ({@link #DEFAULT_EPSILON}), using the default PRNG
     * as source of random data.
     */
    public RandomPercentile() {
        this(DEFAULT_EPSILON, new Well19937c());
    }

    /**
     * Copy constructor, creates a new {@code RandomPercentile} identical
     * to the {@code original}.  Note: the RandomGenerator used by the new
     * instance is referenced, not copied - i.e., the new instance shares
     * a generator with the original.
     *
     * @param original the {@code PSquarePercentile} instance to copy
     */
    public RandomPercentile(RandomPercentile original) {
        super();
        this.h = original.h;
        this.n = original.n;
        this.s = original.s;
        this.epsilon = original.epsilon;
        this.bufferMap = new BufferMap(original.bufferMap);
        this.randomGenerator = original.randomGenerator;
        Iterator<Buffer> iterator = bufferMap.iterator();
        Buffer current = null;
        Buffer curr = null;
        // See if there is a partially filled buffer - that will be currentBuffer
        while (current == null && iterator.hasNext()) {
            curr = iterator.next();
            if (curr.hasCapacity()) {
                current = curr;
            }
        }
        // If there is no partially filled buffer, just assign the last one.
        // Next increment() will find no capacity and create a new one or trigger
        // a merge.
        this.currentBuffer = current == null ? curr : current;
    }

    @Override
    public long getN() {
        return n;
    }

    /**
     * Returns an estimate of the given percentile, computed using the designated
     * array segment as input data.
     *
     * @param values source of input data
     * @param begin position of the first element of the values array to include
     * @param length number of array elements to include
     * @param percentile desired percentile (scaled 0 - 100)
     *
     * @return estimated percentile
     * @throws MathIllegalArgumentException if percentile is out of the range [0, 100]
     */
    public double evaluate(final double percentile, final double[] values, final int begin, final int length)
        throws MathIllegalArgumentException {
        if (MathArrays.verifyValues(values, begin, length)) {
            RandomPercentile randomPercentile = new RandomPercentile(this.epsilon,
                                                                     this.randomGenerator);
            randomPercentile.incrementAll(values, begin, length);
            return randomPercentile.getResult(percentile);
        }
        return Double.NaN;
    }

    /**
     * Returns an estimate of the median, computed using the designated
     * array segment as input data.
     *
     * @param values source of input data
     * @param begin position of the first element of the values array to include
     * @param length number of array elements to include
     *
     * @return estimated percentile
     * @throws MathIllegalArgumentException if percentile is out of the range [0, 100]
     */
    @Override
    public double evaluate(final double[] values, final int begin, final int length) {
        return evaluate(50d, values, begin, length);
    }

    /**
     * Returns an estimate of percentile over the given array.
     *
     * @param values source of input data
     * @param percentile desired percentile (scaled 0 - 100)
     *
     * @return estimated percentile
     * @throws MathIllegalArgumentException if percentile is out of the range [0, 100]
     */
    public double evaluate(final double percentile, final double[] values) {
        return evaluate(percentile, values, 0, values.length);
    }

    @Override
    public RandomPercentile copy() {
        return new RandomPercentile(this);
    }

    @Override
    public void clear() {
        n = 0;
        bufferMap.clear();
        currentBuffer = bufferMap.create(0);
    }

    /**
     * Returns an estimate of the median.
     */
    @Override
    public double getResult() {
        return getResult(50d);

    }

    /**
     * Returns an estimate of the given percentile.
     *
     * @param percentile desired percentile (scaled 0 - 100)
     * @return estimated percentile
     * @throws MathIllegalArgumentException if percentile is out of the range [0, 100]
     */
    public double getResult(double percentile) {
        if (percentile > 100 || percentile < 0) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE,
                                                   percentile, 0, 100);
        }
        // Convert to internal quantile scale
        final double q = percentile / 100;
        // First get global min and max to bound search.
        double min = Double.POSITIVE_INFINITY;
        double max = Double.NEGATIVE_INFINITY;
        double bMin;
        double bMax;
        Iterator<Buffer> bufferIterator = bufferMap.iterator();
        while (bufferIterator.hasNext()) {
            Buffer buffer = bufferIterator.next();
            bMin = buffer.min();
            if (bMin < min) {
                min = bMin;
            }
            bMax = buffer.max();
            if (bMax > max) {
                max = bMax;
            }
        }

        // Handle degenerate cases
        if (Double.compare(q, 0d) == 0 || n == 1) {
            return min;
        }
        if (Double.compare(q, 1) == 0) {
            return max;
        }
        if (n == 0) {
            return Double.NaN;
        }

        // See if we have all data in memory and enough free memory to copy.
        // If so, use Percentile to perform exact computation.
        if (bufferMap.halfEmpty()) {
            return new Percentile(percentile).evaluate(bufferMap.levelZeroData());
        }

        // Compute target rank
        final double targetRank = q * n;

        // Start with initial guess min + quantile * (max - min).
        double estimate = min + q * (max - min);
        double estimateRank = getRank(estimate);
        double lower;
        double upper;
        if (estimateRank > targetRank) {
            upper = estimate;
            lower = min;
        } else if (estimateRank < targetRank) {
            lower = estimate;
            upper = max;
        } else {
            return estimate;
        }
        final double eps = epsilon / 2;
        final double rankTolerance = eps * n;
        final double minWidth = eps / n;
        double intervalWidth = FastMath.abs(upper - lower);
        while (FastMath.abs(estimateRank - targetRank) > rankTolerance && intervalWidth > minWidth) {
            if (estimateRank > targetRank) {
                upper = estimate;
            } else {
                lower = estimate;
            }
            intervalWidth = upper - lower;
            estimate = lower + intervalWidth / 2;
            estimateRank = getRank(estimate);
        }
        return estimate;
    }

    /**
     * Gets the estimated rank of {@code value}, i.e.  \(|\{x \in X : x &lt; value\}|\)
     * where \(X\) is the set of values that have been consumed from the stream.
     *
     * @param value value whose overall rank is sought
     * @return estimated number of sample values that are strictly less than {@code value}
     */
    public double getRank(double value) {
        double rankSum = 0;
        Iterator<Buffer> bufferIterator = bufferMap.iterator();
        while (bufferIterator.hasNext()) {
            Buffer buffer = bufferIterator.next();
            rankSum += buffer.rankOf(value) * FastMath.pow(2, buffer.level);
        }
        return rankSum;
    }

    /**
     * Returns the estimated quantile position of value in the dataset.
     * Specifically, what is returned is an estimate of \(|\{x \in X : x &lt; value\}| / |X|\)
     * where \(X\) is the set of values that have been consumed from the stream.
     *
     * @param value value whose quantile rank is sought.
     * @return estimated proportion of sample values that are strictly less than {@code value}
     */
    public double getQuantileRank(double value) {
        return getRank(value) / getN();
    }

    @Override
    public void increment(double d) {
        n++;
        if (!currentBuffer.hasCapacity()) { // Need to get a new buffer to fill
            // First see if we have not yet created all the buffers
            if (bufferMap.canCreate()) {
                final int level = (int) Math.ceil(Math.max(0, log2(n/(s * FastMath.pow(2, h - 1)))));
                currentBuffer = bufferMap.create(level);
            } else { // All buffers have been created - need to merge to free one
                currentBuffer = bufferMap.merge();
            }
        }
        currentBuffer.consume(d);
    }

    /**
     * Maintains a buffer of values sampled from the input data stream.
     * <p>
     * The {@link #level} of a buffer determines its sampling frequency.
     * The {@link #consume(double)} method retains 1 out of every 2^level values
     * read in from the stream.
     * <p>
     * The {@link #size} of the buffer is the number of values that it can store
     * The buffer is considered full when it has consumed 2^level * size values.
     * <p>
     * The {@link #blockSize} of a buffer is 2^level.
     * The consume method starts each block by generating a random integer in
     * [0, blockSize - 1].  It then skips over all but the element with that offset
     * in the block, retaining only the selected value. So after 2^level * size
     * elements have been consumed, it will have retained size elements - one
     * from each 2^level block.
     * <p>
     * The {@link #mergeWith(Buffer)} method merges this buffer with another one,
     * The merge operation merges the data from the other buffer into this and clears
     * the other buffer (so it can consume data). Both buffers have their level
     * incremented by the merge. This operation is only used on full buffers.
     */
    private static class Buffer implements Serializable {
        /** Serialization version id */
        private static final long serialVersionUID = 1L;
        /** Number of values actually stored in the buffer */
        private final int size;
        /** Data sampled from the stream */
        private final double[] data;
        /** PRNG used for merges and stream sampling */
        private final RandomGenerator randomGenerator;
        /** Level of the buffer */
        private int level;
        /** Block size  = 2^level */
        private long blockSize;
        /** Next location in backing array for stored (taken) value */
        private int next;
        /** Number of values consumed in current 2^level block of values from the stream */
        private long consumed;
        /** Index of next value to take in current 2^level block */
        private long nextToTake;
        /** ID */
        private final UUID id;

        /**
         * Creates a new buffer capable of retaining size values with the given level.
         *
         * @param size number of values the buffer can retain
         * @param level the base 2 log of the sampling frequency
         *        (one out of every 2^level values is retained)
         * @param randomGenerator PRNG used for sampling and merge operations
         */
        Buffer(int size, int level, RandomGenerator randomGenerator) {
            this.size = size;
            data = new double[size];
            this.level = level;
            this.randomGenerator = randomGenerator;
            this.id = UUID.randomUUID();
            computeBlockSize();
        }

        /**
         * Sets blockSize and nextToTake based on level.
         */
        private void computeBlockSize() {
            if (level == 0) {
                blockSize = 1;
            } else {
                long product = 1;
                for (int i = 0; i < level; i++) {
                    product *= 2;
                }
                blockSize = product;
            }
            if (blockSize > 1) {
                nextToTake = randomGenerator.nextLong(blockSize);
            }
        }

        /**
         * Consumes a value from the input stream.
         * <p>
         * For each 2^level values consumed, one is added to the buffer.
         * The buffer is not considered full until 2^level * size values
         * have been consumed.
         * <p>
         * Sorts the data array if the consumption renders the buffer full.
         * <p>
         * There is no capacity check in this method.  Clients are expected
         * to use {@link #hasCapacity()} before invoking this method.  If
         * it is invoked on a full buffer, an ArrayIndexOutOfBounds exception
         * will result.
         *
         * @param value value to consume from the stream
         */
        public void consume(double value) {
            if (consumed == nextToTake) {
                data[next] = value;
                next++;
            }
            consumed++;
            if (consumed == blockSize) {
                if (next == size) {   // Buffer is full
                    Arrays.sort(data);
                } else {              // Reset in-block counter and nextToTake
                    consumed = 0;
                    if (blockSize > 1) {
                        nextToTake = randomGenerator.nextLong(blockSize);
                    }
                }
            }
        }

        /**
         * Merges this with other.
         * <p>
         * After the merge, this will be the merged buffer and other will be free.
         * Both will have level+1.
         * Post-merge, other can be used to accept new data.
         * <p>
         * The contents of the merged buffer (this after the merge) are determined
         * by randomly choosing one of the two retained elements in each of the
         * [0...size - 1] positions in the two buffers.
         * <p>
         * This and other must have the same level and both must be full.
         *
         * @param other initially full other buffer at the same level as this.
         * @throws MathIllegalArgumentException if either buffer is not full or they
         * have different levels
         */
        public void mergeWith(Buffer other) {
            // Make sure both this and other are full and have the same level
            if (this.hasCapacity() || other.hasCapacity() || other.level != this.level) {
                throw new MathIllegalArgumentException(LocalizedCoreFormats.INTERNAL_ERROR);
            }
            // Randomly select one of the two entries for each slot
            for (int i = 0; i < size; i++) {
                if (randomGenerator.nextBoolean()) {
                    data[i] = other.data[i];
                }
            }
            // Re-sort data
            Arrays.sort(data);
            // Bump level of both buffers
            other.setLevel(level + 1);
            this.setLevel(level + 1);
            // Clear the free one (and compute new blocksize)
            other.clear();
        }

        /**
         * Merge this into a higher-level buffer.
         * <p>
         * Does not alter this; but after the merge, higher may have some of its
         * data replaced by data from this.  Levels are not changed for either buffer.
         * <p>
         * Probability of selection into the newly constituted higher buffer is weighted
         * according to level. So for example, if this has level 0 and higher has level
         * 2, the ith element of higher is 4 times more likely than the corresponding
         * element of this to retain its spot.
         * <p>
         * This method is only used when aggregating RandomPercentile instances.
         * <p>
         * Preconditions:
         * <ol><li> this.level < higher.level </li>
         *     <li> this.size = higher.size </li>
         *     <li> Both buffers are full </li>
         * </ol>
         *
         * @param higher higher-level buffer to merge this into
         * @throws MathIllegalArgumentException if the buffers have different sizes,
         * either buffer is not full or this has level greater than or equal to higher
         */
        public void mergeInto(Buffer higher) {
            // Check preconditions
            if (this.size != higher.size || this.hasCapacity() || higher.hasCapacity() ||
                    this.level >= higher.level) {
                throw new MathIllegalArgumentException(LocalizedCoreFormats.INTERNAL_ERROR);
            }
            final int levelDifference = higher.level - this.level;
            int m = 1;
            for (int i = 0; i < levelDifference; i++) {
                m *= 2;
            }
            // Randomly select one of the two entries for each slot in higher, giving
            // m-times higher weight to the entries of higher.
            for (int i = 0; i < size; i++) {
                // data[i] <-> {0}, higher.data[i] <-> {1, ..., m}
                if (randomGenerator.nextInt(m + 1) == 0) {
                    higher.data[i] = data[i];
                }
            }
            // Resort higher's data
            Arrays.sort(higher.data);
        }

        /**
         * @return true if the buffer has capacity; false if it is full
         */
        public boolean hasCapacity() {
            // Buffer has capacity if it has not yet set all of its data
            // values or if it has but still has not finished its last block
            return next < size || consumed < blockSize;
        }

        /**
         * Sets the level of the buffer.
         *
         * @param level new level value
         */
        public void setLevel(int level) {
            this.level = level;
        }

        /**
         * Clears data, recomputes blockSize and resets consumed and nextToTake.
         */
        public void clear() {
            consumed = 0;
            next = 0;
            computeBlockSize();
        }

        /**
         * Returns a copy of the data that has been added to the buffer
         *
         * @return possibly unsorted copy of the portion of the buffer that has been filled
         */
        public double[] getData() {
            final double[] out = new double[next];
            System.arraycopy(data, 0, out, 0, next);
            return out;
        }

        /**
         * Returns the ordinal rank of value among the sampled values in this buffer.
         *
         * @param value value whose rank is sought
         * @return |{v in data : v < value}|
         */
        public int rankOf(double value) {
            int ret = 0;
            if (!hasCapacity()) { // Full sorted buffer, can do binary search
                ret = Arrays.binarySearch(data, value);
                if (ret < 0) {
                    return -ret - 1;
                } else {
                    return ret;
                }
            } else { // have to count - not sorted yet and can't sort yet
                for (int i = 0; i < next; i++) {
                    if (data[i] < value) {
                        ret++;
                    }
                }
                return ret;
            }
        }

        /**
         * @return the smallest value held in this buffer
         */
        public double min() {
            if (!hasCapacity()) {
                return data[0];
            } else {
                return StatUtils.min(getData());
            }
        }

        /**
         * @return the largest value held in this buffer
         */
        public double max() {
            if (!hasCapacity()) {
                return data[data.length - 1];
            } else {
                return StatUtils.max(getData());
            }
        }

        /**
         * @return the level of this buffer
         */
        public int getLevel() {
            return level;
        }

        /**
         * @return the id
         */
        public UUID getId() {
            return id;
        }
    }

    /**
     * Computes base 2 log of the argument.
     *
     * @param x input value
     * @return the value y such that 2^y = x
     */
    private static double log2(double x) {
        return Math.log(x) / Math.log(2);
    }

    /**
     * A map structure to hold the buffers.
     * Keys are levels and values are lists of buffers at the given level.
     * Overall capacity is limited by the total number of buffers.
     */
    private static class BufferMap implements Iterable<Buffer>, Serializable {
        /** Serialization version ID */
        private static final long serialVersionUID = 1L;
        /** Total number of buffers that can be created - cap for count */
        private final int capacity;
        /** PRNG used in merges */
        private final RandomGenerator randomGenerator;
        /** Total count of all buffers */
        private int count;
        /** Uniform buffer size */
        private final int bufferSize;
        /** Backing store for the buffer map. Keys are levels, values are lists of registered buffers. */
        private final Map<Integer,List<Buffer>> registry;
        /** Maximum buffer level */
        private int maxLevel;

        /**
         * Creates a BufferMap that can manage up to capacity buffers.
         * Buffers created by the pool with have size = buffersize.
         *
         * @param capacity cap on the number of buffers
         * @param bufferSize size of each buffer
         * @param randomGenerator RandomGenerator to use in merges
         */
        BufferMap(int capacity, int bufferSize, RandomGenerator randomGenerator) {
            this.bufferSize = bufferSize;
            this.capacity = capacity;
            this.randomGenerator = randomGenerator;
            this.registry = new HashMap<>();
        }

        /**
         * Copy constructor.
         *
         * @param original BufferMap to copy
         */
        BufferMap(BufferMap original) {
            super();
            this.bufferSize = original.bufferSize;
            this.capacity = original.capacity;
            this.count = 0;
            this.randomGenerator = original.randomGenerator;
            this.registry = new HashMap<>();
            Iterator<Buffer> iterator = original.iterator();
            while (iterator.hasNext()) {
                final Buffer current = iterator.next();
                // Create and register a new buffer at the same level
                final Buffer newCopy = create(current.getLevel());
                // Consume the data
                final double[] data = current.getData();
                for (double value : data) {
                    newCopy.consume(value);
                }
            }
        }

        /**
         * Tries to create a buffer with the given level.
         * <p>
         * If there is capacity to create a new buffer (i.e., fewer than
         * count have been created), a new buffer is created with the given
         * level, registered and returned.  If capacity has been reached,
         * null is returned.
         *
         * @param level level of the new buffer
         * @return an empty buffer or null if a buffer can't be provided
         */
        public Buffer create(int level) {
            if (!canCreate()) {
                return null;
            }
            count++;
            Buffer buffer = new Buffer(bufferSize, level, randomGenerator);
            List<Buffer> bufferList = registry.get(level);
            if (bufferList == null) {
                bufferList = new ArrayList<>();
                registry.put(level, bufferList);
            }
            bufferList.add(buffer);
            if (level > maxLevel) {
                maxLevel = level;
            }
            return buffer;
        }

        /**
         * Returns true if there is capacity to create a new buffer.
         *
         * @return true if fewer than capacity buffers have been created.
         */
        public boolean canCreate() {
            return count < capacity;
        }

        /**
         * Returns true if we have used less than half of the allocated storage.
         * <p>
         * Includes a check to make sure all buffers have level 0;
         * but this should always be the case.
         * <p>
         * When this method returns true, we have all consumed data in storage
         * and enough space to make a copy of the combined dataset.
         *
         * @return true if all buffers have level 0 and less than half of the
         * available storage has been used
         */
        public boolean halfEmpty() {
            return count * 2 < capacity &&
                    registry.size() == 1 &&
                    registry.containsKey(0);
        }

        /**
         * Returns a fresh copy of all data from level 0 buffers.
         *
         * @return combined data stored in all level 0 buffers
         */
        public double[] levelZeroData() {
            List<Buffer> levelZeroBuffers = registry.get(0);
            // First determine the combined size of the data
            int length = 0;
            for (Buffer buffer : levelZeroBuffers) {
                if (!buffer.hasCapacity()) { // full buffer
                    length += buffer.size;
                } else {
                    length += buffer.next;  // filled amount
                }
            }
            // Copy the data
            int pos = 0;
            int currLen;
            final double[] out = new double[length];
            for (Buffer buffer : levelZeroBuffers) {
                if (!buffer.hasCapacity()) {
                    currLen = buffer.size;
                } else {
                    currLen =  buffer.next;
                }
                System.arraycopy(buffer.data, 0, out, pos, currLen);
                pos += currLen;
            }
            return out;
        }

        /**
         * Finds the lowest level l where there exist at least two buffers,
         * merges them to create a new buffer with level l+1 and returns
         * a free buffer with level l+1.
         *
         * @return free buffer that can accept data
         */
        public Buffer merge() {
            int l = 0;
            List<Buffer> mergeCandidates = null;
            // Find the lowest level containing at least two buffers
            while (mergeCandidates == null && l <= maxLevel) {
                final List<Buffer> bufferList = registry.get(l);
                if (bufferList != null && bufferList.size() > 1) {
                    mergeCandidates = bufferList;
                } else {
                    l++;
                }
            }
            if (mergeCandidates == null) {
                // Should never happen
                throw new MathIllegalStateException(LocalizedCoreFormats.INTERNAL_ERROR);
            }
            Buffer buffer1 = mergeCandidates.get(0);
            Buffer buffer2 = mergeCandidates.get(1);
            // Remove buffers to be merged
            mergeCandidates.remove(0);
            mergeCandidates.remove(0);
            // If these are the last level-l buffers, remove the empty list
            if (registry.get(l).isEmpty()) {
                registry.remove(l);
            }
            // Merge the buffers
            buffer1.mergeWith(buffer2);
            // Now both buffers have level l+1; buffer1 is full and buffer2 is free.
            // Register both buffers
            register(buffer1);
            register(buffer2);

            // Return the free one
            return buffer2;
        }

        /**
         * Clears the buffer map.
         */
        public void clear() {
            for (List<Buffer> bufferList : registry.values()) {
                bufferList.clear();
            }
            registry.clear();
            count = 0;
        }

        /**
         * Registers a buffer.
         *
         * @param buffer Buffer to be registered.
         */
        public void register(Buffer buffer) {
            final int level = buffer.getLevel();
            List<Buffer> list = registry.get(level);
            if (list == null) {
                list = new ArrayList<>();
                registry.put(level, list);
                if (level > maxLevel) {
                    maxLevel = level;
                }
            }
            list.add(buffer);
        }

        /**
         * De-register a buffer, without clearing it.
         *
         * @param buffer Buffer to be de-registered
         * @throws IllegalStateException if the buffer is not registered
         */
        public void deRegister(Buffer buffer) {
            final Iterator<Buffer> iterator = registry.get(buffer.getLevel()).iterator();
            while (iterator.hasNext()) {
                if (iterator.next().getId().equals(buffer.getId())) {
                    iterator.remove();
                    return;
                }
            }
            throw new MathIllegalStateException(LocalizedCoreFormats.INTERNAL_ERROR);
        }

        /**
         * Returns an iterator over all of the buffers. Iteration goes by level, with
         * level 0 first.  Assumes there are no empty buffer lists.
         */
        @Override
        public Iterator<Buffer> iterator() {
            return new Iterator<Buffer>() {

                /** Outer loop iterator, from level to level. */
                private final Iterator<Integer> levelIterator = registry.keySet().iterator();

                /** List of buffers at current level. */
                private List<Buffer> currentList = registry.get(levelIterator.next());

                /** Inner loop iterator, from buffer to buffer. */
                private Iterator<Buffer> bufferIterator =
                        currentList == null ? null : currentList.iterator();

                @Override
                public boolean hasNext() {
                    if (bufferIterator == null) {
                        return false;
                    }
                    if (bufferIterator.hasNext()) {
                        return true;
                    }
                    // The current level iterator has just finished, try to bump level
                    if (levelIterator.hasNext()) {
                        List<Buffer> currentList = registry.get(levelIterator.next());
                        bufferIterator = currentList.iterator();
                        return true;
                    } else {
                        // Nothing left, signal this by nulling bufferIterator
                        bufferIterator = null;
                        return false;
                    }
                }

                @Override
                public Buffer next() {
                     if (hasNext()) {
                         return bufferIterator.next();
                     }
                     throw new NoSuchElementException();
                }

                @Override
                public void remove() {
                    throw new UnsupportedOperationException();
                }
            };
        }

        /**
         * Absorbs the data in other into this, merging buffers as necessary to trim
         * the aggregate down to capacity. This method is only used when aggregating
         * RandomPercentile instances.
         *
         * @param other other BufferMap to merge in
         */
        public void absorb(BufferMap other) {
            // Add all of other's buffers to the map - possibly exceeding cap
            boolean full = true;
            Iterator<Buffer> otherIterator = other.iterator();
            while (otherIterator.hasNext()) {
                Buffer buffer = otherIterator.next();
                if (buffer.hasCapacity()) {
                    full = false;
                }
                register(buffer);
                count++;
            }
            // Now eliminate the excess by merging
            while (count > capacity || (count == capacity && full)) {
                mergeUp();
                count--;
            }
        }

        /**
         * Find two buffers, first and second, of minimal level. Then merge
         * first into second and discard first.
         * <p>
         * If the buffers have different levels, make second the higher level
         * buffer and make probability of selection in the merge proportional
         * to level weight ratio.
         * <p>
         * This method is only used when aggregating RandomPercentile instances.
         */
        public void mergeUp() {
            // Find two minimum-level buffers to merge
            // Loop depends on two invariants:
            //   0) iterator goes in level order
            //   1) there are no empty lists in the registry
            Iterator<Buffer> bufferIterator = iterator();
            Buffer first = null;
            Buffer second = null;
            while ((first == null || second == null) && bufferIterator.hasNext()) {
                Buffer buffer = bufferIterator.next();
                if (!buffer.hasCapacity()) { // Skip not full buffers
                    if (first == null) {
                        first = buffer;
                    } else {
                        second = buffer;
                    }
                }
            }
            if (first == null || second == null || first.level > second.level) {
                throw new MathIllegalStateException(LocalizedCoreFormats.INTERNAL_ERROR);
            }
            // Merge first into second and deregister first.
            // Assumes that first has level <= second (checked above).
            if (first.getLevel() == second.getLevel()) {
                deRegister(first);
                deRegister(second);
                second.mergeWith(first);
                register(second);
            } else {
                deRegister(first);
                first.mergeInto(second);
            }
        }
    }

    /**
     * Computes the given percentile by combining the data from the collection
     * of aggregates. The result describes the combined sample of all data added
     * to any of the aggregates.
     *
     * @param percentile desired percentile (scaled 0-100)
     * @param aggregates RandomPercentile instances to combine data from
     * @return estimate of the given percentile using combined data from the aggregates
     * @throws MathIllegalArgumentException if percentile is out of the range [0, 100]
     */
    public double reduce(double percentile, Collection<RandomPercentile> aggregates) {
        if (percentile > 100 || percentile < 0) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE,
                                                   percentile, 0, 100);
        }

        // First see if we can copy all data and just compute exactly.
        // The following could be improved to verify that all have only level 0 buffers
        // and the sum of the data sizes is less than 1/2 total capacity.  Here we
        // just check that each of the aggregates is less than half full.
        Iterator<RandomPercentile> iterator = aggregates.iterator();
        boolean small = true;
        while (small && iterator.hasNext()) {
            small = iterator.next().bufferMap.halfEmpty();
        }
        if (small) {
            iterator = aggregates.iterator();
            double[] combined = {};
            while (iterator.hasNext()) {
               combined = MathArrays.concatenate(combined, iterator.next().bufferMap.levelZeroData());
            }
            final Percentile exactP = new Percentile(percentile);
            return exactP.evaluate(combined);
        }

        // Below largely duplicates code in getResult(percentile).
        // Common binary search code with function parameter should be factored out.

        // Get global max and min to bound binary search and total N
        double min = Double.POSITIVE_INFINITY;
        double max = Double.NEGATIVE_INFINITY;
        double combinedN = 0;
        iterator = aggregates.iterator();
        while (iterator.hasNext()) {
            final RandomPercentile curr = iterator.next();
            final double curMin = curr.getResult(0);
            final double curMax = curr.getResult(100);
            if (curMin < min) {
                min = curMin;
            }
            if (curMax > max) {
                max = curMax;
            }
            combinedN += curr.getN();
        }

        final double q = percentile / 100;
        // Handle degenerate cases
        if (Double.compare(q, 0d) == 0) {
            return min;
        }
        if (Double.compare(q, 1) == 0) {
            return max;
        }

        // Compute target rank
        final double targetRank = q * combinedN;

        // Perform binary search using aggregated rank computation
        // Start with initial guess min + quantile * (max - min).
        double estimate = min + q * (max - min);
        double estimateRank = getAggregateRank(estimate, aggregates);
        double lower;
        double upper;
        if (estimateRank > targetRank) {
            upper = estimate;
            lower = min;
        } else if (estimateRank < targetRank) {
            lower = estimate;
            upper = max;
        } else {
            return estimate;
        }
        final double eps = epsilon / 2;
        double intervalWidth = FastMath.abs(upper - lower);
        while (FastMath.abs(estimateRank / combinedN - q) > eps && intervalWidth > eps / combinedN) {
            if (estimateRank == targetRank) {
                return estimate;
            }
            if (estimateRank > targetRank) {
                upper = estimate;
            } else {
                lower = estimate;
            }
            intervalWidth = FastMath.abs(upper - lower);
            estimate = lower + intervalWidth / 2;
            estimateRank = getAggregateRank(estimate, aggregates);
        }
        return estimate;
    }

    /**
     * Computes the estimated rank of value in the combined dataset of the aggregates.
     * Sums the values from {@link #getRank(double)}.
     *
     * @param value value whose rank is sought
     * @param aggregates collection to aggregate rank over
     * @return estimated number of elements in the combined dataset that are less than value
     */
    public double getAggregateRank(double value, Collection<RandomPercentile> aggregates) {
        double result = 0;
        final Iterator<RandomPercentile> iterator = aggregates.iterator();
        while (iterator.hasNext()) {
            result += iterator.next().getRank(value);
        }
        return result;
    }

    /**
     * Returns the estimated quantile position of value in the combined dataset of the aggregates.
     * Specifically, what is returned is an estimate of \(|\{x \in X : x &lt; value\}| / |X|\)
     * where \(X\) is the set of values that have been consumed from all of the datastreams
     * feeding the aggregates.
     *
     * @param value value whose quantile rank is sought.
     * @param aggregates collection of RandomPercentile instances being combined
     * @return estimated proportion of combined sample values that are strictly less than {@code value}
     */
    public double getAggregateQuantileRank(double value, Collection<RandomPercentile> aggregates) {
        return getAggregateRank(value, aggregates) / getAggregateN(aggregates);
    }

    /**
     * Returns the total number of values that have been consumed by the aggregates.
     *
     * @param aggregates collection of RandomPercentile instances whose combined sample size is sought
     * @return total number of values that have been consumed by the aggregates
     */
    public double getAggregateN(Collection<RandomPercentile> aggregates) {
        double result = 0;
        final Iterator<RandomPercentile> iterator = aggregates.iterator();
        while (iterator.hasNext()) {
            result += iterator.next().getN();
        }
        return result;
    }

    /**
     * Aggregates the provided instance into this instance.
     * <p>
     * Other must have the same buffer size as this. If the combined data size
     * exceeds the maximum storage configured for this instance, buffers are
     * merged to create capacity. If all that is needed is computation of
     * aggregate results, {@link #reduce(double, Collection)} is faster,
     * may be more accurate and does not require the buffer sizes to be the same.
     *
     * @param other the instance to aggregate into this instance
     * @throws NullArgumentException if the input is null
     * @throws IllegalArgumentException if other has different buffer size than this
     */
    @Override
    public void aggregate(RandomPercentile other)
        throws NullArgumentException {
        if (other == null) {
            throw new NullArgumentException();
        }
        if (other.s != s) {
            throw new MathIllegalArgumentException(LocalizedCoreFormats.INTERNAL_ERROR);
        }
        bufferMap.absorb(other.bufferMap);
        n += other.n;
    }

    /**
     * Returns the maximum number of {@code double} values that a {@code RandomPercentile}
     * instance created with the given {@code epsilon} value will retain in memory.
     * <p>
     * If the number of values that have been consumed from the stream is less than 1/2
     * of this value, reported statistics are exact.
     *
     * @param epsilon bound on the relative quantile error (see class javadoc)
     * @return upper bound on the total number of primitive double values retained in memory
     * @throws MathIllegalArgumentException if epsilon is not in the interval (0,1)
     */
    public static long maxValuesRetained(double epsilon) {
        if (epsilon >= 1) {
            throw new MathIllegalArgumentException(
                    LocalizedCoreFormats.NUMBER_TOO_LARGE_BOUND_EXCLUDED, epsilon, 1);
        }
        if (epsilon <= 0) {
            throw new MathIllegalArgumentException(
                    LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED, epsilon, 0);
        }
        final long h = (long) FastMath.ceil(log2(1/epsilon));
        final long s = (long) FastMath.ceil(FastMath.sqrt(log2(1/epsilon)) / epsilon);
        return (h+1) * s;
    }
}