package umontreal.iro.lecuyer.probdistmulti;

import umontreal.iro.lecuyer.util.PrintfFormat;
import umontreal.iro.lecuyer.util.Num;

/**
 * Classes implementing 2-dimensional continuous distributions should inherit 
 * from this class.
 * Such distributions are characterized by a density function f (x, y);
 * thus the signature of a density method is supplied here.
 * This class also provides a default implementation of 
 * bar(F)(x, y), 
 * the upper CDF. The inverse function F-1(u) represents a curve
 * y = h(x) of constant u and it is not implemented.
 * 
 */
public abstract class ContinuousDistribution2Dim
                          extends ContinuousDistributionMulti {

   /**
    * Defines the target number of decimals of accuracy when
    *  approximating a distribution function, but there is no guarantee that
    *  this target is always attained.
    * 
    */
   public int decPrec = 15;
 

    // x infinity for some distributions
     protected static final double XINF = Double.MAX_VALUE;  

    // x infinity for some distributions                                       
    protected static final double XBIG = 1000.0;  

    // EPSARRAY[j]: Epsilon required for j decimal degits of precision
    protected static final double[] EPSARRAY = {
    0.5, 0.5E-1, 0.5E-2, 0.5E-3, 0.5E-4, 0.5E-5, 0.5E-6, 0.5E-7, 0.5E-8,
    0.5E-9, 0.5E-10, 0.5E-11, 0.5E-12, 0.5E-13, 0.5E-14, 0.5E-15, 0.5E-16,
    0.5E-17, 0.5E-18, 0.5E-19, 0.5E-20, 0.5E-21, 0.5E-22, 0.5E-23, 0.5E-24,
    0.5E-25, 0.5E-26, 0.5E-27, 0.5E-28, 0.5E-29, 0.5E-30, 0.5E-31, 0.5E-32,
    0.5E-33, 0.5E-34, 0.5E-35
    };

   /**
    * Returns f (x, y), the density of (X, Y) evaluated at (x, y).
    * 
    * @param x value x at which the density is evaluated
    * 
    *    @param y value y at which the density is evaluated
    * 
    *    @return density function evaluated at (x, y)
    * 
    */
   public abstract double density (double x, double y);


   /**
    * Simply calls density (x[0], x[1]).
    * 
    * @param x point 
    * (x[0], x[1]) at which the density is evaluated
    * 
    *    @return density function evaluated at 
    * (x[0], x[1])
    * 
    */
   public double density (double[] x) {
      if (x.length != 2)
         throw new IllegalArgumentException("x must be in dimension 2");

      return density (x[0], x[1]);
   }


   /**
    * .
    * 
    * Computes the distribution function F(x, y): 
    *  
    * 

    * 

    * F(x, y) = P[X <= x, Y <= y] = ∫-∞xds∫-∞ydt f (s, t).
    * 

    *  
    *  @param x value x at which the distribution function is evaluated
    * 
    *     @param y value y at which the distribution function is evaluated
    * 
    *     @return distribution function evaluated at  (x, y)
    *  
    */
   public abstract double cdf (double x, double y); 

   /**
    * .
    * 
    * Computes  the upper cumulative distribution function
    *   
    * bar(F)(x, y):
    *  
    * 

    * 

    * bar(F)(x, y) = P[X >= x, Y >= y] = ∫∞xds∫∞ydt f (s, t).
    * 

    * 
    * @param x value x at which the upper distribution is evaluated
    * 
    *    @param y value y at which the upper distribution is evaluated
    * 
    *    @return upper distribution function evaluated at  (x, y)
    * 
    */
   public double barF (double x, double y)  {
      double u = 1.0 + cdf (x, y) - cdf (XINF, y) - cdf (x, XINF);
      if (u <= 0.0) return 0.0;
      if (u >= 1.0) return 1.0;
      return u;
   }


   /**
    * .
    * 
    * Computes the cumulative probability in the square region
    *  
    * 

    * 

    * P[a1 <= X <= b1, a2 <= Y <= b2] = ∫a1b1dx∫a2b2dy f (x, y).
    * 
 
    *  
    * @param a1 x lower limit of the square
    * 
    *    @param a2 y lower limit of the square
    * 
    *    @param b1 x upper limit of the square
    * 
    *    @param b2 y upper limit of the square
    * 
    *    @return the cumulative probability in the square region
    * 
    */
   public double cdf (double a1, double a2, double b1, double b2)  {
      if (a1 >= b1 || a2 >= b2) return 0.0;
      return cdf (b1, b2) - cdf (a1, b2) - cdf (b1, a2) + cdf(a1, a2);
   }


}