//Source file: D:\\Projects\\src\\yi\\stat\\em\\AlignedEMCumulator.java

package yi.stat.em;

import org.jmat.data.*;
import java.util.*;

public class AlignedGaussianComponent
    extends GaussianComponent {

  public AlignedGaussianComponent() {
    super();
  }

  public double[][] getCovariance() throws GaussianModelSingularException {
    if (m_covarianceIsFreash)
      return m_covariance;

    int data_length = getDataLength();
    double[] means = getMeans();

    // accumulator
    double prob_sum = 0;

    double[][] covariance = new double[data_length][data_length];

    // loop through each vector
    Iterator iter = m_data.iterator();
    while (iter.hasNext()) {
      ProbabilityData pdata = (ProbabilityData) iter.next();
      // the data
      double[] data = pdata.data;
      double p = pdata.probability;
      prob_sum += p;
      // data - mean
      double[] diff = new double[data_length];
      for (int i = 0; i < data_length; i++) {
        diff[i] = data[i] - means[i];
        // (data - mean)^2 * p

        covariance[i][i] += diff[i] * diff[i] * p;
      }
    }

    if (prob_sum == 0)
      throw new GaussianModelSingularException("probability sum = 0");
    // sum [(data - mean)^2 * p] / sum (p)
    for (int i = 0; i < data_length; i++) {
      covariance[i][i] /= prob_sum;
    }

    m_covariance = covariance;

    // fix zero entries
    if (!fixCovarinceMatrix()) {
      throw new GaussianModelSingularException(
          "zero entry in covariance matrix can't be fixed");
    }

    m_covarianceIsFreash = true;
    return m_covariance;
  }

  AbstractMatrix getCovarianceInverse() throws GaussianModelSingularException {
    if (m_covarianceInverseIsFresh)
      return m_covarianceInverse;

    int data_length = getDataLength();

    double[][] covariance = getCovariance();
    double[][] covarianceInverse = new double[data_length][data_length];

    for (int i = 0; i < data_length; i++) {
      if (covariance[i][i] == 0) {
        throw new GaussianModelSingularException("covariance matrix [" + i + "][" + i +
                                           "] = 0");
      }
      covarianceInverse[i][i] = 1 / covariance[i][i];
    }

    m_covarianceInverse = new Matrix(covarianceInverse);
    m_covarianceInverseIsFresh = true;
    return m_covarianceInverse;
  }


  /**
   * Get sqrt(Sigma) which is useful for standandlize the distribution
   * @return
   * @throws GaussianModelSingularException
   */
  public double[][] getSquareRootOfCovariance() throws
      GaussianModelSingularException {

    double[][] cov = getCovariance();
    int data_length = getDataLength();
    double[][] root_covariance = new double[data_length][data_length];

    for (int i = 0; i < data_length; i++) {
      root_covariance[i][i] = Math.sqrt(cov[i][i]);
    }
    return root_covariance;
  }

  /**
   * Standardize the data (so that the mean of distribution is 0 and the
   * dev is 1)
   * @param data
   * @return
   */
  public double[] standardize(double[] data) throws GaussianModelSingularException{
    double[] mean = getMeans();
    double[][] cov = getCovariance();

    // make zero centered
    double[] norm_data = new double[data.length];
    for (int i=0; i<data.length; i++) {
      norm_data[i] = (data[i] - mean[i])/Math.sqrt(cov[i][i]);
    }

    return norm_data;
  }


  public double getDensityDenominator() throws GaussianModelSingularException {
    if (m_densityDenominatorIsFresh)
      return m_densityDenominator;
    int data_length = getDataLength();
    double[][] covariance = getCovariance();
    double covariance_determinant = 1;
    for (int i=0;i<data_length ; i++) {
      covariance_determinant *= covariance[i][i];
    }
    if (covariance_determinant == 0) {
      fixCovarinceMatrix();
      throw new GaussianModelSingularException("Covariance Matrix Determinant == 0");
    }
    if (covariance_determinant < 0) {
      System.err.println("covariance_determinant < 0, use the inverse one");
      covariance_determinant = -covariance_determinant;
    }

    m_densityDenominator = Math.sqrt(Math.pow(2 * Math.PI, data_length) *
                                     covariance_determinant);
    m_densityDenominatorIsFresh = true;
    return m_densityDenominator;
  }

  /**
   * Return the probability of the data given the model represented by the
   * EMCumulator.
   * @param data
   * @return
   */
  public double getDensity(double[] data) throws GaussianModelSingularException {
    int data_length = data.length;

    double[] mean = getMeans();

    double[][] covariance = getCovariance();

    double up = 0;
    for (int i = 0; i < data_length; i++) {
      double diff = data[i] - mean[i];
      up += diff * diff / covariance[i][i];
    }

    up/=2;


    double densityDenominator = getDensityDenominator();
    return Math.exp(-up) / densityDenominator;
  }


}