--- /dev/null
+#include "muscle.h"\r
+#include <math.h>\r
+#include "pwpath.h"\r
+#include "profile.h"\r
+#include <stdio.h>\r
+\r
+// Textbook Smith-Waterman affine gap implementation.\r
+\r
+#define TRACE 0\r
+\r
+static const char *LocalScoreToStr(SCORE s)\r
+ {\r
+ static char str[16];\r
+ if (MINUS_INFINITY == s)\r
+ return " *";\r
+ sprintf(str, "%6.2f", s);\r
+ return str;\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 SW(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 *DPM_ = new SCORE[LM];\r
+ SCORE *DPD_ = new SCORE[LM];\r
+ SCORE *DPI_ = new SCORE[LM];\r
+\r
+ DPM(0, 0) = 0;\r
+ DPD(0, 0) = MINUS_INFINITY;\r
+ DPI(0, 0) = MINUS_INFINITY;\r
+\r
+ DPM(1, 0) = MINUS_INFINITY;\r
+ DPD(1, 0) = MINUS_INFINITY;\r
+ DPI(1, 0) = MINUS_INFINITY;\r
+\r
+ DPM(0, 1) = MINUS_INFINITY;\r
+ DPD(0, 1) = MINUS_INFINITY;\r
+ DPI(0, 1) = MINUS_INFINITY;\r
+\r
+// Empty prefix of B is special case\r
+ for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
+ {\r
+ // M=LetterA+LetterB, impossible with empty prefix\r
+ DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+\r
+ // D=LetterA+GapB, never optimal in local alignment with gap penalties\r
+ DPD(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+\r
+ // I=GapA+LetterB, impossible with empty prefix\r
+ DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
+ }\r
+\r
+// Empty prefix of A is special case\r
+ for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
+ {\r
+ // M=LetterA+LetterB, impossible with empty prefix\r
+ DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
+\r
+ // D=LetterA+GapB, impossible with empty prefix\r
+ DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
+\r
+ // I=GapA+LetterB, never optimal in local alignment with gap penalties\r
+ DPI(0, uPrefixLengthB) = MINUS_INFINITY;\r
+ }\r
+\r
+ SCORE scoreMax = MINUS_INFINITY;\r
+ unsigned uPrefixLengthAMax = uInsane;\r
+ unsigned uPrefixLengthBMax = uInsane;\r
+\r
+// ============\r
+// Main DP loop\r
+// ============\r
+ SCORE scoreGapCloseB = 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
+ 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
+\r
+ SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
+ SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
+ SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMM >= scoreDM && scoreMM >= scoreIM)\r
+ scoreBest = scoreMM;\r
+ else if (scoreDM >= scoreMM && scoreDM >= scoreIM)\r
+ scoreBest = scoreDM;\r
+ else \r
+ {\r
+ assert(scoreIM >= scoreMM && scoreIM >= scoreDM);\r
+ scoreBest = scoreIM;\r
+ }\r
+ if (scoreBest < 0)\r
+ scoreBest = 0;\r
+ scoreBest += scoreLL;\r
+ if (scoreBest > scoreMax)\r
+ {\r
+ scoreMax = scoreBest;\r
+ uPrefixLengthAMax = uPrefixLengthA;\r
+ uPrefixLengthBMax = uPrefixLengthB;\r
+ }\r
+ DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest;\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);\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMD >= scoreDD)\r
+ scoreBest = scoreMD;\r
+ else\r
+ {\r
+ assert(scoreDD >= scoreMD);\r
+ scoreBest = scoreDD;\r
+ }\r
+ DPD(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);\r
+\r
+ SCORE scoreBest;\r
+ if (scoreMI >= scoreII)\r
+ scoreBest = scoreMI;\r
+ else \r
+ {\r
+ assert(scoreII > scoreMI);\r
+ scoreBest = scoreII;\r
+ }\r
+ DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
+ }\r
+\r
+ scoreGapCloseA = PPA.m_scoreGapClose;\r
+ }\r
+ scoreGapCloseB = PPB.m_scoreGapClose;\r
+ }\r
+\r
+#if TRACE\r
+ Log("DPM:\n");\r
+ ListDP(DPM_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
+ Log("DPD:\n");\r
+ ListDP(DPD_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
+ Log("DPI:\n");\r
+ ListDP(DPI_, PA, PB, uPrefixLengthA, uPrefixLengthB);\r
+#endif\r
+\r
+ assert(scoreMax == DPM(uPrefixLengthAMax, uPrefixLengthBMax));\r
+ TraceBackSW(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, \r
+ uPrefixLengthAMax, uPrefixLengthBMax, 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(scoreMax), LocalScoreToStr(scorePath));\r
+#endif\r
+\r
+ delete[] DPM_;\r
+ delete[] DPD_;\r
+ delete[] DPI_;\r
+\r
+ return scoreMax;\r
+ }\r
git://source.jalview.org / jabaws.git / blobdiff