--- /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