7 void ClusterByHeight(const Tree &tree, double dMaxHeight, unsigned Subtrees[],
\r
8 unsigned *ptruSubtreeCount)
\r
10 if (!tree.IsRooted())
\r
11 Quit("ClusterByHeight: requires rooted tree");
\r
14 Log("ClusterByHeight, max height=%g\n", dMaxHeight);
\r
17 unsigned uSubtreeCount = 0;
\r
18 const unsigned uNodeCount = tree.GetNodeCount();
\r
19 for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)
\r
21 if (tree.IsRoot(uNodeIndex))
\r
23 unsigned uParent = tree.GetParent(uNodeIndex);
\r
24 double dHeight = tree.GetNodeHeight(uNodeIndex);
\r
25 double dParentHeight = tree.GetNodeHeight(uParent);
\r
28 Log("Node %3u Height %5.2f ParentHeight %5.2f\n",
\r
29 uNodeIndex, dHeight, dParentHeight);
\r
31 if (dParentHeight > dMaxHeight && dHeight <= dMaxHeight)
\r
33 Subtrees[uSubtreeCount] = uNodeIndex;
\r
35 Log("Subtree[%u]=%u\n", uSubtreeCount, uNodeIndex);
\r
40 *ptruSubtreeCount = uSubtreeCount;
\r
43 static void ClusterBySubfamCount_Iteration(const Tree &tree, unsigned Subfams[],
\r
46 // Find highest child node of current set of subfamilies.
\r
47 double dHighestHeight = -1e20;
\r
48 int iParentSubscript = -1;
\r
50 for (int n = 0; n < (int) uCount; ++n)
\r
52 const unsigned uNodeIndex = Subfams[n];
\r
53 if (tree.IsLeaf(uNodeIndex))
\r
56 const unsigned uLeft = tree.GetLeft(uNodeIndex);
\r
57 const double dHeightLeft = tree.GetNodeHeight(uLeft);
\r
58 if (dHeightLeft > dHighestHeight)
\r
60 dHighestHeight = dHeightLeft;
\r
61 iParentSubscript = n;
\r
64 const unsigned uRight = tree.GetRight(uNodeIndex);
\r
65 const double dHeightRight = tree.GetNodeHeight(uRight);
\r
66 if (dHeightRight > dHighestHeight)
\r
68 dHighestHeight = dHeightRight;
\r
69 iParentSubscript = n;
\r
73 if (-1 == iParentSubscript)
\r
74 Quit("CBSFCIter: failed to find highest child");
\r
76 const unsigned uNodeIndex = Subfams[iParentSubscript];
\r
77 const unsigned uLeft = tree.GetLeft(uNodeIndex);
\r
78 const unsigned uRight = tree.GetRight(uNodeIndex);
\r
80 // Delete parent by replacing with left child
\r
81 Subfams[iParentSubscript] = uLeft;
\r
83 // Append right child to list
\r
84 Subfams[uCount] = uRight;
\r
88 Log("Iter %3u:", uCount);
\r
89 for (unsigned n = 0; n < uCount; ++n)
\r
90 Log(" %u", Subfams[n]);
\r
96 // Divide a tree containing N leaves into k families by
\r
97 // cutting the tree at a horizontal line at some height.
\r
98 // Each internal node defines a height for the cut,
\r
99 // considering all internal nodes enumerates all distinct
\r
100 // cuts. Visit internal nodes in decreasing order of height.
\r
101 // Visiting the node corresponds to moving the horizontal
\r
102 // line down to cut the tree at the height of that node.
\r
103 // We consider the cut to be "infinitestimally below"
\r
104 // the node, so the effect is to remove the current node
\r
105 // from the list of subfamilies and add its two children.
\r
106 // We must visit a parent before its children (so care may
\r
107 // be needed to handle zero edge lengths properly).
\r
108 // We assume that N is small, and write dumb O(N^2) code.
\r
109 // More efficient strategies are possible for large N
\r
110 // by maintaining a list of nodes sorted by height.
\r
111 void ClusterBySubfamCount(const Tree &tree, unsigned uSubfamCount,
\r
112 unsigned Subfams[], unsigned *ptruSubfamCount)
\r
114 const unsigned uNodeCount = tree.GetNodeCount();
\r
115 const unsigned uLeafCount = (uNodeCount + 1)/2;
\r
117 // Special case: empty tree
\r
118 if (0 == uNodeCount)
\r
120 *ptruSubfamCount = 0;
\r
124 // Special case: more subfamilies than leaves
\r
125 if (uSubfamCount >= uLeafCount)
\r
127 for (unsigned n = 0; n < uLeafCount; ++n)
\r
129 *ptruSubfamCount = uLeafCount;
\r
133 // Initialize list of subfamilies to be root
\r
134 Subfams[0] = tree.GetRootNodeIndex();
\r
137 for (unsigned i = 1; i < uSubfamCount; ++i)
\r
138 ClusterBySubfamCount_Iteration(tree, Subfams, i);
\r
140 *ptruSubfamCount = uSubfamCount;
\r
143 static void GetLeavesRecurse(const Tree &tree, unsigned uNodeIndex,
\r
144 unsigned Leaves[], unsigned &uLeafCount /* in-out */)
\r
146 if (tree.IsLeaf(uNodeIndex))
\r
148 Leaves[uLeafCount] = uNodeIndex;
\r
153 const unsigned uLeft = tree.GetLeft(uNodeIndex);
\r
154 const unsigned uRight = tree.GetRight(uNodeIndex);
\r
156 GetLeavesRecurse(tree, uLeft, Leaves, uLeafCount);
\r
157 GetLeavesRecurse(tree, uRight, Leaves, uLeafCount);
\r
160 void GetLeaves(const Tree &tree, unsigned uNodeIndex, unsigned Leaves[],
\r
161 unsigned *ptruLeafCount)
\r
163 unsigned uLeafCount = 0;
\r
164 GetLeavesRecurse(tree, uNodeIndex, Leaves, uLeafCount);
\r
165 *ptruLeafCount = uLeafCount;
\r
168 void Tree::PruneTree(const Tree &tree, unsigned Subfams[],
\r
169 unsigned uSubfamCount)
\r
171 if (!tree.IsRooted())
\r
172 Quit("Tree::PruneTree: requires rooted tree");
\r
176 m_uNodeCount = 2*uSubfamCount - 1;
\r
177 InitCache(m_uNodeCount);
\r
179 const unsigned uUnprunedNodeCount = tree.GetNodeCount();
\r
181 unsigned *uUnprunedToPrunedIndex = new unsigned[uUnprunedNodeCount];
\r
182 unsigned *uPrunedToUnprunedIndex = new unsigned[m_uNodeCount];
\r
184 for (unsigned n = 0; n < uUnprunedNodeCount; ++n)
\r
185 uUnprunedToPrunedIndex[n] = NULL_NEIGHBOR;
\r
187 for (unsigned n = 0; n < m_uNodeCount; ++n)
\r
188 uPrunedToUnprunedIndex[n] = NULL_NEIGHBOR;
\r
190 // Create mapping between unpruned and pruned node indexes
\r
191 unsigned uInternalNodeIndex = uSubfamCount;
\r
192 for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
\r
194 unsigned uUnprunedNodeIndex = Subfams[uSubfamIndex];
\r
195 uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uSubfamIndex;
\r
196 uPrunedToUnprunedIndex[uSubfamIndex] = uUnprunedNodeIndex;
\r
199 uUnprunedNodeIndex = tree.GetParent(uUnprunedNodeIndex);
\r
200 if (tree.IsRoot(uUnprunedNodeIndex))
\r
203 // Already visited this node?
\r
204 if (NULL_NEIGHBOR != uUnprunedToPrunedIndex[uUnprunedNodeIndex])
\r
207 uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uInternalNodeIndex;
\r
208 uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedNodeIndex;
\r
210 ++uInternalNodeIndex;
\r
214 const unsigned uUnprunedRootIndex = tree.GetRootNodeIndex();
\r
215 uUnprunedToPrunedIndex[uUnprunedRootIndex] = uInternalNodeIndex;
\r
216 uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedRootIndex;
\r
220 Log("Pruned to unpruned:\n");
\r
221 for (unsigned i = 0; i < m_uNodeCount; ++i)
\r
222 Log(" [%u]=%u", i, uPrunedToUnprunedIndex[i]);
\r
224 Log("Unpruned to pruned:\n");
\r
225 for (unsigned i = 0; i < uUnprunedNodeCount; ++i)
\r
227 unsigned n = uUnprunedToPrunedIndex[i];
\r
228 if (n != NULL_NEIGHBOR)
\r
229 Log(" [%u]=%u", i, n);
\r
235 if (uInternalNodeIndex != m_uNodeCount - 1)
\r
236 Quit("Tree::PruneTree, Internal error");
\r
238 // Nodes 0, 1 ... are the leaves
\r
239 for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)
\r
242 sprintf(szName, "Subfam_%u", uSubfamIndex + 1);
\r
243 m_ptrName[uSubfamIndex] = strsave(szName);
\r
246 for (unsigned uPrunedNodeIndex = uSubfamCount; uPrunedNodeIndex < m_uNodeCount;
\r
247 ++uPrunedNodeIndex)
\r
249 unsigned uUnprunedNodeIndex = uPrunedToUnprunedIndex[uPrunedNodeIndex];
\r
251 const unsigned uUnprunedLeft = tree.GetLeft(uUnprunedNodeIndex);
\r
252 const unsigned uUnprunedRight = tree.GetRight(uUnprunedNodeIndex);
\r
254 const unsigned uPrunedLeft = uUnprunedToPrunedIndex[uUnprunedLeft];
\r
255 const unsigned uPrunedRight = uUnprunedToPrunedIndex[uUnprunedRight];
\r
257 const double dLeftLength =
\r
258 tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedLeft);
\r
259 const double dRightLength =
\r
260 tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedRight);
\r
262 m_uNeighbor2[uPrunedNodeIndex] = uPrunedLeft;
\r
263 m_uNeighbor3[uPrunedNodeIndex] = uPrunedRight;
\r
265 m_dEdgeLength1[uPrunedLeft] = dLeftLength;
\r
266 m_dEdgeLength1[uPrunedRight] = dRightLength;
\r
268 m_uNeighbor1[uPrunedLeft] = uPrunedNodeIndex;
\r
269 m_uNeighbor1[uPrunedRight] = uPrunedNodeIndex;
\r
271 m_bHasEdgeLength1[uPrunedLeft] = true;
\r
272 m_bHasEdgeLength1[uPrunedRight] = true;
\r
274 m_dEdgeLength2[uPrunedNodeIndex] = dLeftLength;
\r
275 m_dEdgeLength3[uPrunedNodeIndex] = dRightLength;
\r
277 m_bHasEdgeLength2[uPrunedNodeIndex] = true;
\r
278 m_bHasEdgeLength3[uPrunedNodeIndex] = true;
\r
281 m_uRootNodeIndex = uUnprunedToPrunedIndex[uUnprunedRootIndex];
\r
287 delete[] uUnprunedToPrunedIndex;
\r
290 void LeafIndexesToIds(const Tree &tree, const unsigned Leaves[], unsigned uCount,
\r
293 for (unsigned n = 0; n < uCount; ++n)
\r
294 Ids[n] = tree.GetLeafId(Leaves[n]);
\r