package org.opentrafficsim.core.dsol;
   
   import static nl.tudelft.simulation.jstats.distributions.DistNormal.CUMULATIVE_NORMAL_PROBABILITIES;
   
   import java.util.Locale;
   
   import nl.tudelft.simulation.jstats.distributions.DistContinuous;
   import nl.tudelft.simulation.jstats.streams.MersenneTwister;
   import nl.tudelft.simulation.jstats.streams.StreamInterface;
  
  /**
   * The Normal Truncated distribution. For more information on the normal distribution see
   * <a href="http://mathworld.wolfram.com/NormalDistribution.html"> http://mathworld.wolfram.com/NormalDistribution.html </a>
   * <p>
   * This version of the normal distribution uses the numerically approached inverse cumulative distribution.
   * <p>
   * (c) copyright 2002-2018 <a href="http://www.simulation.tudelft.nl">Delft University of Technology </a>, the Netherlands. <br>
   * See for project information <a href="http://www.simulation.tudelft.nl"> www.simulation.tudelft.nl </a> <br>
   * License of use: <a href="http://www.gnu.org/copyleft/lesser.html">Lesser General Public License (LGPL) </a>, no warranty.
   * @author <a href="mailto:a.verbraeck@tudelft.nl"> Alexander Verbraeck </a> <br>
   *         <a href="http://www.transport.citg.tudelft.nl">Wouter Schakel</a>
   */
  public class DistNormalTrunc extends DistContinuous
  {
      /** */
      private static final long serialVersionUID = 1L;
  
      /** mu refers to the mean of the normal distribution. */
      private final double mu;
  
      /** mu refers to the mean of the normal distribution. */
      private final double sigma;
  
      /** minimum x-value of the distribution. */
      private final double min;
  
      /** maximum x-value of the distribution. */
      private final double max;
  
      /** cumulative distribution value of the minimum. */
      private final double cumulMin;
  
      /** cumulative distribution value difference between max and min. */
      private final double cumulDiff;
  
      /** factor on probability density to normalize to 1. */
      private final double probDensFactor;
  
      /**
       * constructs a normal distribution with mu=0 and sigma=1. Errors of various types, e.g., in the impact point of a bomb;
       * quantities that are the sum of a large number of other quantities by the virtue of the central limit theorem.
       * @param stream StreamInterface; the numberstream
       * @param min double; minimum x-value of the distribution
       * @param max double; maximum x-value of the distribution
       */
      public DistNormalTrunc(final StreamInterface stream, final double min, final double max)
      {
          this(stream, 0.0, 1.0, min, max);
      }
  
      /**
       * constructs a normal distribution with mu and sigma.
       * @param stream StreamInterface; the numberstream
       * @param mu double; the medium
       * @param sigma double; the standard deviation
       * @param min double; minimum x-value of the distribution
       * @param max double; maximum x-value of the distribution
       */
      public DistNormalTrunc(final StreamInterface stream, final double mu, final double sigma, final double min,
              final double max)
      {
          super(stream);
          if (max < min)
          {
              throw new IllegalArgumentException("Error Normal Truncated - max < min");
          }
          this.mu = mu;
          this.sigma = sigma;
          this.min = min;
          this.max = max;
          this.cumulMin = getCumulativeProbabilityNotTruncated(min);
          this.cumulDiff = getCumulativeProbabilityNotTruncated(max) - this.cumulMin;
          this.probDensFactor = 1.0 / this.cumulDiff;
      }
  
      /** {@inheritDoc} */
      @Override
      public double draw()
      {
          return getInverseCumulativeProbabilityNotTruncated(this.cumulMin + this.cumulDiff * this.stream.nextDouble());
      }
  
      /**
       * returns the cumulative probability of the x-value.
       * @param x double; the observation x
       * @return double the cumulative probability
       */
      public double getCumulativeProbability(final double x)
      {
         if (x < this.min)
         {
             return 0.0;
         }
         if (x > this.max)
         {
             return 1.0;
         }
         return (getCumulativeProbabilityNotTruncated(x) - this.cumulMin) * this.probDensFactor;
     }
 
     /**
      * returns the cumulative probability of the x-value.
      * @param x double; the observation x
      * @return double the cumulative probability
      */
     private double getCumulativeProbabilityNotTruncated(final double x)
     {
         double z = (x - this.mu) / this.sigma * 100;
         double absZ = Math.abs(z);
         int intZ = (int) absZ;
         double f = 0.0;
         if (intZ >= 1000)
         {
             intZ = 999;
             f = 1.0;
         }
         else
         {
             f = absZ - intZ;
         }
         if (z >= 0)
         {
             return (1 - f) * CUMULATIVE_NORMAL_PROBABILITIES[intZ] + f * CUMULATIVE_NORMAL_PROBABILITIES[intZ + 1];
         }
         return 1 - ((1 - f) * CUMULATIVE_NORMAL_PROBABILITIES[intZ] + f * CUMULATIVE_NORMAL_PROBABILITIES[intZ + 1]);
     }
 
     /**
      * returns the x-value of the given cumulativePropability.
      * @param cumulativeProbability double; reflects cum prob
      * @return double the inverse cumulative probability
      */
     public double getInverseCumulativeProbability(final double cumulativeProbability)
     {
         if (cumulativeProbability < 0 || cumulativeProbability > 1)
         {
             throw new IllegalArgumentException("1<cumulativeProbability<0 ?");
         }
         /*
          * For extreme cases we return the min and max directly. The method getInverseCumulativeProbabilityNotTruncated() can
          * only return values from "mu - 10*sigma" to "mu + 10*sigma". If min or max is beyond these values, those values would
          * result erroneously. For any cumulative probability that is slightly above 0.0 or slightly below 1.0, values in the
          * range from "mu - 10*sigma" to "mu + 10*sigma" will always result.
          */
         if (cumulativeProbability == 0.0)
         {
             return this.min;
         }
         if (cumulativeProbability == 1.0)
         {
             return this.max;
         }
         return getInverseCumulativeProbabilityNotTruncated(this.cumulMin + cumulativeProbability * this.cumulDiff);
     }
 
     /** {@inheritDoc} */
     @Override
     public double probDensity(final double x)
     {
         if (x < this.min || x > this.max)
         {
             return 0.0;
         }
         return this.probDensFactor / (Math.sqrt(2 * Math.PI * Math.pow(this.sigma, 2)))
                 * Math.exp(-1 * Math.pow(x - this.mu, 2) / (2 * Math.pow(this.sigma, 2)));
     }
 
     /** {@inheritDoc} */
     @Override
     public String toString()
     {
         return "NormalTrunc(" + this.mu + "," + this.sigma + "," + this.min + "," + this.max + ")";
     }
 
     /**
      * Test.
      * @param args String[]; args
      */
     public static void main(final String[] args)
     {
         StreamInterface stream = new MersenneTwister();
         double mu = 2.0;
         double sigma = 3.0;
         double min = -5.0;
         double max = 4.0;
         DistNormalTrunc dist = new DistNormalTrunc(stream, mu, sigma, min, max);
 
         System.out.println("<< probability density >>");
         double sum = 0.0;
         double step = (max - min) / 96;
         for (double x = min - 2 * step; x <= max + 2 * step; x += step)
         {
             double p = dist.probDensity(x);
             System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", x, p));
             sum += p * step;
         }
         System.out.println(String.format(Locale.GERMAN, "Approx. sum = %.8f", sum));
         System.out.println("");
 
         System.out.println("<< cumulative density >>");
         for (double x = min - 2 * step; x <= max + 2 * step; x += step)
         {
             double c = dist.getCumulativeProbability(x);
             System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", x, c));
         }
         System.out.println("");
 
         System.out.println("<< inverse cumulative density >>");
         for (double c = 0.0; c < 1.005; c += 0.01)
         {
             double x = dist.getInverseCumulativeProbability(Math.min(c, 1.0)); // want to include 1.0, also if 1.0000000000001
             System.out.println(String.format(Locale.GERMAN, "%.8f;%.8f", c, x));
         }
         System.out.println("");
 
         System.out.println("<< 10000 random numbers. >>");
         for (int i = 1; i < 10000; i++)
         {
             System.out.println(String.format(Locale.GERMAN, "%,8f", dist.draw()));
         }
 
     }
 
     /**
      * returns the x-value of the given cumulativePropability.
      * @param cumulativeProbability double; reflects cum prob
      * @return double the inverse cumulative probability
      */
     private double getInverseCumulativeProbabilityNotTruncated(final double cumulativeProbability)
     {
         if (cumulativeProbability < 0 || cumulativeProbability > 1)
         {
             throw new IllegalArgumentException("1<cumulativeProbability<0 ?");
         }
         boolean located = false;
         double prob = cumulativeProbability;
         if (cumulativeProbability < 0.5)
         {
             prob = 1 - cumulativeProbability;
         }
         int i = 0;
         double f = 0.0;
         while (!located)
         {
             if (CUMULATIVE_NORMAL_PROBABILITIES[i] < prob && CUMULATIVE_NORMAL_PROBABILITIES[i + 1] >= prob)
             {
                 located = true;
                 if (CUMULATIVE_NORMAL_PROBABILITIES[i] < CUMULATIVE_NORMAL_PROBABILITIES[i + 1])
                 {
                     f = (prob - CUMULATIVE_NORMAL_PROBABILITIES[i])
                             / (CUMULATIVE_NORMAL_PROBABILITIES[i + 1] - CUMULATIVE_NORMAL_PROBABILITIES[i]);
                 }
             }
             else
             {
                 i++;
             }
         }
         if (cumulativeProbability < 0.5)
         {
             return this.mu - ((f + i) / 100.0) * this.sigma;
         }
         return ((f + i) / 100.0) * this.sigma + this.mu;
     }
 
 }