/*
 * 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.math3.distribution;

import java.io.Serializable;

import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.random.RandomGenerator;
import org.apache.commons.math3.util.FastMath;

Base class for probability distributions on the reals.
Default implementations are provided for some of the methods
that do not vary from distribution to distribution.
Since:
3.0/**
 * Base class for probability distributions on the reals.
 * Default implementations are provided for some of the methods
 * that do not vary from distribution to distribution.
 *
 * @since 3.0
 */
public abstract class AbstractRealDistribution
implements RealDistribution, Serializable {
    
Default accuracy. /** Default accuracy. */
    public static final double SOLVER_DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
    
Serializable version identifier /** Serializable version identifier */
    private static final long serialVersionUID = -38038050983108802L;
     
RandomData instance used to generate samples from the distribution.
Deprecated:
As of 3.1, to be removed in 4.0. Please use the random instance variable instead./**
      * RandomData instance used to generate samples from the distribution.
      * @deprecated As of 3.1, to be removed in 4.0. Please use the
      * {@link #random} instance variable instead.
      */
    @Deprecated
    protected org.apache.commons.math3.random.RandomDataImpl randomData =
        new org.apache.commons.math3.random.RandomDataImpl();

    
RNG instance used to generate samples from the distribution.
Since:
3.1/**
     * RNG instance used to generate samples from the distribution.
     * @since 3.1
     */
    protected final RandomGenerator random;

    
Solver absolute accuracy for inverse cumulative computation /** Solver absolute accuracy for inverse cumulative computation */
    private double solverAbsoluteAccuracy = SOLVER_DEFAULT_ABSOLUTE_ACCURACY;

    
Deprecated:
As of 3.1, to be removed in 4.0. Please use AbstractRealDistribution(RandomGenerator) instead./**
     * @deprecated As of 3.1, to be removed in 4.0. Please use
     * {@link #AbstractRealDistribution(RandomGenerator)} instead.
     */
    @Deprecated
    protected AbstractRealDistribution() {
        // Legacy users are only allowed to access the deprecated "randomData".
        // New users are forbidden to use this constructor.
        random = null;
    }
    
Params:
rng – Random number generator.
Since:
3.1/**
     * @param rng Random number generator.
     * @since 3.1
     */
    protected AbstractRealDistribution(RandomGenerator rng) {
        random = rng;
    }

    
{@inheritDoc}
The default implementation uses the identity
P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)
Deprecated:
As of 3.1 (to be removed in 4.0). Please use probability(double, double) instead./**
     * {@inheritDoc}
     *
     * The default implementation uses the identity
     * <p>{@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}</p>
     *
     * @deprecated As of 3.1 (to be removed in 4.0). Please use
     * {@link #probability(double,double)} instead.
     */
    @Deprecated
    public double cumulativeProbability(double x0, double x1) throws NumberIsTooLargeException {
        return probability(x0, x1);
    }

    
For a random variable X whose values are distributed according to this distribution, this method returns P(x0 < X <= x1). 
Params:
x0 – Lower bound (excluded).
x1 – Upper bound (included).
Throws:
NumberIsTooLargeException – if x0 > x1. The default implementation uses the identity P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)
Returns:
the probability that a random variable with this distribution takes a value between x0 and x1, excluding the lower and including the upper endpoint.
Since:
3.1/**
     * For a random variable {@code X} whose values are distributed according
     * to this distribution, this method returns {@code P(x0 < X <= x1)}.
     *
     * @param x0 Lower bound (excluded).
     * @param x1 Upper bound (included).
     * @return the probability that a random variable with this distribution
     * takes a value between {@code x0} and {@code x1}, excluding the lower
     * and including the upper endpoint.
     * @throws NumberIsTooLargeException if {@code x0 > x1}.
     *
     * The default implementation uses the identity
     * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)}
     *
     * @since 3.1
     */
    public double probability(double x0,
                              double x1) {
        if (x0 > x1) {
            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
                                                x0, x1, true);
        }
        return cumulativeProbability(x1) - cumulativeProbability(x0);
    }

    
{@inheritDoc}
The default implementation returns

RealDistribution.getSupportLowerBound() for p = 0,
RealDistribution.getSupportUpperBound() for p = 1.

/**
     * {@inheritDoc}
     *
     * The default implementation returns
     * <ul>
     * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li>
     * <li>{@link #getSupportUpperBound()} for {@code p = 1}.</li>
     * </ul>
     */
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if (p < 0.0 || p > 1.0) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();
        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();
        if (p == 1.0) {
            return upperBound;
        }

        final double mu = getNumericalMean();
        final double sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;
        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) ||
                             Double.isInfinite(sig) || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {
                lowerBound = -1.0;
                while (cumulativeProbability(lowerBound) >= p) {
                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {
                upperBound = 1.0;
                while (cumulativeProbability(upperBound) < p) {
                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {
            
{@inheritDoc} /** {@inheritDoc} */
            public double value(final double x) {
                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve,
                                                   lowerBound,
                                                   upperBound,
                                                   getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {
            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();
            if (x - dx >= getSupportLowerBound()) {
                double px = cumulativeProbability(x);
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x;
                    while (upperBound - lowerBound > dx) {
                        final double midPoint = 0.5 * (lowerBound + upperBound);
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }
                    return upperBound;
                }
            }
        }
        return x;
    }

    
Returns the solver absolute accuracy for inverse cumulative computation.
You can override this method in order to use a Brent solver with an
absolute accuracy different from the default.
Returns:
the maximum absolute error in inverse cumulative probability estimates/**
     * Returns the solver absolute accuracy for inverse cumulative computation.
     * You can override this method in order to use a Brent solver with an
     * absolute accuracy different from the default.
     *
     * @return the maximum absolute error in inverse cumulative probability estimates
     */
    protected double getSolverAbsoluteAccuracy() {
        return solverAbsoluteAccuracy;
    }

    
{@inheritDoc} /** {@inheritDoc} */
    public void reseedRandomGenerator(long seed) {
        random.setSeed(seed);
        randomData.reSeed(seed);
    }

    
{@inheritDoc}
The default implementation uses the

inversion method.

/**
     * {@inheritDoc}
     *
     * The default implementation uses the
     * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling">
     * inversion method.
     * </a>
     */
    public double sample() {
        return inverseCumulativeProbability(random.nextDouble());
    }

    
{@inheritDoc} The default implementation generates the sample by calling sample() in a loop. /**
     * {@inheritDoc}
     *
     * The default implementation generates the sample by calling
     * {@link #sample()} in a loop.
     */
    public double[] sample(int sampleSize) {
        if (sampleSize <= 0) {
            throw new NotStrictlyPositiveException(LocalizedFormats.NUMBER_OF_SAMPLES,
                    sampleSize);
        }
        double[] out = new double[sampleSize];
        for (int i = 0; i < sampleSize; i++) {
            out[i] = sample();
        }
        return out;
    }

    
{@inheritDoc}
Returns:
zero.
Since:
3.1/**
     * {@inheritDoc}
     *
     * @return zero.
     * @since 3.1
     */
    public double probability(double x) {
        return 0d;
    }

    
Returns the natural logarithm of the probability density function (PDF) of this distribution evaluated at the specified point x. In general, the PDF is the derivative of the CDF. If the derivative does not exist at x, then an appropriate replacement should be returned, e.g. Double.POSITIVE_INFINITY, Double.NaN, or the limit inferior or limit superior of the difference quotient. Note that due to the floating point precision and under/overflow issues, this method will for some distributions be more precise and faster than computing the logarithm of RealDistribution.density(double). The default implementation simply computes the logarithm of density(x). 
Params:
x – the point at which the PDF is evaluated
Returns:
the logarithm of the value of the probability density function at point x/**
     * Returns the natural logarithm of the probability density function (PDF) of this distribution
     * evaluated at the specified point {@code x}. In general, the PDF is the derivative of the
     * {@link #cumulativeProbability(double) CDF}. If the derivative does not exist at {@code x},
     * then an appropriate replacement should be returned, e.g. {@code Double.POSITIVE_INFINITY},
     * {@code Double.NaN}, or the limit inferior or limit superior of the difference quotient. Note
     * that due to the floating point precision and under/overflow issues, this method will for some
     * distributions be more precise and faster than computing the logarithm of
     * {@link #density(double)}. The default implementation simply computes the logarithm of
     * {@code density(x)}.
     *
     * @param x the point at which the PDF is evaluated
     * @return the logarithm of the value of the probability density function at point {@code x}
     */
    public double logDensity(double x) {
        return FastMath.log(density(x));
    }
}