+++ /dev/null
-#include "muscle.h"\r
-#include "tree.h"\r
-#include "distcalc.h"\r
-\r
-// UPGMA clustering in O(N^2) time and space.\r
-\r
-#define TRACE 0\r
-\r
-#define MIN(x, y) ((x) < (y) ? (x) : (y))\r
-#define MAX(x, y) ((x) > (y) ? (x) : (y))\r
-#define AVG(x, y) (((x) + (y))/2)\r
-\r
-static unsigned g_uLeafCount;\r
-static unsigned g_uTriangleSize;\r
-static unsigned g_uInternalNodeCount;\r
-static unsigned g_uInternalNodeIndex;\r
-\r
-// Triangular distance matrix is g_Dist, which is allocated\r
-// as a one-dimensional vector of length g_uTriangleSize.\r
-// TriangleSubscript(i,j) maps row,column=i,j to the subscript\r
-// into this vector.\r
-// Row / column coordinates are a bit messy.\r
-// Initially they are leaf indexes 0..N-1.\r
-// But each time we create a new node (=new cluster, new subtree),\r
-// we re-use one of the two rows that become available (the children\r
-// of the new node). This saves memory.\r
-// We keep track of this through the g_uNodeIndex vector.\r
-static dist_t *g_Dist;\r
-\r
-// Distance to nearest neighbor in row i of distance matrix.\r
-// Subscript is distance matrix row.\r
-static dist_t *g_MinDist;\r
-\r
-// Nearest neighbor to row i of distance matrix.\r
-// Subscript is distance matrix row.\r
-static unsigned *g_uNearestNeighbor;\r
-\r
-// Node index of row i in distance matrix.\r
-// Node indexes are 0..N-1 for leaves, N..2N-2 for internal nodes.\r
-// Subscript is distance matrix row.\r
-static unsigned *g_uNodeIndex;\r
-\r
-// The following vectors are defined on internal nodes,\r
-// subscripts are internal node index 0..N-2.\r
-// For g_uLeft/Right, value is the node index 0 .. 2N-2\r
-// because a child can be internal or leaf.\r
-static unsigned *g_uLeft;\r
-static unsigned *g_uRight;\r
-static dist_t *g_Height;\r
-static dist_t *g_LeftLength;\r
-static dist_t *g_RightLength;\r
-\r
-static inline unsigned TriangleSubscript(unsigned uIndex1, unsigned uIndex2)\r
- {\r
-#if DEBUG\r
- if (uIndex1 >= g_uLeafCount || uIndex2 >= g_uLeafCount)\r
- Quit("TriangleSubscript(%u,%u) %u", uIndex1, uIndex2, g_uLeafCount);\r
-#endif\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 < (g_uLeafCount*(g_uLeafCount - 1))/2);\r
- return v;\r
- }\r
-\r
-static void ListState()\r
- {\r
- Log("Dist matrix\n");\r
- Log(" ");\r
- for (unsigned i = 0; i < g_uLeafCount; ++i)\r
- {\r
- if (uInsane == g_uNodeIndex[i])\r
- continue;\r
- Log(" %5u", g_uNodeIndex[i]);\r
- }\r
- Log("\n");\r
-\r
- for (unsigned i = 0; i < g_uLeafCount; ++i)\r
- {\r
- if (uInsane == g_uNodeIndex[i])\r
- continue;\r
- Log("%5u ", g_uNodeIndex[i]);\r
- for (unsigned j = 0; j < g_uLeafCount; ++j)\r
- {\r
- if (uInsane == g_uNodeIndex[j])\r
- continue;\r
- if (i == j)\r
- Log(" ");\r
- else\r
- {\r
- unsigned v = TriangleSubscript(i, j);\r
- Log("%5.2g ", g_Dist[v]);\r
- }\r
- }\r
- Log("\n");\r
- }\r
-\r
- Log("\n");\r
- Log(" i Node NrNb Dist\n");\r
- Log("----- ----- ----- --------\n");\r
- for (unsigned i = 0; i < g_uLeafCount; ++i)\r
- {\r
- if (uInsane == g_uNodeIndex[i])\r
- continue;\r
- Log("%5u %5u %5u %8.3f\n",\r
- i,\r
- g_uNodeIndex[i],\r
- g_uNearestNeighbor[i],\r
- g_MinDist[i]);\r
- }\r
-\r
- Log("\n");\r
- Log(" Node L R Height LLength RLength\n");\r
- Log("----- ----- ----- ------ ------- -------\n");\r
- for (unsigned i = 0; i <= g_uInternalNodeIndex; ++i)\r
- Log("%5u %5u %5u %6.2g %6.2g %6.2g\n",\r
- i,\r
- g_uLeft[i],\r
- g_uRight[i],\r
- g_Height[i],\r
- g_LeftLength[i],\r
- g_RightLength[i]);\r
- }\r
-\r
-void UPGMA2(const DistCalc &DC, Tree &tree, LINKAGE Linkage)\r
- {\r
- g_uLeafCount = DC.GetCount();\r
-\r
- g_uTriangleSize = (g_uLeafCount*(g_uLeafCount - 1))/2;\r
- g_uInternalNodeCount = g_uLeafCount - 1;\r
-\r
- g_Dist = new dist_t[g_uTriangleSize];\r
-\r
- g_uNodeIndex = new unsigned[g_uLeafCount];\r
- g_uNearestNeighbor = new unsigned[g_uLeafCount];\r
- g_MinDist = new dist_t[g_uLeafCount];\r
- unsigned *Ids = new unsigned [g_uLeafCount];\r
- char **Names = new char *[g_uLeafCount];\r
-\r
- g_uLeft = new unsigned[g_uInternalNodeCount];\r
- g_uRight = new unsigned[g_uInternalNodeCount];\r
- g_Height = new dist_t[g_uInternalNodeCount];\r
- g_LeftLength = new dist_t[g_uInternalNodeCount];\r
- g_RightLength = new dist_t[g_uInternalNodeCount];\r
-\r
- for (unsigned i = 0; i < g_uLeafCount; ++i)\r
- {\r
- g_MinDist[i] = BIG_DIST;\r
- g_uNodeIndex[i] = i;\r
- g_uNearestNeighbor[i] = uInsane;\r
- Ids[i] = DC.GetId(i);\r
- Names[i] = strsave(DC.GetName(i));\r
- }\r
-\r
- for (unsigned i = 0; i < g_uInternalNodeCount; ++i)\r
- {\r
- g_uLeft[i] = uInsane;\r
- g_uRight[i] = uInsane;\r
- g_LeftLength[i] = BIG_DIST;\r
- g_RightLength[i] = BIG_DIST;\r
- g_Height[i] = BIG_DIST;\r
- }\r
-\r
-// Compute initial NxN triangular distance matrix.\r
-// Store minimum distance for each full (not triangular) row.\r
-// Loop from 1, not 0, because "row" is 0, 1 ... i-1,\r
-// so nothing to do when i=0.\r
- for (unsigned i = 1; i < g_uLeafCount; ++i)\r
- {\r
- dist_t *Row = g_Dist + TriangleSubscript(i, 0);\r
- DC.CalcDistRange(i, Row);\r
- for (unsigned j = 0; j < i; ++j)\r
- {\r
- const dist_t d = Row[j];\r
- if (d < g_MinDist[i])\r
- {\r
- g_MinDist[i] = d;\r
- g_uNearestNeighbor[i] = j;\r
- }\r
- if (d < g_MinDist[j])\r
- {\r
- g_MinDist[j] = d;\r
- g_uNearestNeighbor[j] = i;\r
- }\r
- }\r
- }\r
-\r
-#if TRACE\r
- Log("Initial state:\n");\r
- ListState();\r
-#endif\r
-\r
- for (g_uInternalNodeIndex = 0; g_uInternalNodeIndex < g_uLeafCount - 1;\r
- ++g_uInternalNodeIndex)\r
- {\r
-#if TRACE\r
- Log("\n");\r
- Log("Internal node index %5u\n", g_uInternalNodeIndex);\r
- Log("-------------------------\n");\r
-#endif\r
-\r
- // Find nearest neighbors\r
- unsigned Lmin = uInsane;\r
- unsigned Rmin = uInsane;\r
- dist_t dtMinDist = BIG_DIST;\r
- for (unsigned j = 0; j < g_uLeafCount; ++j)\r
- {\r
- if (uInsane == g_uNodeIndex[j])\r
- continue;\r
-\r
- dist_t d = g_MinDist[j];\r
- if (d < dtMinDist)\r
- {\r
- dtMinDist = d;\r
- Lmin = j;\r
- Rmin = g_uNearestNeighbor[j];\r
- assert(uInsane != Rmin);\r
- assert(uInsane != g_uNodeIndex[Rmin]);\r
- }\r
- }\r
-\r
- assert(Lmin != uInsane);\r
- assert(Rmin != uInsane);\r
- assert(dtMinDist != BIG_DIST);\r
-\r
-#if TRACE\r
- Log("Nearest neighbors Lmin %u[=%u] Rmin %u[=%u] dist %.3g\n",\r
- Lmin,\r
- g_uNodeIndex[Lmin],\r
- Rmin,\r
- g_uNodeIndex[Rmin],\r
- dtMinDist);\r
-#endif\r
-\r
- // Compute distances to new node\r
- // New node overwrites row currently assigned to Lmin\r
- dist_t dtNewMinDist = BIG_DIST;\r
- unsigned uNewNearestNeighbor = uInsane;\r
- for (unsigned j = 0; j < g_uLeafCount; ++j)\r
- {\r
- if (j == Lmin || j == Rmin)\r
- continue;\r
- if (uInsane == g_uNodeIndex[j])\r
- continue;\r
-\r
- const unsigned vL = TriangleSubscript(Lmin, j);\r
- const unsigned vR = TriangleSubscript(Rmin, j);\r
- const dist_t dL = g_Dist[vL];\r
- const dist_t dR = g_Dist[vR];\r
- dist_t dtNewDist;\r
-\r
- switch (Linkage)\r
- {\r
- case LINKAGE_Avg:\r
- dtNewDist = AVG(dL, dR);\r
- break;\r
-\r
- case LINKAGE_Min:\r
- dtNewDist = MIN(dL, dR);\r
- break;\r
-\r
- case LINKAGE_Max:\r
- dtNewDist = MAX(dL, dR);\r
- break;\r
-\r
- case LINKAGE_Biased:\r
- dtNewDist = g_dSUEFF*AVG(dL, dR) + (1 - g_dSUEFF)*MIN(dL, dR);\r
- break;\r
-\r
- default:\r
- Quit("UPGMA2: Invalid LINKAGE_%u", Linkage);\r
- }\r
-\r
- // Nasty special case.\r
- // If nearest neighbor of j is Lmin or Rmin, then make the new\r
- // node (which overwrites the row currently occupied by Lmin)\r
- // the nearest neighbor. This situation can occur when there are\r
- // equal distances in the matrix. If we don't make this fix,\r
- // the nearest neighbor pointer for j would become invalid.\r
- // (We don't need to test for == Lmin, because in that case\r
- // the net change needed is zero due to the change in row\r
- // numbering).\r
- if (g_uNearestNeighbor[j] == Rmin)\r
- g_uNearestNeighbor[j] = Lmin;\r
-\r
-#if TRACE\r
- Log("New dist to %u = (%u/%.3g + %u/%.3g)/2 = %.3g\n",\r
- j, Lmin, dL, Rmin, dR, dtNewDist);\r
-#endif\r
- g_Dist[vL] = dtNewDist;\r
- if (dtNewDist < dtNewMinDist)\r
- {\r
- dtNewMinDist = dtNewDist;\r
- uNewNearestNeighbor = j;\r
- }\r
- }\r
-\r
- assert(g_uInternalNodeIndex < g_uLeafCount - 1 || BIG_DIST != dtNewMinDist);\r
- assert(g_uInternalNodeIndex < g_uLeafCount - 1 || uInsane != uNewNearestNeighbor);\r
-\r
- const unsigned v = TriangleSubscript(Lmin, Rmin);\r
- const dist_t dLR = g_Dist[v];\r
- const dist_t dHeightNew = dLR/2;\r
- const unsigned uLeft = g_uNodeIndex[Lmin];\r
- const unsigned uRight = g_uNodeIndex[Rmin];\r
- const dist_t HeightLeft =\r
- uLeft < g_uLeafCount ? 0 : g_Height[uLeft - g_uLeafCount];\r
- const dist_t HeightRight =\r
- uRight < g_uLeafCount ? 0 : g_Height[uRight - g_uLeafCount];\r
-\r
- g_uLeft[g_uInternalNodeIndex] = uLeft;\r
- g_uRight[g_uInternalNodeIndex] = uRight;\r
- g_LeftLength[g_uInternalNodeIndex] = dHeightNew - HeightLeft;\r
- g_RightLength[g_uInternalNodeIndex] = dHeightNew - HeightRight;\r
- g_Height[g_uInternalNodeIndex] = dHeightNew;\r
-\r
- // Row for left child overwritten by row for new node\r
- g_uNodeIndex[Lmin] = g_uLeafCount + g_uInternalNodeIndex;\r
- g_uNearestNeighbor[Lmin] = uNewNearestNeighbor;\r
- g_MinDist[Lmin] = dtNewMinDist;\r
-\r
- // Delete row for right child\r
- g_uNodeIndex[Rmin] = uInsane;\r
-\r
-#if TRACE\r
- Log("\nInternalNodeIndex=%u Lmin=%u Rmin=%u\n",\r
- g_uInternalNodeIndex, Lmin, Rmin);\r
- ListState();\r
-#endif\r
- }\r
-\r
- unsigned uRoot = g_uLeafCount - 2;\r
- tree.Create(g_uLeafCount, uRoot, g_uLeft, g_uRight, g_LeftLength, g_RightLength,\r
- Ids, Names);\r
-\r
-#if TRACE\r
- tree.LogMe();\r
-#endif\r
-\r
- delete[] g_Dist;\r
-\r
- delete[] g_uNodeIndex;\r
- delete[] g_uNearestNeighbor;\r
- delete[] g_MinDist;\r
- delete[] g_Height;\r
-\r
- delete[] g_uLeft;\r
- delete[] g_uRight;\r
- delete[] g_LeftLength;\r
- delete[] g_RightLength;\r
- \r
- for (unsigned i = 0; i < g_uLeafCount; ++i)\r
- free(Names[i]);\r
- delete[] Names;\r
- delete[] Ids;\r
- }\r
-\r
-class DistCalcTest : public DistCalc\r
- {\r
- virtual void CalcDistRange(unsigned i, dist_t Dist[]) const\r
- {\r
- static dist_t TestDist[5][5] =\r
- {\r
- 0, 2, 14, 14, 20,\r
- 2, 0, 14, 14, 20,\r
- 14, 14, 0, 4, 20,\r
- 14, 14, 4, 0, 20,\r
- 20, 20, 20, 20, 0,\r
- };\r
- for (unsigned j = 0; j < i; ++j)\r
- Dist[j] = TestDist[i][j];\r
- }\r
- virtual unsigned GetCount() const\r
- {\r
- return 5;\r
- }\r
- virtual unsigned GetId(unsigned i) const\r
- {\r
- return i;\r
- }\r
- virtual const char *GetName(unsigned i) const\r
- {\r
- return "name";\r
- }\r
- };\r
-\r
-void Test()\r
- {\r
- SetListFileName("c:\\tmp\\lobster.log", false);\r
- DistCalcTest DC;\r
- Tree tree;\r
- UPGMA2(DC, tree, LINKAGE_Avg);\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff