/* *
 * 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 GammaProcessPCA}, but the generated uniforms
 * correspond to a bridge transformation of the {@link BrownianMotionPCA}
 * instead of a sequential transformation.
 * 
 */
public class GammaProcessPCABridge extends GammaProcessPCA  {
    protected BrownianMotionBridge BMBridge;
    protected double       mu2OverNu,
                           mu2dTOverNu;
    protected double[]     bMu2dtOverNuL,  // For precomputations for G Bridge
                           bMu2dtOverNuR;
    protected int[]        wIndexList;



   /**
    * Constructs a new GammaProcessPCABridge with parameters
    * 
    * μ = mu, 
    * ν = nu and initial value 
    * S(t0) = s0.
    * The same stream is used to generate the gamma and beta random numbers.  All
    * these numbers are generated by inversion in the following order:
    * the first uniform random number generated
    * is used for the gamma and the other d - 1 for the beta's.
    * 
    */
   public GammaProcessPCABridge (double s0, double mu, double nu, 
                                 RandomStream stream) {
        super (s0, mu, nu,  stream);
        this.BMBridge = new BrownianMotionBridge(0.0, 0.0, Math.sqrt(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.
        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 =  BetaDist.inverseF(bMu2dtOverNuL[newIndex],  bMu2dtOverNuR[newIndex], 8, uniformsV[newIndex]);

            path[newIndex] = path[oldIndexL] +
		(path[oldIndexR] - path[oldIndexL]) * y;
        }
        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);
    }


    public void setParams (double s0, double mu, double nu) {
        super.setParams(s0, mu, nu);
        BMBridge.setParams(0.0, 0.0, Math.sqrt(nu));
        BMPCA.setParams(0.0, 0.0, Math.sqrt(nu));
    }


    public void setObservationTimes (double[] t, int d) {
        super.setObservationTimes(t, d);
        BMBridge.setObservationTimes(t, d);
    }

   /**
    * Returns the inner {@link BrownianMotionPCA}.
    * 
    */
    public BrownianMotionPCA getBMPCA()  {
        return BMPCA;
    }


    protected void init() {
        super.init();
        if (observationTimesSet) {

        // Quantities for gamma bridge process
        bMu2dtOverNuL = new double[d+1];
        bMu2dtOverNuR = new double[d+1];
        wIndexList  = new int[3*d];

        int[] ptIndex = new int[d+1];
        int   indexCounter = 0;
        int   newIndex, oldLeft, oldRight;

        ptIndex[0] = 0;
        ptIndex[1] = d;

        mu2OverNu   = mu * mu / nu;
        mu2dTOverNu = mu2OverNu * (t[d] - t[0]);

        for (int powOfTwo = 1; powOfTwo <= d/2; powOfTwo *= 2) {
            /* Make room in the indexing array "ptIndex" */
            for (int j = powOfTwo; j >= 1; j--) { ptIndex[2*j] = ptIndex[j]; }

            /* Insert new indices and Calculate constants */
            for (int j = 1; j <= powOfTwo; j++) {
                oldLeft  = 2*j - 2;
                oldRight = 2*j;
                newIndex = (int) (0.5*(ptIndex[oldLeft] + ptIndex[oldRight]));

                bMu2dtOverNuL[newIndex] = mu * mu
                                   * (t[newIndex] - t[ptIndex[oldLeft]]) / nu;
                bMu2dtOverNuR[newIndex] = mu * mu 
                                  * (t[ptIndex[oldRight]] - t[newIndex]) / nu;

                ptIndex[oldLeft + 1]       = newIndex;
                wIndexList[indexCounter]   = ptIndex[oldLeft];
                wIndexList[indexCounter+1] = newIndex;
                wIndexList[indexCounter+2] = ptIndex[oldRight];

                indexCounter += 3;
            }
        }
        /* Check if there are holes remaining and fill them */
        for (int k = 1; k < d; k++) {
            if (ptIndex[k-1] + 1 < ptIndex[k]) {
            // there is a hole between (k-1) and k.

                bMu2dtOverNuL[ptIndex[k-1]+1] = mu * mu 
                                  * (t[ptIndex[k-1]+1] - t[ptIndex[k-1]]) / nu;
                bMu2dtOverNuR[ptIndex[k-1]+1] = mu * mu 
                                  * (t[ptIndex[k]] - t[ptIndex[k-1]+1]) / nu;

                wIndexList[indexCounter]   = ptIndex[k]-2;
                wIndexList[indexCounter+1] = ptIndex[k]-1;
                wIndexList[indexCounter+2] = ptIndex[k];
                indexCounter += 3;
            }
        }
        }
    }

}