binaries/src/muscle/phy.cpp

   1 #include "muscle.h"\r
   2 #include "tree.h"\r
   3 #include <math.h>\r
   4 \r
   5 #define TRACE 0\r
   6 \r
   7 /***\r
   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
  13 \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
  18 ***/\r
  19 \r
  20 void Tree::AssertAreNeighbors(unsigned uNodeIndex1, unsigned uNodeIndex2) const\r
  21         {\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
  25 \r
  26         if (m_uNeighbor1[uNodeIndex1] != uNodeIndex2 &&\r
  27           m_uNeighbor2[uNodeIndex1] != uNodeIndex2 &&\r
  28           m_uNeighbor3[uNodeIndex1] != uNodeIndex2)\r
  29                 {\r
  30                 LogMe();\r
  31                 Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);\r
  32                 }\r
  33 \r
  34         if (m_uNeighbor1[uNodeIndex2] != uNodeIndex1 &&\r
  35           m_uNeighbor2[uNodeIndex2] != uNodeIndex1 &&\r
  36           m_uNeighbor3[uNodeIndex2] != uNodeIndex1)\r
  37                 {\r
  38                 LogMe();\r
  39                 Quit("AssertAreNeighbors(%u,%u) failed", uNodeIndex1, uNodeIndex2);\r
  40                 }\r
  41 \r
  42         bool Has12 = HasEdgeLength(uNodeIndex1, uNodeIndex2);\r
  43         bool Has21 = HasEdgeLength(uNodeIndex2, uNodeIndex1);\r
  44         if (Has12 != Has21)\r
  45                 {\r
  46                 HasEdgeLength(uNodeIndex1, uNodeIndex2);\r
  47                 HasEdgeLength(uNodeIndex2, uNodeIndex1);\r
  48                 LogMe();\r
  49                 Log("HasEdgeLength(%u, %u)=%c HasEdgeLength(%u, %u)=%c\n",\r
  50                   uNodeIndex1,\r
  51                   uNodeIndex2,\r
  52                   Has12 ? 'T' : 'F',\r
  53                   uNodeIndex2,\r
  54                   uNodeIndex1,\r
  55                   Has21 ? 'T' : 'F');\r
  56 \r
  57                 Quit("Tree::AssertAreNeighbors, HasEdgeLength not symmetric");\r
  58                 }\r
  59 \r
  60         if (Has12)\r
  61                 {\r
  62                 double d12 = GetEdgeLength(uNodeIndex1, uNodeIndex2);\r
  63                 double d21 = GetEdgeLength(uNodeIndex2, uNodeIndex1);\r
  64                 if (d12 != d21)\r
  65                         {\r
  66                         LogMe();\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
  70                         }\r
  71                 }\r
  72         }\r
  73 \r
  74 void Tree::ValidateNode(unsigned uNodeIndex) const\r
  75         {\r
  76         if (uNodeIndex >= m_uNodeCount)\r
  77                 Quit("ValidateNode(%u), %u nodes", uNodeIndex, m_uNodeCount);\r
  78 \r
  79         const unsigned uNeighborCount = GetNeighborCount(uNodeIndex);\r
  80 \r
  81         if (2 == uNeighborCount)\r
  82                 {\r
  83                 if (!m_bRooted)\r
  84                         {\r
  85                         LogMe();\r
  86                         Quit("Tree::ValidateNode: Node %u has two neighbors, tree is not rooted",\r
  87                          uNodeIndex);\r
  88                         }\r
  89                 if (uNodeIndex != m_uRootNodeIndex)\r
  90                         {\r
  91                         LogMe();\r
  92                         Quit("Tree::ValidateNode: Node %u has two neighbors, but not root node=%u",\r
  93                          uNodeIndex, m_uRootNodeIndex);\r
  94                         }\r
  95                 }\r
  96 \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
 100 \r
 101         if (NULL_NEIGHBOR == n2 && NULL_NEIGHBOR != n3)\r
 102                 {\r
 103                 LogMe();\r
 104                 Quit("Tree::ValidateNode, n2=null, n3!=null", uNodeIndex);\r
 105                 }\r
 106         if (NULL_NEIGHBOR == n3 && NULL_NEIGHBOR != n2)\r
 107                 {\r
 108                 LogMe();\r
 109                 Quit("Tree::ValidateNode, n3=null, n2!=null", uNodeIndex);\r
 110                 }\r
 111 \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
 118 \r
 119         if (n1 != NULL_NEIGHBOR && (n1 == n2 || n1 == n3))\r
 120                 {\r
 121                 LogMe();\r
 122                 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);\r
 123                 }\r
 124         if (n2 != NULL_NEIGHBOR && (n2 == n1 || n2 == n3))\r
 125                 {\r
 126                 LogMe();\r
 127                 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);\r
 128                 }\r
 129         if (n3 != NULL_NEIGHBOR && (n3 == n1 || n3 == n2))\r
 130                 {\r
 131                 LogMe();\r
 132                 Quit("Tree::ValidateNode, duplicate neighbors in node %u", uNodeIndex);\r
 133                 }\r
 134 \r
 135         if (IsRooted())\r
 136                 {\r
 137                 if (NULL_NEIGHBOR == GetParent(uNodeIndex))\r
 138                         {\r
 139                         if (uNodeIndex != m_uRootNodeIndex)\r
 140                                 {\r
 141                                 LogMe();\r
 142                                 Quit("Tree::ValiateNode(%u), no parent", uNodeIndex);\r
 143                                 }\r
 144                         }\r
 145                 else if (GetLeft(GetParent(uNodeIndex)) != uNodeIndex &&\r
 146                   GetRight(GetParent(uNodeIndex)) != uNodeIndex)\r
 147                         {\r
 148                         LogMe();\r
 149                         Quit("Tree::ValidateNode(%u), parent / child mismatch", uNodeIndex);\r
 150                         }\r
 151                 }\r
 152         }\r
 153 \r
 154 void Tree::Validate() const\r
 155         {\r
 156         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
 157                 ValidateNode(uNodeIndex);\r
 158         }\r
 159 \r
 160 bool Tree::IsEdge(unsigned uNodeIndex1, unsigned uNodeIndex2) const\r
 161         {\r
 162         assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);\r
 163 \r
 164         return m_uNeighbor1[uNodeIndex1] == uNodeIndex2 ||\r
 165           m_uNeighbor2[uNodeIndex1] == uNodeIndex2 ||\r
 166           m_uNeighbor3[uNodeIndex1] == uNodeIndex2;\r
 167         }\r
 168 \r
 169 double Tree::GetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const\r
 170         {\r
 171         assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);\r
 172         if (!HasEdgeLength(uNodeIndex1, uNodeIndex2))\r
 173                 {\r
 174                 LogMe();\r
 175                 Quit("Missing edge length in tree %u-%u", uNodeIndex1, uNodeIndex2);\r
 176                 }\r
 177 \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
 184         }\r
 185 \r
 186 void Tree::ExpandCache()\r
 187         {\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
 193 \r
 194         unsigned *uNewIds = new unsigned[uNewCacheCount];\r
 195         memset(uNewIds, 0xff, uNewCacheCount*sizeof(unsigned));\r
 196 \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
 201 \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
 206 \r
 207         char **ptrNewName = new char *[uNewCacheCount];\r
 208         memset(ptrNewName, 0, uNewCacheCount*sizeof(char *));\r
 209 \r
 210         if (m_uCacheCount > 0)\r
 211                 {\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
 216 \r
 217                 memcpy(uNewIds, m_Ids, uUnsignedBytes);\r
 218 \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
 224 \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
 230 \r
 231                 const unsigned uNameBytes = m_uCacheCount*sizeof(char *);\r
 232                 memcpy(ptrNewName, m_ptrName, uNameBytes);\r
 233 \r
 234                 delete[] m_uNeighbor1;\r
 235                 delete[] m_uNeighbor2;\r
 236                 delete[] m_uNeighbor3;\r
 237 \r
 238                 delete[] m_Ids;\r
 239 \r
 240                 delete[] m_dEdgeLength1;\r
 241                 delete[] m_dEdgeLength2;\r
 242                 delete[] m_dEdgeLength3;\r
 243 \r
 244                 delete[] m_bHasEdgeLength1;\r
 245                 delete[] m_bHasEdgeLength2;\r
 246                 delete[] m_bHasEdgeLength3;\r
 247                 delete[] m_bHasHeight;\r
 248 \r
 249                 delete[] m_ptrName;\r
 250                 }\r
 251         m_uCacheCount = uNewCacheCount;\r
 252         m_uNeighbor1 = uNewNeighbor1;\r
 253         m_uNeighbor2 = uNewNeighbor2;\r
 254         m_uNeighbor3 = uNewNeighbor3;\r
 255         m_Ids = uNewIds;\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
 265         }\r
 266 \r
 267 // Creates tree with single node, no edges.\r
 268 // Root node always has index 0.\r
 269 void Tree::CreateRooted()\r
 270         {\r
 271         Clear();\r
 272         ExpandCache();\r
 273         m_uNodeCount = 1;\r
 274 \r
 275         m_uNeighbor1[0] = NULL_NEIGHBOR;\r
 276         m_uNeighbor2[0] = NULL_NEIGHBOR;\r
 277         m_uNeighbor3[0] = NULL_NEIGHBOR;\r
 278 \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
 283 \r
 284         m_uRootNodeIndex = 0;\r
 285         m_bRooted = true;\r
 286 \r
 287 #if     DEBUG\r
 288         Validate();\r
 289 #endif\r
 290         }\r
 291 \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
 295         {\r
 296         Clear();\r
 297         ExpandCache();\r
 298 \r
 299         m_uNeighbor1[0] = 1;\r
 300         m_uNeighbor2[0] = NULL_NEIGHBOR;\r
 301         m_uNeighbor3[0] = NULL_NEIGHBOR;\r
 302 \r
 303         m_uNeighbor1[1] = 0;\r
 304         m_uNeighbor2[1] = NULL_NEIGHBOR;\r
 305         m_uNeighbor3[1] = NULL_NEIGHBOR;\r
 306 \r
 307         m_dEdgeLength1[0] = dEdgeLength;\r
 308         m_dEdgeLength1[1] = dEdgeLength;\r
 309 \r
 310         m_bHasEdgeLength1[0] = true;\r
 311         m_bHasEdgeLength1[1] = true;\r
 312 \r
 313         m_bRooted = false;\r
 314 \r
 315 #if     DEBUG\r
 316         Validate();\r
 317 #endif\r
 318         }\r
 319 \r
 320 void Tree::SetLeafName(unsigned uNodeIndex, const char *ptrName)\r
 321         {\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
 326         }\r
 327 \r
 328 void Tree::SetLeafId(unsigned uNodeIndex, unsigned uId)\r
 329         {\r
 330         assert(uNodeIndex < m_uNodeCount);\r
 331         assert(IsLeaf(uNodeIndex));\r
 332         m_Ids[uNodeIndex] = uId;\r
 333         }\r
 334 \r
 335 const char *Tree::GetLeafName(unsigned uNodeIndex) const\r
 336         {\r
 337         assert(uNodeIndex < m_uNodeCount);\r
 338         assert(IsLeaf(uNodeIndex));\r
 339         return m_ptrName[uNodeIndex];\r
 340         }\r
 341 \r
 342 unsigned Tree::GetLeafId(unsigned uNodeIndex) const\r
 343         {\r
 344         assert(uNodeIndex < m_uNodeCount);\r
 345         assert(IsLeaf(uNodeIndex));\r
 346         return m_Ids[uNodeIndex];\r
 347         }\r
 348 \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
 353         {\r
 354         if (0 == m_uNodeCount)\r
 355                 Quit("Tree::AppendBranch: tree has not been created");\r
 356 \r
 357 #if     DEBUG\r
 358         assert(uExistingLeafIndex < m_uNodeCount);\r
 359         if (!IsLeaf(uExistingLeafIndex))\r
 360                 {\r
 361                 LogMe();\r
 362                 Quit("AppendBranch(%u): not leaf", uExistingLeafIndex);\r
 363                 }\r
 364 #endif\r
 365 \r
 366         if (m_uNodeCount >= m_uCacheCount - 2)\r
 367                 ExpandCache();\r
 368 \r
 369         const unsigned uNewLeaf1 = m_uNodeCount;\r
 370         const unsigned uNewLeaf2 = m_uNodeCount + 1;\r
 371 \r
 372         m_uNodeCount += 2;\r
 373 \r
 374         assert(m_uNeighbor2[uExistingLeafIndex] == NULL_NEIGHBOR);\r
 375         assert(m_uNeighbor3[uExistingLeafIndex] == NULL_NEIGHBOR);\r
 376 \r
 377         m_uNeighbor2[uExistingLeafIndex] = uNewLeaf1;\r
 378         m_uNeighbor3[uExistingLeafIndex] = uNewLeaf2;\r
 379 \r
 380         m_uNeighbor1[uNewLeaf1] = uExistingLeafIndex;\r
 381         m_uNeighbor1[uNewLeaf2] = uExistingLeafIndex;\r
 382 \r
 383         m_uNeighbor2[uNewLeaf1] = NULL_NEIGHBOR;\r
 384         m_uNeighbor2[uNewLeaf2] = NULL_NEIGHBOR;\r
 385 \r
 386         m_uNeighbor3[uNewLeaf1] = NULL_NEIGHBOR;\r
 387         m_uNeighbor3[uNewLeaf2] = NULL_NEIGHBOR;\r
 388 \r
 389         m_dEdgeLength2[uExistingLeafIndex] = 0;\r
 390         m_dEdgeLength3[uExistingLeafIndex] = 0;\r
 391 \r
 392         m_dEdgeLength1[uNewLeaf1] = 0;\r
 393         m_dEdgeLength2[uNewLeaf1] = 0;\r
 394         m_dEdgeLength3[uNewLeaf1] = 0;\r
 395 \r
 396         m_dEdgeLength1[uNewLeaf2] = 0;\r
 397         m_dEdgeLength2[uNewLeaf2] = 0;\r
 398         m_dEdgeLength3[uNewLeaf2] = 0;\r
 399 \r
 400         m_bHasEdgeLength1[uNewLeaf1] = false;\r
 401         m_bHasEdgeLength2[uNewLeaf1] = false;\r
 402         m_bHasEdgeLength3[uNewLeaf1] = false;\r
 403 \r
 404         m_bHasEdgeLength1[uNewLeaf2] = false;\r
 405         m_bHasEdgeLength2[uNewLeaf2] = false;\r
 406         m_bHasEdgeLength3[uNewLeaf2] = false;\r
 407 \r
 408         m_bHasHeight[uNewLeaf1] = false;\r
 409         m_bHasHeight[uNewLeaf2] = false;\r
 410 \r
 411         m_Ids[uNewLeaf1] = uInsane;\r
 412         m_Ids[uNewLeaf2] = uInsane;\r
 413         return uNewLeaf1;\r
 414         }\r
 415 \r
 416 void Tree::LogMe() const\r
 417         {\r
 418         Log("Tree::LogMe %u nodes, ", m_uNodeCount);\r
 419 \r
 420         if (IsRooted())\r
 421                 {\r
 422                 Log("rooted.\n");\r
 423                 Log("\n");\r
 424                 Log("Index  Parnt  LengthP  Left   LengthL  Right  LengthR     Id  Name\n");\r
 425                 Log("-----  -----  -------  ----   -------  -----  -------  -----  ----\n");\r
 426                 }\r
 427         else\r
 428                 {\r
 429                 Log("unrooted.\n");\r
 430                 Log("\n");\r
 431                 Log("Index  Nbr_1  Length1  Nbr_2  Length2  Nbr_3  Length3     Id  Name\n");\r
 432                 Log("-----  -----  -------  -----  -------  -----  -------  -----  ----\n");\r
 433                 }\r
 434 \r
 435         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
 436                 {\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
 442                         {\r
 443                         Log("%5u  ", n1);\r
 444                         if (m_bHasEdgeLength1[uNodeIndex])\r
 445                                 Log("%7.3g  ", m_dEdgeLength1[uNodeIndex]);\r
 446                         else\r
 447                                 Log("      *  ");\r
 448                         }\r
 449                 else\r
 450                         Log("                ");\r
 451 \r
 452                 if (NULL_NEIGHBOR != n2)\r
 453                         {\r
 454                         Log("%5u  ", n2);\r
 455                         if (m_bHasEdgeLength2[uNodeIndex])\r
 456                                 Log("%7.3g  ", m_dEdgeLength2[uNodeIndex]);\r
 457                         else\r
 458                                 Log("      *  ");\r
 459                         }\r
 460                 else\r
 461                         Log("                ");\r
 462 \r
 463                 if (NULL_NEIGHBOR != n3)\r
 464                         {\r
 465                         Log("%5u  ", n3);\r
 466                         if (m_bHasEdgeLength3[uNodeIndex])\r
 467                                 Log("%7.3g  ", m_dEdgeLength3[uNodeIndex]);\r
 468                         else\r
 469                                 Log("      *  ");\r
 470                         }\r
 471                 else\r
 472                         Log("                ");\r
 473 \r
 474                 if (m_Ids != 0 && IsLeaf(uNodeIndex))\r
 475                         {\r
 476                         unsigned uId = m_Ids[uNodeIndex];\r
 477                         if (uId == uInsane)\r
 478                                 Log("    *");\r
 479                         else\r
 480                                 Log("%5u", uId);\r
 481                         }\r
 482                 else\r
 483                         Log("     ");\r
 484 \r
 485                 if (m_bRooted && uNodeIndex == m_uRootNodeIndex)\r
 486                         Log("  [ROOT] ");\r
 487                 const char *ptrName = m_ptrName[uNodeIndex];\r
 488                 if (ptrName != 0)\r
 489                         Log("  %s", ptrName);\r
 490                 Log("\n");\r
 491                 }\r
 492         }\r
 493 \r
 494 void Tree::SetEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2,\r
 495   double dLength)\r
 496         {\r
 497         assert(uNodeIndex1 < m_uNodeCount && uNodeIndex2 < m_uNodeCount);\r
 498         assert(IsEdge(uNodeIndex1, uNodeIndex2));\r
 499 \r
 500         if (m_uNeighbor1[uNodeIndex1] == uNodeIndex2)\r
 501                 {\r
 502                 m_dEdgeLength1[uNodeIndex1] = dLength;\r
 503                 m_bHasEdgeLength1[uNodeIndex1] = true;\r
 504                 }\r
 505         else if (m_uNeighbor2[uNodeIndex1] == uNodeIndex2)\r
 506                 {\r
 507                 m_dEdgeLength2[uNodeIndex1] = dLength;\r
 508                 m_bHasEdgeLength2[uNodeIndex1] = true;\r
 509 \r
 510                 }\r
 511         else\r
 512                 {\r
 513                 assert(m_uNeighbor3[uNodeIndex1] == uNodeIndex2);\r
 514                 m_dEdgeLength3[uNodeIndex1] = dLength;\r
 515                 m_bHasEdgeLength3[uNodeIndex1] = true;\r
 516                 }\r
 517 \r
 518         if (m_uNeighbor1[uNodeIndex2] == uNodeIndex1)\r
 519                 {\r
 520                 m_dEdgeLength1[uNodeIndex2] = dLength;\r
 521                 m_bHasEdgeLength1[uNodeIndex2] = true;\r
 522                 }\r
 523         else if (m_uNeighbor2[uNodeIndex2] == uNodeIndex1)\r
 524                 {\r
 525                 m_dEdgeLength2[uNodeIndex2] = dLength;\r
 526                 m_bHasEdgeLength2[uNodeIndex2] = true;\r
 527                 }\r
 528         else\r
 529                 {\r
 530                 assert(m_uNeighbor3[uNodeIndex2] == uNodeIndex1);\r
 531                 m_dEdgeLength3[uNodeIndex2] = dLength;\r
 532                 m_bHasEdgeLength3[uNodeIndex2] = true;\r
 533                 }\r
 534         }\r
 535 \r
 536 unsigned Tree::UnrootFromFile()\r
 537         {\r
 538 #if     TRACE\r
 539         Log("Before unroot:\n");\r
 540         LogMe();\r
 541 #endif\r
 542 \r
 543         if (!m_bRooted)\r
 544                 Quit("Tree::Unroot, not rooted");\r
 545 \r
 546 // Convention: root node is always node zero\r
 547         assert(IsRoot(0));\r
 548         assert(NULL_NEIGHBOR == m_uNeighbor1[0]);\r
 549 \r
 550         const unsigned uThirdNode = m_uNodeCount++;\r
 551 \r
 552         m_uNeighbor1[0] = uThirdNode;\r
 553         m_uNeighbor1[uThirdNode] = 0;\r
 554 \r
 555         m_uNeighbor2[uThirdNode] = NULL_NEIGHBOR;\r
 556         m_uNeighbor3[uThirdNode] = NULL_NEIGHBOR;\r
 557 \r
 558         m_dEdgeLength1[0] = 0;\r
 559         m_dEdgeLength1[uThirdNode] = 0;\r
 560         m_bHasEdgeLength1[uThirdNode] = true;\r
 561 \r
 562         m_bRooted = false;\r
 563 \r
 564 #if     TRACE\r
 565         Log("After unroot:\n");\r
 566         LogMe();\r
 567 #endif\r
 568 \r
 569         return uThirdNode;\r
 570         }\r
 571 \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
 581         {\r
 582         assert(uNodeIndex < m_uNodeCount);\r
 583         assert(uNeighborIndex < m_uNodeCount);\r
 584         assert(IsEdge(uNodeIndex, uNeighborIndex));\r
 585 \r
 586         for (unsigned n = 0; n < 3; ++n)\r
 587                 {\r
 588                 unsigned uNeighbor = GetNeighbor(uNodeIndex, n);\r
 589                 if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)\r
 590                         return uNeighbor;\r
 591                 }\r
 592         return NULL_NEIGHBOR;\r
 593         }\r
 594 \r
 595 unsigned Tree::GetSecondNeighbor(unsigned uNodeIndex, unsigned uNeighborIndex) const\r
 596         {\r
 597         assert(uNodeIndex < m_uNodeCount);\r
 598         assert(uNeighborIndex < m_uNodeCount);\r
 599         assert(IsEdge(uNodeIndex, uNeighborIndex));\r
 600 \r
 601         bool bFoundOne = false;\r
 602         for (unsigned n = 0; n < 3; ++n)\r
 603                 {\r
 604                 unsigned uNeighbor = GetNeighbor(uNodeIndex, n);\r
 605                 if (NULL_NEIGHBOR != uNeighbor && uNeighborIndex != uNeighbor)\r
 606                         {\r
 607                         if (bFoundOne)\r
 608                                 return uNeighbor;\r
 609                         else\r
 610                                 bFoundOne = true;\r
 611                         }\r
 612                 }\r
 613         return NULL_NEIGHBOR;\r
 614         }\r
 615 \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
 621         {\r
 622         assert(!IsRooted());\r
 623 \r
 624         if (IsLeaf(uNodeIndex2))\r
 625                 {\r
 626                 *ptrdTotalDistance = GetEdgeLength(uNodeIndex1, uNodeIndex2);\r
 627                 return 1;\r
 628                 }\r
 629 \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
 634 \r
 635         double dLeftDistance;\r
 636         double dRightDistance;\r
 637 \r
 638         const unsigned uLeftCount = GetLeafCountUnrooted(uNodeIndex2, uLeft,\r
 639           &dLeftDistance);\r
 640         const unsigned uRightCount = GetLeafCountUnrooted(uNodeIndex2, uRight,\r
 641           &dRightDistance);\r
 642 \r
 643         *ptrdTotalDistance = dLeftDistance + dRightDistance;\r
 644         return uLeftCount + uRightCount;\r
 645         }\r
 646 \r
 647 void Tree::RootUnrootedTree(ROOT Method)\r
 648         {\r
 649         assert(!IsRooted());\r
 650 #if TRACE\r
 651         Log("Tree::RootUnrootedTree, before:");\r
 652         LogMe();\r
 653 #endif\r
 654 \r
 655         unsigned uNode1;\r
 656         unsigned uNode2;\r
 657         double dLength1;\r
 658         double dLength2;\r
 659         FindRoot(*this, &uNode1, &uNode2, &dLength1, &dLength2, Method);\r
 660 \r
 661         if (m_uNodeCount == m_uCacheCount)\r
 662                 ExpandCache();\r
 663         m_uRootNodeIndex = m_uNodeCount++;\r
 664 \r
 665         double dEdgeLength = GetEdgeLength(uNode1, uNode2);\r
 666 \r
 667         m_uNeighbor1[m_uRootNodeIndex] = NULL_NEIGHBOR;\r
 668         m_uNeighbor2[m_uRootNodeIndex] = uNode1;\r
 669         m_uNeighbor3[m_uRootNodeIndex] = uNode2;\r
 670 \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
 675         else\r
 676                 {\r
 677                 assert(m_uNeighbor3[uNode1] == uNode2);\r
 678                 m_uNeighbor3[uNode1] = m_uRootNodeIndex;\r
 679                 }\r
 680 \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
 685         else\r
 686                 {\r
 687                 assert(m_uNeighbor3[uNode2] == uNode1);\r
 688                 m_uNeighbor3[uNode2] = m_uRootNodeIndex;\r
 689                 }\r
 690 \r
 691         OrientParent(uNode1, m_uRootNodeIndex);\r
 692         OrientParent(uNode2, m_uRootNodeIndex);\r
 693 \r
 694         SetEdgeLength(m_uRootNodeIndex, uNode1, dLength1);\r
 695         SetEdgeLength(m_uRootNodeIndex, uNode2, dLength2);\r
 696 \r
 697         m_bHasHeight[m_uRootNodeIndex] = false;\r
 698 \r
 699         m_ptrName[m_uRootNodeIndex] = 0;\r
 700 \r
 701         m_bRooted = true;\r
 702 \r
 703 #if     TRACE\r
 704         Log("\nPhy::RootUnrootedTree, after:");\r
 705         LogMe();\r
 706 #endif\r
 707 \r
 708         Validate();\r
 709         }\r
 710 \r
 711 bool Tree::HasEdgeLength(unsigned uNodeIndex1, unsigned uNodeIndex2) const\r
 712         {\r
 713         assert(uNodeIndex1 < m_uNodeCount);\r
 714         assert(uNodeIndex2 < m_uNodeCount);\r
 715         assert(IsEdge(uNodeIndex1, uNodeIndex2));\r
 716 \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
 723         }\r
 724 \r
 725 void Tree::OrientParent(unsigned uNodeIndex, unsigned uParentNodeIndex)\r
 726         {\r
 727         if (NULL_NEIGHBOR == uNodeIndex)\r
 728                 return;\r
 729 \r
 730         if (m_uNeighbor1[uNodeIndex] == uParentNodeIndex)\r
 731                 ;\r
 732         else if (m_uNeighbor2[uNodeIndex] == uParentNodeIndex)\r
 733                 {\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
 739                 }\r
 740         else\r
 741                 {\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
 748                 }\r
 749 \r
 750         OrientParent(m_uNeighbor2[uNodeIndex], uNodeIndex);\r
 751         OrientParent(m_uNeighbor3[uNodeIndex], uNodeIndex);\r
 752         }\r
 753 \r
 754 unsigned Tree::FirstDepthFirstNode() const\r
 755         {\r
 756         assert(IsRooted());\r
 757 \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
 762         return uNodeIndex;\r
 763         }\r
 764 \r
 765 unsigned Tree::FirstDepthFirstNodeR() const\r
 766         {\r
 767         assert(IsRooted());\r
 768 \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
 773         return uNodeIndex;\r
 774         }\r
 775 \r
 776 unsigned Tree::NextDepthFirstNode(unsigned uNodeIndex) const\r
 777         {\r
 778 #if     TRACE\r
 779         Log("NextDepthFirstNode(%3u) ", uNodeIndex);\r
 780 #endif\r
 781 \r
 782         assert(IsRooted());\r
 783         assert(uNodeIndex < m_uNodeCount);\r
 784 \r
 785         if (IsRoot(uNodeIndex))\r
 786                 {\r
 787 #if     TRACE\r
 788                 Log(">> Node %u is root, end of traversal\n", uNodeIndex);\r
 789 #endif\r
 790                 return NULL_NEIGHBOR;\r
 791                 }\r
 792 \r
 793         unsigned uParent = GetParent(uNodeIndex);\r
 794         if (GetRight(uParent) == uNodeIndex)\r
 795                 {\r
 796 #if     TRACE\r
 797                 Log(">> Is right branch, return parent=%u\n", uParent);\r
 798 #endif\r
 799                 return uParent;\r
 800                 }\r
 801 \r
 802         uNodeIndex = GetRight(uParent);\r
 803 #if     TRACE\r
 804                 Log(">> Descend left from right sibling=%u ... ", uNodeIndex);\r
 805 #endif\r
 806         while (!IsLeaf(uNodeIndex))\r
 807                 uNodeIndex = GetLeft(uNodeIndex);\r
 808 \r
 809 #if     TRACE\r
 810         Log("bottom out at leaf=%u\n", uNodeIndex);\r
 811 #endif\r
 812         return uNodeIndex;\r
 813         }\r
 814 \r
 815 unsigned Tree::NextDepthFirstNodeR(unsigned uNodeIndex) const\r
 816         {\r
 817 #if     TRACE\r
 818         Log("NextDepthFirstNode(%3u) ", uNodeIndex);\r
 819 #endif\r
 820 \r
 821         assert(IsRooted());\r
 822         assert(uNodeIndex < m_uNodeCount);\r
 823 \r
 824         if (IsRoot(uNodeIndex))\r
 825                 {\r
 826 #if     TRACE\r
 827                 Log(">> Node %u is root, end of traversal\n", uNodeIndex);\r
 828 #endif\r
 829                 return NULL_NEIGHBOR;\r
 830                 }\r
 831 \r
 832         unsigned uParent = GetParent(uNodeIndex);\r
 833         if (GetLeft(uParent) == uNodeIndex)\r
 834                 {\r
 835 #if     TRACE\r
 836                 Log(">> Is left branch, return parent=%u\n", uParent);\r
 837 #endif\r
 838                 return uParent;\r
 839                 }\r
 840 \r
 841         uNodeIndex = GetLeft(uParent);\r
 842 #if     TRACE\r
 843                 Log(">> Descend right from left sibling=%u ... ", uNodeIndex);\r
 844 #endif\r
 845         while (!IsLeaf(uNodeIndex))\r
 846                 uNodeIndex = GetRight(uNodeIndex);\r
 847 \r
 848 #if     TRACE\r
 849         Log("bottom out at leaf=%u\n", uNodeIndex);\r
 850 #endif\r
 851         return uNodeIndex;\r
 852         }\r
 853 \r
 854 void Tree::UnrootByDeletingRoot()\r
 855         {\r
 856         assert(IsRooted());\r
 857         assert(m_uNodeCount >= 3);\r
 858 \r
 859         const unsigned uLeft = GetLeft(m_uRootNodeIndex);\r
 860         const unsigned uRight = GetRight(m_uRootNodeIndex);\r
 861 \r
 862         m_uNeighbor1[uLeft] = uRight;\r
 863         m_uNeighbor1[uRight] = uLeft;\r
 864 \r
 865         bool bHasEdgeLength = HasEdgeLength(m_uRootNodeIndex, uLeft) &&\r
 866           HasEdgeLength(m_uRootNodeIndex, uRight);\r
 867         if (bHasEdgeLength)\r
 868                 {\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
 873                 }\r
 874 \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
 879           uUnsBytes);\r
 880         memmove(m_uNeighbor2 + m_uRootNodeIndex, m_uNeighbor2 + m_uRootNodeIndex + 1,\r
 881           uUnsBytes);\r
 882         memmove(m_uNeighbor3 + m_uRootNodeIndex, m_uNeighbor3 + m_uRootNodeIndex + 1,\r
 883           uUnsBytes);\r
 884 \r
 885         const unsigned uDoubleBytes = uMoveCount*sizeof(double);\r
 886         memmove(m_dEdgeLength1 + m_uRootNodeIndex, m_dEdgeLength1 + m_uRootNodeIndex + 1,\r
 887           uDoubleBytes);\r
 888         memmove(m_dEdgeLength2 + m_uRootNodeIndex, m_dEdgeLength2 + m_uRootNodeIndex + 1,\r
 889           uDoubleBytes);\r
 890         memmove(m_dEdgeLength3 + m_uRootNodeIndex, m_dEdgeLength3 + m_uRootNodeIndex + 1,\r
 891           uDoubleBytes);\r
 892 \r
 893         const unsigned uBoolBytes = uMoveCount*sizeof(bool);\r
 894         memmove(m_bHasEdgeLength1 + m_uRootNodeIndex, m_bHasEdgeLength1 + m_uRootNodeIndex + 1,\r
 895           uBoolBytes);\r
 896         memmove(m_bHasEdgeLength2 + m_uRootNodeIndex, m_bHasEdgeLength2 + m_uRootNodeIndex + 1,\r
 897           uBoolBytes);\r
 898         memmove(m_bHasEdgeLength3 + m_uRootNodeIndex, m_bHasEdgeLength3 + m_uRootNodeIndex + 1,\r
 899           uBoolBytes);\r
 900 \r
 901         const unsigned uPtrBytes = uMoveCount*sizeof(char *);\r
 902         memmove(m_ptrName + m_uRootNodeIndex, m_ptrName + m_uRootNodeIndex + 1, uPtrBytes);\r
 903 \r
 904         --m_uNodeCount;\r
 905         m_bRooted = false;\r
 906 \r
 907 // Fix up table entries\r
 908         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
 909                 {\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
 914 #undef  DEC\r
 915                 }\r
 916 \r
 917         Validate();\r
 918         }\r
 919 \r
 920 unsigned Tree::GetLeafParent(unsigned uNodeIndex) const\r
 921         {\r
 922         assert(IsLeaf(uNodeIndex));\r
 923 \r
 924         if (IsRooted())\r
 925                 return GetParent(uNodeIndex);\r
 926 \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
 932         }\r
 933 \r
 934 // TODO: This is not efficient for large trees, should cache.\r
 935 double Tree::GetNodeHeight(unsigned uNodeIndex) const\r
 936         {\r
 937         if (!IsRooted())\r
 938                 Quit("Tree::GetNodeHeight: undefined unless rooted tree");\r
 939         \r
 940         if (IsLeaf(uNodeIndex))\r
 941                 return 0.0;\r
 942 \r
 943         if (m_bHasHeight[uNodeIndex])\r
 944                 return m_dHeight[uNodeIndex];\r
 945 \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
 950 \r
 951         if (dLeftLength < 0)\r
 952                 dLeftLength = 0;\r
 953         if (dRightLength < 0)\r
 954                 dRightLength = 0;\r
 955 \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
 961         return dHeight;\r
 962         }\r
 963 \r
 964 unsigned Tree::GetNeighborSubscript(unsigned uNodeIndex, unsigned uNeighborIndex) const\r
 965         {\r
 966         assert(uNodeIndex < m_uNodeCount);\r
 967         assert(uNeighborIndex < m_uNodeCount);\r
 968         if (uNeighborIndex == m_uNeighbor1[uNodeIndex])\r
 969                 return 0;\r
 970         if (uNeighborIndex == m_uNeighbor2[uNodeIndex])\r
 971                 return 1;\r
 972         if (uNeighborIndex == m_uNeighbor3[uNodeIndex])\r
 973                 return 2;\r
 974         return NULL_NEIGHBOR;\r
 975         }\r
 976 \r
 977 unsigned Tree::GetNeighbor(unsigned uNodeIndex, unsigned uNeighborSubscript) const\r
 978         {\r
 979         switch (uNeighborSubscript)\r
 980                 {\r
 981         case 0:\r
 982                 return m_uNeighbor1[uNodeIndex];\r
 983         case 1:\r
 984                 return m_uNeighbor2[uNodeIndex];\r
 985         case 2:\r
 986                 return m_uNeighbor3[uNodeIndex];\r
 987                 }\r
 988         Quit("Tree::GetNeighbor, sub=%u", uNeighborSubscript);\r
 989         return NULL_NEIGHBOR;\r
 990         }\r
 991 \r
 992 // TODO: check if this is a performance issue, could cache a lookup table\r
 993 unsigned Tree::LeafIndexToNodeIndex(unsigned uLeafIndex) const\r
 994         {\r
 995         const unsigned uNodeCount = GetNodeCount();\r
 996         unsigned uLeafCount = 0;\r
 997         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
 998                 {\r
 999                 if (IsLeaf(uNodeIndex))\r
1000                         {\r
1001                         if (uLeafCount == uLeafIndex)\r
1002                                 return uNodeIndex;\r
1003                         else\r
1004                                 ++uLeafCount;\r
1005                         }\r
1006                 }\r
1007         Quit("LeafIndexToNodeIndex: out of range");\r
1008         return 0;\r
1009         }\r
1010 \r
1011 unsigned Tree::GetLeafNodeIndex(const char *ptrName) const\r
1012         {\r
1013         const unsigned uNodeCount = GetNodeCount();\r
1014         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
1015                 {\r
1016                 if (!IsLeaf(uNodeIndex))\r
1017                         continue;\r
1018                 const char *ptrLeafName = GetLeafName(uNodeIndex);\r
1019                 if (0 == strcmp(ptrName, ptrLeafName))\r
1020                         return uNodeIndex;\r
1021                 }\r
1022         Quit("Tree::GetLeafNodeIndex, name not found");\r
1023         return 0;\r
1024         }\r
1025 \r
1026 void Tree::Copy(const Tree &tree)\r
1027         {\r
1028         const unsigned uNodeCount = tree.GetNodeCount();\r
1029         InitCache(uNodeCount);\r
1030 \r
1031         m_uNodeCount = uNodeCount;\r
1032 \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
1036 \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
1040 \r
1041         memcpy(m_Ids, tree.m_Ids, UnsignedBytes);\r
1042 \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
1046 \r
1047         memcpy(m_dHeight, tree.m_dHeight, DoubleBytes);\r
1048 \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
1052 \r
1053         memcpy(m_bHasHeight, tree.m_bHasHeight, BoolBytes);\r
1054 \r
1055         m_uRootNodeIndex = tree.m_uRootNodeIndex;\r
1056         m_bRooted = tree.m_bRooted;\r
1057 \r
1058         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
1059                 {\r
1060                 if (tree.IsLeaf(uNodeIndex))\r
1061                         {\r
1062                         const char *ptrName = tree.GetLeafName(uNodeIndex);\r
1063                         m_ptrName[uNodeIndex] = strsave(ptrName);\r
1064                         }\r
1065                 else\r
1066                         m_ptrName[uNodeIndex] = 0;\r
1067                 }\r
1068 \r
1069 #if     DEBUG\r
1070         Validate();\r
1071 #endif\r
1072         }\r
1073 \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
1087         {\r
1088         Clear();\r
1089 \r
1090         m_uNodeCount = 2*uLeafCount - 1;\r
1091         InitCache(m_uNodeCount);\r
1092 \r
1093         for (unsigned uNodeIndex = 0; uNodeIndex < uLeafCount; ++uNodeIndex)\r
1094                 {\r
1095                 m_Ids[uNodeIndex] = LeafIds[uNodeIndex];\r
1096                 m_ptrName[uNodeIndex] = strsave(LeafNames[uNodeIndex]);\r
1097                 }\r
1098 \r
1099         for (unsigned uNodeIndex = uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
1100                 {\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
1106 \r
1107                 m_uNeighbor2[uNodeIndex] = uLeft;\r
1108                 m_uNeighbor3[uNodeIndex] = uRight;\r
1109 \r
1110                 m_bHasEdgeLength2[uNodeIndex] = true;\r
1111                 m_bHasEdgeLength3[uNodeIndex] = true;\r
1112 \r
1113                 m_dEdgeLength2[uNodeIndex] = fLeft;\r
1114                 m_dEdgeLength3[uNodeIndex] = fRight;\r
1115 \r
1116                 m_uNeighbor1[uLeft] = uNodeIndex;\r
1117                 m_uNeighbor1[uRight] = uNodeIndex;\r
1118 \r
1119                 m_dEdgeLength1[uLeft] = fLeft;\r
1120                 m_dEdgeLength1[uRight] = fRight;\r
1121 \r
1122                 m_bHasEdgeLength1[uLeft] = true;\r
1123                 m_bHasEdgeLength1[uRight] = true;\r
1124                 }\r
1125 \r
1126         m_bRooted = true;\r
1127         m_uRootNodeIndex = uRoot + uLeafCount;\r
1128 \r
1129         Validate();\r
1130         }\r