/*
    * 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());
     }
 }