+++ /dev/null
-#include "muscle.h"\r
-#include "tree.h"\r
-#include <stdio.h>\r
-\r
-#define TRACE 0\r
-\r
-void ClusterByHeight(const Tree &tree, double dMaxHeight, unsigned Subtrees[],\r
- unsigned *ptruSubtreeCount)\r
- {\r
- if (!tree.IsRooted())\r
- Quit("ClusterByHeight: requires rooted tree");\r
-\r
-#if TRACE\r
- Log("ClusterByHeight, max height=%g\n", dMaxHeight);\r
-#endif\r
-\r
- unsigned uSubtreeCount = 0;\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
- {\r
- if (tree.IsRoot(uNodeIndex))\r
- continue;\r
- unsigned uParent = tree.GetParent(uNodeIndex);\r
- double dHeight = tree.GetNodeHeight(uNodeIndex);\r
- double dParentHeight = tree.GetNodeHeight(uParent);\r
-\r
-#if TRACE\r
- Log("Node %3u Height %5.2f ParentHeight %5.2f\n",\r
- uNodeIndex, dHeight, dParentHeight);\r
-#endif\r
- if (dParentHeight > dMaxHeight && dHeight <= dMaxHeight)\r
- {\r
- Subtrees[uSubtreeCount] = uNodeIndex;\r
-#if TRACE\r
- Log("Subtree[%u]=%u\n", uSubtreeCount, uNodeIndex);\r
-#endif\r
- ++uSubtreeCount;\r
- }\r
- }\r
- *ptruSubtreeCount = uSubtreeCount;\r
- }\r
-\r
-static void ClusterBySubfamCount_Iteration(const Tree &tree, unsigned Subfams[],\r
- unsigned uCount)\r
- {\r
-// Find highest child node of current set of subfamilies.\r
- double dHighestHeight = -1e20;\r
- int iParentSubscript = -1;\r
-\r
- for (int n = 0; n < (int) uCount; ++n)\r
- {\r
- const unsigned uNodeIndex = Subfams[n];\r
- if (tree.IsLeaf(uNodeIndex))\r
- continue;\r
-\r
- const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
- const double dHeightLeft = tree.GetNodeHeight(uLeft);\r
- if (dHeightLeft > dHighestHeight)\r
- {\r
- dHighestHeight = dHeightLeft;\r
- iParentSubscript = n;\r
- }\r
-\r
- const unsigned uRight = tree.GetRight(uNodeIndex);\r
- const double dHeightRight = tree.GetNodeHeight(uRight);\r
- if (dHeightRight > dHighestHeight)\r
- {\r
- dHighestHeight = dHeightRight;\r
- iParentSubscript = n;\r
- }\r
- }\r
-\r
- if (-1 == iParentSubscript)\r
- Quit("CBSFCIter: failed to find highest child");\r
-\r
- const unsigned uNodeIndex = Subfams[iParentSubscript];\r
- const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
- const unsigned uRight = tree.GetRight(uNodeIndex);\r
-\r
-// Delete parent by replacing with left child\r
- Subfams[iParentSubscript] = uLeft;\r
-\r
-// Append right child to list\r
- Subfams[uCount] = uRight;\r
-\r
-#if TRACE\r
- {\r
- Log("Iter %3u:", uCount);\r
- for (unsigned n = 0; n < uCount; ++n)\r
- Log(" %u", Subfams[n]);\r
- Log("\n");\r
- }\r
-#endif\r
- }\r
-\r
-// Divide a tree containing N leaves into k families by\r
-// cutting the tree at a horizontal line at some height.\r
-// Each internal node defines a height for the cut, \r
-// considering all internal nodes enumerates all distinct\r
-// cuts. Visit internal nodes in decreasing order of height.\r
-// Visiting the node corresponds to moving the horizontal\r
-// line down to cut the tree at the height of that node.\r
-// We consider the cut to be "infinitestimally below"\r
-// the node, so the effect is to remove the current node \r
-// from the list of subfamilies and add its two children.\r
-// We must visit a parent before its children (so care may\r
-// be needed to handle zero edge lengths properly).\r
-// We assume that N is small, and write dumb O(N^2) code.\r
-// More efficient strategies are possible for large N\r
-// by maintaining a list of nodes sorted by height.\r
-void ClusterBySubfamCount(const Tree &tree, unsigned uSubfamCount,\r
- unsigned Subfams[], unsigned *ptruSubfamCount)\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- const unsigned uLeafCount = (uNodeCount + 1)/2;\r
-\r
-// Special case: empty tree\r
- if (0 == uNodeCount)\r
- {\r
- *ptruSubfamCount = 0;\r
- return;\r
- }\r
-\r
-// Special case: more subfamilies than leaves\r
- if (uSubfamCount >= uLeafCount)\r
- {\r
- for (unsigned n = 0; n < uLeafCount; ++n)\r
- Subfams[n] = n;\r
- *ptruSubfamCount = uLeafCount;\r
- return;\r
- }\r
-\r
-// Initialize list of subfamilies to be root\r
- Subfams[0] = tree.GetRootNodeIndex();\r
-\r
-// Iterate\r
- for (unsigned i = 1; i < uSubfamCount; ++i)\r
- ClusterBySubfamCount_Iteration(tree, Subfams, i);\r
- \r
- *ptruSubfamCount = uSubfamCount;\r
- }\r
-\r
-static void GetLeavesRecurse(const Tree &tree, unsigned uNodeIndex,\r
- unsigned Leaves[], unsigned &uLeafCount /* in-out */)\r
- {\r
- if (tree.IsLeaf(uNodeIndex))\r
- {\r
- Leaves[uLeafCount] = uNodeIndex;\r
- ++uLeafCount;\r
- return;\r
- }\r
-\r
- const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
- const unsigned uRight = tree.GetRight(uNodeIndex);\r
-\r
- GetLeavesRecurse(tree, uLeft, Leaves, uLeafCount);\r
- GetLeavesRecurse(tree, uRight, Leaves, uLeafCount);\r
- }\r
-\r
-void GetLeaves(const Tree &tree, unsigned uNodeIndex, unsigned Leaves[],\r
- unsigned *ptruLeafCount)\r
- {\r
- unsigned uLeafCount = 0;\r
- GetLeavesRecurse(tree, uNodeIndex, Leaves, uLeafCount);\r
- *ptruLeafCount = uLeafCount;\r
- }\r
-\r
-void Tree::PruneTree(const Tree &tree, unsigned Subfams[],\r
- unsigned uSubfamCount)\r
- {\r
- if (!tree.IsRooted())\r
- Quit("Tree::PruneTree: requires rooted tree");\r
-\r
- Clear();\r
-\r
- m_uNodeCount = 2*uSubfamCount - 1;\r
- InitCache(m_uNodeCount);\r
-\r
- const unsigned uUnprunedNodeCount = tree.GetNodeCount();\r
-\r
- unsigned *uUnprunedToPrunedIndex = new unsigned[uUnprunedNodeCount];\r
- unsigned *uPrunedToUnprunedIndex = new unsigned[m_uNodeCount];\r
-\r
- for (unsigned n = 0; n < uUnprunedNodeCount; ++n)\r
- uUnprunedToPrunedIndex[n] = NULL_NEIGHBOR;\r
-\r
- for (unsigned n = 0; n < m_uNodeCount; ++n)\r
- uPrunedToUnprunedIndex[n] = NULL_NEIGHBOR;\r
-\r
-// Create mapping between unpruned and pruned node indexes\r
- unsigned uInternalNodeIndex = uSubfamCount;\r
- for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
- {\r
- unsigned uUnprunedNodeIndex = Subfams[uSubfamIndex];\r
- uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uSubfamIndex;\r
- uPrunedToUnprunedIndex[uSubfamIndex] = uUnprunedNodeIndex;\r
- for (;;)\r
- {\r
- uUnprunedNodeIndex = tree.GetParent(uUnprunedNodeIndex);\r
- if (tree.IsRoot(uUnprunedNodeIndex))\r
- break;\r
-\r
- // Already visited this node?\r
- if (NULL_NEIGHBOR != uUnprunedToPrunedIndex[uUnprunedNodeIndex])\r
- break;\r
-\r
- uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uInternalNodeIndex;\r
- uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedNodeIndex;\r
-\r
- ++uInternalNodeIndex;\r
- }\r
- }\r
-\r
- const unsigned uUnprunedRootIndex = tree.GetRootNodeIndex();\r
- uUnprunedToPrunedIndex[uUnprunedRootIndex] = uInternalNodeIndex;\r
- uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedRootIndex;\r
-\r
-#if TRACE\r
- {\r
- Log("Pruned to unpruned:\n");\r
- for (unsigned i = 0; i < m_uNodeCount; ++i)\r
- Log(" [%u]=%u", i, uPrunedToUnprunedIndex[i]);\r
- Log("\n");\r
- Log("Unpruned to pruned:\n");\r
- for (unsigned i = 0; i < uUnprunedNodeCount; ++i)\r
- {\r
- unsigned n = uUnprunedToPrunedIndex[i];\r
- if (n != NULL_NEIGHBOR)\r
- Log(" [%u]=%u", i, n);\r
- }\r
- Log("\n");\r
- }\r
-#endif\r
-\r
- if (uInternalNodeIndex != m_uNodeCount - 1)\r
- Quit("Tree::PruneTree, Internal error");\r
-\r
-// Nodes 0, 1 ... are the leaves\r
- for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
- {\r
- char szName[32];\r
- sprintf(szName, "Subfam_%u", uSubfamIndex + 1);\r
- m_ptrName[uSubfamIndex] = strsave(szName);\r
- }\r
-\r
- for (unsigned uPrunedNodeIndex = uSubfamCount; uPrunedNodeIndex < m_uNodeCount;\r
- ++uPrunedNodeIndex)\r
- {\r
- unsigned uUnprunedNodeIndex = uPrunedToUnprunedIndex[uPrunedNodeIndex];\r
-\r
- const unsigned uUnprunedLeft = tree.GetLeft(uUnprunedNodeIndex);\r
- const unsigned uUnprunedRight = tree.GetRight(uUnprunedNodeIndex);\r
-\r
- const unsigned uPrunedLeft = uUnprunedToPrunedIndex[uUnprunedLeft];\r
- const unsigned uPrunedRight = uUnprunedToPrunedIndex[uUnprunedRight];\r
-\r
- const double dLeftLength =\r
- tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedLeft);\r
- const double dRightLength =\r
- tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedRight);\r
-\r
- m_uNeighbor2[uPrunedNodeIndex] = uPrunedLeft;\r
- m_uNeighbor3[uPrunedNodeIndex] = uPrunedRight;\r
-\r
- m_dEdgeLength1[uPrunedLeft] = dLeftLength;\r
- m_dEdgeLength1[uPrunedRight] = dRightLength;\r
-\r
- m_uNeighbor1[uPrunedLeft] = uPrunedNodeIndex;\r
- m_uNeighbor1[uPrunedRight] = uPrunedNodeIndex;\r
-\r
- m_bHasEdgeLength1[uPrunedLeft] = true;\r
- m_bHasEdgeLength1[uPrunedRight] = true;\r
-\r
- m_dEdgeLength2[uPrunedNodeIndex] = dLeftLength;\r
- m_dEdgeLength3[uPrunedNodeIndex] = dRightLength;\r
-\r
- m_bHasEdgeLength2[uPrunedNodeIndex] = true;\r
- m_bHasEdgeLength3[uPrunedNodeIndex] = true;\r
- }\r
-\r
- m_uRootNodeIndex = uUnprunedToPrunedIndex[uUnprunedRootIndex];\r
-\r
- m_bRooted = true;\r
-\r
- Validate();\r
-\r
- delete[] uUnprunedToPrunedIndex;\r
- }\r
-\r
-void LeafIndexesToIds(const Tree &tree, const unsigned Leaves[], unsigned uCount,\r
- unsigned Ids[])\r
- {\r
- for (unsigned n = 0; n < uCount; ++n)\r
- Ids[n] = tree.GetLeafId(Leaves[n]);\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff