// NeighborJoiningTree.java
//
// (c) 1999-2001 PAL Development Core Team
//
// This package may be distributed under the
// terms of the Lesser GNU General Public License (LGPL)


// computational complexity O(numSeqs^3)


package de.ugoe.cs.autoquest.tasktrees.alignment.pal.tree;

import de.ugoe.cs.autoquest.tasktrees.alignment.matrix.UPGMAMatrix;


/**
 * constructs a neighbor-joining tree from pairwise distances
 * 
 * @version $Id: NeighborJoiningTree.java,v 1.9 2001/07/13 14:39:13 korbinian Exp $
 *
 * @author Korbinian Strimmer
 * @author Alexei Drummond
 */
public class NeighborJoiningTree extends SimpleTree
{
	//
	// Public stuff
	//	

	/**
	 * construct NJ tree
	 *
	 * @param m distance matrix
	 */
	public NeighborJoiningTree(UPGMAMatrix m)
	{
		if (m.size() < 3)
		{
			new IllegalArgumentException("LESS THAN 3 TAXA IN DISTANCE MATRIX");
		}
				
		init(m);

		//while (numClusters > 3)
		while (true)
		{
			findNextPair();
			newBranchLengths();
			if (numClusters == 3)
			{
				break;
			}
			newCluster();
		}
		
		finish();
	}


	//
	// Private stuff
	//
	
	private int numClusters;
	private Node newCluster;
	private int besti, abi;
	private int bestj, abj;
	private int[] alias;
	private double[][] distance;
	private double[] r;
	private double scale;	

	private double getDist(int a, int b)
	{
		return distance[alias[a]][alias[b]];
	}
	
	private void init(UPGMAMatrix m)
	{
		numClusters = m.size();

		distance = new double[numClusters][numClusters];
		for (int i = 0; i < numClusters; i++)
		{
			for (int j = 0; j < numClusters; j++)
			{
				distance[i][j] = m.get(i,j);
			}
		}

		for (int i = 0; i < numClusters; i++)
		{
			Node tmp = NodeFactory.createNode();
			//tmp.setIdentifier(m.getIdentifier(i));
			getRoot().addChild(tmp);
		}
		
		alias = new int[numClusters];
		for (int i = 0; i < numClusters; i++)
		{
			alias[i] = i;
		}
		
		r = new double[numClusters];
	}

	private void finish()
	{
		if (besti != 0 && bestj != 0)
		{
			getRoot().getChild(0).setBranchLength(updatedDistance(besti, bestj, 0));
		}
		else if (besti != 1 && bestj != 1)
		{
			getRoot().getChild(1).setBranchLength(updatedDistance(besti, bestj, 1));
		}
		else
		{
			getRoot().getChild(2).setBranchLength(updatedDistance(besti, bestj, 2));
		}
		distance = null;

		// make node heights available also
		NodeUtils.lengths2Heights(getRoot());
	}

	private void findNextPair()
	{
		for (int i = 0; i < numClusters; i++)
		{
			r[i] = 0;
			for (int j = 0; j < numClusters; j++)
			{
				r[i] += getDist(i,j);
			}
		}

		besti = 0;
		bestj = 1;
		double smax = -1.0;
		scale = 1.0/(numClusters-2);
		for (int i = 0; i < numClusters-1; i++)
		{
			for (int j = i+1; j < numClusters; j++)
			{
				double sij = (r[i] + r[j] ) * scale - getDist(i, j);
			
				if (sij > smax)
				{
					smax = sij;
					besti = i;
					bestj = j;
				}
			}
		}
		abi = alias[besti];
		abj = alias[bestj];
	}

	private void newBranchLengths()
	{
		double dij = getDist(besti, bestj);
		double li = (dij + (r[besti]-r[bestj])*scale)*0.5;
		double lj = dij - li; // = (dij + (r[bestj]-r[besti])*scale)*0.5

		getRoot().getChild(besti).setBranchLength(li);
		getRoot().getChild(bestj).setBranchLength(lj);
	}

	private void newCluster()
	{
		// Update distances
		for (int k = 0; k < numClusters; k++)
		{
			if (k != besti && k != bestj)
			{
				int ak = alias[k];	
				distance[ak][abi] = distance[abi][ak] = updatedDistance(besti, bestj, k);
			}
		}
		distance[abi][abi] = 0.0;
		
		// Replace besti with new cluster
		getRoot().joinChildren(besti, bestj);
		
		// Update alias
		for (int i = bestj; i < numClusters-1; i++)
		{
			alias[i] = alias[i+1];
		}
		
		numClusters--;
	}
	
	/**
	 * compute updated distance between the new cluster (i,j)
	 * to any other cluster k
	 */
	private double updatedDistance(int i, int j, int k)
	{
		return (getDist(k, i) + getDist(k, j) - getDist(i, j))*0.5;
	}
}