source: trunk/GDE/MUSCLE/src/phy.cpp

Last change on this file was 10390, checked in by aboeckma, 11 years ago

added muscle sourcles amd makefile
File size: 30.2 KB
Line 
1#include "muscle.h"
2#include "tree.h"
3#include <math.h>
4
5#define TRACE 0
6
7/***
8Node has 0 to 3 neighbors:
9        0 neighbors:    singleton root
10        1 neighbor:             leaf, neighbor is parent
11        2 neigbors:             non-singleton root
12        3 neighbors:    internal node (other than root)
13
14Minimal rooted tree is single node.
15Minimal unrooted tree is single edge.
16Leaf node always has nulls in neighbors 2 and 3, neighbor 1 is parent.
17When tree is rooted, neighbor 1=parent, 2=left, 3=right.
18***/
19
20void Tree::AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2) const
21        {
22        if (uNodeIndex1 >= m_uNodeCount || uNodeIndex2 >= m_uNodeCount)
23                Quit("AssertAreNeighbors(%u,%u), are %u nodes",
24                  uNodeIndex1, uNodeIndex2, m_uNodeCount);
25
26        if (m_uNeighbor1[uNodeIndex1] != uNodeIndex2 &&
27          m_uNeighbor2[uNodeIndex1] != uNodeIndex2 &&
28          m_uNeighbor3[uNodeIndex1] != uNodeIndex2)
29                {
30                LogMe();
31                Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
32                }
33
34        if (m_uNeighbor1[uNodeIndex2] != uNodeIndex1 &&
35          m_uNeighbor2[uNodeIndex2] != uNodeIndex1 &&
36          m_uNeighbor3[uNodeIndex2] != uNodeIndex1)
37                {
38                LogMe();
39                Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);
40                }
41
42        bool Has12 = HasEdgeLength(uNodeIndex1, uNodeIndex2);
43        bool Has21 = HasEdgeLength(uNodeIndex2, uNodeIndex1);
44        if (Has12 != Has21)
45                {
46                HasEdgeLength(uNodeIndex1, uNodeIndex2);
47                HasEdgeLength(uNodeIndex2, uNodeIndex1);
48                LogMe();
49                Log("HasEdgeLength(%u, %u)=%c HasEdgeLength(%u, %u)=%c\n",
50                  uNodeIndex1,
51                  uNodeIndex2,
52                  Has12 ? 'T' : 'F',
53                  uNodeIndex2,
54                  uNodeIndex1,
55                  Has21 ? 'T' : 'F');
56
57                Quit("Tree::AssertAreNeighbors, HasEdgeLength not symmetric");
58                }
59
60        if (Has12)
61                {
62                double d12 = GetEdgeLength(uNodeIndex1, uNodeIndex2);
63                double d21 = GetEdgeLength(uNodeIndex2, uNodeIndex1);
64                if (d12 != d21)
65                        {
66                        LogMe();
67                        Quit("Tree::AssertAreNeighbors, Edge length disagrees %u-%u=%.3g, %u-%u=%.3g",
68                          uNodeIndex1, uNodeIndex2, d12,
69                          uNodeIndex2, uNodeIndex1, d21);
70                        }
71                }
72        }
73
74void Tree::ValidateNode(unsigned uNodeIndex) const
75        {
76        if (uNodeIndex >= m_uNodeCount)
77                Quit("ValidateNode(%u), %u nodes", uNodeIndex, m_uNodeCount);
78
79        const unsigned uNeighborCount = GetNeighborCount(uNodeIndex);
80
81        if (2 == uNeighborCount)
82                {
83                if (!m_bRooted)
84                        {
85                        LogMe();
86                        Quit("Tree::ValidateNode: Node %u has two neighbors, tree is not rooted",
87                         uNodeIndex);
88                        }
89                if (uNodeIndex != m_uRootNodeIndex)
90                        {
91                        LogMe();
92                        Quit("Tree::ValidateNode: Node %u has two neighbors, but not root node=%u",
93                         uNodeIndex, m_uRootNodeIndex);
94                        }
95                }
96
97        const unsigned n1 = m_uNeighbor1[uNodeIndex];
98        const unsigned n2 = m_uNeighbor2[uNodeIndex];
99        const unsigned n3 = m_uNeighbor3[uNodeIndex];
100
101        if (NULL_NEIGHBOR == n2 && NULL_NEIGHBOR != n3)
102                {
103                LogMe();
104                Quit("Tree::ValidateNode, n2=null, n3!=null", uNodeIndex);
105                }
106        if (NULL_NEIGHBOR == n3 && NULL_NEIGHBOR != n2)
107                {
108                LogMe();
109                Quit("Tree::ValidateNode, n3=null, n2!=null", uNodeIndex);
110                }
111
112        if (n1 != NULL_NEIGHBOR)
113                AssertAreNeighbors(uNodeIndex, n1);
114        if (n2 != NULL_NEIGHBOR)
115                AssertAreNeighbors(uNodeIndex, n2);
116        if (n3 != NULL_NEIGHBOR)
117                AssertAreNeighbors(uNodeIndex, n3);
118
119        if (n1 != NULL_NEIGHBOR && (n1 == n2 || n1 == n3))
120                {
121                LogMe();
122                Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
123                }
124        if (n2 != NULL_NEIGHBOR && (n2 == n1 || n2 == n3))
125                {
126                LogMe();
127                Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
128                }
129        if (n3 != NULL_NEIGHBOR && (n3 == n1 || n3 == n2))
130                {
131                LogMe();
132                Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);
133                }
134
135        if (IsRooted())
136                {
137                if (NULL_NEIGHBOR == GetParent(uNodeIndex))
138                        {
139                        if (uNodeIndex != m_uRootNodeIndex)
140                                {
141                                LogMe();
142                                Quit("Tree::ValiateNode(%u), no parent", uNodeIndex);
143                                }
144                        }
145                else if (GetLeft(GetParent(uNodeIndex)) != uNodeIndex &&
146                  GetRight(GetParent(uNodeIndex)) != uNodeIndex)
147                        {
148                        LogMe();
149                        Quit("Tree::ValidateNode(%u), parent / child mismatch", uNodeIndex);
150                        }
151                }
152        }
153
154void Tree::Validate() const
155        {
156        for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
157                ValidateNode(uNodeIndex);
158        }
159
160bool Tree::IsEdge(unsigned uNodeIndex1, unsigned uNodeIndex2) const
161        {
162        assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
163
164        return m_uNeighbor1[uNodeIndex1] == uNodeIndex2 ||
165          m_uNeighbor2[uNodeIndex1] == uNodeIndex2 ||
166          m_uNeighbor3[uNodeIndex1] == uNodeIndex2;
167        }
168
169double Tree::GetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
170        {
171        assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
172        if (!HasEdgeLength(uNodeIndex1, uNodeIndex2))
173                {
174                LogMe();
175                Quit("Missing edge length in tree %u-%u", uNodeIndex1, uNodeIndex2);
176                }
177
178        if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
179                return m_dEdgeLength1[uNodeIndex1];
180        else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
181                return m_dEdgeLength2[uNodeIndex1];
182        assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
183        return m_dEdgeLength3[uNodeIndex1];
184        }
185
186void Tree::ExpandCache()
187        {
188        const unsigned uNodeCount = 100;
189        unsigned uNewCacheCount = m_uCacheCount + uNodeCount;
190        unsigned *uNewNeighbor1 = new unsigned[uNewCacheCount];
191        unsigned *uNewNeighbor2 = new unsigned[uNewCacheCount];
192        unsigned *uNewNeighbor3 = new unsigned[uNewCacheCount];
193
194        unsigned *uNewIds = new unsigned[uNewCacheCount];
195        memset(uNewIds, 0xff, uNewCacheCount*sizeof(unsigned));
196
197        double *dNewEdgeLength1 = new double[uNewCacheCount];
198        double *dNewEdgeLength2 = new double[uNewCacheCount];
199        double *dNewEdgeLength3 = new double[uNewCacheCount];
200        double *dNewHeight = new double[uNewCacheCount];
201
202        bool *bNewHasEdgeLength1 = new bool[uNewCacheCount];
203        bool *bNewHasEdgeLength2 = new bool[uNewCacheCount];
204        bool *bNewHasEdgeLength3 = new bool[uNewCacheCount];
205        bool *bNewHasHeight = new bool[uNewCacheCount];
206
207        char **ptrNewName = new char *[uNewCacheCount];
208        memset(ptrNewName, 0, uNewCacheCount*sizeof(char *));
209
210        if (m_uCacheCount > 0)
211                {
212                const unsigned uUnsignedBytes = m_uCacheCount*sizeof(unsigned);
213                memcpy(uNewNeighbor1, m_uNeighbor1, uUnsignedBytes);
214                memcpy(uNewNeighbor2, m_uNeighbor2, uUnsignedBytes);
215                memcpy(uNewNeighbor3, m_uNeighbor3, uUnsignedBytes);
216
217                memcpy(uNewIds, m_Ids, uUnsignedBytes);
218
219                const unsigned uEdgeBytes = m_uCacheCount*sizeof(double);
220                memcpy(dNewEdgeLength1, m_dEdgeLength1, uEdgeBytes);
221                memcpy(dNewEdgeLength2, m_dEdgeLength2, uEdgeBytes);
222                memcpy(dNewEdgeLength3, m_dEdgeLength3, uEdgeBytes);
223                memcpy(dNewHeight, m_dHeight, uEdgeBytes);
224
225                const unsigned uBoolBytes = m_uCacheCount*sizeof(bool);
226                memcpy(bNewHasEdgeLength1, m_bHasEdgeLength1, uBoolBytes);
227                memcpy(bNewHasEdgeLength2, m_bHasEdgeLength2, uBoolBytes);
228                memcpy(bNewHasEdgeLength3, m_bHasEdgeLength3, uBoolBytes);
229                memcpy(bNewHasHeight, m_bHasHeight, uBoolBytes);
230
231                const unsigned uNameBytes = m_uCacheCount*sizeof(char *);
232                memcpy(ptrNewName, m_ptrName, uNameBytes);
233
234                delete[] m_uNeighbor1;
235                delete[] m_uNeighbor2;
236                delete[] m_uNeighbor3;
237
238                delete[] m_Ids;
239
240                delete[] m_dEdgeLength1;
241                delete[] m_dEdgeLength2;
242                delete[] m_dEdgeLength3;
243
244                delete[] m_bHasEdgeLength1;
245                delete[] m_bHasEdgeLength2;
246                delete[] m_bHasEdgeLength3;
247                delete[] m_bHasHeight;
248
249                delete[] m_ptrName;
250                }
251        m_uCacheCount = uNewCacheCount;
252        m_uNeighbor1 = uNewNeighbor1;
253        m_uNeighbor2 = uNewNeighbor2;
254        m_uNeighbor3 = uNewNeighbor3;
255        m_Ids = uNewIds;
256        m_dEdgeLength1 = dNewEdgeLength1;
257        m_dEdgeLength2 = dNewEdgeLength2;
258        m_dEdgeLength3 = dNewEdgeLength3;
259        m_dHeight = dNewHeight;
260        m_bHasEdgeLength1 = bNewHasEdgeLength1;
261        m_bHasEdgeLength2 = bNewHasEdgeLength2;
262        m_bHasEdgeLength3 = bNewHasEdgeLength3;
263        m_bHasHeight = bNewHasHeight;
264        m_ptrName = ptrNewName;
265        }
266
267// Creates tree with single node, no edges.
268// Root node always has index 0.
269void Tree::CreateRooted()
270        {
271        Clear();
272        ExpandCache();
273        m_uNodeCount = 1;
274
275        m_uNeighbor1[0] = NULL_NEIGHBOR;
276        m_uNeighbor2[0] = NULL_NEIGHBOR;
277        m_uNeighbor3[0] = NULL_NEIGHBOR;
278
279        m_bHasEdgeLength1[0] = false;
280        m_bHasEdgeLength2[0] = false;
281        m_bHasEdgeLength3[0] = false;
282        m_bHasHeight[0] = false;
283
284        m_uRootNodeIndex = 0;
285        m_bRooted = true;
286
287#if     DEBUG
288        Validate();
289#endif
290        }
291
292// Creates unrooted tree with single edge.
293// Nodes for that edge are always 0 and 1.
294void Tree::CreateUnrooted(double dEdgeLength)
295        {
296        Clear();
297        ExpandCache();
298
299        m_uNeighbor1[0] = 1;
300        m_uNeighbor2[0] = NULL_NEIGHBOR;
301        m_uNeighbor3[0] = NULL_NEIGHBOR;
302
303        m_uNeighbor1[1] = 0;
304        m_uNeighbor2[1] = NULL_NEIGHBOR;
305        m_uNeighbor3[1] = NULL_NEIGHBOR;
306
307        m_dEdgeLength1[0] = dEdgeLength;
308        m_dEdgeLength1[1] = dEdgeLength;
309
310        m_bHasEdgeLength1[0] = true;
311        m_bHasEdgeLength1[1] = true;
312
313        m_bRooted = false;
314
315#if     DEBUG
316        Validate();
317#endif
318        }
319
320void Tree::SetLeafName(unsigned uNodeIndex, const char *ptrName)
321        {
322        assert(uNodeIndex < m_uNodeCount);
323        assert(IsLeaf(uNodeIndex));
324        free(m_ptrName[uNodeIndex]);
325        m_ptrName[uNodeIndex] = strsave(ptrName);
326        }
327
328void Tree::SetLeafId(unsigned uNodeIndex, unsigned uId)
329        {
330        assert(uNodeIndex < m_uNodeCount);
331        assert(IsLeaf(uNodeIndex));
332        m_Ids[uNodeIndex] = uId;
333        }
334
335const char *Tree::GetLeafName(unsigned uNodeIndex) const
336        {
337        assert(uNodeIndex < m_uNodeCount);
338        assert(IsLeaf(uNodeIndex));
339        return m_ptrName[uNodeIndex];
340        }
341
342unsigned Tree::GetLeafId(unsigned uNodeIndex) const
343        {
344        assert(uNodeIndex < m_uNodeCount);
345        assert(IsLeaf(uNodeIndex));
346        return m_Ids[uNodeIndex];
347        }
348
349// Append a new branch.
350// This adds two new nodes and joins them to an existing leaf node.
351// Return value is k, new nodes have indexes k and k+1 respectively.
352unsigned Tree::AppendBranch(unsigned uExistingLeafIndex)
353        {
354        if (0 == m_uNodeCount)
355                Quit("Tree::AppendBranch: tree has not been created");
356
357#if     DEBUG
358        assert(uExistingLeafIndex < m_uNodeCount);
359        if (!IsLeaf(uExistingLeafIndex))
360                {
361                LogMe();
362                Quit("AppendBranch(%u): not leaf", uExistingLeafIndex);
363                }
364#endif
365
366        if (m_uNodeCount >= m_uCacheCount - 2)
367                ExpandCache();
368
369        const unsigned uNewLeaf1 = m_uNodeCount;
370        const unsigned uNewLeaf2 = m_uNodeCount + 1;
371
372        m_uNodeCount += 2;
373
374        assert(m_uNeighbor2[uExistingLeafIndex] == NULL_NEIGHBOR);
375        assert(m_uNeighbor3[uExistingLeafIndex] == NULL_NEIGHBOR);
376
377        m_uNeighbor2[uExistingLeafIndex] = uNewLeaf1;
378        m_uNeighbor3[uExistingLeafIndex] = uNewLeaf2;
379
380        m_uNeighbor1[uNewLeaf1] = uExistingLeafIndex;
381        m_uNeighbor1[uNewLeaf2] = uExistingLeafIndex;
382
383        m_uNeighbor2[uNewLeaf1] = NULL_NEIGHBOR;
384        m_uNeighbor2[uNewLeaf2] = NULL_NEIGHBOR;
385
386        m_uNeighbor3[uNewLeaf1] = NULL_NEIGHBOR;
387        m_uNeighbor3[uNewLeaf2] = NULL_NEIGHBOR;
388
389        m_dEdgeLength2[uExistingLeafIndex] = 0;
390        m_dEdgeLength3[uExistingLeafIndex] = 0;
391
392        m_dEdgeLength1[uNewLeaf1] = 0;
393        m_dEdgeLength2[uNewLeaf1] = 0;
394        m_dEdgeLength3[uNewLeaf1] = 0;
395
396        m_dEdgeLength1[uNewLeaf2] = 0;
397        m_dEdgeLength2[uNewLeaf2] = 0;
398        m_dEdgeLength3[uNewLeaf2] = 0;
399
400        m_bHasEdgeLength1[uNewLeaf1] = false;
401        m_bHasEdgeLength2[uNewLeaf1] = false;
402        m_bHasEdgeLength3[uNewLeaf1] = false;
403
404        m_bHasEdgeLength1[uNewLeaf2] = false;
405        m_bHasEdgeLength2[uNewLeaf2] = false;
406        m_bHasEdgeLength3[uNewLeaf2] = false;
407
408        m_bHasHeight[uNewLeaf1] = false;
409        m_bHasHeight[uNewLeaf2] = false;
410
411        m_Ids[uNewLeaf1] = uInsane;
412        m_Ids[uNewLeaf2] = uInsane;
413        return uNewLeaf1;
414        }
415
416void Tree::LogMe() const
417        {
418        Log("Tree::LogMe %u nodes, ", m_uNodeCount);
419
420        if (IsRooted())
421                {
422                Log("rooted.\n");
423                Log("\n");
424                Log("Index  Parnt  LengthP  Left   LengthL  Right  LengthR     Id  Name\n");
425                Log("-----  -----  -------  ----   -------  -----  -------  -----  ----\n");
426                }
427        else
428                {
429                Log("unrooted.\n");
430                Log("\n");
431                Log("Index  Nbr_1  Length1  Nbr_2  Length2  Nbr_3  Length3     Id  Name\n");
432                Log("-----  -----  -------  -----  -------  -----  -------  -----  ----\n");
433                }
434
435        for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
436                {
437                Log("%5u  ", uNodeIndex);
438                const unsigned n1 = m_uNeighbor1[uNodeIndex];
439                const unsigned n2 = m_uNeighbor2[uNodeIndex];
440                const unsigned n3 = m_uNeighbor3[uNodeIndex];
441                if (NULL_NEIGHBOR != n1)
442                        {
443                        Log("%5u  ", n1);
444                        if (m_bHasEdgeLength1[uNodeIndex])
445                                Log("%7.4f  ", m_dEdgeLength1[uNodeIndex]);
446                        else
447                                Log("      *  ");
448                        }
449                else
450                        Log("                ");
451
452                if (NULL_NEIGHBOR != n2)
453                        {
454                        Log("%5u  ", n2);
455                        if (m_bHasEdgeLength2[uNodeIndex])
456                                Log("%7.4f  ", m_dEdgeLength2[uNodeIndex]);
457                        else
458                                Log("      *  ");
459                        }
460                else
461                        Log("                ");
462
463                if (NULL_NEIGHBOR != n3)
464                        {
465                        Log("%5u  ", n3);
466                        if (m_bHasEdgeLength3[uNodeIndex])
467                                Log("%7.4f  ", m_dEdgeLength3[uNodeIndex]);
468                        else
469                                Log("      *  ");
470                        }
471                else
472                        Log("                ");
473
474                if (m_Ids != 0 && IsLeaf(uNodeIndex))
475                        {
476                        unsigned uId = m_Ids[uNodeIndex];
477                        if (uId == uInsane)
478                                Log("    *");
479                        else
480                                Log("%5u", uId);
481                        }
482                else
483                        Log("     ");
484
485                if (m_bRooted && uNodeIndex == m_uRootNodeIndex)
486                        Log("  [ROOT] ");
487                const char *ptrName = m_ptrName[uNodeIndex];
488                if (ptrName != 0)
489                        Log("  %s", ptrName);
490                Log("\n");
491                }
492        }
493
494void Tree::SetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2,
495  double dLength)
496        {
497        assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);
498        assert(IsEdge(uNodeIndex1, uNodeIndex2));
499
500        if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
501                {
502                m_dEdgeLength1[uNodeIndex1] = dLength;
503                m_bHasEdgeLength1[uNodeIndex1] = true;
504                }
505        else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
506                {
507                m_dEdgeLength2[uNodeIndex1] = dLength;
508                m_bHasEdgeLength2[uNodeIndex1] = true;
509
510                }
511        else
512                {
513                assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
514                m_dEdgeLength3[uNodeIndex1] = dLength;
515                m_bHasEdgeLength3[uNodeIndex1] = true;
516                }
517
518        if (m_uNeighbor1[uNodeIndex2] == uNodeIndex1)
519                {
520                m_dEdgeLength1[uNodeIndex2] = dLength;
521                m_bHasEdgeLength1[uNodeIndex2] = true;
522                }
523        else if (m_uNeighbor2[uNodeIndex2] == uNodeIndex1)
524                {
525                m_dEdgeLength2[uNodeIndex2] = dLength;
526                m_bHasEdgeLength2[uNodeIndex2] = true;
527                }
528        else
529                {
530                assert(m_uNeighbor3[uNodeIndex2] == uNodeIndex1);
531                m_dEdgeLength3[uNodeIndex2] = dLength;
532                m_bHasEdgeLength3[uNodeIndex2] = true;
533                }
534        }
535
536unsigned Tree::UnrootFromFile()
537        {
538#if     TRACE
539        Log("Before unroot:\n");
540        LogMe();
541#endif
542
543        if (!m_bRooted)
544                Quit("Tree::Unroot, not rooted");
545
546// Convention: root node is always node zero
547        assert(IsRoot(0));
548        assert(NULL_NEIGHBOR == m_uNeighbor1[0]);
549
550        const unsigned uThirdNode = m_uNodeCount++;
551
552        m_uNeighbor1[0] = uThirdNode;
553        m_uNeighbor1[uThirdNode] = 0;
554
555        m_uNeighbor2[uThirdNode] = NULL_NEIGHBOR;
556        m_uNeighbor3[uThirdNode] = NULL_NEIGHBOR;
557
558        m_dEdgeLength1[0] = 0;
559        m_dEdgeLength1[uThirdNode] = 0;
560        m_bHasEdgeLength1[uThirdNode] = true;
561
562        m_bRooted = false;
563
564#if     TRACE
565        Log("After unroot:\n");
566        LogMe();
567#endif
568
569        return uThirdNode;
570        }
571
572// In an unrooted tree, equivalent of GetLeft/Right is
573// GetFirst/SecondNeighbor.
574// uNeighborIndex must be a known neighbor of uNodeIndex.
575// This is the way to find the other two neighbor nodes of
576// an internal node.
577// The labeling as "First" and "Second" neighbor is arbitrary.
578// Calling these functions on a leaf returns NULL_NEIGHBOR, as
579// for GetLeft/Right.
580unsigned Tree::GetFirstNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
581        {
582        assert(uNodeIndex < m_uNodeCount);
583        assert(uNeighborIndex < m_uNodeCount);
584        assert(IsEdge(uNodeIndex, uNeighborIndex));
585
586        for (unsigned n = 0; n < 3; ++n)
587                {
588                unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
589                if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
590                        return uNeighbor;
591                }
592        return NULL_NEIGHBOR;
593        }
594
595unsigned Tree::GetSecondNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const
596        {
597        assert(uNodeIndex < m_uNodeCount);
598        assert(uNeighborIndex < m_uNodeCount);
599        assert(IsEdge(uNodeIndex, uNeighborIndex));
600
601        bool bFoundOne = false;
602        for (unsigned n = 0; n < 3; ++n)
603                {
604                unsigned uNeighbor = GetNeighbor(uNodeIndex, n);
605                if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)
606                        {
607                        if (bFoundOne)
608                                return uNeighbor;
609                        else
610                                bFoundOne = true;
611                        }
612                }
613        return NULL_NEIGHBOR;
614        }
615
616// Compute the number of leaves in the sub-tree defined by an edge
617// in an unrooted tree. Conceptually, the tree is cut at this edge,
618// and uNodeIndex2 considered the root of the sub-tree.
619unsigned Tree::GetLeafCountUnrooted(unsigned uNodeIndex1, unsigned uNodeIndex2,
620  double *ptrdTotalDistance) const
621        {
622        assert(!IsRooted());
623
624        if (IsLeaf(uNodeIndex2))
625                {
626                *ptrdTotalDistance = GetEdgeLength(uNodeIndex1, uNodeIndex2);
627                return 1;
628                }
629
630// Recurse down the rooted sub-tree defined by cutting the edge
631// and considering uNodeIndex2 as the root.
632        const unsigned uLeft = GetFirstNeighbor(uNodeIndex2, uNodeIndex1);
633        const unsigned uRight = GetSecondNeighbor(uNodeIndex2, uNodeIndex1);
634
635        double dLeftDistance;
636        double dRightDistance;
637
638        const unsigned uLeftCount = GetLeafCountUnrooted(uNodeIndex2, uLeft,
639          &dLeftDistance);
640        const unsigned uRightCount = GetLeafCountUnrooted(uNodeIndex2, uRight,
641          &dRightDistance);
642
643        *ptrdTotalDistance = dLeftDistance + dRightDistance;
644        return uLeftCount + uRightCount;
645        }
646
647void Tree::RootUnrootedTree(ROOT Method)
648        {
649        assert(!IsRooted());
650#if TRACE
651        Log("Tree::RootUnrootedTree, before:");
652        LogMe();
653#endif
654
655        unsigned uNode1;
656        unsigned uNode2;
657        double dLength1;
658        double dLength2;
659        FindRoot(*this, &uNode1, &uNode2, &dLength1, &dLength2, Method);
660
661        if (m_uNodeCount == m_uCacheCount)
662                ExpandCache();
663        m_uRootNodeIndex = m_uNodeCount++;
664
665        double dEdgeLength = GetEdgeLength(uNode1, uNode2);
666
667        m_uNeighbor1[m_uRootNodeIndex] = NULL_NEIGHBOR;
668        m_uNeighbor2[m_uRootNodeIndex] = uNode1;
669        m_uNeighbor3[m_uRootNodeIndex] = uNode2;
670
671        if (m_uNeighbor1[uNode1] == uNode2)
672                m_uNeighbor1[uNode1] = m_uRootNodeIndex;
673        else if (m_uNeighbor2[uNode1] == uNode2)
674                m_uNeighbor2[uNode1] = m_uRootNodeIndex;
675        else
676                {
677                assert(m_uNeighbor3[uNode1] == uNode2);
678                m_uNeighbor3[uNode1] = m_uRootNodeIndex;
679                }
680
681        if (m_uNeighbor1[uNode2] == uNode1)
682                m_uNeighbor1[uNode2] = m_uRootNodeIndex;
683        else if (m_uNeighbor2[uNode2] == uNode1)
684                m_uNeighbor2[uNode2] = m_uRootNodeIndex;
685        else
686                {
687                assert(m_uNeighbor3[uNode2] == uNode1);
688                m_uNeighbor3[uNode2] = m_uRootNodeIndex;
689                }
690
691        OrientParent(uNode1, m_uRootNodeIndex);
692        OrientParent(uNode2, m_uRootNodeIndex);
693
694        SetEdgeLength(m_uRootNodeIndex, uNode1, dLength1);
695        SetEdgeLength(m_uRootNodeIndex, uNode2, dLength2);
696
697        m_bHasHeight[m_uRootNodeIndex] = false;
698
699        m_ptrName[m_uRootNodeIndex] = 0;
700
701        m_bRooted = true;
702
703#if     TRACE
704        Log("\nPhy::RootUnrootedTree, after:");
705        LogMe();
706#endif
707
708        Validate();
709        }
710
711bool Tree::HasEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const
712        {
713        assert(uNodeIndex1 < m_uNodeCount);
714        assert(uNodeIndex2 < m_uNodeCount);
715        assert(IsEdge(uNodeIndex1, uNodeIndex2));
716
717        if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)
718                return m_bHasEdgeLength1[uNodeIndex1];
719        else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)
720                return m_bHasEdgeLength2[uNodeIndex1];
721        assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);
722        return m_bHasEdgeLength3[uNodeIndex1];
723        }
724
725void Tree::OrientParent(unsigned uNodeIndex, unsigned uParentNodeIndex)
726        {
727        if (NULL_NEIGHBOR == uNodeIndex)
728                return;
729
730        if (m_uNeighbor1[uNodeIndex] == uParentNodeIndex)
731                ;
732        else if (m_uNeighbor2[uNodeIndex] == uParentNodeIndex)
733                {
734                double dEdgeLength2 = m_dEdgeLength2[uNodeIndex];
735                m_uNeighbor2[uNodeIndex] = m_uNeighbor1[uNodeIndex];
736                m_dEdgeLength2[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
737                m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
738                m_dEdgeLength1[uNodeIndex] = dEdgeLength2;
739                }
740        else
741                {
742                assert(m_uNeighbor3[uNodeIndex] == uParentNodeIndex);
743                double dEdgeLength3 = m_dEdgeLength3[uNodeIndex];
744                m_uNeighbor3[uNodeIndex] = m_uNeighbor1[uNodeIndex];
745                m_dEdgeLength3[uNodeIndex] = m_dEdgeLength1[uNodeIndex];
746                m_uNeighbor1[uNodeIndex] = uParentNodeIndex;
747                m_dEdgeLength1[uNodeIndex] = dEdgeLength3;
748                }
749
750        OrientParent(m_uNeighbor2[uNodeIndex], uNodeIndex);
751        OrientParent(m_uNeighbor3[uNodeIndex], uNodeIndex);
752        }
753
754unsigned Tree::FirstDepthFirstNode() const
755        {
756        assert(IsRooted());
757
758// Descend via left branches until we hit a leaf
759        unsigned uNodeIndex = m_uRootNodeIndex;
760        while (!IsLeaf(uNodeIndex))
761                uNodeIndex = GetLeft(uNodeIndex);
762        return uNodeIndex;
763        }
764
765unsigned Tree::FirstDepthFirstNodeR() const
766        {
767        assert(IsRooted());
768
769// Descend via left branches until we hit a leaf
770        unsigned uNodeIndex = m_uRootNodeIndex;
771        while (!IsLeaf(uNodeIndex))
772                uNodeIndex = GetRight(uNodeIndex);
773        return uNodeIndex;
774        }
775
776unsigned Tree::NextDepthFirstNode(unsigned uNodeIndex) const
777        {
778#if     TRACE
779        Log("NextDepthFirstNode(%3u) ", uNodeIndex);
780#endif
781
782        assert(IsRooted());
783        assert(uNodeIndex < m_uNodeCount);
784
785        if (IsRoot(uNodeIndex))
786                {
787#if     TRACE
788                Log(">> Node %u is root, end of traversal\n", uNodeIndex);
789#endif
790                return NULL_NEIGHBOR;
791                }
792
793        unsigned uParent = GetParent(uNodeIndex);
794        if (GetRight(uParent) == uNodeIndex)
795                {
796#if     TRACE
797                Log(">> Is right branch, return parent=%u\n", uParent);
798#endif
799                return uParent;
800                }
801
802        uNodeIndex = GetRight(uParent);
803#if     TRACE
804                Log(">> Descend left from right sibling=%u ... ", uNodeIndex);
805#endif
806        while (!IsLeaf(uNodeIndex))
807                uNodeIndex = GetLeft(uNodeIndex);
808
809#if     TRACE
810        Log("bottom out at leaf=%u\n", uNodeIndex);
811#endif
812        return uNodeIndex;
813        }
814
815unsigned Tree::NextDepthFirstNodeR(unsigned uNodeIndex) const
816        {
817#if     TRACE
818        Log("NextDepthFirstNode(%3u) ", uNodeIndex);
819#endif
820
821        assert(IsRooted());
822        assert(uNodeIndex < m_uNodeCount);
823
824        if (IsRoot(uNodeIndex))
825                {
826#if     TRACE
827                Log(">> Node %u is root, end of traversal\n", uNodeIndex);
828#endif
829                return NULL_NEIGHBOR;
830                }
831
832        unsigned uParent = GetParent(uNodeIndex);
833        if (GetLeft(uParent) == uNodeIndex)
834                {
835#if     TRACE
836                Log(">> Is left branch, return parent=%u\n", uParent);
837#endif
838                return uParent;
839                }
840
841        uNodeIndex = GetLeft(uParent);
842#if     TRACE
843                Log(">> Descend right from left sibling=%u ... ", uNodeIndex);
844#endif
845        while (!IsLeaf(uNodeIndex))
846                uNodeIndex = GetRight(uNodeIndex);
847
848#if     TRACE
849        Log("bottom out at leaf=%u\n", uNodeIndex);
850#endif
851        return uNodeIndex;
852        }
853
854void Tree::UnrootByDeletingRoot()
855        {
856        assert(IsRooted());
857        assert(m_uNodeCount >= 3);
858
859        const unsigned uLeft = GetLeft(m_uRootNodeIndex);
860        const unsigned uRight = GetRight(m_uRootNodeIndex);
861
862        m_uNeighbor1[uLeft] = uRight;
863        m_uNeighbor1[uRight] = uLeft;
864
865        bool bHasEdgeLength = HasEdgeLength(m_uRootNodeIndex, uLeft) &&
866          HasEdgeLength(m_uRootNodeIndex, uRight);
867        if (bHasEdgeLength)
868                {
869                double dEdgeLength = GetEdgeLength(m_uRootNodeIndex, uLeft) +
870                  GetEdgeLength(m_uRootNodeIndex, uRight);
871                m_dEdgeLength1[uLeft] = dEdgeLength;
872                m_dEdgeLength1[uRight] = dEdgeLength;
873                }
874
875// Remove root node entry from arrays
876        const unsigned uMoveCount = m_uNodeCount - m_uRootNodeIndex;
877        const unsigned uUnsBytes = uMoveCount*sizeof(unsigned);
878        memmove(m_uNeighbor1 + m_uRootNodeIndex, m_uNeighbor1 + m_uRootNodeIndex + 1,
879          uUnsBytes);
880        memmove(m_uNeighbor2 + m_uRootNodeIndex, m_uNeighbor2 + m_uRootNodeIndex + 1,
881          uUnsBytes);
882        memmove(m_uNeighbor3 + m_uRootNodeIndex, m_uNeighbor3 + m_uRootNodeIndex + 1,
883          uUnsBytes);
884
885        const unsigned uDoubleBytes = uMoveCount*sizeof(double);
886        memmove(m_dEdgeLength1 + m_uRootNodeIndex, m_dEdgeLength1 + m_uRootNodeIndex + 1,
887          uDoubleBytes);
888        memmove(m_dEdgeLength2 + m_uRootNodeIndex, m_dEdgeLength2 + m_uRootNodeIndex + 1,
889          uDoubleBytes);
890        memmove(m_dEdgeLength3 + m_uRootNodeIndex, m_dEdgeLength3 + m_uRootNodeIndex + 1,
891          uDoubleBytes);
892
893        const unsigned uBoolBytes = uMoveCount*sizeof(bool);
894        memmove(m_bHasEdgeLength1 + m_uRootNodeIndex, m_bHasEdgeLength1 + m_uRootNodeIndex + 1,
895          uBoolBytes);
896        memmove(m_bHasEdgeLength2 + m_uRootNodeIndex, m_bHasEdgeLength2 + m_uRootNodeIndex + 1,
897          uBoolBytes);
898        memmove(m_bHasEdgeLength3 + m_uRootNodeIndex, m_bHasEdgeLength3 + m_uRootNodeIndex + 1,
899          uBoolBytes);
900
901        const unsigned uPtrBytes = uMoveCount*sizeof(char *);
902        memmove(m_ptrName + m_uRootNodeIndex, m_ptrName + m_uRootNodeIndex + 1, uPtrBytes);
903
904        --m_uNodeCount;
905        m_bRooted = false;
906
907// Fix up table entries
908        for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
909                {
910#define DEC(x)  if (x != NULL_NEIGHBOR && x > m_uRootNodeIndex) --x;
911                DEC(m_uNeighbor1[uNodeIndex])
912                DEC(m_uNeighbor2[uNodeIndex])
913                DEC(m_uNeighbor3[uNodeIndex])
914#undef  DEC
915                }
916
917        Validate();
918        }
919
920unsigned Tree::GetLeafParent(unsigned uNodeIndex) const
921        {
922        assert(IsLeaf(uNodeIndex));
923
924        if (IsRooted())
925                return GetParent(uNodeIndex);
926
927        if (m_uNeighbor1[uNodeIndex] != NULL_NEIGHBOR)
928                return m_uNeighbor1[uNodeIndex];
929        if (m_uNeighbor2[uNodeIndex] != NULL_NEIGHBOR)
930                return m_uNeighbor2[uNodeIndex];
931        return m_uNeighbor3[uNodeIndex];
932        }
933
934// TODO: This is not efficient for large trees, should cache.
935double Tree::GetNodeHeight(unsigned uNodeIndex) const
936        {
937        if (!IsRooted())
938                Quit("Tree::GetNodeHeight: undefined unless rooted tree");
939       
940        if (IsLeaf(uNodeIndex))
941                return 0.0;
942
943        if (m_bHasHeight[uNodeIndex])
944                return m_dHeight[uNodeIndex];
945
946        const unsigned uLeft = GetLeft(uNodeIndex);
947        const unsigned uRight = GetRight(uNodeIndex);
948        double dLeftLength = GetEdgeLength(uNodeIndex, uLeft);
949        double dRightLength = GetEdgeLength(uNodeIndex, uRight);
950
951        if (dLeftLength < 0)
952                dLeftLength = 0;
953        if (dRightLength < 0)
954                dRightLength = 0;
955
956        const double dLeftHeight = dLeftLength + GetNodeHeight(uLeft);
957        const double dRightHeight = dRightLength + GetNodeHeight(uRight);
958        const double dHeight = (dLeftHeight + dRightHeight)/2;
959        m_bHasHeight[uNodeIndex] = true;
960        m_dHeight[uNodeIndex] = dHeight;
961        return dHeight;
962        }
963
964unsigned Tree::GetNeighborSubscript(unsigned uNodeIndex, unsigned uNeighborIndex) const
965        {
966        assert(uNodeIndex < m_uNodeCount);
967        assert(uNeighborIndex < m_uNodeCount);
968        if (uNeighborIndex == m_uNeighbor1[uNodeIndex])
969                return 0;
970        if (uNeighborIndex == m_uNeighbor2[uNodeIndex])
971                return 1;
972        if (uNeighborIndex == m_uNeighbor3[uNodeIndex])
973                return 2;
974        return NULL_NEIGHBOR;
975        }
976
977unsigned Tree::GetNeighbor(unsigned uNodeIndex, unsigned uNeighborSubscript) const
978        {
979        switch (uNeighborSubscript)
980                {
981        case 0:
982                return m_uNeighbor1[uNodeIndex];
983        case 1:
984                return m_uNeighbor2[uNodeIndex];
985        case 2:
986                return m_uNeighbor3[uNodeIndex];
987                }
988        Quit("Tree::GetNeighbor, sub=%u", uNeighborSubscript);
989        return NULL_NEIGHBOR;
990        }
991
992// TODO: check if this is a performance issue, could cache a lookup table
993unsigned Tree::LeafIndexToNodeIndex(unsigned uLeafIndex) const
994        {
995        const unsigned uNodeCount = GetNodeCount();
996        unsigned uLeafCount = 0;
997        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
998                {
999                if (IsLeaf(uNodeIndex))
1000                        {
1001                        if (uLeafCount == uLeafIndex)
1002                                return uNodeIndex;
1003                        else
1004                                ++uLeafCount;
1005                        }
1006                }
1007        Quit("LeafIndexToNodeIndex: out of range");
1008        return 0;
1009        }
1010
1011unsigned Tree::GetLeafNodeIndex(const char *ptrName) const
1012        {
1013        const unsigned uNodeCount = GetNodeCount();
1014        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
1015                {
1016                if (!IsLeaf(uNodeIndex))
1017                        continue;
1018                const char *ptrLeafName = GetLeafName(uNodeIndex);
1019                if (0 == strcmp(ptrName, ptrLeafName))
1020                        return uNodeIndex;
1021                }
1022        Quit("Tree::GetLeafNodeIndex, name not found");
1023        return 0;
1024        }
1025
1026void Tree::Copy(const Tree &tree)
1027        {
1028        const unsigned uNodeCount = tree.GetNodeCount();
1029        InitCache(uNodeCount);
1030
1031        m_uNodeCount = uNodeCount;
1032
1033        const size_t UnsignedBytes = uNodeCount*sizeof(unsigned);
1034        const size_t DoubleBytes = uNodeCount*sizeof(double);
1035        const size_t BoolBytes = uNodeCount*sizeof(bool);
1036
1037        memcpy(m_uNeighbor1, tree.m_uNeighbor1, UnsignedBytes);
1038        memcpy(m_uNeighbor2, tree.m_uNeighbor2, UnsignedBytes);
1039        memcpy(m_uNeighbor3, tree.m_uNeighbor3, UnsignedBytes);
1040
1041        memcpy(m_Ids, tree.m_Ids, UnsignedBytes);
1042
1043        memcpy(m_dEdgeLength1, tree.m_dEdgeLength1, DoubleBytes);
1044        memcpy(m_dEdgeLength2, tree.m_dEdgeLength2, DoubleBytes);
1045        memcpy(m_dEdgeLength3, tree.m_dEdgeLength3, DoubleBytes);
1046
1047        memcpy(m_dHeight, tree.m_dHeight, DoubleBytes);
1048
1049        memcpy(m_bHasEdgeLength1, tree.m_bHasEdgeLength1, BoolBytes);
1050        memcpy(m_bHasEdgeLength2, tree.m_bHasEdgeLength2, BoolBytes);
1051        memcpy(m_bHasEdgeLength3, tree.m_bHasEdgeLength3, BoolBytes);
1052
1053        memcpy(m_bHasHeight, tree.m_bHasHeight, BoolBytes);
1054
1055        m_uRootNodeIndex = tree.m_uRootNodeIndex;
1056        m_bRooted = tree.m_bRooted;
1057
1058        for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)
1059                {
1060                if (tree.IsLeaf(uNodeIndex))
1061                        {
1062                        const char *ptrName = tree.GetLeafName(uNodeIndex);
1063                        m_ptrName[uNodeIndex] = strsave(ptrName);
1064                        }
1065                else
1066                        m_ptrName[uNodeIndex] = 0;
1067                }
1068
1069#if     DEBUG
1070        Validate();
1071#endif
1072        }
1073
1074// Create rooted tree from a vector description.
1075// Node indexes are 0..N-1 for leaves, N..2N-2 for
1076// internal nodes.
1077// Vector subscripts are i-N and have values for
1078// internal nodes only, but those values are node
1079// indexes 0..2N-2. So e.g. if N=6 and Left[2]=1,
1080// this means that the third internal node (node index 8)
1081// has the second leaf (node index 1) as its left child.
1082// uRoot gives the vector subscript of the root, so add N
1083// to get the node index.
1084void Tree::Create(unsigned uLeafCount, unsigned uRoot, const unsigned Left[],
1085  const unsigned Right[], const float LeftLength[], const float RightLength[],
1086  const unsigned LeafIds[], char **LeafNames)
1087        {
1088        Clear();
1089
1090        m_uNodeCount = 2*uLeafCount - 1;
1091        InitCache(m_uNodeCount);
1092
1093        for (unsigned uNodeIndex = 0; uNodeIndex < uLeafCount; ++uNodeIndex)
1094                {
1095                m_Ids[uNodeIndex] = LeafIds[uNodeIndex];
1096                m_ptrName[uNodeIndex] = strsave(LeafNames[uNodeIndex]);
1097                }
1098
1099        for (unsigned uNodeIndex = uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)
1100                {
1101                unsigned v = uNodeIndex - uLeafCount;
1102                unsigned uLeft = Left[v];
1103                unsigned uRight = Right[v];
1104                float fLeft = LeftLength[v];
1105                float fRight = RightLength[v];
1106
1107                m_uNeighbor2[uNodeIndex] = uLeft;
1108                m_uNeighbor3[uNodeIndex] = uRight;
1109
1110                m_bHasEdgeLength2[uNodeIndex] = true;
1111                m_bHasEdgeLength3[uNodeIndex] = true;
1112
1113                m_dEdgeLength2[uNodeIndex] = fLeft;
1114                m_dEdgeLength3[uNodeIndex] = fRight;
1115
1116                m_uNeighbor1[uLeft] = uNodeIndex;
1117                m_uNeighbor1[uRight] = uNodeIndex;
1118
1119                m_dEdgeLength1[uLeft] = fLeft;
1120                m_dEdgeLength1[uRight] = fRight;
1121
1122                m_bHasEdgeLength1[uLeft] = true;
1123                m_bHasEdgeLength1[uRight] = true;
1124                }
1125
1126        m_bRooted = true;
1127        m_uRootNodeIndex = uRoot + uLeafCount;
1128
1129        Validate();
1130        }
Note: See TracBrowser for help on using the repository browser.