--- /dev/null
+#include "muscle.h"\r
+#include "pwpath.h"\r
+#include "estring.h"\r
+#include "seq.h"\r
+#include "msa.h"\r
+\r
+/***\r
+An "estring" is an edit string that operates on a sequence.\r
+An estring is represented as a vector of integers.\r
+It is interpreted in order of increasing suffix.\r
+A positive value n means copy n letters.\r
+A negative value -n means insert n indels.\r
+Zero marks the end of the vector.\r
+Consecutive entries must have opposite sign, i.e. the\r
+shortest possible representation must be used.\r
+\r
+A "tpair" is a traceback path for a pairwise alignment\r
+represented as two estrings, one for each sequence.\r
+***/\r
+\r
+#define c2(c,d) (((unsigned char) c) << 8 | (unsigned char) d)\r
+\r
+unsigned LengthEstring(const short es[])\r
+ {\r
+ unsigned i = 0;\r
+ while (*es++ != 0)\r
+ ++i;\r
+ return i;\r
+ }\r
+\r
+short *EstringNewCopy(const short es[])\r
+ {\r
+ unsigned n = LengthEstring(es) + 1;\r
+ short *esNew = new short[n];\r
+ memcpy(esNew, es, n*sizeof(short));\r
+ return esNew;\r
+ }\r
+\r
+void LogEstring(const short es[])\r
+ {\r
+ Log("<");\r
+ for (unsigned i = 0; es[i] != 0; ++i)\r
+ {\r
+ if (i > 0)\r
+ Log(" ");\r
+ Log("%d", es[i]);\r
+ }\r
+ Log(">");\r
+ }\r
+\r
+static bool EstringsEq(const short es1[], const short es2[])\r
+ {\r
+ for (;;)\r
+ {\r
+ if (*es1 != *es2)\r
+ return false;\r
+ if (0 == *es1)\r
+ break;\r
+ ++es1;\r
+ ++es2;\r
+ }\r
+ return true;\r
+ }\r
+\r
+static void EstringCounts(const short es[], unsigned *ptruSymbols,\r
+ unsigned *ptruIndels)\r
+ {\r
+ unsigned uSymbols = 0;\r
+ unsigned uIndels = 0;\r
+ for (unsigned i = 0; es[i] != 0; ++i)\r
+ {\r
+ short n = es[i];\r
+ if (n > 0)\r
+ uSymbols += n;\r
+ else if (n < 0)\r
+ uIndels += -n;\r
+ }\r
+ *ptruSymbols = uSymbols;\r
+ *ptruIndels = uIndels;\r
+ }\r
+\r
+static char *EstringOp(const short es[], const char s[])\r
+ {\r
+ unsigned uSymbols;\r
+ unsigned uIndels;\r
+ EstringCounts(es, &uSymbols, &uIndels);\r
+ assert((unsigned) strlen(s) == uSymbols);\r
+ char *sout = new char[uSymbols + uIndels + 1];\r
+ char *psout = sout;\r
+ for (;;)\r
+ {\r
+ int n = *es++;\r
+ if (0 == n)\r
+ break;\r
+ if (n > 0)\r
+ for (int i = 0; i < n; ++i)\r
+ *psout++ = *s++;\r
+ else\r
+ for (int i = 0; i < -n; ++i)\r
+ *psout++ = '-';\r
+ }\r
+ assert(0 == *s);\r
+ *psout = 0;\r
+ return sout;\r
+ }\r
+\r
+void EstringOp(const short es[], const Seq &sIn, Seq &sOut)\r
+ {\r
+#if DEBUG\r
+ unsigned uSymbols;\r
+ unsigned uIndels;\r
+ EstringCounts(es, &uSymbols, &uIndels);\r
+ assert(sIn.Length() == uSymbols);\r
+#endif\r
+ sOut.Clear();\r
+ sOut.SetName(sIn.GetName());\r
+ int p = 0;\r
+ for (;;)\r
+ {\r
+ int n = *es++;\r
+ if (0 == n)\r
+ break;\r
+ if (n > 0)\r
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ char c = sIn[p++];\r
+ sOut.push_back(c);\r
+ }\r
+ else\r
+ for (int i = 0; i < -n; ++i)\r
+ sOut.push_back('-');\r
+ }\r
+ }\r
+\r
+unsigned EstringOp(const short es[], const Seq &sIn, MSA &a)\r
+ {\r
+ unsigned uSymbols;\r
+ unsigned uIndels;\r
+ EstringCounts(es, &uSymbols, &uIndels);\r
+ assert(sIn.Length() == uSymbols);\r
+\r
+ unsigned uColCount = uSymbols + uIndels;\r
+\r
+ a.Clear();\r
+ a.SetSize(1, uColCount);\r
+\r
+ a.SetSeqName(0, sIn.GetName());\r
+ a.SetSeqId(0, sIn.GetId());\r
+\r
+ unsigned p = 0;\r
+ unsigned uColIndex = 0;\r
+ for (;;)\r
+ {\r
+ int n = *es++;\r
+ if (0 == n)\r
+ break;\r
+ if (n > 0)\r
+ for (int i = 0; i < n; ++i)\r
+ {\r
+ char c = sIn[p++];\r
+ a.SetChar(0, uColIndex++, c);\r
+ }\r
+ else\r
+ for (int i = 0; i < -n; ++i)\r
+ a.SetChar(0, uColIndex++, '-');\r
+ }\r
+ assert(uColIndex == uColCount);\r
+ return uColCount;\r
+ }\r
+\r
+void PathToEstrings(const PWPath &Path, short **ptresA, short **ptresB)\r
+ {\r
+// First pass to determine size of estrings esA and esB\r
+ const unsigned uEdgeCount = Path.GetEdgeCount();\r
+ if (0 == uEdgeCount)\r
+ {\r
+ short *esA = new short[1];\r
+ short *esB = new short[1];\r
+ esA[0] = 0;\r
+ esB[0] = 0;\r
+ *ptresA = esA;\r
+ *ptresB = esB;\r
+ return;\r
+ }\r
+\r
+ unsigned iLengthA = 1;\r
+ unsigned iLengthB = 1;\r
+ const char cFirstEdgeType = Path.GetEdge(0).cType;\r
+ char cPrevEdgeType = cFirstEdgeType;\r
+ for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ char cEdgeType = Edge.cType;\r
+\r
+ switch (c2(cPrevEdgeType, cEdgeType))\r
+ {\r
+ case c2('M', 'M'):\r
+ case c2('D', 'D'):\r
+ case c2('I', 'I'):\r
+ break;\r
+\r
+ case c2('D', 'M'):\r
+ case c2('M', 'D'):\r
+ ++iLengthB;\r
+ break;\r
+\r
+ case c2('I', 'M'):\r
+ case c2('M', 'I'):\r
+ ++iLengthA;\r
+ break;\r
+\r
+ case c2('I', 'D'):\r
+ case c2('D', 'I'):\r
+ ++iLengthB;\r
+ ++iLengthA;\r
+ break;\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+ cPrevEdgeType = cEdgeType;\r
+ }\r
+\r
+// Pass2 for seq A\r
+ {\r
+ short *esA = new short[iLengthA+1];\r
+ unsigned iA = 0;\r
+ switch (Path.GetEdge(0).cType)\r
+ {\r
+ case 'M':\r
+ case 'D':\r
+ esA[0] = 1;\r
+ break;\r
+\r
+ case 'I':\r
+ esA[0] = -1;\r
+ break;\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+\r
+ char cPrevEdgeType = cFirstEdgeType;\r
+ for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ char cEdgeType = Edge.cType;\r
+\r
+ switch (c2(cPrevEdgeType, cEdgeType))\r
+ {\r
+ case c2('M', 'M'):\r
+ case c2('D', 'D'):\r
+ case c2('D', 'M'):\r
+ case c2('M', 'D'):\r
+ ++(esA[iA]);\r
+ break;\r
+\r
+ case c2('I', 'D'):\r
+ case c2('I', 'M'):\r
+ ++iA;\r
+ esA[iA] = 1;\r
+ break;\r
+\r
+ case c2('M', 'I'):\r
+ case c2('D', 'I'):\r
+ ++iA;\r
+ esA[iA] = -1;\r
+ break;\r
+\r
+ case c2('I', 'I'):\r
+ --(esA[iA]);\r
+ break;\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+\r
+ cPrevEdgeType = cEdgeType;\r
+ }\r
+ assert(iA == iLengthA - 1);\r
+ esA[iLengthA] = 0;\r
+ *ptresA = esA;\r
+ }\r
+\r
+ {\r
+// Pass2 for seq B\r
+ short *esB = new short[iLengthB+1];\r
+ unsigned iB = 0;\r
+ switch (Path.GetEdge(0).cType)\r
+ {\r
+ case 'M':\r
+ case 'I':\r
+ esB[0] = 1;\r
+ break;\r
+\r
+ case 'D':\r
+ esB[0] = -1;\r
+ break;\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+\r
+ char cPrevEdgeType = cFirstEdgeType;\r
+ for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
+ {\r
+ const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
+ char cEdgeType = Edge.cType;\r
+\r
+ switch (c2(cPrevEdgeType, cEdgeType))\r
+ {\r
+ case c2('M', 'M'):\r
+ case c2('I', 'I'):\r
+ case c2('I', 'M'):\r
+ case c2('M', 'I'):\r
+ ++(esB[iB]);\r
+ break;\r
+\r
+ case c2('D', 'I'):\r
+ case c2('D', 'M'):\r
+ ++iB;\r
+ esB[iB] = 1;\r
+ break;\r
+\r
+ case c2('M', 'D'):\r
+ case c2('I', 'D'):\r
+ ++iB;\r
+ esB[iB] = -1;\r
+ break;\r
+\r
+ case c2('D', 'D'):\r
+ --(esB[iB]);\r
+ break;\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+\r
+ cPrevEdgeType = cEdgeType;\r
+ }\r
+ assert(iB == iLengthB - 1);\r
+ esB[iLengthB] = 0;\r
+ *ptresB = esB;\r
+ }\r
+\r
+#if DEBUG\r
+ {\r
+ const PWEdge &LastEdge = Path.GetEdge(uEdgeCount - 1);\r
+ unsigned uSymbols;\r
+ unsigned uIndels;\r
+ EstringCounts(*ptresA, &uSymbols, &uIndels);\r
+ assert(uSymbols == LastEdge.uPrefixLengthA);\r
+ assert(uSymbols + uIndels == uEdgeCount);\r
+\r
+ EstringCounts(*ptresB, &uSymbols, &uIndels);\r
+ assert(uSymbols == LastEdge.uPrefixLengthB);\r
+ assert(uSymbols + uIndels == uEdgeCount);\r
+\r
+ PWPath TmpPath;\r
+ EstringsToPath(*ptresA, *ptresB, TmpPath);\r
+ TmpPath.AssertEqual(Path);\r
+ }\r
+#endif\r
+ }\r
+\r
+void EstringsToPath(const short esA[], const short esB[], PWPath &Path)\r
+ {\r
+ Path.Clear();\r
+ unsigned iA = 0;\r
+ unsigned iB = 0;\r
+ int nA = esA[iA++];\r
+ int nB = esB[iB++];\r
+ unsigned uPrefixLengthA = 0;\r
+ unsigned uPrefixLengthB = 0;\r
+ for (;;)\r
+ {\r
+ char cType;\r
+ if (nA > 0)\r
+ {\r
+ if (nB > 0)\r
+ {\r
+ cType = 'M';\r
+ --nA;\r
+ --nB;\r
+ }\r
+ else if (nB < 0)\r
+ {\r
+ cType = 'D';\r
+ --nA;\r
+ ++nB;\r
+ }\r
+ else\r
+ assert(false);\r
+ }\r
+ else if (nA < 0)\r
+ {\r
+ if (nB > 0)\r
+ {\r
+ cType = 'I';\r
+ ++nA;\r
+ --nB;\r
+ }\r
+ else\r
+ assert(false);\r
+ }\r
+ else\r
+ assert(false);\r
+\r
+ switch (cType)\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
+ }\r
+\r
+ PWEdge Edge;\r
+ Edge.cType = cType;\r
+ Edge.uPrefixLengthA = uPrefixLengthA;\r
+ Edge.uPrefixLengthB = uPrefixLengthB;\r
+ Path.AppendEdge(Edge);\r
+\r
+ if (nA == 0)\r
+ {\r
+ if (0 == esA[iA])\r
+ {\r
+ assert(0 == esB[iB]);\r
+ break;\r
+ }\r
+ nA = esA[iA++];\r
+ }\r
+ if (nB == 0)\r
+ nB = esB[iB++];\r
+ }\r
+ }\r
+\r
+/***\r
+Multiply two estrings to make a third estring.\r
+The product of two estrings e1*e2 is defined to be\r
+the estring that produces the same result as applying\r
+e1 then e2. Multiplication is not commutative. In fact,\r
+the reversed order is undefined unless both estrings\r
+consist of a single, identical, positive entry.\r
+A primary motivation for using estrings is that\r
+multiplication is very fast, reducing the time\r
+needed to construct the root alignment.\r
+\r
+Example\r
+\r
+ <-1,3>(XXX) = -XXX\r
+ <2,-1,2>(-XXX) = -X-XX\r
+\r
+Therefore,\r
+\r
+ <-1,3>*<2,-1,2> = <-1,1,-1,2>\r
+***/\r
+\r
+static bool CanMultiplyEstrings(const short es1[], const short es2[])\r
+ {\r
+ unsigned uSymbols1;\r
+ unsigned uSymbols2;\r
+ unsigned uIndels1;\r
+ unsigned uIndels2;\r
+ EstringCounts(es1, &uSymbols1, &uIndels1);\r
+ EstringCounts(es2, &uSymbols2, &uIndels2);\r
+ return uSymbols1 + uIndels1 == uSymbols2;\r
+ }\r
+\r
+static inline void AppendGaps(short esp[], int &ip, int n)\r
+ {\r
+ if (-1 == ip)\r
+ esp[++ip] = n;\r
+ else if (esp[ip] < 0)\r
+ esp[ip] += n;\r
+ else\r
+ esp[++ip] = n;\r
+ }\r
+\r
+static inline void AppendSymbols(short esp[], int &ip, int n)\r
+ {\r
+ if (-1 == ip)\r
+ esp[++ip] = n;\r
+ else if (esp[ip] > 0)\r
+ esp[ip] += n;\r
+ else\r
+ esp[++ip] = n;\r
+ }\r
+\r
+void MulEstrings(const short es1[], const short es2[], short esp[])\r
+ {\r
+ assert(CanMultiplyEstrings(es1, es2));\r
+\r
+ unsigned i1 = 0;\r
+ int ip = -1;\r
+ int n1 = es1[i1++];\r
+ for (unsigned i2 = 0; ; ++i2)\r
+ {\r
+ int n2 = es2[i2];\r
+ if (0 == n2)\r
+ break;\r
+ if (n2 > 0)\r
+ {\r
+ for (;;)\r
+ {\r
+ if (n1 < 0)\r
+ {\r
+ if (n2 > -n1)\r
+ {\r
+ AppendGaps(esp, ip, n1);\r
+ n2 += n1;\r
+ n1 = es1[i1++];\r
+ }\r
+ else if (n2 == -n1)\r
+ {\r
+ AppendGaps(esp, ip, n1);\r
+ n1 = es1[i1++];\r
+ break;\r
+ }\r
+ else\r
+ {\r
+ assert(n2 < -n1);\r
+ AppendGaps(esp, ip, -n2);\r
+ n1 += n2;\r
+ break;\r
+ }\r
+ }\r
+ else\r
+ {\r
+ assert(n1 > 0);\r
+ if (n2 > n1)\r
+ {\r
+ AppendSymbols(esp, ip, n1);\r
+ n2 -= n1;\r
+ n1 = es1[i1++];\r
+ }\r
+ else if (n2 == n1)\r
+ {\r
+ AppendSymbols(esp, ip, n1);\r
+ n1 = es1[i1++];\r
+ break;\r
+ }\r
+ else\r
+ {\r
+ assert(n2 < n1);\r
+ AppendSymbols(esp, ip, n2);\r
+ n1 -= n2;\r
+ break;\r
+ }\r
+ }\r
+ }\r
+ }\r
+ else\r
+ {\r
+ assert(n2 < 0);\r
+ AppendGaps(esp, ip, n2);\r
+ }\r
+ }\r
+ esp[++ip] = 0;\r
+\r
+#if DEBUG\r
+ {\r
+ int MaxLen = (int) (LengthEstring(es1) + LengthEstring(es2) + 1);\r
+ assert(ip < MaxLen);\r
+ if (ip >= 2)\r
+ for (int i = 0; i < ip - 2; ++i)\r
+ {\r
+ if (!(esp[i] > 0 && esp[i+1] < 0 || esp[i] < 0 && esp[i+1] > 0))\r
+ {\r
+ Log("Bad result of MulEstring: ");\r
+ LogEstring(esp);\r
+ Quit("Assert failed (alternating signs)");\r
+ }\r
+ }\r
+ unsigned uSymbols1;\r
+ unsigned uSymbols2;\r
+ unsigned uSymbolsp;\r
+ unsigned uIndels1;\r
+ unsigned uIndels2;\r
+ unsigned uIndelsp;\r
+ EstringCounts(es1, &uSymbols1, &uIndels1);\r
+ EstringCounts(es2, &uSymbols2, &uIndels2);\r
+ EstringCounts(esp, &uSymbolsp, &uIndelsp);\r
+ if (uSymbols1 + uIndels1 != uSymbols2)\r
+ {\r
+ Log("Bad result of MulEstring: ");\r
+ LogEstring(esp);\r
+ Quit("Assert failed (counts1 %u %u %u)",\r
+ uSymbols1, uIndels1, uSymbols2);\r
+ }\r
+ }\r
+#endif\r
+ }\r
+\r
+static void test(const short es1[], const short es2[], const short esa[])\r
+ {\r
+ unsigned uSymbols1;\r
+ unsigned uSymbols2;\r
+ unsigned uIndels1;\r
+ unsigned uIndels2;\r
+ EstringCounts(es1, &uSymbols1, &uIndels1);\r
+ EstringCounts(es2, &uSymbols2, &uIndels2);\r
+\r
+ char s[4096];\r
+ memset(s, 'X', sizeof(s));\r
+ s[uSymbols1] = 0;\r
+\r
+ char *s1 = EstringOp(es1, s);\r
+ char *s12 = EstringOp(es2, s1);\r
+\r
+ memset(s, 'X', sizeof(s));\r
+ s[uSymbols2] = 0;\r
+ char *s2 = EstringOp(es2, s);\r
+\r
+ Log("%s * %s = %s\n", s1, s2, s12);\r
+\r
+ LogEstring(es1);\r
+ Log(" * ");\r
+ LogEstring(es2);\r
+ Log(" = ");\r
+ LogEstring(esa);\r
+ Log("\n");\r
+\r
+ short esp[4096];\r
+ MulEstrings(es1, es2, esp);\r
+ LogEstring(esp);\r
+ if (!EstringsEq(esp, esa))\r
+ Log(" *ERROR* ");\r
+ Log("\n");\r
+\r
+ memset(s, 'X', sizeof(s));\r
+ s[uSymbols1] = 0;\r
+ char *sp = EstringOp(esp, s);\r
+ Log("%s\n", sp);\r
+ Log("\n==========\n\n");\r
+ }\r
+\r
+void TestEstrings()\r
+ {\r
+ SetListFileName("c:\\tmp\\muscle.log", false);\r
+ //{\r
+ //short es1[] = { -1, 1, -1, 0 };\r
+ //short es2[] = { 1, -1, 2, 0 };\r
+ //short esa[] = { -2, 1, -1, 0 };\r
+ //test(es1, es2, esa);\r
+ //}\r
+ //{\r
+ //short es1[] = { 2, -1, 2, 0 };\r
+ //short es2[] = { 1, -1, 3, -1, 1, 0 };\r
+ //short esa[] = { 1, -1, 1, -1, 1, -1, 1, 0 };\r
+ //test(es1, es2, esa);\r
+ //}\r
+ //{\r
+ //short es1[] = { -1, 3, 0 };\r
+ //short es2[] = { 2, -1, 2, 0 };\r
+ //short esa[] = { -1, 1, -1, 2, 0 };\r
+ //test(es1, es2, esa);\r
+ //}\r
+ //{\r
+ //short es1[] = { -1, 1, -1, 1, 0};\r
+ //short es2[] = { 4, 0 };\r
+ //short esa[] = { -1, 1, -1, 1, 0};\r
+ //test(es1, es2, esa);\r
+ //}\r
+ //{\r
+ //short es1[] = { 1, -1, 1, -1, 0};\r
+ //short es2[] = { 4, 0 };\r
+ //short esa[] = { 1, -1, 1, -1, 0};\r
+ //test(es1, es2, esa);\r
+ //}\r
+ //{\r
+ //short es1[] = { 1, -1, 1, -1, 0};\r
+ //short es2[] = { -1, 4, -1, 0 };\r
+ //short esa[] = { -1, 1, -1, 1, -2, 0};\r
+ //test(es1, es2, esa);\r
+ //}\r
+ {\r
+ short es1[] = { 106, -77, 56, -2, 155, -3, 123, -2, 0};\r
+ short es2[] = { 50, -36, 34, -3, 12, -6, 1, -6, 18, -17, 60, -5, 349, -56, 0 };\r
+ short esa[] = { 0 };\r
+ test(es1, es2, esa);\r
+ }\r
+ exit(0);\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff