package umontreal.iro.lecuyer.randvarmulti;

import umontreal.iro.lecuyer.probdist.GammaDist;
import umontreal.iro.lecuyer.randvar.RandomVariateGen;
import umontreal.iro.lecuyer.randvar.GammaAcceptanceRejectionGen;
import umontreal.iro.lecuyer.rng.RandomStream;


/**
 * Extends {@link RandomMultivariateGen} for a
 * Dirichlet distribution.  This distribution uses the
 * parameters
 * 
 * 
 * α1,..., αk, and has density
 * 

 * 
 * 

 * f (x1,…, xk) = Γ(α0)∏i=1kxiαi-1 / ∏i=1kΓ(αi))
 * 

 * where 
 * α0 = ∑i=1kαi.
 * 
 * 

 * Here, the successive coordinates of the Dirichlet vector are generated
 * 
 * via the class
 * {@link umontreal.iro.lecuyer.randvar.GammaAcceptanceRejectionGen GammaAcceptanceRejectionGen}
 * in package randvar, using the same stream for all the uniforms.
 * 
 */
public class DirichletGen extends RandomMultivariateGen {
   private GammaAcceptanceRejectionGen[] ggens;


   /**
    * Constructs a new Dirichlet
    *  generator with parameters 
    * αi+1 = alphas[i],
    *  for 
    * i = 0,…, k - 1, and the stream stream.
    * 
    * @param stream the random number stream used to generate uniforms.
    * 
    *    @param alphas the αi parameters of the generated distribution.
    * 
    *    @exception IllegalArgumentException if one αk is negative or 0.
    * 
    *    @exception NullPointerException if any argument is null.
    * 
    */
   public DirichletGen (RandomStream stream, double[] alphas) {
      if (stream == null)
         throw new NullPointerException ("stream is null");
      this.stream = stream;
      dimension = alphas.length;
      ggens = new GammaAcceptanceRejectionGen[alphas.length];
      for (int k = 0; k < alphas.length; k++)
         ggens[k] = new GammaAcceptanceRejectionGen
            (stream, new GammaDist (alphas[k], 1.0/2.0));
   }


   /**
    * Returns the 
    * αi+1 parameter for this
    *  Dirichlet generator.
    * 
    * @param i the index of the parameter.
    * 
    *    @return the value of the parameter.
    *    @exception ArrayIndexOutOfBoundsException if i is
    *     negative or greater than or equal to {@link #getDimension(()) getDimension}.
    * 
    * 
    */
   public double getAlpha (int i) {
      return ((GammaDist)ggens[i].getDistribution()).getAlpha();
   }


   /**
    * Generates a new point from the Dirichlet distribution with
    *  parameters alphas, using the stream stream.
    *  The generated values are placed into p.
    * 
    * @param stream the random number stream used to generate the uniforms.
    * 
    *    @param alphas the αi parameters of the distribution, for
    *     
    * i = 1,…, k.
    * 
    *    @param p the array to be filled with the generated point.
    * 
    * 
    */
   public static void nextPoint (RandomStream stream, double[] alphas,
                                 double[] p) {
      double total = 0;
      for (int i = 0; i < alphas.length; i++) {
         p[i] = GammaAcceptanceRejectionGen.nextDouble
            (stream, stream, alphas[i], 1.0/2.0);
         total += p[i];
      }
      for (int i = 0; i < alphas.length; i++)
         p[i] /= total;
   }


   /**
    * Generates a point from the Dirichlet distribution.
    * 
    * @param p the array to be filled with the generated point.
    * 
    */
   public void nextPoint (double[] p) {
      int n = ggens.length;
      double total = 0;
      for (int i = 0; i < n; i++) {
         p[i] = ggens[i].nextDouble();
         total += p[i];
      }
      for (int i = 0; i < n; i++)
         p[i] /= total;
   }
}