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