+++ /dev/null
-#include "muscle.h"\r
-#include <math.h>\r
-#include "pwpath.h"\r
-#include "profile.h"\r
-#include <stdio.h>\r
-\r
-// NW small memory\r
-\r
-#define TRACE 0\r
-\r
-#if TRACE\r
-extern bool g_bKeepSimpleDP;\r
-extern SCORE *g_DPM;\r
-extern SCORE *g_DPD;\r
-extern SCORE *g_DPI;\r
-extern char *g_TBM;\r
-extern char *g_TBD;\r
-extern char *g_TBI;\r
-#endif\r
-\r
-#if TRACE\r
-#define ALLOC_TRACE() \\r
- const SCORE UNINIT = MINUS_INFINITY; \\r
- const size_t LM = uPrefixCountA*uPrefixCountB; \\r
- \\r
- SCORE *DPM_ = new SCORE[LM]; \\r
- SCORE *DPD_ = new SCORE[LM]; \\r
- SCORE *DPI_ = new SCORE[LM]; \\r
- \\r
- char *TBM_ = new char[LM]; \\r
- char *TBD_ = new char[LM]; \\r
- char *TBI_ = new char[LM]; \\r
- \\r
- memset(TBM_, '?', LM); \\r
- memset(TBD_, '?', LM); \\r
- memset(TBI_, '?', LM); \\r
- \\r
- for (unsigned i = 0; i <= uLengthA; ++i) \\r
- for (unsigned j = 0; j <= uLengthB; ++j) \\r
- { \\r
- DPM(i, j) = UNINIT; \\r
- DPD(i, j) = UNINIT; \\r
- DPI(i, j) = UNINIT; \\r
- }\r
-#else\r
-#define ALLOC_TRACE()\r
-#endif\r
-\r
-#if TRACE\r
-#define SetDPM(i, j, x) DPM(i, j) = x\r
-#define SetDPD(i, j, x) DPD(i, j) = x\r
-#define SetDPI(i, j, x) DPI(i, j) = x\r
-#define SetTBM(i, j, x) TBM(i, j) = x\r
-#define SetTBD(i, j, x) TBD(i, j) = x\r
-#define SetTBI(i, j, x) TBI(i, j) = x\r
-#else\r
-#define SetDPM(i, j, x) /* empty */\r
-#define SetDPD(i, j, x) /* empty */\r
-#define SetDPI(i, j, x) /* empty */\r
-#define SetTBM(i, j, x) /* empty */\r
-#define SetTBD(i, j, x) /* empty */\r
-#define SetTBI(i, j, x) /* empty */\r
-#endif\r
-\r
-#define RECURSE_D(i, j) \\r
- { \\r
- SCORE DD = DRow[j] + e; \\r
- SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\\r
- if (DD > MD) \\r
- { \\r
- DRow[j] = DD; \\r
- SetTBD(i, j, 'D'); \\r
- } \\r
- else \\r
- { \\r
- DRow[j] = MD; \\r
- /* SetBitTBD(TB, i, j, 'M'); */ \\r
- TBRow[j] &= ~BIT_xD; \\r
- TBRow[j] |= BIT_MD; \\r
- SetTBD(i, j, 'M'); \\r
- } \\r
- SetDPD(i, j, DRow[j]); \\r
- }\r
-\r
-#define RECURSE_D_ATerm(j) RECURSE_D(uLengthA, j)\r
-#define RECURSE_D_BTerm(j) RECURSE_D(i, uLengthB)\r
-\r
-#define RECURSE_I(i, j) \\r
- { \\r
- Iij += e; \\r
- SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\\r
- if (MI >= Iij) \\r
- { \\r
- Iij = MI; \\r
- /* SetBitTBI(TB, i, j, 'M'); */ \\r
- TBRow[j] &= ~BIT_xI; \\r
- TBRow[j] |= BIT_MI; \\r
- SetTBI(i, j, 'M'); \\r
- } \\r
- else \\r
- SetTBI(i, j, 'I'); \\r
- SetDPI(i, j, Iij); \\r
- }\r
-\r
-#define RECURSE_I_ATerm(j) RECURSE_I(uLengthA, j)\r
-#define RECURSE_I_BTerm(j) RECURSE_I(i, uLengthB)\r
-\r
-#define RECURSE_M(i, j) \\r
- { \\r
- SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose; \\r
- SCORE IM = Iij + PB[j-1].m_scoreGapClose; \\r
- SCORE MM = MCurr[j]; \\r
- TB[i+1][j+1] &= ~BIT_xM; \\r
- if (MM >= DM && MM >= IM) \\r
- { \\r
- MNext[j+1] += MM; \\r
- SetDPM(i+1, j+1, MNext[j+1]); \\r
- SetTBM(i+1, j+1, 'M'); \\r
- /* SetBitTBM(TB, i+1, j+1, 'M'); */ \\r
- TB[i+1][j+1] |= BIT_MM; \\r
- } \\r
- else if (DM >= MM && DM >= IM) \\r
- { \\r
- MNext[j+1] += DM; \\r
- SetDPM(i+1, j+1, MNext[j+1]); \\r
- SetTBM(i+1, j+1, 'D'); \\r
- /* SetBitTBM(TB, i+1, j+1, 'D'); */ \\r
- TB[i+1][j+1] |= BIT_DM; \\r
- } \\r
- else \\r
- { \\r
- assert(IM >= MM && IM >= DM); \\r
- MNext[j+1] += IM; \\r
- SetDPM(i+1, j+1, MNext[j+1]); \\r
- SetTBM(i+1, j+1, 'I'); \\r
- /* SetBitTBM(TB, i+1, j+1, 'I'); */ \\r
- TB[i+1][j+1] |= BIT_IM; \\r
- } \\r
- }\r
-\r
-#if TRACE\r
-static bool LocalEq(BASETYPE b1, BASETYPE b2)\r
- {\r
- if (b1 < -100000 && b2 < -100000)\r
- return true;\r
- double diff = fabs(b1 - b2);\r
- if (diff < 0.0001)\r
- return true;\r
- double sum = fabs(b1) + fabs(b2);\r
- return diff/sum < 0.005;\r
- }\r
-\r
-static char Get_M_Char(char Bits)\r
- {\r
- switch (Bits & BIT_xM)\r
- {\r
- case BIT_MM:\r
- return 'M';\r
- case BIT_DM:\r
- return 'D';\r
- case BIT_IM:\r
- return 'I';\r
- }\r
- Quit("Huh?");\r
- return '?';\r
- }\r
-\r
-static char Get_D_Char(char Bits)\r
- {\r
- return (Bits & BIT_xD) ? 'M' : 'D';\r
- }\r
-\r
-static char Get_I_Char(char Bits)\r
- {\r
- return (Bits & BIT_xI) ? 'M' : 'I';\r
- }\r
-\r
-static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- SCORE *DPM_ = g_DP;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- if (!LocalEq(DPM(i, j), DPD(i, j)))\r
- {\r
- Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Fast = %.2g\n",\r
- c, i, j, DPM(i, j), DPD(i, j));\r
- return false;\r
- }\r
- return true;\r
- }\r
-\r
-static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBI_, \r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- SCORE *DPM_ = g_DPM;\r
- bool Eq = true;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- {\r
- char c1 = TBM(i, j);\r
- char c2 = Get_M_Char(TB[i][j]);\r
- if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)\r
- {\r
- Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);\r
- Eq = false;\r
- goto D;\r
- }\r
- }\r
-\r
-D:\r
- SCORE *DPD_ = g_DPD;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- {\r
- char c1 = TBD(i, j);\r
- char c2 = Get_D_Char(TB[i][j]);\r
- if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)\r
- {\r
- Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);\r
- Eq = false;\r
- goto I;\r
- }\r
- }\r
-I:\r
- SCORE *DPI_ = g_DPI;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- {\r
- char c1 = TBI(i, j);\r
- char c2 = Get_I_Char(TB[i][j]);\r
- if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)\r
- {\r
- Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);\r
- Eq = false;\r
- goto Done;\r
- }\r
- }\r
-Done:\r
- if (Eq)\r
- Log("TB success\n");\r
- return Eq;\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 LogDP(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
-static void LogBitTB(char **TB, 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
- Log("Bit TBM:\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
- {\r
- char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);\r
- Log(" %6c", c);\r
- }\r
- Log("\n");\r
- }\r
-\r
- Log("\n");\r
- Log("Bit TBD:\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
- {\r
- char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);\r
- Log(" %6c", c);\r
- }\r
- Log("\n");\r
- }\r
-\r
- Log("\n");\r
- Log("Bit TBI:\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
- {\r
- char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);\r
- Log(" %6c", c);\r
- }\r
- Log("\n");\r
- }\r
- }\r
-\r
-static void ListTB(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
- {\r
- char c = TBM(uPrefixLengthA, uPrefixLengthB);\r
- Log(" %6c", c);\r
- }\r
- Log("\n");\r
- }\r
- }\r
-\r
-static const char *BitsToStr(char Bits)\r
- {\r
- static char Str[9];\r
-\r
- sprintf(Str, "%cM %cD %cI",\r
- Get_M_Char(Bits),\r
- Get_D_Char(Bits),\r
- Get_I_Char(Bits));\r
- }\r
-#endif // TRACE\r
-\r
-static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)\r
- {\r
- char Bit;\r
- switch (c)\r
- {\r
- case 'M':\r
- Bit = BIT_MM;\r
- break;\r
- case 'D':\r
- Bit = BIT_DM;\r
- break;\r
- case 'I':\r
- Bit = BIT_IM;\r
- break;\r
- default:\r
- Quit("Huh?!");\r
- }\r
- TB[i][j] &= ~BIT_xM;\r
- TB[i][j] |= Bit;\r
- }\r
-\r
-static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)\r
- {\r
- char Bit;\r
- switch (c)\r
- {\r
- case 'M':\r
- Bit = BIT_MD;\r
- break;\r
- case 'D':\r
- Bit = BIT_DD;\r
- break;\r
- default:\r
- Quit("Huh?!");\r
- }\r
- TB[i][j] &= ~BIT_xD;\r
- TB[i][j] |= Bit;\r
- }\r
-\r
-static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)\r
- {\r
- char Bit;\r
- switch (c)\r
- {\r
- case 'M':\r
- Bit = BIT_MI;\r
- break;\r
- case 'I':\r
- Bit = BIT_II;\r
- break;\r
- default:\r
- Quit("Huh?!");\r
- }\r
- TB[i][j] &= ~BIT_xI;\r
- TB[i][j] |= Bit;\r
- }\r
-\r
-#if TRACE\r
-#define LogMatrices() \\r
- { \\r
- Log("Bit DPM:\n"); \\r
- LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- Log("Bit DPD:\n"); \\r
- LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- Log("Bit DPI:\n"); \\r
- LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- Log("Bit TB:\n"); \\r
- LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- bool Same; \\r
- Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPM success\n"); \\r
- Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPD success\n"); \\r
- Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPI success\n"); \\r
- CompareTB(TB, g_TBM, g_TBD, g_TBI, uPrefixCountA, uPrefixCountB);\\r
- }\r
-#else\r
-#define LogMatrices() /* empty */\r
-#endif\r
-\r
-static unsigned uCachePrefixCountB;\r
-static unsigned uCachePrefixCountA;\r
-static SCORE *CacheMCurr;\r
-static SCORE *CacheMNext;\r
-static SCORE *CacheMPrev;\r
-static SCORE *CacheDRow;\r
-static char **CacheTB;\r
-\r
-static void AllocCache(unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- if (uPrefixCountA <= uCachePrefixCountA && uPrefixCountB <= uCachePrefixCountB)\r
- return;\r
-\r
- delete[] CacheMCurr;\r
- delete[] CacheMNext;\r
- delete[] CacheMPrev;\r
- delete[] CacheDRow;\r
- for (unsigned i = 0; i < uCachePrefixCountA; ++i)\r
- delete[] CacheTB[i];\r
- delete[] CacheTB;\r
-\r
- uCachePrefixCountA = uPrefixCountA + 1024;\r
- uCachePrefixCountB = uPrefixCountB + 1024;\r
-\r
- CacheMCurr = new SCORE[uCachePrefixCountB];\r
- CacheMNext = new SCORE[uCachePrefixCountB];\r
- CacheMPrev = new SCORE[uCachePrefixCountB];\r
- CacheDRow = new SCORE[uCachePrefixCountB];\r
-\r
- CacheTB = new char *[uCachePrefixCountA];\r
- for (unsigned i = 0; i < uCachePrefixCountA; ++i)\r
- CacheTB[i] = new char [uCachePrefixCountB];\r
- }\r
-\r
-SCORE NWSmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
- unsigned uLengthB, PWPath &Path)\r
- {\r
- if (0 == uLengthB || 0 == uLengthA )\r
- Quit("Internal error, NWSmall: length=0");\r
-\r
- SetTermGaps(PA, uLengthA);\r
- SetTermGaps(PB, uLengthB);\r
-\r
- const unsigned uPrefixCountA = uLengthA + 1;\r
- const unsigned uPrefixCountB = uLengthB + 1;\r
- const SCORE e = g_scoreGapExtend;\r
-\r
- ALLOC_TRACE()\r
-\r
- AllocCache(uPrefixCountA, uPrefixCountB);\r
-\r
- SCORE *MCurr = CacheMCurr;\r
- SCORE *MNext = CacheMNext;\r
- SCORE *MPrev = CacheMPrev;\r
- SCORE *DRow = CacheDRow;\r
-\r
- char **TB = CacheTB;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- memset(TB[i], 0, uPrefixCountB);\r
-\r
- SCORE Iij = MINUS_INFINITY;\r
- SetDPI(0, 0, Iij);\r
-\r
- Iij = PB[0].m_scoreGapOpen;\r
- SetDPI(0, 1, Iij);\r
-\r
- for (unsigned j = 2; j <= uLengthB; ++j)\r
- {\r
- Iij += e;\r
- SetDPI(0, j, Iij);\r
- SetTBI(0, j, 'I');\r
- }\r
-\r
- for (unsigned j = 0; j <= uLengthB; ++j)\r
- {\r
- DRow[j] = MINUS_INFINITY;\r
- SetDPD(0, j, DRow[j]);\r
- SetTBD(0, j, 'D');\r
- }\r
-\r
- MPrev[0] = 0;\r
- SetDPM(0, 0, MPrev[0]);\r
- for (unsigned j = 1; j <= uLengthB; ++j)\r
- {\r
- MPrev[j] = MINUS_INFINITY;\r
- SetDPM(0, j, MPrev[j]);\r
- }\r
-\r
- MCurr[0] = MINUS_INFINITY;\r
- SetDPM(1, 0, MCurr[0]);\r
-\r
- MCurr[1] = ScoreProfPos2(PA[0], PB[0]);\r
- SetDPM(1, 1, MCurr[1]);\r
- SetBitTBM(TB, 1, 1, 'M');\r
- SetTBM(1, 1, 'M');\r
-\r
- for (unsigned j = 2; j <= uLengthB; ++j)\r
- {\r
- MCurr[j] = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +\r
- (j - 2)*e + PB[j-2].m_scoreGapClose;\r
- SetDPM(1, j, MCurr[j]);\r
- SetBitTBM(TB, 1, j, 'I');\r
- SetTBM(1, j, 'I');\r
- }\r
-\r
-// Main DP loop\r
- for (unsigned i = 1; i < uLengthA; ++i)\r
- {\r
- char *TBRow = TB[i];\r
-\r
- Iij = MINUS_INFINITY;\r
- SetDPI(i, 0, Iij);\r
-\r
- DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;\r
- SetDPD(i, 0, DRow[0]);\r
-\r
- MCurr[0] = MINUS_INFINITY; \r
- if (i == 1)\r
- {\r
- MCurr[1] = ScoreProfPos2(PA[0], PB[0]);\r
- SetBitTBM(TB, i, 1, 'M');\r
- SetTBM(i, 1, 'M');\r
- }\r
- else\r
- {\r
- MCurr[1] = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +\r
- (i - 2)*e + PA[i-2].m_scoreGapClose;\r
- SetBitTBM(TB, i, 1, 'D');\r
- SetTBM(i, 1, 'D');\r
- }\r
- SetDPM(i, 0, MCurr[0]);\r
- SetDPM(i, 1, MCurr[1]);\r
-\r
- for (unsigned j = 1; j < uLengthB; ++j)\r
- MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);\r
-\r
- for (unsigned j = 1; j < uLengthB; ++j)\r
- {\r
- RECURSE_D(i, j)\r
- RECURSE_I(i, j)\r
- RECURSE_M(i, j)\r
- }\r
- // Special case for j=uLengthB\r
- RECURSE_D_BTerm(i)\r
- RECURSE_I_BTerm(i)\r
-\r
- // Prev := Curr, Curr := Next, Next := Prev\r
- Rotate(MPrev, MCurr, MNext);\r
- }\r
-\r
-// Special case for i=uLengthA\r
- char *TBRow = TB[uLengthA];\r
- MCurr[0] = MINUS_INFINITY;\r
- if (uLengthA > 1)\r
- MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +\r
- PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;\r
- else\r
- MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + PA[0].m_scoreGapOpen +\r
- PA[0].m_scoreGapClose;\r
- SetBitTBM(TB, uLengthA, 1, 'D');\r
- SetTBM(uLengthA, 1, 'D');\r
- SetDPM(uLengthA, 0, MCurr[0]);\r
- SetDPM(uLengthA, 1, MCurr[1]);\r
-\r
- DRow[0] = MINUS_INFINITY;\r
- SetDPD(uLengthA, 0, DRow[0]);\r
- for (unsigned j = 1; j <= uLengthB; ++j)\r
- RECURSE_D_ATerm(j);\r
-\r
- Iij = MINUS_INFINITY;\r
- for (unsigned j = 1; j <= uLengthB; ++j)\r
- RECURSE_I_ATerm(j)\r
-\r
- LogMatrices();\r
-\r
- SCORE MAB = MCurr[uLengthB];\r
- SCORE DAB = DRow[uLengthB];\r
- SCORE IAB = Iij;\r
-\r
- SCORE Score = MAB;\r
- char cEdgeType = 'M';\r
- if (DAB > Score)\r
- {\r
- Score = DAB;\r
- cEdgeType = 'D';\r
- }\r
- if (IAB > Score)\r
- {\r
- Score = IAB;\r
- cEdgeType = 'I';\r
- }\r
-\r
-#if TRACE\r
- Log(" Fast: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n",\r
- MAB, DAB, IAB, cEdgeType);\r
-#endif\r
-\r
- BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);\r
-\r
-#if DBEUG\r
- Path.Validate();\r
-#endif\r
-\r
- return 0;\r
- }\r