/*
    * 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;
  
  /**
   * Implementation of the trapezoidal distribution.
   *
   * <p>The probability density function of \( X \) is:
   *
   * <p>\[ f(x; a, b, c, d) = \begin{cases}
   *       \frac{2}{d+c-a-b}\frac{x-a}{b-a} &amp; \text{for } a\le x \lt b \\
   *       \frac{2}{d+c-a-b}                &amp; \text{for } b\le x \lt c \\
   *       \frac{2}{d+c-a-b}\frac{d-x}{d-c} &amp; \text{for } c\le x \le d
   *       \end{cases} \]
   *
   * <p>for \( -\infty \lt a \le b \le c \le d \lt \infty \) and
   * \( x \in [a, d] \).
   *
   * <p>Note the special cases:
   * <ul>
   * <li>\( b = c \) is the triangular distribution
   * <li>\( a = b \) and \( c = d \) is the uniform distribution
   * </ul>
   *
   * @see <a href="https://en.wikipedia.org/wiki/Trapezoidal_distribution">Trapezoidal distribution (Wikipedia)</a>
   */
  public abstract class TrapezoidalDistribution extends AbstractContinuousDistribution {
      /** Lower limit of this distribution (inclusive). */
      protected final double a;
      /** Start of the trapezoid constant density. */
      protected final double b;
      /** End of the trapezoid constant density. */
      protected final double c;
      /** Upper limit of this distribution (inclusive). */
      protected final double d;
  
      /**
       * Specialisation of the trapezoidal distribution used when the distribution simplifies
       * to an alternative distribution.
       */
      private static class DelegatedTrapezoidalDistribution extends TrapezoidalDistribution {
          /** Distribution delegate. */
          private final ContinuousDistribution delegate;
  
          /**
           * @param a Lower limit of this distribution (inclusive).
           * @param b Start of the trapezoid constant density.
           * @param c End of the trapezoid constant density.
           * @param d Upper limit of this distribution (inclusive).
           * @param delegate Distribution delegate.
           */
          DelegatedTrapezoidalDistribution(double a, double b, double c, double d,
                                           ContinuousDistribution delegate) {
              super(a, b, c, d);
              this.delegate = delegate;
          }
  
          @Override
          public double density(double x) {
              return delegate.density(x);
          }
  
          @Override
          public double probability(double x0, double x1) {
              return delegate.probability(x0, x1);
          }
  
          @Override
          public double logDensity(double x) {
              return delegate.logDensity(x);
          }
  
          @Override
          public double cumulativeProbability(double x) {
              return delegate.cumulativeProbability(x);
          }
  
          @Override
          public double inverseCumulativeProbability(double p) {
              return delegate.inverseCumulativeProbability(p);
          }
  
         @Override
         public double survivalProbability(double x) {
             return delegate.survivalProbability(x);
         }
 
         @Override
         public double inverseSurvivalProbability(double p) {
             return delegate.inverseSurvivalProbability(p);
         }
 
         @Override
         public double getMean() {
             return delegate.getMean();
         }
 
         @Override
         public double getVariance() {
             return delegate.getVariance();
         }
 
         @Override
         public Sampler createSampler(UniformRandomProvider rng) {
             return delegate.createSampler(rng);
         }
     }
 
     /**
      * Specialisation of the trapezoidal distribution used when {@code b == c}.
      *
      * <p>This delegates all methods to the triangular distribution.
      */
     private static class TriangularTrapezoidalDistribution extends DelegatedTrapezoidalDistribution {
         /**
          * @param a Lower limit of this distribution (inclusive).
          * @param b Start/end of the trapezoid constant density (mode).
          * @param d Upper limit of this distribution (inclusive).
          */
         TriangularTrapezoidalDistribution(double a, double b, double d) {
             super(a, b, b, d, TriangularDistribution.of(a, b, d));
         }
     }
 
     /**
      * Specialisation of the trapezoidal distribution used when {@code a == b} and {@code c == d}.
      *
      * <p>This delegates all methods to the uniform distribution.
      */
     private static class UniformTrapezoidalDistribution extends DelegatedTrapezoidalDistribution {
         /**
          * @param a Lower limit of this distribution (inclusive).
          * @param d Upper limit of this distribution (inclusive).
          */
         UniformTrapezoidalDistribution(double a, double d) {
             super(a, a, d, d, UniformContinuousDistribution.of(a, d));
         }
     }
 
     /**
      * Regular implementation of the trapezoidal distribution.
      */
     private static class RegularTrapezoidalDistribution extends TrapezoidalDistribution {
         /** Cached value (d + c - a - b). */
         private final double divisor;
         /** Cached value (b - a). */
         private final double bma;
         /** Cached value (d - c). */
         private final double dmc;
         /** Cumulative probability at b. */
         private final double cdfB;
         /** Cumulative probability at c. */
         private final double cdfC;
         /** Survival probability at b. */
         private final double sfB;
         /** Survival probability at c. */
         private final double sfC;
 
         /**
          * @param a Lower limit of this distribution (inclusive).
          * @param b Start of the trapezoid constant density.
          * @param c End of the trapezoid constant density.
          * @param d Upper limit of this distribution (inclusive).
          */
         RegularTrapezoidalDistribution(double a, double b, double c, double d) {
             super(a, b, c, d);
 
             // Sum positive terms
             divisor = (d - a) + (c - b);
             bma = b - a;
             dmc = d - c;
 
             cdfB = bma / divisor;
             sfB = 1 - cdfB;
             sfC = dmc / divisor;
             cdfC = 1 - sfC;
         }
 
         @Override
         public double density(double x) {
             // Note: x < a allows correct density where a == b
             if (x < a) {
                 return 0;
             }
             if (x < b) {
                 final double divident = (x - a) / bma;
                 return 2 * (divident / divisor);
             }
             if (x < c) {
                 return 2 / divisor;
             }
             if (x < d) {
                 final double divident = (d - x) / dmc;
                 return 2 * (divident / divisor);
             }
             return 0;
         }
 
         @Override
         public double cumulativeProbability(double x)  {
             if (x <= a) {
                 return 0;
             }
             if (x < b) {
                 final double divident = (x - a) * (x - a) / bma;
                 return divident / divisor;
             }
             if (x < c) {
                 final double divident = 2 * x - b - a;
                 return divident / divisor;
             }
             if (x < d) {
                 final double divident = (d - x) * (d - x) / dmc;
                 return 1 - divident / divisor;
             }
             return 1;
         }
 
         @Override
         public double survivalProbability(double x)  {
             // By symmetry:
             if (x <= a) {
                 return 1;
             }
             if (x < b) {
                 final double divident = (x - a) * (x - a) / bma;
                 return 1 - divident / divisor;
             }
             if (x < c) {
                 final double divident = 2 * x - b - a;
                 return 1 - divident / divisor;
             }
             if (x < d) {
                 final double divident = (d - x) * (d - x) / dmc;
                 return divident / divisor;
             }
             return 0;
         }
 
         @Override
         public double inverseCumulativeProbability(double p) {
             ArgumentUtils.checkProbability(p);
             if (p == 0) {
                 return a;
             }
             if (p == 1) {
                 return d;
             }
             if (p < cdfB) {
                 return a + Math.sqrt(p * divisor * bma);
             }
             if (p < cdfC) {
                 return 0.5 * ((p * divisor) + a + b);
             }
             return d - Math.sqrt((1 - p) * divisor * dmc);
         }
 
         @Override
         public double inverseSurvivalProbability(double p) {
             // By symmetry:
             ArgumentUtils.checkProbability(p);
             if (p == 1) {
                 return a;
             }
             if (p == 0) {
                 return d;
             }
             if (p > sfB) {
                 return a + Math.sqrt((1 - p) * divisor * bma);
             }
             if (p > sfC) {
                 return 0.5 * (((1 - p) * divisor) + a + b);
             }
             return d - Math.sqrt(p * divisor * dmc);
         }
 
         @Override
         public double getMean() {
             // Compute using a standardized distribution
             // b' = (b-a) / (d-a)
             // c' = (c-a) / (d-a)
             final double scale = d - a;
             final double bp = bma / scale;
             final double cp = (c - a) / scale;
             return nonCentralMoment(1, bp, cp) * scale + a;
         }
 
         @Override
         public double getVariance() {
             // Compute using a standardized distribution
             // b' = (b-a) / (d-a)
             // c' = (c-a) / (d-a)
             final double scale = d - a;
             final double bp = bma / scale;
             final double cp = (c - a) / scale;
             final double mu = nonCentralMoment(1, bp, cp);
             return (nonCentralMoment(2, bp, cp) - mu * mu) * scale * scale;
         }
 
         /**
          * Compute the {@code k}-th non-central moment of the standardized trapezoidal
          * distribution.
          *
          * <p>Shifting the distribution by scale {@code (d - a)} and location {@code a}
          * creates a standardized trapezoidal distribution. This has a simplified
          * non-central moment as {@code a = 0, d = 1, 0 <= b < c <= 1}.
          * <pre>
          *               2             1       ( 1 - c^(k+2)           )
          * E[X^k] = ----------- -------------- ( ----------- - b^(k+1) )
          *          (1 + c - b) (k + 1)(k + 2) (    1 - c              )
          * </pre>
          *
          * <p>Simplification eliminates issues computing the moments when {@code a == b}
          * or {@code c == d} in the original (non-standardized) distribution.
          *
          * @param k Moment to compute
          * @param b Start of the trapezoid constant density (shape parameter in [0, 1]).
          * @param c End of the trapezoid constant density (shape parameter in [0, 1]).
          * @return the moment
          */
         private static double nonCentralMoment(int k, double b, double c) {
             // As c -> 1 then (1 - c^(k+2)) loses precision
             // 1 - x^y == -(x^y - 1)    [high precision powm1]
             //         == -(exp(y * log(x)) - 1)
             // Note: avoid log(1) using the limit:
             // (1 - c^(k+2)) / (1-c) -> (k+2) as c -> 1
             final double term1 = c == 1 ? k + 2 : Math.expm1((k + 2) * Math.log(c)) / (c - 1);
             final double term2 = Math.pow(b, k + 1);
             return 2 * ((term1 - term2) / (c - b + 1) / ((k + 1) * (k + 2)));
         }
     }
 
     /**
      * @param a Lower limit of this distribution (inclusive).
      * @param b Start of the trapezoid constant density.
      * @param c End of the trapezoid constant density.
      * @param d Upper limit of this distribution (inclusive).
      */
     TrapezoidalDistribution(double a, double b, double c, double d) {
         this.a = a;
         this.b = b;
         this.c = c;
         this.d = d;
     }
 
     /**
      * Creates a trapezoidal distribution.
      *
      * <p>The distribution density is represented as an up sloping line from
      * {@code a} to {@code b}, constant from {@code b} to {@code c}, and then a down
      * sloping line from {@code c} to {@code d}.
      *
      * @param a Lower limit of this distribution (inclusive).
      * @param b Start of the trapezoid constant density (first shape parameter).
      * @param c End of the trapezoid constant density (second shape parameter).
      * @param d Upper limit of this distribution (inclusive).
      * @return the distribution
      * @throws IllegalArgumentException if {@code a >= d}, if {@code b < a}, if
      * {@code c < b} or if {@code c > d}.
      */
     public static TrapezoidalDistribution of(double a, double b, double c, double d) {
         if (a >= d) {
             throw new DistributionException(DistributionException.INVALID_RANGE_LOW_GTE_HIGH,
                                             a, d);
         }
         if (b < a) {
             throw new DistributionException(DistributionException.TOO_SMALL,
                                             b, a);
         }
         if (c < b) {
             throw new DistributionException(DistributionException.TOO_SMALL,
                                             c, b);
         }
         if (c > d) {
             throw new DistributionException(DistributionException.TOO_LARGE,
                                             c, d);
         }
         // For consistency, delegate to the appropriate simplified distribution.
         // Note: Floating-point equality comparison is intentional.
         if (b == c) {
             return new TriangularTrapezoidalDistribution(a, b, d);
         }
         if (d - a == c - b) {
             return new UniformTrapezoidalDistribution(a, d);
         }
         return new RegularTrapezoidalDistribution(a, b, c, d);
     }
 
     /**
      * {@inheritDoc}
      *
      * <p>For lower limit \( a \), start of the density constant region \( b \),
      * end of the density constant region \( c \) and upper limit \( d \), the
      * mean is:
      *
      * <p>\[ \frac{1}{3(d+c-b-a)}\left(\frac{d^3-c^3}{d-c}-\frac{b^3-a^3}{b-a}\right) \]
      */
     @Override
     public abstract double getMean();
 
     /**
      * {@inheritDoc}
      *
      * <p>For lower limit \( a \), start of the density constant region \( b \),
      * end of the density constant region \( c \) and upper limit \( d \), the
      * variance is:
      *
      * <p>\[ \frac{1}{6(d+c-b-a)}\left(\frac{d^4-c^4}{d-c}-\frac{b^4-a^4}{b-a}\right) - \mu^2 \]
      *
      * <p>where \( \mu \) is the mean.
      */
     @Override
     public abstract double getVariance();
 
     /**
      * Gets the start of the constant region of the density function.
      *
      * <p>This is the first shape parameter {@code b} of the distribution.
      *
      * @return the first shape parameter {@code b}
      */
     public double getB() {
         return b;
     }
 
     /**
      * Gets the end of the constant region of the density function.
      *
      * <p>This is the second shape parameter {@code c} of the distribution.
      *
      * @return the second shape parameter {@code c}
      */
     public double getC() {
         return c;
     }
 
     /**
      * {@inheritDoc}
      *
      * <p>The lower bound of the support is equal to the lower limit parameter
      * {@code a} of the distribution.
      */
     @Override
     public double getSupportLowerBound() {
         return a;
     }
 
     /**
      * {@inheritDoc}
      *
      * <p>The upper bound of the support is equal to the upper limit parameter
      * {@code d} of the distribution.
      */
     @Override
     public double getSupportUpperBound() {
         return d;
     }
 }