source: branches/lib/GDE/MUSCLE/src/nwdasmall.cpp

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

added muscle sourcles amd makefile
File size: 21.9 KB
Line 
1#include "muscle.h"
2#include <math.h>
3#include "pwpath.h"
4#include "profile.h"
5#include <stdio.h>
6
7#if     DOUBLE_AFFINE
8
9// NW double affine small memory, term gaps fully penalized
10// (so up to caller to adjust in profile if desired).
11
12#define TRACE   0
13
14#define MIN(x, y)       ((x) < (y) ? (x) : (y))
15
16#if     TRACE
17extern bool g_bKeepSimpleDP;
18extern SCORE *g_DPM;
19extern SCORE *g_DPD;
20extern SCORE *g_DPE;
21extern SCORE *g_DPI;
22extern SCORE *g_DPJ;
23extern char *g_TBM;
24extern char *g_TBD;
25extern char *g_TBE;
26extern char *g_TBI;
27extern char *g_TBJ;
28#endif
29
30#if     TRACE
31#define ALLOC_TRACE()                                                           \
32        const SCORE UNINIT = MINUS_INFINITY;                    \
33        const size_t LM = uPrefixCountA*uPrefixCountB;  \
34                                                                                                        \
35        SCORE *DPM_ = new SCORE[LM];                                    \
36        SCORE *DPD_ = new SCORE[LM];                                    \
37        SCORE *DPE_ = new SCORE[LM];                                    \
38        SCORE *DPI_ = new SCORE[LM];                                    \
39        SCORE *DPJ_ = new SCORE[LM];                                    \
40                                                                                                        \
41        char *TBM_ = new char[LM];                                              \
42        char *TBD_ = new char[LM];                                              \
43        char *TBE_ = new char[LM];                                              \
44        char *TBI_ = new char[LM];                                              \
45        char *TBJ_ = new char[LM];                                              \
46                                                                                                        \
47        memset(TBM_, '?', LM);                                                  \
48        memset(TBD_, '?', LM);                                                  \
49        memset(TBE_, '?', LM);                                                  \
50        memset(TBI_, '?', LM);                                                  \
51        memset(TBJ_, '?', LM);                                                  \
52                                                                                                        \
53        for (unsigned i = 0; i <= uLengthA; ++i)                \
54                for (unsigned j = 0; j <= uLengthB; ++j)        \
55                        {                                                                               \
56                        DPM(i, j) = UNINIT;                                             \
57                        DPD(i, j) = UNINIT;                                             \
58                        DPE(i, j) = UNINIT;                                             \
59                        DPI(i, j) = UNINIT;                                             \
60                        DPJ(i, j) = UNINIT;                                             \
61                        }
62#else
63#define ALLOC_TRACE()
64#endif
65
66#if     TRACE
67#define SetDPM(i, j, x)         DPM(i, j) = x
68#define SetDPD(i, j, x)         DPD(i, j) = x
69#define SetDPE(i, j, x)         DPE(i, j) = x
70#define SetDPI(i, j, x)         DPI(i, j) = x
71#define SetDPJ(i, j, x)         DPJ(i, j) = x
72#define SetTBM(i, j, x)         TBM(i, j) = x
73#define SetTBD(i, j, x)         TBD(i, j) = x
74#define SetTBE(i, j, x)         TBE(i, j) = x
75#define SetTBI(i, j, x)         TBI(i, j) = x
76#define SetTBJ(i, j, x)         TBJ(i, j) = x
77#else
78#define SetDPM(i, j, x)         /* empty  */
79#define SetDPD(i, j, x)         /* empty  */
80#define SetDPE(i, j, x)         /* empty  */
81#define SetDPI(i, j, x)         /* empty  */
82#define SetDPJ(i, j, x)         /* empty  */
83#define SetTBM(i, j, x)         /* empty  */
84#define SetTBD(i, j, x)         /* empty  */
85#define SetTBE(i, j, x)         /* empty  */
86#define SetTBI(i, j, x)         /* empty  */
87#define SetTBJ(i, j, x)         /* empty  */
88#endif
89
90#define RECURSE_D(i, j)                         \
91        {                                                               \
92        SCORE DD = DRow[j] + e;                 \
93        SCORE MD = MPrev[j] + PA[i-1].m_scoreGapOpen;\
94        if (DD > MD)                                    \
95                {                                                       \
96                DRow[j] = DD;                           \
97                SetTBD(i, j, 'D');                      \
98                }                                                       \
99        else                                                    \
100                {                                                       \
101                DRow[j] = MD;                           \
102                SetBitTBD(TB, i, j, 'M');       \
103                SetTBD(i, j, 'M');                      \
104                }                                                       \
105        SetDPD(i, j, DRow[j]);                  \
106        }
107
108#define RECURSE_E(i, j)                         \
109        {                                                               \
110        SCORE EE = ERow[j] + e2;                \
111        SCORE ME = MPrev[j] + PA[i-1].m_scoreGapOpen2;\
112        if (EE > ME)                                    \
113                {                                                       \
114                ERow[j] = EE;                           \
115                SetTBE(i, j, 'E');                      \
116                }                                                       \
117        else                                                    \
118                {                                                       \
119                ERow[j] = ME;                           \
120                SetBitTBE(TB, i, j, 'M');       \
121                SetTBE(i, j, 'M');                      \
122                }                                                       \
123        SetDPE(i, j, ERow[j]);                  \
124        }
125
126#define RECURSE_D_ATerm(j)      RECURSE_D(uLengthA, j)
127#define RECURSE_E_ATerm(j)      RECURSE_E(uLengthA, j)
128
129#define RECURSE_D_BTerm(j)      RECURSE_D(i, uLengthB)
130#define RECURSE_E_BTerm(j)      RECURSE_E(i, uLengthB)
131
132#define RECURSE_I(i, j)                         \
133        {                                                               \
134        Iij += e;                                               \
135        SCORE MI = MCurr[j-1] + PB[j-1].m_scoreGapOpen;\
136        if (MI >= Iij)                                  \
137                {                                                       \
138                Iij = MI;                                       \
139                SetBitTBI(TB, i, j, 'M');       \
140                SetTBI(i, j, 'M');                      \
141                }                                                       \
142        else                                                    \
143                SetTBI(i, j, 'I');                      \
144        SetDPI(i, j, Iij);                              \
145        }
146
147#define RECURSE_J(i, j)                         \
148        {                                                               \
149        Jij += e2;                                              \
150        SCORE MJ = MCurr[j-1] + PB[j-1].m_scoreGapOpen2;\
151        if (MJ >= Jij)                                  \
152                {                                                       \
153                Jij = MJ;                                       \
154                SetBitTBJ(TB, i, j, 'M');       \
155                SetTBJ(i, j, 'M');                      \
156                }                                                       \
157        else                                                    \
158                SetTBJ(i, j, 'I');                      \
159        SetDPJ(i, j, Jij);                              \
160        }
161
162#define RECURSE_I_ATerm(j)      RECURSE_I(uLengthA, j)
163#define RECURSE_J_ATerm(j)      RECURSE_J(uLengthA, j)
164
165#define RECURSE_I_BTerm(j)      RECURSE_I(i, uLengthB)
166#define RECURSE_J_BTerm(j)      RECURSE_J(i, uLengthB)
167
168#define RECURSE_M(i, j)                                                                 \
169        {                                                                                                       \
170        SCORE Best = MCurr[j];  /*  MM  */                                      \
171        SetTBM(i+1, j+1, 'M');                                                          \
172        SetBitTBM(TB, i+1, j+1, 'M');                                           \
173                                                                                                                \
174        SCORE DM = DRow[j] + PA[i-1].m_scoreGapClose;           \
175        if (DM > Best)                                                                          \
176                {                                                                                               \
177                Best = DM;                                                                              \
178                SetTBM(i+1, j+1, 'D');                                                  \
179                SetBitTBM(TB, i+1, j+1, 'D');                                   \
180                }                                                                                               \
181                                                                                                                \
182        SCORE EM = ERow[j] + PA[i-1].m_scoreGapClose2;          \
183        if (EM > Best)                                                                          \
184                {                                                                                               \
185                Best = EM;                                                                              \
186                SetTBM(i+1, j+1, 'E');                                                  \
187                SetBitTBM(TB, i+1, j+1, 'E');                                   \
188                }                                                                                               \
189                                                                                                                \
190        SCORE IM = Iij + PB[j-1].m_scoreGapClose;                       \
191        if (IM > Best)                                                                          \
192                {                                                                                               \
193                Best = IM;                                                                              \
194                SetTBM(i+1, j+1, 'I');                                                  \
195                SetBitTBM(TB, i+1, j+1, 'I');                                   \
196                }                                                                                               \
197                                                                                                                \
198        SCORE JM = Jij + PB[j-1].m_scoreGapClose2;                      \
199        if (JM > Best)                                                                          \
200                {                                                                                               \
201                Best = JM;                                                                              \
202                SetTBM(i+1, j+1, 'J');                                                  \
203                SetBitTBM(TB, i+1, j+1, 'J');                                   \
204                }                                                                                               \
205        MNext[j+1] += Best;                                                                     \
206        SetDPM(i+1, j+1, MNext[j+1]);                                           \
207        }
208
209#if     TRACE
210static bool LocalEq(BASETYPE b1, BASETYPE b2)
211        {
212        if (b1 < -100000 && b2 < -100000)
213                return true;
214        double diff = fabs(b1 - b2);
215        if (diff < 0.0001)
216                return true;
217        double sum = fabs(b1) + fabs(b2);
218        return diff/sum < 0.005;
219        }
220
221static char Get_M_Char(char Bits)
222        {
223        switch (Bits & BIT_xM)
224                {
225        case BIT_MM:
226                return 'M';
227        case BIT_DM:
228                return 'D';
229        case BIT_EM:
230                return 'E';
231        case BIT_IM:
232                return 'I';
233        case BIT_JM:
234                return 'J';
235                }
236        Quit("Huh?");
237        return '?';
238        }
239
240static char Get_D_Char(char Bits)
241        {
242        return (Bits & BIT_xD) ? 'M' : 'D';
243        }
244
245static char Get_E_Char(char Bits)
246        {
247        return (Bits & BIT_xE) ? 'M' : 'E';
248        }
249
250static char Get_I_Char(char Bits)
251        {
252        return (Bits & BIT_xI) ? 'M' : 'I';
253        }
254
255static char Get_J_Char(char Bits)
256        {
257        return (Bits & BIT_xJ) ? 'M' : 'J';
258        }
259
260static bool DPEq(char c, SCORE *g_DP, SCORE *DPD_,
261  unsigned uPrefixCountA, unsigned uPrefixCountB)
262        {
263        if (0 == g_DP)
264                {
265                Log("***DPDIFF*** DP%c=NULL\n", c);
266                return true;
267                }
268
269        SCORE *DPM_ = g_DP;
270        for (unsigned i = 0; i < uPrefixCountA; ++i)
271                for (unsigned j = 0; j < uPrefixCountB; ++j)
272                        if (!LocalEq(DPM(i, j), DPD(i, j)))
273                                {
274                                Log("***DPDIFF*** DP%c(%d, %d) Simple = %.2g, Small = %.2g\n",
275                                  c, i, j, DPM(i, j), DPD(i, j));
276                                return false;
277                                }
278        return true;
279        }
280
281static bool CompareTB(char **TB, char *TBM_, char *TBD_, char *TBE_, char *TBI_, char *TBJ_,
282  unsigned uPrefixCountA, unsigned uPrefixCountB)
283        {
284        if (!g_bKeepSimpleDP)
285                return true;
286        SCORE *DPM_ = g_DPM;
287        bool Eq = true;
288        for (unsigned i = 0; i < uPrefixCountA; ++i)
289                for (unsigned j = 0; j < uPrefixCountB; ++j)
290                        {
291                        char c1 = TBM(i, j);
292                        char c2 = Get_M_Char(TB[i][j]);
293                        if (c1 != '?' && c1 != c2 && DPM(i, j) > -100000)
294                                {
295                                Log("TBM(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
296                                Eq = false;
297                                goto D;
298                                }
299                        }
300
301D:
302        SCORE *DPD_ = g_DPD;
303        for (unsigned i = 0; i < uPrefixCountA; ++i)
304                for (unsigned j = 0; j < uPrefixCountB; ++j)
305                        {
306                        char c1 = TBD(i, j);
307                        char c2 = Get_D_Char(TB[i][j]);
308                        if (c1 != '?' && c1 != c2 && DPD(i, j) > -100000)
309                                {
310                                Log("TBD(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
311                                Eq = false;
312                                goto E;
313                                }
314                        }
315E:
316        SCORE *DPE_ = g_DPE;
317        if (0 == TBE_)
318                goto I;
319        for (unsigned i = 0; i < uPrefixCountA; ++i)
320                for (unsigned j = 0; j < uPrefixCountB; ++j)
321                        {
322                        char c1 = TBE(i, j);
323                        char c2 = Get_E_Char(TB[i][j]);
324                        if (c1 != '?' && c1 != c2 && DPE(i, j) > -100000)
325                                {
326                                Log("TBE(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
327                                Eq = false;
328                                goto I;
329                                }
330                        }
331I:
332        SCORE *DPI_ = g_DPI;
333        for (unsigned i = 0; i < uPrefixCountA; ++i)
334                for (unsigned j = 0; j < uPrefixCountB; ++j)
335                        {
336                        char c1 = TBI(i, j);
337                        char c2 = Get_I_Char(TB[i][j]);
338                        if (c1 != '?' && c1 != c2 && DPI(i, j) > -100000)
339                                {
340                                Log("TBI(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
341                                Eq = false;
342                                goto J;
343                                }
344                        }
345J:
346        SCORE *DPJ_ = g_DPJ;
347        if (0 == DPJ_)
348                goto Done;
349        for (unsigned i = 0; i < uPrefixCountA; ++i)
350                for (unsigned j = 0; j < uPrefixCountB; ++j)
351                        {
352                        char c1 = TBJ(i, j);
353                        char c2 = Get_J_Char(TB[i][j]);
354                        if (c1 != '?' && c1 != c2 && DPJ(i, j) > -100000)
355                                {
356                                Log("TBJ(%d, %d) Simple = %c, NW = %c\n", i, j, c1, c2);
357                                Eq = false;
358                                goto Done;
359                                }
360                        }
361Done:
362        if (Eq)
363                Log("TB success\n");
364        return Eq;
365        }
366
367static const char *LocalScoreToStr(SCORE s)
368        {
369        static char str[16];
370        if (s < -100000)
371                return "     *";
372        sprintf(str, "%6.1f", s);
373        return str;
374        }
375
376static void LogDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,
377  unsigned uPrefixCountA, unsigned uPrefixCountB)
378        {
379        Log("        ");
380        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
381                {
382                char c = ' ';
383                if (uPrefixLengthB > 0)
384                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
385                Log(" %4u:%c", uPrefixLengthB, c);
386                }
387        Log("\n");
388        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
389                {
390                char c = ' ';
391                if (uPrefixLengthA > 0)
392                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
393                Log("%4u:%c  ", uPrefixLengthA, c);
394                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
395                        Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));
396                Log("\n");
397                }
398        }
399
400static void LogBitTB(char **TB, const ProfPos *PA, const ProfPos *PB,
401  unsigned uPrefixCountA, unsigned uPrefixCountB)
402        {
403        Log("        ");
404        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
405                {
406                char c = ' ';
407                if (uPrefixLengthB > 0)
408                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
409                Log(" %4u:%c", uPrefixLengthB, c);
410                }
411        Log("\n");
412        Log("Bit TBM:\n");
413        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
414                {
415                char c = ' ';
416                if (uPrefixLengthA > 0)
417                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
418                Log("%4u:%c  ", uPrefixLengthA, c);
419                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
420                        {
421                        char c = Get_M_Char(TB[uPrefixLengthA][uPrefixLengthB]);
422                        Log(" %6c", c);
423                        }
424                Log("\n");
425                }
426
427        Log("\n");
428        Log("Bit TBD:\n");
429        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
430                {
431                char c = ' ';
432                if (uPrefixLengthA > 0)
433                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
434                Log("%4u:%c  ", uPrefixLengthA, c);
435                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
436                        {
437                        char c = Get_D_Char(TB[uPrefixLengthA][uPrefixLengthB]);
438                        Log(" %6c", c);
439                        }
440                Log("\n");
441                }
442
443        Log("\n");
444        Log("Bit TBE:\n");
445        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
446                {
447                char c = ' ';
448                if (uPrefixLengthA > 0)
449                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
450                Log("%4u:%c  ", uPrefixLengthA, c);
451                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
452                        {
453                        char c = Get_E_Char(TB[uPrefixLengthA][uPrefixLengthB]);
454                        Log(" %6c", c);
455                        }
456                Log("\n");
457                }
458
459        Log("\n");
460        Log("Bit TBI:\n");
461        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
462                {
463                char c = ' ';
464                if (uPrefixLengthA > 0)
465                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
466                Log("%4u:%c  ", uPrefixLengthA, c);
467                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
468                        {
469                        char c = Get_I_Char(TB[uPrefixLengthA][uPrefixLengthB]);
470                        Log(" %6c", c);
471                        }
472                Log("\n");
473                }
474
475        Log("\n");
476        Log("Bit TBJ:\n");
477        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
478                {
479                char c = ' ';
480                if (uPrefixLengthA > 0)
481                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
482                Log("%4u:%c  ", uPrefixLengthA, c);
483                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
484                        {
485                        char c = Get_J_Char(TB[uPrefixLengthA][uPrefixLengthB]);
486                        Log(" %6c", c);
487                        }
488                Log("\n");
489                }
490        }
491
492static void ListTB(char *TBM_, const ProfPos *PA, const ProfPos *PB,
493  unsigned uPrefixCountA, unsigned uPrefixCountB)
494        {
495        Log("        ");
496        for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
497                {
498                char c = ' ';
499                if (uPrefixLengthB > 0)
500                        c = ConsensusChar(PB[uPrefixLengthB - 1]);
501                Log(" %4u:%c", uPrefixLengthB, c);
502                }
503        Log("\n");
504        for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)
505                {
506                char c = ' ';
507                if (uPrefixLengthA > 0)
508                        c = ConsensusChar(PA[uPrefixLengthA - 1]);
509                Log("%4u:%c  ", uPrefixLengthA, c);
510                for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)
511                        {
512                        char c = TBM(uPrefixLengthA, uPrefixLengthB);
513                        Log(" %6c", c);
514                        }
515                Log("\n");
516                }
517        }
518
519static const char *BitsToStr(char Bits)
520        {
521        static char Str[32];
522
523        sprintf(Str, "%cM %cD %cE %cI %cJ",
524          Get_M_Char(Bits),
525          Get_D_Char(Bits),
526          Get_E_Char(Bits),
527          Get_I_Char(Bits),
528          Get_J_Char(Bits));
529        }
530#endif  // TRACE
531
532static inline void SetBitTBM(char **TB, unsigned i, unsigned j, char c)
533        {
534        char Bit;
535        switch (c)
536                {
537        case 'M':
538                Bit = BIT_MM;
539                break;
540        case 'D':
541                Bit = BIT_DM;
542                break;
543#if     DOUBLE_AFFINE
544        case 'E':
545                Bit = BIT_EM;
546                break;
547        case 'I':
548                Bit = BIT_IM;
549                break;
550        case 'J':
551                Bit = BIT_JM;
552                break;
553#endif
554        default:
555                Quit("Huh?!");
556                }
557        TB[i][j] &= ~BIT_xM;
558        TB[i][j] |= Bit;
559        }
560
561static inline void SetBitTBD(char **TB, unsigned i, unsigned j, char c)
562        {
563        char Bit;
564        switch (c)
565                {
566        case 'M':
567                Bit = BIT_MD;
568                break;
569        case 'D':
570                Bit = BIT_DD;
571                break;
572        default:
573                Quit("Huh?!");
574                }
575        TB[i][j] &= ~BIT_xD;
576        TB[i][j] |= Bit;
577        }
578
579static inline void SetBitTBI(char **TB, unsigned i, unsigned j, char c)
580        {
581        char Bit;
582        switch (c)
583                {
584        case 'M':
585                Bit = BIT_MI;
586                break;
587        case 'I':
588                Bit = BIT_II;
589                break;
590        default:
591                Quit("Huh?!");
592                }
593        TB[i][j] &= ~BIT_xI;
594        TB[i][j] |= Bit;
595        }
596
597#if     DOUBLE_AFFINE
598static inline void SetBitTBE(char **TB, unsigned i, unsigned j, char c)
599        {
600        char Bit;
601        switch (c)
602                {
603        case 'M':
604                Bit = BIT_ME;
605                break;
606        case 'E':
607                Bit = BIT_EE;
608                break;
609        default:
610                Quit("Huh?!");
611                }
612        TB[i][j] &= ~BIT_xE;
613        TB[i][j] |= Bit;
614        }
615
616static inline void SetBitTBJ(char **TB, unsigned i, unsigned j, char c)
617        {
618        char Bit;
619        switch (c)
620                {
621        case 'M':
622                Bit = BIT_MJ;
623                break;
624        case 'J':
625                Bit = BIT_JJ;
626                break;
627        default:
628                Quit("Huh?!");
629                }
630        TB[i][j] &= ~BIT_xJ;
631        TB[i][j] |= Bit;
632        }
633#endif
634
635#if     TRACE
636#define LogMatrices()                                                                                   \
637        {                                                                                                                       \
638        Log("Bit DPM:\n");                                                                                      \
639        LogDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
640        Log("Bit DPD:\n");                                                                                      \
641        LogDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
642        Log("Bit DPE:\n");                                                                                      \
643        LogDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
644        Log("Bit DPI:\n");                                                                                      \
645        LogDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
646        Log("Bit DPJ:\n");                                                                                      \
647        LogDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);                      \
648        Log("Bit TB:\n");                                                                                       \
649        LogBitTB(TB, PA, PB, uPrefixCountA, uPrefixCountB);                     \
650        bool Same;                                                                                                      \
651        Same = DPEq('M', g_DPM, DPM_, uPrefixCountA, uPrefixCountB);\
652        if (Same)                                                                                                       \
653                Log("DPM success\n");                                                                   \
654        Same = DPEq('D', g_DPD, DPD_, uPrefixCountA, uPrefixCountB);\
655        if (Same)                                                                                                       \
656                Log("DPD success\n");                                                                   \
657        Same = DPEq('E', g_DPE, DPE_, uPrefixCountA, uPrefixCountB);\
658        if (Same)                                                                                                       \
659                Log("DPE success\n");                                                                   \
660        Same = DPEq('I', g_DPI, DPI_, uPrefixCountA, uPrefixCountB);\
661        if (Same)                                                                                                       \
662                Log("DPI success\n");                                                                   \
663        Same = DPEq('J', g_DPJ, DPJ_, uPrefixCountA, uPrefixCountB);\
664        if (Same)                                                                                                       \
665                Log("DPJ success\n");                                                                   \
666        CompareTB(TB, g_TBM, g_TBD, g_TBE, g_TBI, g_TBJ, uPrefixCountA, uPrefixCountB);\
667        }
668#else
669#define LogMatrices()   /* empty */
670#endif
671
672SCORE NWDASmall(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
673  unsigned uLengthB, PWPath &Path)
674        {
675        assert(uLengthB > 0 && uLengthA > 0);
676
677        ProfPos *pa0 = (ProfPos *) PA;
678        ProfPos *pb0 = (ProfPos *) PB;
679        ProfPos *paa = (ProfPos *) (PA + uLengthA - 1);
680        ProfPos *pbb = (ProfPos *) (PB + uLengthB - 1);
681
682        pa0->m_scoreGapOpen *= -1;
683        pb0->m_scoreGapOpen *= -1;
684
685        paa->m_scoreGapClose *= -1;
686        pbb->m_scoreGapClose *= -1;
687
688        pa0->m_scoreGapOpen2 *= -1;
689        pb0->m_scoreGapOpen2 *= -1;
690        paa->m_scoreGapClose2 *= -1;
691        pbb->m_scoreGapClose2 *= -1;
692
693        const unsigned uPrefixCountA = uLengthA + 1;
694        const unsigned uPrefixCountB = uLengthB + 1;
695        const SCORE e = g_scoreGapExtend;
696
697        const SCORE e2 = g_scoreGapExtend2;
698        const SCORE min_e = MIN(g_scoreGapExtend, g_scoreGapExtend2);
699
700        ALLOC_TRACE()
701
702        SCORE *MCurr = new SCORE[uPrefixCountB];
703        SCORE *MNext = new SCORE[uPrefixCountB];
704        SCORE *MPrev = new SCORE[uPrefixCountB];
705        SCORE *DRow = new SCORE[uPrefixCountB];
706        SCORE *ERow = new SCORE[uPrefixCountB];
707
708        char **TB = new char *[uPrefixCountA];
709        for (unsigned i = 0; i < uPrefixCountA; ++i)
710                {
711                TB[i] = new char [uPrefixCountB];
712                memset(TB[i], 0, uPrefixCountB);
713                }
714
715        SCORE Iij = MINUS_INFINITY;
716        SetDPI(0, 0, Iij);
717
718        SCORE Jij = MINUS_INFINITY;
719        SetDPJ(0, 0, Jij);
720
721        Iij = PB[0].m_scoreGapOpen;
722        SetDPI(0, 1, Iij);
723
724        Jij = PB[0].m_scoreGapOpen2;
725        SetDPJ(0, 1, Jij);
726
727        for (unsigned j = 2; j <= uLengthB; ++j)
728                {
729                Iij += e;
730                Jij += e2;
731
732                SetDPI(0, j, Iij);
733                SetDPJ(0, j, Jij);
734
735                SetTBI(0, j, 'I');
736                SetTBJ(0, j, 'J');
737                }
738
739        for (unsigned j = 0; j <= uLengthB; ++j)
740                {
741                DRow[j] = MINUS_INFINITY;
742                ERow[j] = MINUS_INFINITY;
743
744                SetDPD(0, j, DRow[j]);
745                SetDPE(0, j, ERow[j]);
746
747                SetTBD(0, j, 'D');
748                SetTBE(0, j, 'E');
749                }
750
751        MPrev[0] = 0;
752        SetDPM(0, 0, MPrev[0]);
753        for (unsigned j = 1; j <= uLengthB; ++j)
754                {
755                MPrev[j] = MINUS_INFINITY;
756                SetDPM(0, j, MPrev[j]);
757                }
758
759        MCurr[0] = MINUS_INFINITY;
760        SetDPM(1, 0, MCurr[0]);
761
762        MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
763        SetDPM(1, 1, MCurr[1]);
764        SetBitTBM(TB, 1, 1, 'M');
765        SetTBM(1, 1, 'M');
766
767        for (unsigned j = 2; j <= uLengthB; ++j)
768                {
769                SCORE M = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen +
770                  (j - 2)*e + PB[j-2].m_scoreGapClose;
771                SCORE M2 = ScoreProfPos2(PA[0], PB[j-1]) + PB[0].m_scoreGapOpen2 +
772                  (j - 2)*e2 + PB[j-2].m_scoreGapClose2;
773               
774                if (M >= M2)
775                        {
776                        MCurr[j] = M;
777                        SetBitTBM(TB, 1, j, 'I');
778                        SetTBM(1, j, 'I');
779                        }
780                else
781                        {
782                        MCurr[j] = M2;
783                        SetBitTBM(TB, 1, j, 'J');
784                        SetTBM(1, j, 'J');
785                        }
786                SetDPM(1, j, MCurr[j]);
787                }
788
789// Main DP loop
790        for (unsigned i = 1; i < uLengthA; ++i)
791                {
792                Iij = MINUS_INFINITY;
793                Jij = MINUS_INFINITY;
794                SetDPI(i, 0, Iij);
795                SetDPJ(i, 0, Jij);
796
797                DRow[0] = PA[0].m_scoreGapOpen + (i - 1)*e;
798                ERow[0] = PA[0].m_scoreGapOpen2 + (i - 1)*e2;
799                SetDPD(i, 0, DRow[0]);
800                SetDPE(i, 0, ERow[0]);
801
802                MCurr[0] = MINUS_INFINITY; 
803                if (i == 1)
804                        {
805                        MCurr[1] = ScoreProfPos2(PA[0], PB[0]);
806                        SetBitTBM(TB, i, 1, 'M');
807                        SetTBM(i, 1, 'M');
808                        }
809                else
810                        {
811                        SCORE M = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen +
812                          (i - 2)*e + PA[i-2].m_scoreGapClose;
813                        SCORE M2 = ScoreProfPos2(PA[i-1], PB[0]) + PA[0].m_scoreGapOpen2 +
814                          (i - 2)*e2 + PA[i-2].m_scoreGapClose2;
815                        if (M >= M2)
816                                {
817                                MCurr[1] = M;
818                                SetBitTBM(TB, i, 1, 'D');
819                                SetTBM(i, 1, 'D');
820                                }
821                        else
822                                {
823                                MCurr[1] = M2;
824                                SetBitTBM(TB, i, 1, 'E');
825                                SetTBM(i, 1, 'E');
826                                }
827                        }
828                SetDPM(i, 0, MCurr[0]);
829                SetDPM(i, 1, MCurr[1]);
830
831                for (unsigned j = 1; j < uLengthB; ++j)
832                        MNext[j+1] = ScoreProfPos2(PA[i], PB[j]);
833
834                for (unsigned j = 1; j < uLengthB; ++j)
835                        {
836                        RECURSE_D(i, j)
837                        RECURSE_E(i, j)
838                        RECURSE_I(i, j)
839                        RECURSE_J(i, j)
840                        RECURSE_M(i, j)
841                        }
842        // Special case for j=uLengthB
843                RECURSE_D_BTerm(i)
844                RECURSE_E_BTerm(i)
845                RECURSE_I_BTerm(i)
846                RECURSE_J_BTerm(i)
847
848        // Prev := Curr, Curr := Next, Next := Prev
849                Rotate(MPrev, MCurr, MNext);
850                }
851
852// Special case for i=uLengthA
853        MCurr[0] = MINUS_INFINITY;
854        SCORE M = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
855          PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
856        SCORE M2 = ScoreProfPos2(PA[uLengthA-1], PB[0]) + (uLengthA - 2)*e +
857          PA[0].m_scoreGapOpen + PA[uLengthA-2].m_scoreGapClose;
858        if (M >= M2)
859                {
860                MCurr[1] = M;
861                SetBitTBM(TB, uLengthA, 1, 'D');
862                SetTBM(uLengthA, 1, 'D');
863                }
864        else
865                {
866                MCurr[1] = M2;
867                SetBitTBM(TB, uLengthA, 1, 'E');
868                SetTBM(uLengthA, 1, 'D');
869                }
870        SetDPM(uLengthA, 0, MCurr[0]);
871        SetDPM(uLengthA, 1, MCurr[1]);
872
873        DRow[0] = MINUS_INFINITY;
874        ERow[0] = MINUS_INFINITY;
875
876        SetDPD(uLengthA, 0, DRow[0]);
877        SetDPE(uLengthA, 0, ERow[0]);
878
879        for (unsigned j = 1; j <= uLengthB; ++j)
880                {
881                RECURSE_D_ATerm(j);
882                RECURSE_E_ATerm(j);
883                }
884
885        Iij = MINUS_INFINITY;
886        Jij = MINUS_INFINITY;
887
888        for (unsigned j = 1; j <= uLengthB; ++j)
889                {
890                RECURSE_I_ATerm(j)
891                RECURSE_J_ATerm(j)
892                }
893
894        LogMatrices();
895
896        SCORE MAB = MCurr[uLengthB];
897        SCORE DAB = DRow[uLengthB] + PA[uLengthA-1].m_scoreGapClose;
898        SCORE EAB = ERow[uLengthB] + PA[uLengthA-1].m_scoreGapClose2;
899        SCORE IAB = Iij + PB[uLengthB-1].m_scoreGapClose;
900        SCORE JAB = Jij + PB[uLengthB-1].m_scoreGapClose2;
901
902        SCORE Score = MAB;
903        char cEdgeType = 'M';
904        if (DAB > Score)
905                {
906                Score = DAB;
907                cEdgeType = 'D';
908                }
909        if (EAB > Score)
910                {
911                Score = EAB;
912                cEdgeType = 'E';
913                }
914        if (IAB > Score)
915                {
916                Score = IAB;
917                cEdgeType = 'I';
918                }
919        if (JAB > Score)
920                {
921                Score = JAB;
922                cEdgeType = 'J';
923                }
924
925#if TRACE
926        Log("    Small: MAB=%.4g DAB=%.4g EAB=%.4g IAB=%.4g JAB=%.4g best=%c\n",
927          MAB, DAB, EAB, IAB, JAB, cEdgeType);
928#endif
929
930        BitTraceBack(TB, uLengthA, uLengthB, cEdgeType, Path);
931
932#if     DBEUG
933        Path.Validate();
934#endif
935
936        delete[] MCurr;
937        delete[] MNext;
938        delete[] MPrev;
939        delete[] DRow;
940        delete[] ERow;
941        for (unsigned i = 0; i < uPrefixCountA; ++i)
942                delete[] TB[i];
943        delete[] TB;
944
945        return 0;
946        }
947#endif  // DOUBLE_AFFINE
Note: See TracBrowser for help on using the repository browser.