9 bool g_bKeepSimpleDP;
\r
23 static char XlatEdgeType(char c)
\r
32 static const char *LocalScoreToStr(SCORE s)
\r
34 static char str[16];
\r
37 sprintf(str, "%6.1f", s);
\r
41 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,
\r
42 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
45 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
48 if (uPrefixLengthB > 0)
\r
49 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
50 Log(" %4u:%c", uPrefixLengthB, c);
\r
53 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
56 if (uPrefixLengthA > 0)
\r
57 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
58 Log("%4u:%c ", uPrefixLengthA, c);
\r
59 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
60 Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));
\r
65 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
\r
66 unsigned uPrefixCountA, unsigned uPrefixCountB)
\r
69 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
72 if (uPrefixLengthB > 0)
\r
73 c = ConsensusChar(PB[uPrefixLengthB - 1]);
\r
74 Log(" %4u:%c", uPrefixLengthB, c);
\r
77 for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
80 if (uPrefixLengthA > 0)
\r
81 c = ConsensusChar(PA[uPrefixLengthA - 1]);
\r
82 Log("%4u:%c ", uPrefixLengthA, c);
\r
83 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
84 Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
\r
89 SCORE NWDASimple(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
90 unsigned uLengthB, PWPath &Path)
\r
92 assert(uLengthB > 0 && uLengthA > 0);
\r
94 const unsigned uPrefixCountA = uLengthA + 1;
\r
95 const unsigned uPrefixCountB = uLengthB + 1;
\r
97 // Allocate DP matrices
\r
98 const size_t LM = uPrefixCountA*uPrefixCountB;
\r
99 SCORE *DPL_ = new SCORE[LM];
\r
100 SCORE *DPM_ = new SCORE[LM];
\r
101 SCORE *DPD_ = new SCORE[LM];
\r
102 SCORE *DPE_ = new SCORE[LM];
\r
103 SCORE *DPI_ = new SCORE[LM];
\r
104 SCORE *DPJ_ = new SCORE[LM];
\r
106 char *TBM_ = new char[LM];
\r
107 char *TBD_ = new char[LM];
\r
108 char *TBE_ = new char[LM];
\r
109 char *TBI_ = new char[LM];
\r
110 char *TBJ_ = new char[LM];
\r
112 memset(TBM_, '?', LM);
\r
113 memset(TBD_, '?', LM);
\r
114 memset(TBE_, '?', LM);
\r
115 memset(TBI_, '?', LM);
\r
116 memset(TBJ_, '?', LM);
\r
119 DPD(0, 0) = MINUS_INFINITY;
\r
120 DPE(0, 0) = MINUS_INFINITY;
\r
121 DPI(0, 0) = MINUS_INFINITY;
\r
122 DPJ(0, 0) = MINUS_INFINITY;
\r
124 DPM(1, 0) = MINUS_INFINITY;
\r
125 DPD(1, 0) = PA[0].m_scoreGapOpen;
\r
126 DPE(1, 0) = PA[0].m_scoreGapOpen2;
\r
129 DPI(1, 0) = MINUS_INFINITY;
\r
130 DPJ(1, 0) = MINUS_INFINITY;
\r
132 DPM(0, 1) = MINUS_INFINITY;
\r
133 DPD(0, 1) = MINUS_INFINITY;
\r
134 DPE(0, 1) = MINUS_INFINITY;
\r
135 DPI(0, 1) = PB[0].m_scoreGapOpen;
\r
136 DPJ(0, 1) = PB[0].m_scoreGapOpen2;
\r
140 // Empty prefix of B is special case
\r
141 for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
143 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
145 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;
\r
146 DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;
\r
148 TBD(uPrefixLengthA, 0) = 'D';
\r
149 TBE(uPrefixLengthA, 0) = 'E';
\r
151 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
152 DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;
\r
155 // Empty prefix of A is special case
\r
156 for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
158 DPM(0, uPrefixLengthB) = MINUS_INFINITY;
\r
160 DPD(0, uPrefixLengthB) = MINUS_INFINITY;
\r
161 DPE(0, uPrefixLengthB) = MINUS_INFINITY;
\r
163 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;
\r
164 DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;
\r
166 TBI(0, uPrefixLengthB) = 'I';
\r
167 TBJ(0, uPrefixLengthB) = 'J';
\r
170 // Special case to agree with NWFast, no D-I transitions so...
\r
171 DPD(uLengthA, 0) = MINUS_INFINITY;
\r
172 DPE(uLengthA, 0) = MINUS_INFINITY;
\r
173 // DPI(0, uLengthB) = MINUS_INFINITY;
\r
174 // DPJ(0, uLengthB) = MINUS_INFINITY;
\r
179 SCORE scoreGapCloseB = MINUS_INFINITY;
\r
180 SCORE scoreGapClose2B = MINUS_INFINITY;
\r
181 for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
\r
183 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
185 SCORE scoreGapCloseA = MINUS_INFINITY;
\r
186 SCORE scoreGapClose2A = MINUS_INFINITY;
\r
187 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
\r
189 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
192 // Match M=LetterA+LetterB
\r
193 SCORE scoreLL = ScoreProfPos2(PPA, PPB);
\r
194 DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;
\r
196 SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
\r
197 SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;
\r
198 SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;
\r
199 SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;
\r
200 SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;
\r
203 if (scoreMM >= scoreDM && scoreMM >= scoreEM && scoreMM >= scoreIM && scoreMM >= scoreJM)
\r
205 scoreBest = scoreMM;
\r
206 TBM(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
208 else if (scoreDM >= scoreMM && scoreDM >= scoreEM && scoreDM >= scoreIM && scoreDM >= scoreJM)
\r
210 scoreBest = scoreDM;
\r
211 TBM(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
213 else if (scoreEM >= scoreMM && scoreEM >= scoreDM && scoreEM >= scoreIM && scoreEM >= scoreJM)
\r
215 scoreBest = scoreEM;
\r
216 TBM(uPrefixLengthA, uPrefixLengthB) = 'E';
\r
218 else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)
\r
220 scoreBest = scoreIM;
\r
221 TBM(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
225 assert(scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM);
\r
226 scoreBest = scoreJM;
\r
227 TBM(uPrefixLengthA, uPrefixLengthB) = 'J';
\r
229 DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;
\r
233 // Delete D=LetterA+GapB
\r
234 SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
235 PA[uPrefixLengthA-1].m_scoreGapOpen;
\r
236 SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend;
\r
239 if (scoreMD >= scoreDD)
\r
241 scoreBest = scoreMD;
\r
242 TBD(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
246 assert(scoreDD >= scoreMD);
\r
247 scoreBest = scoreDD;
\r
248 TBD(uPrefixLengthA, uPrefixLengthB) = 'D';
\r
250 DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
254 // Delete E=LetterA+GapB
\r
255 SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +
\r
256 PA[uPrefixLengthA-1].m_scoreGapOpen2;
\r
257 SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) + g_scoreGapExtend2;
\r
260 if (scoreME >= scoreEE)
\r
262 scoreBest = scoreME;
\r
263 TBE(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
267 assert(scoreEE >= scoreME);
\r
268 scoreBest = scoreEE;
\r
269 TBE(uPrefixLengthA, uPrefixLengthB) = 'E';
\r
271 DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
274 // Insert I=GapA+LetterB
\r
276 SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
277 PB[uPrefixLengthB - 1].m_scoreGapOpen;
\r
278 SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend;
\r
281 if (scoreMI >= scoreII)
\r
283 scoreBest = scoreMI;
\r
284 TBI(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
288 assert(scoreII > scoreMI);
\r
289 scoreBest = scoreII;
\r
290 TBI(uPrefixLengthA, uPrefixLengthB) = 'I';
\r
292 DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
295 // Insert J=GapA+LetterB
\r
297 SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +
\r
298 PB[uPrefixLengthB - 1].m_scoreGapOpen2;
\r
299 SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) + g_scoreGapExtend2;
\r
302 if (scoreMJ >= scoreJJ)
\r
304 scoreBest = scoreMJ;
\r
305 TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';
\r
309 assert(scoreJJ > scoreMJ);
\r
310 scoreBest = scoreJJ;
\r
311 TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';
\r
313 DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;
\r
316 scoreGapCloseA = PPA.m_scoreGapClose;
\r
317 scoreGapClose2A = PPA.m_scoreGapClose2;
\r
319 scoreGapCloseB = PPB.m_scoreGapClose;
\r
320 scoreGapClose2B = PPB.m_scoreGapClose2;
\r
325 Log("DA Simple DPL:\n");
\r
326 ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
328 Log("DA Simple DPM:\n");
\r
329 ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
331 Log("DA Simple DPD:\n");
\r
332 ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
334 Log("DA Simple DPE:\n");
\r
335 ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
337 Log("DA Simple DPI:\n");
\r
338 ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
340 Log("DA Simple DPJ:\n");
\r
341 ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
343 Log("DA Simple TBM:\n");
\r
344 ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
346 Log("DA Simple TBD:\n");
\r
347 ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
349 Log("DA Simple TBE:\n");
\r
350 ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
352 Log("DA Simple TBI:\n");
\r
353 ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
355 Log("DA Simple TBJ:\n");
\r
356 ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);
\r
364 SCORE M = DPM(uLengthA, uLengthB);
\r
365 SCORE D = DPD(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose;
\r
366 SCORE E = DPE(uLengthA, uLengthB) + PA[uLengthA-1].m_scoreGapClose2;
\r
367 SCORE I = DPI(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose;
\r
368 SCORE J = DPJ(uLengthA, uLengthB) + PB[uLengthB-1].m_scoreGapClose2;
\r
369 char cEdgeType = '?';
\r
371 SCORE BestScore = M;
\r
395 Log("DA Simple: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",
\r
396 M, D, E, I, J, cEdgeType);
\r
399 unsigned PLA = uLengthA;
\r
400 unsigned PLB = uLengthB;
\r
404 Edge.cType = XlatEdgeType(cEdgeType);
\r
405 Edge.uPrefixLengthA = PLA;
\r
406 Edge.uPrefixLengthB = PLB;
\r
408 Log("Prepend %c%d.%d\n", Edge.cType, PLA, PLB);
\r
410 Path.PrependEdge(Edge);
\r
417 cEdgeType = TBM(PLA, PLB);
\r
424 cEdgeType = TBD(PLA, PLB);
\r
430 cEdgeType = TBE(PLA, PLB);
\r
436 cEdgeType = TBI(PLA, PLB);
\r
442 cEdgeType = TBJ(PLA, PLB);
\r
447 Quit("Invalid edge %c", cEdgeType);
\r
449 if (0 == PLA && 0 == PLB)
\r
454 // SCORE Score = TraceBack(PA, uLengthA, PB, uLengthB, DPM_, DPD_, DPI_, Path);
\r
457 SCORE scorePath = FastScorePath2(PA, uLengthA, PB, uLengthB, Path);
\r
459 Log("Score = %s Path = %s\n", LocalScoreToStr(BestScore), LocalScoreToStr(scorePath));
\r
462 if (g_bKeepSimpleDP)
\r
494 #endif // DOUBLE_AFFINE
\r