/*
 * Licensed to the Apache Software Foundation (ASF) 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 org.apache.commons.statistics.descriptive;

import java.util.Arrays;
import java.util.Objects;
import java.util.function.IntToDoubleFunction;
import org.apache.commons.numbers.arrays.Selection;

/**
 * Provides quantile computation.
 *
 * <p>For values of length {@code n}:
 * <ul>
 * <li>The result is {@code NaN} if {@code n = 0}.
 * <li>The result is {@code values[0]} if {@code n = 1}.
 * <li>Otherwise the result is computed using the {@link EstimationMethod}.
 * </ul>
 *
 * <p>Computation of multiple quantiles and will handle duplicate and unordered
 * probabilities. Passing ordered probabilities is recommended if the order is already
 * known as this can improve efficiency; for example using uniform spacing through the
 * array data, or to identify extreme values from the data such as {@code [0.001, 0.999]}.
 *
 * <p>This implementation respects the ordering imposed by
 * {@link Double#compare(double, double)} for {@code NaN} values. If a {@code NaN} occurs
 * in the selected positions in the fully sorted values then the result is {@code NaN}.
 *
 * <p>The {@link NaNPolicy} can be used to change the behaviour on {@code NaN} values.
 *
 * <p>Instances of this class are immutable and thread-safe.
 *
 * @see #with(NaNPolicy)
 * @see <a href="http://en.wikipedia.org/wiki/Quantile">Quantile (Wikipedia)</a>
 * @since 1.1
 */
public final class Quantile {
    /** Message when the probability is not in the range {@code [0, 1]}. */
    private static final String INVALID_PROBABILITY = "Invalid probability: ";
    /** Message when no probabilities are provided for the varargs method. */
    private static final String NO_PROBABILITIES_SPECIFIED = "No probabilities specified";
    /** Message when the size is not valid. */
    private static final String INVALID_SIZE = "Invalid size: ";
    /** Message when the number of probabilities in a range is not valid. */
    private static final String INVALID_NUMBER_OF_PROBABILITIES = "Invalid number of probabilities: ";

    /** Default instance. Method 8 is recommended by Hyndman and Fan. */
    private static final Quantile DEFAULT = new Quantile(false, NaNPolicy.INCLUDE, EstimationMethod.HF8);

    /** Flag to indicate if the data should be copied. */
    private final boolean copy;
    /** NaN policy for floating point data. */
    private final NaNPolicy nanPolicy;
    /** Transformer for NaN data. */
    private final NaNTransformer nanTransformer;
    /** Estimation type used to determine the value from the quantile. */
    private final EstimationMethod estimationType;

    /**
     * @param copy Flag to indicate if the data should be copied.
     * @param nanPolicy NaN policy.
     * @param estimationType Estimation type used to determine the value from the quantile.
     */
    private Quantile(boolean copy, NaNPolicy nanPolicy, EstimationMethod estimationType) {
        this.copy = copy;
        this.nanPolicy = nanPolicy;
        this.estimationType = estimationType;
        nanTransformer = NaNTransformers.createNaNTransformer(nanPolicy, copy);
    }

    /**
     * Return a new instance with the default options.
     *
     * <ul>
     * <li>{@linkplain #withCopy(boolean) Copy = false}
     * <li>{@linkplain #with(NaNPolicy) NaN policy = include}
     * <li>{@linkplain #with(EstimationMethod) Estimation method = HF8}
     * </ul>
     *
     * <p>Note: The default options configure for processing in-place and including
     * {@code NaN} values in the data. This is the most efficient mode and has the
     * smallest memory consumption.
     *
     * @return the quantile implementation
     * @see #withCopy(boolean)
     * @see #with(NaNPolicy)
     * @see #with(EstimationMethod)
     */
    public static Quantile withDefaults() {
        return DEFAULT;
    }

    /**
     * Return an instance with the configured copy behaviour. If {@code false} then
     * the input array will be modified by the call to evaluate the quantiles; otherwise
     * the computation uses a copy of the data.
     *
     * @param v Value.
     * @return an instance
     */
    public Quantile withCopy(boolean v) {
        return new Quantile(v, nanPolicy, estimationType);
    }

    /**
     * Return an instance with the configured {@link NaNPolicy}.
     *
     * <p>Note: This implementation respects the ordering imposed by
     * {@link Double#compare(double, double)} for {@code NaN} values: {@code NaN} is
     * considered greater than all other values, and all {@code NaN} values are equal. The
     * {@link NaNPolicy} changes the computation of the statistic in the presence of
     * {@code NaN} values.
     *
     * <ul>
     * <li>{@link NaNPolicy#INCLUDE}: {@code NaN} values are moved to the end of the data;
     * the size of the data <em>includes</em> the {@code NaN} values and the quantile will be
     * {@code NaN} if any value used for quantile interpolation is {@code NaN}.
     * <li>{@link NaNPolicy#EXCLUDE}: {@code NaN} values are moved to the end of the data;
     * the size of the data <em>excludes</em> the {@code NaN} values and the quantile will
     * never be {@code NaN} for non-zero size. If all data are {@code NaN} then the size is zero
     * and the result is {@code NaN}.
     * <li>{@link NaNPolicy#ERROR}: An exception is raised if the data contains {@code NaN}
     * values.
     * </ul>
     *
     * <p>Note that the result is identical for all policies if no {@code NaN} values are present.
     *
     * @param v Value.
     * @return an instance
     */
    public Quantile with(NaNPolicy v) {
        return new Quantile(copy, Objects.requireNonNull(v), estimationType);
    }

    /**
     * Return an instance with the configured {@link EstimationMethod}.
     *
     * @param v Value.
     * @return an instance
     */
    public Quantile with(EstimationMethod v) {
        return new Quantile(copy, nanPolicy, Objects.requireNonNull(v));
    }

    /**
     * Generate {@code n} evenly spaced probabilities in the range {@code [0, 1]}.
     *
     * <pre>
     * 1/(n + 1), 2/(n + 1), ..., n/(n + 1)
     * </pre>
     *
     * @param n Number of probabilities.
     * @return the probabilities
     * @throws IllegalArgumentException if {@code n < 1}
     */
    public static double[] probabilities(int n) {
        checkNumberOfProbabilities(n);
        final double c1 = n + 1.0;
        final double[] p = new double[n];
        for (int i = 0; i < n; i++) {
            p[i] = (i + 1.0) / c1;
        }
        return p;
    }

    /**
     * Generate {@code n} evenly spaced probabilities in the range {@code [p1, p2]}.
     *
     * <pre>
     * w = p2 - p1
     * p1 + w/(n + 1), p1 + 2w/(n + 1), ..., p1 + nw/(n + 1)
     * </pre>
     *
     * @param n Number of probabilities.
     * @param p1 Lower probability.
     * @param p2 Upper probability.
     * @return the probabilities
     * @throws IllegalArgumentException if {@code n < 1}; if the probabilities are not in the
     * range {@code [0, 1]}; or {@code p2 <= p1}.
     */
    public static double[] probabilities(int n, double p1, double p2) {
        checkProbability(p1);
        checkProbability(p2);
        if (p2 <= p1) {
            throw new IllegalArgumentException("Invalid range: [" + p1 + ", " + p2 + "]");
        }
        final double[] p = probabilities(n);
        for (int i = 0; i < n; i++) {
            p[i] = (1 - p[i]) * p1 + p[i] * p2;
        }
        return p;
    }

    /**
     * Evaluate the {@code p}-th quantile of the values.
     *
     * <p>Note: This method may partially sort the input values if not configured to
     * {@link #withCopy(boolean) copy} the input data.
     *
     * <p><strong>Performance</strong>
     *
     * <p>It is not recommended to use this method for repeat calls for different quantiles
     * within the same values. The {@link #evaluate(double[], double...)} method should be used
     * which provides better performance.
     *
     * @param values Values.
     * @param p Probability for the quantile to compute.
     * @return the quantile
     * @throws IllegalArgumentException if the probability {@code p} is not in the range {@code [0, 1]}
     * @see #evaluate(double[], double...)
     */
    public double evaluate(double[] values, double p) {
        checkProbability(p);
        // Floating-point data handling
        final int[] bounds = new int[1];
        final double[] x = nanTransformer.apply(values, bounds);
        final int n = bounds[0];
        // Special cases
        if (n <= 1) {
            return n == 0 ? Double.NaN : x[0];
        }
        final double pos = estimationType.index(p, n);
        final int i = (int) pos;

        // Partition and compute
        if (pos > i) {
            Selection.select(x, 0, n, new int[] {i, i + 1});
            return Interpolation.interpolate(x[i], x[i + 1], pos - i);
        }
        Selection.select(x, 0, n, i);
        return x[i];
    }

    /**
     * Evaluate the {@code p}-th quantiles of the values.
     *
     * <p>Note: This method may partially sort the input values if not configured to
     * {@link #withCopy(boolean) copy} the input data.
     *
     * @param values Values.
     * @param p Probabilities for the quantiles to compute.
     * @return the quantiles
     * @throws IllegalArgumentException if any probability {@code p} is not in the range {@code [0, 1]};
     * or no probabilities are specified.
     */
    public double[] evaluate(double[] values, double... p) {
        checkProbabilities(p);
        // Floating-point data handling
        final int[] bounds = new int[1];
        final double[] x = nanTransformer.apply(values, bounds);
        final int n = bounds[0];
        // Special cases
        final double[] q = new double[p.length];
        if (n <= 1) {
            Arrays.fill(q, n == 0 ? Double.NaN : x[0]);
            return q;
        }

        // Collect interpolation positions. We use the output q as storage.
        final int[] indices = computeIndices(n, p, q);

        // Partition
        Selection.select(x, 0, n, indices);

        // Compute
        for (int k = 0; k < p.length; k++) {
            final int i = (int) q[k];
            if (q[k] > i) {
                q[k] = Interpolation.interpolate(x[i], x[i + 1], q[k] - i);
            } else {
                q[k] = x[i];
            }
        }
        return q;
    }

    /**
     * Evaluate the {@code p}-th quantile of the values.
     *
     * <p>Note: This method may partially sort the input values if not configured to
     * {@link #withCopy(boolean) copy} the input data.
     *
     * <p><strong>Performance</strong>
     *
     * <p>It is not recommended to use this method for repeat calls for different quantiles
     * within the same values. The {@link #evaluate(int[], double...)} method should be used
     * which provides better performance.
     *
     * @param values Values.
     * @param p Probability for the quantile to compute.
     * @return the quantile
     * @throws IllegalArgumentException if the probability {@code p} is not in the range {@code [0, 1]}
     * @see #evaluate(int[], double...)
     */
    public double evaluate(int[] values, double p) {
        checkProbability(p);
        final int n = values.length;
        // Special cases
        if (n <= 1) {
            return n == 0 ? Double.NaN : values[0];
        }
        final double pos = estimationType.index(p, n);
        final int i = (int) pos;

        // Partition and compute
        final int[] x = copy ? values.clone() : values;
        if (pos > i) {
            Selection.select(x, 0, n, new int[] {i, i + 1});
            return Interpolation.interpolate(x[i], x[i + 1], pos - i);
        }
        Selection.select(x, 0, n, i);
        return x[i];
    }

    /**
     * Evaluate the {@code p}-th quantiles of the values.
     *
     * <p>Note: This method may partially sort the input values if not configured to
     * {@link #withCopy(boolean) copy} the input data.
     *
     * @param values Values.
     * @param p Probabilities for the quantiles to compute.
     * @return the quantiles
     * @throws IllegalArgumentException if any probability {@code p} is not in the range {@code [0, 1]};
     * or no probabilities are specified.
     */
    public double[] evaluate(int[] values, double... p) {
        checkProbabilities(p);
        final int n = values.length;
        // Special cases
        final double[] q = new double[p.length];
        if (n <= 1) {
            Arrays.fill(q, n == 0 ? Double.NaN : values[0]);
            return q;
        }

        // Collect interpolation positions. We use the output q as storage.
        final int[] indices = computeIndices(n, p, q);

        // Partition
        final int[] x = copy ? values.clone() : values;
        Selection.select(x, 0, n, indices);

        // Compute
        for (int k = 0; k < p.length; k++) {
            final int i = (int) q[k];
            if (q[k] > i) {
                q[k] = Interpolation.interpolate(x[i], x[i + 1], q[k] - i);
            } else {
                q[k] = x[i];
            }
        }
        return q;
    }

    /**
     * Evaluate the {@code p}-th quantile of the values.
     *
     * <p>This method can be used when the values of known size are already sorted.
     *
     * <pre>{@code
     * short[] x = ...
     * Arrays.sort(x);
     * double q = Quantile.withDefaults().evaluate(x.length, i -> x[i], 0.05);
     * }</pre>
     *
     * @param n Size of the values.
     * @param values Values function.
     * @param p Probability for the quantile to compute.
     * @return the quantile
     * @throws IllegalArgumentException if {@code size < 0}; or if the probability {@code p} is
     * not in the range {@code [0, 1]}.
     */
    public double evaluate(int n, IntToDoubleFunction values, double p) {
        checkSize(n);
        checkProbability(p);
        // Special case
        if (n <= 1) {
            return n == 0 ? Double.NaN : values.applyAsDouble(0);
        }
        final double pos = estimationType.index(p, n);
        final int i = (int) pos;
        final double v1 = values.applyAsDouble(i);
        if (pos > i) {
            final double v2 = values.applyAsDouble(i + 1);
            return Interpolation.interpolate(v1, v2, pos - i);
        }
        return v1;
    }

    /**
     * Evaluate the {@code p}-th quantiles of the values.
     *
     * <p>This method can be used when the values of known size are already sorted.
     *
     * <pre>{@code
     * short[] x = ...
     * Arrays.sort(x);
     * double[] q = Quantile.withDefaults().evaluate(x.length, i -> x[i], 0.25, 0.5, 0.75);
     * }</pre>
     *
     * @param n Size of the values.
     * @param values Values function.
     * @param p Probabilities for the quantiles to compute.
     * @return the quantiles
     * @throws IllegalArgumentException if {@code size < 0}; if any probability {@code p} is
     * not in the range {@code [0, 1]}; or no probabilities are specified.
     */
    public double[] evaluate(int n, IntToDoubleFunction values, double... p) {
        checkSize(n);
        checkProbabilities(p);
        // Special case
        final double[] q = new double[p.length];
        if (n <= 1) {
            Arrays.fill(q, n == 0 ? Double.NaN : values.applyAsDouble(0));
            return q;
        }
        for (int k = 0; k < p.length; k++) {
            final double pos = estimationType.index(p[k], n);
            final int i = (int) pos;
            final double v1 = values.applyAsDouble(i);
            if (pos > i) {
                final double v2 = values.applyAsDouble(i + 1);
                q[k] = Interpolation.interpolate(v1, v2, pos - i);
            } else {
                q[k] = v1;
            }
        }
        return q;
    }

    /**
     * Check the probability {@code p} is in the range {@code [0, 1]}.
     *
     * @param p Probability for the quantile to compute.
     * @throws IllegalArgumentException if the probability is not in the range {@code [0, 1]}
     */
    private static void checkProbability(double p) {
        // Logic negation will detect NaN
        if (!(p >= 0 && p <= 1)) {
            throw new IllegalArgumentException(INVALID_PROBABILITY + p);
        }
    }

    /**
     * Check the probabilities {@code p} are in the range {@code [0, 1]}.
     *
     * @param p Probabilities for the quantiles to compute.
     * @throws IllegalArgumentException if any probabilities {@code p} is not in the range {@code [0, 1]};
     * or no probabilities are specified.
     */
    private static void checkProbabilities(double... p) {
        if (p.length == 0) {
            throw new IllegalArgumentException(NO_PROBABILITIES_SPECIFIED);
        }
        for (final double pp : p) {
            checkProbability(pp);
        }
    }

    /**
     * Check the {@code size} is positive.
     *
     * @param n Size of the values.
     * @throws IllegalArgumentException if {@code size < 0}
     */
    private static void checkSize(int n) {
        if (n < 0) {
            throw new IllegalArgumentException(INVALID_SIZE + n);
        }
    }

    /**
     * Check the number of probabilities {@code n} is strictly positive.
     *
     * @param n Number of probabilities.
     * @throws IllegalArgumentException if {@code c < 1}
     */
    private static void checkNumberOfProbabilities(int n) {
        if (n < 1) {
            throw new IllegalArgumentException(INVALID_NUMBER_OF_PROBABILITIES + n);
        }
    }

    /**
     * Compute the indices required for quantile interpolation.
     *
     * <p>The zero-based interpolation index in {@code [0, n)} is
     * saved into the working array {@code q} for each {@code p}.
     *
     * @param n Size of the data.
     * @param p Probabilities for the quantiles to compute.
     * @param q Working array for quantiles.
     * @return the indices
     */
    private int[] computeIndices(int n, double[] p, double[] q) {
        final int[] indices = new int[p.length << 1];
        int count = 0;
        for (int k = 0; k < p.length; k++) {
            final double pos = estimationType.index(p[k], n);
            q[k] = pos;
            final int i = (int) pos;
            indices[count++] = i;
            if (pos > i) {
                // Require the next index for interpolation
                indices[count++] = i + 1;
            }
        }
        if (count < indices.length) {
            return Arrays.copyOf(indices, count);
        }
        return indices;
    }

    /**
     * Estimation methods for a quantile. Provides the nine quantile algorithms
     * defined in Hyndman and Fan (1996)[1] as {@code HF1 - HF9}.
     *
     * <p>Samples quantiles are defined by:
     *
     * <p>\[ Q(p) = (1 - \gamma) x_j + \gamma x_{j+1} \]
     *
     * <p>where \( \frac{j-m}{n} \leq p \le \frac{j-m+1}{n} \), \( x_j \) is the \( j \)th
     * order statistic, \( n \) is the sample size, the value of \( \gamma \) is a function
     * of \( j = \lfloor np+m \rfloor \) and \( g = np + m - j \), and \( m \) is a constant
     * determined by the sample quantile type.
     *
     * <p>Note that the real-valued position \( np + m \) is a 1-based index and
     * \( j \in [1, n] \). If the real valued position is computed as beyond the lowest or
     * highest values in the sample, this implementation will return the minimum or maximum
     * observation respectively.
     *
     * <p>Types 1, 2, and 3 are discontinuous functions of \( p \); types 4 to 9 are continuous
     * functions of \( p \).
     *
     * <p>For the continuous functions, the probability \( p_k \) is provided for the \( k \)-th order
     * statistic in size \( n \). Samples quantiles are equivalently obtained to \( Q(p) \) by
     * linear interpolation between points \( (p_k, x_k) \) and \( (p_{k+1}, x_{k+1}) \) for
     * any \( p_k \leq p \leq p_{k+1} \).
     *
     * <ol>
     * <li>Hyndman and Fan (1996)
     *     <i>Sample Quantiles in Statistical Packages.</i>
     *     The American Statistician, 50, 361-365.
     *     <a href="https://www.jstor.org/stable/2684934">doi.org/10.2307/2684934</a>
     * <li><a href="http://en.wikipedia.org/wiki/Quantile">Quantile (Wikipedia)</a>
     * </ol>
     */
    public enum EstimationMethod {
        /**
         * Inverse of the empirical distribution function.
         *
         * <p>\( m = 0 \). \( \gamma = 0 \) if \( g = 0 \), and 1 otherwise.
         */
        HF1 {
            @Override
            double position0(double p, int n) {
                // position = np + 0. This is 1-based so adjust to 0-based.
                return Math.ceil(n * p) - 1;
            }
        },
        /**
         * Similar to {@link #HF1} with averaging at discontinuities.
         *
         * <p>\( m = 0 \). \( \gamma = 0.5 \) if \( g = 0 \), and 1 otherwise.
         */
        HF2 {
            @Override
            double position0(double p, int n) {
                final double pos = n * p;
                // Average at discontinuities
                final int j = (int) pos;
                final double g = pos - j;
                if (g == 0) {
                    return j - 0.5;
                }
                // As HF1 : ceil(j + g) - 1
                return j;
            }
        },
        /**
         * The observation closest to \( np \). Ties are resolved to the nearest even order statistic.
         *
         * <p>\( m = -1/2 \). \( \gamma = 0 \) if \( g = 0 \) and \( j \) is even, and 1 otherwise.
         */
        HF3 {
            @Override
            double position0(double p, int n) {
                // Let rint do the work for ties to even
                return Math.rint(n * p) - 1;
            }
        },
        /**
         * Linear interpolation of the inverse of the empirical CDF.
         *
         * <p>\( m = 0 \). \( p_k = \frac{k}{n} \).
         */
        HF4 {
            @Override
            double position0(double p, int n) {
                // np + 0 - 1
                return n * p - 1;
            }
        },
        /**
         * A piecewise linear function where the knots are the values midway through the steps of
         * the empirical CDF. Proposed by Hazen (1914) and popular amongst hydrologists.
         *
         * <p>\( m = 1/2 \). \( p_k = \frac{k - 1/2}{n} \).
         */
        HF5 {
            @Override
            double position0(double p, int n) {
                // np + 0.5 - 1
                return n * p - 0.5;
            }
        },
        /**
         * Linear interpolation of the expectations for the order statistics for the uniform
         * distribution on [0,1]. Proposed by Weibull (1939).
         *
         * <p>\( m = p \). \( p_k = \frac{k}{n + 1} \).
         *
         * <p>This method computes the quantile as per the Apache Commons Math Percentile
         * legacy implementation.
         */
        HF6 {
            @Override
            double position0(double p, int n) {
                // np + p - 1
                return (n + 1) * p - 1;
            }
        },
        /**
         * Linear interpolation of the modes for the order statistics for the uniform
         * distribution on [0,1]. Proposed by Gumbull (1939).
         *
         * <p>\( m = 1 - p \). \( p_k = \frac{k - 1}{n - 1} \).
         */
        HF7 {
            @Override
            double position0(double p, int n) {
                // np + 1-p - 1
                return (n - 1) * p;
            }
        },
        /**
         * Linear interpolation of the approximate medians for order statistics.
         *
         * <p>\( m = (p + 1)/3 \). \( p_k = \frac{k - 1/3}{n + 1/3} \).
         *
         * <p>As per Hyndman and Fan (1996) this approach is most recommended as it provides
         * an approximate median-unbiased estimate regardless of distribution.
         */
        HF8 {
            @Override
            double position0(double p, int n) {
                return n * p + (p + 1) / 3 - 1;
            }
        },
        /**
         * Quantile estimates are approximately unbiased for the expected order statistics if
         * \( x \) is normally distributed.
         *
         * <p>\( m = p/4 + 3/8 \). \( p_k = \frac{k - 3/8}{n + 1/4} \).
         */
        HF9 {
            @Override
            double position0(double p, int n) {
                // np + p/4 + 3/8 - 1
                return (n + 0.25) * p - 0.625;
            }
        };

        /**
         * Finds the real-valued position for calculation of the quantile.
         *
         * <p>Return {@code i + g} such that the quantile value from sorted data is:
         *
         * <p>value = data[i] + g * (data[i+1] - data[i])
         *
         * <p>Warning: Interpolation should not use {@code data[i+1]} unless {@code g != 0}.
         *
         * <p>Note: In contrast to the definition of Hyndman and Fan in the class header
         * which uses a 1-based position, this is a zero based index. This change is for
         * convenience when addressing array positions.
         *
         * @param p p<sup>th</sup> quantile.
         * @param n Size.
         * @return a real-valued position (0-based) into the range {@code [0, n)}
         */
        abstract double position0(double p, int n);

        /**
         * Finds the index {@code i} and fractional part {@code g} of a real-valued position
         * to interpolate the quantile.
         *
         * <p>Return {@code i + g} such that the quantile value from sorted data is:
         *
         * <p>value = data[i] + g * (data[i+1] - data[i])
         *
         * <p>Note: Interpolation should not use {@code data[i+1]} unless {@code g != 0}.
         *
         * @param p p<sup>th</sup> quantile.
         * @param n Size.
         * @return index (in [0, n-1])
         */
        final double index(double p, int n) {
            final double pos = position0(p, n);
            // Bounds check in [0, n-1]
            if (pos < 0) {
                return 0;
            }
            if (pos > n - 1.0) {
                return n - 1.0;
            }
            return pos;
        }
    }
}
