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

   1 #include "muscle.h"\r
   2 #include "clust.h"\r
   3 #include "clustset.h"\r
   4 #include <stdio.h>\r
   5 \r
   6 #define TRACE           0\r
   7 \r
   8 Clust::Clust()\r
   9         {\r
  10         m_Nodes = 0;\r
  11         m_uNodeCount = 0;\r
  12         m_uLeafCount = 0;\r
  13         m_uClusterCount = 0;\r
  14         m_JoinStyle = JOIN_Undefined;\r
  15         m_dDist = 0;\r
  16         m_uLeafCount = 0;\r
  17         m_ptrSet = 0;\r
  18         }\r
  19 \r
  20 Clust::~Clust()\r
  21         {\r
  22         delete[] m_Nodes;\r
  23         delete[] m_dDist;\r
  24         delete[] m_ClusterIndexToNodeIndex;\r
  25         }\r
  26 \r
  27 void Clust::Create(ClustSet &Set, CLUSTER Method)\r
  28         {\r
  29         m_ptrSet = &Set;\r
  30 \r
  31         SetLeafCount(Set.GetLeafCount());\r
  32 \r
  33         switch (Method)\r
  34                 {\r
  35         case CLUSTER_UPGMA:\r
  36                 m_JoinStyle = JOIN_NearestNeighbor;\r
  37                 m_CentroidStyle = LINKAGE_Avg;\r
  38                 break;\r
  39 \r
  40         case CLUSTER_UPGMAMax:\r
  41                 m_JoinStyle = JOIN_NearestNeighbor;\r
  42                 m_CentroidStyle = LINKAGE_Max;\r
  43                 break;\r
  44 \r
  45         case CLUSTER_UPGMAMin:\r
  46                 m_JoinStyle = JOIN_NearestNeighbor;\r
  47                 m_CentroidStyle = LINKAGE_Min;\r
  48                 break;\r
  49 \r
  50         case CLUSTER_UPGMB:\r
  51                 m_JoinStyle = JOIN_NearestNeighbor;\r
  52                 m_CentroidStyle = LINKAGE_Biased;\r
  53                 break;\r
  54 \r
  55         case CLUSTER_NeighborJoining:\r
  56                 m_JoinStyle = JOIN_NeighborJoining;\r
  57                 m_CentroidStyle = LINKAGE_NeighborJoining;\r
  58                 break;\r
  59 \r
  60         default:\r
  61                 Quit("Clust::Create, invalid method %d", Method);\r
  62                 }\r
  63 \r
  64         if (m_uLeafCount <= 1)\r
  65                 Quit("Clust::Create: no leaves");\r
  66 \r
  67         m_uNodeCount = 2*m_uLeafCount - 1;\r
  68         m_Nodes = new ClustNode[m_uNodeCount];\r
  69         m_ClusterIndexToNodeIndex = new unsigned[m_uLeafCount];\r
  70 \r
  71         m_ptrClusterList = 0;\r
  72         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
  73                 {\r
  74                 ClustNode &Node = m_Nodes[uNodeIndex];\r
  75                 Node.m_uIndex = uNodeIndex;\r
  76                 if (uNodeIndex < m_uLeafCount)\r
  77                         {\r
  78                         Node.m_uSize = 1;\r
  79                         Node.m_uLeafIndexes = new unsigned[1];\r
  80                         Node.m_uLeafIndexes[0] = uNodeIndex;\r
  81                         AddToClusterList(uNodeIndex);\r
  82                         }\r
  83                 else\r
  84                         Node.m_uSize = 0;\r
  85                 }\r
  86 \r
  87 // Compute initial distance matrix between leaves\r
  88         SetProgressDesc("Build dist matrix");\r
  89         unsigned uPairIndex = 0;\r
  90         const unsigned uPairCount = (m_uLeafCount*(m_uLeafCount - 1))/2;\r
  91         for (unsigned i = 0; i < m_uLeafCount; ++i)\r
  92                 for (unsigned j = 0; j < i; ++j)\r
  93                         {\r
  94                         const float dDist = (float) m_ptrSet->ComputeDist(*this, i, j);\r
  95                         SetDist(i, j, dDist);\r
  96                         if (0 == uPairIndex%10000)\r
  97                                 Progress(uPairIndex, uPairCount);\r
  98                         ++uPairIndex;\r
  99                         }\r
 100         ProgressStepsDone();\r
 101 \r
 102 // Call CreateCluster once for each internal node in the tree\r
 103         SetProgressDesc("Build guide tree");\r
 104         m_uClusterCount = m_uLeafCount;\r
 105         const unsigned uInternalNodeCount = m_uNodeCount - m_uLeafCount;\r
 106         for (unsigned uNodeIndex = m_uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
 107                 {\r
 108                 unsigned i = uNodeIndex + 1 - m_uLeafCount;\r
 109                 Progress(i, uInternalNodeCount);\r
 110                 CreateCluster();\r
 111                 }\r
 112         ProgressStepsDone();\r
 113         }\r
 114 \r
 115 void Clust::CreateCluster()\r
 116         {\r
 117         unsigned uLeftNodeIndex;\r
 118         unsigned uRightNodeIndex;\r
 119         float dLeftLength;\r
 120         float dRightLength;\r
 121         ChooseJoin(&uLeftNodeIndex, &uRightNodeIndex, &dLeftLength, &dRightLength);\r
 122 \r
 123         const unsigned uNewNodeIndex = m_uNodeCount - m_uClusterCount + 1;\r
 124 \r
 125         JoinNodes(uLeftNodeIndex, uRightNodeIndex, dLeftLength, dRightLength,\r
 126           uNewNodeIndex);\r
 127 \r
 128 #if     TRACE\r
 129         Log("Merge New=%u L=%u R=%u Ld=%7.2g Rd=%7.2g\n",\r
 130           uNewNodeIndex, uLeftNodeIndex, uRightNodeIndex, dLeftLength, dRightLength);\r
 131 #endif\r
 132 \r
 133 // Compute distances to other clusters\r
 134         --m_uClusterCount;\r
 135         for (unsigned uNodeIndex = GetFirstCluster(); uNodeIndex != uInsane;\r
 136           uNodeIndex = GetNextCluster(uNodeIndex))\r
 137                 {\r
 138                 if (uNodeIndex == uLeftNodeIndex || uNodeIndex == uRightNodeIndex)\r
 139                         continue;\r
 140 \r
 141                 if (uNewNodeIndex == uNodeIndex)\r
 142                         continue;\r
 143 \r
 144                 const float dDist = ComputeDist(uNewNodeIndex, uNodeIndex);\r
 145                 SetDist(uNewNodeIndex, uNodeIndex, dDist);\r
 146                 }\r
 147 \r
 148         for (unsigned uNodeIndex = GetFirstCluster(); uNodeIndex != uInsane;\r
 149           uNodeIndex = GetNextCluster(uNodeIndex))\r
 150                 {\r
 151                 if (uNodeIndex == uLeftNodeIndex || uNodeIndex == uRightNodeIndex)\r
 152                         continue;\r
 153 \r
 154                 if (uNewNodeIndex == uNodeIndex)\r
 155                         continue;\r
 156 \r
 157 #if     REDLACK\r
 158                 const float dMetric = ComputeMetric(uNewNodeIndex, uNodeIndex);\r
 159                 InsertMetric(uNewNodeIndex, uNodeIndex, dMetric);\r
 160 #endif\r
 161                 }\r
 162         }\r
 163 \r
 164 void Clust::ChooseJoin(unsigned *ptruLeftIndex, unsigned *ptruRightIndex,\r
 165   float *ptrdLeftLength, float *ptrdRightLength)\r
 166         {\r
 167         switch (m_JoinStyle)\r
 168                 {\r
 169         case JOIN_NearestNeighbor:\r
 170                 ChooseJoinNearestNeighbor(ptruLeftIndex, ptruRightIndex, ptrdLeftLength,\r
 171                   ptrdRightLength);\r
 172                 return;\r
 173         case JOIN_NeighborJoining:\r
 174                 ChooseJoinNeighborJoining(ptruLeftIndex, ptruRightIndex, ptrdLeftLength,\r
 175                   ptrdRightLength);\r
 176                 return;\r
 177                 }\r
 178         Quit("Clust::ChooseJoin, Invalid join style %u", m_JoinStyle);\r
 179         }\r
 180 \r
 181 void Clust::ChooseJoinNearestNeighbor(unsigned *ptruLeftIndex,\r
 182   unsigned *ptruRightIndex, float *ptrdLeftLength, float *ptrdRightLength)\r
 183         {\r
 184         const unsigned uClusterCount = GetClusterCount();\r
 185 \r
 186         unsigned uMinLeftNodeIndex;\r
 187         unsigned uMinRightNodeIndex;\r
 188         GetMinMetric(&uMinLeftNodeIndex, &uMinRightNodeIndex);\r
 189 \r
 190         float dMinDist = GetDist(uMinLeftNodeIndex, uMinRightNodeIndex);\r
 191 \r
 192         const float dLeftHeight = GetHeight(uMinLeftNodeIndex);\r
 193         const float dRightHeight = GetHeight(uMinRightNodeIndex);\r
 194 \r
 195         *ptruLeftIndex = uMinLeftNodeIndex;\r
 196         *ptruRightIndex = uMinRightNodeIndex;\r
 197         *ptrdLeftLength = dMinDist/2 - dLeftHeight;\r
 198         *ptrdRightLength = dMinDist/2 - dRightHeight;\r
 199         }\r
 200 \r
 201 void Clust::ChooseJoinNeighborJoining(unsigned *ptruLeftIndex,\r
 202   unsigned *ptruRightIndex, float *ptrdLeftLength, float *ptrdRightLength)\r
 203         {\r
 204         const unsigned uClusterCount = GetClusterCount();\r
 205 \r
 206         //unsigned uMinLeftNodeIndex = uInsane;\r
 207         //unsigned uMinRightNodeIndex = uInsane;\r
 208         //float dMinD = PLUS_INFINITY;\r
 209         //for (unsigned i = GetFirstCluster(); i != uInsane; i = GetNextCluster(i))\r
 210         //      {\r
 211         //      const float ri = Calc_r(i);\r
 212         //      for (unsigned j = GetNextCluster(i); j != uInsane; j = GetNextCluster(j))\r
 213         //              {\r
 214         //              const float rj = Calc_r(j);\r
 215         //              const float dij = GetDist(i, j);\r
 216         //              const float Dij = dij - (ri + rj);\r
 217         //              if (Dij < dMinD)\r
 218         //                      {\r
 219         //                      dMinD = Dij;\r
 220         //                      uMinLeftNodeIndex = i;\r
 221         //                      uMinRightNodeIndex = j;\r
 222         //                      }\r
 223         //              }\r
 224         //      }\r
 225 \r
 226         unsigned uMinLeftNodeIndex;\r
 227         unsigned uMinRightNodeIndex;\r
 228         GetMinMetric(&uMinLeftNodeIndex, &uMinRightNodeIndex);\r
 229 \r
 230         const float dDistLR = GetDist(uMinLeftNodeIndex, uMinRightNodeIndex);\r
 231         const float rL = Calc_r(uMinLeftNodeIndex);\r
 232         const float rR = Calc_r(uMinRightNodeIndex);\r
 233 \r
 234         const float dLeftLength = (dDistLR + rL - rR)/2;\r
 235         const float dRightLength = (dDistLR - rL + rR)/2;\r
 236 \r
 237         *ptruLeftIndex = uMinLeftNodeIndex;\r
 238         *ptruRightIndex = uMinRightNodeIndex;\r
 239         *ptrdLeftLength = dLeftLength;\r
 240         *ptrdRightLength = dRightLength;\r
 241         }\r
 242 \r
 243 void Clust::JoinNodes(unsigned uLeftIndex, unsigned uRightIndex, float dLeftLength,\r
 244   float dRightLength, unsigned uNodeIndex)\r
 245         {\r
 246         ClustNode &Parent = m_Nodes[uNodeIndex];\r
 247         ClustNode &Left = m_Nodes[uLeftIndex];\r
 248         ClustNode &Right = m_Nodes[uRightIndex];\r
 249 \r
 250         Left.m_dLength = dLeftLength;\r
 251         Right.m_dLength = dRightLength;\r
 252 \r
 253         Parent.m_ptrLeft = &Left;\r
 254         Parent.m_ptrRight = &Right;\r
 255 \r
 256         Left.m_ptrParent = &Parent;\r
 257         Right.m_ptrParent = &Parent;\r
 258 \r
 259         const unsigned uLeftSize = Left.m_uSize;\r
 260         const unsigned uRightSize = Right.m_uSize;\r
 261         const unsigned uParentSize = uLeftSize + uRightSize;\r
 262         Parent.m_uSize = uParentSize;\r
 263 \r
 264         assert(0 == Parent.m_uLeafIndexes);\r
 265         Parent.m_uLeafIndexes = new unsigned[uParentSize];\r
 266 \r
 267         const unsigned uLeftBytes = uLeftSize*sizeof(unsigned);\r
 268         const unsigned uRightBytes = uRightSize*sizeof(unsigned);\r
 269         memcpy(Parent.m_uLeafIndexes, Left.m_uLeafIndexes, uLeftBytes);\r
 270         memcpy(Parent.m_uLeafIndexes + uLeftSize, Right.m_uLeafIndexes, uRightBytes);\r
 271 \r
 272         DeleteFromClusterList(uLeftIndex);\r
 273         DeleteFromClusterList(uRightIndex);\r
 274         AddToClusterList(uNodeIndex);\r
 275         }\r
 276 \r
 277 float Clust::Calc_r(unsigned uNodeIndex) const\r
 278         {\r
 279         const unsigned uClusterCount = GetClusterCount();\r
 280         if (2 == uClusterCount)\r
 281                 return 0;\r
 282 \r
 283         float dSum = 0;\r
 284         for (unsigned i = GetFirstCluster(); i != uInsane; i = GetNextCluster(i))\r
 285                 {\r
 286                 if (i == uNodeIndex)\r
 287                         continue;\r
 288                 dSum += GetDist(uNodeIndex, i);\r
 289                 }\r
 290         return dSum/(uClusterCount - 2);\r
 291         }\r
 292 \r
 293 float Clust::ComputeDist(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 294         {\r
 295         switch (m_CentroidStyle)\r
 296                 {\r
 297         case LINKAGE_Avg:\r
 298                 return ComputeDistAverageLinkage(uNewNodeIndex, uNodeIndex);\r
 299 \r
 300         case LINKAGE_Min:\r
 301                 return ComputeDistMinLinkage(uNewNodeIndex, uNodeIndex);\r
 302 \r
 303         case LINKAGE_Max:\r
 304                 return ComputeDistMaxLinkage(uNewNodeIndex, uNodeIndex);\r
 305 \r
 306         case LINKAGE_Biased:\r
 307                 return ComputeDistMAFFT(uNewNodeIndex, uNodeIndex);\r
 308 \r
 309         case LINKAGE_NeighborJoining:\r
 310                 return ComputeDistNeighborJoining(uNewNodeIndex, uNodeIndex);\r
 311                 }\r
 312         Quit("Clust::ComputeDist, invalid centroid style %u", m_CentroidStyle);\r
 313         return (float) g_dNAN;\r
 314         }\r
 315 \r
 316 float Clust::ComputeDistMinLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 317         {\r
 318         const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
 319         const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
 320         const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
 321         const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
 322         return (dDistL < dDistR ? dDistL : dDistR);\r
 323         }\r
 324 \r
 325 float Clust::ComputeDistMaxLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 326         {\r
 327         const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
 328         const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
 329         const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
 330         const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
 331         return (dDistL > dDistR ? dDistL : dDistR);\r
 332         }\r
 333 \r
 334 float Clust::ComputeDistAverageLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 335         {\r
 336         const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
 337         const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
 338         const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
 339         const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
 340         return (dDistL + dDistR)/2;\r
 341         }\r
 342 \r
 343 float Clust::ComputeDistNeighborJoining(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 344         {\r
 345         const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
 346         const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
 347         const float dDistLR = GetDist(uLeftNodeIndex, uRightNodeIndex);\r
 348         const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
 349         const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
 350         const float dDist = (dDistL + dDistR - dDistLR)/2;\r
 351         return dDist;\r
 352         }\r
 353 \r
 354 // This is a mysterious variant of UPGMA reverse-engineered from MAFFT source.\r
 355 float Clust::ComputeDistMAFFT(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
 356         {\r
 357         const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
 358         const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
 359 \r
 360         const float dDistLR = GetDist(uLeftNodeIndex, uRightNodeIndex);\r
 361         const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
 362         const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
 363         const float dMinDistLR = (dDistL < dDistR ? dDistL : dDistR);\r
 364         const float dSumDistLR = dDistL + dDistR;\r
 365         const float dDist = dMinDistLR*(1 - g_dSUEFF) + dSumDistLR*g_dSUEFF/2;\r
 366         return dDist;\r
 367         }\r
 368 \r
 369 unsigned Clust::GetClusterCount() const\r
 370         {\r
 371         return m_uClusterCount;\r
 372         }\r
 373 \r
 374 void Clust::LogMe() const\r
 375         {\r
 376         Log("Clust %u leaves, %u nodes, %u clusters.\n",\r
 377           m_uLeafCount, m_uNodeCount, m_uClusterCount);\r
 378 \r
 379         Log("Distance matrix\n");\r
 380         const unsigned uNodeCount = GetNodeCount();\r
 381         Log("       ");\r
 382         for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
 383                 Log(" %7u", i);\r
 384         Log("\n");\r
 385 \r
 386         Log("       ");\r
 387         for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
 388                 Log("  ------");\r
 389         Log("\n");\r
 390 \r
 391         for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
 392                 {\r
 393                 Log("%4u:  ", i);\r
 394                 for (unsigned j = 0; j < i; ++j)\r
 395                         Log(" %7.2g", GetDist(i, j));\r
 396                 Log("\n");\r
 397                 }\r
 398 \r
 399         Log("\n");\r
 400         Log("Node  Size  Prnt  Left  Rght   Length  Name\n");\r
 401         Log("----  ----  ----  ----  ----   ------  ----\n");\r
 402         for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
 403                 {\r
 404                 const ClustNode &Node = m_Nodes[uNodeIndex];\r
 405                 Log("%4u  %4u", uNodeIndex, Node.m_uSize);\r
 406                 if (0 != Node.m_ptrParent)\r
 407                         Log("  %4u", Node.m_ptrParent->m_uIndex);\r
 408                 else\r
 409                         Log("      ");\r
 410 \r
 411                 if (0 != Node.m_ptrLeft)\r
 412                         Log("  %4u", Node.m_ptrLeft->m_uIndex);\r
 413                 else\r
 414                         Log("      ");\r
 415 \r
 416                 if (0 != Node.m_ptrRight)\r
 417                         Log("  %4u", Node.m_ptrRight->m_uIndex);\r
 418                 else\r
 419                         Log("      ");\r
 420 \r
 421                 if (uNodeIndex != m_uNodeCount - 1)\r
 422                         Log("  %7.3g", Node.m_dLength);\r
 423                 if (IsLeaf(uNodeIndex))\r
 424                         {\r
 425                         const char *ptrName = GetNodeName(uNodeIndex);\r
 426                         if (0 != ptrName)\r
 427                                 Log("  %s", ptrName);\r
 428                         }\r
 429                 if (GetRootNodeIndex() == uNodeIndex)\r
 430                         Log("    [ROOT]");\r
 431                 Log("\n");\r
 432                 }\r
 433         }\r
 434 \r
 435 const ClustNode &Clust::GetNode(unsigned uNodeIndex) const\r
 436         {\r
 437         if (uNodeIndex >= m_uNodeCount)\r
 438                 Quit("ClustNode::GetNode(%u) %u", uNodeIndex, m_uNodeCount);\r
 439         return m_Nodes[uNodeIndex];\r
 440         }\r
 441 \r
 442 bool Clust::IsLeaf(unsigned uNodeIndex) const\r
 443         {\r
 444         return uNodeIndex < m_uLeafCount;\r
 445         }\r
 446 \r
 447 unsigned Clust::GetClusterSize(unsigned uNodeIndex) const\r
 448         {\r
 449         const ClustNode &Node = GetNode(uNodeIndex);\r
 450         return Node.m_uSize;\r
 451         }\r
 452 \r
 453 unsigned Clust::GetLeftIndex(unsigned uNodeIndex) const\r
 454         {\r
 455         const ClustNode &Node = GetNode(uNodeIndex);\r
 456         if (0 == Node.m_ptrLeft)\r
 457                 Quit("Clust::GetLeftIndex: leaf");\r
 458         return Node.m_ptrLeft->m_uIndex;\r
 459         }\r
 460 \r
 461 unsigned Clust::GetRightIndex(unsigned uNodeIndex) const\r
 462         {\r
 463         const ClustNode &Node = GetNode(uNodeIndex);\r
 464         if (0 == Node.m_ptrRight)\r
 465                 Quit("Clust::GetRightIndex: leaf");\r
 466         return Node.m_ptrRight->m_uIndex;\r
 467         }\r
 468 \r
 469 float Clust::GetLength(unsigned uNodeIndex) const\r
 470         {\r
 471         const ClustNode &Node = GetNode(uNodeIndex);\r
 472         return Node.m_dLength;\r
 473         }\r
 474 \r
 475 void Clust::SetLeafCount(unsigned uLeafCount)\r
 476         {\r
 477         if (uLeafCount <= 1)\r
 478                 Quit("Clust::SetLeafCount(%u)", uLeafCount);\r
 479 \r
 480         m_uLeafCount = uLeafCount;\r
 481         const unsigned uNodeCount = GetNodeCount();\r
 482 \r
 483 // Triangular matrix size excluding diagonal (all zeros in our case).\r
 484         m_uTriangularMatrixSize = (uNodeCount*(uNodeCount - 1))/2;\r
 485         m_dDist = new float[m_uTriangularMatrixSize];\r
 486         }\r
 487 \r
 488 unsigned Clust::GetLeafCount() const\r
 489         {\r
 490         return m_uLeafCount;\r
 491         }\r
 492 \r
 493 unsigned Clust::VectorIndex(unsigned uIndex1, unsigned uIndex2) const\r
 494         {\r
 495         const unsigned uNodeCount = GetNodeCount();\r
 496         if (uIndex1 >= uNodeCount || uIndex2 >= uNodeCount)\r
 497                 Quit("DistVectorIndex(%u,%u) %u", uIndex1, uIndex2, uNodeCount);\r
 498         unsigned v;\r
 499         if (uIndex1 >= uIndex2)\r
 500                 v = uIndex2 + (uIndex1*(uIndex1 - 1))/2;\r
 501         else\r
 502                 v = uIndex1 + (uIndex2*(uIndex2 - 1))/2;\r
 503         assert(v < m_uTriangularMatrixSize);\r
 504         return v;\r
 505         }\r
 506 \r
 507 float Clust::GetDist(unsigned uIndex1, unsigned uIndex2) const\r
 508         {\r
 509         unsigned v = VectorIndex(uIndex1, uIndex2);\r
 510         return m_dDist[v];\r
 511         }\r
 512 \r
 513 void Clust::SetDist(unsigned uIndex1, unsigned uIndex2, float dDist)\r
 514         {\r
 515         unsigned v = VectorIndex(uIndex1, uIndex2);\r
 516         m_dDist[v] = dDist;\r
 517         }\r
 518 \r
 519 float Clust::GetHeight(unsigned uNodeIndex) const\r
 520         {\r
 521         if (IsLeaf(uNodeIndex))\r
 522                 return 0;\r
 523 \r
 524         const unsigned uLeftIndex = GetLeftIndex(uNodeIndex);\r
 525         const unsigned uRightIndex = GetRightIndex(uNodeIndex);\r
 526         const float dLeftLength = GetLength(uLeftIndex);\r
 527         const float dRightLength = GetLength(uRightIndex);\r
 528         const float dLeftHeight = dLeftLength + GetHeight(uLeftIndex);\r
 529         const float dRightHeight = dRightLength + GetHeight(uRightIndex);\r
 530         return (dLeftHeight + dRightHeight)/2;\r
 531         }\r
 532 \r
 533 const char *Clust::GetNodeName(unsigned uNodeIndex) const\r
 534         {\r
 535         if (!IsLeaf(uNodeIndex))\r
 536                 Quit("Clust::GetNodeName, is not leaf");\r
 537         return m_ptrSet->GetLeafName(uNodeIndex);\r
 538         }\r
 539 \r
 540 unsigned Clust::GetNodeId(unsigned uNodeIndex) const\r
 541         {\r
 542         if (uNodeIndex >= GetLeafCount())\r
 543                 return 0;\r
 544         return m_ptrSet->GetLeafId(uNodeIndex);\r
 545         }\r
 546 \r
 547 unsigned Clust::GetLeaf(unsigned uNodeIndex, unsigned uLeafIndex) const\r
 548         {\r
 549         const ClustNode &Node = GetNode(uNodeIndex);\r
 550         const unsigned uLeafCount = Node.m_uSize;\r
 551         if (uLeafIndex >= uLeafCount)\r
 552                 Quit("Clust::GetLeaf, invalid index");\r
 553         const unsigned uIndex = Node.m_uLeafIndexes[uLeafIndex];\r
 554         if (uIndex >= m_uNodeCount)\r
 555                 Quit("Clust::GetLeaf, index out of range");\r
 556         return uIndex;\r
 557         }\r
 558 \r
 559 unsigned Clust::GetFirstCluster() const\r
 560         {\r
 561         if (0 == m_ptrClusterList)\r
 562                 return uInsane;\r
 563         return m_ptrClusterList->m_uIndex;\r
 564         }\r
 565 \r
 566 unsigned Clust::GetNextCluster(unsigned uIndex) const\r
 567         {\r
 568         ClustNode *ptrNode = &m_Nodes[uIndex];\r
 569         if (0 == ptrNode->m_ptrNextCluster)\r
 570                 return uInsane;\r
 571         return ptrNode->m_ptrNextCluster->m_uIndex;\r
 572         }\r
 573 \r
 574 void Clust::DeleteFromClusterList(unsigned uNodeIndex)\r
 575         {\r
 576         assert(uNodeIndex < m_uNodeCount);\r
 577         ClustNode *ptrNode = &m_Nodes[uNodeIndex];\r
 578         ClustNode *ptrPrev = ptrNode->m_ptrPrevCluster;\r
 579         ClustNode *ptrNext = ptrNode->m_ptrNextCluster;\r
 580 \r
 581         if (0 != ptrNext)\r
 582                 ptrNext->m_ptrPrevCluster = ptrPrev;\r
 583         if (0 == ptrPrev)\r
 584                 {\r
 585                 assert(m_ptrClusterList == ptrNode);\r
 586                 m_ptrClusterList = ptrNext;\r
 587                 }\r
 588         else\r
 589                 ptrPrev->m_ptrNextCluster = ptrNext;\r
 590 \r
 591         ptrNode->m_ptrNextCluster = 0;\r
 592         ptrNode->m_ptrPrevCluster = 0;\r
 593         }\r
 594 \r
 595 void Clust::AddToClusterList(unsigned uNodeIndex)\r
 596         {\r
 597         assert(uNodeIndex < m_uNodeCount);\r
 598         ClustNode *ptrNode = &m_Nodes[uNodeIndex];\r
 599 \r
 600         if (0 != m_ptrClusterList)\r
 601                 m_ptrClusterList->m_ptrPrevCluster = ptrNode;\r
 602 \r
 603         ptrNode->m_ptrNextCluster = m_ptrClusterList;\r
 604         ptrNode->m_ptrPrevCluster = 0;\r
 605 \r
 606         m_ptrClusterList = ptrNode;\r
 607         }\r
 608 \r
 609 float Clust::ComputeMetric(unsigned uIndex1, unsigned uIndex2) const\r
 610         {\r
 611         switch (m_JoinStyle)\r
 612                 {\r
 613         case JOIN_NearestNeighbor:\r
 614                 return ComputeMetricNearestNeighbor(uIndex1, uIndex2);\r
 615 \r
 616         case JOIN_NeighborJoining:\r
 617                 return ComputeMetricNeighborJoining(uIndex1, uIndex2);\r
 618                 }\r
 619         Quit("Clust::ComputeMetric");\r
 620         return 0;\r
 621         }\r
 622 \r
 623 float Clust::ComputeMetricNeighborJoining(unsigned i, unsigned j) const\r
 624         {\r
 625         float ri = Calc_r(i);\r
 626         float rj = Calc_r(j);\r
 627         float dij = GetDist(i, j);\r
 628         float dMetric = dij - (ri + rj);\r
 629         return (float) dMetric;\r
 630         }\r
 631 \r
 632 float Clust::ComputeMetricNearestNeighbor(unsigned i, unsigned j) const\r
 633         {\r
 634         return (float) GetDist(i, j);\r
 635         }\r
 636 \r
 637 float Clust::GetMinMetricBruteForce(unsigned *ptruIndex1, unsigned *ptruIndex2) const\r
 638         {\r
 639         unsigned uMinLeftNodeIndex = uInsane;\r
 640         unsigned uMinRightNodeIndex = uInsane;\r
 641         float dMinMetric = PLUS_INFINITY;\r
 642         for (unsigned uLeftNodeIndex = GetFirstCluster(); uLeftNodeIndex != uInsane;\r
 643           uLeftNodeIndex = GetNextCluster(uLeftNodeIndex))\r
 644                 {\r
 645                 for (unsigned uRightNodeIndex = GetNextCluster(uLeftNodeIndex);\r
 646                   uRightNodeIndex != uInsane;\r
 647                   uRightNodeIndex = GetNextCluster(uRightNodeIndex))\r
 648                         {\r
 649                         float dMetric = ComputeMetric(uLeftNodeIndex, uRightNodeIndex);\r
 650                         if (dMetric < dMinMetric)\r
 651                                 {\r
 652                                 dMinMetric = dMetric;\r
 653                                 uMinLeftNodeIndex = uLeftNodeIndex;\r
 654                                 uMinRightNodeIndex = uRightNodeIndex;\r
 655                                 }\r
 656                         }\r
 657                 }\r
 658         *ptruIndex1 = uMinLeftNodeIndex;\r
 659         *ptruIndex2 = uMinRightNodeIndex;\r
 660         return dMinMetric;\r
 661         }\r
 662 \r
 663 float Clust::GetMinMetric(unsigned *ptruIndex1, unsigned *ptruIndex2) const\r
 664         {\r
 665         return GetMinMetricBruteForce(ptruIndex1, ptruIndex2);\r
 666         }\r