+++ /dev/null
-#include "muscle.h"\r
-#include <math.h>\r
-#include "pwpath.h"\r
-#include "profile.h"\r
-#include <stdio.h>\r
-\r
-#if DOUBLE_AFFINE\r
-\r
-// NW double affine small memory, term gaps fully penalized\r
-// (so up to caller to adjust in profile if desired).\r
-\r
-#define TRACE 0\r
-\r
-#define MIN(x, y) ((x) < (y) ? (x) : (y))\r
-\r
-#if TRACE\r
-extern bool g_bKeepSimpleDP;\r
-extern SCORE *g_DPM;\r
-extern SCORE *g_DPD;\r
-extern SCORE *g_DPE;\r
-extern SCORE *g_DPI;\r
-extern SCORE *g_DPJ;\r
-extern char *g_TBM;\r
-extern char *g_TBD;\r
-extern char *g_TBE;\r
-extern char *g_TBI;\r
-extern char *g_TBJ;\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 *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
- 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
- DPE(i, j) = UNINIT; \\r
- DPI(i, j) = UNINIT; \\r
- DPJ(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 SetDPE(i, j, x) DPE(i, j) = x\r
-#define SetDPI(i, j, x) DPI(i, j) = x\r
-#define SetDPJ(i, j, x) DPJ(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 SetTBE(i, j, x) TBE(i, j) = x\r
-#define SetTBI(i, j, x) TBI(i, j) = x\r
-#define SetTBJ(i, j, x) TBJ(i, j) = x\r
-#else\r
-#define SetDPM(i, j, x) /* empty */\r
-#define SetDPD(i, j, x) /* empty */\r
-#define SetDPE(i, j, x) /* empty */\r
-#define SetDPI(i, j, x) /* empty */\r
-#define SetDPJ(i, j, x) /* empty */\r
-#define SetTBM(i, j, x) /* empty */\r
-#define SetTBD(i, j, x) /* empty */\r
-#define SetTBE(i, j, x) /* empty */\r
-#define SetTBI(i, j, x) /* empty */\r
-#define SetTBJ(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
- SetTBD(i, j, 'M'); \\r
- } \\r
- SetDPD(i, j, DRow[j]); \\r
- }\r
-\r
-#define RECURSE_E(i, j) \\r
- { \\r
- SCORE EE = ERow[j] + e2; \\r
- SCORE ME = MPrev[j] + PA[i-1].m_scoreGapOpen2;\\r
- if (EE > ME) \\r
- { \\r
- ERow[j] = EE; \\r
- SetTBE(i, j, 'E'); \\r
- } \\r
- else \\r
- { \\r
- ERow[j] = ME; \\r
- SetBitTBE(TB, i, j, 'M'); \\r
- SetTBE(i, j, 'M'); \\r
- } \\r
- SetDPE(i, j, ERow[j]); \\r
- }\r
-\r
-#define RECURSE_D_ATerm(j) RECURSE_D(uLengthA, j)\r
-#define RECURSE_E_ATerm(j) RECURSE_E(uLengthA, j)\r
-\r
-#define RECURSE_D_BTerm(j) RECURSE_D(i, uLengthB)\r
-#define RECURSE_E_BTerm(j) RECURSE_E(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
- SetTBI(i, j, 'M'); \\r
- } \\r
- else \\r
- SetTBI(i, j, 'I'); \\r
- SetDPI(i, j, Iij); \\r
- }\r
-\r
-#define RECURSE_J(i, j) \\r
- { \\r
- Jij += e2; \\r
- SCORE MJ = MCurr[j-1] + PB[j-1].m_scoreGapOpen2;\\r
- if (MJ >= Jij) \\r
- { \\r
- Jij = MJ; \\r
- SetBitTBJ(TB, i, j, 'M'); \\r
- SetTBJ(i, j, 'M'); \\r
- } \\r
- else \\r
- SetTBJ(i, j, 'I'); \\r
- SetDPJ(i, j, Jij); \\r
- }\r
-\r
-#define RECURSE_I_ATerm(j) RECURSE_I(uLengthA, j)\r
-#define RECURSE_J_ATerm(j) RECURSE_J(uLengthA, j)\r
-\r
-#define RECURSE_I_BTerm(j) RECURSE_I(i, uLengthB)\r
-#define RECURSE_J_BTerm(j) RECURSE_J(i, uLengthB)\r
-\r
-#define RECURSE_M(i, j) \\r
- { \\r
- SCORE Best = MCurr[j]; /* MM */ \\r
- SetTBM(i+1, j+1, 'M'); \\r
- SetBitTBM(TB, i+1, j+1, 'M'); \\r
- \\r
- SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose; \\r
- if (DM > Best) \\r
- { \\r
- Best = DM; \\r
- SetTBM(i+1, j+1, 'D'); \\r
- SetBitTBM(TB, i+1, j+1, 'D'); \\r
- } \\r
- \\r
- SCORE EM = ERow[j] + PA[i-1].m_scoreGapClose2; \\r
- if (EM > Best) \\r
- { \\r
- Best = EM; \\r
- SetTBM(i+1, j+1, 'E'); \\r
- SetBitTBM(TB, i+1, j+1, 'E'); \\r
- } \\r
- \\r
- SCORE IM = Iij + PB[j-1].m_scoreGapClose; \\r
- if (IM > Best) \\r
- { \\r
- Best = IM; \\r
- SetTBM(i+1, j+1, 'I'); \\r
- SetBitTBM(TB, i+1, j+1, 'I'); \\r
- } \\r
- \\r
- SCORE JM = Jij + PB[j-1].m_scoreGapClose2; \\r
- if (JM > Best) \\r
- { \\r
- Best = JM; \\r
- SetTBM(i+1, j+1, 'J'); \\r
- SetBitTBM(TB, i+1, j+1, 'J'); \\r
- } \\r
- MNext[j+1] += Best; \\r
- SetDPM(i+1, j+1, MNext[j+1]); \\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_EM:\r
- return 'E';\r
- case BIT_IM:\r
- return 'I';\r
- case BIT_JM:\r
- return 'J';\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_E_Char(char Bits)\r
- {\r
- return (Bits & BIT_xE) ? 'M' : 'E';\r
- }\r
-\r
-static char Get_I_Char(char Bits)\r
- {\r
- return (Bits & BIT_xI) ? 'M' : 'I';\r
- }\r
-\r
-static char Get_J_Char(char Bits)\r
- {\r
- return (Bits & BIT_xJ) ? 'M' : 'J';\r
- }\r
-\r
-static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- if (0 == g_DP)\r
- {\r
- Log("***DPDIFF*** DP%c=NULL\n", c);\r
- return true;\r
- }\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, Small = %.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 *TBE_, char *TBI_, char *TBJ_,\r
- unsigned uPrefixCountA, unsigned uPrefixCountB)\r
- {\r
- if (!g_bKeepSimpleDP)\r
- return true;\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 E;\r
- }\r
- }\r
-E:\r
- SCORE *DPE_ = g_DPE;\r
- if (0 == TBE_)\r
- goto I;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- {\r
- char c1 = TBE(i, j);\r
- char c2 = Get_E_Char(TB[i][j]);\r
- if (c1 != '?' && c1 != c2 && DPE(i, j) > -100000)\r
- {\r
- Log("TBE(%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 J;\r
- }\r
- }\r
-J:\r
- SCORE *DPJ_ = g_DPJ;\r
- if (0 == DPJ_)\r
- goto Done;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- for (unsigned j = 0; j < uPrefixCountB; ++j)\r
- {\r
- char c1 = TBJ(i, j);\r
- char c2 = Get_J_Char(TB[i][j]);\r
- if (c1 != '?' && c1 != c2 && DPJ(i, j) > -100000)\r
- {\r
- Log("TBJ(%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 TBE:\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_E_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
- Log("\n");\r
- Log("Bit TBJ:\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_J_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[32];\r
-\r
- sprintf(Str, "%cM %cD %cE %cI %cJ",\r
- Get_M_Char(Bits),\r
- Get_D_Char(Bits),\r
- Get_E_Char(Bits),\r
- Get_I_Char(Bits),\r
- Get_J_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
-#if DOUBLE_AFFINE\r
- case 'E':\r
- Bit = BIT_EM;\r
- break;\r
- case 'I':\r
- Bit = BIT_IM;\r
- break;\r
- case 'J':\r
- Bit = BIT_JM;\r
- break;\r
-#endif\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 DOUBLE_AFFINE\r
-static inline void SetBitTBE(char **TB, unsigned i, unsigned j, char c)\r
- {\r
- char Bit;\r
- switch (c)\r
- {\r
- case 'M':\r
- Bit = BIT_ME;\r
- break;\r
- case 'E':\r
- Bit = BIT_EE;\r
- break;\r
- default:\r
- Quit("Huh?!");\r
- }\r
- TB[i][j] &= ~BIT_xE;\r
- TB[i][j] |= Bit;\r
- }\r
-\r
-static inline void SetBitTBJ(char **TB, unsigned i, unsigned j, char c)\r
- {\r
- char Bit;\r
- switch (c)\r
- {\r
- case 'M':\r
- Bit = BIT_MJ;\r
- break;\r
- case 'J':\r
- Bit = BIT_JJ;\r
- break;\r
- default:\r
- Quit("Huh?!");\r
- }\r
- TB[i][j] &= ~BIT_xJ;\r
- TB[i][j] |= Bit;\r
- }\r
-#endif\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 DPE:\n"); \\r
- LogDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- Log("Bit DPI:\n"); \\r
- LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB); \\r
- Log("Bit DPJ:\n"); \\r
- LogDP(DPJ_, 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('E', g_DPE, DPE_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPE success\n"); \\r
- Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPI success\n"); \\r
- Same = DPEq('J', g_DPJ, DPJ_, uPrefixCountA, uPrefixCountB);\\r
- if (Same) \\r
- Log("DPJ success\n"); \\r
- CompareTB(TB, g_TBM, g_TBD, g_TBE, g_TBI, g_TBJ, uPrefixCountA, uPrefixCountB);\\r
- }\r
-#else\r
-#define LogMatrices() /* empty */\r
-#endif\r
-\r
-SCORE NWDASmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
- unsigned uLengthB, PWPath &Path)\r
- {\r
- assert(uLengthB > 0 && uLengthA > 0);\r
-\r
- ProfPos *pa0 = (ProfPos *) PA;\r
- ProfPos *pb0 = (ProfPos *) PB;\r
- ProfPos *paa = (ProfPos *) (PA + uLengthA - 1);\r
- ProfPos *pbb = (ProfPos *) (PB + uLengthB - 1);\r
-\r
- pa0->m_scoreGapOpen *= -1;\r
- pb0->m_scoreGapOpen *= -1;\r
-\r
- paa->m_scoreGapClose *= -1;\r
- pbb->m_scoreGapClose *= -1;\r
-\r
- pa0->m_scoreGapOpen2 *= -1;\r
- pb0->m_scoreGapOpen2 *= -1;\r
- paa->m_scoreGapClose2 *= -1;\r
- pbb->m_scoreGapClose2 *= -1;\r
-\r
- const unsigned uPrefixCountA = uLengthA + 1;\r
- const unsigned uPrefixCountB = uLengthB + 1;\r
- const SCORE e = g_scoreGapExtend;\r
-\r
- const SCORE e2 = g_scoreGapExtend2;\r
- const SCORE min_e = MIN(g_scoreGapExtend, g_scoreGapExtend2);\r
-\r
- ALLOC_TRACE()\r
-\r
- SCORE *MCurr = new SCORE[uPrefixCountB];\r
- SCORE *MNext = new SCORE[uPrefixCountB];\r
- SCORE *MPrev = new SCORE[uPrefixCountB];\r
- SCORE *DRow = new SCORE[uPrefixCountB];\r
- SCORE *ERow = new SCORE[uPrefixCountB];\r
-\r
- char **TB = new char *[uPrefixCountA];\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- {\r
- TB[i] = new char [uPrefixCountB];\r
- memset(TB[i], 0, uPrefixCountB);\r
- }\r
-\r
- SCORE Iij = MINUS_INFINITY;\r
- SetDPI(0, 0, Iij);\r
-\r
- SCORE Jij = MINUS_INFINITY;\r
- SetDPJ(0, 0, Jij);\r
-\r
- Iij = PB[0].m_scoreGapOpen;\r
- SetDPI(0, 1, Iij);\r
-\r
- Jij = PB[0].m_scoreGapOpen2;\r
- SetDPJ(0, 1, Jij);\r
-\r
- for (unsigned j = 2; j <= uLengthB; ++j)\r
- {\r
- Iij += e;\r
- Jij += e2;\r
-\r
- SetDPI(0, j, Iij);\r
- SetDPJ(0, j, Jij);\r
-\r
- SetTBI(0, j, 'I');\r
- SetTBJ(0, j, 'J');\r
- }\r
-\r
- for (unsigned j = 0; j <= uLengthB; ++j)\r
- {\r
- DRow[j] = MINUS_INFINITY;\r
- ERow[j] = MINUS_INFINITY;\r
-\r
- SetDPD(0, j, DRow[j]);\r
- SetDPE(0, j, ERow[j]);\r
-\r
- SetTBD(0, j, 'D');\r
- SetTBE(0, j, 'E');\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
- SCORE M = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +\r
- (j - 2)*e + PB[j-2].m_scoreGapClose;\r
- SCORE M2 = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen2 +\r
- (j - 2)*e2 + PB[j-2].m_scoreGapClose2;\r
- \r
- if (M >= M2)\r
- {\r
- MCurr[j] = M;\r
- SetBitTBM(TB, 1, j, 'I');\r
- SetTBM(1, j, 'I');\r
- }\r
- else\r
- {\r
- MCurr[j] = M2;\r
- SetBitTBM(TB, 1, j, 'J');\r
- SetTBM(1, j, 'J');\r
- }\r
- SetDPM(1, j, MCurr[j]);\r
- }\r
-\r
-// Main DP loop\r
- for (unsigned i = 1; i < uLengthA; ++i)\r
- {\r
- Iij = MINUS_INFINITY;\r
- Jij = MINUS_INFINITY;\r
- SetDPI(i, 0, Iij);\r
- SetDPJ(i, 0, Jij);\r
-\r
- DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;\r
- ERow[0] = PA[0].m_scoreGapOpen2 + (i - 1)*e2;\r
- SetDPD(i, 0, DRow[0]);\r
- SetDPE(i, 0, ERow[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
- SCORE M = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +\r
- (i - 2)*e + PA[i-2].m_scoreGapClose;\r
- SCORE M2 = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen2 +\r
- (i - 2)*e2 + PA[i-2].m_scoreGapClose2;\r
- if (M >= M2)\r
- {\r
- MCurr[1] = M;\r
- SetBitTBM(TB, i, 1, 'D');\r
- SetTBM(i, 1, 'D');\r
- }\r
- else\r
- {\r
- MCurr[1] = M2;\r
- SetBitTBM(TB, i, 1, 'E');\r
- SetTBM(i, 1, 'E');\r
- }\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_E(i, j)\r
- RECURSE_I(i, j)\r
- RECURSE_J(i, j)\r
- RECURSE_M(i, j)\r
- }\r
- // Special case for j=uLengthB\r
- RECURSE_D_BTerm(i)\r
- RECURSE_E_BTerm(i)\r
- RECURSE_I_BTerm(i)\r
- RECURSE_J_BTerm(i)\r
-\r
- // Prev := Curr, Curr := Next, Next := Prev\r
- Rotate(MPrev, MCurr, MNext);\r
- }\r
-\r
-// Special case for i=uLengthA\r
- MCurr[0] = MINUS_INFINITY;\r
- SCORE M = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +\r
- PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;\r
- SCORE M2 = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +\r
- PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;\r
- if (M >= M2)\r
- {\r
- MCurr[1] = M;\r
- SetBitTBM(TB, uLengthA, 1, 'D');\r
- SetTBM(uLengthA, 1, 'D');\r
- }\r
- else\r
- {\r
- MCurr[1] = M2;\r
- SetBitTBM(TB, uLengthA, 1, 'E');\r
- SetTBM(uLengthA, 1, 'D');\r
- }\r
- SetDPM(uLengthA, 0, MCurr[0]);\r
- SetDPM(uLengthA, 1, MCurr[1]);\r
-\r
- DRow[0] = MINUS_INFINITY;\r
- ERow[0] = MINUS_INFINITY;\r
-\r
- SetDPD(uLengthA, 0, DRow[0]);\r
- SetDPE(uLengthA, 0, ERow[0]);\r
-\r
- for (unsigned j = 1; j <= uLengthB; ++j)\r
- {\r
- RECURSE_D_ATerm(j);\r
- RECURSE_E_ATerm(j);\r
- }\r
-\r
- Iij = MINUS_INFINITY;\r
- Jij = MINUS_INFINITY;\r
-\r
- for (unsigned j = 1; j <= uLengthB; ++j)\r
- {\r
- RECURSE_I_ATerm(j)\r
- RECURSE_J_ATerm(j)\r
- }\r
-\r
- LogMatrices();\r
-\r
- SCORE MAB = MCurr[uLengthB];\r
- SCORE DAB = DRow[uLengthB] + PA[uLengthA-1].m_scoreGapClose;\r
- SCORE EAB = ERow[uLengthB] + PA[uLengthA-1].m_scoreGapClose2;\r
- SCORE IAB = Iij + PB[uLengthB-1].m_scoreGapClose;\r
- SCORE JAB = Jij + PB[uLengthB-1].m_scoreGapClose2;\r
-\r
- SCORE Score = MAB;\r
- char cEdgeType = 'M';\r
- if (DAB > Score)\r
- {\r
- Score = DAB;\r
- cEdgeType = 'D';\r
- }\r
- if (EAB > Score)\r
- {\r
- Score = EAB;\r
- cEdgeType = 'E';\r
- }\r
- if (IAB > Score)\r
- {\r
- Score = IAB;\r
- cEdgeType = 'I';\r
- }\r
- if (JAB > Score)\r
- {\r
- Score = JAB;\r
- cEdgeType = 'J';\r
- }\r
-\r
-#if TRACE\r
- Log(" Small: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",\r
- MAB, DAB, EAB, IAB, JAB, cEdgeType);\r
-#endif\r
-\r
- BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);\r
-\r
-#if DBEUG\r
- Path.Validate();\r
-#endif\r
-\r
- delete[] MCurr;\r
- delete[] MNext;\r
- delete[] MPrev;\r
- delete[] DRow;\r
- delete[] ERow;\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- delete[] TB[i];\r
- delete[] TB;\r
-\r
- return 0;\r
- }\r
-#endif // DOUBLE_AFFINE\r
git://source.jalview.org / jabaws.git / blobdiff