+++ /dev/null
-#include "muscle.h"\r
-#include "clust.h"\r
-#include "clustset.h"\r
-#include <stdio.h>\r
-\r
-#define TRACE 0\r
-\r
-Clust::Clust()\r
- {\r
- m_Nodes = 0;\r
- m_uNodeCount = 0;\r
- m_uLeafCount = 0;\r
- m_uClusterCount = 0;\r
- m_JoinStyle = JOIN_Undefined;\r
- m_dDist = 0;\r
- m_uLeafCount = 0;\r
- m_ptrSet = 0;\r
- }\r
-\r
-Clust::~Clust()\r
- {\r
- delete[] m_Nodes;\r
- delete[] m_dDist;\r
- delete[] m_ClusterIndexToNodeIndex;\r
- }\r
-\r
-void Clust::Create(ClustSet &Set, CLUSTER Method)\r
- {\r
- m_ptrSet = &Set;\r
-\r
- SetLeafCount(Set.GetLeafCount());\r
-\r
- switch (Method)\r
- {\r
- case CLUSTER_UPGMA:\r
- m_JoinStyle = JOIN_NearestNeighbor;\r
- m_CentroidStyle = LINKAGE_Avg;\r
- break;\r
-\r
- case CLUSTER_UPGMAMax:\r
- m_JoinStyle = JOIN_NearestNeighbor;\r
- m_CentroidStyle = LINKAGE_Max;\r
- break;\r
-\r
- case CLUSTER_UPGMAMin:\r
- m_JoinStyle = JOIN_NearestNeighbor;\r
- m_CentroidStyle = LINKAGE_Min;\r
- break;\r
-\r
- case CLUSTER_UPGMB:\r
- m_JoinStyle = JOIN_NearestNeighbor;\r
- m_CentroidStyle = LINKAGE_Biased;\r
- break;\r
-\r
- case CLUSTER_NeighborJoining:\r
- m_JoinStyle = JOIN_NeighborJoining;\r
- m_CentroidStyle = LINKAGE_NeighborJoining;\r
- break;\r
-\r
- default:\r
- Quit("Clust::Create, invalid method %d", Method);\r
- }\r
-\r
- if (m_uLeafCount <= 1)\r
- Quit("Clust::Create: no leaves");\r
-\r
- m_uNodeCount = 2*m_uLeafCount - 1;\r
- m_Nodes = new ClustNode[m_uNodeCount];\r
- m_ClusterIndexToNodeIndex = new unsigned[m_uLeafCount];\r
-\r
- m_ptrClusterList = 0;\r
- for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
- {\r
- ClustNode &Node = m_Nodes[uNodeIndex];\r
- Node.m_uIndex = uNodeIndex;\r
- if (uNodeIndex < m_uLeafCount)\r
- {\r
- Node.m_uSize = 1;\r
- Node.m_uLeafIndexes = new unsigned[1];\r
- Node.m_uLeafIndexes[0] = uNodeIndex;\r
- AddToClusterList(uNodeIndex);\r
- }\r
- else\r
- Node.m_uSize = 0;\r
- }\r
-\r
-// Compute initial distance matrix between leaves\r
- SetProgressDesc("Build dist matrix");\r
- unsigned uPairIndex = 0;\r
- const unsigned uPairCount = (m_uLeafCount*(m_uLeafCount - 1))/2;\r
- for (unsigned i = 0; i < m_uLeafCount; ++i)\r
- for (unsigned j = 0; j < i; ++j)\r
- {\r
- const float dDist = (float) m_ptrSet->ComputeDist(*this, i, j);\r
- SetDist(i, j, dDist);\r
- if (0 == uPairIndex%10000)\r
- Progress(uPairIndex, uPairCount);\r
- ++uPairIndex;\r
- }\r
- ProgressStepsDone();\r
-\r
-// Call CreateCluster once for each internal node in the tree\r
- SetProgressDesc("Build guide tree");\r
- m_uClusterCount = m_uLeafCount;\r
- const unsigned uInternalNodeCount = m_uNodeCount - m_uLeafCount;\r
- for (unsigned uNodeIndex = m_uLeafCount; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
- {\r
- unsigned i = uNodeIndex + 1 - m_uLeafCount;\r
- Progress(i, uInternalNodeCount);\r
- CreateCluster();\r
- }\r
- ProgressStepsDone();\r
- }\r
-\r
-void Clust::CreateCluster()\r
- {\r
- unsigned uLeftNodeIndex;\r
- unsigned uRightNodeIndex;\r
- float dLeftLength;\r
- float dRightLength;\r
- ChooseJoin(&uLeftNodeIndex, &uRightNodeIndex, &dLeftLength, &dRightLength);\r
-\r
- const unsigned uNewNodeIndex = m_uNodeCount - m_uClusterCount + 1;\r
-\r
- JoinNodes(uLeftNodeIndex, uRightNodeIndex, dLeftLength, dRightLength,\r
- uNewNodeIndex);\r
-\r
-#if TRACE\r
- Log("Merge New=%u L=%u R=%u Ld=%7.2g Rd=%7.2g\n",\r
- uNewNodeIndex, uLeftNodeIndex, uRightNodeIndex, dLeftLength, dRightLength);\r
-#endif\r
-\r
-// Compute distances to other clusters\r
- --m_uClusterCount;\r
- for (unsigned uNodeIndex = GetFirstCluster(); uNodeIndex != uInsane;\r
- uNodeIndex = GetNextCluster(uNodeIndex))\r
- {\r
- if (uNodeIndex == uLeftNodeIndex || uNodeIndex == uRightNodeIndex)\r
- continue;\r
-\r
- if (uNewNodeIndex == uNodeIndex)\r
- continue;\r
-\r
- const float dDist = ComputeDist(uNewNodeIndex, uNodeIndex);\r
- SetDist(uNewNodeIndex, uNodeIndex, dDist);\r
- }\r
-\r
- for (unsigned uNodeIndex = GetFirstCluster(); uNodeIndex != uInsane;\r
- uNodeIndex = GetNextCluster(uNodeIndex))\r
- {\r
- if (uNodeIndex == uLeftNodeIndex || uNodeIndex == uRightNodeIndex)\r
- continue;\r
-\r
- if (uNewNodeIndex == uNodeIndex)\r
- continue;\r
-\r
-#if REDLACK\r
- const float dMetric = ComputeMetric(uNewNodeIndex, uNodeIndex);\r
- InsertMetric(uNewNodeIndex, uNodeIndex, dMetric);\r
-#endif\r
- }\r
- }\r
-\r
-void Clust::ChooseJoin(unsigned *ptruLeftIndex, unsigned *ptruRightIndex,\r
- float *ptrdLeftLength, float *ptrdRightLength)\r
- {\r
- switch (m_JoinStyle)\r
- {\r
- case JOIN_NearestNeighbor:\r
- ChooseJoinNearestNeighbor(ptruLeftIndex, ptruRightIndex, ptrdLeftLength,\r
- ptrdRightLength);\r
- return;\r
- case JOIN_NeighborJoining:\r
- ChooseJoinNeighborJoining(ptruLeftIndex, ptruRightIndex, ptrdLeftLength,\r
- ptrdRightLength);\r
- return;\r
- }\r
- Quit("Clust::ChooseJoin, Invalid join style %u", m_JoinStyle);\r
- }\r
-\r
-void Clust::ChooseJoinNearestNeighbor(unsigned *ptruLeftIndex,\r
- unsigned *ptruRightIndex, float *ptrdLeftLength, float *ptrdRightLength)\r
- {\r
- const unsigned uClusterCount = GetClusterCount();\r
-\r
- unsigned uMinLeftNodeIndex;\r
- unsigned uMinRightNodeIndex;\r
- GetMinMetric(&uMinLeftNodeIndex, &uMinRightNodeIndex);\r
-\r
- float dMinDist = GetDist(uMinLeftNodeIndex, uMinRightNodeIndex);\r
-\r
- const float dLeftHeight = GetHeight(uMinLeftNodeIndex);\r
- const float dRightHeight = GetHeight(uMinRightNodeIndex);\r
-\r
- *ptruLeftIndex = uMinLeftNodeIndex;\r
- *ptruRightIndex = uMinRightNodeIndex;\r
- *ptrdLeftLength = dMinDist/2 - dLeftHeight;\r
- *ptrdRightLength = dMinDist/2 - dRightHeight;\r
- }\r
-\r
-void Clust::ChooseJoinNeighborJoining(unsigned *ptruLeftIndex,\r
- unsigned *ptruRightIndex, float *ptrdLeftLength, float *ptrdRightLength)\r
- {\r
- const unsigned uClusterCount = GetClusterCount();\r
-\r
- //unsigned uMinLeftNodeIndex = uInsane;\r
- //unsigned uMinRightNodeIndex = uInsane;\r
- //float dMinD = PLUS_INFINITY;\r
- //for (unsigned i = GetFirstCluster(); i != uInsane; i = GetNextCluster(i))\r
- // {\r
- // const float ri = Calc_r(i);\r
- // for (unsigned j = GetNextCluster(i); j != uInsane; j = GetNextCluster(j))\r
- // {\r
- // const float rj = Calc_r(j);\r
- // const float dij = GetDist(i, j);\r
- // const float Dij = dij - (ri + rj);\r
- // if (Dij < dMinD)\r
- // {\r
- // dMinD = Dij;\r
- // uMinLeftNodeIndex = i;\r
- // uMinRightNodeIndex = j;\r
- // }\r
- // }\r
- // }\r
-\r
- unsigned uMinLeftNodeIndex;\r
- unsigned uMinRightNodeIndex;\r
- GetMinMetric(&uMinLeftNodeIndex, &uMinRightNodeIndex);\r
-\r
- const float dDistLR = GetDist(uMinLeftNodeIndex, uMinRightNodeIndex);\r
- const float rL = Calc_r(uMinLeftNodeIndex);\r
- const float rR = Calc_r(uMinRightNodeIndex);\r
-\r
- const float dLeftLength = (dDistLR + rL - rR)/2;\r
- const float dRightLength = (dDistLR - rL + rR)/2;\r
-\r
- *ptruLeftIndex = uMinLeftNodeIndex;\r
- *ptruRightIndex = uMinRightNodeIndex;\r
- *ptrdLeftLength = dLeftLength;\r
- *ptrdRightLength = dRightLength;\r
- }\r
-\r
-void Clust::JoinNodes(unsigned uLeftIndex, unsigned uRightIndex, float dLeftLength,\r
- float dRightLength, unsigned uNodeIndex)\r
- {\r
- ClustNode &Parent = m_Nodes[uNodeIndex];\r
- ClustNode &Left = m_Nodes[uLeftIndex];\r
- ClustNode &Right = m_Nodes[uRightIndex];\r
-\r
- Left.m_dLength = dLeftLength;\r
- Right.m_dLength = dRightLength;\r
-\r
- Parent.m_ptrLeft = &Left;\r
- Parent.m_ptrRight = &Right;\r
-\r
- Left.m_ptrParent = &Parent;\r
- Right.m_ptrParent = &Parent;\r
-\r
- const unsigned uLeftSize = Left.m_uSize;\r
- const unsigned uRightSize = Right.m_uSize;\r
- const unsigned uParentSize = uLeftSize + uRightSize;\r
- Parent.m_uSize = uParentSize;\r
-\r
- assert(0 == Parent.m_uLeafIndexes);\r
- Parent.m_uLeafIndexes = new unsigned[uParentSize];\r
-\r
- const unsigned uLeftBytes = uLeftSize*sizeof(unsigned);\r
- const unsigned uRightBytes = uRightSize*sizeof(unsigned);\r
- memcpy(Parent.m_uLeafIndexes, Left.m_uLeafIndexes, uLeftBytes);\r
- memcpy(Parent.m_uLeafIndexes + uLeftSize, Right.m_uLeafIndexes, uRightBytes);\r
-\r
- DeleteFromClusterList(uLeftIndex);\r
- DeleteFromClusterList(uRightIndex);\r
- AddToClusterList(uNodeIndex);\r
- }\r
-\r
-float Clust::Calc_r(unsigned uNodeIndex) const\r
- {\r
- const unsigned uClusterCount = GetClusterCount();\r
- if (2 == uClusterCount)\r
- return 0;\r
-\r
- float dSum = 0;\r
- for (unsigned i = GetFirstCluster(); i != uInsane; i = GetNextCluster(i))\r
- {\r
- if (i == uNodeIndex)\r
- continue;\r
- dSum += GetDist(uNodeIndex, i);\r
- }\r
- return dSum/(uClusterCount - 2);\r
- }\r
-\r
-float Clust::ComputeDist(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- switch (m_CentroidStyle)\r
- {\r
- case LINKAGE_Avg:\r
- return ComputeDistAverageLinkage(uNewNodeIndex, uNodeIndex);\r
-\r
- case LINKAGE_Min:\r
- return ComputeDistMinLinkage(uNewNodeIndex, uNodeIndex);\r
-\r
- case LINKAGE_Max:\r
- return ComputeDistMaxLinkage(uNewNodeIndex, uNodeIndex);\r
-\r
- case LINKAGE_Biased:\r
- return ComputeDistMAFFT(uNewNodeIndex, uNodeIndex);\r
-\r
- case LINKAGE_NeighborJoining:\r
- return ComputeDistNeighborJoining(uNewNodeIndex, uNodeIndex);\r
- }\r
- Quit("Clust::ComputeDist, invalid centroid style %u", m_CentroidStyle);\r
- return (float) g_dNAN;\r
- }\r
-\r
-float Clust::ComputeDistMinLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
- const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
- const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
- const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
- return (dDistL < dDistR ? dDistL : dDistR);\r
- }\r
-\r
-float Clust::ComputeDistMaxLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
- const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
- const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
- const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
- return (dDistL > dDistR ? dDistL : dDistR);\r
- }\r
-\r
-float Clust::ComputeDistAverageLinkage(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
- const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
- const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
- const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
- return (dDistL + dDistR)/2;\r
- }\r
-\r
-float Clust::ComputeDistNeighborJoining(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
- const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
- const float dDistLR = GetDist(uLeftNodeIndex, uRightNodeIndex);\r
- const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
- const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
- const float dDist = (dDistL + dDistR - dDistLR)/2;\r
- return dDist;\r
- }\r
-\r
-// This is a mysterious variant of UPGMA reverse-engineered from MAFFT source.\r
-float Clust::ComputeDistMAFFT(unsigned uNewNodeIndex, unsigned uNodeIndex)\r
- {\r
- const unsigned uLeftNodeIndex = GetLeftIndex(uNewNodeIndex);\r
- const unsigned uRightNodeIndex = GetRightIndex(uNewNodeIndex);\r
-\r
- const float dDistLR = GetDist(uLeftNodeIndex, uRightNodeIndex);\r
- const float dDistL = GetDist(uLeftNodeIndex, uNodeIndex);\r
- const float dDistR = GetDist(uRightNodeIndex, uNodeIndex);\r
- const float dMinDistLR = (dDistL < dDistR ? dDistL : dDistR);\r
- const float dSumDistLR = dDistL + dDistR;\r
- const float dDist = dMinDistLR*(1 - g_dSUEFF) + dSumDistLR*g_dSUEFF/2;\r
- return dDist;\r
- }\r
-\r
-unsigned Clust::GetClusterCount() const\r
- {\r
- return m_uClusterCount;\r
- }\r
-\r
-void Clust::LogMe() const\r
- {\r
- Log("Clust %u leaves, %u nodes, %u clusters.\n",\r
- m_uLeafCount, m_uNodeCount, m_uClusterCount);\r
-\r
- Log("Distance matrix\n");\r
- const unsigned uNodeCount = GetNodeCount();\r
- Log(" ");\r
- for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
- Log(" %7u", i);\r
- Log("\n");\r
-\r
- Log(" ");\r
- for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
- Log(" ------");\r
- Log("\n");\r
-\r
- for (unsigned i = 0; i < uNodeCount - 1; ++i)\r
- {\r
- Log("%4u: ", i);\r
- for (unsigned j = 0; j < i; ++j)\r
- Log(" %7.2g", GetDist(i, j));\r
- Log("\n");\r
- }\r
-\r
- Log("\n");\r
- Log("Node Size Prnt Left Rght Length Name\n");\r
- Log("---- ---- ---- ---- ---- ------ ----\n");\r
- for (unsigned uNodeIndex = 0; uNodeIndex < m_uNodeCount; ++uNodeIndex)\r
- {\r
- const ClustNode &Node = m_Nodes[uNodeIndex];\r
- Log("%4u %4u", uNodeIndex, Node.m_uSize);\r
- if (0 != Node.m_ptrParent)\r
- Log(" %4u", Node.m_ptrParent->m_uIndex);\r
- else\r
- Log(" ");\r
-\r
- if (0 != Node.m_ptrLeft)\r
- Log(" %4u", Node.m_ptrLeft->m_uIndex);\r
- else\r
- Log(" ");\r
-\r
- if (0 != Node.m_ptrRight)\r
- Log(" %4u", Node.m_ptrRight->m_uIndex);\r
- else\r
- Log(" ");\r
-\r
- if (uNodeIndex != m_uNodeCount - 1)\r
- Log(" %7.3g", Node.m_dLength);\r
- if (IsLeaf(uNodeIndex))\r
- {\r
- const char *ptrName = GetNodeName(uNodeIndex);\r
- if (0 != ptrName)\r
- Log(" %s", ptrName);\r
- }\r
- if (GetRootNodeIndex() == uNodeIndex)\r
- Log(" [ROOT]");\r
- Log("\n");\r
- }\r
- }\r
-\r
-const ClustNode &Clust::GetNode(unsigned uNodeIndex) const\r
- {\r
- if (uNodeIndex >= m_uNodeCount)\r
- Quit("ClustNode::GetNode(%u) %u", uNodeIndex, m_uNodeCount);\r
- return m_Nodes[uNodeIndex];\r
- }\r
-\r
-bool Clust::IsLeaf(unsigned uNodeIndex) const\r
- {\r
- return uNodeIndex < m_uLeafCount;\r
- }\r
-\r
-unsigned Clust::GetClusterSize(unsigned uNodeIndex) const\r
- {\r
- const ClustNode &Node = GetNode(uNodeIndex);\r
- return Node.m_uSize;\r
- }\r
-\r
-unsigned Clust::GetLeftIndex(unsigned uNodeIndex) const\r
- {\r
- const ClustNode &Node = GetNode(uNodeIndex);\r
- if (0 == Node.m_ptrLeft)\r
- Quit("Clust::GetLeftIndex: leaf");\r
- return Node.m_ptrLeft->m_uIndex;\r
- }\r
-\r
-unsigned Clust::GetRightIndex(unsigned uNodeIndex) const\r
- {\r
- const ClustNode &Node = GetNode(uNodeIndex);\r
- if (0 == Node.m_ptrRight)\r
- Quit("Clust::GetRightIndex: leaf");\r
- return Node.m_ptrRight->m_uIndex;\r
- }\r
-\r
-float Clust::GetLength(unsigned uNodeIndex) const\r
- {\r
- const ClustNode &Node = GetNode(uNodeIndex);\r
- return Node.m_dLength;\r
- }\r
-\r
-void Clust::SetLeafCount(unsigned uLeafCount)\r
- {\r
- if (uLeafCount <= 1)\r
- Quit("Clust::SetLeafCount(%u)", uLeafCount);\r
-\r
- m_uLeafCount = uLeafCount;\r
- const unsigned uNodeCount = GetNodeCount();\r
-\r
-// Triangular matrix size excluding diagonal (all zeros in our case).\r
- m_uTriangularMatrixSize = (uNodeCount*(uNodeCount - 1))/2;\r
- m_dDist = new float[m_uTriangularMatrixSize];\r
- }\r
-\r
-unsigned Clust::GetLeafCount() const\r
- {\r
- return m_uLeafCount;\r
- }\r
-\r
-unsigned Clust::VectorIndex(unsigned uIndex1, unsigned uIndex2) const\r
- {\r
- const unsigned uNodeCount = GetNodeCount();\r
- if (uIndex1 >= uNodeCount || uIndex2 >= uNodeCount)\r
- Quit("DistVectorIndex(%u,%u) %u", uIndex1, uIndex2, uNodeCount);\r
- unsigned v;\r
- if (uIndex1 >= uIndex2)\r
- v = uIndex2 + (uIndex1*(uIndex1 - 1))/2;\r
- else\r
- v = uIndex1 + (uIndex2*(uIndex2 - 1))/2;\r
- assert(v < m_uTriangularMatrixSize);\r
- return v;\r
- }\r
-\r
-float Clust::GetDist(unsigned uIndex1, unsigned uIndex2) const\r
- {\r
- unsigned v = VectorIndex(uIndex1, uIndex2);\r
- return m_dDist[v];\r
- }\r
-\r
-void Clust::SetDist(unsigned uIndex1, unsigned uIndex2, float dDist)\r
- {\r
- unsigned v = VectorIndex(uIndex1, uIndex2);\r
- m_dDist[v] = dDist;\r
- }\r
-\r
-float Clust::GetHeight(unsigned uNodeIndex) const\r
- {\r
- if (IsLeaf(uNodeIndex))\r
- return 0;\r
-\r
- const unsigned uLeftIndex = GetLeftIndex(uNodeIndex);\r
- const unsigned uRightIndex = GetRightIndex(uNodeIndex);\r
- const float dLeftLength = GetLength(uLeftIndex);\r
- const float dRightLength = GetLength(uRightIndex);\r
- const float dLeftHeight = dLeftLength + GetHeight(uLeftIndex);\r
- const float dRightHeight = dRightLength + GetHeight(uRightIndex);\r
- return (dLeftHeight + dRightHeight)/2;\r
- }\r
-\r
-const char *Clust::GetNodeName(unsigned uNodeIndex) const\r
- {\r
- if (!IsLeaf(uNodeIndex))\r
- Quit("Clust::GetNodeName, is not leaf");\r
- return m_ptrSet->GetLeafName(uNodeIndex);\r
- }\r
-\r
-unsigned Clust::GetNodeId(unsigned uNodeIndex) const\r
- {\r
- if (uNodeIndex >= GetLeafCount())\r
- return 0;\r
- return m_ptrSet->GetLeafId(uNodeIndex);\r
- }\r
-\r
-unsigned Clust::GetLeaf(unsigned uNodeIndex, unsigned uLeafIndex) const\r
- {\r
- const ClustNode &Node = GetNode(uNodeIndex);\r
- const unsigned uLeafCount = Node.m_uSize;\r
- if (uLeafIndex >= uLeafCount)\r
- Quit("Clust::GetLeaf, invalid index");\r
- const unsigned uIndex = Node.m_uLeafIndexes[uLeafIndex];\r
- if (uIndex >= m_uNodeCount)\r
- Quit("Clust::GetLeaf, index out of range");\r
- return uIndex;\r
- }\r
-\r
-unsigned Clust::GetFirstCluster() const\r
- {\r
- if (0 == m_ptrClusterList)\r
- return uInsane;\r
- return m_ptrClusterList->m_uIndex;\r
- }\r
-\r
-unsigned Clust::GetNextCluster(unsigned uIndex) const\r
- {\r
- ClustNode *ptrNode = &m_Nodes[uIndex];\r
- if (0 == ptrNode->m_ptrNextCluster)\r
- return uInsane;\r
- return ptrNode->m_ptrNextCluster->m_uIndex;\r
- }\r
-\r
-void Clust::DeleteFromClusterList(unsigned uNodeIndex)\r
- {\r
- assert(uNodeIndex < m_uNodeCount);\r
- ClustNode *ptrNode = &m_Nodes[uNodeIndex];\r
- ClustNode *ptrPrev = ptrNode->m_ptrPrevCluster;\r
- ClustNode *ptrNext = ptrNode->m_ptrNextCluster;\r
-\r
- if (0 != ptrNext)\r
- ptrNext->m_ptrPrevCluster = ptrPrev;\r
- if (0 == ptrPrev)\r
- {\r
- assert(m_ptrClusterList == ptrNode);\r
- m_ptrClusterList = ptrNext;\r
- }\r
- else\r
- ptrPrev->m_ptrNextCluster = ptrNext;\r
-\r
- ptrNode->m_ptrNextCluster = 0;\r
- ptrNode->m_ptrPrevCluster = 0;\r
- }\r
-\r
-void Clust::AddToClusterList(unsigned uNodeIndex)\r
- {\r
- assert(uNodeIndex < m_uNodeCount);\r
- ClustNode *ptrNode = &m_Nodes[uNodeIndex];\r
-\r
- if (0 != m_ptrClusterList)\r
- m_ptrClusterList->m_ptrPrevCluster = ptrNode;\r
-\r
- ptrNode->m_ptrNextCluster = m_ptrClusterList;\r
- ptrNode->m_ptrPrevCluster = 0;\r
-\r
- m_ptrClusterList = ptrNode;\r
- }\r
-\r
-float Clust::ComputeMetric(unsigned uIndex1, unsigned uIndex2) const\r
- {\r
- switch (m_JoinStyle)\r
- {\r
- case JOIN_NearestNeighbor:\r
- return ComputeMetricNearestNeighbor(uIndex1, uIndex2);\r
-\r
- case JOIN_NeighborJoining:\r
- return ComputeMetricNeighborJoining(uIndex1, uIndex2);\r
- }\r
- Quit("Clust::ComputeMetric");\r
- return 0;\r
- }\r
-\r
-float Clust::ComputeMetricNeighborJoining(unsigned i, unsigned j) const\r
- {\r
- float ri = Calc_r(i);\r
- float rj = Calc_r(j);\r
- float dij = GetDist(i, j);\r
- float dMetric = dij - (ri + rj);\r
- return (float) dMetric;\r
- }\r
-\r
-float Clust::ComputeMetricNearestNeighbor(unsigned i, unsigned j) const\r
- {\r
- return (float) GetDist(i, j);\r
- }\r
-\r
-float Clust::GetMinMetricBruteForce(unsigned *ptruIndex1, unsigned *ptruIndex2) const\r
- {\r
- unsigned uMinLeftNodeIndex = uInsane;\r
- unsigned uMinRightNodeIndex = uInsane;\r
- float dMinMetric = PLUS_INFINITY;\r
- for (unsigned uLeftNodeIndex = GetFirstCluster(); uLeftNodeIndex != uInsane;\r
- uLeftNodeIndex = GetNextCluster(uLeftNodeIndex))\r
- {\r
- for (unsigned uRightNodeIndex = GetNextCluster(uLeftNodeIndex);\r
- uRightNodeIndex != uInsane;\r
- uRightNodeIndex = GetNextCluster(uRightNodeIndex))\r
- {\r
- float dMetric = ComputeMetric(uLeftNodeIndex, uRightNodeIndex);\r
- if (dMetric < dMinMetric)\r
- {\r
- dMinMetric = dMetric;\r
- uMinLeftNodeIndex = uLeftNodeIndex;\r
- uMinRightNodeIndex = uRightNodeIndex;\r
- }\r
- }\r
- }\r
- *ptruIndex1 = uMinLeftNodeIndex;\r
- *ptruIndex2 = uMinRightNodeIndex;\r
- return dMinMetric;\r
- }\r
-\r
-float Clust::GetMinMetric(unsigned *ptruIndex1, unsigned *ptruIndex2) const\r
- {\r
- return GetMinMetricBruteForce(ptruIndex1, ptruIndex2);\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff