binaries/src/muscle/phy2.cpp

   1 #include "muscle.h"\r
   2 #include "tree.h"\r
   3 \r
   4 #define TRACE   0\r
   5 \r
   6 // Return false when done\r
   7 bool PhyEnumEdges(const Tree &tree, PhyEnumEdgeState &ES)\r
   8         {\r
   9         unsigned uNode1 = uInsane;\r
  10 \r
  11         if (!ES.m_bInit)\r
  12                 {\r
  13                 if (tree.GetNodeCount() <= 1)\r
  14                         {\r
  15                         ES.m_uNodeIndex1 = NULL_NEIGHBOR;\r
  16                         ES.m_uNodeIndex2 = NULL_NEIGHBOR;\r
  17                         return false;\r
  18                         }\r
  19                 uNode1 = tree.FirstDepthFirstNode();\r
  20                 ES.m_bInit = true;\r
  21                 }\r
  22         else\r
  23                 {\r
  24                 uNode1 = tree.NextDepthFirstNode(ES.m_uNodeIndex1);\r
  25                 if (NULL_NEIGHBOR == uNode1)\r
  26                         return false;\r
  27                 if (tree.IsRooted() && tree.IsRoot(uNode1))\r
  28                         {\r
  29                         uNode1 = tree.NextDepthFirstNode(uNode1);\r
  30                         if (NULL_NEIGHBOR == uNode1)\r
  31                                 return false;\r
  32                         }\r
  33                 }\r
  34         unsigned uNode2 = tree.GetParent(uNode1);\r
  35 \r
  36         ES.m_uNodeIndex1 = uNode1;\r
  37         ES.m_uNodeIndex2 = uNode2;\r
  38         return true;\r
  39         }\r
  40 \r
  41 bool PhyEnumEdgesR(const Tree &tree, PhyEnumEdgeState &ES)\r
  42         {\r
  43         unsigned uNode1 = uInsane;\r
  44 \r
  45         if (!ES.m_bInit)\r
  46                 {\r
  47                 if (tree.GetNodeCount() <= 1)\r
  48                         {\r
  49                         ES.m_uNodeIndex1 = NULL_NEIGHBOR;\r
  50                         ES.m_uNodeIndex2 = NULL_NEIGHBOR;\r
  51                         return false;\r
  52                         }\r
  53                 uNode1 = tree.FirstDepthFirstNodeR();\r
  54                 ES.m_bInit = true;\r
  55                 }\r
  56         else\r
  57                 {\r
  58                 uNode1 = tree.NextDepthFirstNodeR(ES.m_uNodeIndex1);\r
  59                 if (NULL_NEIGHBOR == uNode1)\r
  60                         return false;\r
  61                 if (tree.IsRooted() && tree.IsRoot(uNode1))\r
  62                         {\r
  63                         uNode1 = tree.NextDepthFirstNode(uNode1);\r
  64                         if (NULL_NEIGHBOR == uNode1)\r
  65                                 return false;\r
  66                         }\r
  67                 }\r
  68         unsigned uNode2 = tree.GetParent(uNode1);\r
  69 \r
  70         ES.m_uNodeIndex1 = uNode1;\r
  71         ES.m_uNodeIndex2 = uNode2;\r
  72         return true;\r
  73         }\r
  74 \r
  75 static void GetLeavesSubtree(const Tree &tree, unsigned uNodeIndex1,\r
  76   const unsigned uNodeIndex2, unsigned Leaves[], unsigned *ptruCount)\r
  77         {\r
  78         if (tree.IsLeaf(uNodeIndex1))\r
  79                 {\r
  80                 Leaves[*ptruCount] = uNodeIndex1;\r
  81                 ++(*ptruCount);\r
  82                 return;\r
  83                 }\r
  84 \r
  85         const unsigned uLeft = tree.GetFirstNeighbor(uNodeIndex1, uNodeIndex2);\r
  86         const unsigned uRight = tree.GetSecondNeighbor(uNodeIndex1, uNodeIndex2);\r
  87         if (NULL_NEIGHBOR != uLeft)\r
  88                 GetLeavesSubtree(tree, uLeft, uNodeIndex1, Leaves, ptruCount);\r
  89         if (NULL_NEIGHBOR != uRight)\r
  90                 GetLeavesSubtree(tree, uRight, uNodeIndex1, Leaves, ptruCount);\r
  91         }\r
  92 \r
  93 static void PhyGetLeaves(const Tree &tree, unsigned uNodeIndex1, unsigned uNodeIndex2,\r
  94   unsigned Leaves[], unsigned *ptruCount)\r
  95         {\r
  96         *ptruCount = 0;\r
  97         GetLeavesSubtree(tree, uNodeIndex1, uNodeIndex2, Leaves, ptruCount);\r
  98         }\r
  99 \r
 100 bool PhyEnumBiParts(const Tree &tree, PhyEnumEdgeState &ES,\r
 101   unsigned Leaves1[], unsigned *ptruCount1,\r
 102   unsigned Leaves2[], unsigned *ptruCount2)\r
 103         {\r
 104         bool bOk = PhyEnumEdges(tree, ES);\r
 105         if (!bOk)\r
 106                 {\r
 107                 *ptruCount1 = 0;\r
 108                 *ptruCount2 = 0;\r
 109                 return false;\r
 110                 }\r
 111 \r
 112 // Special case: in a rooted tree, both edges from the root\r
 113 // give the same bipartition, so skip one of them.\r
 114         if (tree.IsRooted() && tree.IsRoot(ES.m_uNodeIndex2)\r
 115           && tree.GetRight(ES.m_uNodeIndex2) == ES.m_uNodeIndex1)\r
 116                 {\r
 117                 bOk = PhyEnumEdges(tree, ES);\r
 118                 if (!bOk)\r
 119                         return false;\r
 120                 }\r
 121 \r
 122         PhyGetLeaves(tree, ES.m_uNodeIndex1, ES.m_uNodeIndex2, Leaves1, ptruCount1);\r
 123         PhyGetLeaves(tree, ES.m_uNodeIndex2, ES.m_uNodeIndex1, Leaves2, ptruCount2);\r
 124 \r
 125         if (*ptruCount1 + *ptruCount2 != tree.GetLeafCount())\r
 126                 Quit("PhyEnumBiParts %u + %u != %u",\r
 127                   *ptruCount1, *ptruCount2, tree.GetLeafCount());\r
 128 #if     DEBUG\r
 129         {\r
 130         for (unsigned i = 0; i < *ptruCount1; ++i)\r
 131                 {\r
 132                 if (!tree.IsLeaf(Leaves1[i]))\r
 133                         Quit("PhyEnumByParts: not leaf");\r
 134                 for (unsigned j = 0; j < *ptruCount2; ++j)\r
 135                         {\r
 136                         if (!tree.IsLeaf(Leaves2[j]))\r
 137                                 Quit("PhyEnumByParts: not leaf");\r
 138                         if (Leaves1[i] == Leaves2[j])\r
 139                                 Quit("PhyEnumByParts: dupe");\r
 140                         }\r
 141                 }\r
 142         }\r
 143 #endif\r
 144 \r
 145         return true;\r
 146         }\r
 147 \r
 148 #if     0\r
 149 void TestBiPart()\r
 150         {\r
 151         SetListFileName("c:\\tmp\\lobster.log", false);\r
 152         Tree tree;\r
 153         TextFile fileIn("c:\\tmp\\test.phy");\r
 154         tree.FromFile(fileIn);\r
 155         tree.LogMe();\r
 156 \r
 157         const unsigned uNodeCount = tree.GetNodeCount();\r
 158         unsigned *Leaves1 = new unsigned[uNodeCount];\r
 159         unsigned *Leaves2 = new unsigned[uNodeCount];\r
 160 \r
 161         PhyEnumEdgeState ES;\r
 162         bool bDone = false;\r
 163         for (;;)\r
 164                 {\r
 165                 unsigned uCount1 = uInsane;\r
 166                 unsigned uCount2 = uInsane;\r
 167                 bool bOk = PhyEnumBiParts(tree, ES, Leaves1, &uCount1, Leaves2, &uCount2);\r
 168                 Log("PEBP=%d ES.Init=%d ES.ni1=%d ES.ni2=%d\n",\r
 169                   bOk,\r
 170                   ES.m_bInit,\r
 171                   ES.m_uNodeIndex1,\r
 172                   ES.m_uNodeIndex2);\r
 173                 if (!bOk)\r
 174                         break;\r
 175                 Log("\n");\r
 176                 Log("Part1: ");\r
 177                 for (unsigned n = 0; n < uCount1; ++n)\r
 178                         Log(" %d(%s)", Leaves1[n], tree.GetLeafName(Leaves1[n]));\r
 179                 Log("\n");\r
 180                 Log("Part2: ");\r
 181                 for (unsigned n = 0; n < uCount2; ++n)\r
 182                         Log(" %d(%s)", Leaves2[n], tree.GetLeafName(Leaves2[n]));\r
 183                 Log("\n");\r
 184                 }\r
 185         }\r
 186 #endif\r
 187 \r
 188 static void GetLeavesSubtreeExcluding(const Tree &tree, unsigned uNodeIndex,\r
 189   unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)\r
 190         {\r
 191         if (uNodeIndex == uExclude)\r
 192                 return;\r
 193 \r
 194         if (tree.IsLeaf(uNodeIndex))\r
 195                 {\r
 196                 Leaves[*ptruCount] = uNodeIndex;\r
 197                 ++(*ptruCount);\r
 198                 return;\r
 199                 }\r
 200 \r
 201         const unsigned uLeft = tree.GetLeft(uNodeIndex);\r
 202         const unsigned uRight = tree.GetRight(uNodeIndex);\r
 203         if (NULL_NEIGHBOR != uLeft)\r
 204                 GetLeavesSubtreeExcluding(tree, uLeft, uExclude, Leaves, ptruCount);\r
 205         if (NULL_NEIGHBOR != uRight)\r
 206                 GetLeavesSubtreeExcluding(tree, uRight, uExclude, Leaves, ptruCount);\r
 207         }\r
 208 \r
 209 void GetLeavesExcluding(const Tree &tree, unsigned uNodeIndex,\r
 210   unsigned uExclude, unsigned Leaves[], unsigned *ptruCount)\r
 211         {\r
 212         *ptruCount = 0;\r
 213         GetLeavesSubtreeExcluding(tree, uNodeIndex, uExclude, Leaves, ptruCount);\r
 214         }\r
 215 \r
 216 void GetInternalNodesInHeightOrder(const Tree &tree, unsigned NodeIndexes[])\r
 217         {\r
 218         const unsigned uNodeCount = tree.GetNodeCount();\r
 219         if (uNodeCount < 3)\r
 220                 Quit("GetInternalNodesInHeightOrder: %u nodes, none are internal",\r
 221                   uNodeCount);\r
 222         const unsigned uInternalNodeCount = (uNodeCount - 1)/2;\r
 223         double *Heights = new double[uInternalNodeCount];\r
 224 \r
 225         unsigned uIndex = 0;\r
 226         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
 227                 {\r
 228                 if (tree.IsLeaf(uNodeIndex))\r
 229                         continue;\r
 230                 NodeIndexes[uIndex] = uNodeIndex;\r
 231                 Heights[uIndex] = tree.GetNodeHeight(uNodeIndex);\r
 232                 ++uIndex;\r
 233                 }\r
 234         if (uIndex != uInternalNodeCount)\r
 235                 Quit("Internal error: GetInternalNodesInHeightOrder");\r
 236 \r
 237 // Simple but slow bubble sort (probably don't care about speed here)\r
 238         bool bDone = false;\r
 239         while (!bDone)\r
 240                 {\r
 241                 bDone = true;\r
 242                 for (unsigned i = 0; i < uInternalNodeCount - 1; ++i)\r
 243                         {\r
 244                         if (Heights[i] > Heights[i+1])\r
 245                                 {\r
 246                                 double dTmp = Heights[i];\r
 247                                 Heights[i] = Heights[i+1];\r
 248                                 Heights[i+1] = dTmp;\r
 249 \r
 250                                 unsigned uTmp = NodeIndexes[i];\r
 251                                 NodeIndexes[i] = NodeIndexes[i+1];\r
 252                                 NodeIndexes[i+1] = uTmp;\r
 253                                 bDone = false;\r
 254                                 }\r
 255                         }\r
 256                 }\r
 257 #if     TRACE\r
 258         Log("Internal node index     Height\n");\r
 259         Log("-------------------   --------\n");\r
 260         //    1234567890123456789  123456789\r
 261         for (unsigned n = 0; n < uInternalNodeCount; ++n)\r
 262                 Log("%19u  %9.3f\n", NodeIndexes[n], Heights[n]);\r
 263 #endif\r
 264         delete[] Heights;\r
 265         }\r
 266 \r
 267 void ApplyMinEdgeLength(Tree &tree, double dMinEdgeLength)\r
 268         {\r
 269         const unsigned uNodeCount = tree.GetNodeCount();\r
 270         for (unsigned uNodeIndex = 0; uNodeIndex < uNodeCount; ++uNodeIndex)\r
 271                 {\r
 272                 const unsigned uNeighborCount = tree.GetNeighborCount(uNodeIndex);\r
 273                 for (unsigned n = 0; n < uNeighborCount; ++n)\r
 274                         {\r
 275                         const unsigned uNeighborNodeIndex = tree.GetNeighbor(uNodeIndex, n);\r
 276                         if (!tree.HasEdgeLength(uNodeIndex, uNeighborNodeIndex))\r
 277                                 continue;\r
 278                         if (tree.GetEdgeLength(uNodeIndex, uNeighborNodeIndex) < dMinEdgeLength)\r
 279                                 tree.SetEdgeLength(uNodeIndex, uNeighborNodeIndex, dMinEdgeLength);\r
 280                         }\r
 281                 }\r
 282         }\r