/* *
 * Title:        

 * Description:  

 * Copyright:    Copyright (c) 

 * Company:      

 * @author
 * @version 1.0
 */


package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * The inverse Gaussian process is a non-decreasing process
 * where the increments are additive and are given by the
 * inverse gaussian distribution,
 * {@link umontreal.iro.lecuyer.probdist.InverseGaussianDist InverseGaussianDist}.
 * With parameters δ and γ, the
 * time increments are given by
 * {@link umontreal.iro.lecuyer.probdist.InverseGaussianDist InverseGaussianDist}
 * (δdt/γ, δ2dt2).
 * 
 * 

 * [We here use the inverse gaussian distribution
 * parametrized with IGDist
 * (μ, λ), where 
 * μ = δ/γ
 * and 
 * λ = δ2.  If we instead used the parametrization
 * 
 * IGDist(δ, γ),
 * then the increment distribution of our process would have been written
 * more simply as 
 * IGDist[tex2html_wrap_inline158](δdt, γ).]
 * 
 * 

 * The increments are generated by using
 * the inversion of the cumulative distribution function.
 * It therefore uses only one {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}.
 * Subclasses of this class use different generating methods and some need
 * two {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}'s.
 * 
 * 

 * The initial value of this process is the initial observation time.
 * 
 */
public class InverseGaussianProcess extends StochasticProcess  {

    protected RandomStream       stream;

    protected double   delta;
    protected double   gamma;

    protected double   deltaOverGamma;
    protected double   deltaSquare;
    // mu and lambda are the common names of the params for InverseGaussianGen.
    protected double[] imu;
    protected double[] ilam;

    // Number of random streams needed by the current class
    // to generate an IG process.  For this class, = 1, for subclasses
    // will sometimes be, = 2.
    int numberOfRandomStreams;



    // needed by InverseGaussianProcessMSH
   protected InverseGaussianProcess()  { }


   /**
    * Constructs a new InverseGaussianProcess.
    * The initial value s0 will be overridden by t[0] when
    * the observation times are set.
    * 
    */
   public InverseGaussianProcess (double s0, double delta, double gamma,
                                  RandomStream stream)  {
        this.x0 = s0;
        setParams(delta, gamma);
        this.stream = stream;
        numberOfRandomStreams = 1;
    }


   public double[] generatePath() {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianGen.nextDouble(stream, imu[i], ilam[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * Instead of using the internal stream to generate the path,
    * uses an array of uniforms U[0, 1).  The array should be
    * of the length of the number of periods in the observation
    * times. This method is useful for {@link NormalInverseGaussianProcess}.
    * 
    */
   public double[] generatePath (double[] uniforms01)  {
        double s = x0;
        for (int i = 0; i < d; i++) {
            s += InverseGaussianDist.inverseF(imu[i], ilam[i], uniforms01[i]);
            path[i+1] = s;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   /**
    * This method does not work for this class, but will
    * be useful for the subclasses that require two streams.
    * 
    */
   public double[] generatePath (double[] uniforms01, double[] uniforms01b)  {
       throw new UnsupportedOperationException("Use generatePath with 1 stream.");
    }


   public double nextObservation () {
        double s = path[observationIndex];
        s += InverseGaussianGen.nextDouble(stream,
                  imu[observationIndex], ilam[observationIndex]);
        observationIndex++;
        observationCounter = observationIndex;
        path[observationIndex] = s;
        return s;
    }

   /**
    * Sets the parameters.
    * 
    */
   public void setParams (double delta, double gamma)  {
        this.delta = delta;
        this.gamma = gamma;
        deltaOverGamma = delta/gamma;
        deltaSquare    = delta*delta;
        init();
    }


   /**
    * Returns δ.
    * 
    */
   public double getDelta()  {
        return delta;
    }


   /**
    * Returns γ.
    * 
    */
   public double getGamma()  {
        return gamma;
    }


   /**
    * Returns the analytic average which is 
    * δt/γ,
    *    with t = time.
    * 
    */
   public double getAnalyticAverage (double time)  {
        return delta*time/gamma;
    }


   /**
    * Returns the analytic variance which is 
    * (δt)2,
    *    with t = time.
    * 
    */
   public double getAnalyticVariance (double time)  {
        return delta*delta*time*time;
    }


    protected void init () {
        super.init(); // set path[0] to s0
        if(observationTimesSet){
            x0 = t[0];
            path[0] = x0;
            imu  = new double[d];
            ilam = new double[d];
            for (int j = 0; j < d; j++)
            {
                double temp = delta * (t[j+1] - t[j]);
                imu[j]  = temp / gamma;
                ilam[j] = temp * temp;
            }
        }
    }

   public RandomStream getStream () {
        // It is assumed that stream is always the same
        // as the stream in the InverseGaussianGen.
        return stream;
    }

   public void setStream (RandomStream stream) {
       this.stream = stream;
    }

   /**
    * Returns the number of random streams of this process.
    * It is useful because some subclasses use different number
    * of streams.  It returns 1 for {@link InverseGaussianProcess}.
    * 
    */
   public int getNumberOfRandomStreams()  {
        return numberOfRandomStreams;
    }


}