6 // Return false when done
\r
7 bool PhyEnumEdges(const Tree &tree, PhyEnumEdgeState &ES)
\r
9 unsigned uNode1 = uInsane;
\r
13 if (tree.GetNodeCount() <= 1)
\r
15 ES.m_uNodeIndex1 = NULL_NEIGHBOR;
\r
16 ES.m_uNodeIndex2 = NULL_NEIGHBOR;
\r
19 uNode1 = tree.FirstDepthFirstNode();
\r
24 uNode1 = tree.NextDepthFirstNode(ES.m_uNodeIndex1);
\r
25 if (NULL_NEIGHBOR == uNode1)
\r
27 if (tree.IsRooted() && tree.IsRoot(uNode1))
\r
29 uNode1 = tree.NextDepthFirstNode(uNode1);
\r
30 if (NULL_NEIGHBOR == uNode1)
\r
34 unsigned uNode2 = tree.GetParent(uNode1);
\r
36 ES.m_uNodeIndex1 = uNode1;
\r
37 ES.m_uNodeIndex2 = uNode2;
\r
41 bool PhyEnumEdgesR(const Tree &tree, PhyEnumEdgeState &ES)
\r
43 unsigned uNode1 = uInsane;
\r
47 if (tree.GetNodeCount() <= 1)
\r
49 ES.m_uNodeIndex1 = NULL_NEIGHBOR;
\r
50 ES.m_uNodeIndex2 = NULL_NEIGHBOR;
\r
53 uNode1 = tree.FirstDepthFirstNodeR();
\r
58 uNode1 = tree.NextDepthFirstNodeR(ES.m_uNodeIndex1);
\r
59 if (NULL_NEIGHBOR == uNode1)
\r
61 if (tree.IsRooted() && tree.IsRoot(uNode1))
\r
63 uNode1 = tree.NextDepthFirstNode(uNode1);
\r
64 if (NULL_NEIGHBOR == uNode1)
\r
68 unsigned uNode2 = tree.GetParent(uNode1);
\r
70 ES.m_uNodeIndex1 = uNode1;
\r
71 ES.m_uNodeIndex2 = uNode2;
\r
75 static void GetLeavesSubtree(const Tree &tree, unsigned uNodeIndex1,
\r
76 const unsigned uNodeIndex2, unsigned Leaves[], unsigned *ptruCount)
\r
78 if (tree.IsLeaf(uNodeIndex1))
\r
80 Leaves[*ptruCount] = uNodeIndex1;
\r
85 const unsigned uLeft = tree.GetFirstNeighbor(uNodeIndex1, uNodeIndex2);
\r
86 const unsigned uRight = tree.GetSecondNeighbor(uNodeIndex1, uNodeIndex2);
\r
87 if (NULL_NEIGHBOR != uLeft)
\r
88 GetLeavesSubtree(tree, uLeft, uNodeIndex1, Leaves, ptruCount);
\r
89 if (NULL_NEIGHBOR != uRight)
\r
90 GetLeavesSubtree(tree, uRight, uNodeIndex1, Leaves, ptruCount);
\r
93 static void PhyGetLeaves(const Tree &tree, unsigned uNodeIndex1, unsigned uNodeIndex2,
\r
94 unsigned Leaves[], unsigned *ptruCount)
\r
97 GetLeavesSubtree(tree, uNodeIndex1, uNodeIndex2, Leaves, ptruCount);
\r
100 bool PhyEnumBiParts(const Tree &tree, PhyEnumEdgeState &ES,
\r
101 unsigned Leaves1[], unsigned *ptruCount1,
\r
102 unsigned Leaves2[], unsigned *ptruCount2)
\r
104 bool bOk = PhyEnumEdges(tree, ES);
\r
112 // Special case: in a rooted tree, both edges from the root
\r
113 // give the same bipartition, so skip one of them.
\r
114 if (tree.IsRooted() && tree.IsRoot(ES.m_uNodeIndex2)
\r
115 && tree.GetRight(ES.m_uNodeIndex2) == ES.m_uNodeIndex1)
\r
117 bOk = PhyEnumEdges(tree, ES);
\r
122 PhyGetLeaves(tree, ES.m_uNodeIndex1, ES.m_uNodeIndex2, Leaves1, ptruCount1);
\r
123 PhyGetLeaves(tree, ES.m_uNodeIndex2, ES.m_uNodeIndex1, Leaves2, ptruCount2);
\r
125 if (*ptruCount1 + *ptruCount2 != tree.GetLeafCount())
\r
126 Quit("PhyEnumBiParts %u + %u != %u",
\r
127 *ptruCount1, *ptruCount2, tree.GetLeafCount());
\r
130 for (unsigned i = 0; i < *ptruCount1; ++i)
\r
132 if (!tree.IsLeaf(Leaves1[i]))
\r
133 Quit("PhyEnumByParts: not leaf");
\r
134 for (unsigned j = 0; j < *ptruCount2; ++j)
\r
136 if (!tree.IsLeaf(Leaves2[j]))
\r
137 Quit("PhyEnumByParts: not leaf");
\r
138 if (Leaves1[i] == Leaves2[j])
\r
139 Quit("PhyEnumByParts: dupe");
\r
151 SetListFileName("c:\\tmp\\lobster.log", false);
\r
153 TextFile fileIn("c:\\tmp\\test.phy");
\r
154 tree.FromFile(fileIn);
\r
157 const unsigned uNodeCount = tree.GetNodeCount();
\r
158 unsigned *Leaves1 = new unsigned[uNodeCount];
\r
159 unsigned *Leaves2 = new unsigned[uNodeCount];
\r
161 PhyEnumEdgeState ES;
\r
162 bool bDone = false;
\r
165 unsigned uCount1 = uInsane;
\r
166 unsigned uCount2 = uInsane;
\r
167 bool bOk = PhyEnumBiParts(tree, ES, Leaves1, &uCount1, Leaves2, &uCount2);
\r
168 Log("PEBP=%d ES.Init=%d ES.ni1=%d ES.ni2=%d\n",
\r
177 for (unsigned n = 0; n < uCount1; ++n)
\r
178 Log(" %d(%s)", Leaves1[n], tree.GetLeafName(Leaves1[n]));
\r
181 for (unsigned n = 0; n < uCount2; ++n)
\r
182 Log(" %d(%s)", Leaves2[n], tree.GetLeafName(Leaves2[n]));
\r
188 static void GetLeavesSubtreeExcluding(const Tree &tree, unsigned uNodeIndex,
\r
189 unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
\r
191 if (uNodeIndex == uExclude)
\r
194 if (tree.IsLeaf(uNodeIndex))
\r
196 Leaves[*ptruCount] = uNodeIndex;
\r
201 const unsigned uLeft = tree.GetLeft(uNodeIndex);
\r
202 const unsigned uRight = tree.GetRight(uNodeIndex);
\r
203 if (NULL_NEIGHBOR != uLeft)
\r
204 GetLeavesSubtreeExcluding(tree, uLeft, uExclude, Leaves, ptruCount);
\r
205 if (NULL_NEIGHBOR != uRight)
\r
206 GetLeavesSubtreeExcluding(tree, uRight, uExclude, Leaves, ptruCount);
\r
209 void GetLeavesExcluding(const Tree &tree, unsigned uNodeIndex,
\r
210 unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
\r
213 GetLeavesSubtreeExcluding(tree, uNodeIndex, uExclude, Leaves, ptruCount);
\r
216 void GetInternalNodesInHeightOrder(const Tree &tree, unsigned NodeIndexes[])
\r
218 const unsigned uNodeCount = tree.GetNodeCount();
\r
219 if (uNodeCount < 3)
\r
220 Quit("GetInternalNodesInHeightOrder: %u nodes, none are internal",
\r
222 const unsigned uInternalNodeCount = (uNodeCount - 1)/2;
\r
223 double *Heights = new double[uInternalNodeCount];
\r
225 unsigned uIndex = 0;
\r
226 for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
\r
228 if (tree.IsLeaf(uNodeIndex))
\r
230 NodeIndexes[uIndex] = uNodeIndex;
\r
231 Heights[uIndex] = tree.GetNodeHeight(uNodeIndex);
\r
234 if (uIndex != uInternalNodeCount)
\r
235 Quit("Internal error: GetInternalNodesInHeightOrder");
\r
237 // Simple but slow bubble sort (probably don't care about speed here)
\r
238 bool bDone = false;
\r
242 for (unsigned i = 0; i < uInternalNodeCount - 1; ++i)
\r
244 if (Heights[i] > Heights[i+1])
\r
246 double dTmp = Heights[i];
\r
247 Heights[i] = Heights[i+1];
\r
248 Heights[i+1] = dTmp;
\r
250 unsigned uTmp = NodeIndexes[i];
\r
251 NodeIndexes[i] = NodeIndexes[i+1];
\r
252 NodeIndexes[i+1] = uTmp;
\r
258 Log("Internal node index Height\n");
\r
259 Log("------------------- --------\n");
\r
260 // 1234567890123456789 123456789
\r
261 for (unsigned n = 0; n < uInternalNodeCount; ++n)
\r
262 Log("%19u %9.3f\n", NodeIndexes[n], Heights[n]);
\r
267 void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)
\r
269 const unsigned uNodeCount = tree.GetNodeCount();
\r
270 for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
\r
272 const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);
\r
273 for (unsigned n = 0; n < uNeighborCount; ++n)
\r
275 const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);
\r
276 if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))
\r
278 if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)
\r
279 tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);
\r