8 Node has 0 to 3 neighbors:
\r
9 0 neighbors: singleton root
\r
10 1 neighbor: leaf, neighbor is parent
\r
11 2 neigbors: non-singleton root
\r
12 3 neighbors: internal node (other than root)
\r
14 Minimal rooted tree is single node.
\r
15 Minimal unrooted tree is single edge.
\r
16 Leaf node always has nulls in neighbors 2 and 3, neighbor 1 is parent.
\r
17 When tree is rooted, neighbor 1=parent, 2=left, 3=right.
\r
20 void Tree::AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2) const
\r
22 if (uNodeIndex1 >= m_uNodeCount || uNodeIndex2 >= m_uNodeCount)
\r
23 Quit("AssertAreNeighbors(%u,%u), are %u nodes",
\r
24 uNodeIndex1, uNodeIndex2, m_uNodeCount);
\r
26 if (m_uNeighbor1[uNodeIndex1] != uNodeIndex2 &&
\r
27 m_uNeighbor2[uNodeIndex1] != uNodeIndex2 &&
\r
28 m_uNeighbor3[uNodeIndex1] != uNodeIndex2)
\r
31 Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
\r
34 if (m_uNeighbor1[uNodeIndex2] != uNodeIndex1 &&
\r
35 m_uNeighbor2[uNodeIndex2] != uNodeIndex1 &&
\r
36 m_uNeighbor3[uNodeIndex2] != uNodeIndex1)
\r
39 Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
\r
42 bool Has12 = HasEdgeLength(uNodeIndex1, uNodeIndex2);
\r
43 bool Has21 = HasEdgeLength(uNodeIndex2, uNodeIndex1);
\r
46 HasEdgeLength(uNodeIndex1, uNodeIndex2);
\r
47 HasEdgeLength(uNodeIndex2, uNodeIndex1);
\r
49 Log("HasEdgeLength(%u, %u)=%c HasEdgeLength(%u, %u)=%c\n",
\r
57 Quit("Tree::AssertAreNeighbors, HasEdgeLength not symmetric");
\r
62 double d12 = GetEdgeLength(uNodeIndex1, uNodeIndex2);
\r
63 double d21 = GetEdgeLength(uNodeIndex2, uNodeIndex1);
\r
67 Quit("Tree::AssertAreNeighbors, Edge length disagrees %u-%u=%.3g, %u-%u=%.3g",
\r
68 uNodeIndex1, uNodeIndex2, d12,
\r
69 uNodeIndex2, uNodeIndex1, d21);
\r
74 void Tree::ValidateNode(unsigned uNodeIndex) const
\r
76 if (uNodeIndex >= m_uNodeCount)
\r
77 Quit("ValidateNode(%u), %u nodes", uNodeIndex, m_uNodeCount);
\r
79 const unsigned uNeighborCount = GetNeighborCount(uNodeIndex);
\r
81 if (2 == uNeighborCount)
\r
86 Quit("Tree::ValidateNode: Node %u has two neighbors, tree is not rooted",
\r
89 if (uNodeIndex != m_uRootNodeIndex)
\r
92 Quit("Tree::ValidateNode: Node %u has two neighbors, but not root node=%u",
\r
93 uNodeIndex, m_uRootNodeIndex);
\r
97 const unsigned n1 = m_uNeighbor1[uNodeIndex];
\r
98 const unsigned n2 = m_uNeighbor2[uNodeIndex];
\r
99 const unsigned n3 = m_uNeighbor3[uNodeIndex];
\r
101 if (NULL_NEIGHBOR == n2 && NULL_NEIGHBOR != n3)
\r
104 Quit("Tree::ValidateNode, n2=null, n3!=null", uNodeIndex);
\r
106 if (NULL_NEIGHBOR == n3 && NULL_NEIGHBOR != n2)
\r
109 Quit("Tree::ValidateNode, n3=null, n2!=null", uNodeIndex);
\r
112 if (n1 != NULL_NEIGHBOR)
\r
113 AssertAreNeighbors(uNodeIndex, n1);
\r
114 if (n2 != NULL_NEIGHBOR)
\r
115 AssertAreNeighbors(uNodeIndex, n2);
\r
116 if (n3 != NULL_NEIGHBOR)
\r
117 AssertAreNeighbors(uNodeIndex, n3);
\r
119 if (n1 != NULL_NEIGHBOR && (n1 == n2 || n1 == n3))
\r
122 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
\r
124 if (n2 != NULL_NEIGHBOR && (n2 == n1 || n2 == n3))
\r
127 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
\r
129 if (n3 != NULL_NEIGHBOR && (n3 == n1 || n3 == n2))
\r
132 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
\r
137 if (NULL_NEIGHBOR == GetParent(uNodeIndex))
\r
139 if (uNodeIndex != m_uRootNodeIndex)
\r
142 Quit("Tree::ValiateNode(%u), no parent", uNodeIndex);
\r
145 else if (GetLeft(GetParent(uNodeIndex)) != uNodeIndex &&
\r
146 GetRight(GetParent(uNodeIndex)) != uNodeIndex)
\r
149 Quit("Tree::ValidateNode(%u), parent / child mismatch", uNodeIndex);
\r
154 void Tree::Validate() const
\r
156 for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
\r
157 ValidateNode(uNodeIndex);
\r
160 bool Tree::IsEdge(unsigned uNodeIndex1, unsigned uNodeIndex2) const
\r
162 assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
\r
164 return m_uNeighbor1[uNodeIndex1] == uNodeIndex2 ||
\r
165 m_uNeighbor2[uNodeIndex1] == uNodeIndex2 ||
\r
166 m_uNeighbor3[uNodeIndex1] == uNodeIndex2;
\r
169 double Tree::GetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
\r
171 assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
\r
172 if (!HasEdgeLength(uNodeIndex1, uNodeIndex2))
\r
175 Quit("Missing edge length in tree %u-%u", uNodeIndex1, uNodeIndex2);
\r
178 if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
\r
179 return m_dEdgeLength1[uNodeIndex1];
\r
180 else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
\r
181 return m_dEdgeLength2[uNodeIndex1];
\r
182 assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
\r
183 return m_dEdgeLength3[uNodeIndex1];
\r
186 void Tree::ExpandCache()
\r
188 const unsigned uNodeCount = 100;
\r
189 unsigned uNewCacheCount = m_uCacheCount + uNodeCount;
\r
190 unsigned *uNewNeighbor1 = new unsigned[uNewCacheCount];
\r
191 unsigned *uNewNeighbor2 = new unsigned[uNewCacheCount];
\r
192 unsigned *uNewNeighbor3 = new unsigned[uNewCacheCount];
\r
194 unsigned *uNewIds = new unsigned[uNewCacheCount];
\r
195 memset(uNewIds, 0xff, uNewCacheCount*sizeof(unsigned));
\r
197 double *dNewEdgeLength1 = new double[uNewCacheCount];
\r
198 double *dNewEdgeLength2 = new double[uNewCacheCount];
\r
199 double *dNewEdgeLength3 = new double[uNewCacheCount];
\r
200 double *dNewHeight = new double[uNewCacheCount];
\r
202 bool *bNewHasEdgeLength1 = new bool[uNewCacheCount];
\r
203 bool *bNewHasEdgeLength2 = new bool[uNewCacheCount];
\r
204 bool *bNewHasEdgeLength3 = new bool[uNewCacheCount];
\r
205 bool *bNewHasHeight = new bool[uNewCacheCount];
\r
207 char **ptrNewName = new char *[uNewCacheCount];
\r
208 memset(ptrNewName, 0, uNewCacheCount*sizeof(char *));
\r
210 if (m_uCacheCount > 0)
\r
212 const unsigned uUnsignedBytes = m_uCacheCount*sizeof(unsigned);
\r
213 memcpy(uNewNeighbor1, m_uNeighbor1, uUnsignedBytes);
\r
214 memcpy(uNewNeighbor2, m_uNeighbor2, uUnsignedBytes);
\r
215 memcpy(uNewNeighbor3, m_uNeighbor3, uUnsignedBytes);
\r
217 memcpy(uNewIds, m_Ids, uUnsignedBytes);
\r
219 const unsigned uEdgeBytes = m_uCacheCount*sizeof(double);
\r
220 memcpy(dNewEdgeLength1, m_dEdgeLength1, uEdgeBytes);
\r
221 memcpy(dNewEdgeLength2, m_dEdgeLength2, uEdgeBytes);
\r
222 memcpy(dNewEdgeLength3, m_dEdgeLength3, uEdgeBytes);
\r
223 memcpy(dNewHeight, m_dHeight, uEdgeBytes);
\r
225 const unsigned uBoolBytes = m_uCacheCount*sizeof(bool);
\r
226 memcpy(bNewHasEdgeLength1, m_bHasEdgeLength1, uBoolBytes);
\r
227 memcpy(bNewHasEdgeLength2, m_bHasEdgeLength2, uBoolBytes);
\r
228 memcpy(bNewHasEdgeLength3, m_bHasEdgeLength3, uBoolBytes);
\r
229 memcpy(bNewHasHeight, m_bHasHeight, uBoolBytes);
\r
231 const unsigned uNameBytes = m_uCacheCount*sizeof(char *);
\r
232 memcpy(ptrNewName, m_ptrName, uNameBytes);
\r
234 delete[] m_uNeighbor1;
\r
235 delete[] m_uNeighbor2;
\r
236 delete[] m_uNeighbor3;
\r
240 delete[] m_dEdgeLength1;
\r
241 delete[] m_dEdgeLength2;
\r
242 delete[] m_dEdgeLength3;
\r
244 delete[] m_bHasEdgeLength1;
\r
245 delete[] m_bHasEdgeLength2;
\r
246 delete[] m_bHasEdgeLength3;
\r
247 delete[] m_bHasHeight;
\r
249 delete[] m_ptrName;
\r
251 m_uCacheCount = uNewCacheCount;
\r
252 m_uNeighbor1 = uNewNeighbor1;
\r
253 m_uNeighbor2 = uNewNeighbor2;
\r
254 m_uNeighbor3 = uNewNeighbor3;
\r
256 m_dEdgeLength1 = dNewEdgeLength1;
\r
257 m_dEdgeLength2 = dNewEdgeLength2;
\r
258 m_dEdgeLength3 = dNewEdgeLength3;
\r
259 m_dHeight = dNewHeight;
\r
260 m_bHasEdgeLength1 = bNewHasEdgeLength1;
\r
261 m_bHasEdgeLength2 = bNewHasEdgeLength2;
\r
262 m_bHasEdgeLength3 = bNewHasEdgeLength3;
\r
263 m_bHasHeight = bNewHasHeight;
\r
264 m_ptrName = ptrNewName;
\r
267 // Creates tree with single node, no edges.
\r
268 // Root node always has index 0.
\r
269 void Tree::CreateRooted()
\r
275 m_uNeighbor1[0] = NULL_NEIGHBOR;
\r
276 m_uNeighbor2[0] = NULL_NEIGHBOR;
\r
277 m_uNeighbor3[0] = NULL_NEIGHBOR;
\r
279 m_bHasEdgeLength1[0] = false;
\r
280 m_bHasEdgeLength2[0] = false;
\r
281 m_bHasEdgeLength3[0] = false;
\r
282 m_bHasHeight[0] = false;
\r
284 m_uRootNodeIndex = 0;
\r
292 // Creates unrooted tree with single edge.
\r
293 // Nodes for that edge are always 0 and 1.
\r
294 void Tree::CreateUnrooted(double dEdgeLength)
\r
299 m_uNeighbor1[0] = 1;
\r
300 m_uNeighbor2[0] = NULL_NEIGHBOR;
\r
301 m_uNeighbor3[0] = NULL_NEIGHBOR;
\r
303 m_uNeighbor1[1] = 0;
\r
304 m_uNeighbor2[1] = NULL_NEIGHBOR;
\r
305 m_uNeighbor3[1] = NULL_NEIGHBOR;
\r
307 m_dEdgeLength1[0] = dEdgeLength;
\r
308 m_dEdgeLength1[1] = dEdgeLength;
\r
310 m_bHasEdgeLength1[0] = true;
\r
311 m_bHasEdgeLength1[1] = true;
\r
320 void Tree::SetLeafName(unsigned uNodeIndex, const char *ptrName)
\r
322 assert(uNodeIndex < m_uNodeCount);
\r
323 assert(IsLeaf(uNodeIndex));
\r
324 free(m_ptrName[uNodeIndex]);
\r
325 m_ptrName[uNodeIndex] = strsave(ptrName);
\r
328 void Tree::SetLeafId(unsigned uNodeIndex, unsigned uId)
\r
330 assert(uNodeIndex < m_uNodeCount);
\r
331 assert(IsLeaf(uNodeIndex));
\r
332 m_Ids[uNodeIndex] = uId;
\r
335 const char *Tree::GetLeafName(unsigned uNodeIndex) const
\r
337 assert(uNodeIndex < m_uNodeCount);
\r
338 assert(IsLeaf(uNodeIndex));
\r
339 return m_ptrName[uNodeIndex];
\r
342 unsigned Tree::GetLeafId(unsigned uNodeIndex) const
\r
344 assert(uNodeIndex < m_uNodeCount);
\r
345 assert(IsLeaf(uNodeIndex));
\r
346 return m_Ids[uNodeIndex];
\r
349 // Append a new branch.
\r
350 // This adds two new nodes and joins them to an existing leaf node.
\r
351 // Return value is k, new nodes have indexes k and k+1 respectively.
\r
352 unsigned Tree::AppendBranch(unsigned uExistingLeafIndex)
\r
354 if (0 == m_uNodeCount)
\r
355 Quit("Tree::AppendBranch: tree has not been created");
\r
358 assert(uExistingLeafIndex < m_uNodeCount);
\r
359 if (!IsLeaf(uExistingLeafIndex))
\r
362 Quit("AppendBranch(%u): not leaf", uExistingLeafIndex);
\r
366 if (m_uNodeCount >= m_uCacheCount - 2)
\r
369 const unsigned uNewLeaf1 = m_uNodeCount;
\r
370 const unsigned uNewLeaf2 = m_uNodeCount + 1;
\r
374 assert(m_uNeighbor2[uExistingLeafIndex] == NULL_NEIGHBOR);
\r
375 assert(m_uNeighbor3[uExistingLeafIndex] == NULL_NEIGHBOR);
\r
377 m_uNeighbor2[uExistingLeafIndex] = uNewLeaf1;
\r
378 m_uNeighbor3[uExistingLeafIndex] = uNewLeaf2;
\r
380 m_uNeighbor1[uNewLeaf1] = uExistingLeafIndex;
\r
381 m_uNeighbor1[uNewLeaf2] = uExistingLeafIndex;
\r
383 m_uNeighbor2[uNewLeaf1] = NULL_NEIGHBOR;
\r
384 m_uNeighbor2[uNewLeaf2] = NULL_NEIGHBOR;
\r
386 m_uNeighbor3[uNewLeaf1] = NULL_NEIGHBOR;
\r
387 m_uNeighbor3[uNewLeaf2] = NULL_NEIGHBOR;
\r
389 m_dEdgeLength2[uExistingLeafIndex] = 0;
\r
390 m_dEdgeLength3[uExistingLeafIndex] = 0;
\r
392 m_dEdgeLength1[uNewLeaf1] = 0;
\r
393 m_dEdgeLength2[uNewLeaf1] = 0;
\r
394 m_dEdgeLength3[uNewLeaf1] = 0;
\r
396 m_dEdgeLength1[uNewLeaf2] = 0;
\r
397 m_dEdgeLength2[uNewLeaf2] = 0;
\r
398 m_dEdgeLength3[uNewLeaf2] = 0;
\r
400 m_bHasEdgeLength1[uNewLeaf1] = false;
\r
401 m_bHasEdgeLength2[uNewLeaf1] = false;
\r
402 m_bHasEdgeLength3[uNewLeaf1] = false;
\r
404 m_bHasEdgeLength1[uNewLeaf2] = false;
\r
405 m_bHasEdgeLength2[uNewLeaf2] = false;
\r
406 m_bHasEdgeLength3[uNewLeaf2] = false;
\r
408 m_bHasHeight[uNewLeaf1] = false;
\r
409 m_bHasHeight[uNewLeaf2] = false;
\r
411 m_Ids[uNewLeaf1] = uInsane;
\r
412 m_Ids[uNewLeaf2] = uInsane;
\r
416 void Tree::LogMe() const
\r
418 Log("Tree::LogMe %u nodes, ", m_uNodeCount);
\r
424 Log("Index Parnt LengthP Left LengthL Right LengthR Id Name\n");
\r
425 Log("----- ----- ------- ---- ------- ----- ------- ----- ----\n");
\r
429 Log("unrooted.\n");
\r
431 Log("Index Nbr_1 Length1 Nbr_2 Length2 Nbr_3 Length3 Id Name\n");
\r
432 Log("----- ----- ------- ----- ------- ----- ------- ----- ----\n");
\r
435 for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
\r
437 Log("%5u ", uNodeIndex);
\r
438 const unsigned n1 = m_uNeighbor1[uNodeIndex];
\r
439 const unsigned n2 = m_uNeighbor2[uNodeIndex];
\r
440 const unsigned n3 = m_uNeighbor3[uNodeIndex];
\r
441 if (NULL_NEIGHBOR != n1)
\r
444 if (m_bHasEdgeLength1[uNodeIndex])
\r
445 Log("%7.3g ", m_dEdgeLength1[uNodeIndex]);
\r
452 if (NULL_NEIGHBOR != n2)
\r
455 if (m_bHasEdgeLength2[uNodeIndex])
\r
456 Log("%7.3g ", m_dEdgeLength2[uNodeIndex]);
\r
463 if (NULL_NEIGHBOR != n3)
\r
466 if (m_bHasEdgeLength3[uNodeIndex])
\r
467 Log("%7.3g ", m_dEdgeLength3[uNodeIndex]);
\r
474 if (m_Ids != 0 && IsLeaf(uNodeIndex))
\r
476 unsigned uId = m_Ids[uNodeIndex];
\r
477 if (uId == uInsane)
\r
485 if (m_bRooted && uNodeIndex == m_uRootNodeIndex)
\r
487 const char *ptrName = m_ptrName[uNodeIndex];
\r
489 Log(" %s", ptrName);
\r
494 void Tree::SetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2,
\r
497 assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
\r
498 assert(IsEdge(uNodeIndex1, uNodeIndex2));
\r
500 if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
\r
502 m_dEdgeLength1[uNodeIndex1] = dLength;
\r
503 m_bHasEdgeLength1[uNodeIndex1] = true;
\r
505 else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
\r
507 m_dEdgeLength2[uNodeIndex1] = dLength;
\r
508 m_bHasEdgeLength2[uNodeIndex1] = true;
\r
513 assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
\r
514 m_dEdgeLength3[uNodeIndex1] = dLength;
\r
515 m_bHasEdgeLength3[uNodeIndex1] = true;
\r
518 if (m_uNeighbor1[uNodeIndex2] == uNodeIndex1)
\r
520 m_dEdgeLength1[uNodeIndex2] = dLength;
\r
521 m_bHasEdgeLength1[uNodeIndex2] = true;
\r
523 else if (m_uNeighbor2[uNodeIndex2] == uNodeIndex1)
\r
525 m_dEdgeLength2[uNodeIndex2] = dLength;
\r
526 m_bHasEdgeLength2[uNodeIndex2] = true;
\r
530 assert(m_uNeighbor3[uNodeIndex2] == uNodeIndex1);
\r
531 m_dEdgeLength3[uNodeIndex2] = dLength;
\r
532 m_bHasEdgeLength3[uNodeIndex2] = true;
\r
536 unsigned Tree::UnrootFromFile()
\r
539 Log("Before unroot:\n");
\r
544 Quit("Tree::Unroot, not rooted");
\r
546 // Convention: root node is always node zero
\r
548 assert(NULL_NEIGHBOR == m_uNeighbor1[0]);
\r
550 const unsigned uThirdNode = m_uNodeCount++;
\r
552 m_uNeighbor1[0] = uThirdNode;
\r
553 m_uNeighbor1[uThirdNode] = 0;
\r
555 m_uNeighbor2[uThirdNode] = NULL_NEIGHBOR;
\r
556 m_uNeighbor3[uThirdNode] = NULL_NEIGHBOR;
\r
558 m_dEdgeLength1[0] = 0;
\r
559 m_dEdgeLength1[uThirdNode] = 0;
\r
560 m_bHasEdgeLength1[uThirdNode] = true;
\r
565 Log("After unroot:\n");
\r
572 // In an unrooted tree, equivalent of GetLeft/Right is
\r
573 // GetFirst/SecondNeighbor.
\r
574 // uNeighborIndex must be a known neighbor of uNodeIndex.
\r
575 // This is the way to find the other two neighbor nodes of
\r
576 // an internal node.
\r
577 // The labeling as "First" and "Second" neighbor is arbitrary.
\r
578 // Calling these functions on a leaf returns NULL_NEIGHBOR, as
\r
579 // for GetLeft/Right.
\r
580 unsigned Tree::GetFirstNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
\r
582 assert(uNodeIndex < m_uNodeCount);
\r
583 assert(uNeighborIndex < m_uNodeCount);
\r
584 assert(IsEdge(uNodeIndex, uNeighborIndex));
\r
586 for (unsigned n = 0; n < 3; ++n)
\r
588 unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
\r
589 if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
\r
592 return NULL_NEIGHBOR;
\r
595 unsigned Tree::GetSecondNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
\r
597 assert(uNodeIndex < m_uNodeCount);
\r
598 assert(uNeighborIndex < m_uNodeCount);
\r
599 assert(IsEdge(uNodeIndex, uNeighborIndex));
\r
601 bool bFoundOne = false;
\r
602 for (unsigned n = 0; n < 3; ++n)
\r
604 unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
\r
605 if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
\r
613 return NULL_NEIGHBOR;
\r
616 // Compute the number of leaves in the sub-tree defined by an edge
\r
617 // in an unrooted tree. Conceptually, the tree is cut at this edge,
\r
618 // and uNodeIndex2 considered the root of the sub-tree.
\r
619 unsigned Tree::GetLeafCountUnrooted(unsigned uNodeIndex1, unsigned uNodeIndex2,
\r
620 double *ptrdTotalDistance) const
\r
622 assert(!IsRooted());
\r
624 if (IsLeaf(uNodeIndex2))
\r
626 *ptrdTotalDistance = GetEdgeLength(uNodeIndex1, uNodeIndex2);
\r
630 // Recurse down the rooted sub-tree defined by cutting the edge
\r
631 // and considering uNodeIndex2 as the root.
\r
632 const unsigned uLeft = GetFirstNeighbor(uNodeIndex2, uNodeIndex1);
\r
633 const unsigned uRight = GetSecondNeighbor(uNodeIndex2, uNodeIndex1);
\r
635 double dLeftDistance;
\r
636 double dRightDistance;
\r
638 const unsigned uLeftCount = GetLeafCountUnrooted(uNodeIndex2, uLeft,
\r
640 const unsigned uRightCount = GetLeafCountUnrooted(uNodeIndex2, uRight,
\r
643 *ptrdTotalDistance = dLeftDistance + dRightDistance;
\r
644 return uLeftCount + uRightCount;
\r
647 void Tree::RootUnrootedTree(ROOT Method)
\r
649 assert(!IsRooted());
\r
651 Log("Tree::RootUnrootedTree, before:");
\r
659 FindRoot(*this, &uNode1, &uNode2, &dLength1, &dLength2, Method);
\r
661 if (m_uNodeCount == m_uCacheCount)
\r
663 m_uRootNodeIndex = m_uNodeCount++;
\r
665 double dEdgeLength = GetEdgeLength(uNode1, uNode2);
\r
667 m_uNeighbor1[m_uRootNodeIndex] = NULL_NEIGHBOR;
\r
668 m_uNeighbor2[m_uRootNodeIndex] = uNode1;
\r
669 m_uNeighbor3[m_uRootNodeIndex] = uNode2;
\r
671 if (m_uNeighbor1[uNode1] == uNode2)
\r
672 m_uNeighbor1[uNode1] = m_uRootNodeIndex;
\r
673 else if (m_uNeighbor2[uNode1] == uNode2)
\r
674 m_uNeighbor2[uNode1] = m_uRootNodeIndex;
\r
677 assert(m_uNeighbor3[uNode1] == uNode2);
\r
678 m_uNeighbor3[uNode1] = m_uRootNodeIndex;
\r
681 if (m_uNeighbor1[uNode2] == uNode1)
\r
682 m_uNeighbor1[uNode2] = m_uRootNodeIndex;
\r
683 else if (m_uNeighbor2[uNode2] == uNode1)
\r
684 m_uNeighbor2[uNode2] = m_uRootNodeIndex;
\r
687 assert(m_uNeighbor3[uNode2] == uNode1);
\r
688 m_uNeighbor3[uNode2] = m_uRootNodeIndex;
\r
691 OrientParent(uNode1, m_uRootNodeIndex);
\r
692 OrientParent(uNode2, m_uRootNodeIndex);
\r
694 SetEdgeLength(m_uRootNodeIndex, uNode1, dLength1);
\r
695 SetEdgeLength(m_uRootNodeIndex, uNode2, dLength2);
\r
697 m_bHasHeight[m_uRootNodeIndex] = false;
\r
699 m_ptrName[m_uRootNodeIndex] = 0;
\r
704 Log("\nPhy::RootUnrootedTree, after:");
\r
711 bool Tree::HasEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
\r
713 assert(uNodeIndex1 < m_uNodeCount);
\r
714 assert(uNodeIndex2 < m_uNodeCount);
\r
715 assert(IsEdge(uNodeIndex1, uNodeIndex2));
\r
717 if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
\r
718 return m_bHasEdgeLength1[uNodeIndex1];
\r
719 else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
\r
720 return m_bHasEdgeLength2[uNodeIndex1];
\r
721 assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
\r
722 return m_bHasEdgeLength3[uNodeIndex1];
\r
725 void Tree::OrientParent(unsigned uNodeIndex, unsigned uParentNodeIndex)
\r
727 if (NULL_NEIGHBOR == uNodeIndex)
\r
730 if (m_uNeighbor1[uNodeIndex] == uParentNodeIndex)
\r
732 else if (m_uNeighbor2[uNodeIndex] == uParentNodeIndex)
\r
734 double dEdgeLength2 = m_dEdgeLength2[uNodeIndex];
\r
735 m_uNeighbor2[uNodeIndex] = m_uNeighbor1[uNodeIndex];
\r
736 m_dEdgeLength2[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
\r
737 m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
\r
738 m_dEdgeLength1[uNodeIndex] = dEdgeLength2;
\r
742 assert(m_uNeighbor3[uNodeIndex] == uParentNodeIndex);
\r
743 double dEdgeLength3 = m_dEdgeLength3[uNodeIndex];
\r
744 m_uNeighbor3[uNodeIndex] = m_uNeighbor1[uNodeIndex];
\r
745 m_dEdgeLength3[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
\r
746 m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
\r
747 m_dEdgeLength1[uNodeIndex] = dEdgeLength3;
\r
750 OrientParent(m_uNeighbor2[uNodeIndex], uNodeIndex);
\r
751 OrientParent(m_uNeighbor3[uNodeIndex], uNodeIndex);
\r
754 unsigned Tree::FirstDepthFirstNode() const
\r
756 assert(IsRooted());
\r
758 // Descend via left branches until we hit a leaf
\r
759 unsigned uNodeIndex = m_uRootNodeIndex;
\r
760 while (!IsLeaf(uNodeIndex))
\r
761 uNodeIndex = GetLeft(uNodeIndex);
\r
765 unsigned Tree::FirstDepthFirstNodeR() const
\r
767 assert(IsRooted());
\r
769 // Descend via left branches until we hit a leaf
\r
770 unsigned uNodeIndex = m_uRootNodeIndex;
\r
771 while (!IsLeaf(uNodeIndex))
\r
772 uNodeIndex = GetRight(uNodeIndex);
\r
776 unsigned Tree::NextDepthFirstNode(unsigned uNodeIndex) const
\r
779 Log("NextDepthFirstNode(%3u) ", uNodeIndex);
\r
782 assert(IsRooted());
\r
783 assert(uNodeIndex < m_uNodeCount);
\r
785 if (IsRoot(uNodeIndex))
\r
788 Log(">> Node %u is root, end of traversal\n", uNodeIndex);
\r
790 return NULL_NEIGHBOR;
\r
793 unsigned uParent = GetParent(uNodeIndex);
\r
794 if (GetRight(uParent) == uNodeIndex)
\r
797 Log(">> Is right branch, return parent=%u\n", uParent);
\r
802 uNodeIndex = GetRight(uParent);
\r
804 Log(">> Descend left from right sibling=%u ... ", uNodeIndex);
\r
806 while (!IsLeaf(uNodeIndex))
\r
807 uNodeIndex = GetLeft(uNodeIndex);
\r
810 Log("bottom out at leaf=%u\n", uNodeIndex);
\r
815 unsigned Tree::NextDepthFirstNodeR(unsigned uNodeIndex) const
\r
818 Log("NextDepthFirstNode(%3u) ", uNodeIndex);
\r
821 assert(IsRooted());
\r
822 assert(uNodeIndex < m_uNodeCount);
\r
824 if (IsRoot(uNodeIndex))
\r
827 Log(">> Node %u is root, end of traversal\n", uNodeIndex);
\r
829 return NULL_NEIGHBOR;
\r
832 unsigned uParent = GetParent(uNodeIndex);
\r
833 if (GetLeft(uParent) == uNodeIndex)
\r
836 Log(">> Is left branch, return parent=%u\n", uParent);
\r
841 uNodeIndex = GetLeft(uParent);
\r
843 Log(">> Descend right from left sibling=%u ... ", uNodeIndex);
\r
845 while (!IsLeaf(uNodeIndex))
\r
846 uNodeIndex = GetRight(uNodeIndex);
\r
849 Log("bottom out at leaf=%u\n", uNodeIndex);
\r
854 void Tree::UnrootByDeletingRoot()
\r
856 assert(IsRooted());
\r
857 assert(m_uNodeCount >= 3);
\r
859 const unsigned uLeft = GetLeft(m_uRootNodeIndex);
\r
860 const unsigned uRight = GetRight(m_uRootNodeIndex);
\r
862 m_uNeighbor1[uLeft] = uRight;
\r
863 m_uNeighbor1[uRight] = uLeft;
\r
865 bool bHasEdgeLength = HasEdgeLength(m_uRootNodeIndex, uLeft) &&
\r
866 HasEdgeLength(m_uRootNodeIndex, uRight);
\r
867 if (bHasEdgeLength)
\r
869 double dEdgeLength = GetEdgeLength(m_uRootNodeIndex, uLeft) +
\r
870 GetEdgeLength(m_uRootNodeIndex, uRight);
\r
871 m_dEdgeLength1[uLeft] = dEdgeLength;
\r
872 m_dEdgeLength1[uRight] = dEdgeLength;
\r
875 // Remove root node entry from arrays
\r
876 const unsigned uMoveCount = m_uNodeCount - m_uRootNodeIndex;
\r
877 const unsigned uUnsBytes = uMoveCount*sizeof(unsigned);
\r
878 memmove(m_uNeighbor1 + m_uRootNodeIndex, m_uNeighbor1 + m_uRootNodeIndex + 1,
\r
880 memmove(m_uNeighbor2 + m_uRootNodeIndex, m_uNeighbor2 + m_uRootNodeIndex + 1,
\r
882 memmove(m_uNeighbor3 + m_uRootNodeIndex, m_uNeighbor3 + m_uRootNodeIndex + 1,
\r
885 const unsigned uDoubleBytes = uMoveCount*sizeof(double);
\r
886 memmove(m_dEdgeLength1 + m_uRootNodeIndex, m_dEdgeLength1 + m_uRootNodeIndex + 1,
\r
888 memmove(m_dEdgeLength2 + m_uRootNodeIndex, m_dEdgeLength2 + m_uRootNodeIndex + 1,
\r
890 memmove(m_dEdgeLength3 + m_uRootNodeIndex, m_dEdgeLength3 + m_uRootNodeIndex + 1,
\r
893 const unsigned uBoolBytes = uMoveCount*sizeof(bool);
\r
894 memmove(m_bHasEdgeLength1 + m_uRootNodeIndex, m_bHasEdgeLength1 + m_uRootNodeIndex + 1,
\r
896 memmove(m_bHasEdgeLength2 + m_uRootNodeIndex, m_bHasEdgeLength2 + m_uRootNodeIndex + 1,
\r
898 memmove(m_bHasEdgeLength3 + m_uRootNodeIndex, m_bHasEdgeLength3 + m_uRootNodeIndex + 1,
\r
901 const unsigned uPtrBytes = uMoveCount*sizeof(char *);
\r
902 memmove(m_ptrName + m_uRootNodeIndex, m_ptrName + m_uRootNodeIndex + 1, uPtrBytes);
\r
907 // Fix up table entries
\r
908 for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
\r
910 #define DEC(x) if (x != NULL_NEIGHBOR && x > m_uRootNodeIndex) --x;
\r
911 DEC(m_uNeighbor1[uNodeIndex])
\r
912 DEC(m_uNeighbor2[uNodeIndex])
\r
913 DEC(m_uNeighbor3[uNodeIndex])
\r
920 unsigned Tree::GetLeafParent(unsigned uNodeIndex) const
\r
922 assert(IsLeaf(uNodeIndex));
\r
925 return GetParent(uNodeIndex);
\r
927 if (m_uNeighbor1[uNodeIndex] != NULL_NEIGHBOR)
\r
928 return m_uNeighbor1[uNodeIndex];
\r
929 if (m_uNeighbor2[uNodeIndex] != NULL_NEIGHBOR)
\r
930 return m_uNeighbor2[uNodeIndex];
\r
931 return m_uNeighbor3[uNodeIndex];
\r
934 // TODO: This is not efficient for large trees, should cache.
\r
935 double Tree::GetNodeHeight(unsigned uNodeIndex) const
\r
938 Quit("Tree::GetNodeHeight: undefined unless rooted tree");
\r
940 if (IsLeaf(uNodeIndex))
\r
943 if (m_bHasHeight[uNodeIndex])
\r
944 return m_dHeight[uNodeIndex];
\r
946 const unsigned uLeft = GetLeft(uNodeIndex);
\r
947 const unsigned uRight = GetRight(uNodeIndex);
\r
948 double dLeftLength = GetEdgeLength(uNodeIndex, uLeft);
\r
949 double dRightLength = GetEdgeLength(uNodeIndex, uRight);
\r
951 if (dLeftLength < 0)
\r
953 if (dRightLength < 0)
\r
956 const double dLeftHeight = dLeftLength + GetNodeHeight(uLeft);
\r
957 const double dRightHeight = dRightLength + GetNodeHeight(uRight);
\r
958 const double dHeight = (dLeftHeight + dRightHeight)/2;
\r
959 m_bHasHeight[uNodeIndex] = true;
\r
960 m_dHeight[uNodeIndex] = dHeight;
\r
964 unsigned Tree::GetNeighborSubscript(unsigned uNodeIndex, unsigned uNeighborIndex) const
\r
966 assert(uNodeIndex < m_uNodeCount);
\r
967 assert(uNeighborIndex < m_uNodeCount);
\r
968 if (uNeighborIndex == m_uNeighbor1[uNodeIndex])
\r
970 if (uNeighborIndex == m_uNeighbor2[uNodeIndex])
\r
972 if (uNeighborIndex == m_uNeighbor3[uNodeIndex])
\r
974 return NULL_NEIGHBOR;
\r
977 unsigned Tree::GetNeighbor(unsigned uNodeIndex, unsigned uNeighborSubscript) const
\r
979 switch (uNeighborSubscript)
\r
982 return m_uNeighbor1[uNodeIndex];
\r
984 return m_uNeighbor2[uNodeIndex];
\r
986 return m_uNeighbor3[uNodeIndex];
\r
988 Quit("Tree::GetNeighbor, sub=%u", uNeighborSubscript);
\r
989 return NULL_NEIGHBOR;
\r
992 // TODO: check if this is a performance issue, could cache a lookup table
\r
993 unsigned Tree::LeafIndexToNodeIndex(unsigned uLeafIndex) const
\r
995 const unsigned uNodeCount = GetNodeCount();
\r
996 unsigned uLeafCount = 0;
\r
997 for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
\r
999 if (IsLeaf(uNodeIndex))
\r
1001 if (uLeafCount == uLeafIndex)
\r
1002 return uNodeIndex;
\r
1007 Quit("LeafIndexToNodeIndex: out of range");
\r
1011 unsigned Tree::GetLeafNodeIndex(const char *ptrName) const
\r
1013 const unsigned uNodeCount = GetNodeCount();
\r
1014 for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
\r
1016 if (!IsLeaf(uNodeIndex))
\r
1018 const char *ptrLeafName = GetLeafName(uNodeIndex);
\r
1019 if (0 == strcmp(ptrName, ptrLeafName))
\r
1020 return uNodeIndex;
\r
1022 Quit("Tree::GetLeafNodeIndex, name not found");
\r
1026 void Tree::Copy(const Tree &tree)
\r
1028 const unsigned uNodeCount = tree.GetNodeCount();
\r
1029 InitCache(uNodeCount);
\r
1031 m_uNodeCount = uNodeCount;
\r
1033 const size_t UnsignedBytes = uNodeCount*sizeof(unsigned);
\r
1034 const size_t DoubleBytes = uNodeCount*sizeof(double);
\r
1035 const size_t BoolBytes = uNodeCount*sizeof(bool);
\r
1037 memcpy(m_uNeighbor1, tree.m_uNeighbor1, UnsignedBytes);
\r
1038 memcpy(m_uNeighbor2, tree.m_uNeighbor2, UnsignedBytes);
\r
1039 memcpy(m_uNeighbor3, tree.m_uNeighbor3, UnsignedBytes);
\r
1041 memcpy(m_Ids, tree.m_Ids, UnsignedBytes);
\r
1043 memcpy(m_dEdgeLength1, tree.m_dEdgeLength1, DoubleBytes);
\r
1044 memcpy(m_dEdgeLength2, tree.m_dEdgeLength2, DoubleBytes);
\r
1045 memcpy(m_dEdgeLength3, tree.m_dEdgeLength3, DoubleBytes);
\r
1047 memcpy(m_dHeight, tree.m_dHeight, DoubleBytes);
\r
1049 memcpy(m_bHasEdgeLength1, tree.m_bHasEdgeLength1, BoolBytes);
\r
1050 memcpy(m_bHasEdgeLength2, tree.m_bHasEdgeLength2, BoolBytes);
\r
1051 memcpy(m_bHasEdgeLength3, tree.m_bHasEdgeLength3, BoolBytes);
\r
1053 memcpy(m_bHasHeight, tree.m_bHasHeight, BoolBytes);
\r
1055 m_uRootNodeIndex = tree.m_uRootNodeIndex;
\r
1056 m_bRooted = tree.m_bRooted;
\r
1058 for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
\r
1060 if (tree.IsLeaf(uNodeIndex))
\r
1062 const char *ptrName = tree.GetLeafName(uNodeIndex);
\r
1063 m_ptrName[uNodeIndex] = strsave(ptrName);
\r
1066 m_ptrName[uNodeIndex] = 0;
\r
1074 // Create rooted tree from a vector description.
\r
1075 // Node indexes are 0..N-1 for leaves, N..2N-2 for
\r
1076 // internal nodes.
\r
1077 // Vector subscripts are i-N and have values for
\r
1078 // internal nodes only, but those values are node
\r
1079 // indexes 0..2N-2. So e.g. if N=6 and Left[2]=1,
\r
1080 // this means that the third internal node (node index 8)
\r
1081 // has the second leaf (node index 1) as its left child.
\r
1082 // uRoot gives the vector subscript of the root, so add N
\r
1083 // to get the node index.
\r
1084 void Tree::Create(unsigned uLeafCount, unsigned uRoot, const unsigned Left[],
\r
1085 const unsigned Right[], const float LeftLength[], const float RightLength[],
\r
1086 const unsigned LeafIds[], char **LeafNames)
\r
1090 m_uNodeCount = 2*uLeafCount - 1;
\r
1091 InitCache(m_uNodeCount);
\r
1093 for (unsigned uNodeIndex = 0; uNodeIndex < uLeafCount; ++uNodeIndex)
\r
1095 m_Ids[uNodeIndex] = LeafIds[uNodeIndex];
\r
1096 m_ptrName[uNodeIndex] = strsave(LeafNames[uNodeIndex]);
\r
1099 for (unsigned uNodeIndex = uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)
\r
1101 unsigned v = uNodeIndex - uLeafCount;
\r
1102 unsigned uLeft = Left[v];
\r
1103 unsigned uRight = Right[v];
\r
1104 float fLeft = LeftLength[v];
\r
1105 float fRight = RightLength[v];
\r
1107 m_uNeighbor2[uNodeIndex] = uLeft;
\r
1108 m_uNeighbor3[uNodeIndex] = uRight;
\r
1110 m_bHasEdgeLength2[uNodeIndex] = true;
\r
1111 m_bHasEdgeLength3[uNodeIndex] = true;
\r
1113 m_dEdgeLength2[uNodeIndex] = fLeft;
\r
1114 m_dEdgeLength3[uNodeIndex] = fRight;
\r
1116 m_uNeighbor1[uLeft] = uNodeIndex;
\r
1117 m_uNeighbor1[uRight] = uNodeIndex;
\r
1119 m_dEdgeLength1[uLeft] = fLeft;
\r
1120 m_dEdgeLength1[uRight] = fRight;
\r
1122 m_bHasEdgeLength1[uLeft] = true;
\r
1123 m_bHasEdgeLength1[uRight] = true;
\r
1127 m_uRootNodeIndex = uRoot + uLeafCount;
\r