← Back to branch summary

~ubuntu-branches/ubuntu/karmic/ugene/karmic

« back to all changes in this revision

Viewing changes to src/plugins_3rdparty/umuscle/src/muscle/phy2.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Morten Kjeldgaard
  • Date: 2009-06-16 13:20:00 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20090616132000-a0eezpre9uvx91a5
Tags: 1.4.2+repack-0ubuntu1
New upstream release

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/plugins_3rdparty/umuscle/src/muscle/phy2.cpp
1
 
/*****************************************************************
2
 
* Unipro UGENE - Integrated Bioinformatics Suite
3
 
* Copyright (C) 2008 Unipro, Russia (http://ugene.unipro.ru)
4
 
* All Rights Reserved
5
 
* 
6
 
*     This source code is distributed under the terms of the
7
 
*     GNU General Public License. See the files COPYING and LICENSE
8
 
*     for details.
9
 
*****************************************************************/
10
 
 
11
 
#include "muscle.h"
12
 
#include "tree.h"
13
 
 
14
 
#define TRACE   0
15
 
 
16
 
// Return false when done
17
 
bool PhyEnumEdges(const Tree &tree, PhyEnumEdgeState &ES)
18
 
        {
19
 
        unsigned uNode1 = uInsane;
20
 
 
21
 
        if (!ES.m_bInit)
22
 
                {
23
 
                if (tree.GetNodeCount() <= 1)
24
 
                        {
25
 
                        ES.m_uNodeIndex1 = NULL_NEIGHBOR;
26
 
                        ES.m_uNodeIndex2 = NULL_NEIGHBOR;
27
 
                        return false;
28
 
                        }
29
 
                uNode1 = tree.FirstDepthFirstNode();
30
 
                ES.m_bInit = true;
31
 
                }
32
 
        else
33
 
                {
34
 
                uNode1 = tree.NextDepthFirstNode(ES.m_uNodeIndex1);
35
 
                if (NULL_NEIGHBOR == uNode1)
36
 
                        return false;
37
 
                if (tree.IsRooted() && tree.IsRoot(uNode1))
38
 
                        {
39
 
                        uNode1 = tree.NextDepthFirstNode(uNode1);
40
 
                        if (NULL_NEIGHBOR == uNode1)
41
 
                                return false;
42
 
                        }
43
 
                }
44
 
        unsigned uNode2 = tree.GetParent(uNode1);
45
 
 
46
 
        ES.m_uNodeIndex1 = uNode1;
47
 
        ES.m_uNodeIndex2 = uNode2;
48
 
        return true;
49
 
        }
50
 
 
51
 
bool PhyEnumEdgesR(const Tree &tree, PhyEnumEdgeState &ES)
52
 
        {
53
 
        unsigned uNode1 = uInsane;
54
 
 
55
 
        if (!ES.m_bInit)
56
 
                {
57
 
                if (tree.GetNodeCount() <= 1)
58
 
                        {
59
 
                        ES.m_uNodeIndex1 = NULL_NEIGHBOR;
60
 
                        ES.m_uNodeIndex2 = NULL_NEIGHBOR;
61
 
                        return false;
62
 
                        }
63
 
                uNode1 = tree.FirstDepthFirstNodeR();
64
 
                ES.m_bInit = true;
65
 
                }
66
 
        else
67
 
                {
68
 
                uNode1 = tree.NextDepthFirstNodeR(ES.m_uNodeIndex1);
69
 
                if (NULL_NEIGHBOR == uNode1)
70
 
                        return false;
71
 
                if (tree.IsRooted() && tree.IsRoot(uNode1))
72
 
                        {
73
 
                        uNode1 = tree.NextDepthFirstNode(uNode1);
74
 
                        if (NULL_NEIGHBOR == uNode1)
75
 
                                return false;
76
 
                        }
77
 
                }
78
 
        unsigned uNode2 = tree.GetParent(uNode1);
79
 
 
80
 
        ES.m_uNodeIndex1 = uNode1;
81
 
        ES.m_uNodeIndex2 = uNode2;
82
 
        return true;
83
 
        }
84
 
 
85
 
static void GetLeavesSubtree(const Tree &tree, unsigned uNodeIndex1,
86
 
  const unsigned uNodeIndex2, unsigned Leaves[], unsigned *ptruCount)
87
 
        {
88
 
        if (tree.IsLeaf(uNodeIndex1))
89
 
                {
90
 
                Leaves[*ptruCount] = uNodeIndex1;
91
 
                ++(*ptruCount);
92
 
                return;
93
 
                }
94
 
 
95
 
        const unsigned uLeft = tree.GetFirstNeighbor(uNodeIndex1, uNodeIndex2);
96
 
        const unsigned uRight = tree.GetSecondNeighbor(uNodeIndex1, uNodeIndex2);
97
 
        if (NULL_NEIGHBOR != uLeft)
98
 
                GetLeavesSubtree(tree, uLeft, uNodeIndex1, Leaves, ptruCount);
99
 
        if (NULL_NEIGHBOR != uRight)
100
 
                GetLeavesSubtree(tree, uRight, uNodeIndex1, Leaves, ptruCount);
101
 
        }
102
 
 
103
 
static void PhyGetLeaves(const Tree &tree, unsigned uNodeIndex1, unsigned uNodeIndex2,
104
 
  unsigned Leaves[], unsigned *ptruCount)
105
 
        {
106
 
        *ptruCount = 0;
107
 
        GetLeavesSubtree(tree, uNodeIndex1, uNodeIndex2, Leaves, ptruCount);
108
 
        }
109
 
 
110
 
bool PhyEnumBiParts(const Tree &tree, PhyEnumEdgeState &ES,
111
 
  unsigned Leaves1[], unsigned *ptruCount1,
112
 
  unsigned Leaves2[], unsigned *ptruCount2)
113
 
        {
114
 
        bool bOk = PhyEnumEdges(tree, ES);
115
 
        if (!bOk)
116
 
                {
117
 
                *ptruCount1 = 0;
118
 
                *ptruCount2 = 0;
119
 
                return false;
120
 
                }
121
 
 
122
 
// Special case: in a rooted tree, both edges from the root
123
 
// give the same bipartition, so skip one of them.
124
 
        if (tree.IsRooted() && tree.IsRoot(ES.m_uNodeIndex2)
125
 
          && tree.GetRight(ES.m_uNodeIndex2) == ES.m_uNodeIndex1)
126
 
                {
127
 
                bOk = PhyEnumEdges(tree, ES);
128
 
                if (!bOk)
129
 
                        return false;
130
 
                }
131
 
 
132
 
        PhyGetLeaves(tree, ES.m_uNodeIndex1, ES.m_uNodeIndex2, Leaves1, ptruCount1);
133
 
        PhyGetLeaves(tree, ES.m_uNodeIndex2, ES.m_uNodeIndex1, Leaves2, ptruCount2);
134
 
 
135
 
        if (*ptruCount1 + *ptruCount2 != tree.GetLeafCount())
136
 
                Quit("PhyEnumBiParts %u + %u != %u",
137
 
                  *ptruCount1, *ptruCount2, tree.GetLeafCount());
138
 
#if     DEBUG
139
 
        {
140
 
        for (unsigned i = 0; i < *ptruCount1; ++i)
141
 
                {
142
 
                if (!tree.IsLeaf(Leaves1[i]))
143
 
                        Quit("PhyEnumByParts: not leaf");
144
 
                for (unsigned j = 0; j < *ptruCount2; ++j)
145
 
                        {
146
 
                        if (!tree.IsLeaf(Leaves2[j]))
147
 
                                Quit("PhyEnumByParts: not leaf");
148
 
                        if (Leaves1[i] == Leaves2[j])
149
 
                                Quit("PhyEnumByParts: dupe");
150
 
                        }
151
 
                }
152
 
        }
153
 
#endif
154
 
 
155
 
        return true;
156
 
        }
157
 
 
158
 
#if     0
159
 
void TestBiPart()
160
 
        {
161
 
        SetListFileName("c:\\tmp\\lobster.log", false);
162
 
        Tree tree;
163
 
        TextFile fileIn("c:\\tmp\\test.phy");
164
 
        tree.FromFile(fileIn);
165
 
        tree.LogMe();
166
 
 
167
 
        const unsigned uNodeCount = tree.GetNodeCount();
168
 
        unsigned *Leaves1 = new unsigned[uNodeCount];
169
 
        unsigned *Leaves2 = new unsigned[uNodeCount];
170
 
 
171
 
        PhyEnumEdgeState ES;
172
 
        bool bDone = false;
173
 
        for (;;)
174
 
                {
175
 
                unsigned uCount1 = uInsane;
176
 
                unsigned uCount2 = uInsane;
177
 
                bool bOk = PhyEnumBiParts(tree, ES, Leaves1, &uCount1, Leaves2, &uCount2);
178
 
                Log("PEBP=%d ES.Init=%d ES.ni1=%d ES.ni2=%d\n",
179
 
                  bOk,
180
 
                  ES.m_bInit,
181
 
                  ES.m_uNodeIndex1,
182
 
                  ES.m_uNodeIndex2);
183
 
                if (!bOk)
184
 
                        break;
185
 
                Log("\n");
186
 
                Log("Part1: ");
187
 
                for (unsigned n = 0; n < uCount1; ++n)
188
 
                        Log(" %d(%s)", Leaves1[n], tree.GetLeafName(Leaves1[n]));
189
 
                Log("\n");
190
 
                Log("Part2: ");
191
 
                for (unsigned n = 0; n < uCount2; ++n)
192
 
                        Log(" %d(%s)", Leaves2[n], tree.GetLeafName(Leaves2[n]));
193
 
                Log("\n");
194
 
                }
195
 
        }
196
 
#endif
197
 
 
198
 
static void GetLeavesSubtreeExcluding(const Tree &tree, unsigned uNodeIndex,
199
 
  unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
200
 
        {
201
 
        if (uNodeIndex == uExclude)
202
 
                return;
203
 
 
204
 
        if (tree.IsLeaf(uNodeIndex))
205
 
                {
206
 
                Leaves[*ptruCount] = uNodeIndex;
207
 
                ++(*ptruCount);
208
 
                return;
209
 
                }
210
 
 
211
 
        const unsigned uLeft = tree.GetLeft(uNodeIndex);
212
 
        const unsigned uRight = tree.GetRight(uNodeIndex);
213
 
        if (NULL_NEIGHBOR != uLeft)
214
 
                GetLeavesSubtreeExcluding(tree, uLeft, uExclude, Leaves, ptruCount);
215
 
        if (NULL_NEIGHBOR != uRight)
216
 
                GetLeavesSubtreeExcluding(tree, uRight, uExclude, Leaves, ptruCount);
217
 
        }
218
 
 
219
 
void GetLeavesExcluding(const Tree &tree, unsigned uNodeIndex,
220
 
  unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
221
 
        {
222
 
        *ptruCount = 0;
223
 
        GetLeavesSubtreeExcluding(tree, uNodeIndex, uExclude, Leaves, ptruCount);
224
 
        }
225
 
 
226
 
void GetInternalNodesInHeightOrder(const Tree &tree, unsigned NodeIndexes[])
227
 
        {
228
 
        const unsigned uNodeCount = tree.GetNodeCount();
229
 
        if (uNodeCount < 3)
230
 
                Quit("GetInternalNodesInHeightOrder: %u nodes, none are internal",
231
 
                  uNodeCount);
232
 
        const unsigned uInternalNodeCount = (uNodeCount - 1)/2;
233
 
        double *Heights = new double[uInternalNodeCount];
234
 
 
235
 
        unsigned uIndex = 0;
236
 
        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
237
 
                {
238
 
                if (tree.IsLeaf(uNodeIndex))
239
 
                        continue;
240
 
                NodeIndexes[uIndex] = uNodeIndex;
241
 
                Heights[uIndex] = tree.GetNodeHeight(uNodeIndex);
242
 
                ++uIndex;
243
 
                }
244
 
        if (uIndex != uInternalNodeCount)
245
 
                Quit("Internal error: GetInternalNodesInHeightOrder");
246
 
 
247
 
// Simple but slow bubble sort (probably don't care about speed here)
248
 
        bool bDone = false;
249
 
        while (!bDone)
250
 
                {
251
 
                bDone = true;
252
 
                for (unsigned i = 0; i < uInternalNodeCount - 1; ++i)
253
 
                        {
254
 
                        if (Heights[i] > Heights[i+1])
255
 
                                {
256
 
                                double dTmp = Heights[i];
257
 
                                Heights[i] = Heights[i+1];
258
 
                                Heights[i+1] = dTmp;
259
 
 
260
 
                                unsigned uTmp = NodeIndexes[i];
261
 
                                NodeIndexes[i] = NodeIndexes[i+1];
262
 
                                NodeIndexes[i+1] = uTmp;
263
 
                                bDone = false;
264
 
                                }
265
 
                        }
266
 
                }
267
 
#if     TRACE
268
 
        Log("Internal node index     Height\n");
269
 
        Log("-------------------   --------\n");
270
 
        //    1234567890123456789  123456789
271
 
        for (unsigned n = 0; n < uInternalNodeCount; ++n)
272
 
                Log("%19u  %9.3f\n", NodeIndexes[n], Heights[n]);
273
 
#endif
274
 
        delete[] Heights;
275
 
        }
276
 
 
277
 
void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)
278
 
        {
279
 
        const unsigned uNodeCount = tree.GetNodeCount();
280
 
        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
281
 
                {
282
 
                const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);
283
 
                for (unsigned n = 0; n < uNeighborCount; ++n)
284
 
                        {
285
 
                        const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);
286
 
                        if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))
287
 
                                continue;
288
 
                        if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)
289
 
                                tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);
290
 
                        }
291
 
                }
292
 
        }
 
1
/*****************************************************************
 
2
* Unipro UGENE - Integrated Bioinformatics Suite
 
3
* Copyright (C) 2008 Unipro, Russia (http://ugene.unipro.ru)
 
4
* All Rights Reserved
 
5
* 
 
6
*     This source code is distributed under the terms of the
 
7
*     GNU General Public License. See the files COPYING and LICENSE
 
8
*     for details.
 
9
*****************************************************************/
 
10
 
 
11
#include "muscle.h"
 
12
#include "tree.h"
 
13
 
 
14
#define TRACE   0
 
15
 
 
16
// Return false when done
 
17
bool PhyEnumEdges(const Tree &tree, PhyEnumEdgeState &ES)
 
18
        {
 
19
        unsigned uNode1 = uInsane;
 
20
 
 
21
        if (!ES.m_bInit)
 
22
                {
 
23
                if (tree.GetNodeCount() <= 1)
 
24
                        {
 
25
                        ES.m_uNodeIndex1 = NULL_NEIGHBOR;
 
26
                        ES.m_uNodeIndex2 = NULL_NEIGHBOR;
 
27
                        return false;
 
28
                        }
 
29
                uNode1 = tree.FirstDepthFirstNode();
 
30
                ES.m_bInit = true;
 
31
                }
 
32
        else
 
33
                {
 
34
                uNode1 = tree.NextDepthFirstNode(ES.m_uNodeIndex1);
 
35
                if (NULL_NEIGHBOR == uNode1)
 
36
                        return false;
 
37
                if (tree.IsRooted() && tree.IsRoot(uNode1))
 
38
                        {
 
39
                        uNode1 = tree.NextDepthFirstNode(uNode1);
 
40
                        if (NULL_NEIGHBOR == uNode1)
 
41
                                return false;
 
42
                        }
 
43
                }
 
44
        unsigned uNode2 = tree.GetParent(uNode1);
 
45
 
 
46
        ES.m_uNodeIndex1 = uNode1;
 
47
        ES.m_uNodeIndex2 = uNode2;
 
48
        return true;
 
49
        }
 
50
 
 
51
bool PhyEnumEdgesR(const Tree &tree, PhyEnumEdgeState &ES)
 
52
        {
 
53
        unsigned uNode1 = uInsane;
 
54
 
 
55
        if (!ES.m_bInit)
 
56
                {
 
57
                if (tree.GetNodeCount() <= 1)
 
58
                        {
 
59
                        ES.m_uNodeIndex1 = NULL_NEIGHBOR;
 
60
                        ES.m_uNodeIndex2 = NULL_NEIGHBOR;
 
61
                        return false;
 
62
                        }
 
63
                uNode1 = tree.FirstDepthFirstNodeR();
 
64
                ES.m_bInit = true;
 
65
                }
 
66
        else
 
67
                {
 
68
                uNode1 = tree.NextDepthFirstNodeR(ES.m_uNodeIndex1);
 
69
                if (NULL_NEIGHBOR == uNode1)
 
70
                        return false;
 
71
                if (tree.IsRooted() && tree.IsRoot(uNode1))
 
72
                        {
 
73
                        uNode1 = tree.NextDepthFirstNode(uNode1);
 
74
                        if (NULL_NEIGHBOR == uNode1)
 
75
                                return false;
 
76
                        }
 
77
                }
 
78
        unsigned uNode2 = tree.GetParent(uNode1);
 
79
 
 
80
        ES.m_uNodeIndex1 = uNode1;
 
81
        ES.m_uNodeIndex2 = uNode2;
 
82
        return true;
 
83
        }
 
84
 
 
85
static void GetLeavesSubtree(const Tree &tree, unsigned uNodeIndex1,
 
86
  const unsigned uNodeIndex2, unsigned Leaves[], unsigned *ptruCount)
 
87
        {
 
88
        if (tree.IsLeaf(uNodeIndex1))
 
89
                {
 
90
                Leaves[*ptruCount] = uNodeIndex1;
 
91
                ++(*ptruCount);
 
92
                return;
 
93
                }
 
94
 
 
95
        const unsigned uLeft = tree.GetFirstNeighbor(uNodeIndex1, uNodeIndex2);
 
96
        const unsigned uRight = tree.GetSecondNeighbor(uNodeIndex1, uNodeIndex2);
 
97
        if (NULL_NEIGHBOR != uLeft)
 
98
                GetLeavesSubtree(tree, uLeft, uNodeIndex1, Leaves, ptruCount);
 
99
        if (NULL_NEIGHBOR != uRight)
 
100
                GetLeavesSubtree(tree, uRight, uNodeIndex1, Leaves, ptruCount);
 
101
        }
 
102
 
 
103
static void PhyGetLeaves(const Tree &tree, unsigned uNodeIndex1, unsigned uNodeIndex2,
 
104
  unsigned Leaves[], unsigned *ptruCount)
 
105
        {
 
106
        *ptruCount = 0;
 
107
        GetLeavesSubtree(tree, uNodeIndex1, uNodeIndex2, Leaves, ptruCount);
 
108
        }
 
109
 
 
110
bool PhyEnumBiParts(const Tree &tree, PhyEnumEdgeState &ES,
 
111
  unsigned Leaves1[], unsigned *ptruCount1,
 
112
  unsigned Leaves2[], unsigned *ptruCount2)
 
113
        {
 
114
        bool bOk = PhyEnumEdges(tree, ES);
 
115
        if (!bOk)
 
116
                {
 
117
                *ptruCount1 = 0;
 
118
                *ptruCount2 = 0;
 
119
                return false;
 
120
                }
 
121
 
 
122
// Special case: in a rooted tree, both edges from the root
 
123
// give the same bipartition, so skip one of them.
 
124
        if (tree.IsRooted() && tree.IsRoot(ES.m_uNodeIndex2)
 
125
          && tree.GetRight(ES.m_uNodeIndex2) == ES.m_uNodeIndex1)
 
126
                {
 
127
                bOk = PhyEnumEdges(tree, ES);
 
128
                if (!bOk)
 
129
                        return false;
 
130
                }
 
131
 
 
132
        PhyGetLeaves(tree, ES.m_uNodeIndex1, ES.m_uNodeIndex2, Leaves1, ptruCount1);
 
133
        PhyGetLeaves(tree, ES.m_uNodeIndex2, ES.m_uNodeIndex1, Leaves2, ptruCount2);
 
134
 
 
135
        if (*ptruCount1 + *ptruCount2 != tree.GetLeafCount())
 
136
                Quit("PhyEnumBiParts %u + %u != %u",
 
137
                  *ptruCount1, *ptruCount2, tree.GetLeafCount());
 
138
#if     DEBUG
 
139
        {
 
140
        for (unsigned i = 0; i < *ptruCount1; ++i)
 
141
                {
 
142
                if (!tree.IsLeaf(Leaves1[i]))
 
143
                        Quit("PhyEnumByParts: not leaf");
 
144
                for (unsigned j = 0; j < *ptruCount2; ++j)
 
145
                        {
 
146
                        if (!tree.IsLeaf(Leaves2[j]))
 
147
                                Quit("PhyEnumByParts: not leaf");
 
148
                        if (Leaves1[i] == Leaves2[j])
 
149
                                Quit("PhyEnumByParts: dupe");
 
150
                        }
 
151
                }
 
152
        }
 
153
#endif
 
154
 
 
155
        return true;
 
156
        }
 
157
 
 
158
#if     0
 
159
void TestBiPart()
 
160
        {
 
161
        SetListFileName("c:\\tmp\\lobster.log", false);
 
162
        Tree tree;
 
163
        TextFile fileIn("c:\\tmp\\test.phy");
 
164
        tree.FromFile(fileIn);
 
165
        tree.LogMe();
 
166
 
 
167
        const unsigned uNodeCount = tree.GetNodeCount();
 
168
        unsigned *Leaves1 = new unsigned[uNodeCount];
 
169
        unsigned *Leaves2 = new unsigned[uNodeCount];
 
170
 
 
171
        PhyEnumEdgeState ES;
 
172
        bool bDone = false;
 
173
        for (;;)
 
174
                {
 
175
                unsigned uCount1 = uInsane;
 
176
                unsigned uCount2 = uInsane;
 
177
                bool bOk = PhyEnumBiParts(tree, ES, Leaves1, &uCount1, Leaves2, &uCount2);
 
178
                Log("PEBP=%d ES.Init=%d ES.ni1=%d ES.ni2=%d\n",
 
179
                  bOk,
 
180
                  ES.m_bInit,
 
181
                  ES.m_uNodeIndex1,
 
182
                  ES.m_uNodeIndex2);
 
183
                if (!bOk)
 
184
                        break;
 
185
                Log("\n");
 
186
                Log("Part1: ");
 
187
                for (unsigned n = 0; n < uCount1; ++n)
 
188
                        Log(" %d(%s)", Leaves1[n], tree.GetLeafName(Leaves1[n]));
 
189
                Log("\n");
 
190
                Log("Part2: ");
 
191
                for (unsigned n = 0; n < uCount2; ++n)
 
192
                        Log(" %d(%s)", Leaves2[n], tree.GetLeafName(Leaves2[n]));
 
193
                Log("\n");
 
194
                }
 
195
        }
 
196
#endif
 
197
 
 
198
static void GetLeavesSubtreeExcluding(const Tree &tree, unsigned uNodeIndex,
 
199
  unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
 
200
        {
 
201
        if (uNodeIndex == uExclude)
 
202
                return;
 
203
 
 
204
        if (tree.IsLeaf(uNodeIndex))
 
205
                {
 
206
                Leaves[*ptruCount] = uNodeIndex;
 
207
                ++(*ptruCount);
 
208
                return;
 
209
                }
 
210
 
 
211
        const unsigned uLeft = tree.GetLeft(uNodeIndex);
 
212
        const unsigned uRight = tree.GetRight(uNodeIndex);
 
213
        if (NULL_NEIGHBOR != uLeft)
 
214
                GetLeavesSubtreeExcluding(tree, uLeft, uExclude, Leaves, ptruCount);
 
215
        if (NULL_NEIGHBOR != uRight)
 
216
                GetLeavesSubtreeExcluding(tree, uRight, uExclude, Leaves, ptruCount);
 
217
        }
 
218
 
 
219
void GetLeavesExcluding(const Tree &tree, unsigned uNodeIndex,
 
220
  unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)
 
221
        {
 
222
        *ptruCount = 0;
 
223
        GetLeavesSubtreeExcluding(tree, uNodeIndex, uExclude, Leaves, ptruCount);
 
224
        }
 
225
 
 
226
void GetInternalNodesInHeightOrder(const Tree &tree, unsigned NodeIndexes[])
 
227
        {
 
228
        const unsigned uNodeCount = tree.GetNodeCount();
 
229
        if (uNodeCount < 3)
 
230
                Quit("GetInternalNodesInHeightOrder: %u nodes, none are internal",
 
231
                  uNodeCount);
 
232
        const unsigned uInternalNodeCount = (uNodeCount - 1)/2;
 
233
        double *Heights = new double[uInternalNodeCount];
 
234
 
 
235
        unsigned uIndex = 0;
 
236
        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
 
237
                {
 
238
                if (tree.IsLeaf(uNodeIndex))
 
239
                        continue;
 
240
                NodeIndexes[uIndex] = uNodeIndex;
 
241
                Heights[uIndex] = tree.GetNodeHeight(uNodeIndex);
 
242
                ++uIndex;
 
243
                }
 
244
        if (uIndex != uInternalNodeCount)
 
245
                Quit("Internal error: GetInternalNodesInHeightOrder");
 
246
 
 
247
// Simple but slow bubble sort (probably don't care about speed here)
 
248
        bool bDone = false;
 
249
        while (!bDone)
 
250
                {
 
251
                bDone = true;
 
252
                for (unsigned i = 0; i < uInternalNodeCount - 1; ++i)
 
253
                        {
 
254
                        if (Heights[i] > Heights[i+1])
 
255
                                {
 
256
                                double dTmp = Heights[i];
 
257
                                Heights[i] = Heights[i+1];
 
258
                                Heights[i+1] = dTmp;
 
259
 
 
260
                                unsigned uTmp = NodeIndexes[i];
 
261
                                NodeIndexes[i] = NodeIndexes[i+1];
 
262
                                NodeIndexes[i+1] = uTmp;
 
263
                                bDone = false;
 
264
                                }
 
265
                        }
 
266
                }
 
267
#if     TRACE
 
268
        Log("Internal node index     Height\n");
 
269
        Log("-------------------   --------\n");
 
270
        //    1234567890123456789  123456789
 
271
        for (unsigned n = 0; n < uInternalNodeCount; ++n)
 
272
                Log("%19u  %9.3f\n", NodeIndexes[n], Heights[n]);
 
273
#endif
 
274
        delete[] Heights;
 
275
        }
 
276
 
 
277
void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)
 
278
        {
 
279
        const unsigned uNodeCount = tree.GetNodeCount();
 
280
        for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
 
281
                {
 
282
                const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);
 
283
                for (unsigned n = 0; n < uNeighborCount; ++n)
 
284
                        {
 
285
                        const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);
 
286
                        if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))
 
287
                                continue;
 
288
                        if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)
 
289
                                tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);
 
290
                        }
 
291
                }
 
292
        }