/* *
 * 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.*;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.*;



/**
 * A Brownian motion process 
 * {X(t) : t >= 0} sampled using the
 * principal component decomposition (PCA).
 * 
 */
public class BrownianMotionPCA extends BrownianMotion  {

    protected double[][]  sigmaCov; // Matrice de covariance du vecteur des observ.
                                 // sigmaCov [i][j] = Cov[X(t_{i+1}),X(t_{j+1})].
    protected double[][]  A;     // sigmaCov = AA' (PCA decomposition).
    protected double[]    z;     // vector of standard normals.
    protected double[]    sortedEigenvalues;
    protected boolean     isDecompPCA;




   /**
    * Constructs a new BrownianMotionBridge with
    * parameters 
    * μ = mu, 
    * σ = sigma and initial value
    * 
    * X(t0) = x0.
    * The normal variates will be generated by inversion using stream.
    * 
    */
   public BrownianMotionPCA (double x0, double mu, double sigma,
                             RandomStream stream)  {
        super (x0, mu, sigma, stream);
        isDecompPCA = false;
    }


   /**
    * Constructs a new BrownianMotionBridge with
    * parameters 
    * μ = mu, 
    * σ = sigma and initial value
    * 
    * X(t0) = x0.
    * The normal variates will be generated by gen.
    * 
    */
   public BrownianMotionPCA (double x0, double mu, double sigma,
                             NormalGen gen)  {
        super (x0, mu, sigma, gen);
        isDecompPCA = false;
    }


   public double nextObservation() {
        throw new UnsupportedOperationException ("nextObservation() is not defined in BrownianMotionPCA");
    }

   public void setParams (double x0, double mu, double sigma) {
        super.setParams(x0, mu, sigma);
        isDecompPCA = true;  // the setParams method defined in BrownianMotion calls init() which does
                             // the PCA decompostion
    }

   public double[] generatePath() {
       if(!isDecompPCA) {init();}  // if the decomposition is not done, do it...
       for (int j = 0; j < d; j++)
           z[j] = gen.nextDouble ();
       for (int j = 0; j < d; j++) {
           double sum = 0.0;
           for (int k = 0; k < d; k++)
               sum += A[j][k] * z[k];
           path[j+1] = x0 + mu * t[j+1] + sum;
       }
       observationIndex   = d;
       observationCounter = d;
       return path;
    }

   public double[] generatePath (double[] uniform01) {
       if(!isDecompPCA) {init();}  // if the decomposition is not done, do it...
       for (int j = 0; j < d; j++)
           z[j] = NormalDist.inverseF01(uniform01[j]);
       for (int j = 0; j < d; j++) {
           double sum = 0.0;
           for (int k = 0; k < d; k++)
               sum += A[j][k] * z[k];
           path[j+1] = x0 + mu * t[j+1] + sum;
       }
       observationIndex   = d;
       observationCounter = d;
       return path;
    }

   public double[][] decompPCA (double[][] sigma){
      // L'objet SingularValueDecomposition permet de recuperer la matrice
      // des valeurs propres en ordre decroissant et celle des vecteurs propres de
      // sigma (pour une matrice symetrique et definie-positive seulement).
      SingularValueDecomposition sv =
         new SingularValueDecomposition (new DenseDoubleMatrix2D (sigma));
      DoubleMatrix2D D = sv.getS ();    // diagonal
      // Calculer la racine carree des valeurs propres
      for (int i = 0; i < D.rows(); i++){
         sortedEigenvalues[i] = D.getQuick (i, i);
         D.setQuick (i, i, Math.sqrt (D.getQuick (i, i)));
      }
      DoubleMatrix2D P = sv.getV();   // right factor matrix
      return P.zMult (D, null).toArray ();
   }

   /**
    * Returns the sorted eigenvalues obtained in the PCA decomposition.
    * 
    */
   public double[] getSortedEigenvalues()  {
        return sortedEigenvalues;
   }



   protected void init() {
        super.init();
        z = new double[d];
        sortedEigenvalues = new double[d];

        // Initialize sigmaCov, based on the observation times.
        sigmaCov = new double[d][d];
        for (int i = 0; i < d; i++) {
           for (int j = i; j < d; j++) {
              sigmaCov[i][j] = sigma * sigma * t[i+1];
              sigmaCov[j][i] = sigmaCov[i][j];
           }
        }
        A = decompPCA (sigmaCov);
        isDecompPCA = true;
    }
}