/*  
 *  src/org/metaqtl/algo/HClustAlgorithm.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.Tree;
import org.metaqtl.TreeNode;

/**
 * This class implements stantard agglomerative hierarchical
 * clustering methods.
 */
public class HClustAlgorithm {	  	
	 /**
	  * Centroid method
	  */
	 public static final int CENTROID_METHOD = 1;
	 /**
	  * Ward method
	  */
	 public static final int WARD_METHOD     = 2;
	 /**
	  * A global variable to control the custering
	  * method.
	  */
	 public static int METHOD = CENTROID_METHOD;
	 /**
	  * This method performs a agglomerative clustering
	  * of the given data points where <code>x</code> is the 
	  * vector of observed values and <code>sd</code> the
	  * vector of their standard deviations. It returns an
	  * array of {@link TreeNode} which represents the 
	  * hierarchical clustering.  
	  * 
	  * @param x  the vector of data points.
	  * @param sd the vector of standard deviations.
	  * @return the tree which represents the hierarchical clustering.
	  */
	  public static Tree run(double[] x, double[] sd) {
	  	switch(METHOD) {
	  		case CENTROID_METHOD :
	  			return CentroidMethod(x,sd);
	  		case WARD_METHOD :
	  			return WardMethod(x,sd);
	  		default :
	  			return CentroidMethod(x,sd);
	  	}
	  }
	  /**
	   * This method performs a hierarchical clustering
	   * of the given data points using a centroid method
	   * based on the mahalanobis distance.
	   * @param x  the vector of data points.
	  *  @param sd the vector of standard deviations.
	  *  @return the tree which represents the hierarchical clustering.
	   */
	  public static Tree CentroidMethod(double[] x, double[] sd) {
	  	int i,j,ii,n,k,mini,minj,iter;
	  	double t,s,wi,wj,dmin;
	  	double[] xx,var2,var;
	  	TreeNode[] nodes=null;
	  	 
	  	k = n = x.length;
	  	
	  	nodes = new TreeNode[2*n-1];
	  	
	  	xx    = new double[n];
	  	var2   = new double[n];
	  	var   = new double[n];
	  	
	  	for(i=0;i<n;i++) {
	  		nodes[i] = new TreeNode(i);
	  		xx[i]    = x[i];
	  		var2[i]   = var[i] = sd[i]*sd[i];
	  		//System.out.println(i+" "+x[i]+" "+sd[i]);
	  	}
		
	  	for(iter=0;iter<n-1;iter++) {
		  	/* 
		  	* First compute the pairwise
		  	* mahalanobis distance between 
		  	* the k elements and find the min
		  	*/
		  	mini=minj=-1;
		  	dmin=Double.POSITIVE_INFINITY;
			for(i=0;i<k;i++) {
				for(j=i+1;j<k;j++) {
					t=(xx[i]-xx[j])*(xx[i]-xx[j]);
					t/=(var2[i]+var2[j]);
					if (Math.sqrt(t) < dmin) {
						dmin=Math.sqrt(t);
						mini=i;
						minj=j;
					}					
				}
			}
			/*
			 * Then compute the weighted mean
			 * for this cluster. 
			 */
			wj=wi=s=t=0.0;
			for(i=0;i<nodes[mini].card;i++) {
				ii=nodes[mini].nidx[i];
				t +=x[ii]/var[ii];
				s +=1.0/var[ii];
				wi+=var[ii];
			}
			for(i=0;i<nodes[minj].card;i++) {
				ii=nodes[minj].nidx[i];
				t +=x[ii]/var[ii];
				s +=1.0/var[ii];
				wj+=var[ii];
			}
			//System.out.println("Join "+mini+" and "+minj+" "+dmin+" centroid " + t/s + " " + 1.0/s);
			nodes[mini].dist = dmin*wi/(wi+wj);
			nodes[minj].dist = dmin*wj/(wi+wj);
			TreeNode[] children = new TreeNode[2];
			children[0] = nodes[mini];
			children[1] = nodes[minj];
			nodes[n+iter] 				= new TreeNode(n+iter, children);
			nodes[Math.min(mini,minj)]  = nodes[n+iter];
			xx[Math.min(mini,minj)]=t/s;
			var2[Math.min(mini,minj)]=1.0/s;
			for(i=Math.max(mini,minj);i<k-1;i++) {
				xx[i]=xx[i+1];
				var2[i]=var2[i+1];
				nodes[i]=nodes[i+1];
			}
			k--;
	  	}
	  	Tree tree = new Tree();
	  	tree.nodes = nodes;
	  	tree.root  = nodes[nodes.length-1];
		 
	  	return tree; 
	  }
	  /**
	   * This method performs a hierarchical clustering
	   * of the given data points using the Ward's method.
	   * 
	   * @param x  the vector of data points.
	  *  @param sd the vector of standard deviations.
	  *  @return the tree which represents the hierarchical clustering.
	   */
	  public static Tree WardMethod(double[] x, double[] sd) {
	  	int i,j,l,ii,n,k,mini,minj,iter;
	  	double t,s,chi,min;
	  	double[] var;
	  	TreeNode[] nodes=null;
	  	 
	  	k = n = x.length;
	  	
	  	nodes = new TreeNode[2*n-1];
	  	var   = new double[n];
	  	
	  	for(i=0;i<n;i++) {
	  		nodes[i] = new TreeNode(i);
	  		var[i]   = sd[i]*sd[i];
	  		//System.out.println(i+" "+x[i]+" "+sd[i]);
	  	}
		
	  	for(iter=0;iter<n-1;iter++) {
		  	/* 
		  	* For the k(k-1)/2 pairwise combinations
		  	* computes the objective function and find
		  	* which pair minimizes it.
		  	*/
		  	mini=minj=-1;
		  	min=Double.POSITIVE_INFINITY;
			for(i=0;i<k;i++) {
				for(j=i+1;j<k;j++) {
					t=s=0.0;
					for(ii=0;ii<nodes[i].card;ii++) {
						l=nodes[i].nidx[ii];
						t +=x[l]/var[l];
						s +=1.0/var[l];
					}
					for(ii=0;ii<nodes[j].card;ii++) {
						l=nodes[j].nidx[ii];
						t +=x[l]/var[l];
						s +=1.0/var[l];
					}
					t /= s;
					/* Compute the chi-square */
					chi=0.0;
					for(ii=0;ii<nodes[i].card;ii++) {
						l=nodes[i].nidx[ii];
						chi+=(x[l]-t)*(x[l]-t)/var[l];
					}
					for(ii=0;ii<nodes[j].card;ii++) {
						l=nodes[j].nidx[ii];
						chi+=(x[l]-t)*(x[l]-t)/var[l];
					}
					chi /= s;
					//System.out.println("Cluster " +i+" "+j+" "+chi);
					if (chi < min) {
						min=chi;
						mini=i;
						minj=j;
					}					
				}
			}
			//System.out.println("Join "+mini+" and "+minj+" "+min);
			TreeNode[] children = new TreeNode[2];
			nodes[mini].dist = min*0.5;
			nodes[minj].dist = min*0.5;
			children[0] = nodes[mini];
			children[1] = nodes[minj];
			nodes[n+iter] 				= new TreeNode(n+iter, children);
			nodes[Math.min(mini,minj)]  = nodes[n+iter];
			for(i=Math.max(mini,minj);i<k-1;i++)
				nodes[i]=nodes[i+1];
			k--;
	  	}
		 
	  	Tree tree = new Tree();
	  	tree.nodes = nodes;
	  	tree.root  = nodes[nodes.length-1]; 
		 
	  	return tree; 
	  }
}
