9 extern bool g_bKeepSimpleDP;
\r
10 extern SCORE *g_DPM;
\r
11 extern SCORE *g_DPD;
\r
12 extern SCORE *g_DPE;
\r
13 extern SCORE *g_DPI;
\r
14 extern SCORE *g_DPJ;
\r
21 static char XlatEdgeType(char c)
\r
30 static const char *LocalScoreToStr(SCORE s)
\r
32 static char str[16];
\r
35 sprintf(str, "%6.1f", s);
\r
39 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
40 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
43 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
46 if (uPrefixLengthB > 0)
\r
47 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
48 Log(" %4u:%c", uPrefixLengthB, c);
\r
51 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
54 if (uPrefixLengthA > 0)
\r
55 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
56 Log("%4u:%c ", uPrefixLengthA, c);
\r
57 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
58 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
63 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
64 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
67 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
70 if (uPrefixLengthB > 0)
\r
71 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
72 Log(" %4u:%c", uPrefixLengthB, c);
\r
75 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
78 if (uPrefixLengthA > 0)
\r
79 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
80 Log("%4u:%c ", uPrefixLengthA, c);
\r
81 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
82 Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));
\r
87 static void ListDPM(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
88 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
91 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
94 if (uPrefixLengthB > 0)
\r
95 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
96 Log(" %4u:%c", uPrefixLengthB, c);
\r
99 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
102 if (uPrefixLengthA > 0)
\r
103 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
104 Log("%4u:%c ", uPrefixLengthA, c);
\r
105 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
107 SCORE x = (uPrefixLengthA + uPrefixLengthB)*g_scoreGapExtend;
\r
108 SCORE s = DPM(uPrefixLengthA, uPrefixLengthB) - x;
\r
109 Log(" %s", LocalScoreToStr(s));
\r
115 extern SCORE ScoreProfPos2(const ProfPos &PP, const ProfPos &PPB);
\r
117 SCORE NWDASimple2(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
118 unsigned uLengthB, PWPath &Path)
\r
120 assert(uLengthB > 0 && uLengthA > 0);
\r
122 const unsigned uPrefixCountA = uLengthA + 1;
\r
123 const unsigned uPrefixCountB = uLengthB + 1;
\r
125 // Allocate DP matrices
\r
126 const size_t LM = uPrefixCountA*uPrefixCountB;
\r
127 SCORE *DPM_ = new SCORE[LM];
\r
128 SCORE *DPD_ = new SCORE[LM];
\r
129 SCORE *DPE_ = new SCORE[LM];
\r
130 SCORE *DPI_ = new SCORE[LM];
\r
131 SCORE *DPJ_ = new SCORE[LM];
\r
132 SCORE *DPL_ = new SCORE[LM];
\r
134 char *TBM_ = new char[LM];
\r
135 char *TBD_ = new char[LM];
\r
136 char *TBE_ = new char[LM];
\r
137 char *TBI_ = new char[LM];
\r
138 char *TBJ_ = new char[LM];
\r
140 memset(DPM_, 0, LM*sizeof(SCORE));
\r
141 memset(DPD_, 0, LM*sizeof(SCORE));
\r
142 memset(DPE_, 0, LM*sizeof(SCORE));
\r
143 memset(DPI_, 0, LM*sizeof(SCORE));
\r
144 memset(DPJ_, 0, LM*sizeof(SCORE));
\r
146 // memset(DPL_, 0, LM*sizeof(SCORE));
\r
148 memset(TBM_, '?', LM);
\r
149 memset(TBD_, '?', LM);
\r
150 memset(TBE_, '?', LM);
\r
151 memset(TBI_, '?', LM);
\r
152 memset(TBJ_, '?', LM);
\r
155 DPD(0, 0) = MINUS_INFINITY;
\r
156 DPE(0, 0) = MINUS_INFINITY;
\r
157 DPI(0, 0) = MINUS_INFINITY;
\r
158 DPJ(0, 0) = MINUS_INFINITY;
\r
160 DPM(1, 0) = MINUS_INFINITY;
\r
161 DPD(1, 0) = PA[0].m_scoreGapOpen;
\r
162 DPE(1, 0) = PA[0].m_scoreGapOpen2;
\r
163 DPI(1, 0) = MINUS_INFINITY;
\r
164 DPJ(1, 0) = MINUS_INFINITY;
\r
166 DPM(0, 1) = MINUS_INFINITY;
\r
167 DPD(0, 1) = MINUS_INFINITY;
\r
168 DPE(0, 1) = MINUS_INFINITY;
\r
169 DPI(0, 1) = PB[0].m_scoreGapOpen;
\r
170 DPJ(0, 1) = PB[0].m_scoreGapOpen2;
\r
172 // Empty prefix of B is special case
\r
173 for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
175 // M=LetterA+LetterB, impossible with empty prefix
\r
176 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
179 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;
\r
180 TBD(uPrefixLengthA, 0) = 'D';
\r
182 DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;
\r
183 TBE(uPrefixLengthA, 0) = 'E';
\r
185 // I=GapA+LetterB, impossible with empty prefix
\r
186 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
187 DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
190 // Empty prefix of A is special case
\r
191 for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
193 // M=LetterA+LetterB, impossible with empty prefix
\r
194 DPM(0, uPrefixLengthB) = MINUS_INFINITY;
\r
196 // D=LetterA+GapB, impossible with empty prefix
\r
197 DPD(0, uPrefixLengthB) = MINUS_INFINITY;
\r
198 DPE(0, uPrefixLengthB) = MINUS_INFINITY;
\r
201 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;
\r
202 TBI(0, uPrefixLengthB) = 'I';
\r
204 DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;
\r
205 TBJ(0, uPrefixLengthB) = 'J';
\r
211 for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
213 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
214 SCORE scoreGapCloseB;
\r
215 if (uPrefixLengthB == 1)
\r
216 scoreGapCloseB = MINUS_INFINITY;
\r
218 scoreGapCloseB = PB[uPrefixLengthB-2].m_scoreGapClose;
\r
220 SCORE scoreGapClose2B;
\r
221 if (uPrefixLengthB == 1)
\r
222 scoreGapClose2B = MINUS_INFINITY;
\r
224 scoreGapClose2B = PB[uPrefixLengthB-2].m_scoreGapClose2;
\r
226 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
228 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
231 // Match M=LetterA+LetterB
\r
232 SCORE scoreLL = ScoreProfPos2(PPA, PPB);
\r
233 DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;
\r
235 SCORE scoreGapCloseA;
\r
236 if (uPrefixLengthA == 1)
\r
237 scoreGapCloseA = MINUS_INFINITY;
\r
239 scoreGapCloseA = PA[uPrefixLengthA-2].m_scoreGapClose;
\r
241 SCORE scoreGapClose2A;
\r
242 if (uPrefixLengthA == 1)
\r
243 scoreGapClose2A = MINUS_INFINITY;
\r
245 scoreGapClose2A = PA[uPrefixLengthA-2].m_scoreGapClose2;
\r
247 SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
\r
248 SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
\r
249 SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;
\r
250 SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;
\r
251 SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;
\r
253 if (scoreMM >= scoreDM && scoreMM >= scoreIM && scoreMM >= scoreEM && scoreMM >= scoreJM)
\r
255 scoreBest = scoreMM;
\r
256 TBM(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
258 else if (scoreDM >= scoreMM && scoreDM >= scoreIM && scoreDM >= scoreEM && scoreDM >= scoreJM)
\r
260 scoreBest = scoreDM;
\r
261 TBM(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
263 else if (scoreEM >= scoreMM && scoreEM >= scoreIM && scoreEM >= scoreDM && scoreEM >= scoreJM)
\r
265 scoreBest = scoreEM;
\r
266 TBM(uPrefixLengthA, uPrefixLengthB) = 'E';
\r
268 else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)
\r
270 scoreBest = scoreIM;
\r
271 TBM(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
273 else if (scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM)
\r
275 scoreBest = scoreJM;
\r
276 TBM(uPrefixLengthA, uPrefixLengthB) = 'J';
\r
279 Quit("Max failed (M)");
\r
281 DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;
\r
285 // Delete D=LetterA+GapB
\r
286 SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
287 PA[uPrefixLengthA-1].m_scoreGapOpen;
\r
288 SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) +
\r
292 if (scoreMD >= scoreDD)
\r
294 scoreBest = scoreMD;
\r
295 TBD(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
299 assert(scoreDD >= scoreMD);
\r
300 scoreBest = scoreDD;
\r
301 TBD(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
303 DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
307 // Delete E=LetterA+GapB
\r
308 SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
309 PA[uPrefixLengthA-1].m_scoreGapOpen2;
\r
310 SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) +
\r
314 if (scoreME >= scoreEE)
\r
316 scoreBest = scoreME;
\r
317 TBE(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
321 assert(scoreEE >= scoreME);
\r
322 scoreBest = scoreEE;
\r
323 TBE(uPrefixLengthA, uPrefixLengthB) = 'E';
\r
325 DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
328 // Insert I=GapA+LetterB
\r
330 SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
331 PB[uPrefixLengthB-1].m_scoreGapOpen;
\r
332 SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) +
\r
336 if (scoreMI >= scoreII)
\r
338 scoreBest = scoreMI;
\r
339 TBI(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
343 assert(scoreII > scoreMI);
\r
344 scoreBest = scoreII;
\r
345 TBI(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
347 DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
350 // Insert J=GapA+LetterB
\r
352 SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
353 PB[uPrefixLengthB-1].m_scoreGapOpen2;
\r
354 SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) +
\r
358 if (scoreMJ > scoreJJ)
\r
360 scoreBest = scoreMJ;
\r
361 TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
365 assert(scoreJJ >= scoreMJ);
\r
366 scoreBest = scoreJJ;
\r
367 TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';
\r
369 DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
374 // Special case: close gaps at end of alignment
\r
375 DPD(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose;
\r
376 DPE(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose2;
\r
378 DPI(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose;
\r
379 DPJ(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose2;
\r
383 ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
385 ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
387 ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
389 ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
391 ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
393 ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
395 ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
397 ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
399 ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
401 ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
403 ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
413 char cEdgeType = '?';
\r
414 SCORE BestScore = MINUS_INFINITY;
\r
415 SCORE M = DPM(uLengthA, uLengthB);
\r
416 SCORE D = DPD(uLengthA, uLengthB);
\r
417 SCORE E = DPE(uLengthA, uLengthB);
\r
418 SCORE I = DPI(uLengthA, uLengthB);
\r
419 SCORE J = DPJ(uLengthA, uLengthB);
\r
421 if (M >= D && M >= E && M >= I && M >= J)
\r
426 else if (D >= M && D >= E && D >= I && D >= J)
\r
431 else if (E >= M && E >= D && E >= I && E >= J)
\r
436 else if (I >= M && I >= D && I >= E && I >= J)
\r
441 else if (J >= M && J >= D && J >= E && J >= I)
\r
449 unsigned PLA = uLengthA;
\r
450 unsigned PLB = uLengthB;
\r
451 unsigned ECount = 0;
\r
452 unsigned JCount = 0;
\r
456 Log("TraceBack: %c%u.%u\n", cEdgeType, PLA, PLB);
\r
459 Edge.cType = XlatEdgeType(cEdgeType);
\r
460 Edge.uPrefixLengthA = PLA;
\r
461 Edge.uPrefixLengthB = PLB;
\r
462 Path.PrependEdge(Edge);
\r
469 cEdgeType = TBM(PLA, PLB);
\r
476 cEdgeType = TBD(PLA, PLB);
\r
483 cEdgeType = TBE(PLA, PLB);
\r
489 cEdgeType = TBI(PLA, PLB);
\r
496 cEdgeType = TBJ(PLA, PLB);
\r
501 Quit("Invalid edge %c", cEdgeType);
\r
503 if (0 == PLA && 0 == PLB)
\r
506 //if (ECount > 0 || JCount > 0)
\r
507 // fprintf(stderr, "E=%d J=%d\n", ECount, JCount);
\r
509 if (Path.GetMatchCount() + Path.GetDeleteCount() != uLengthA)
\r
510 Quit("Path count A");
\r
511 if (Path.GetMatchCount() + Path.GetInsertCount() != uLengthB)
\r
512 Quit("Path count B");
\r
514 if (g_bKeepSimpleDP)
\r
544 Log("BestScore=%.6g\n", BestScore);
\r
549 #endif // DOUBLE_AFFINE
\r