Next version of JABA
[jabaws.git] / binaries / src / muscle / estring.cpp
1 #include "muscle.h"\r
2 #include "pwpath.h"\r
3 #include "estring.h"\r
4 #include "seq.h"\r
5 #include "msa.h"\r
6 \r
7 /***\r
8 An "estring" is an edit string that operates on a sequence.\r
9 An estring is represented as a vector of integers.\r
10 It is interpreted in order of increasing suffix.\r
11 A positive value n means copy n letters.\r
12 A negative value -n means insert n indels.\r
13 Zero marks the end of the vector.\r
14 Consecutive entries must have opposite sign, i.e. the\r
15 shortest possible representation must be used.\r
16 \r
17 A "tpair" is a traceback path for a pairwise alignment\r
18 represented as two estrings, one for each sequence.\r
19 ***/\r
20 \r
21 #define c2(c,d) (((unsigned char) c) << 8 | (unsigned char) d)\r
22 \r
23 unsigned LengthEstring(const short es[])\r
24         {\r
25         unsigned i = 0;\r
26         while (*es++ != 0)\r
27                 ++i;\r
28         return i;\r
29         }\r
30 \r
31 short *EstringNewCopy(const short es[])\r
32         {\r
33         unsigned n = LengthEstring(es) + 1;\r
34         short *esNew = new short[n];\r
35         memcpy(esNew, es, n*sizeof(short));\r
36         return esNew;\r
37         }\r
38 \r
39 void LogEstring(const short es[])\r
40         {\r
41         Log("<");\r
42         for (unsigned i = 0; es[i] != 0; ++i)\r
43                 {\r
44                 if (i > 0)\r
45                         Log(" ");\r
46                 Log("%d", es[i]);\r
47                 }\r
48         Log(">");\r
49         }\r
50 \r
51 static bool EstringsEq(const short es1[], const short es2[])\r
52         {\r
53         for (;;)\r
54                 {\r
55                 if (*es1 != *es2)\r
56                         return false;\r
57                 if (0 == *es1)\r
58                         break;\r
59                 ++es1;\r
60                 ++es2;\r
61                 }\r
62         return true;\r
63         }\r
64 \r
65 static void EstringCounts(const short es[], unsigned *ptruSymbols,\r
66   unsigned *ptruIndels)\r
67         {\r
68         unsigned uSymbols = 0;\r
69         unsigned uIndels = 0;\r
70         for (unsigned i = 0; es[i] != 0; ++i)\r
71                 {\r
72                 short n = es[i];\r
73                 if (n > 0)\r
74                         uSymbols += n;\r
75                 else if (n < 0)\r
76                         uIndels += -n;\r
77                 }\r
78         *ptruSymbols = uSymbols;\r
79         *ptruIndels = uIndels;\r
80         }\r
81 \r
82 static char *EstringOp(const short es[], const char s[])\r
83         {\r
84         unsigned uSymbols;\r
85         unsigned uIndels;\r
86         EstringCounts(es, &uSymbols, &uIndels);\r
87         assert((unsigned) strlen(s) == uSymbols);\r
88         char *sout = new char[uSymbols + uIndels + 1];\r
89         char *psout = sout;\r
90         for (;;)\r
91                 {\r
92                 int n = *es++;\r
93                 if (0 == n)\r
94                         break;\r
95                 if (n > 0)\r
96                         for (int i = 0; i < n; ++i)\r
97                                 *psout++ = *s++;\r
98                 else\r
99                         for (int i = 0; i < -n; ++i)\r
100                                 *psout++ = '-';\r
101                 }\r
102         assert(0 == *s);\r
103         *psout = 0;\r
104         return sout;\r
105         }\r
106 \r
107 void EstringOp(const short es[], const Seq &sIn, Seq &sOut)\r
108         {\r
109 #if     DEBUG\r
110         unsigned uSymbols;\r
111         unsigned uIndels;\r
112         EstringCounts(es, &uSymbols, &uIndels);\r
113         assert(sIn.Length() == uSymbols);\r
114 #endif\r
115         sOut.Clear();\r
116         sOut.SetName(sIn.GetName());\r
117         int p = 0;\r
118         for (;;)\r
119                 {\r
120                 int n = *es++;\r
121                 if (0 == n)\r
122                         break;\r
123                 if (n > 0)\r
124                         for (int i = 0; i < n; ++i)\r
125                                 {\r
126                                 char c = sIn[p++];\r
127                                 sOut.push_back(c);\r
128                                 }\r
129                 else\r
130                         for (int i = 0; i < -n; ++i)\r
131                                 sOut.push_back('-');\r
132                 }\r
133         }\r
134 \r
135 unsigned EstringOp(const short es[], const Seq &sIn, MSA &a)\r
136         {\r
137         unsigned uSymbols;\r
138         unsigned uIndels;\r
139         EstringCounts(es, &uSymbols, &uIndels);\r
140         assert(sIn.Length() == uSymbols);\r
141 \r
142         unsigned uColCount = uSymbols + uIndels;\r
143 \r
144         a.Clear();\r
145         a.SetSize(1, uColCount);\r
146 \r
147         a.SetSeqName(0, sIn.GetName());\r
148         a.SetSeqId(0, sIn.GetId());\r
149 \r
150         unsigned p = 0;\r
151         unsigned uColIndex = 0;\r
152         for (;;)\r
153                 {\r
154                 int n = *es++;\r
155                 if (0 == n)\r
156                         break;\r
157                 if (n > 0)\r
158                         for (int i = 0; i < n; ++i)\r
159                                 {\r
160                                 char c = sIn[p++];\r
161                                 a.SetChar(0, uColIndex++, c);\r
162                                 }\r
163                 else\r
164                         for (int i = 0; i < -n; ++i)\r
165                                 a.SetChar(0, uColIndex++, '-');\r
166                 }\r
167         assert(uColIndex == uColCount);\r
168         return uColCount;\r
169         }\r
170 \r
171 void PathToEstrings(const PWPath &Path, short **ptresA, short **ptresB)\r
172         {\r
173 // First pass to determine size of estrings esA and esB\r
174         const unsigned uEdgeCount = Path.GetEdgeCount();\r
175         if (0 == uEdgeCount)\r
176                 {\r
177                 short *esA = new short[1];\r
178                 short *esB = new short[1];\r
179                 esA[0] = 0;\r
180                 esB[0] = 0;\r
181                 *ptresA = esA;\r
182                 *ptresB = esB;\r
183                 return;\r
184                 }\r
185 \r
186         unsigned iLengthA = 1;\r
187         unsigned iLengthB = 1;\r
188         const char cFirstEdgeType = Path.GetEdge(0).cType;\r
189         char cPrevEdgeType = cFirstEdgeType;\r
190         for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
191                 {\r
192                 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
193                 char cEdgeType = Edge.cType;\r
194 \r
195                 switch (c2(cPrevEdgeType, cEdgeType))\r
196                         {\r
197                 case c2('M', 'M'):\r
198                 case c2('D', 'D'):\r
199                 case c2('I', 'I'):\r
200                         break;\r
201 \r
202                 case c2('D', 'M'):\r
203                 case c2('M', 'D'):\r
204                         ++iLengthB;\r
205                         break;\r
206 \r
207                 case c2('I', 'M'):\r
208                 case c2('M', 'I'):\r
209                         ++iLengthA;\r
210                         break;\r
211 \r
212                 case c2('I', 'D'):\r
213                 case c2('D', 'I'):\r
214                         ++iLengthB;\r
215                         ++iLengthA;\r
216                         break;\r
217 \r
218                 default:\r
219                         assert(false);\r
220                         }\r
221                 cPrevEdgeType = cEdgeType;\r
222                 }\r
223 \r
224 // Pass2 for seq A\r
225         {\r
226         short *esA = new short[iLengthA+1];\r
227         unsigned iA = 0;\r
228         switch (Path.GetEdge(0).cType)\r
229                 {\r
230         case 'M':\r
231         case 'D':\r
232                 esA[0] = 1;\r
233                 break;\r
234 \r
235         case 'I':\r
236                 esA[0] = -1;\r
237                 break;\r
238 \r
239         default:\r
240                 assert(false);\r
241                 }\r
242 \r
243         char cPrevEdgeType = cFirstEdgeType;\r
244         for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
245                 {\r
246                 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
247                 char cEdgeType = Edge.cType;\r
248 \r
249                 switch (c2(cPrevEdgeType, cEdgeType))\r
250                         {\r
251                 case c2('M', 'M'):\r
252                 case c2('D', 'D'):\r
253                 case c2('D', 'M'):\r
254                 case c2('M', 'D'):\r
255                         ++(esA[iA]);\r
256                         break;\r
257 \r
258                 case c2('I', 'D'):\r
259                 case c2('I', 'M'):\r
260                         ++iA;\r
261                         esA[iA] = 1;\r
262                         break;\r
263 \r
264                 case c2('M', 'I'):\r
265                 case c2('D', 'I'):\r
266                         ++iA;\r
267                         esA[iA] = -1;\r
268                         break;\r
269 \r
270                 case c2('I', 'I'):\r
271                         --(esA[iA]);\r
272                         break;\r
273 \r
274                 default:\r
275                         assert(false);\r
276                         }\r
277 \r
278                 cPrevEdgeType = cEdgeType;\r
279                 }\r
280         assert(iA == iLengthA - 1);\r
281         esA[iLengthA] = 0;\r
282         *ptresA = esA;\r
283         }\r
284 \r
285         {\r
286 // Pass2 for seq B\r
287         short *esB = new short[iLengthB+1];\r
288         unsigned iB = 0;\r
289         switch (Path.GetEdge(0).cType)\r
290                 {\r
291         case 'M':\r
292         case 'I':\r
293                 esB[0] = 1;\r
294                 break;\r
295 \r
296         case 'D':\r
297                 esB[0] = -1;\r
298                 break;\r
299 \r
300         default:\r
301                 assert(false);\r
302                 }\r
303 \r
304         char cPrevEdgeType = cFirstEdgeType;\r
305         for (unsigned uEdgeIndex = 1; uEdgeIndex < uEdgeCount; ++uEdgeIndex)\r
306                 {\r
307                 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);\r
308                 char cEdgeType = Edge.cType;\r
309 \r
310                 switch (c2(cPrevEdgeType, cEdgeType))\r
311                         {\r
312                 case c2('M', 'M'):\r
313                 case c2('I', 'I'):\r
314                 case c2('I', 'M'):\r
315                 case c2('M', 'I'):\r
316                         ++(esB[iB]);\r
317                         break;\r
318 \r
319                 case c2('D', 'I'):\r
320                 case c2('D', 'M'):\r
321                         ++iB;\r
322                         esB[iB] = 1;\r
323                         break;\r
324 \r
325                 case c2('M', 'D'):\r
326                 case c2('I', 'D'):\r
327                         ++iB;\r
328                         esB[iB] = -1;\r
329                         break;\r
330 \r
331                 case c2('D', 'D'):\r
332                         --(esB[iB]);\r
333                         break;\r
334 \r
335                 default:\r
336                         assert(false);\r
337                         }\r
338 \r
339                 cPrevEdgeType = cEdgeType;\r
340                 }\r
341         assert(iB == iLengthB - 1);\r
342         esB[iLengthB] = 0;\r
343         *ptresB = esB;\r
344         }\r
345 \r
346 #if     DEBUG\r
347         {\r
348         const PWEdge &LastEdge = Path.GetEdge(uEdgeCount - 1);\r
349         unsigned uSymbols;\r
350         unsigned uIndels;\r
351         EstringCounts(*ptresA, &uSymbols, &uIndels);\r
352         assert(uSymbols == LastEdge.uPrefixLengthA);\r
353         assert(uSymbols + uIndels == uEdgeCount);\r
354 \r
355         EstringCounts(*ptresB, &uSymbols, &uIndels);\r
356         assert(uSymbols == LastEdge.uPrefixLengthB);\r
357         assert(uSymbols + uIndels == uEdgeCount);\r
358 \r
359         PWPath TmpPath;\r
360         EstringsToPath(*ptresA, *ptresB, TmpPath);\r
361         TmpPath.AssertEqual(Path);\r
362         }\r
363 #endif\r
364         }\r
365 \r
366 void EstringsToPath(const short esA[], const short esB[], PWPath &Path)\r
367         {\r
368         Path.Clear();\r
369         unsigned iA = 0;\r
370         unsigned iB = 0;\r
371         int nA = esA[iA++];\r
372         int nB = esB[iB++];\r
373         unsigned uPrefixLengthA = 0;\r
374         unsigned uPrefixLengthB = 0;\r
375         for (;;)\r
376                 {\r
377                 char cType;\r
378                 if (nA > 0)\r
379                         {\r
380                         if (nB > 0)\r
381                                 {\r
382                                 cType = 'M';\r
383                                 --nA;\r
384                                 --nB;\r
385                                 }\r
386                         else if (nB < 0)\r
387                                 {\r
388                                 cType = 'D';\r
389                                 --nA;\r
390                                 ++nB;\r
391                                 }\r
392                         else\r
393                                 assert(false);\r
394                         }\r
395                 else if (nA < 0)\r
396                         {\r
397                         if (nB > 0)\r
398                                 {\r
399                                 cType = 'I';\r
400                                 ++nA;\r
401                                 --nB;\r
402                                 }\r
403                         else\r
404                                 assert(false);\r
405                         }\r
406                 else\r
407                         assert(false);\r
408 \r
409                 switch (cType)\r
410                         {\r
411                 case 'M':\r
412                         ++uPrefixLengthA;\r
413                         ++uPrefixLengthB;\r
414                         break;\r
415                 case 'D':\r
416                         ++uPrefixLengthA;\r
417                         break;\r
418                 case 'I':\r
419                         ++uPrefixLengthB;\r
420                         break;\r
421                         }\r
422 \r
423                 PWEdge Edge;\r
424                 Edge.cType = cType;\r
425                 Edge.uPrefixLengthA = uPrefixLengthA;\r
426                 Edge.uPrefixLengthB = uPrefixLengthB;\r
427                 Path.AppendEdge(Edge);\r
428 \r
429                 if (nA == 0)\r
430                         {\r
431                         if (0 == esA[iA])\r
432                                 {\r
433                                 assert(0 == esB[iB]);\r
434                                 break;\r
435                                 }\r
436                         nA = esA[iA++];\r
437                         }\r
438                 if (nB == 0)\r
439                         nB = esB[iB++];\r
440                 }\r
441         }\r
442 \r
443 /***\r
444 Multiply two estrings to make a third estring.\r
445 The product of two estrings e1*e2 is defined to be\r
446 the estring that produces the same result as applying\r
447 e1 then e2. Multiplication is not commutative. In fact,\r
448 the reversed order is undefined unless both estrings\r
449 consist of a single, identical, positive entry.\r
450 A primary motivation for using estrings is that\r
451 multiplication is very fast, reducing the time\r
452 needed to construct the root alignment.\r
453 \r
454 Example\r
455 \r
456         <-1,3>(XXX)     = -XXX\r
457         <2,-1,2>(-XXX) = -X-XX\r
458 \r
459 Therefore,\r
460 \r
461         <-1,3>*<2,-1,2> = <-1,1,-1,2>\r
462 ***/\r
463 \r
464 static bool CanMultiplyEstrings(const short es1[], const short es2[])\r
465         {\r
466         unsigned uSymbols1;\r
467         unsigned uSymbols2;\r
468         unsigned uIndels1;\r
469         unsigned uIndels2;\r
470         EstringCounts(es1, &uSymbols1, &uIndels1);\r
471         EstringCounts(es2, &uSymbols2, &uIndels2);\r
472         return uSymbols1 + uIndels1 == uSymbols2;\r
473         }\r
474 \r
475 static inline void AppendGaps(short esp[], int &ip, int n)\r
476         {\r
477         if (-1 == ip)\r
478                 esp[++ip] = n;\r
479         else if (esp[ip] < 0)\r
480                 esp[ip] += n;\r
481         else\r
482                 esp[++ip] = n;\r
483         }\r
484 \r
485 static inline void AppendSymbols(short esp[], int &ip, int n)\r
486         {\r
487         if (-1 == ip)\r
488                 esp[++ip] = n;\r
489         else if (esp[ip] > 0)\r
490                 esp[ip] += n;\r
491         else\r
492                 esp[++ip] = n;\r
493         }\r
494 \r
495 void MulEstrings(const short es1[], const short es2[], short esp[])\r
496         {\r
497         assert(CanMultiplyEstrings(es1, es2));\r
498 \r
499         unsigned i1 = 0;\r
500         int ip = -1;\r
501         int n1 = es1[i1++];\r
502         for (unsigned i2 = 0; ; ++i2)\r
503                 {\r
504                 int n2 = es2[i2];\r
505                 if (0 == n2)\r
506                         break;\r
507                 if (n2 > 0)\r
508                         {\r
509                         for (;;)\r
510                                 {\r
511                                 if (n1 < 0)\r
512                                         {\r
513                                         if (n2 > -n1)\r
514                                                 {\r
515                                                 AppendGaps(esp, ip, n1);\r
516                                                 n2 += n1;\r
517                                                 n1 = es1[i1++];\r
518                                                 }\r
519                                         else if (n2 == -n1)\r
520                                                 {\r
521                                                 AppendGaps(esp, ip, n1);\r
522                                                 n1 = es1[i1++];\r
523                                                 break;\r
524                                                 }\r
525                                         else\r
526                                                 {\r
527                                                 assert(n2 < -n1);\r
528                                                 AppendGaps(esp, ip, -n2);\r
529                                                 n1 += n2;\r
530                                                 break;\r
531                                                 }\r
532                                         }\r
533                                 else\r
534                                         {\r
535                                         assert(n1 > 0);\r
536                                         if (n2 > n1)\r
537                                                 {\r
538                                                 AppendSymbols(esp, ip, n1);\r
539                                                 n2 -= n1;\r
540                                                 n1 = es1[i1++];\r
541                                                 }\r
542                                         else if (n2 == n1)\r
543                                                 {\r
544                                                 AppendSymbols(esp, ip, n1);\r
545                                                 n1 = es1[i1++];\r
546                                                 break;\r
547                                                 }\r
548                                         else\r
549                                                 {\r
550                                                 assert(n2 < n1);\r
551                                                 AppendSymbols(esp, ip, n2);\r
552                                                 n1 -= n2;\r
553                                                 break;\r
554                                                 }\r
555                                         }\r
556                                 }\r
557                         }\r
558                 else\r
559                         {\r
560                         assert(n2 < 0);\r
561                         AppendGaps(esp, ip, n2);\r
562                         }\r
563                 }\r
564         esp[++ip] = 0;\r
565 \r
566 #if     DEBUG\r
567         {\r
568         int MaxLen = (int) (LengthEstring(es1) + LengthEstring(es2) + 1);\r
569         assert(ip < MaxLen);\r
570         if (ip >= 2)\r
571                 for (int i = 0; i < ip - 2; ++i)\r
572                         {\r
573                         if (!(esp[i] > 0 && esp[i+1] < 0 || esp[i] < 0 && esp[i+1] > 0))\r
574                                 {\r
575                                 Log("Bad result of MulEstring: ");\r
576                                 LogEstring(esp);\r
577                                 Quit("Assert failed (alternating signs)");\r
578                                 }\r
579                         }\r
580         unsigned uSymbols1;\r
581         unsigned uSymbols2;\r
582         unsigned uSymbolsp;\r
583         unsigned uIndels1;\r
584         unsigned uIndels2;\r
585         unsigned uIndelsp;\r
586         EstringCounts(es1, &uSymbols1, &uIndels1);\r
587         EstringCounts(es2, &uSymbols2, &uIndels2);\r
588         EstringCounts(esp, &uSymbolsp, &uIndelsp);\r
589         if (uSymbols1 + uIndels1 != uSymbols2)\r
590                 {\r
591                 Log("Bad result of MulEstring: ");\r
592                 LogEstring(esp);\r
593                 Quit("Assert failed (counts1 %u %u %u)",\r
594                   uSymbols1, uIndels1, uSymbols2);\r
595                 }\r
596         }\r
597 #endif\r
598         }\r
599 \r
600 static void test(const short es1[], const short es2[], const short esa[])\r
601         {\r
602         unsigned uSymbols1;\r
603         unsigned uSymbols2;\r
604         unsigned uIndels1;\r
605         unsigned uIndels2;\r
606         EstringCounts(es1, &uSymbols1, &uIndels1);\r
607         EstringCounts(es2, &uSymbols2, &uIndels2);\r
608 \r
609         char s[4096];\r
610         memset(s, 'X', sizeof(s));\r
611         s[uSymbols1] = 0;\r
612 \r
613         char *s1 = EstringOp(es1, s);\r
614         char *s12 = EstringOp(es2, s1);\r
615 \r
616         memset(s, 'X', sizeof(s));\r
617         s[uSymbols2] = 0;\r
618         char *s2 = EstringOp(es2, s);\r
619 \r
620         Log("%s * %s = %s\n", s1, s2, s12);\r
621 \r
622         LogEstring(es1);\r
623         Log(" * ");\r
624         LogEstring(es2);\r
625         Log(" = ");\r
626         LogEstring(esa);\r
627         Log("\n");\r
628 \r
629         short esp[4096];\r
630         MulEstrings(es1, es2, esp);\r
631         LogEstring(esp);\r
632         if (!EstringsEq(esp, esa))\r
633                 Log(" *ERROR* ");\r
634         Log("\n");\r
635 \r
636         memset(s, 'X', sizeof(s));\r
637         s[uSymbols1] = 0;\r
638         char *sp = EstringOp(esp, s);\r
639         Log("%s\n", sp);\r
640         Log("\n==========\n\n");\r
641         }\r
642 \r
643 void TestEstrings()\r
644         {\r
645         SetListFileName("c:\\tmp\\muscle.log", false);\r
646         //{\r
647         //short es1[] = { -1, 1, -1, 0 };\r
648         //short es2[] = { 1, -1, 2, 0 };\r
649         //short esa[] = { -2, 1, -1, 0 };\r
650         //test(es1, es2, esa);\r
651         //}\r
652         //{\r
653         //short es1[] = { 2, -1, 2, 0 };\r
654         //short es2[] = { 1, -1, 3, -1, 1, 0 };\r
655         //short esa[] = { 1, -1, 1, -1, 1, -1, 1, 0 };\r
656         //test(es1, es2, esa);\r
657         //}\r
658         //{\r
659         //short es1[] = { -1, 3, 0 };\r
660         //short es2[] = { 2, -1, 2, 0 };\r
661         //short esa[] = { -1, 1, -1, 2, 0 };\r
662         //test(es1, es2, esa);\r
663         //}\r
664         //{\r
665         //short es1[] = { -1, 1, -1, 1, 0};\r
666         //short es2[] = { 4, 0 };\r
667         //short esa[] = { -1, 1, -1, 1, 0};\r
668         //test(es1, es2, esa);\r
669         //}\r
670         //{\r
671         //short es1[] = { 1, -1, 1, -1, 0};\r
672         //short es2[] = { 4, 0 };\r
673         //short esa[] = { 1, -1, 1, -1, 0};\r
674         //test(es1, es2, esa);\r
675         //}\r
676         //{\r
677         //short es1[] = { 1, -1, 1, -1, 0};\r
678         //short es2[] = { -1, 4, -1, 0 };\r
679         //short esa[] = { -1, 1, -1, 1, -2, 0};\r
680         //test(es1, es2, esa);\r
681         //}\r
682         {\r
683         short es1[] = { 106, -77, 56, -2, 155, -3, 123, -2, 0};\r
684         short es2[] = { 50, -36, 34, -3, 12, -6, 1, -6, 18, -17, 60, -5, 349, -56, 0 };\r
685         short esa[] = { 0 };\r
686         test(es1, es2, esa);\r
687         }\r
688         exit(0);\r
689         }\r