website/archive/binaries/mac/src/muscle/nwsmall.cpp

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