/*  
 *  src/org/metaqtl/algo/QtlClustAlgorithm.java
 * 
 *  $Author: veyrieras $, $Date: 2006-10-17 15:41:55 $, $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.algo;

import org.metaqtl.EMResult;
import org.metaqtl.MetaQtlData;
import org.metaqtl.MetaQtlResult;
import org.metaqtl.adapter.MetaQtlAnalysisAdapter;
import org.metaqtl.bio.IBioGenome;
import org.metaqtl.bio.IBioOntology;
/**
 * This class defines methods to perform a QTL meta-analysis
 * on a genetic map by applying a clustering of QTL
 * based on a gaussian mixture model and the EM-algorihtm. 
 */
public class QtlClustAlgorithm extends ClustAlgorithm {
	
	/**
	 * The maximum number of clusters.
	 */
	private int nmax_cluster	  = 0;
	
	private int n_em_start        = EMAlgorithm.EM_START;
	
	private double em_eps         = EMAlgorithm.EM_ERR;
	/**
	 * @param map
	 * @param ontology
	 */
	public QtlClustAlgorithm(IBioGenome map, IBioOntology ontology) {
		super(map, ontology);
	}
	
	public QtlClustAlgorithm(IBioGenome map, String chrom, IBioOntology ontology) {
		super(map, chrom, ontology);
	}
	/**
	 * This method performs a QTL meta-analysis on each chromosome 
	 * of the given map. For each chromosome all the possible
	 * QTL cluster combinations are tested and the results are
	 * added to the chromosome. Then the method returns the array
	 * of chromosomes for which the meta-analysis has been performed.
	 * 
	 */
	public void run() {
		int i;
		MetaQtlData data;
		
		/* Loop on linkage groups */
		for(i=0;i<chromosomes.length;i++) {
			
			if (chromosomes[i].hasQTL()) {
				
				boolean doClust = (chrName==null || (chrName!=null && chrName.equals(chromosomes[i].getName())));
				
				if (doClust) {
				
					workProgress++;
					/* First compute the standard deviation of the 
					 * QTL positions */  
					chromosomes[i].computeQtlSD(sd_mode);
					/* Then adapt the QTL in the chromosome into a metaQTLData object */
					data      = new MetaQtlData(chromosomes[i].getQtls());
					
					if (ontology == null)
						data.doTraitGroupClustering();
					else
						data.doTraitOntologyClustering(ontology);
					
					/* Do imputation for qtl with missing CI */
					data.manageMissingData(missing_sd_mode);
					/* Do the clustering */
					
					if (isLoggerEnable()) {
						getLogger().println("Start Clustering on chromosome " + chromosomes[i].getName());
						getLogger().flush();
					}
					
					chromosomes[i].metaQtlAnalysis = doQtlClustering(data);
					
				}
				
			}
		}
		
		result = MetaQtlAnalysisAdapter.adapt(chromosomes);
	}
	
	/**
	 * For a given data set <code>data</code> performs the clustering
	 * on each trait group and returns the results as an array of 
	 * {@link MetaQtlResult} objects.
	 * 
	 * @param data the data set.
	 * @return the results of teh clustering.
	 * 
	 * @see MetaQtlData
	 * @see MetaQtlResult
	 */
	public MetaQtlResult[] doQtlClustering(MetaQtlData data) {
		int i,j,ntg;
		Double[][] X;
		MetaQtlResult[] results = null;
		
		if (data==null) return null;
		
		ntg = data.getTraitClusterNumber();
		
		results = new MetaQtlResult[ntg];
		
		for(i=0;i<ntg;i++) {
			
			X                  = data.getDataPoints(i, true);
			
			if (X != null) {
				
				String trait = data.getTraitClusterName(i);
				
				if (isLoggerEnable()) {
					getLogger().println("Clustering for trait " + trait);
					getLogger().flush();
				}
				
				results[i]         = doDataPointClustering(trait, X[0], X[1]);
				
				if (results[i] != null) {
					
					results[i].setQtlIdx(data.getQtlIdx(i));
					
						
				} else {
					
					if (isLoggerEnable()) {
						getLogger().print("[ WARNING ] : ");
						getLogger().println("Any result for trait cluster " + trait);
						getLogger().flush();
					}
				}
				
			}
		}
		
		return results;
	}
	/**
	 * For the given data point where <code>xx</code> is the vector of
	 * observed values and <code>ssd</code> their standard deviations this
	 * method performs the EM algorithm for each possible value of the number
	 * of cluster, i.e from one to the number of data points and returns the
	 * results as a {@link MetaQtlResult} object.
	 * @param xx the vector of observed values.
	 * @param sdd the vector of standard deviations.
	 * @return the result of the clustering.
	 * 
	 * @see MetaQtlResult
	 * @see EMAlgorithm
	 */
	private MetaQtlResult doDataPointClustering(String trait, Double[] x, Double[] sd) {
		int i,n;
		MetaQtlResult result = null;
		EMResult[] emResults = null;
		
		
		if (x == null || sd == null) return null;
		if (x.length < 2) return null;
		
		n = x.length; /* the number of positions */
		
		result = new MetaQtlResult(trait,n);
		
		result.setX(x);
		result.setSD(sd);
		
		if (nmax_cluster > 0 && nmax_cluster < n) { 
			n=nmax_cluster;
			emResults = new EMResult[n + 1];
		} else emResults = new EMResult[n];
		
		
		EMAlgorithm.EM_START = n_em_start;
		EMAlgorithm.EM_ERR   = em_eps;
		
		for(i=0;i<n;i++) {
			
			emResults[i] = EMAlgorithm.doEM(result.x, result.sd, i+1, null);
			
			if (emResults[i] != null)
				emResults[i].computeCriteria();
			
			if (isLoggerEnable()) {
				getLogger().println("EM [ rs=" + n_em_start +", eps="+em_eps+" ]");
				getLogger().println("Cluster K=" + (i+1) +" : [done]");
				getLogger().flush();
			}
			
		}
		
		if (n < x.length) {
			
			emResults[i] = EMAlgorithm.doEM(result.x, result.sd, x.length, null);
			emResults[i].computeCriteria();
			
			if (isLoggerEnable()) {
				getLogger().println("Cluster K=" + x.length +" : [done]");
				getLogger().flush();
			}
		}
		
		result.setClusterings(emResults);
		
		return result;
	}
	/**
	 * 
	 * 
	 * @param nmax the maximum number of clusters to test.
	 */
	public void setMaxClusterNumber(int nmax) {
		this.nmax_cluster = nmax;
	}
	
	/**
	 * @return Returns the number of random starting points.
	 */
	public int getEMRanStartNumber() {
		return n_em_start;
	}
	/**
	 * @param n_em_start The n_em_start to set.
	 */
	public void setEMRanStartNumber(int n) {
		this.n_em_start = n;
	}
	/**
	 * @return Returns the convergence error for the EM.
	 */
	public double getEMEps() {
		return em_eps;
	}
	/**
	 * @param em_err The em_err to set.
	 */
	public void setEMEps(double epsilon) {
		this.em_eps = epsilon;
	}
}
