--- /dev/null
+#include "muscle.h"\r
+#include <math.h>\r
+#include "pwpath.h"\r
+#include "profile.h"\r
+#include <stdio.h>\r
+\r
+#define TRACE 0\r
+\r
+bool g_bKeepSimpleDP;\r
+SCORE *g_DPM;\r
+SCORE *g_DPD;\r
+SCORE *g_DPE;\r
+SCORE *g_DPI;\r
+SCORE *g_DPJ;\r
+char *g_TBM;\r
+char *g_TBD;\r
+char *g_TBE;\r
+char *g_TBI;\r
+char *g_TBJ;\r
+\r
+#if DOUBLE_AFFINE\r
+\r
+static char XlatEdgeType(char c)\r
+ {\r
+ if ('E' == c)\r
+ return 'D';\r
+ if ('J' == c)\r
+ return 'I';\r
+ return c;\r
+ }\r
+\r
+static const char *LocalScoreToStr(SCORE s)\r
+ {\r
+ static char str[16];\r
+ if (s < -100000)\r
+ return " *";\r
+ sprintf(str, "%6.1f", s);\r
+ return str;\r
+ }\r
+\r
+static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,\r
+ unsigned uPrefixCountA, unsigned uPrefixCountB)\r
+ {\r
+ Log(" ");\r
+ for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ {\r
+ char c = ' ';\r
+ if (uPrefixLengthB > 0)\r
+ c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
+ Log(" %4u:%c", uPrefixLengthB, c);\r
+ }\r
+ Log("\n");\r
+ for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
+ {\r
+ char c = ' ';\r
+ if (uPrefixLengthA > 0)\r
+ c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
+ Log("%4u:%c ", uPrefixLengthA, c);\r
+ for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));\r
+ Log("\n");\r
+ }\r
+ }\r
+\r
+static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
+ unsigned uPrefixCountA, unsigned uPrefixCountB)\r
+ {\r
+ Log(" ");\r
+ for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ {\r
+ char c = ' ';\r
+ if (uPrefixLengthB > 0)\r
+ c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
+ Log(" %4u:%c", uPrefixLengthB, c);\r
+ }\r
+ Log("\n");\r
+ for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
+ {\r
+ char c = ' ';\r
+ if (uPrefixLengthA > 0)\r
+ c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
+ Log("%4u:%c ", uPrefixLengthA, c);\r
+ for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));\r
+ Log("\n");\r
+ }\r
+ }\r
+\r
+SCORE NWDASimple(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
+ unsigned uLengthB, PWPath &Path)\r
+ {\r
+ assert(uLengthB > 0 && uLengthA > 0);\r
+\r
+ const unsigned uPrefixCountA = uLengthA + 1;\r
+ const unsigned uPrefixCountB = uLengthB + 1;\r
+\r
+// Allocate DP matrices\r
+ const size_t LM = uPrefixCountA*uPrefixCountB;\r
+ SCORE *DPL_ = new SCORE[LM];\r
+ SCORE *DPM_ = new SCORE[LM];\r
+ SCORE *DPD_ = new SCORE[LM];\r
+ SCORE *DPE_ = new SCORE[LM];\r
+ SCORE *DPI_ = new SCORE[LM];\r
+ SCORE *DPJ_ = new SCORE[LM];\r
+\r
+ char *TBM_ = new char[LM];\r
+ char *TBD_ = new char[LM];\r
+ char *TBE_ = new char[LM];\r
+ char *TBI_ = new char[LM];\r
+ char *TBJ_ = new char[LM];\r
+\r
+ memset(TBM_, '?', LM);\r
+ memset(TBD_, '?', LM);\r
+ memset(TBE_, '?', LM);\r
+ memset(TBI_, '?', LM);\r
+ memset(TBJ_, '?', LM);\r
+\r
+ DPM(0, 0) = 0;\r
+ DPD(0, 0) = MINUS_INFINITY;\r
+ DPE(0, 0) = MINUS_INFINITY;\r
+ DPI(0, 0) = MINUS_INFINITY;\r
+ DPJ(0, 0) = MINUS_INFINITY;\r
+\r
+ DPM(1, 0) = MINUS_INFINITY;\r
+ DPD(1, 0) = PA[0].m_scoreGapOpen;\r
+ DPE(1, 0) = PA[0].m_scoreGapOpen2;\r
+ TBD(1, 0) = 'D';\r
+ TBE(1, 0) = 'E';\r
+ DPI(1, 0) = MINUS_INFINITY;\r
+ DPJ(1, 0) = MINUS_INFINITY;\r
+\r
+ DPM(0, 1) = MINUS_INFINITY;\r
+ DPD(0, 1) = MINUS_INFINITY;\r
+ DPE(0, 1) = MINUS_INFINITY;\r
+ DPI(0, 1) = PB[0].m_scoreGapOpen;\r
+ DPJ(0, 1) = PB[0].m_scoreGapOpen2;\r
+ TBI(0, 1) = 'I';\r
+ TBJ(0, 1) = 'J';\r
+\r
+// Empty prefix of B is special case\r
+ for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
+ {\r
+ DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+\r
+ DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;\r
+ DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;\r
+\r
+ TBD(uPrefixLengthA, 0) = 'D';\r
+ TBE(uPrefixLengthA, 0) = 'E';\r
+\r
+ DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+ DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+ }\r
+\r
+// Empty prefix of A is special case\r
+ for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ {\r
+ DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
+\r
+ DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
+ DPE(0, uPrefixLengthB) = MINUS_INFINITY;\r
+\r
+ DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;\r
+ DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;\r
+\r
+ TBI(0, uPrefixLengthB) = 'I';\r
+ TBJ(0, uPrefixLengthB) = 'J';\r
+ }\r
+\r
+// Special case to agree with NWFast, no D-I transitions so...\r
+ DPD(uLengthA, 0) = MINUS_INFINITY;\r
+ DPE(uLengthA, 0) = MINUS_INFINITY;\r
+// DPI(0, uLengthB) = MINUS_INFINITY;\r
+// DPJ(0, uLengthB) = MINUS_INFINITY;\r
+\r
+// ============\r
+// Main DP loop\r
+// ============\r
+ SCORE scoreGapCloseB = MINUS_INFINITY;\r
+ SCORE scoreGapClose2B = MINUS_INFINITY;\r
+ for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ {\r
+ const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
+\r
+ SCORE scoreGapCloseA = MINUS_INFINITY;\r
+ SCORE scoreGapClose2A = MINUS_INFINITY;\r
+ for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
+ {\r
+ const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
+\r
+ {\r
+ // Match M=LetterA+LetterB\r
+ SCORE scoreLL = ScoreProfPos2(PPA, PPB);\r
+ DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;\r
+\r
+ SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
+ SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
+ SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;\r
+ SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
+ SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMM >= scoreDM && scoreMM >= scoreEM && scoreMM >= scoreIM && scoreMM >= scoreJM)\r
+ {\r
+ scoreBest = scoreMM;\r
+ TBM(uPrefixLengthA, uPrefixLengthB) = 'M';\r
+ }\r
+ else if (scoreDM >= scoreMM && scoreDM >= scoreEM && scoreDM >= scoreIM && scoreDM >= scoreJM)\r
+ {\r
+ scoreBest = scoreDM;\r
+ TBM(uPrefixLengthA, uPrefixLengthB) = 'D';\r
+ }\r
+ else if (scoreEM >= scoreMM && scoreEM >= scoreDM && scoreEM >= scoreIM && scoreEM >= scoreJM)\r
+ {\r
+ scoreBest = scoreEM;\r
+ TBM(uPrefixLengthA, uPrefixLengthB) = 'E';\r
+ }\r
+ else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)\r
+ {\r
+ scoreBest = scoreIM;\r
+ TBM(uPrefixLengthA, uPrefixLengthB) = 'I';\r
+ }\r
+ else\r
+ {\r
+ assert(scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM);\r
+ scoreBest = scoreJM;\r
+ TBM(uPrefixLengthA, uPrefixLengthB) = 'J';\r
+ }\r
+ DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;\r
+ }\r
+\r
+ {\r
+ // Delete D=LetterA+GapB\r
+ SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
+ PA[uPrefixLengthA-1].m_scoreGapOpen;\r
+ SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMD >= scoreDD)\r
+ {\r
+ scoreBest = scoreMD;\r
+ TBD(uPrefixLengthA, uPrefixLengthB) = 'M';\r
+ }\r
+ else\r
+ {\r
+ assert(scoreDD >= scoreMD);\r
+ scoreBest = scoreDD;\r
+ TBD(uPrefixLengthA, uPrefixLengthB) = 'D';\r
+ }\r
+ DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
+ }\r
+\r
+ {\r
+ // Delete E=LetterA+GapB\r
+ SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
+ PA[uPrefixLengthA-1].m_scoreGapOpen2;\r
+ SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend2;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreME >= scoreEE)\r
+ {\r
+ scoreBest = scoreME;\r
+ TBE(uPrefixLengthA, uPrefixLengthB) = 'M';\r
+ }\r
+ else\r
+ {\r
+ assert(scoreEE >= scoreME);\r
+ scoreBest = scoreEE;\r
+ TBE(uPrefixLengthA, uPrefixLengthB) = 'E';\r
+ }\r
+ DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
+ }\r
+\r
+ // Insert I=GapA+LetterB\r
+ {\r
+ SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
+ PB[uPrefixLengthB - 1].m_scoreGapOpen;\r
+ SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMI >= scoreII)\r
+ {\r
+ scoreBest = scoreMI;\r
+ TBI(uPrefixLengthA, uPrefixLengthB) = 'M';\r
+ }\r
+ else \r
+ {\r
+ assert(scoreII > scoreMI);\r
+ scoreBest = scoreII;\r
+ TBI(uPrefixLengthA, uPrefixLengthB) = 'I';\r
+ }\r
+ DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
+ }\r
+\r
+ // Insert J=GapA+LetterB\r
+ {\r
+ SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
+ PB[uPrefixLengthB - 1].m_scoreGapOpen2;\r
+ SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend2;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMJ >= scoreJJ)\r
+ {\r
+ scoreBest = scoreMJ;\r
+ TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';\r
+ }\r
+ else \r
+ {\r
+ assert(scoreJJ > scoreMJ);\r
+ scoreBest = scoreJJ;\r
+ TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';\r
+ }\r
+ DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
+ }\r
+\r
+ scoreGapCloseA = PPA.m_scoreGapClose;\r
+ scoreGapClose2A = PPA.m_scoreGapClose2;\r
+ }\r
+ scoreGapCloseB = PPB.m_scoreGapClose;\r
+ scoreGapClose2B = PPB.m_scoreGapClose2;\r
+ }\r
+\r
+#if TRACE\r
+ Log("\n");\r
+ Log("DA Simple DPL:\n");\r
+ ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple DPM:\n");\r
+ ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple DPD:\n");\r
+ ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple DPE:\n");\r
+ ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple DPI:\n");\r
+ ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple DPJ:\n");\r
+ ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple TBM:\n");\r
+ ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple TBD:\n");\r
+ ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple TBE:\n");\r
+ ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple TBI:\n");\r
+ ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+ Log("\n");\r
+ Log("DA Simple TBJ:\n");\r
+ ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
+#endif\r
+\r
+// Trace-back\r
+// ==========\r
+ Path.Clear();\r
+\r
+// Find last edge\r
+ SCORE M = DPM(uLengthA, uLengthB);\r
+ SCORE D = DPD(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose;\r
+ SCORE E = DPE(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose2;\r
+ SCORE I = DPI(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose;\r
+ SCORE J = DPJ(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose2;\r
+ char cEdgeType = '?';\r
+\r
+ SCORE BestScore = M;\r
+ cEdgeType = 'M';\r
+ if (D > BestScore)\r
+ {\r
+ cEdgeType = 'D';\r
+ BestScore = D;\r
+ }\r
+ if (E > BestScore)\r
+ {\r
+ cEdgeType = 'E';\r
+ BestScore = E;\r
+ }\r
+ if (I > BestScore)\r
+ {\r
+ cEdgeType = 'I';\r
+ BestScore = I;\r
+ }\r
+ if (J > BestScore)\r
+ {\r
+ cEdgeType = 'J';\r
+ BestScore = J;\r
+ }\r
+\r
+#if TRACE\r
+ Log("DA Simple: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",\r
+ M, D, E, I, J, cEdgeType);\r
+#endif\r
+\r
+ unsigned PLA = uLengthA;\r
+ unsigned PLB = uLengthB;\r
+ for (;;)\r
+ {\r
+ PWEdge Edge;\r
+ Edge.cType = XlatEdgeType(cEdgeType);\r
+ Edge.uPrefixLengthA = PLA;\r
+ Edge.uPrefixLengthB = PLB;\r
+#if TRACE\r
+ Log("Prepend %c%d.%d\n", Edge.cType, PLA, PLB);\r
+#endif\r
+ Path.PrependEdge(Edge);\r
+\r
+ switch (cEdgeType)\r
+ {\r
+ case 'M':\r
+ assert(PLA > 0);\r
+ assert(PLB > 0);\r
+ cEdgeType = TBM(PLA, PLB);\r
+ --PLA;\r
+ --PLB;\r
+ break;\r
+\r
+ case 'D':\r
+ assert(PLA > 0);\r
+ cEdgeType = TBD(PLA, PLB);\r
+ --PLA;\r
+ break;\r
+\r
+ case 'E':\r
+ assert(PLA > 0);\r
+ cEdgeType = TBE(PLA, PLB);\r
+ --PLA;\r
+ break;\r
+\r
+ case 'I':\r
+ assert(PLB > 0);\r
+ cEdgeType = TBI(PLA, PLB);\r
+ --PLB;\r
+ break;\r
+ \r
+ case 'J':\r
+ assert(PLB > 0);\r
+ cEdgeType = TBJ(PLA, PLB);\r
+ --PLB;\r
+ break;\r
+\r
+ default:\r
+ Quit("Invalid edge %c", cEdgeType);\r
+ }\r
+ if (0 == PLA && 0 == PLB)\r
+ break;\r
+ }\r
+ Path.Validate();\r
+\r
+// SCORE Score = TraceBack(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, Path);\r
+\r
+#if TRACE\r
+ SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);\r
+ Path.LogMe();\r
+ Log("Score = %s Path = %s\n", LocalScoreToStr(BestScore), LocalScoreToStr(scorePath));\r
+#endif\r
+\r
+ if (g_bKeepSimpleDP)\r
+ {\r
+ g_DPM = DPM_;\r
+ g_DPD = DPD_;\r
+ g_DPE = DPE_;\r
+ g_DPI = DPI_;\r
+ g_DPJ = DPJ_;\r
+\r
+ g_TBM = TBM_;\r
+ g_TBD = TBD_;\r
+ g_TBE = TBE_;\r
+ g_TBI = TBI_;\r
+ g_TBJ = TBJ_;\r
+ }\r
+ else\r
+ {\r
+ delete[] DPM_;\r
+ delete[] DPD_;\r
+ delete[] DPE_;\r
+ delete[] DPI_;\r
+ delete[] DPJ_;\r
+\r
+ delete[] TBM_;\r
+ delete[] TBD_;\r
+ delete[] TBE_;\r
+ delete[] TBI_;\r
+ delete[] TBJ_;\r
+ }\r
+\r
+ return BestScore;\r
+ }\r
+\r
+#endif // DOUBLE_AFFINE\r
git://source.jalview.org / jabaws.git / blobdiff