/* *
 * 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.*;



/**
 * Same as {@link GammaProcessPCABridge}, but uses the fast inversion method
 * for the symmetrical beta distribution, proposed by L'Ecuyer and Simard, to accelerate the generation of the beta random variables.
 * This class works only in the case where the number of intervals is a power of
 * 2 and all these intervals are of equal size.
 * 
 */
public class GammaProcessPCASymmetricalBridge extends GammaProcessPCABridge  {



   /**
    * Constructs a new GammaProcessPCASymmetricalBridge 
    * with parameters 
    * μ = mu, 
    * ν = nu and initial 
    * value 
    * S(t0) = s0.
    * The {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream} stream
    * is used in the 
    * {@link umontreal.iro.lecuyer.randvar.GammaGen GammaGen}
    * and in the {@link umontreal.iro.lecuyer.randvar.BetaSymmetricalGen BetaSymmetricalGen}.
    * These two generators use inversion to generate random numbers.  The first
    * uniform random number generated by stream is used for the gamma, and the 
    * other d - 1 for the beta's.
    * 
    */
   public GammaProcessPCASymmetricalBridge (double s0, double mu, double nu,
                                            RandomStream stream)  {
        super (s0, mu, nu,  stream);
   }


   public double[] generatePath (double[] uniform01)  {
    // uniformsV[] of size d+1, but element 0 never used.
        double[] uniformsV = new double[d+1];

    // generate BrownianMotion PCA path
        double[] BMPCApath = BMPCA.generatePath(uniform01);
        int oldIndexL;
        int newIndex;
        int oldIndexR;
        double temp, y;
    // Transform BMPCA points to uniforms using an inverse bridge.
        for (int j = 0; j < 3*(d-1); j+=3) {
            oldIndexL   = BMBridge.wIndexList[j];
            newIndex    = BMBridge.wIndexList[j + 1];
            oldIndexR   = BMBridge.wIndexList[j + 2];

            temp = BMPCApath[newIndex] - BMPCApath[oldIndexL];
            temp -= (BMPCApath[oldIndexR] - BMPCApath[oldIndexL]) *
                                          BMBridge.wMuDt[newIndex];
            temp /= BMBridge.wSqrtDt[newIndex];
            uniformsV[newIndex] = NormalDist.cdf01(temp);
        }
    double dT = BMPCA.t[d] - BMPCA.t[0];
    uniformsV[d] = NormalDist.cdf01( ( BMPCApath[d] - BMPCApath[0] - BMPCA.mu*dT )/
                     ( BMPCA.sigma * Math.sqrt(dT) ) );


    // Reconstruct the bridge for the GammaProcess from the Brownian uniforms.
    // Here we have to hope that the bridge is implemented in the
    // same order for the Brownian and Gamma processes.

        path[0] = x0;
        path[d] = x0 + GammaDist.inverseF(mu2dTOverNu, muOverNu, 10, uniformsV[d]);
        for (int j = 0; j < 3*(d-1); j+=3) {
            oldIndexL   = wIndexList[j];
            newIndex    = wIndexList[j + 1];
            oldIndexR   = wIndexList[j + 2];

            y = BetaSymmetricalDist.inverseF(bMu2dtOverNuL[newIndex], uniformsV[newIndex]);

            path[newIndex] = path[oldIndexL] +
        (path[oldIndexR] - path[oldIndexL]) * y;
        }
        //resetStartProcess();
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   public double[] generatePath()  {
        double[] u = new double[d];
        for(int i =0; i < d; i++)
            u[i] = stream.nextDouble();
        return generatePath(u);
    }


    // code taken from GammaSymmetricalBridge to check time is power of 2,
    // as it is required for the symmetrical bridge.
   protected void init () {
        super.init ();
        if (observationTimesSet) {

            /* Testing to make sure number of observations n = 2^k */
                int x = d;
            int y = 1;
            while (x>1) {
            x = x / 2; 
            y = y * 2;
            }
            if (y != d) throw new IllegalArgumentException
            ( "GammaSymmetricalBridgeProcess:"
                +"Number 'n' of observation times is not a power of 2" );
       }
    }
}