source: trunk/GDE/MUSCLE/src/phy4.cpp

#include "muscle.h"
#include "tree.h"
#include <stdio.h>

#define TRACE   0

void ClusterByHeight(const Tree &tree, double dMaxHeight, unsigned Subtrees[],
  unsigned *ptruSubtreeCount)
        {
        if (!tree.IsRooted())
                Quit("ClusterByHeight: requires rooted tree");

#if     TRACE
        Log("ClusterByHeight, max height=%g\n", dMaxHeight);
#endif

        unsigned uSubtreeCount = 0;
        const unsigned uNodeCount = tree.GetNodeCount();
        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
                {
                if (tree.IsRoot(uNodeIndex))
                        continue;
                unsigned uParent = tree.GetParent(uNodeIndex);
                double dHeight = tree.GetNodeHeight(uNodeIndex);
                double dParentHeight = tree.GetNodeHeight(uParent);

#if     TRACE
                Log("Node %3u  Height %5.2f  ParentHeight %5.2f\n",
                  uNodeIndex, dHeight, dParentHeight);
#endif
                if (dParentHeight > dMaxHeight && dHeight <= dMaxHeight)
                        {
                        Subtrees[uSubtreeCount] = uNodeIndex;
#if     TRACE
                        Log("Subtree[%u]=%u\n", uSubtreeCount, uNodeIndex);
#endif
                        ++uSubtreeCount;
                        }
                }
        *ptruSubtreeCount = uSubtreeCount;
        }

static void ClusterBySubfamCount_Iteration(const Tree &tree, unsigned Subfams[],
  unsigned uCount)
        {
// Find highest child node of current set of subfamilies.
        double dHighestHeight = -1e20;
        int iParentSubscript = -1;

        for (int n = 0; n < (int) uCount; ++n)
                {
                const unsigned uNodeIndex = Subfams[n];
                if (tree.IsLeaf(uNodeIndex))
                        continue;

                const unsigned uLeft = tree.GetLeft(uNodeIndex);
                const double dHeightLeft = tree.GetNodeHeight(uLeft);
                if (dHeightLeft > dHighestHeight)
                        {
                        dHighestHeight = dHeightLeft;
                        iParentSubscript = n;
                        }

                const unsigned uRight = tree.GetRight(uNodeIndex);
                const double dHeightRight = tree.GetNodeHeight(uRight);
                if (dHeightRight > dHighestHeight)
                        {
                        dHighestHeight = dHeightRight;
                        iParentSubscript = n;
                        }
                }

        if (-1 == iParentSubscript)
                Quit("CBSFCIter: failed to find highest child");

        const unsigned uNodeIndex = Subfams[iParentSubscript];
        const unsigned uLeft = tree.GetLeft(uNodeIndex);
        const unsigned uRight = tree.GetRight(uNodeIndex);

// Delete parent by replacing with left child
        Subfams[iParentSubscript] = uLeft;

// Append right child to list
        Subfams[uCount] = uRight;

#if     TRACE
        {
        Log("Iter %3u:", uCount);
        for (unsigned n = 0; n < uCount; ++n)
                Log(" %u", Subfams[n]);
        Log("\n");
        }
#endif
        }

// Divide a tree containing N leaves into k families by
// cutting the tree at a horizontal line at some height.
// Each internal node defines a height for the cut, 
// considering all internal nodes enumerates all distinct
// cuts. Visit internal nodes in decreasing order of height.
// Visiting the node corresponds to moving the horizontal
// line down to cut the tree at the height of that node.
// We consider the cut to be "infinitestimally below"
// the node, so the effect is to remove the current node 
// from the list of subfamilies and add its two children.
// We must visit a parent before its children (so care may
// be needed to handle zero edge lengths properly).
// We assume that N is small, and write dumb O(N^2) code.
// More efficient strategies are possible for large N
// by maintaining a list of nodes sorted by height.
void ClusterBySubfamCount(const Tree &tree, unsigned uSubfamCount,
  unsigned Subfams[], unsigned *ptruSubfamCount)
        {
        const unsigned uNodeCount = tree.GetNodeCount();
        const unsigned uLeafCount = (uNodeCount + 1)/2;

// Special case: empty tree
        if (0 == uNodeCount)
                {
                *ptruSubfamCount = 0;
                return;
                }

// Special case: more subfamilies than leaves
        if (uSubfamCount >= uLeafCount)
                {
                for (unsigned n = 0; n < uLeafCount; ++n)
                        Subfams[n] = n;
                *ptruSubfamCount = uLeafCount;
                return;
                }

// Initialize list of subfamilies to be root
        Subfams[0] = tree.GetRootNodeIndex();

// Iterate
        for (unsigned i = 1; i < uSubfamCount; ++i)
                ClusterBySubfamCount_Iteration(tree, Subfams, i);
        
        *ptruSubfamCount = uSubfamCount;
        }

static void GetLeavesRecurse(const Tree &tree, unsigned uNodeIndex,
  unsigned Leaves[], unsigned &uLeafCount /* in-out */)
        {
        if (tree.IsLeaf(uNodeIndex))
                {
                Leaves[uLeafCount] = uNodeIndex;
                ++uLeafCount;
                return;
                }

        const unsigned uLeft = tree.GetLeft(uNodeIndex);
        const unsigned uRight = tree.GetRight(uNodeIndex);

        GetLeavesRecurse(tree, uLeft, Leaves, uLeafCount);
        GetLeavesRecurse(tree, uRight, Leaves, uLeafCount);
        }

void GetLeaves(const Tree &tree, unsigned uNodeIndex, unsigned Leaves[],
  unsigned *ptruLeafCount)
        {
        unsigned uLeafCount = 0;
        GetLeavesRecurse(tree, uNodeIndex, Leaves, uLeafCount);
        *ptruLeafCount = uLeafCount;
        }

void Tree::PruneTree(const Tree &tree, unsigned Subfams[],
  unsigned uSubfamCount)
        {
        if (!tree.IsRooted())
                Quit("Tree::PruneTree: requires rooted tree");

        Clear();

        m_uNodeCount = 2*uSubfamCount - 1;
        InitCache(m_uNodeCount);

        const unsigned uUnprunedNodeCount = tree.GetNodeCount();

        unsigned *uUnprunedToPrunedIndex = new unsigned[uUnprunedNodeCount];
        unsigned *uPrunedToUnprunedIndex = new unsigned[m_uNodeCount];

        for (unsigned n = 0; n < uUnprunedNodeCount; ++n)
                uUnprunedToPrunedIndex[n] = NULL_NEIGHBOR;

        for (unsigned n = 0; n < m_uNodeCount; ++n)
                uPrunedToUnprunedIndex[n] = NULL_NEIGHBOR;

// Create mapping between unpruned and pruned node indexes
        unsigned uInternalNodeIndex = uSubfamCount;
        for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
                {
                unsigned uUnprunedNodeIndex = Subfams[uSubfamIndex];
                uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uSubfamIndex;
                uPrunedToUnprunedIndex[uSubfamIndex] = uUnprunedNodeIndex;
                for (;;)
                        {
                        uUnprunedNodeIndex = tree.GetParent(uUnprunedNodeIndex);
                        if (tree.IsRoot(uUnprunedNodeIndex))
                                break;

                // Already visited this node?
                        if (NULL_NEIGHBOR != uUnprunedToPrunedIndex[uUnprunedNodeIndex])
                                break;

                        uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uInternalNodeIndex;
                        uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedNodeIndex;

                        ++uInternalNodeIndex;
                        }
                }

        const unsigned uUnprunedRootIndex = tree.GetRootNodeIndex();
        uUnprunedToPrunedIndex[uUnprunedRootIndex] = uInternalNodeIndex;
        uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedRootIndex;

#if     TRACE
        {
        Log("Pruned to unpruned:\n");
        for (unsigned i = 0; i < m_uNodeCount; ++i)
                Log(" [%u]=%u", i, uPrunedToUnprunedIndex[i]);
        Log("\n");
        Log("Unpruned to pruned:\n");
        for (unsigned i = 0; i < uUnprunedNodeCount; ++i)
                {
                unsigned n = uUnprunedToPrunedIndex[i];
                if (n != NULL_NEIGHBOR)
                        Log(" [%u]=%u", i, n);
                }
        Log("\n");
        }
#endif

        if (uInternalNodeIndex != m_uNodeCount - 1)
                Quit("Tree::PruneTree, Internal error");

// Nodes 0, 1 ... are the leaves
        for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
                {
                char szName[32];
                sprintf(szName, "Subfam_%u", uSubfamIndex + 1);
                m_ptrName[uSubfamIndex] = strsave(szName);
                }

        for (unsigned uPrunedNodeIndex = uSubfamCount; uPrunedNodeIndex < m_uNodeCount;
          ++uPrunedNodeIndex)
                {
                unsigned uUnprunedNodeIndex = uPrunedToUnprunedIndex[uPrunedNodeIndex];

                const unsigned uUnprunedLeft = tree.GetLeft(uUnprunedNodeIndex);
                const unsigned uUnprunedRight = tree.GetRight(uUnprunedNodeIndex);

                const unsigned uPrunedLeft = uUnprunedToPrunedIndex[uUnprunedLeft];
                const unsigned uPrunedRight = uUnprunedToPrunedIndex[uUnprunedRight];

                const double dLeftLength =
                  tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedLeft);
                const double dRightLength =
                  tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedRight);

                m_uNeighbor2[uPrunedNodeIndex] = uPrunedLeft;
                m_uNeighbor3[uPrunedNodeIndex] = uPrunedRight;

                m_dEdgeLength1[uPrunedLeft] = dLeftLength;
                m_dEdgeLength1[uPrunedRight] = dRightLength;

                m_uNeighbor1[uPrunedLeft] = uPrunedNodeIndex;
                m_uNeighbor1[uPrunedRight] = uPrunedNodeIndex;

                m_bHasEdgeLength1[uPrunedLeft] = true;
                m_bHasEdgeLength1[uPrunedRight] = true;

                m_dEdgeLength2[uPrunedNodeIndex] = dLeftLength;
                m_dEdgeLength3[uPrunedNodeIndex] = dRightLength;

                m_bHasEdgeLength2[uPrunedNodeIndex] = true;
                m_bHasEdgeLength3[uPrunedNodeIndex] = true;
                }

        m_uRootNodeIndex = uUnprunedToPrunedIndex[uUnprunedRootIndex];

        m_bRooted = true;

        Validate();

        delete[] uUnprunedToPrunedIndex;
        }

void LeafIndexesToIds(const Tree &tree, const unsigned Leaves[], unsigned uCount,
  unsigned Ids[])
        {
        for (unsigned n = 0; n < uCount; ++n)
                Ids[n] = tree.GetLeafId(Leaves[n]);
        }