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

import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.rng.sampling.distribution.ContinuousUniformSampler;
import org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler;

/**
 * Implementation of the log-uniform distribution. This is also known as the reciprocal distribution.
 *
 * <p>The probability density function of \( X \) is:
 *
 * <p>\[ f(x; a, b) = \frac{1}{x \ln \frac b a} \]
 *
 * <p>for \( 0 \lt a \lt b \lt \infty \) and
 * \( x \in [a, b] \).
 *
 * @see <a href="https://en.wikipedia.org/wiki/Reciprocal_distribution">Reciprocal distribution (Wikipedia)</a>
 * @since 1.1
 */
public final class LogUniformDistribution extends AbstractContinuousDistribution {
    /** Lower bound (a) of this distribution (inclusive). */
    private final double lower;
    /** Upper bound (b) of this distribution (exclusive). */
    private final double upper;
    /** log(a). */
    private final double logA;
    /** log(b). */
    private final double logB;
    /** log(b) - log(a). */
    private final double logBmLogA;
    /** log(log(b) - log(a)). */
    private final double logLogBmLogA;

    /**
     * @param lower Lower bound of this distribution (inclusive).
     * @param upper Upper bound of this distribution (inclusive).
     */
    private LogUniformDistribution(double lower,
                                   double upper) {
        this.lower = lower;
        this.upper = upper;
        logA = Math.log(lower);
        logB = Math.log(upper);
        logBmLogA = logB - logA;
        logLogBmLogA = Math.log(logBmLogA);
    }

    /**
     * Creates a log-uniform distribution.
     *
     * @param lower Lower bound of this distribution (inclusive).
     * @param upper Upper bound of this distribution (inclusive).
     * @return the distribution
     * @throws IllegalArgumentException if {@code lower >= upper}; the range between the bounds
     * is not finite; or {@code lower <= 0}
     */
    public static LogUniformDistribution of(double lower,
                                            double upper) {
        if (lower >= upper) {
            throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GTE_HIGH,
                                            lower, upper);
        }
        if (!Double.isFinite(upper - lower)) {
            throw new DistributionException("Range %s is not finite", upper - lower);
        }
        if (lower <= 0) {
            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, lower);
        }
        return new LogUniformDistribution(lower, upper);
    }

    /** {@inheritDoc} */
    @Override
    public double density(double x) {
        if (x < lower || x > upper) {
            return 0;
        }
        return Math.exp(logDensity(x));
    }

    /** {@inheritDoc} */
    @Override
    public double logDensity(double x) {
        if (x < lower || x > upper) {
            return Double.NEGATIVE_INFINITY;
        }
        return -Math.log(x) - logLogBmLogA;
    }

    /** {@inheritDoc} */
    @Override
    public double cumulativeProbability(double x)  {
        if (x <= lower) {
            return 0;
        }
        if (x >= upper) {
            return 1;
        }
        return (Math.log(x) - logA) / logBmLogA;
    }

    /** {@inheritDoc} */
    @Override
    public double survivalProbability(double x) {
        if (x <= lower) {
            return 1;
        }
        if (x >= upper) {
            return 0;
        }
        return (logB - Math.log(x)) / logBmLogA;
    }

    /** {@inheritDoc} */
    @Override
    public double inverseCumulativeProbability(double p) {
        ArgumentUtils.checkProbability(p);
        // Avoid floating-point error at the bounds
        return clipToRange(Math.exp(logA + p * logBmLogA));
    }

    @Override
    public double inverseSurvivalProbability(double p) {
        ArgumentUtils.checkProbability(p);
        // Avoid floating-point error at the bounds
        return clipToRange(Math.exp(logB - p * logBmLogA));
    }

    /**
     * {@inheritDoc}
     *
     * <p>For lower bound \( a \) and upper bound \( b \), the mean is:
     *
     * <p>\[ \frac{b - a}{\ln \frac b a} \]
     */
    @Override
    public double getMean() {
        return (upper - lower) / logBmLogA;
    }

    /**
     * {@inheritDoc}
     *
     * <p>For lower bound \( a \) and upper bound \( b \), the variance is:
     *
     * <p>\[ \frac{b^2 - a^2}{2 \ln \frac b a} - \left( \frac{b - a}{\ln \frac b a} \right)^2 \]
     */
    @Override
    public double getVariance() {
        // Compute u_2 via a stabilising rearrangement:
        // https://docs.scipy.org/doc/scipy/tutorial/stats/continuous_loguniform.html
        final double a = lower;
        final double b = upper;
        final double d = -logBmLogA;
        return (a - b) * (a * (d - 2) + b * (d + 2)) / (2 * d * d);
    }

    /**
     * {@inheritDoc}
     *
     * <p>The lower bound of the support is equal to the lower bound parameter
     * of the distribution.
     */
    @Override
    public double getSupportLowerBound() {
        return lower;
    }

    /**
     * {@inheritDoc}
     *
     * <p>The upper bound of the support is equal to the upper bound parameter
     * of the distribution.
     */
    @Override
    public double getSupportUpperBound() {
        return upper;
    }

    /**
     * Clip the value to the range [lower, upper].
     * This is used to handle floating-point error at the support bound.
     *
     * @param x Value x
     * @return x clipped to the range
     */
    private double clipToRange(double x) {
        return clip(x, lower, upper);
    }

    /**
     * Clip the value to the range [lower, upper].
     *
     * @param x Value x
     * @param lower Lower bound (inclusive)
     * @param upper Upper bound (inclusive)
     * @return x clipped to the range
     */
    private static double clip(double x, double lower, double upper) {
        if (x <= lower) {
            return lower;
        }
        return x < upper ? x : upper;
    }

    /** {@inheritDoc} */
    @Override
    double getMedian() {
        // Overridden for the probability(double, double) method.
        // This is intentionally not a public method.
        // sqrt(ab) avoiding overflow
        return Math.exp(0.5 * (logA + logB));
    }

    /** {@inheritDoc} */
    @Override
    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
        // Exponentiate a uniform distribution sampler of the logarithmic range.
        final SharedStateContinuousSampler s = ContinuousUniformSampler.of(rng, logA, logB);
        return () -> Math.exp(s.sample());
    }
}