+++ /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