9 // NW double affine small memory, term gaps fully penalized
\r
10 // (so up to caller to adjust in profile if desired).
\r
14 #define MIN(x, y) ((x) < (y) ? (x) : (y))
\r
17 extern bool g_bKeepSimpleDP;
\r
18 extern SCORE *g_DPM;
\r
19 extern SCORE *g_DPD;
\r
20 extern SCORE *g_DPE;
\r
21 extern SCORE *g_DPI;
\r
22 extern SCORE *g_DPJ;
\r
31 #define ALLOC_TRACE() \
\r
32 const SCORE UNINIT = MINUS_INFINITY; \
\r
33 const size_t LM = uPrefixCountA*uPrefixCountB; \
\r
35 SCORE *DPM_ = new SCORE[LM]; \
\r
36 SCORE *DPD_ = new SCORE[LM]; \
\r
37 SCORE *DPE_ = new SCORE[LM]; \
\r
38 SCORE *DPI_ = new SCORE[LM]; \
\r
39 SCORE *DPJ_ = new SCORE[LM]; \
\r
41 char *TBM_ = new char[LM]; \
\r
42 char *TBD_ = new char[LM]; \
\r
43 char *TBE_ = new char[LM]; \
\r
44 char *TBI_ = new char[LM]; \
\r
45 char *TBJ_ = new char[LM]; \
\r
47 memset(TBM_, '?', LM); \
\r
48 memset(TBD_, '?', LM); \
\r
49 memset(TBE_, '?', LM); \
\r
50 memset(TBI_, '?', LM); \
\r
51 memset(TBJ_, '?', LM); \
\r
53 for (unsigned i = 0; i <= uLengthA; ++i) \
\r
54 for (unsigned j = 0; j <= uLengthB; ++j) \
\r
56 DPM(i, j) = UNINIT; \
\r
57 DPD(i, j) = UNINIT; \
\r
58 DPE(i, j) = UNINIT; \
\r
59 DPI(i, j) = UNINIT; \
\r
60 DPJ(i, j) = UNINIT; \
\r
63 #define ALLOC_TRACE()
\r
67 #define SetDPM(i, j, x) DPM(i, j) = x
\r
68 #define SetDPD(i, j, x) DPD(i, j) = x
\r
69 #define SetDPE(i, j, x) DPE(i, j) = x
\r
70 #define SetDPI(i, j, x) DPI(i, j) = x
\r
71 #define SetDPJ(i, j, x) DPJ(i, j) = x
\r
72 #define SetTBM(i, j, x) TBM(i, j) = x
\r
73 #define SetTBD(i, j, x) TBD(i, j) = x
\r
74 #define SetTBE(i, j, x) TBE(i, j) = x
\r
75 #define SetTBI(i, j, x) TBI(i, j) = x
\r
76 #define SetTBJ(i, j, x) TBJ(i, j) = x
\r
78 #define SetDPM(i, j, x) /* empty */
\r
79 #define SetDPD(i, j, x) /* empty */
\r
80 #define SetDPE(i, j, x) /* empty */
\r
81 #define SetDPI(i, j, x) /* empty */
\r
82 #define SetDPJ(i, j, x) /* empty */
\r
83 #define SetTBM(i, j, x) /* empty */
\r
84 #define SetTBD(i, j, x) /* empty */
\r
85 #define SetTBE(i, j, x) /* empty */
\r
86 #define SetTBI(i, j, x) /* empty */
\r
87 #define SetTBJ(i, j, x) /* empty */
\r
90 #define RECURSE_D(i, j) \
\r
92 SCORE DD = DRow[j] + e; \
\r
93 SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\
\r
97 SetTBD(i, j, 'D'); \
\r
102 SetBitTBD(TB, i, j, 'M'); \
\r
103 SetTBD(i, j, 'M'); \
\r
105 SetDPD(i, j, DRow[j]); \
\r
108 #define RECURSE_E(i, j) \
\r
110 SCORE EE = ERow[j] + e2; \
\r
111 SCORE ME = MPrev[j] + PA[i-1].m_scoreGapOpen2;\
\r
115 SetTBE(i, j, 'E'); \
\r
120 SetBitTBE(TB, i, j, 'M'); \
\r
121 SetTBE(i, j, 'M'); \
\r
123 SetDPE(i, j, ERow[j]); \
\r
126 #define RECURSE_D_ATerm(j) RECURSE_D(uLengthA, j)
\r
127 #define RECURSE_E_ATerm(j) RECURSE_E(uLengthA, j)
\r
129 #define RECURSE_D_BTerm(j) RECURSE_D(i, uLengthB)
\r
130 #define RECURSE_E_BTerm(j) RECURSE_E(i, uLengthB)
\r
132 #define RECURSE_I(i, j) \
\r
135 SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\
\r
139 SetBitTBI(TB, i, j, 'M'); \
\r
140 SetTBI(i, j, 'M'); \
\r
143 SetTBI(i, j, 'I'); \
\r
144 SetDPI(i, j, Iij); \
\r
147 #define RECURSE_J(i, j) \
\r
150 SCORE MJ = MCurr[j-1] + PB[j-1].m_scoreGapOpen2;\
\r
154 SetBitTBJ(TB, i, j, 'M'); \
\r
155 SetTBJ(i, j, 'M'); \
\r
158 SetTBJ(i, j, 'I'); \
\r
159 SetDPJ(i, j, Jij); \
\r
162 #define RECURSE_I_ATerm(j) RECURSE_I(uLengthA, j)
\r
163 #define RECURSE_J_ATerm(j) RECURSE_J(uLengthA, j)
\r
165 #define RECURSE_I_BTerm(j) RECURSE_I(i, uLengthB)
\r
166 #define RECURSE_J_BTerm(j) RECURSE_J(i, uLengthB)
\r
168 #define RECURSE_M(i, j) \
\r
170 SCORE Best = MCurr[j]; /* MM */ \
\r
171 SetTBM(i+1, j+1, 'M'); \
\r
172 SetBitTBM(TB, i+1, j+1, 'M'); \
\r
174 SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose; \
\r
178 SetTBM(i+1, j+1, 'D'); \
\r
179 SetBitTBM(TB, i+1, j+1, 'D'); \
\r
182 SCORE EM = ERow[j] + PA[i-1].m_scoreGapClose2; \
\r
186 SetTBM(i+1, j+1, 'E'); \
\r
187 SetBitTBM(TB, i+1, j+1, 'E'); \
\r
190 SCORE IM = Iij + PB[j-1].m_scoreGapClose; \
\r
194 SetTBM(i+1, j+1, 'I'); \
\r
195 SetBitTBM(TB, i+1, j+1, 'I'); \
\r
198 SCORE JM = Jij + PB[j-1].m_scoreGapClose2; \
\r
202 SetTBM(i+1, j+1, 'J'); \
\r
203 SetBitTBM(TB, i+1, j+1, 'J'); \
\r
205 MNext[j+1] += Best; \
\r
206 SetDPM(i+1, j+1, MNext[j+1]); \
\r
210 static bool LocalEq(BASETYPE b1, BASETYPE b2)
\r
212 if (b1 < -100000 && b2 < -100000)
\r
214 double diff = fabs(b1 - b2);
\r
217 double sum = fabs(b1) + fabs(b2);
\r
218 return diff/sum < 0.005;
\r
221 static char Get_M_Char(char Bits)
\r
223 switch (Bits & BIT_xM)
\r
240 static char Get_D_Char(char Bits)
\r
242 return (Bits & BIT_xD) ? 'M' : 'D';
\r
245 static char Get_E_Char(char Bits)
\r
247 return (Bits & BIT_xE) ? 'M' : 'E';
\r
250 static char Get_I_Char(char Bits)
\r
252 return (Bits & BIT_xI) ? 'M' : 'I';
\r
255 static char Get_J_Char(char Bits)
\r
257 return (Bits & BIT_xJ) ? 'M' : 'J';
\r
260 static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,
\r
261 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
265 Log("***DPDIFF*** DP%c=NULL\n", c);
\r
269 SCORE *DPM_ = g_DP;
\r
270 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
271 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
272 if (!LocalEq(DPM(i, j), DPD(i, j)))
\r
274 Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Small = %.2g\n",
\r
275 c, i, j, DPM(i, j), DPD(i, j));
\r
281 static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBE_, char *TBI_, char *TBJ_,
\r
282 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
284 if (!g_bKeepSimpleDP)
\r
286 SCORE *DPM_ = g_DPM;
\r
288 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
289 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
291 char c1 = TBM(i, j);
\r
292 char c2 = Get_M_Char(TB[i][j]);
\r
293 if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)
\r
295 Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
302 SCORE *DPD_ = g_DPD;
\r
303 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
304 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
306 char c1 = TBD(i, j);
\r
307 char c2 = Get_D_Char(TB[i][j]);
\r
308 if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)
\r
310 Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
316 SCORE *DPE_ = g_DPE;
\r
319 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
320 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
322 char c1 = TBE(i, j);
\r
323 char c2 = Get_E_Char(TB[i][j]);
\r
324 if (c1 != '?' && c1 != c2 && DPE(i, j) > -100000)
\r
326 Log("TBE(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
332 SCORE *DPI_ = g_DPI;
\r
333 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
334 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
336 char c1 = TBI(i, j);
\r
337 char c2 = Get_I_Char(TB[i][j]);
\r
338 if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)
\r
340 Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
346 SCORE *DPJ_ = g_DPJ;
\r
349 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
350 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
352 char c1 = TBJ(i, j);
\r
353 char c2 = Get_J_Char(TB[i][j]);
\r
354 if (c1 != '?' && c1 != c2 && DPJ(i, j) > -100000)
\r
356 Log("TBJ(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
363 Log("TB success\n");
\r
367 static const char *LocalScoreToStr(SCORE s)
\r
369 static char str[16];
\r
372 sprintf(str, "%6.1f", s);
\r
376 static void LogDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
377 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
380 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
383 if (uPrefixLengthB > 0)
\r
384 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
385 Log(" %4u:%c", uPrefixLengthB, c);
\r
388 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
391 if (uPrefixLengthA > 0)
\r
392 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
393 Log("%4u:%c ", uPrefixLengthA, c);
\r
394 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
395 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
400 static void LogBitTB(char **TB, const ProfPos *PA, const ProfPos *PB,
\r
401 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
404 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
407 if (uPrefixLengthB > 0)
\r
408 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
409 Log(" %4u:%c", uPrefixLengthB, c);
\r
413 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
416 if (uPrefixLengthA > 0)
\r
417 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
418 Log("%4u:%c ", uPrefixLengthA, c);
\r
419 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
421 char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
429 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
432 if (uPrefixLengthA > 0)
\r
433 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
434 Log("%4u:%c ", uPrefixLengthA, c);
\r
435 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
437 char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
445 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
448 if (uPrefixLengthA > 0)
\r
449 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
450 Log("%4u:%c ", uPrefixLengthA, c);
\r
451 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
453 char c = Get_E_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
461 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
464 if (uPrefixLengthA > 0)
\r
465 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
466 Log("%4u:%c ", uPrefixLengthA, c);
\r
467 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
469 char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
477 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
480 if (uPrefixLengthA > 0)
\r
481 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
482 Log("%4u:%c ", uPrefixLengthA, c);
\r
483 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
485 char c = Get_J_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
492 static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
493 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
496 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
499 if (uPrefixLengthB > 0)
\r
500 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
501 Log(" %4u:%c", uPrefixLengthB, c);
\r
504 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
507 if (uPrefixLengthA > 0)
\r
508 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
509 Log("%4u:%c ", uPrefixLengthA, c);
\r
510 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
512 char c = TBM(uPrefixLengthA, uPrefixLengthB);
\r
519 static const char *BitsToStr(char Bits)
\r
521 static char Str[32];
\r
523 sprintf(Str, "%cM %cD %cE %cI %cJ",
\r
532 static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
\r
557 TB[i][j] &= ~BIT_xM;
\r
561 static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
\r
575 TB[i][j] &= ~BIT_xD;
\r
579 static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
\r
593 TB[i][j] &= ~BIT_xI;
\r
598 static inline void SetBitTBE(char **TB, unsigned i, unsigned j, char c)
\r
612 TB[i][j] &= ~BIT_xE;
\r
616 static inline void SetBitTBJ(char **TB, unsigned i, unsigned j, char c)
\r
630 TB[i][j] &= ~BIT_xJ;
\r
636 #define LogMatrices() \
\r
638 Log("Bit DPM:\n"); \
\r
639 LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
640 Log("Bit DPD:\n"); \
\r
641 LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
642 Log("Bit DPE:\n"); \
\r
643 LogDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
644 Log("Bit DPI:\n"); \
\r
645 LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
646 Log("Bit DPJ:\n"); \
\r
647 LogDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
648 Log("Bit TB:\n"); \
\r
649 LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
651 Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
\r
653 Log("DPM success\n"); \
\r
654 Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
\r
656 Log("DPD success\n"); \
\r
657 Same = DPEq('E', g_DPE, DPE_, uPrefixCountA, uPrefixCountB);\
\r
659 Log("DPE success\n"); \
\r
660 Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
\r
662 Log("DPI success\n"); \
\r
663 Same = DPEq('J', g_DPJ, DPJ_, uPrefixCountA, uPrefixCountB);\
\r
665 Log("DPJ success\n"); \
\r
666 CompareTB(TB, g_TBM, g_TBD, g_TBE, g_TBI, g_TBJ, uPrefixCountA, uPrefixCountB);\
\r
669 #define LogMatrices() /* empty */
\r
672 SCORE NWDASmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
673 unsigned uLengthB, PWPath &Path)
\r
675 assert(uLengthB > 0 && uLengthA > 0);
\r
677 ProfPos *pa0 = (ProfPos *) PA;
\r
678 ProfPos *pb0 = (ProfPos *) PB;
\r
679 ProfPos *paa = (ProfPos *) (PA + uLengthA - 1);
\r
680 ProfPos *pbb = (ProfPos *) (PB + uLengthB - 1);
\r
682 pa0->m_scoreGapOpen *= -1;
\r
683 pb0->m_scoreGapOpen *= -1;
\r
685 paa->m_scoreGapClose *= -1;
\r
686 pbb->m_scoreGapClose *= -1;
\r
688 pa0->m_scoreGapOpen2 *= -1;
\r
689 pb0->m_scoreGapOpen2 *= -1;
\r
690 paa->m_scoreGapClose2 *= -1;
\r
691 pbb->m_scoreGapClose2 *= -1;
\r
693 const unsigned uPrefixCountA = uLengthA + 1;
\r
694 const unsigned uPrefixCountB = uLengthB + 1;
\r
695 const SCORE e = g_scoreGapExtend;
\r
697 const SCORE e2 = g_scoreGapExtend2;
\r
698 const SCORE min_e = MIN(g_scoreGapExtend, g_scoreGapExtend2);
\r
702 SCORE *MCurr = new SCORE[uPrefixCountB];
\r
703 SCORE *MNext = new SCORE[uPrefixCountB];
\r
704 SCORE *MPrev = new SCORE[uPrefixCountB];
\r
705 SCORE *DRow = new SCORE[uPrefixCountB];
\r
706 SCORE *ERow = new SCORE[uPrefixCountB];
\r
708 char **TB = new char *[uPrefixCountA];
\r
709 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
711 TB[i] = new char [uPrefixCountB];
\r
712 memset(TB[i], 0, uPrefixCountB);
\r
715 SCORE Iij = MINUS_INFINITY;
\r
718 SCORE Jij = MINUS_INFINITY;
\r
721 Iij = PB[0].m_scoreGapOpen;
\r
724 Jij = PB[0].m_scoreGapOpen2;
\r
727 for (unsigned j = 2; j <= uLengthB; ++j)
\r
739 for (unsigned j = 0; j <= uLengthB; ++j)
\r
741 DRow[j] = MINUS_INFINITY;
\r
742 ERow[j] = MINUS_INFINITY;
\r
744 SetDPD(0, j, DRow[j]);
\r
745 SetDPE(0, j, ERow[j]);
\r
752 SetDPM(0, 0, MPrev[0]);
\r
753 for (unsigned j = 1; j <= uLengthB; ++j)
\r
755 MPrev[j] = MINUS_INFINITY;
\r
756 SetDPM(0, j, MPrev[j]);
\r
759 MCurr[0] = MINUS_INFINITY;
\r
760 SetDPM(1, 0, MCurr[0]);
\r
762 MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
\r
763 SetDPM(1, 1, MCurr[1]);
\r
764 SetBitTBM(TB, 1, 1, 'M');
\r
767 for (unsigned j = 2; j <= uLengthB; ++j)
\r
769 SCORE M = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
\r
770 (j - 2)*e + PB[j-2].m_scoreGapClose;
\r
771 SCORE M2 = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen2 +
\r
772 (j - 2)*e2 + PB[j-2].m_scoreGapClose2;
\r
777 SetBitTBM(TB, 1, j, 'I');
\r
783 SetBitTBM(TB, 1, j, 'J');
\r
786 SetDPM(1, j, MCurr[j]);
\r
790 for (unsigned i = 1; i < uLengthA; ++i)
\r
792 Iij = MINUS_INFINITY;
\r
793 Jij = MINUS_INFINITY;
\r
797 DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
\r
798 ERow[0] = PA[0].m_scoreGapOpen2 + (i - 1)*e2;
\r
799 SetDPD(i, 0, DRow[0]);
\r
800 SetDPE(i, 0, ERow[0]);
\r
802 MCurr[0] = MINUS_INFINITY;
\r
805 MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
\r
806 SetBitTBM(TB, i, 1, 'M');
\r
811 SCORE M = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
\r
812 (i - 2)*e + PA[i-2].m_scoreGapClose;
\r
813 SCORE M2 = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen2 +
\r
814 (i - 2)*e2 + PA[i-2].m_scoreGapClose2;
\r
818 SetBitTBM(TB, i, 1, 'D');
\r
824 SetBitTBM(TB, i, 1, 'E');
\r
828 SetDPM(i, 0, MCurr[0]);
\r
829 SetDPM(i, 1, MCurr[1]);
\r
831 for (unsigned j = 1; j < uLengthB; ++j)
\r
832 MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);
\r
834 for (unsigned j = 1; j < uLengthB; ++j)
\r
842 // Special case for j=uLengthB
\r
848 // Prev := Curr, Curr := Next, Next := Prev
\r
849 Rotate(MPrev, MCurr, MNext);
\r
852 // Special case for i=uLengthA
\r
853 MCurr[0] = MINUS_INFINITY;
\r
854 SCORE M = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
\r
855 PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
\r
856 SCORE M2 = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
\r
857 PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
\r
861 SetBitTBM(TB, uLengthA, 1, 'D');
\r
862 SetTBM(uLengthA, 1, 'D');
\r
867 SetBitTBM(TB, uLengthA, 1, 'E');
\r
868 SetTBM(uLengthA, 1, 'D');
\r
870 SetDPM(uLengthA, 0, MCurr[0]);
\r
871 SetDPM(uLengthA, 1, MCurr[1]);
\r
873 DRow[0] = MINUS_INFINITY;
\r
874 ERow[0] = MINUS_INFINITY;
\r
876 SetDPD(uLengthA, 0, DRow[0]);
\r
877 SetDPE(uLengthA, 0, ERow[0]);
\r
879 for (unsigned j = 1; j <= uLengthB; ++j)
\r
881 RECURSE_D_ATerm(j);
\r
882 RECURSE_E_ATerm(j);
\r
885 Iij = MINUS_INFINITY;
\r
886 Jij = MINUS_INFINITY;
\r
888 for (unsigned j = 1; j <= uLengthB; ++j)
\r
896 SCORE MAB = MCurr[uLengthB];
\r
897 SCORE DAB = DRow[uLengthB] + PA[uLengthA-1].m_scoreGapClose;
\r
898 SCORE EAB = ERow[uLengthB] + PA[uLengthA-1].m_scoreGapClose2;
\r
899 SCORE IAB = Iij + PB[uLengthB-1].m_scoreGapClose;
\r
900 SCORE JAB = Jij + PB[uLengthB-1].m_scoreGapClose2;
\r
903 char cEdgeType = 'M';
\r
926 Log(" Small: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",
\r
927 MAB, DAB, EAB, IAB, JAB, cEdgeType);
\r
930 BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);
\r
941 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
947 #endif // DOUBLE_AFFINE
\r