12 extern bool g_bKeepSimpleDP;
\r
13 extern SCORE *g_DPM;
\r
14 extern SCORE *g_DPD;
\r
15 extern SCORE *g_DPI;
\r
22 #define ALLOC_TRACE() \
\r
23 const SCORE UNINIT = MINUS_INFINITY; \
\r
24 const size_t LM = uPrefixCountA*uPrefixCountB; \
\r
26 SCORE *DPM_ = new SCORE[LM]; \
\r
27 SCORE *DPD_ = new SCORE[LM]; \
\r
28 SCORE *DPI_ = new SCORE[LM]; \
\r
30 char *TBM_ = new char[LM]; \
\r
31 char *TBD_ = new char[LM]; \
\r
32 char *TBI_ = new char[LM]; \
\r
34 memset(TBM_, '?', LM); \
\r
35 memset(TBD_, '?', LM); \
\r
36 memset(TBI_, '?', LM); \
\r
38 for (unsigned i = 0; i <= uLengthA; ++i) \
\r
39 for (unsigned j = 0; j <= uLengthB; ++j) \
\r
41 DPM(i, j) = UNINIT; \
\r
42 DPD(i, j) = UNINIT; \
\r
43 DPI(i, j) = UNINIT; \
\r
46 #define ALLOC_TRACE()
\r
50 #define SetDPM(i, j, x) DPM(i, j) = x
\r
51 #define SetDPD(i, j, x) DPD(i, j) = x
\r
52 #define SetDPI(i, j, x) DPI(i, j) = x
\r
53 #define SetTBM(i, j, x) TBM(i, j) = x
\r
54 #define SetTBD(i, j, x) TBD(i, j) = x
\r
55 #define SetTBI(i, j, x) TBI(i, j) = x
\r
57 #define SetDPM(i, j, x) /* empty */
\r
58 #define SetDPD(i, j, x) /* empty */
\r
59 #define SetDPI(i, j, x) /* empty */
\r
60 #define SetTBM(i, j, x) /* empty */
\r
61 #define SetTBD(i, j, x) /* empty */
\r
62 #define SetTBI(i, j, x) /* empty */
\r
65 #define RECURSE_D(i, j) \
\r
67 SCORE DD = DRow[j] + e; \
\r
68 SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\
\r
72 SetTBD(i, j, 'D'); \
\r
77 /* SetBitTBD(TB, i, j, 'M'); */ \
\r
78 TBRow[j] &= ~BIT_xD; \
\r
79 TBRow[j] |= BIT_MD; \
\r
80 SetTBD(i, j, 'M'); \
\r
82 SetDPD(i, j, DRow[j]); \
\r
85 #define RECURSE_D_ATerm(j) RECURSE_D(uLengthA, j)
\r
86 #define RECURSE_D_BTerm(j) RECURSE_D(i, uLengthB)
\r
88 #define RECURSE_I(i, j) \
\r
91 SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\
\r
95 /* SetBitTBI(TB, i, j, 'M'); */ \
\r
96 TBRow[j] &= ~BIT_xI; \
\r
97 TBRow[j] |= BIT_MI; \
\r
98 SetTBI(i, j, 'M'); \
\r
101 SetTBI(i, j, 'I'); \
\r
102 SetDPI(i, j, Iij); \
\r
105 #define RECURSE_I_ATerm(j) RECURSE_I(uLengthA, j)
\r
106 #define RECURSE_I_BTerm(j) RECURSE_I(i, uLengthB)
\r
108 #define RECURSE_M(i, j) \
\r
110 SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose; \
\r
111 SCORE IM = Iij + PB[j-1].m_scoreGapClose; \
\r
112 SCORE MM = MCurr[j]; \
\r
113 TB[i+1][j+1] &= ~BIT_xM; \
\r
114 if (MM >= DM && MM >= IM) \
\r
116 MNext[j+1] += MM; \
\r
117 SetDPM(i+1, j+1, MNext[j+1]); \
\r
118 SetTBM(i+1, j+1, 'M'); \
\r
119 /* SetBitTBM(TB, i+1, j+1, 'M'); */ \
\r
120 TB[i+1][j+1] |= BIT_MM; \
\r
122 else if (DM >= MM && DM >= IM) \
\r
124 MNext[j+1] += DM; \
\r
125 SetDPM(i+1, j+1, MNext[j+1]); \
\r
126 SetTBM(i+1, j+1, 'D'); \
\r
127 /* SetBitTBM(TB, i+1, j+1, 'D'); */ \
\r
128 TB[i+1][j+1] |= BIT_DM; \
\r
132 assert(IM >= MM && IM >= DM); \
\r
133 MNext[j+1] += IM; \
\r
134 SetDPM(i+1, j+1, MNext[j+1]); \
\r
135 SetTBM(i+1, j+1, 'I'); \
\r
136 /* SetBitTBM(TB, i+1, j+1, 'I'); */ \
\r
137 TB[i+1][j+1] |= BIT_IM; \
\r
142 static bool LocalEq(BASETYPE b1, BASETYPE b2)
\r
144 if (b1 < -100000 && b2 < -100000)
\r
146 double diff = fabs(b1 - b2);
\r
149 double sum = fabs(b1) + fabs(b2);
\r
150 return diff/sum < 0.005;
\r
153 static char Get_M_Char(char Bits)
\r
155 switch (Bits & BIT_xM)
\r
168 static char Get_D_Char(char Bits)
\r
170 return (Bits & BIT_xD) ? 'M' : 'D';
\r
173 static char Get_I_Char(char Bits)
\r
175 return (Bits & BIT_xI) ? 'M' : 'I';
\r
178 static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,
\r
179 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
181 SCORE *DPM_ = g_DP;
\r
182 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
183 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
184 if (!LocalEq(DPM(i, j), DPD(i, j)))
\r
186 Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Fast = %.2g\n",
\r
187 c, i, j, DPM(i, j), DPD(i, j));
\r
193 static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBI_,
\r
194 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
196 SCORE *DPM_ = g_DPM;
\r
198 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
199 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
201 char c1 = TBM(i, j);
\r
202 char c2 = Get_M_Char(TB[i][j]);
\r
203 if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)
\r
205 Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
212 SCORE *DPD_ = g_DPD;
\r
213 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
214 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
216 char c1 = TBD(i, j);
\r
217 char c2 = Get_D_Char(TB[i][j]);
\r
218 if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)
\r
220 Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
226 SCORE *DPI_ = g_DPI;
\r
227 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
228 for (unsigned j = 0; j < uPrefixCountB; ++j)
\r
230 char c1 = TBI(i, j);
\r
231 char c2 = Get_I_Char(TB[i][j]);
\r
232 if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)
\r
234 Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
\r
241 Log("TB success\n");
\r
245 static const char *LocalScoreToStr(SCORE s)
\r
247 static char str[16];
\r
250 sprintf(str, "%6.1f", s);
\r
254 static void LogDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
255 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
258 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
261 if (uPrefixLengthB > 0)
\r
262 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
263 Log(" %4u:%c", uPrefixLengthB, c);
\r
266 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
269 if (uPrefixLengthA > 0)
\r
270 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
271 Log("%4u:%c ", uPrefixLengthA, c);
\r
272 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
273 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
278 static void LogBitTB(char **TB, const ProfPos *PA, const ProfPos *PB,
\r
279 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
282 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
285 if (uPrefixLengthB > 0)
\r
286 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
287 Log(" %4u:%c", uPrefixLengthB, c);
\r
291 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
294 if (uPrefixLengthA > 0)
\r
295 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
296 Log("%4u:%c ", uPrefixLengthA, c);
\r
297 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
299 char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
307 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
310 if (uPrefixLengthA > 0)
\r
311 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
312 Log("%4u:%c ", uPrefixLengthA, c);
\r
313 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
315 char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
323 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
326 if (uPrefixLengthA > 0)
\r
327 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
328 Log("%4u:%c ", uPrefixLengthA, c);
\r
329 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
331 char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);
\r
338 static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
339 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
342 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
345 if (uPrefixLengthB > 0)
\r
346 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
347 Log(" %4u:%c", uPrefixLengthB, c);
\r
350 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
353 if (uPrefixLengthA > 0)
\r
354 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
355 Log("%4u:%c ", uPrefixLengthA, c);
\r
356 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
358 char c = TBM(uPrefixLengthA, uPrefixLengthB);
\r
365 static const char *BitsToStr(char Bits)
\r
367 static char Str[9];
\r
369 sprintf(Str, "%cM %cD %cI",
\r
376 static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
\r
393 TB[i][j] &= ~BIT_xM;
\r
397 static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
\r
411 TB[i][j] &= ~BIT_xD;
\r
415 static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
\r
429 TB[i][j] &= ~BIT_xI;
\r
434 #define LogMatrices() \
\r
436 Log("Bit DPM:\n"); \
\r
437 LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
438 Log("Bit DPD:\n"); \
\r
439 LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
440 Log("Bit DPI:\n"); \
\r
441 LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
442 Log("Bit TB:\n"); \
\r
443 LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB); \
\r
445 Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
\r
447 Log("DPM success\n"); \
\r
448 Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
\r
450 Log("DPD success\n"); \
\r
451 Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
\r
453 Log("DPI success\n"); \
\r
454 CompareTB(TB, g_TBM, g_TBD, g_TBI, uPrefixCountA, uPrefixCountB);\
\r
457 #define LogMatrices() /* empty */
\r
460 static unsigned uCachePrefixCountB;
\r
461 static unsigned uCachePrefixCountA;
\r
462 static SCORE *CacheMCurr;
\r
463 static SCORE *CacheMNext;
\r
464 static SCORE *CacheMPrev;
\r
465 static SCORE *CacheDRow;
\r
466 static char **CacheTB;
\r
468 static void AllocCache(unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
470 if (uPrefixCountA <= uCachePrefixCountA && uPrefixCountB <= uCachePrefixCountB)
\r
473 delete[] CacheMCurr;
\r
474 delete[] CacheMNext;
\r
475 delete[] CacheMPrev;
\r
476 delete[] CacheDRow;
\r
477 for (unsigned i = 0; i < uCachePrefixCountA; ++i)
\r
478 delete[] CacheTB[i];
\r
481 uCachePrefixCountA = uPrefixCountA + 1024;
\r
482 uCachePrefixCountB = uPrefixCountB + 1024;
\r
484 CacheMCurr = new SCORE[uCachePrefixCountB];
\r
485 CacheMNext = new SCORE[uCachePrefixCountB];
\r
486 CacheMPrev = new SCORE[uCachePrefixCountB];
\r
487 CacheDRow = new SCORE[uCachePrefixCountB];
\r
489 CacheTB = new char *[uCachePrefixCountA];
\r
490 for (unsigned i = 0; i < uCachePrefixCountA; ++i)
\r
491 CacheTB[i] = new char [uCachePrefixCountB];
\r
494 SCORE NWSmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
495 unsigned uLengthB, PWPath &Path)
\r
497 if (0 == uLengthB || 0 == uLengthA )
\r
498 Quit("Internal error, NWSmall: length=0");
\r
500 SetTermGaps(PA, uLengthA);
\r
501 SetTermGaps(PB, uLengthB);
\r
503 const unsigned uPrefixCountA = uLengthA + 1;
\r
504 const unsigned uPrefixCountB = uLengthB + 1;
\r
505 const SCORE e = g_scoreGapExtend;
\r
509 AllocCache(uPrefixCountA, uPrefixCountB);
\r
511 SCORE *MCurr = CacheMCurr;
\r
512 SCORE *MNext = CacheMNext;
\r
513 SCORE *MPrev = CacheMPrev;
\r
514 SCORE *DRow = CacheDRow;
\r
516 char **TB = CacheTB;
\r
517 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
518 memset(TB[i], 0, uPrefixCountB);
\r
520 SCORE Iij = MINUS_INFINITY;
\r
523 Iij = PB[0].m_scoreGapOpen;
\r
526 for (unsigned j = 2; j <= uLengthB; ++j)
\r
533 for (unsigned j = 0; j <= uLengthB; ++j)
\r
535 DRow[j] = MINUS_INFINITY;
\r
536 SetDPD(0, j, DRow[j]);
\r
541 SetDPM(0, 0, MPrev[0]);
\r
542 for (unsigned j = 1; j <= uLengthB; ++j)
\r
544 MPrev[j] = MINUS_INFINITY;
\r
545 SetDPM(0, j, MPrev[j]);
\r
548 MCurr[0] = MINUS_INFINITY;
\r
549 SetDPM(1, 0, MCurr[0]);
\r
551 MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
\r
552 SetDPM(1, 1, MCurr[1]);
\r
553 SetBitTBM(TB, 1, 1, 'M');
\r
556 for (unsigned j = 2; j <= uLengthB; ++j)
\r
558 MCurr[j] = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
\r
559 (j - 2)*e + PB[j-2].m_scoreGapClose;
\r
560 SetDPM(1, j, MCurr[j]);
\r
561 SetBitTBM(TB, 1, j, 'I');
\r
566 for (unsigned i = 1; i < uLengthA; ++i)
\r
568 char *TBRow = TB[i];
\r
570 Iij = MINUS_INFINITY;
\r
573 DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
\r
574 SetDPD(i, 0, DRow[0]);
\r
576 MCurr[0] = MINUS_INFINITY;
\r
579 MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
\r
580 SetBitTBM(TB, i, 1, 'M');
\r
585 MCurr[1] = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
\r
586 (i - 2)*e + PA[i-2].m_scoreGapClose;
\r
587 SetBitTBM(TB, i, 1, 'D');
\r
590 SetDPM(i, 0, MCurr[0]);
\r
591 SetDPM(i, 1, MCurr[1]);
\r
593 for (unsigned j = 1; j < uLengthB; ++j)
\r
594 MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);
\r
596 for (unsigned j = 1; j < uLengthB; ++j)
\r
602 // Special case for j=uLengthB
\r
606 // Prev := Curr, Curr := Next, Next := Prev
\r
607 Rotate(MPrev, MCurr, MNext);
\r
610 // Special case for i=uLengthA
\r
611 char *TBRow = TB[uLengthA];
\r
612 MCurr[0] = MINUS_INFINITY;
\r
614 MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
\r
615 PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
\r
617 MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + PA[0].m_scoreGapOpen +
\r
618 PA[0].m_scoreGapClose;
\r
619 SetBitTBM(TB, uLengthA, 1, 'D');
\r
620 SetTBM(uLengthA, 1, 'D');
\r
621 SetDPM(uLengthA, 0, MCurr[0]);
\r
622 SetDPM(uLengthA, 1, MCurr[1]);
\r
624 DRow[0] = MINUS_INFINITY;
\r
625 SetDPD(uLengthA, 0, DRow[0]);
\r
626 for (unsigned j = 1; j <= uLengthB; ++j)
\r
627 RECURSE_D_ATerm(j);
\r
629 Iij = MINUS_INFINITY;
\r
630 for (unsigned j = 1; j <= uLengthB; ++j)
\r
635 SCORE MAB = MCurr[uLengthB];
\r
636 SCORE DAB = DRow[uLengthB];
\r
640 char cEdgeType = 'M';
\r
653 Log(" Fast: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n",
\r
654 MAB, DAB, IAB, cEdgeType);
\r
657 BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);
\r