--- /dev/null
+#include "muscle.h"\r
+#include "profile.h"\r
+#include "pwpath.h"\r
+#include <math.h>\r
+\r
+#define TRACE 0\r
+\r
+#define EQ(a, b) (fabs(a-b) < 0.1)\r
+\r
+SCORE TraceBack(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
+ unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,\r
+ PWPath &Path)\r
+ {\r
+#if TRACE\r
+ Log("\n");\r
+ Log("TraceBack LengthA=%u LengthB=%u\n", uLengthA, uLengthB);\r
+#endif\r
+ assert(uLengthB > 0 && uLengthA > 0);\r
+\r
+ const unsigned uPrefixCountA = uLengthA + 1;\r
+ const unsigned uPrefixCountB = uLengthB + 1;\r
+\r
+ Path.Clear();\r
+\r
+ unsigned uPrefixLengthA = uLengthA;\r
+ unsigned uPrefixLengthB = uLengthB;\r
+\r
+ const SCORE scoreM = DPM(uPrefixLengthA, uPrefixLengthB);\r
+ SCORE scoreD = DPD(uPrefixLengthA, uPrefixLengthB);\r
+ SCORE scoreI = DPI(uPrefixLengthA, uPrefixLengthB);\r
+\r
+ const ProfPos &LastPPA = PA[uLengthA - 1];\r
+ const ProfPos &LastPPB = PB[uLengthB - 1];\r
+\r
+ scoreD += LastPPA.m_scoreGapClose;\r
+ scoreI += LastPPB.m_scoreGapClose;\r
+\r
+ char cEdgeType = cInsane;\r
+ SCORE scoreMax;\r
+ if (scoreM >= scoreD && scoreM >= scoreI)\r
+ {\r
+ scoreMax = scoreM;\r
+ cEdgeType = 'M';\r
+ }\r
+ else if (scoreD >= scoreM && scoreD >= scoreI)\r
+ {\r
+ scoreMax = scoreD;\r
+ cEdgeType = 'D';\r
+ }\r
+ else\r
+ {\r
+ assert(scoreI >= scoreM && scoreI >= scoreD);\r
+ scoreMax = scoreI;\r
+ cEdgeType = 'I';\r
+ }\r
+\r
+ for (;;)\r
+ {\r
+ if ('S' == cEdgeType)\r
+ break;\r
+\r
+ PWEdge Edge;\r
+ Edge.cType = cEdgeType;\r
+ Edge.uPrefixLengthA = uPrefixLengthA;\r
+ Edge.uPrefixLengthB = uPrefixLengthB;\r
+ Path.PrependEdge(Edge);\r
+\r
+ char cPrevEdgeType;\r
+ unsigned uPrevPrefixLengthA = uPrefixLengthA;\r
+ unsigned uPrevPrefixLengthB = uPrefixLengthB;\r
+\r
+ switch (cEdgeType)\r
+ {\r
+ case 'M':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ assert(uPrefixLengthB > 0);\r
+ const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
+ const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
+\r
+ const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);\r
+ const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);\r
+\r
+ SCORE scoreSM;\r
+ if (1 == uPrefixLengthA && 1 == uPrefixLengthB)\r
+ scoreSM = scoreMatch;\r
+ else\r
+ scoreSM = MINUS_INFINITY;\r
+\r
+ SCORE scoreMM = MINUS_INFINITY;\r
+ SCORE scoreDM = MINUS_INFINITY;\r
+ SCORE scoreIM = MINUS_INFINITY;\r
+ if (uPrefixLengthA > 1 && uPrefixLengthB > 1)\r
+ scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1) + scoreMatch;\r
+ if (uPrefixLengthA > 1)\r
+ {\r
+ SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;\r
+ scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;\r
+ }\r
+ if (uPrefixLengthB > 1)\r
+ {\r
+ SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;\r
+ scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;\r
+ }\r
+\r
+ if (EQ(scoreMM, Score))\r
+ cPrevEdgeType = 'M';\r
+ else if (EQ(scoreDM, Score))\r
+ cPrevEdgeType = 'D';\r
+ else if (EQ(scoreIM, Score))\r
+ cPrevEdgeType = 'I';\r
+ else if (EQ(scoreSM, Score))\r
+ cPrevEdgeType = 'S';\r
+ else\r
+ Quit("TraceBack: failed to match M score=%g M=%g D=%g I=%g S=%g",\r
+ Score, scoreMM, scoreDM, scoreIM, scoreSM);\r
+\r
+ --uPrevPrefixLengthA;\r
+ --uPrevPrefixLengthB;\r
+ break;\r
+ }\r
+\r
+ case 'D':\r
+ {\r
+ assert(uPrefixLengthA > 0);\r
+ const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);\r
+\r
+ SCORE scoreMD = MINUS_INFINITY;\r
+ SCORE scoreDD = MINUS_INFINITY;\r
+ SCORE scoreSD = MINUS_INFINITY;\r
+ if (uPrefixLengthB == 0)\r
+ {\r
+ if (uPrefixLengthA == 1)\r
+ scoreSD = PA[0].m_scoreGapOpen;\r
+ else\r
+ scoreSD = DPD(uPrefixLengthA - 1, 0);\r
+ }\r
+ if (uPrefixLengthA > 1)\r
+ {\r
+ const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
+ SCORE scoreTransMD = PPA.m_scoreGapOpen;\r
+ scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;\r
+ scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);\r
+ }\r
+\r
+ if (EQ(Score, scoreMD))\r
+ cPrevEdgeType = 'M';\r
+ else if (EQ(Score, scoreDD))\r
+ cPrevEdgeType = 'D';\r
+ else if (EQ(Score, scoreSD))\r
+ cPrevEdgeType = 'S';\r
+ else\r
+ Quit("TraceBack: failed to match D");\r
+\r
+ --uPrevPrefixLengthA;\r
+ break;\r
+ }\r
+\r
+ case 'I':\r
+ {\r
+ assert(uPrefixLengthB > 0);\r
+ const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);\r
+\r
+ SCORE scoreMI = MINUS_INFINITY;\r
+ SCORE scoreII = MINUS_INFINITY;\r
+ SCORE scoreSI = MINUS_INFINITY;\r
+ if (uPrefixLengthA == 0)\r
+ {\r
+ if (uPrefixLengthB == 1)\r
+ scoreSI = PB[0].m_scoreGapOpen;\r
+ else\r
+ scoreSI = DPI(0, uPrefixLengthB - 1);\r
+ }\r
+ if (uPrefixLengthB > 1)\r
+ {\r
+ const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
+ SCORE scoreTransMI = PPB.m_scoreGapOpen;\r
+ scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;\r
+ scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);\r
+ }\r
+\r
+ if (EQ(Score, scoreMI))\r
+ cPrevEdgeType = 'M';\r
+ else if (EQ(Score, scoreII))\r
+ cPrevEdgeType = 'I';\r
+ else if (EQ(Score, scoreSI))\r
+ cPrevEdgeType = 'S';\r
+ else\r
+ Quit("TraceBack: failed to match I");\r
+\r
+ --uPrevPrefixLengthB;\r
+ break;\r
+ }\r
+\r
+ default:\r
+ assert(false);\r
+ }\r
+#if TRACE\r
+ Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);\r
+ Log("\n");\r
+#endif\r
+ cEdgeType = cPrevEdgeType;\r
+ uPrefixLengthA = uPrevPrefixLengthA;\r
+ uPrefixLengthB = uPrevPrefixLengthB;\r
+ }\r
+\r
+ return scoreMax;\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff