--- /dev/null
+#include "muscle.h"\r
+#include "pwpath.h"\r
+#include "seq.h"\r
+#include "textfile.h"\r
+#include "msa.h"\r
+\r
+PWPath::PWPath()\r
+ {\r
+ m_uArraySize = 0;\r
+ m_uEdgeCount = 0;\r
+ m_Edges = 0;\r
+ }\r
+\r
+PWPath::~PWPath()\r
+ {\r
+ Clear();\r
+ }\r
+\r
+void PWPath::Clear()\r
+ {\r
+ delete[] m_Edges;\r
+ m_Edges = 0;\r
+ m_uArraySize = 0;\r
+ m_uEdgeCount = 0;\r
+ }\r
+\r
+void PWPath::ExpandPath(unsigned uAdditionalEdgeCount)\r
+ {\r
+ PWEdge *OldPath = m_Edges;\r
+ unsigned uEdgeCount = m_uArraySize + uAdditionalEdgeCount;\r
+\r
+ m_Edges = new PWEdge[uEdgeCount];\r
+ m_uArraySize = uEdgeCount;\r
+ if (m_uEdgeCount > 0)\r
+ memcpy(m_Edges, OldPath, m_uEdgeCount*sizeof(PWEdge));\r
+ delete[] OldPath;\r
+ }\r
+\r
+void PWPath::AppendEdge(const PWEdge &Edge)\r
+ {\r
+ if (0 == m_uArraySize || m_uEdgeCount + 1 == m_uArraySize)\r
+ ExpandPath(200);\r
+\r
+ m_Edges[m_uEdgeCount] = Edge;\r
+ ++m_uEdgeCount;\r
+ }\r
+\r
+void PWPath::AppendEdge(char cType, unsigned uPrefixLengthA, unsigned uPrefixLengthB)\r
+ {\r
+ PWEdge e;\r
+ e.uPrefixLengthA = uPrefixLengthA;\r
+ e.uPrefixLengthB = uPrefixLengthB;\r
+ e.cType = cType;\r
+ AppendEdge(e);\r
+ }\r
+\r
+void PWPath::PrependEdge(const PWEdge &Edge)\r
+ {\r
+ if (0 == m_uArraySize || m_uEdgeCount + 1 == m_uArraySize)\r
+ ExpandPath(1000);\r
+ if (m_uEdgeCount > 0)\r
+ memmove(m_Edges + 1, m_Edges, sizeof(PWEdge)*m_uEdgeCount);\r
+ m_Edges[0] = Edge;\r
+ ++m_uEdgeCount;\r
+ }\r
+\r
+const PWEdge &PWPath::GetEdge(unsigned uEdgeIndex) const\r
+ {\r
+ assert(uEdgeIndex < m_uEdgeCount);\r
+ return m_Edges[uEdgeIndex];\r
+ }\r
+\r
+void PWPath::Validate() const\r
+ {\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ if (0 == uEdgeCount)\r
+ return;\r
+ const PWEdge &FirstEdge = GetEdge(0);\r
+ const PWEdge &LastEdge = GetEdge(uEdgeCount - 1);\r
+ unsigned uStartA = FirstEdge.uPrefixLengthA;\r
+ unsigned uStartB = FirstEdge.uPrefixLengthB;\r
+ if (FirstEdge.cType != 'I')\r
+ --uStartA;\r
+ if (FirstEdge.cType != 'D')\r
+ --uStartB;\r
+\r
+ unsigned uPrefixLengthA = FirstEdge.uPrefixLengthA;\r
+ unsigned uPrefixLengthB = FirstEdge.uPrefixLengthB;\r
+ for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = GetEdge(uEdgeIndex);\r
+ switch (Edge.cType)\r
+ {\r
+ case 'M':\r
+ if (uPrefixLengthA + 1 != Edge.uPrefixLengthA)\r
+ Quit("PWPath::Validate MA %u", uPrefixLengthA);\r
+ if (uPrefixLengthB + 1 != Edge.uPrefixLengthB)\r
+ Quit("PWPath::Validate MB %u", uPrefixLengthB);\r
+ ++uPrefixLengthA;\r
+ ++uPrefixLengthB;\r
+ break;\r
+ case 'D':\r
+ if (uPrefixLengthA + 1 != Edge.uPrefixLengthA)\r
+ Quit("PWPath::Validate DA %u", uPrefixLengthA);\r
+ if (uPrefixLengthB != Edge.uPrefixLengthB)\r
+ Quit("PWPath::Validate DB %u", uPrefixLengthB);\r
+ ++uPrefixLengthA;\r
+ break;\r
+ case 'I':\r
+ if (uPrefixLengthA != Edge.uPrefixLengthA)\r
+ Quit("PWPath::Validate IA %u", uPrefixLengthA);\r
+ if (uPrefixLengthB + 1 != Edge.uPrefixLengthB)\r
+ Quit("PWPath::Validate IB %u", uPrefixLengthB);\r
+ ++uPrefixLengthB;\r
+ break;\r
+ }\r
+ }\r
+ }\r
+\r
+void PWPath::LogMe() const\r
+ {\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < GetEdgeCount(); ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = GetEdge(uEdgeIndex);\r
+ if (uEdgeIndex > 0)\r
+ Log(" ");\r
+ Log("%c%d.%d",\r
+ Edge.cType,\r
+ Edge.uPrefixLengthA,\r
+ Edge.uPrefixLengthB);\r
+ if ((uEdgeIndex > 0 && uEdgeIndex%10 == 0) ||\r
+ uEdgeIndex == GetEdgeCount() - 1)\r
+ Log("\n");\r
+ }\r
+ }\r
+\r
+void PWPath::Copy(const PWPath &Path)\r
+ {\r
+ Clear();\r
+ const unsigned uEdgeCount = Path.GetEdgeCount();\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ AppendEdge(Edge);\r
+ }\r
+ }\r
+\r
+void PWPath::FromMSAPair(const MSA &msaA, const MSA &msaB)\r
+ {\r
+ const unsigned uColCount = msaA.GetColCount();\r
+ if (uColCount != msaB.GetColCount())\r
+ Quit("PWPath::FromMSAPair, lengths differ");\r
+\r
+ Clear();\r
+\r
+ unsigned uPrefixLengthA = 0;\r
+ unsigned uPrefixLengthB = 0;\r
+ for (unsigned uColIndex = 0; uColIndex < uColCount; ++uColIndex)\r
+ {\r
+ bool bIsGapA = msaA.IsGapColumn(uColIndex);\r
+ bool bIsGapB = msaB.IsGapColumn(uColIndex);\r
+\r
+ PWEdge Edge;\r
+ char cType;\r
+ if (!bIsGapA && !bIsGapB)\r
+ {\r
+ cType = 'M';\r
+ ++uPrefixLengthA;\r
+ ++uPrefixLengthB;\r
+ }\r
+ else if (bIsGapA && !bIsGapB)\r
+ {\r
+ cType = 'I';\r
+ ++uPrefixLengthB;\r
+ }\r
+ else if (!bIsGapA && bIsGapB)\r
+ {\r
+ cType = 'D';\r
+ ++uPrefixLengthA;\r
+ }\r
+ else\r
+ {\r
+ assert(bIsGapB && bIsGapA);\r
+ continue;\r
+ }\r
+\r
+ Edge.cType = cType;\r
+ Edge.uPrefixLengthA = uPrefixLengthA;\r
+ Edge.uPrefixLengthB = uPrefixLengthB;\r
+ AppendEdge(Edge);\r
+ }\r
+ }\r
+\r
+// Very similar to HMMPath::FromFile, should consolidate.\r
+void PWPath::FromFile(TextFile &File)\r
+ {\r
+ Clear();\r
+ char szToken[1024];\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (0 != strcmp(szToken, "Path"))\r
+ Quit("Invalid path file (Path)");\r
+\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (0 != strcmp(szToken, "edges"))\r
+ Quit("Invalid path file (edges)");\r
+\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (!IsValidInteger(szToken))\r
+ Quit("Invalid path file (edges value)");\r
+\r
+ const unsigned uEdgeCount = (unsigned) atoi(szToken);\r
+ unsigned uEdgeIndex = 0;\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ // index\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (!IsValidInteger(szToken))\r
+ Quit("Invalid path file, invalid index '%s'", szToken);\r
+ unsigned n = (unsigned) atoi(szToken);\r
+ if (n != uEdgeIndex)\r
+ Quit("Invalid path file, expecting edge %u got %u", uEdgeIndex, n);\r
+\r
+ // type\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (1 != strlen(szToken))\r
+ Quit("Invalid path file, expecting state, got '%s'", szToken);\r
+ const char cType = szToken[0];\r
+ if ('M' != cType && 'D' != cType && cType != 'I' && 'S' != cType)\r
+ Quit("Invalid path file, expecting state, got '%c'", cType);\r
+\r
+ // prefix length A\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (!IsValidInteger(szToken))\r
+ Quit("Invalid path file, bad prefix length A '%s'", szToken);\r
+ const unsigned uPrefixLengthA = (unsigned) atoi(szToken);\r
+\r
+ // prefix length B\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (!IsValidInteger(szToken))\r
+ Quit("Invalid path file, bad prefix length B '%s'", szToken);\r
+ const unsigned uPrefixLengthB = (unsigned) atoi(szToken);\r
+\r
+ PWEdge Edge;\r
+ Edge.cType = cType;\r
+ Edge.uPrefixLengthA = uPrefixLengthA;\r
+ Edge.uPrefixLengthB = uPrefixLengthB;\r
+ AppendEdge(Edge);\r
+ }\r
+ File.GetTokenX(szToken, sizeof(szToken));\r
+ if (0 != strcmp(szToken, "//"))\r
+ Quit("Invalid path file (//)");\r
+ }\r
+\r
+void PWPath::ToFile(TextFile &File) const\r
+ {\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+\r
+ File.PutString("Path\n");\r
+ File.PutFormat("edges %u\n", uEdgeCount);\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = GetEdge(uEdgeIndex);\r
+ File.PutFormat("%u %c %u %u\n",\r
+ uEdgeIndex,\r
+ Edge.cType,\r
+ Edge.uPrefixLengthA,\r
+ Edge.uPrefixLengthB);\r
+ }\r
+ File.PutString("//\n");\r
+ }\r
+\r
+void PWPath::AssertEqual(const PWPath &Path) const\r
+ {\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ if (uEdgeCount != Path.GetEdgeCount())\r
+ {\r
+ Log("PWPath::AssertEqual, this=\n");\r
+ LogMe();\r
+ Log("\nOther path=\n");\r
+ Path.LogMe();\r
+ Log("\n");\r
+ Quit("PWPath::AssertEqual, Edge count different %u %u\n",\r
+ uEdgeCount, Path.GetEdgeCount());\r
+ }\r
+\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &e1 = GetEdge(uEdgeIndex);\r
+ const PWEdge &e2 = Path.GetEdge(uEdgeIndex);\r
+ if (e1.cType != e2.cType || e1.uPrefixLengthA != e2.uPrefixLengthA ||\r
+ e1.uPrefixLengthB != e2.uPrefixLengthB)\r
+ {\r
+ Log("PWPath::AssertEqual, this=\n");\r
+ LogMe();\r
+ Log("\nOther path=\n");\r
+ Path.LogMe();\r
+ Log("\n");\r
+ Log("This edge %c%u.%u, other edge %c%u.%u\n",\r
+ e1.cType, e1.uPrefixLengthA, e1.uPrefixLengthB,\r
+ e2.cType, e2.uPrefixLengthA, e2.uPrefixLengthB);\r
+ Quit("PWPath::AssertEqual, edge %u different\n", uEdgeIndex);\r
+ }\r
+ }\r
+ }\r
+\r
+bool PWPath::Equal(const PWPath &Path) const\r
+ {\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ if (uEdgeCount != Path.GetEdgeCount())\r
+ return false;\r
+\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &e1 = GetEdge(uEdgeIndex);\r
+ const PWEdge &e2 = Path.GetEdge(uEdgeIndex);\r
+ if (e1.cType != e2.cType || e1.uPrefixLengthA != e2.uPrefixLengthA ||\r
+ e1.uPrefixLengthB != e2.uPrefixLengthB)\r
+ return false;\r
+ }\r
+ return true;\r
+ }\r
+\r
+unsigned PWPath::GetMatchCount() const\r
+ {\r
+ unsigned uMatchCount = 0;\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &e = GetEdge(uEdgeIndex);\r
+ if ('M' == e.cType)\r
+ ++uMatchCount;\r
+ }\r
+ return uMatchCount;\r
+ }\r
+\r
+unsigned PWPath::GetInsertCount() const\r
+ {\r
+ unsigned uInsertCount = 0;\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &e = GetEdge(uEdgeIndex);\r
+ if ('I' == e.cType)\r
+ ++uInsertCount;\r
+ }\r
+ return uInsertCount;\r
+ }\r
+\r
+unsigned PWPath::GetDeleteCount() const\r
+ {\r
+ unsigned uDeleteCount = 0;\r
+ const unsigned uEdgeCount = GetEdgeCount();\r
+ for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &e = GetEdge(uEdgeIndex);\r
+ if ('D' == e.cType)\r
+ ++uDeleteCount;\r
+ }\r
+ return uDeleteCount;\r
+ }\r
+\r
+void PWPath::FromStr(const char Str[])\r
+ {\r
+ Clear();\r
+ unsigned uPrefixLengthA = 0;\r
+ unsigned uPrefixLengthB = 0;\r
+ while (char c = *Str++)\r
+ {\r
+ switch (c)\r
+ {\r
+ case 'M':\r
+ ++uPrefixLengthA;\r
+ ++uPrefixLengthB;\r
+ break;\r
+ case 'D':\r
+ ++uPrefixLengthA;\r
+ break;\r
+ case 'I':\r
+ ++uPrefixLengthB;\r
+ break;\r
+ default:\r
+ Quit("PWPath::FromStr, invalid state %c", c);\r
+ }\r
+ AppendEdge(c, uPrefixLengthA, uPrefixLengthB);\r
+ }\r
+ }\r