website/archive/binaries/mac/src/muscle/phy4.cpp

   1 #include "muscle.h"\r
   2 #include "tree.h"\r
   3 #include <stdio.h>\r
   4 \r
   5 #define TRACE   0\r
   6 \r
   7 void ClusterByHeight(const Tree &tree, double dMaxHeight, unsigned Subtrees[],\r
   8   unsigned *ptruSubtreeCount)\r
   9         {\r
  10         if (!tree.IsRooted())\r
  11                 Quit("ClusterByHeight: requires rooted tree");\r
  12 \r
  13 #if     TRACE\r
  14         Log("ClusterByHeight, max height=%g\n", dMaxHeight);\r
  15 #endif\r
  16 \r
  17         unsigned uSubtreeCount = 0;\r
  18         const unsigned uNodeCount = tree.GetNodeCount();\r
  19         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
  20                 {\r
  21                 if (tree.IsRoot(uNodeIndex))\r
  22                         continue;\r
  23                 unsigned uParent = tree.GetParent(uNodeIndex);\r
  24                 double dHeight = tree.GetNodeHeight(uNodeIndex);\r
  25                 double dParentHeight = tree.GetNodeHeight(uParent);\r
  26 \r
  27 #if     TRACE\r
  28                 Log("Node %3u  Height %5.2f  ParentHeight %5.2f\n",\r
  29                   uNodeIndex, dHeight, dParentHeight);\r
  30 #endif\r
  31                 if (dParentHeight > dMaxHeight && dHeight <= dMaxHeight)\r
  32                         {\r
  33                         Subtrees[uSubtreeCount] = uNodeIndex;\r
  34 #if     TRACE\r
  35                         Log("Subtree[%u]=%u\n", uSubtreeCount, uNodeIndex);\r
  36 #endif\r
  37                         ++uSubtreeCount;\r
  38                         }\r
  39                 }\r
  40         *ptruSubtreeCount = uSubtreeCount;\r
  41         }\r
  42 \r
  43 static void ClusterBySubfamCount_Iteration(const Tree &tree, unsigned Subfams[],\r
  44   unsigned uCount)\r
  45         {\r
  46 // Find highest child node of current set of subfamilies.\r
  47         double dHighestHeight = -1e20;\r
  48         int iParentSubscript = -1;\r
  49 \r
  50         for (int n = 0; n < (int) uCount; ++n)\r
  51                 {\r
  52                 const unsigned uNodeIndex = Subfams[n];\r
  53                 if (tree.IsLeaf(uNodeIndex))\r
  54                         continue;\r
  55 \r
  56                 const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
  57                 const double dHeightLeft = tree.GetNodeHeight(uLeft);\r
  58                 if (dHeightLeft > dHighestHeight)\r
  59                         {\r
  60                         dHighestHeight = dHeightLeft;\r
  61                         iParentSubscript = n;\r
  62                         }\r
  63 \r
  64                 const unsigned uRight = tree.GetRight(uNodeIndex);\r
  65                 const double dHeightRight = tree.GetNodeHeight(uRight);\r
  66                 if (dHeightRight > dHighestHeight)\r
  67                         {\r
  68                         dHighestHeight = dHeightRight;\r
  69                         iParentSubscript = n;\r
  70                         }\r
  71                 }\r
  72 \r
  73         if (-1 == iParentSubscript)\r
  74                 Quit("CBSFCIter: failed to find highest child");\r
  75 \r
  76         const unsigned uNodeIndex = Subfams[iParentSubscript];\r
  77         const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
  78         const unsigned uRight = tree.GetRight(uNodeIndex);\r
  79 \r
  80 // Delete parent by replacing with left child\r
  81         Subfams[iParentSubscript] = uLeft;\r
  82 \r
  83 // Append right child to list\r
  84         Subfams[uCount] = uRight;\r
  85 \r
  86 #if     TRACE\r
  87         {\r
  88         Log("Iter %3u:", uCount);\r
  89         for (unsigned n = 0; n < uCount; ++n)\r
  90                 Log(" %u", Subfams[n]);\r
  91         Log("\n");\r
  92         }\r
  93 #endif\r
  94         }\r
  95 \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
 113         {\r
 114         const unsigned uNodeCount = tree.GetNodeCount();\r
 115         const unsigned uLeafCount = (uNodeCount + 1)/2;\r
 116 \r
 117 // Special case: empty tree\r
 118         if (0 == uNodeCount)\r
 119                 {\r
 120                 *ptruSubfamCount = 0;\r
 121                 return;\r
 122                 }\r
 123 \r
 124 // Special case: more subfamilies than leaves\r
 125         if (uSubfamCount >= uLeafCount)\r
 126                 {\r
 127                 for (unsigned n = 0; n < uLeafCount; ++n)\r
 128                         Subfams[n] = n;\r
 129                 *ptruSubfamCount = uLeafCount;\r
 130                 return;\r
 131                 }\r
 132 \r
 133 // Initialize list of subfamilies to be root\r
 134         Subfams[0] = tree.GetRootNodeIndex();\r
 135 \r
 136 // Iterate\r
 137         for (unsigned i = 1; i < uSubfamCount; ++i)\r
 138                 ClusterBySubfamCount_Iteration(tree, Subfams, i);\r
 139         \r
 140         *ptruSubfamCount = uSubfamCount;\r
 141         }\r
 142 \r
 143 static void GetLeavesRecurse(const Tree &tree, unsigned uNodeIndex,\r
 144   unsigned Leaves[], unsigned &uLeafCount /* in-out */)\r
 145         {\r
 146         if (tree.IsLeaf(uNodeIndex))\r
 147                 {\r
 148                 Leaves[uLeafCount] = uNodeIndex;\r
 149                 ++uLeafCount;\r
 150                 return;\r
 151                 }\r
 152 \r
 153         const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
 154         const unsigned uRight = tree.GetRight(uNodeIndex);\r
 155 \r
 156         GetLeavesRecurse(tree, uLeft, Leaves, uLeafCount);\r
 157         GetLeavesRecurse(tree, uRight, Leaves, uLeafCount);\r
 158         }\r
 159 \r
 160 void GetLeaves(const Tree &tree, unsigned uNodeIndex, unsigned Leaves[],\r
 161   unsigned *ptruLeafCount)\r
 162         {\r
 163         unsigned uLeafCount = 0;\r
 164         GetLeavesRecurse(tree, uNodeIndex, Leaves, uLeafCount);\r
 165         *ptruLeafCount = uLeafCount;\r
 166         }\r
 167 \r
 168 void Tree::PruneTree(const Tree &tree, unsigned Subfams[],\r
 169   unsigned uSubfamCount)\r
 170         {\r
 171         if (!tree.IsRooted())\r
 172                 Quit("Tree::PruneTree: requires rooted tree");\r
 173 \r
 174         Clear();\r
 175 \r
 176         m_uNodeCount = 2*uSubfamCount - 1;\r
 177         InitCache(m_uNodeCount);\r
 178 \r
 179         const unsigned uUnprunedNodeCount = tree.GetNodeCount();\r
 180 \r
 181         unsigned *uUnprunedToPrunedIndex = new unsigned[uUnprunedNodeCount];\r
 182         unsigned *uPrunedToUnprunedIndex = new unsigned[m_uNodeCount];\r
 183 \r
 184         for (unsigned n = 0; n < uUnprunedNodeCount; ++n)\r
 185                 uUnprunedToPrunedIndex[n] = NULL_NEIGHBOR;\r
 186 \r
 187         for (unsigned n = 0; n < m_uNodeCount; ++n)\r
 188                 uPrunedToUnprunedIndex[n] = NULL_NEIGHBOR;\r
 189 \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
 193                 {\r
 194                 unsigned uUnprunedNodeIndex = Subfams[uSubfamIndex];\r
 195                 uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uSubfamIndex;\r
 196                 uPrunedToUnprunedIndex[uSubfamIndex] = uUnprunedNodeIndex;\r
 197                 for (;;)\r
 198                         {\r
 199                         uUnprunedNodeIndex = tree.GetParent(uUnprunedNodeIndex);\r
 200                         if (tree.IsRoot(uUnprunedNodeIndex))\r
 201                                 break;\r
 202 \r
 203                 // Already visited this node?\r
 204                         if (NULL_NEIGHBOR != uUnprunedToPrunedIndex[uUnprunedNodeIndex])\r
 205                                 break;\r
 206 \r
 207                         uUnprunedToPrunedIndex[uUnprunedNodeIndex] = uInternalNodeIndex;\r
 208                         uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedNodeIndex;\r
 209 \r
 210                         ++uInternalNodeIndex;\r
 211                         }\r
 212                 }\r
 213 \r
 214         const unsigned uUnprunedRootIndex = tree.GetRootNodeIndex();\r
 215         uUnprunedToPrunedIndex[uUnprunedRootIndex] = uInternalNodeIndex;\r
 216         uPrunedToUnprunedIndex[uInternalNodeIndex] = uUnprunedRootIndex;\r
 217 \r
 218 #if     TRACE\r
 219         {\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
 223         Log("\n");\r
 224         Log("Unpruned to pruned:\n");\r
 225         for (unsigned i = 0; i < uUnprunedNodeCount; ++i)\r
 226                 {\r
 227                 unsigned n = uUnprunedToPrunedIndex[i];\r
 228                 if (n != NULL_NEIGHBOR)\r
 229                         Log(" [%u]=%u", i, n);\r
 230                 }\r
 231         Log("\n");\r
 232         }\r
 233 #endif\r
 234 \r
 235         if (uInternalNodeIndex != m_uNodeCount - 1)\r
 236                 Quit("Tree::PruneTree, Internal error");\r
 237 \r
 238 // Nodes 0, 1 ... are the leaves\r
 239         for (unsigned uSubfamIndex = 0; uSubfamIndex < uSubfamCount; ++uSubfamIndex)\r
 240                 {\r
 241                 char szName[32];\r
 242                 sprintf(szName, "Subfam_%u", uSubfamIndex + 1);\r
 243                 m_ptrName[uSubfamIndex] = strsave(szName);\r
 244                 }\r
 245 \r
 246         for (unsigned uPrunedNodeIndex = uSubfamCount; uPrunedNodeIndex < m_uNodeCount;\r
 247           ++uPrunedNodeIndex)\r
 248                 {\r
 249                 unsigned uUnprunedNodeIndex = uPrunedToUnprunedIndex[uPrunedNodeIndex];\r
 250 \r
 251                 const unsigned uUnprunedLeft = tree.GetLeft(uUnprunedNodeIndex);\r
 252                 const unsigned uUnprunedRight = tree.GetRight(uUnprunedNodeIndex);\r
 253 \r
 254                 const unsigned uPrunedLeft = uUnprunedToPrunedIndex[uUnprunedLeft];\r
 255                 const unsigned uPrunedRight = uUnprunedToPrunedIndex[uUnprunedRight];\r
 256 \r
 257                 const double dLeftLength =\r
 258                   tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedLeft);\r
 259                 const double dRightLength =\r
 260                   tree.GetEdgeLength(uUnprunedNodeIndex, uUnprunedRight);\r
 261 \r
 262                 m_uNeighbor2[uPrunedNodeIndex] = uPrunedLeft;\r
 263                 m_uNeighbor3[uPrunedNodeIndex] = uPrunedRight;\r
 264 \r
 265                 m_dEdgeLength1[uPrunedLeft] = dLeftLength;\r
 266                 m_dEdgeLength1[uPrunedRight] = dRightLength;\r
 267 \r
 268                 m_uNeighbor1[uPrunedLeft] = uPrunedNodeIndex;\r
 269                 m_uNeighbor1[uPrunedRight] = uPrunedNodeIndex;\r
 270 \r
 271                 m_bHasEdgeLength1[uPrunedLeft] = true;\r
 272                 m_bHasEdgeLength1[uPrunedRight] = true;\r
 273 \r
 274                 m_dEdgeLength2[uPrunedNodeIndex] = dLeftLength;\r
 275                 m_dEdgeLength3[uPrunedNodeIndex] = dRightLength;\r
 276 \r
 277                 m_bHasEdgeLength2[uPrunedNodeIndex] = true;\r
 278                 m_bHasEdgeLength3[uPrunedNodeIndex] = true;\r
 279                 }\r
 280 \r
 281         m_uRootNodeIndex = uUnprunedToPrunedIndex[uUnprunedRootIndex];\r
 282 \r
 283         m_bRooted = true;\r
 284 \r
 285         Validate();\r
 286 \r
 287         delete[] uUnprunedToPrunedIndex;\r
 288         }\r
 289 \r
 290 void LeafIndexesToIds(const Tree &tree, const unsigned Leaves[], unsigned uCount,\r
 291   unsigned Ids[])\r
 292         {\r
 293         for (unsigned n = 0; n < uCount; ++n)\r
 294                 Ids[n] = tree.GetLeafId(Leaves[n]);\r
 295         }\r