source: trunk/GDE/MUSCLE/src/nwsmall.cpp

Last change on this file was 10390, checked in by aboeckma, 11 years ago

added muscle sourcles amd makefile

File size: 15.8 KB
Line 
1#include "muscle.h"
2#include <math.h>
3#include "pwpath.h"
4#include "profile.h"
5#include <stdio.h>
6
7// NW small memory
8
9#define TRACE   0
10
11#if     TRACE
12extern bool g_bKeepSimpleDP;
13extern SCORE *g_DPM;
14extern SCORE *g_DPD;
15extern SCORE *g_DPI;
16extern char *g_TBM;
17extern char *g_TBD;
18extern char *g_TBI;
19#endif
20
21#if     TRACE
22#define ALLOC_TRACE()                                                           \
23        const SCORE UNINIT = MINUS_INFINITY;                    \
24        const size_t LM = uPrefixCountA*uPrefixCountB;  \
25                                                                                                        \
26        SCORE *DPM_ = new SCORE[LM];                                    \
27        SCORE *DPD_ = new SCORE[LM];                                    \
28        SCORE *DPI_ = new SCORE[LM];                                    \
29                                                                                                        \
30        char *TBM_ = new char[LM];                                              \
31        char *TBD_ = new char[LM];                                              \
32        char *TBI_ = new char[LM];                                              \
33                                                                                                        \
34        memset(TBM_, '?', LM);                                                  \
35        memset(TBD_, '?', LM);                                                  \
36        memset(TBI_, '?', LM);                                                  \
37                                                                                                        \
38        for (unsigned i = 0; i <= uLengthA; ++i)                \
39                for (unsigned j = 0; j <= uLengthB; ++j)        \
40                        {                                                                               \
41                        DPM(i, j) = UNINIT;                                             \
42                        DPD(i, j) = UNINIT;                                             \
43                        DPI(i, j) = UNINIT;                                             \
44                        }
45#else
46#define ALLOC_TRACE()
47#endif
48
49#if     TRACE
50#define SetDPM(i, j, x)         DPM(i, j) = x
51#define SetDPD(i, j, x)         DPD(i, j) = x
52#define SetDPI(i, j, x)         DPI(i, j) = x
53#define SetTBM(i, j, x)         TBM(i, j) = x
54#define SetTBD(i, j, x)         TBD(i, j) = x
55#define SetTBI(i, j, x)         TBI(i, j) = x
56#else
57#define SetDPM(i, j, x)         /* empty  */
58#define SetDPD(i, j, x)         /* empty  */
59#define SetDPI(i, j, x)         /* empty  */
60#define SetTBM(i, j, x)         /* empty  */
61#define SetTBD(i, j, x)         /* empty  */
62#define SetTBI(i, j, x)         /* empty  */
63#endif
64
65#define RECURSE_D(i, j)                         \
66        {                                                               \
67        SCORE DD = DRow[j] + e;                 \
68        SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\
69        if (DD > MD)                                    \
70                {                                                       \
71                DRow[j] = DD;                           \
72                SetTBD(i, j, 'D');                      \
73                }                                                       \
74        else                                                    \
75                {                                                       \
76                DRow[j] = MD;                           \
77                /* SetBitTBD(TB, i, j, 'M'); */ \
78                TBRow[j] &= ~BIT_xD;            \
79                TBRow[j] |= BIT_MD;                     \
80                SetTBD(i, j, 'M');                      \
81                }                                                       \
82        SetDPD(i, j, DRow[j]);                  \
83        }
84
85#define RECURSE_D_ATerm(j)      RECURSE_D(uLengthA, j)
86#define RECURSE_D_BTerm(j)      RECURSE_D(i, uLengthB)
87
88#define RECURSE_I(i, j)                         \
89        {                                                               \
90        Iij += e;                                               \
91        SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\
92        if (MI >= Iij)                                  \
93                {                                                       \
94                Iij = MI;                                       \
95                /* SetBitTBI(TB, i, j, 'M'); */ \
96                TBRow[j] &= ~BIT_xI;            \
97                TBRow[j] |= BIT_MI;                     \
98                SetTBI(i, j, 'M');                      \
99                }                                                       \
100        else                                                    \
101                SetTBI(i, j, 'I');                      \
102        SetDPI(i, j, Iij);                              \
103        }
104
105#define RECURSE_I_ATerm(j)      RECURSE_I(uLengthA, j)
106#define RECURSE_I_BTerm(j)      RECURSE_I(i, uLengthB)
107
108#define RECURSE_M(i, j)                                                         \
109        {                                                                                               \
110        SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose;   \
111        SCORE IM = Iij +     PB[j-1].m_scoreGapClose;   \
112        SCORE MM = MCurr[j];                                                    \
113        TB[i+1][j+1] &= ~BIT_xM;                                                        \
114        if (MM >= DM && MM >= IM)                                               \
115                {                                                                                       \
116                MNext[j+1] += MM;                                                       \
117                SetDPM(i+1, j+1, MNext[j+1]);                           \
118                SetTBM(i+1, j+1, 'M');                                          \
119                /* SetBitTBM(TB, i+1, j+1, 'M');        */              \
120                TB[i+1][j+1] |= BIT_MM;                                         \
121                }                                                                                       \
122        else if (DM >= MM && DM >= IM)                                  \
123                {                                                                                       \
124                MNext[j+1] += DM;                                                       \
125                SetDPM(i+1, j+1, MNext[j+1]);                           \
126                SetTBM(i+1, j+1, 'D');                                          \
127                /* SetBitTBM(TB, i+1, j+1, 'D'); */                     \
128                TB[i+1][j+1] |= BIT_DM;                                         \
129                }                                                                                       \
130        else                                                                                    \
131                {                                                                                       \
132                assert(IM >= MM && IM >= DM);                           \
133                MNext[j+1] += IM;                                                       \
134                SetDPM(i+1, j+1, MNext[j+1]);                           \
135                SetTBM(i+1, j+1, 'I');                                          \
136                /* SetBitTBM(TB, i+1, j+1, 'I'); */                     \
137                TB[i+1][j+1] |= BIT_IM;                                         \
138                }                                                                                       \
139        }
140
141#if     TRACE
142static bool LocalEq(BASETYPE b1, BASETYPE b2)
143        {
144        if (b1 < -100000 && b2 < -100000)
145                return true;
146        double diff = fabs(b1 - b2);
147        if (diff < 0.0001)
148                return true;
149        double sum = fabs(b1) + fabs(b2);
150        return diff/sum < 0.005;
151        }
152
153static char Get_M_Char(char Bits)
154        {
155        switch (Bits & BIT_xM)
156                {
157        case BIT_MM:
158                return 'M';
159        case BIT_DM:
160                return 'D';
161        case BIT_IM:
162                return 'I';
163                }
164        Quit("Huh?");
165        return '?';
166        }
167
168static char Get_D_Char(char Bits)
169        {
170        return (Bits & BIT_xD) ? 'M' : 'D';
171        }
172
173static char Get_I_Char(char Bits)
174        {
175        return (Bits & BIT_xI) ? 'M' : 'I';
176        }
177
178static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,
179  unsigned uPrefixCountA, unsigned uPrefixCountB)
180        {
181        SCORE *DPM_ = g_DP;
182        for (unsigned i = 0; i < uPrefixCountA; ++i)
183                for (unsigned j = 0; j < uPrefixCountB; ++j)
184                        if (!LocalEq(DPM(i, j), DPD(i, j)))
185                                {
186                                Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Fast = %.2g\n",
187                                  c, i, j, DPM(i, j), DPD(i, j));
188                                return false;
189                                }
190        return true;
191        }
192
193static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBI_, 
194  unsigned uPrefixCountA, unsigned uPrefixCountB)
195        {
196        SCORE *DPM_ = g_DPM;
197        bool Eq = true;
198        for (unsigned i = 0; i < uPrefixCountA; ++i)
199                for (unsigned j = 0; j < uPrefixCountB; ++j)
200                        {
201                        char c1 = TBM(i, j);
202                        char c2 = Get_M_Char(TB[i][j]);
203                        if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)
204                                {
205                                Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
206                                Eq = false;
207                                goto D;
208                                }
209                        }
210
211D:
212        SCORE *DPD_ = g_DPD;
213        for (unsigned i = 0; i < uPrefixCountA; ++i)
214                for (unsigned j = 0; j < uPrefixCountB; ++j)
215                        {
216                        char c1 = TBD(i, j);
217                        char c2 = Get_D_Char(TB[i][j]);
218                        if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)
219                                {
220                                Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
221                                Eq = false;
222                                goto I;
223                                }
224                        }
225I:
226        SCORE *DPI_ = g_DPI;
227        for (unsigned i = 0; i < uPrefixCountA; ++i)
228                for (unsigned j = 0; j < uPrefixCountB; ++j)
229                        {
230                        char c1 = TBI(i, j);
231                        char c2 = Get_I_Char(TB[i][j]);
232                        if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)
233                                {
234                                Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
235                                Eq = false;
236                                goto Done;
237                                }
238                        }
239Done:
240        if (Eq)
241                Log("TB success\n");
242        return Eq;
243        }
244
245static const char *LocalScoreToStr(SCORE s)
246        {
247        static char str[16];
248        if (s < -100000)
249                return "     *";
250        sprintf(str, "%6.1f", s);
251        return str;
252        }
253
254static void LogDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
255  unsigned uPrefixCountA, unsigned uPrefixCountB)
256        {
257        Log("        ");
258        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
259                {
260                char c = ' ';
261                if (uPrefixLengthB > 0)
262                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
263                Log(" %4u:%c", uPrefixLengthB, c);
264                }
265        Log("\n");
266        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
267                {
268                char c = ' ';
269                if (uPrefixLengthA > 0)
270                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
271                Log("%4u:%c  ", uPrefixLengthA, c);
272                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
273                        Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
274                Log("\n");
275                }
276        }
277
278static void LogBitTB(char **TB, const ProfPos *PA, const ProfPos *PB,
279  unsigned uPrefixCountA, unsigned uPrefixCountB)
280        {
281        Log("        ");
282        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
283                {
284                char c = ' ';
285                if (uPrefixLengthB > 0)
286                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
287                Log(" %4u:%c", uPrefixLengthB, c);
288                }
289        Log("\n");
290        Log("Bit TBM:\n");
291        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
292                {
293                char c = ' ';
294                if (uPrefixLengthA > 0)
295                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
296                Log("%4u:%c  ", uPrefixLengthA, c);
297                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
298                        {
299                        char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);
300                        Log(" %6c", c);
301                        }
302                Log("\n");
303                }
304
305        Log("\n");
306        Log("Bit TBD:\n");
307        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
308                {
309                char c = ' ';
310                if (uPrefixLengthA > 0)
311                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
312                Log("%4u:%c  ", uPrefixLengthA, c);
313                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
314                        {
315                        char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);
316                        Log(" %6c", c);
317                        }
318                Log("\n");
319                }
320
321        Log("\n");
322        Log("Bit TBI:\n");
323        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
324                {
325                char c = ' ';
326                if (uPrefixLengthA > 0)
327                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
328                Log("%4u:%c  ", uPrefixLengthA, c);
329                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
330                        {
331                        char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);
332                        Log(" %6c", c);
333                        }
334                Log("\n");
335                }
336        }
337
338static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
339  unsigned uPrefixCountA, unsigned uPrefixCountB)
340        {
341        Log("        ");
342        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
343                {
344                char c = ' ';
345                if (uPrefixLengthB > 0)
346                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
347                Log(" %4u:%c", uPrefixLengthB, c);
348                }
349        Log("\n");
350        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
351                {
352                char c = ' ';
353                if (uPrefixLengthA > 0)
354                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
355                Log("%4u:%c  ", uPrefixLengthA, c);
356                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
357                        {
358                        char c = TBM(uPrefixLengthA, uPrefixLengthB);
359                        Log(" %6c", c);
360                        }
361                Log("\n");
362                }
363        }
364
365static const char *BitsToStr(char Bits)
366        {
367        static char Str[9];
368
369        sprintf(Str, "%cM %cD %cI",
370          Get_M_Char(Bits),
371          Get_D_Char(Bits),
372          Get_I_Char(Bits));
373        }
374#endif  // TRACE
375
376static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
377        {
378        char Bit;
379        switch (c)
380                {
381        case 'M':
382                Bit = BIT_MM;
383                break;
384        case 'D':
385                Bit = BIT_DM;
386                break;
387        case 'I':
388                Bit = BIT_IM;
389                break;
390        default:
391                Quit("Huh?!");
392                }
393        TB[i][j] &= ~BIT_xM;
394        TB[i][j] |= Bit;
395        }
396
397static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
398        {
399        char Bit;
400        switch (c)
401                {
402        case 'M':
403                Bit = BIT_MD;
404                break;
405        case 'D':
406                Bit = BIT_DD;
407                break;
408        default:
409                Quit("Huh?!");
410                }
411        TB[i][j] &= ~BIT_xD;
412        TB[i][j] |= Bit;
413        }
414
415static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
416        {
417        char Bit;
418        switch (c)
419                {
420        case 'M':
421                Bit = BIT_MI;
422                break;
423        case 'I':
424                Bit = BIT_II;
425                break;
426        default:
427                Quit("Huh?!");
428                }
429        TB[i][j] &= ~BIT_xI;
430        TB[i][j] |= Bit;
431        }
432
433#if     TRACE
434#define LogMatrices()                                                                                   \
435        {                                                                                                                       \
436        Log("Bit DPM:\n");                                                                                      \
437        LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
438        Log("Bit DPD:\n");                                                                                      \
439        LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
440        Log("Bit DPI:\n");                                                                                      \
441        LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
442        Log("Bit TB:\n");                                                                                       \
443        LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB);                     \
444        bool Same;                                                                                                      \
445        Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
446        if (Same)                                                                                                       \
447                Log("DPM success\n");                                                                   \
448        Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
449        if (Same)                                                                                                       \
450                Log("DPD success\n");                                                                   \
451        Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
452        if (Same)                                                                                                       \
453                Log("DPI success\n");                                                                   \
454        CompareTB(TB, g_TBM, g_TBD, g_TBI, uPrefixCountA, uPrefixCountB);\
455        }
456#else
457#define LogMatrices()   /* empty */
458#endif
459
460static unsigned uCachePrefixCountB;
461static unsigned uCachePrefixCountA;
462static SCORE *CacheMCurr;
463static SCORE *CacheMNext;
464static SCORE *CacheMPrev;
465static SCORE *CacheDRow;
466static char **CacheTB;
467
468static void AllocCache(unsigned uPrefixCountA, unsigned uPrefixCountB)
469        {
470        if (uPrefixCountA <= uCachePrefixCountA && uPrefixCountB <= uCachePrefixCountB)
471                return;
472
473        delete[] CacheMCurr;
474        delete[] CacheMNext;
475        delete[] CacheMPrev;
476        delete[] CacheDRow;
477        for (unsigned i = 0; i < uCachePrefixCountA; ++i)
478                delete[] CacheTB[i];
479        delete[] CacheTB;
480
481        uCachePrefixCountA = uPrefixCountA + 1024;
482        uCachePrefixCountB = uPrefixCountB + 1024;
483
484        CacheMCurr = new SCORE[uCachePrefixCountB];
485        CacheMNext = new SCORE[uCachePrefixCountB];
486        CacheMPrev = new SCORE[uCachePrefixCountB];
487        CacheDRow = new SCORE[uCachePrefixCountB];
488
489        CacheTB = new char *[uCachePrefixCountA];
490        for (unsigned i = 0; i < uCachePrefixCountA; ++i)
491                CacheTB[i] = new char [uCachePrefixCountB];
492        }
493
494SCORE NWSmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
495  unsigned uLengthB, PWPath &Path)
496        {
497        if (0 == uLengthB || 0 == uLengthA )
498                Quit("Internal error, NWSmall: length=0");
499
500        SetTermGaps(PA, uLengthA);
501        SetTermGaps(PB, uLengthB);
502
503        const unsigned uPrefixCountA = uLengthA + 1;
504        const unsigned uPrefixCountB = uLengthB + 1;
505        const SCORE e = g_scoreGapExtend;
506
507        ALLOC_TRACE()
508
509        AllocCache(uPrefixCountA, uPrefixCountB);
510
511        SCORE *MCurr = CacheMCurr;
512        SCORE *MNext = CacheMNext;
513        SCORE *MPrev = CacheMPrev;
514        SCORE *DRow = CacheDRow;
515
516        char **TB = CacheTB;
517        for (unsigned i = 0; i < uPrefixCountA; ++i)
518                memset(TB[i], 0, uPrefixCountB);
519
520        SCORE Iij = MINUS_INFINITY;
521        SetDPI(0, 0, Iij);
522
523        Iij = PB[0].m_scoreGapOpen;
524        SetDPI(0, 1, Iij);
525
526        for (unsigned j = 2; j <= uLengthB; ++j)
527                {
528                Iij += e;
529                SetDPI(0, j, Iij);
530                SetTBI(0, j, 'I');
531                }
532
533        for (unsigned j = 0; j <= uLengthB; ++j)
534                {
535                DRow[j] = MINUS_INFINITY;
536                SetDPD(0, j, DRow[j]);
537                SetTBD(0, j, 'D');
538                }
539
540        MPrev[0] = 0;
541        SetDPM(0, 0, MPrev[0]);
542        for (unsigned j = 1; j <= uLengthB; ++j)
543                {
544                MPrev[j] = MINUS_INFINITY;
545                SetDPM(0, j, MPrev[j]);
546                }
547
548        MCurr[0] = MINUS_INFINITY;
549        SetDPM(1, 0, MCurr[0]);
550
551        MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
552        SetDPM(1, 1, MCurr[1]);
553        SetBitTBM(TB, 1, 1, 'M');
554        SetTBM(1, 1, 'M');
555
556        for (unsigned j = 2; j <= uLengthB; ++j)
557                {
558                MCurr[j] = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
559                  (j - 2)*e + PB[j-2].m_scoreGapClose;
560                SetDPM(1, j, MCurr[j]);
561                SetBitTBM(TB, 1, j, 'I');
562                SetTBM(1, j, 'I');
563                }
564
565// Main DP loop
566        for (unsigned i = 1; i < uLengthA; ++i)
567                {
568                char *TBRow = TB[i];
569
570                Iij = MINUS_INFINITY;
571                SetDPI(i, 0, Iij);
572
573                DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
574                SetDPD(i, 0, DRow[0]);
575
576                MCurr[0] = MINUS_INFINITY; 
577                if (i == 1)
578                        {
579                        MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
580                        SetBitTBM(TB, i, 1, 'M');
581                        SetTBM(i, 1, 'M');
582                        }
583                else
584                        {
585                        MCurr[1] = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
586                          (i - 2)*e + PA[i-2].m_scoreGapClose;
587                        SetBitTBM(TB, i, 1, 'D');
588                        SetTBM(i, 1, 'D');
589                        }
590                SetDPM(i, 0, MCurr[0]);
591                SetDPM(i, 1, MCurr[1]);
592
593                for (unsigned j = 1; j < uLengthB; ++j)
594                        MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);
595
596                for (unsigned j = 1; j < uLengthB; ++j)
597                        {
598                        RECURSE_D(i, j)
599                        RECURSE_I(i, j)
600                        RECURSE_M(i, j)
601                        }
602        // Special case for j=uLengthB
603                RECURSE_D_BTerm(i)
604                RECURSE_I_BTerm(i)
605
606        // Prev := Curr, Curr := Next, Next := Prev
607                Rotate(MPrev, MCurr, MNext);
608                }
609
610// Special case for i=uLengthA
611        char *TBRow = TB[uLengthA];
612        MCurr[0] = MINUS_INFINITY;
613        if (uLengthA > 1)
614                MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
615                  PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
616        else
617                MCurr[1] = ScoreProfPos2(PA[uLengthA-1], PB[0]) + PA[0].m_scoreGapOpen +
618                  PA[0].m_scoreGapClose;
619        SetBitTBM(TB, uLengthA, 1, 'D');
620        SetTBM(uLengthA, 1, 'D');
621        SetDPM(uLengthA, 0, MCurr[0]);
622        SetDPM(uLengthA, 1, MCurr[1]);
623
624        DRow[0] = MINUS_INFINITY;
625        SetDPD(uLengthA, 0, DRow[0]);
626        for (unsigned j = 1; j <= uLengthB; ++j)
627                RECURSE_D_ATerm(j);
628
629        Iij = MINUS_INFINITY;
630        for (unsigned j = 1; j <= uLengthB; ++j)
631                RECURSE_I_ATerm(j)
632
633        LogMatrices();
634
635        SCORE MAB = MCurr[uLengthB];
636        SCORE DAB = DRow[uLengthB];
637        SCORE IAB = Iij;
638
639        SCORE Score = MAB;
640        char cEdgeType = 'M';
641        if (DAB > Score)
642                {
643                Score = DAB;
644                cEdgeType = 'D';
645                }
646        if (IAB > Score)
647                {
648                Score = IAB;
649                cEdgeType = 'I';
650                }
651
652#if TRACE
653        Log("    Fast: MAB=%.4g DAB=%.4g IAB=%.4g best=%c\n",
654          MAB, DAB, IAB, cEdgeType);
655#endif
656
657        BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);
658
659#if     DBEUG
660        Path.Validate();
661#endif
662
663        return 0;
664        }
Note: See TracBrowser for help on using the repository browser.