source: branches/ralph/src/main/java/de/ugoe/cs/autoquest/tasktrees/alignment/pal/tree/NeighborJoiningTree.java @ 1573

// 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;
        }
}