--- /dev/null
+#include "muscle.h"\r
+#include "cluster.h"\r
+#include "distfunc.h"\r
+\r
+static inline float Min(float d1, float d2)\r
+ {\r
+ return d1 < d2 ? d1 : d2;\r
+ }\r
+\r
+static inline float Max(float d1, float d2)\r
+ {\r
+ return d1 > d2 ? d1 : d2;\r
+ }\r
+\r
+static inline float Mean(float d1, float d2)\r
+ {\r
+ return (float) ((d1 + d2)/2.0);\r
+ }\r
+\r
+#if _DEBUG\r
+void ClusterTree::Validate(unsigned uNodeCount)\r
+ {\r
+ unsigned n;\r
+ ClusterNode *pNode;\r
+ unsigned uDisjointListCount = 0;\r
+ for (pNode = m_ptrDisjoints; pNode; pNode = pNode->GetNextDisjoint())\r
+ {\r
+ ClusterNode *pPrev = pNode->GetPrevDisjoint();\r
+ ClusterNode *pNext = pNode->GetNextDisjoint();\r
+ if (0 != pPrev)\r
+ {\r
+ if (pPrev->GetNextDisjoint() != pNode)\r
+ {\r
+ Log("Prev->This mismatch, prev=\n");\r
+ pPrev->LogMe();\r
+ Log("This=\n");\r
+ pNode->LogMe();\r
+ Quit("ClusterTree::Validate()");\r
+ }\r
+ }\r
+ else\r
+ {\r
+ if (pNode != m_ptrDisjoints)\r
+ {\r
+ Log("[%u]->prev = 0 but != m_ptrDisjoints=%d\n",\r
+ pNode->GetIndex(),\r
+ m_ptrDisjoints ? m_ptrDisjoints->GetIndex() : 0xffffffff);\r
+ pNode->LogMe();\r
+ Quit("ClusterTree::Validate()");\r
+ }\r
+ }\r
+ if (0 != pNext)\r
+ {\r
+ if (pNext->GetPrevDisjoint() != pNode)\r
+ {\r
+ Log("Next->This mismatch, next=\n");\r
+ pNext->LogMe();\r
+ Log("This=\n");\r
+ pNode->LogMe();\r
+ Quit("ClusterTree::Validate()");\r
+ }\r
+ }\r
+ ++uDisjointListCount;\r
+ if (uDisjointListCount > m_uNodeCount)\r
+ Quit("Loop in disjoint list");\r
+ }\r
+\r
+ unsigned uParentlessNodeCount = 0;\r
+ for (n = 0; n < uNodeCount; ++n)\r
+ if (0 == m_Nodes[n].GetParent())\r
+ ++uParentlessNodeCount;\r
+ \r
+ if (uDisjointListCount != uParentlessNodeCount)\r
+ Quit("Disjoints = %u Parentless = %u\n", uDisjointListCount,\r
+ uParentlessNodeCount);\r
+ }\r
+#else // !_DEBUG\r
+#define Validate(uNodeCount) // empty\r
+#endif\r
+\r
+void ClusterNode::LogMe() const\r
+ {\r
+ unsigned uClusterSize = GetClusterSize();\r
+ Log("[%02u] w=%5.3f CW=%5.3f LBW=%5.3f RBW=%5.3f LWT=%5.3f RWT=%5.3f L=%02d R=%02d P=%02d NxDj=%02d PvDj=%02d Sz=%02d {",\r
+ m_uIndex,\r
+ m_dWeight,\r
+ GetClusterWeight(),\r
+ GetLeftBranchWeight(),\r
+ GetRightBranchWeight(),\r
+ GetLeftWeight(),\r
+ GetRightWeight(),\r
+ m_ptrLeft ? m_ptrLeft->GetIndex() : 0xffffffff,\r
+ m_ptrRight ? m_ptrRight->GetIndex() : 0xffffffff,\r
+ m_ptrParent ? m_ptrParent->GetIndex() : 0xffffffff,\r
+ m_ptrNextDisjoint ? m_ptrNextDisjoint->GetIndex() : 0xffffffff,\r
+ m_ptrPrevDisjoint ? m_ptrPrevDisjoint->GetIndex() : 0xffffffff,\r
+ uClusterSize);\r
+ for (unsigned i = 0; i < uClusterSize; ++i)\r
+ Log(" %u", GetClusterLeaf(i)->GetIndex());\r
+ Log(" }\n");\r
+ }\r
+\r
+// How many leaves in the sub-tree under this node?\r
+unsigned ClusterNode::GetClusterSize() const\r
+ {\r
+ unsigned uLeafCount = 0;\r
+\r
+ if (0 == m_ptrLeft && 0 == m_ptrRight)\r
+ return 1;\r
+\r
+ if (0 != m_ptrLeft)\r
+ uLeafCount += m_ptrLeft->GetClusterSize();\r
+ if (0 != m_ptrRight)\r
+ uLeafCount += m_ptrRight->GetClusterSize();\r
+ assert(uLeafCount > 0);\r
+ return uLeafCount;\r
+ }\r
+\r
+double ClusterNode::GetClusterWeight() const\r
+ {\r
+ double dWeight = 0.0;\r
+ if (0 != m_ptrLeft)\r
+ dWeight += m_ptrLeft->GetClusterWeight();\r
+ if (0 != m_ptrRight)\r
+ dWeight += m_ptrRight->GetClusterWeight();\r
+ return dWeight + GetWeight();\r
+ }\r
+\r
+double ClusterNode::GetLeftBranchWeight() const\r
+ {\r
+ const ClusterNode *ptrLeft = GetLeft();\r
+ if (0 == ptrLeft)\r
+ return 0.0;\r
+\r
+ return GetWeight() - ptrLeft->GetWeight();\r
+ }\r
+\r
+double ClusterNode::GetRightBranchWeight() const\r
+ {\r
+ const ClusterNode *ptrRight = GetRight();\r
+ if (0 == ptrRight)\r
+ return 0.0;\r
+\r
+ return GetWeight() - ptrRight->GetWeight();\r
+ }\r
+\r
+double ClusterNode::GetRightWeight() const\r
+ {\r
+ const ClusterNode *ptrRight = GetRight();\r
+ if (0 == ptrRight)\r
+ return 0.0;\r
+ return ptrRight->GetClusterWeight() + GetWeight();\r
+ }\r
+\r
+double ClusterNode::GetLeftWeight() const\r
+ {\r
+ const ClusterNode *ptrLeft = GetLeft();\r
+ if (0 == ptrLeft)\r
+ return 0.0;\r
+ return ptrLeft->GetClusterWeight() + GetWeight();\r
+ }\r
+\r
+// Return n'th leaf in the sub-tree under this node.\r
+const ClusterNode *ClusterNode::GetClusterLeaf(unsigned uLeafIndex) const\r
+ {\r
+ if (0 != m_ptrLeft)\r
+ {\r
+ if (0 == m_ptrRight)\r
+ return this;\r
+\r
+ unsigned uLeftLeafCount = m_ptrLeft->GetClusterSize();\r
+\r
+ if (uLeafIndex < uLeftLeafCount)\r
+ return m_ptrLeft->GetClusterLeaf(uLeafIndex);\r
+\r
+ assert(uLeafIndex >= uLeftLeafCount);\r
+ return m_ptrRight->GetClusterLeaf(uLeafIndex - uLeftLeafCount);\r
+ }\r
+ if (0 == m_ptrRight)\r
+ return this;\r
+ return m_ptrRight->GetClusterLeaf(uLeafIndex);\r
+ }\r
+\r
+void ClusterTree::DeleteFromDisjoints(ClusterNode *ptrNode)\r
+ {\r
+ ClusterNode *ptrPrev = ptrNode->GetPrevDisjoint();\r
+ ClusterNode *ptrNext = ptrNode->GetNextDisjoint();\r
+\r
+ if (0 != ptrPrev)\r
+ ptrPrev->SetNextDisjoint(ptrNext);\r
+ else\r
+ m_ptrDisjoints = ptrNext;\r
+\r
+ if (0 != ptrNext)\r
+ ptrNext->SetPrevDisjoint(ptrPrev);\r
+\r
+#if _DEBUG\r
+// not algorithmically necessary, but improves clarity\r
+// and supports Validate().\r
+ ptrNode->SetPrevDisjoint(0);\r
+ ptrNode->SetNextDisjoint(0);\r
+#endif\r
+ }\r
+\r
+void ClusterTree::AddToDisjoints(ClusterNode *ptrNode)\r
+ {\r
+ ptrNode->SetNextDisjoint(m_ptrDisjoints);\r
+ ptrNode->SetPrevDisjoint(0);\r
+ if (0 != m_ptrDisjoints)\r
+ m_ptrDisjoints->SetPrevDisjoint(ptrNode);\r
+ m_ptrDisjoints = ptrNode;\r
+ }\r
+\r
+ClusterTree::ClusterTree()\r
+ {\r
+ m_ptrDisjoints = 0;\r
+ m_Nodes = 0;\r
+ m_uNodeCount = 0;\r
+ }\r
+\r
+ClusterTree::~ClusterTree()\r
+ {\r
+ delete[] m_Nodes;\r
+ }\r
+\r
+void ClusterTree::LogMe() const\r
+ {\r
+ Log("Disjoints=%d\n", m_ptrDisjoints ? m_ptrDisjoints->GetIndex() : 0xffffffff);\r
+ for (unsigned i = 0; i < m_uNodeCount; ++i)\r
+ {\r
+ m_Nodes[i].LogMe();\r
+ }\r
+ }\r
+\r
+ClusterNode *ClusterTree::GetRoot() const\r
+ {\r
+ return &m_Nodes[m_uNodeCount - 1];\r
+ }\r
+\r
+// This is the UPGMA algorithm as described in Durbin et al. p166.\r
+void ClusterTree::Create(const DistFunc &Dist)\r
+ {\r
+ unsigned i;\r
+ m_uLeafCount = Dist.GetCount();\r
+ m_uNodeCount = 2*m_uLeafCount - 1;\r
+\r
+ delete[] m_Nodes;\r
+ m_Nodes = new ClusterNode[m_uNodeCount];\r
+\r
+ for (i = 0; i < m_uNodeCount; ++i)\r
+ m_Nodes[i].SetIndex(i);\r
+\r
+ for (i = 0; i < m_uLeafCount - 1; ++i)\r
+ m_Nodes[i].SetNextDisjoint(&m_Nodes[i+1]);\r
+\r
+ for (i = 1; i < m_uLeafCount; ++i)\r
+ m_Nodes[i].SetPrevDisjoint(&m_Nodes[i-1]);\r
+ \r
+ m_ptrDisjoints = &m_Nodes[0];\r
+\r
+// Log("Initial state\n");\r
+// LogMe();\r
+// Log("\n");\r
+\r
+ DistFunc ClusterDist;\r
+ ClusterDist.SetCount(m_uNodeCount);\r
+ double dMaxDist = 0.0;\r
+ for (i = 0; i < m_uLeafCount; ++i)\r
+ for (unsigned j = 0; j < m_uLeafCount; ++j)\r
+ {\r
+ float dDist = Dist.GetDist(i, j);\r
+ ClusterDist.SetDist(i, j, dDist);\r
+ }\r
+\r
+ Validate(m_uLeafCount);\r
+\r
+// Iteration. N-1 joins needed to create a binary tree from N leaves.\r
+ for (unsigned uJoinIndex = m_uLeafCount; uJoinIndex < m_uNodeCount;\r
+ ++uJoinIndex)\r
+ {\r
+ // Find closest pair of clusters\r
+ unsigned uIndexClosest1;\r
+ unsigned uIndexClosest2;\r
+ bool bFound = false;\r
+ double dDistClosest = 9e99;\r
+ for (ClusterNode *ptrNode1 = m_ptrDisjoints; ptrNode1;\r
+ ptrNode1 = ptrNode1->GetNextDisjoint())\r
+ {\r
+ for (ClusterNode *ptrNode2 = ptrNode1->GetNextDisjoint(); ptrNode2;\r
+ ptrNode2 = ptrNode2->GetNextDisjoint())\r
+ {\r
+ unsigned i1 = ptrNode1->GetIndex();\r
+ unsigned i2 = ptrNode2->GetIndex();\r
+ double dDist = ClusterDist.GetDist(i1, i2);\r
+ if (dDist < dDistClosest)\r
+ {\r
+ bFound = true;\r
+ dDistClosest = dDist;\r
+ uIndexClosest1 = i1;\r
+ uIndexClosest2 = i2;\r
+ }\r
+ }\r
+ }\r
+ assert(bFound);\r
+\r
+ ClusterNode &Join = m_Nodes[uJoinIndex];\r
+ ClusterNode &Child1 = m_Nodes[uIndexClosest1];\r
+ ClusterNode &Child2 = m_Nodes[uIndexClosest2];\r
+\r
+ Join.SetLeft(&Child1);\r
+ Join.SetRight(&Child2);\r
+ Join.SetWeight(dDistClosest);\r
+\r
+ Child1.SetParent(&Join);\r
+ Child2.SetParent(&Join);\r
+\r
+ DeleteFromDisjoints(&Child1);\r
+ DeleteFromDisjoints(&Child2);\r
+ AddToDisjoints(&Join);\r
+\r
+// Log("After join %d %d\n", uIndexClosest1, uIndexClosest2);\r
+// LogMe();\r
+\r
+ // Calculate distance of every remaining disjoint cluster to the\r
+ // new cluster created by the join\r
+ for (ClusterNode *ptrNode = m_ptrDisjoints; ptrNode;\r
+ ptrNode = ptrNode->GetNextDisjoint())\r
+ {\r
+ unsigned uNodeIndex = ptrNode->GetIndex();\r
+ float dDist1 = ClusterDist.GetDist(uNodeIndex, uIndexClosest1);\r
+ float dDist2 = ClusterDist.GetDist(uNodeIndex, uIndexClosest2);\r
+ float dDist = Min(dDist1, dDist2);\r
+ ClusterDist.SetDist(uJoinIndex, uNodeIndex, dDist);\r
+ }\r
+ Validate(uJoinIndex+1);\r
+ }\r
+ GetRoot()->GetClusterWeight();\r
+// LogMe();\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff