+++ /dev/null
-#include "muscle.h"\r
-#include "tree.h"\r
-\r
-#define TRACE 0\r
-\r
-// Return false when done\r
-bool PhyEnumEdges(const Tree &tree, PhyEnumEdgeState &ES)\r
- {\r
- unsigned uNode1 = uInsane;\r
-\r
- if (!ES.m_bInit)\r
- {\r
- if (tree.GetNodeCount() <= 1)\r
- {\r
- ES.m_uNodeIndex1 = NULL_NEIGHBOR;\r
- ES.m_uNodeIndex2 = NULL_NEIGHBOR;\r
- return false;\r
- }\r
- uNode1 = tree.FirstDepthFirstNode();\r
- ES.m_bInit = true;\r
- }\r
- else\r
- {\r
- uNode1 = tree.NextDepthFirstNode(ES.m_uNodeIndex1);\r
- if (NULL_NEIGHBOR == uNode1)\r
- return false;\r
- if (tree.IsRooted() && tree.IsRoot(uNode1))\r
- {\r
- uNode1 = tree.NextDepthFirstNode(uNode1);\r
- if (NULL_NEIGHBOR == uNode1)\r
- return false;\r
- }\r
- }\r
- unsigned uNode2 = tree.GetParent(uNode1);\r
-\r
- ES.m_uNodeIndex1 = uNode1;\r
- ES.m_uNodeIndex2 = uNode2;\r
- return true;\r
- }\r
-\r
-bool PhyEnumEdgesR(const Tree &tree, PhyEnumEdgeState &ES)\r
- {\r
- unsigned uNode1 = uInsane;\r
-\r
- if (!ES.m_bInit)\r
- {\r
- if (tree.GetNodeCount() <= 1)\r
- {\r
- ES.m_uNodeIndex1 = NULL_NEIGHBOR;\r
- ES.m_uNodeIndex2 = NULL_NEIGHBOR;\r
- return false;\r
- }\r
- uNode1 = tree.FirstDepthFirstNodeR();\r
- ES.m_bInit = true;\r
- }\r
- else\r
- {\r
- uNode1 = tree.NextDepthFirstNodeR(ES.m_uNodeIndex1);\r
- if (NULL_NEIGHBOR == uNode1)\r
- return false;\r
- if (tree.IsRooted() && tree.IsRoot(uNode1))\r
- {\r
- uNode1 = tree.NextDepthFirstNode(uNode1);\r
- if (NULL_NEIGHBOR == uNode1)\r
- return false;\r
- }\r
- }\r
- unsigned uNode2 = tree.GetParent(uNode1);\r
-\r
- ES.m_uNodeIndex1 = uNode1;\r
- ES.m_uNodeIndex2 = uNode2;\r
- return true;\r
- }\r
-\r
-static void GetLeavesSubtree(const Tree &tree, unsigned uNodeIndex1,\r
- const unsigned uNodeIndex2, unsigned Leaves[], unsigned *ptruCount)\r
- {\r
- if (tree.IsLeaf(uNodeIndex1))\r
- {\r
- Leaves[*ptruCount] = uNodeIndex1;\r
- ++(*ptruCount);\r
- return;\r
- }\r
-\r
- const unsigned uLeft = tree.GetFirstNeighbor(uNodeIndex1, uNodeIndex2);\r
- const unsigned uRight = tree.GetSecondNeighbor(uNodeIndex1, uNodeIndex2);\r
- if (NULL_NEIGHBOR != uLeft)\r
- GetLeavesSubtree(tree, uLeft, uNodeIndex1, Leaves, ptruCount);\r
- if (NULL_NEIGHBOR != uRight)\r
- GetLeavesSubtree(tree, uRight, uNodeIndex1, Leaves, ptruCount);\r
- }\r
-\r
-static void PhyGetLeaves(const Tree &tree, unsigned uNodeIndex1, unsigned uNodeIndex2,\r
- unsigned Leaves[], unsigned *ptruCount)\r
- {\r
- *ptruCount = 0;\r
- GetLeavesSubtree(tree, uNodeIndex1, uNodeIndex2, Leaves, ptruCount);\r
- }\r
-\r
-bool PhyEnumBiParts(const Tree &tree, PhyEnumEdgeState &ES,\r
- unsigned Leaves1[], unsigned *ptruCount1,\r
- unsigned Leaves2[], unsigned *ptruCount2)\r
- {\r
- bool bOk = PhyEnumEdges(tree, ES);\r
- if (!bOk)\r
- {\r
- *ptruCount1 = 0;\r
- *ptruCount2 = 0;\r
- return false;\r
- }\r
-\r
-// Special case: in a rooted tree, both edges from the root\r
-// give the same bipartition, so skip one of them.\r
- if (tree.IsRooted() && tree.IsRoot(ES.m_uNodeIndex2)\r
- && tree.GetRight(ES.m_uNodeIndex2) == ES.m_uNodeIndex1)\r
- {\r
- bOk = PhyEnumEdges(tree, ES);\r
- if (!bOk)\r
- return false;\r
- }\r
-\r
- PhyGetLeaves(tree, ES.m_uNodeIndex1, ES.m_uNodeIndex2, Leaves1, ptruCount1);\r
- PhyGetLeaves(tree, ES.m_uNodeIndex2, ES.m_uNodeIndex1, Leaves2, ptruCount2);\r
-\r
- if (*ptruCount1 + *ptruCount2 != tree.GetLeafCount())\r
- Quit("PhyEnumBiParts %u + %u != %u",\r
- *ptruCount1, *ptruCount2, tree.GetLeafCount());\r
-#if DEBUG\r
- {\r
- for (unsigned i = 0; i < *ptruCount1; ++i)\r
- {\r
- if (!tree.IsLeaf(Leaves1[i]))\r
- Quit("PhyEnumByParts: not leaf");\r
- for (unsigned j = 0; j < *ptruCount2; ++j)\r
- {\r
- if (!tree.IsLeaf(Leaves2[j]))\r
- Quit("PhyEnumByParts: not leaf");\r
- if (Leaves1[i] == Leaves2[j])\r
- Quit("PhyEnumByParts: dupe");\r
- }\r
- }\r
- }\r
-#endif\r
-\r
- return true;\r
- }\r
-\r
-#if 0\r
-void TestBiPart()\r
- {\r
- SetListFileName("c:\\tmp\\lobster.log", false);\r
- Tree tree;\r
- TextFile fileIn("c:\\tmp\\test.phy");\r
- tree.FromFile(fileIn);\r
- tree.LogMe();\r
-\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- unsigned *Leaves1 = new unsigned[uNodeCount];\r
- unsigned *Leaves2 = new unsigned[uNodeCount];\r
-\r
- PhyEnumEdgeState ES;\r
- bool bDone = false;\r
- for (;;)\r
- {\r
- unsigned uCount1 = uInsane;\r
- unsigned uCount2 = uInsane;\r
- bool bOk = PhyEnumBiParts(tree, ES, Leaves1, &uCount1, Leaves2, &uCount2);\r
- Log("PEBP=%d ES.Init=%d ES.ni1=%d ES.ni2=%d\n",\r
- bOk,\r
- ES.m_bInit,\r
- ES.m_uNodeIndex1,\r
- ES.m_uNodeIndex2);\r
- if (!bOk)\r
- break;\r
- Log("\n");\r
- Log("Part1: ");\r
- for (unsigned n = 0; n < uCount1; ++n)\r
- Log(" %d(%s)", Leaves1[n], tree.GetLeafName(Leaves1[n]));\r
- Log("\n");\r
- Log("Part2: ");\r
- for (unsigned n = 0; n < uCount2; ++n)\r
- Log(" %d(%s)", Leaves2[n], tree.GetLeafName(Leaves2[n]));\r
- Log("\n");\r
- }\r
- }\r
-#endif\r
-\r
-static void GetLeavesSubtreeExcluding(const Tree &tree, unsigned uNodeIndex,\r
- unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)\r
- {\r
- if (uNodeIndex == uExclude)\r
- return;\r
-\r
- if (tree.IsLeaf(uNodeIndex))\r
- {\r
- Leaves[*ptruCount] = uNodeIndex;\r
- ++(*ptruCount);\r
- return;\r
- }\r
-\r
- const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
- const unsigned uRight = tree.GetRight(uNodeIndex);\r
- if (NULL_NEIGHBOR != uLeft)\r
- GetLeavesSubtreeExcluding(tree, uLeft, uExclude, Leaves, ptruCount);\r
- if (NULL_NEIGHBOR != uRight)\r
- GetLeavesSubtreeExcluding(tree, uRight, uExclude, Leaves, ptruCount);\r
- }\r
-\r
-void GetLeavesExcluding(const Tree &tree, unsigned uNodeIndex,\r
- unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)\r
- {\r
- *ptruCount = 0;\r
- GetLeavesSubtreeExcluding(tree, uNodeIndex, uExclude, Leaves, ptruCount);\r
- }\r
-\r
-void GetInternalNodesInHeightOrder(const Tree &tree, unsigned NodeIndexes[])\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- if (uNodeCount < 3)\r
- Quit("GetInternalNodesInHeightOrder: %u nodes, none are internal",\r
- uNodeCount);\r
- const unsigned uInternalNodeCount = (uNodeCount - 1)/2;\r
- double *Heights = new double[uInternalNodeCount];\r
-\r
- unsigned uIndex = 0;\r
- for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
- {\r
- if (tree.IsLeaf(uNodeIndex))\r
- continue;\r
- NodeIndexes[uIndex] = uNodeIndex;\r
- Heights[uIndex] = tree.GetNodeHeight(uNodeIndex);\r
- ++uIndex;\r
- }\r
- if (uIndex != uInternalNodeCount)\r
- Quit("Internal error: GetInternalNodesInHeightOrder");\r
-\r
-// Simple but slow bubble sort (probably don't care about speed here)\r
- bool bDone = false;\r
- while (!bDone)\r
- {\r
- bDone = true;\r
- for (unsigned i = 0; i < uInternalNodeCount - 1; ++i)\r
- {\r
- if (Heights[i] > Heights[i+1])\r
- {\r
- double dTmp = Heights[i];\r
- Heights[i] = Heights[i+1];\r
- Heights[i+1] = dTmp;\r
-\r
- unsigned uTmp = NodeIndexes[i];\r
- NodeIndexes[i] = NodeIndexes[i+1];\r
- NodeIndexes[i+1] = uTmp;\r
- bDone = false;\r
- }\r
- }\r
- }\r
-#if TRACE\r
- Log("Internal node index Height\n");\r
- Log("------------------- --------\n");\r
- // 1234567890123456789 123456789\r
- for (unsigned n = 0; n < uInternalNodeCount; ++n)\r
- Log("%19u %9.3f\n", NodeIndexes[n], Heights[n]);\r
-#endif\r
- delete[] Heights;\r
- }\r
-\r
-void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)\r
- {\r
- const unsigned uNodeCount = tree.GetNodeCount();\r
- for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
- {\r
- const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);\r
- for (unsigned n = 0; n < uNeighborCount; ++n)\r
- {\r
- const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);\r
- if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))\r
- continue;\r
- if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)\r
- tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);\r
- }\r
- }\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff