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