/*  
 *  src/org/metaqtl/EMResult.java
 * 
 *  $Author: veyrieras $, $Date: 2006-10-17 15:41:54 $, $Version$
 *
 *  MetaQTL : a Java package for QTL Meta-Analysis
 * 
 *  Copyright (C) 2003-2006  Jean-Baptiste Veyrieras, INRA, France.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA * 
 */
package org.metaqtl;

import org.metaqtl.numrec.NumericalUtilities;

/**
 * 
 */
public class EMResult {
	/**
	 * The number of clusters.
	 */
	public int k;
	/**
	 * The number of data points.
	 */
	public int n;
	/**
	 * The quadratric rate of convergence
	 * of the algorithm.
	 */
	public double rate;
	/**
	 * The euclidean distance between the 
	 * current and the previous vector of
	 * parameters (mu + pi) 
	 */
	public double edist;
	/**
	 * The means of the clusters.
	 */
	public double[] mu;
	/**
	 * The mixing proportions.
	 */
	public double[] pi;
	/**
	 * The cluster membership probabilties.
	 */
	public double[][] z;
	/**
	 * The observed loglikelihood.
	 */
	public double olog;
	/**
	 * The complete loglikelihood.
	 */
	public double clog;
	/**
	 * The matrix of convergence rate
	 * for the EM-Algorithm (just for means).
	 */
	public double[][] dm;
	/**
	 * The observed information (just for means). 
	 */
	public double[][] ocov;
	/**
	 * The complete information (just for means).
	 */
	public double[] ccov;
	/**
	 * The model choice criteria.
	 */
	public EMCriteria criteria;
	
	/**
	 * 
	 */
	public EMResult() {}
	/**
	 * @param n the number of data points.
	 * @param k the number of components.
	 */
	public EMResult(int n, int k) {
		this.n = n;
		this.k = k;
		mu = new double[k];
		pi = new double[k];
		z  = new double[k][n];
		ocov = new double[k][k];
		dm   = new double[k][k];
		ccov = new double[k];
		olog = Double.NEGATIVE_INFINITY;
		clog = Double.NEGATIVE_INFINITY;
		criteria = new EMCriteria();
		rate  = 0.0;
		edist = 1.0;
	}

	

	/**
	 * Sort the components in increasing order.
	 */
	public void sortCluster() {
		if (k > 1) {
			int i;
			/* First test if it is necessary to
			 * sort the means 
			 */
			for(i=0;i<k-1;i++)
				if (mu[i] > mu[i+1])
					break;
				
			if (i<k-1) {
				
				int[] idx;
				double[] muu, pii;
				double[][] zz;
				
				idx = new int[k+1];
				muu = new double[k+1];
				pii = new double[k+1];
				zz  = new double[k+1][];
				
				for(i=0;i<k;i++) {
					muu[i+1]=mu[i];
					pii[i+1]=pi[i];
					zz[i+1] =z[i];
				}
				
				NumericalUtilities.indexx(k, muu, idx);
				
				for(i=0;i<k;i++) {
					mu[i] = muu[idx[i+1]];
					pi[i] = pii[idx[i+1]];
					z[i]  = zz[idx[i+1]];
				}
				
			}
		}
	}

	/**
	 * @param theta
	 * @param spoint
	 */
	public static void copy(EMResult dest, EMResult src) {
		if (src == null) return;
		if (dest == null) return;
		if (dest.k != src.k) return;
		if (dest.n != src.n) return;
		
		int i,j;
		
		for(i=0;i<src.k;i++) {
			dest.mu[i] = src.mu[i];
			dest.pi[i] = src.pi[i];
			for(j=0;j<src.n;j++)
				dest.z[i][j]=src.z[i][j];
		}
		
		dest.olog = src.olog;
		dest.clog = src.clog;
		
	}

	/**
	 * This methos computes the following criteria :
	 * <ul>
	 * <li> AIC </li>
	 * <li> AIC3 </li>
	 * <li> ICOMP </li>
	 * <li> MIR   <li>
	 * <li> BIC </li>
	 * <li> AWE </li>
	 * </ul>
	 */
	public void computeCriteria() {
		int i;
		double t,d;
		double nfp = (2.0*k-1); /* the number of free parameters */
		
		if (k==n) nfp=n;
		
		criteria.aic   = -2.0*olog + 2.0*nfp;
		if (n-nfp > 1)
			criteria.aicc  = -2.0*olog + 2.0*nfp + 2.0*(double)(nfp*(nfp+1))/(double)(n-nfp-1);
		else criteria.aicc  = -2.0*olog + 2.0*nfp;
		criteria.aic3  = -2.0*olog + 3.0*nfp;
		criteria.bic   = -2.0*olog + Math.log(n)*nfp;
		criteria.awe   = -2.0*clog + 2.0*(3.0/2.0 + Math.log(n))*nfp;
		t=0.0;d=1.0;
		for(i=0;i<k;i++) {
			t+=ocov[i][i];
			d*=ocov[i][i];
		}
		criteria.icomp = -2.0*olog + 0.5*nfp*Math.log(t/nfp) - 0.5*Math.log(d);
		criteria.mir   = 1.0-rate;
	}
	/**
	 * @param critIdx
	 */
	public double getCriterion(int critIdx) {
		return criteria.getCriterion(critIdx);
	}
	/**
	 * @param sd
	 */
	public double[] getMahalanobis(double[] sd) {
		if (sd == null) return null;
		if (k == 1) return null; 
		double[] d = new double[k-1];
		double[] a = new double[k];
		double[] s = new double[k];
		for(int i=0;i<n;i++) {
			for(int j=0;j<k;j++) {
				a[j]+=z[j][i]*sd[i];
				s[j]+=z[j][i];
			}			
		}
		for(int i=0;i<k;i++) a[i] /= s[i];
		for(int i=0;i<k-1;i++) {
			d[i]=(mu[i]-mu[i+1])*(mu[i]-mu[i+1]);
			d[i]/=(a[i]*a[i]+a[i+1]*a[i+1]);
			d[i]=Math.sqrt(d[i]);
		}
		return d;
	}
	/**
	 * @return
	 */
	public double[] getXPred() {
		double[] xp = new double[n];
		for(int i=0;i<n;i++) {
			for(int j=0;j<k;j++) {
				xp[i]+=z[j][i]*mu[j];
			}
		}
		return xp;
	}
	/**
	 * @param i
	 * @return
	 */
	public double getSD(int i) {
		return Math.sqrt(ocov[i][i]);
	}
	/**
	 * @param i
	 * @return
	 */
	public double getMu(int i) {
		return mu[i];
	}
	/**
	 * @return
	 */
	public int getK() {
		return k;
	}
	
	public double[] getZ(int i) {
		double[] zi = new double[k];
		for(int j=0;j<k;j++)
			zi[j]=z[j][i];
		return zi;
	}
	/**
	 * @param sd
	 * @return
	 */
	public double[] getEuclidean() {
		double[] d= new double[k-1];
		for(int i=0;i<k-1;i++) {
			d[i]=(mu[i]-mu[i+1])*(mu[i]-mu[i+1]);
			d[i]=Math.sqrt(d[i]);
		}
		return d;
	}
	/**
	 * @param i
	 * @return
	 */
	public double getPi(int i) {
		return pi[i];
	}
	/**
	 * @param i
	 * @return
	 */
	public int getXPredIdx(int i) {
		double max = 0.0;
		int maxk   = 0;
		for(int j=0;j<k;j++) {
			if (z[j][i] > max) {
				max = z[j][i];
				maxk = j;
			}
		}
		return maxk;
	}
	/**
	 * @return
	 */
	public double[] getXBestPred() {
		double[] xp = new double[n];
		for(int i=0;i<n;i++) {
			xp[i] = getMu(getXPredIdx(i)); 
		}
		return xp;
	}
}
