binaries/src/muscle/nwdasimple2.cpp

   1 #include "muscle.h"\r
   2 #include "pwpath.h"\r
   3 #include "profile.h"\r
   4 \r
   5 #if     DOUBLE_AFFINE\r
   6 \r
   7 #define TRACE   0\r
   8 \r
   9 extern bool g_bKeepSimpleDP;\r
  10 extern SCORE *g_DPM;\r
  11 extern SCORE *g_DPD;\r
  12 extern SCORE *g_DPE;\r
  13 extern SCORE *g_DPI;\r
  14 extern SCORE *g_DPJ;\r
  15 extern char *g_TBM;\r
  16 extern char *g_TBD;\r
  17 extern char *g_TBE;\r
  18 extern char *g_TBI;\r
  19 extern char *g_TBJ;\r
  20 \r
  21 static char XlatEdgeType(char c)\r
  22         {\r
  23         if ('E' == c)\r
  24                 return 'D';\r
  25         if ('J' == c)\r
  26                 return 'I';\r
  27         return c;\r
  28         }\r
  29 \r
  30 static const char *LocalScoreToStr(SCORE s)\r
  31         {\r
  32         static char str[16];\r
  33         if (s < -100000)\r
  34                 return "     *";\r
  35         sprintf(str, "%6.1f", s);\r
  36         return str;\r
  37         }\r
  38 \r
  39 static void ListDP(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
  40   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  41         {\r
  42         Log("        ");\r
  43         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  44                 {\r
  45                 char c = ' ';\r
  46                 if (uPrefixLengthB > 0)\r
  47                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  48                 Log(" %4u:%c", uPrefixLengthB, c);\r
  49                 }\r
  50         Log("\n");\r
  51         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  52                 {\r
  53                 char c = ' ';\r
  54                 if (uPrefixLengthA > 0)\r
  55                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
  56                 Log("%4u:%c  ", uPrefixLengthA, c);\r
  57                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  58                         Log(" %s", LocalScoreToStr(DPM(uPrefixLengthA, uPrefixLengthB)));\r
  59                 Log("\n");\r
  60                 }\r
  61         }\r
  62 \r
  63 static void ListTB(const char *TBM_, const ProfPos *PA, const ProfPos *PB,\r
  64   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  65         {\r
  66         Log("        ");\r
  67         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  68                 {\r
  69                 char c = ' ';\r
  70                 if (uPrefixLengthB > 0)\r
  71                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  72                 Log(" %4u:%c", uPrefixLengthB, c);\r
  73                 }\r
  74         Log("\n");\r
  75         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
  76                 {\r
  77                 char c = ' ';\r
  78                 if (uPrefixLengthA > 0)\r
  79                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
  80                 Log("%4u:%c  ", uPrefixLengthA, c);\r
  81                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  82                         Log(" %6c", TBM(uPrefixLengthA, uPrefixLengthB));\r
  83                 Log("\n");\r
  84                 }\r
  85         }\r
  86 \r
  87 static void ListDPM(const SCORE *DPM_, const ProfPos *PA, const ProfPos *PB,\r
  88   unsigned uPrefixCountA, unsigned uPrefixCountB)\r
  89         {\r
  90         Log("        ");\r
  91         for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
  92                 {\r
  93                 char c = ' ';\r
  94                 if (uPrefixLengthB > 0)\r
  95                         c = ConsensusChar(PB[uPrefixLengthB - 1]);\r
  96                 Log(" %4u:%c", uPrefixLengthB, c);\r
  97                 }\r
  98         Log("\n");\r
  99         for (unsigned uPrefixLengthA = 0; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 100                 {\r
 101                 char c = ' ';\r
 102                 if (uPrefixLengthA > 0)\r
 103                         c = ConsensusChar(PA[uPrefixLengthA - 1]);\r
 104                 Log("%4u:%c  ", uPrefixLengthA, c);\r
 105                 for (unsigned uPrefixLengthB = 0; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 106                         {\r
 107                         SCORE x = (uPrefixLengthA + uPrefixLengthB)*g_scoreGapExtend;\r
 108                         SCORE s = DPM(uPrefixLengthA, uPrefixLengthB) - x;\r
 109                         Log(" %s", LocalScoreToStr(s));\r
 110                         }\r
 111                 Log("\n");\r
 112                 }\r
 113         }\r
 114 \r
 115 extern SCORE ScoreProfPos2(const ProfPos &PP, const ProfPos &PPB);\r
 116 \r
 117 SCORE NWDASimple2(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
 118   unsigned uLengthB, PWPath &Path)\r
 119         {\r
 120         assert(uLengthB > 0 && uLengthA > 0);\r
 121 \r
 122         const unsigned uPrefixCountA = uLengthA + 1;\r
 123         const unsigned uPrefixCountB = uLengthB + 1;\r
 124 \r
 125 // Allocate DP matrices\r
 126         const size_t LM = uPrefixCountA*uPrefixCountB;\r
 127         SCORE *DPM_ = new SCORE[LM];\r
 128         SCORE *DPD_ = new SCORE[LM];\r
 129         SCORE *DPE_ = new SCORE[LM];\r
 130         SCORE *DPI_ = new SCORE[LM];\r
 131         SCORE *DPJ_ = new SCORE[LM];\r
 132         SCORE *DPL_ = new SCORE[LM];\r
 133 \r
 134         char *TBM_ = new char[LM];\r
 135         char *TBD_ = new char[LM];\r
 136         char *TBE_ = new char[LM];\r
 137         char *TBI_ = new char[LM];\r
 138         char *TBJ_ = new char[LM];\r
 139 \r
 140         memset(DPM_, 0, LM*sizeof(SCORE));\r
 141         memset(DPD_, 0, LM*sizeof(SCORE));\r
 142         memset(DPE_, 0, LM*sizeof(SCORE));\r
 143         memset(DPI_, 0, LM*sizeof(SCORE));\r
 144         memset(DPJ_, 0, LM*sizeof(SCORE));\r
 145 \r
 146 //      memset(DPL_, 0, LM*sizeof(SCORE));\r
 147 \r
 148         memset(TBM_, '?', LM);\r
 149         memset(TBD_, '?', LM);\r
 150         memset(TBE_, '?', LM);\r
 151         memset(TBI_, '?', LM);\r
 152         memset(TBJ_, '?', LM);\r
 153 \r
 154         DPM(0, 0) = 0;\r
 155         DPD(0, 0) = MINUS_INFINITY;\r
 156         DPE(0, 0) = MINUS_INFINITY;\r
 157         DPI(0, 0) = MINUS_INFINITY;\r
 158         DPJ(0, 0) = MINUS_INFINITY;\r
 159 \r
 160         DPM(1, 0) = MINUS_INFINITY;\r
 161         DPD(1, 0) = PA[0].m_scoreGapOpen;\r
 162         DPE(1, 0) = PA[0].m_scoreGapOpen2;\r
 163         DPI(1, 0) = MINUS_INFINITY;\r
 164         DPJ(1, 0) = MINUS_INFINITY;\r
 165 \r
 166         DPM(0, 1) = MINUS_INFINITY;\r
 167         DPD(0, 1) = MINUS_INFINITY;\r
 168         DPE(0, 1) = MINUS_INFINITY;\r
 169         DPI(0, 1) = PB[0].m_scoreGapOpen;\r
 170         DPJ(0, 1) = PB[0].m_scoreGapOpen2;\r
 171 \r
 172 // Empty prefix of B is special case\r
 173         for (unsigned uPrefixLengthA = 2; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 174                 {\r
 175         // M=LetterA+LetterB, impossible with empty prefix\r
 176                 DPM(uPrefixLengthA, 0) = MINUS_INFINITY;\r
 177 \r
 178         // D=LetterA+GapB\r
 179                 DPD(uPrefixLengthA, 0) = DPD(uPrefixLengthA - 1, 0) + g_scoreGapExtend;\r
 180                 TBD(uPrefixLengthA, 0) = 'D';\r
 181 \r
 182                 DPE(uPrefixLengthA, 0) = DPE(uPrefixLengthA - 1, 0) + g_scoreGapExtend2;\r
 183                 TBE(uPrefixLengthA, 0) = 'E';\r
 184 \r
 185         // I=GapA+LetterB, impossible with empty prefix\r
 186                 DPI(uPrefixLengthA, 0) = MINUS_INFINITY;\r
 187                 DPJ(uPrefixLengthA, 0) = MINUS_INFINITY;\r
 188                 }\r
 189 \r
 190 // Empty prefix of A is special case\r
 191         for (unsigned uPrefixLengthB = 2; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 192                 {\r
 193         // M=LetterA+LetterB, impossible with empty prefix\r
 194                 DPM(0, uPrefixLengthB) = MINUS_INFINITY;\r
 195 \r
 196         // D=LetterA+GapB, impossible with empty prefix\r
 197                 DPD(0, uPrefixLengthB) = MINUS_INFINITY;\r
 198                 DPE(0, uPrefixLengthB) = MINUS_INFINITY;\r
 199 \r
 200         // I=GapA+LetterB\r
 201                 DPI(0, uPrefixLengthB) = DPI(0, uPrefixLengthB - 1) + g_scoreGapExtend;\r
 202                 TBI(0, uPrefixLengthB) = 'I';\r
 203 \r
 204                 DPJ(0, uPrefixLengthB) = DPJ(0, uPrefixLengthB - 1) + g_scoreGapExtend2;\r
 205                 TBJ(0, uPrefixLengthB) = 'J';\r
 206                 }\r
 207 \r
 208 // ============\r
 209 // Main DP loop\r
 210 // ============\r
 211         for (unsigned uPrefixLengthB = 1; uPrefixLengthB < uPrefixCountB; ++uPrefixLengthB)\r
 212                 {\r
 213                 const ProfPos &PPB = PB[uPrefixLengthB - 1];\r
 214                 SCORE scoreGapCloseB;\r
 215                 if (uPrefixLengthB == 1)\r
 216                         scoreGapCloseB = MINUS_INFINITY;\r
 217                 else\r
 218                         scoreGapCloseB = PB[uPrefixLengthB-2].m_scoreGapClose;\r
 219 \r
 220                 SCORE scoreGapClose2B;\r
 221                 if (uPrefixLengthB == 1)\r
 222                         scoreGapClose2B = MINUS_INFINITY;\r
 223                 else\r
 224                         scoreGapClose2B = PB[uPrefixLengthB-2].m_scoreGapClose2;\r
 225 \r
 226                 for (unsigned uPrefixLengthA = 1; uPrefixLengthA < uPrefixCountA; ++uPrefixLengthA)\r
 227                         {\r
 228                         const ProfPos &PPA = PA[uPrefixLengthA - 1];\r
 229 \r
 230                         {\r
 231                 // Match M=LetterA+LetterB\r
 232                         SCORE scoreLL = ScoreProfPos2(PPA, PPB);\r
 233                         DPL(uPrefixLengthA, uPrefixLengthB) = scoreLL;\r
 234 \r
 235                         SCORE scoreGapCloseA;\r
 236                         if (uPrefixLengthA == 1)\r
 237                                 scoreGapCloseA = MINUS_INFINITY;\r
 238                         else\r
 239                                 scoreGapCloseA = PA[uPrefixLengthA-2].m_scoreGapClose;\r
 240 \r
 241                         SCORE scoreGapClose2A;\r
 242                         if (uPrefixLengthA == 1)\r
 243                                 scoreGapClose2A = MINUS_INFINITY;\r
 244                         else\r
 245                                 scoreGapClose2A = PA[uPrefixLengthA-2].m_scoreGapClose2;\r
 246 \r
 247                         SCORE scoreMM = DPM(uPrefixLengthA-1, uPrefixLengthB-1);\r
 248                         SCORE scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseA;\r
 249                         SCORE scoreEM = DPE(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2A;\r
 250                         SCORE scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapCloseB;\r
 251                         SCORE scoreJM = DPJ(uPrefixLengthA-1, uPrefixLengthB-1) + scoreGapClose2B;\r
 252                         SCORE scoreBest;\r
 253                         if (scoreMM >= scoreDM && scoreMM >= scoreIM && scoreMM >= scoreEM && scoreMM >= scoreJM)\r
 254                                 {\r
 255                                 scoreBest = scoreMM;\r
 256                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 257                                 }\r
 258                         else if (scoreDM >= scoreMM && scoreDM >= scoreIM && scoreDM >= scoreEM && scoreDM >= scoreJM)\r
 259                                 {\r
 260                                 scoreBest = scoreDM;\r
 261                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'D';\r
 262                                 }\r
 263                         else if (scoreEM >= scoreMM && scoreEM >= scoreIM && scoreEM >= scoreDM && scoreEM >= scoreJM)\r
 264                                 {\r
 265                                 scoreBest = scoreEM;\r
 266                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'E';\r
 267                                 }\r
 268                         else if (scoreIM >= scoreMM && scoreIM >= scoreDM && scoreIM >= scoreEM && scoreIM >= scoreJM)\r
 269                                 {\r
 270                                 scoreBest = scoreIM;\r
 271                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'I';\r
 272                                 }\r
 273                         else if (scoreJM >= scoreMM && scoreJM >= scoreDM && scoreJM >= scoreEM && scoreJM >= scoreIM)\r
 274                                 {\r
 275                                 scoreBest = scoreJM;\r
 276                                 TBM(uPrefixLengthA, uPrefixLengthB) = 'J';\r
 277                                 }\r
 278                         else\r
 279                                 Quit("Max failed (M)");\r
 280 \r
 281                         DPM(uPrefixLengthA, uPrefixLengthB) = scoreBest + scoreLL;\r
 282                         }\r
 283 \r
 284                         {\r
 285                 // Delete D=LetterA+GapB\r
 286                         SCORE scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
 287                           PA[uPrefixLengthA-1].m_scoreGapOpen;\r
 288                         SCORE scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB) +\r
 289                           g_scoreGapExtend;\r
 290 \r
 291                         SCORE scoreBest;\r
 292                         if (scoreMD >= scoreDD)\r
 293                                 {\r
 294                                 scoreBest = scoreMD;\r
 295                                 TBD(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 296                                 }\r
 297                         else\r
 298                                 {\r
 299                                 assert(scoreDD >= scoreMD);\r
 300                                 scoreBest = scoreDD;\r
 301                                 TBD(uPrefixLengthA, uPrefixLengthB) = 'D';\r
 302                                 }\r
 303                         DPD(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 304                         }\r
 305 \r
 306                         {\r
 307                 // Delete E=LetterA+GapB\r
 308                         SCORE scoreME = DPM(uPrefixLengthA-1, uPrefixLengthB) +\r
 309                           PA[uPrefixLengthA-1].m_scoreGapOpen2;\r
 310                         SCORE scoreEE = DPE(uPrefixLengthA-1, uPrefixLengthB) +\r
 311                           g_scoreGapExtend2;\r
 312 \r
 313                         SCORE scoreBest;\r
 314                         if (scoreME >= scoreEE)\r
 315                                 {\r
 316                                 scoreBest = scoreME;\r
 317                                 TBE(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 318                                 }\r
 319                         else\r
 320                                 {\r
 321                                 assert(scoreEE >= scoreME);\r
 322                                 scoreBest = scoreEE;\r
 323                                 TBE(uPrefixLengthA, uPrefixLengthB) = 'E';\r
 324                                 }\r
 325                         DPE(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 326                         }\r
 327 \r
 328                 // Insert I=GapA+LetterB\r
 329                         {\r
 330                         SCORE scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
 331                           PB[uPrefixLengthB-1].m_scoreGapOpen;\r
 332                         SCORE scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1) +\r
 333                           g_scoreGapExtend;\r
 334 \r
 335                         SCORE scoreBest;\r
 336                         if (scoreMI >= scoreII)\r
 337                                 {\r
 338                                 scoreBest = scoreMI;\r
 339                                 TBI(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 340                                 }\r
 341                         else \r
 342                                 {\r
 343                                 assert(scoreII > scoreMI);\r
 344                                 scoreBest = scoreII;\r
 345                                 TBI(uPrefixLengthA, uPrefixLengthB) = 'I';\r
 346                                 }\r
 347                         DPI(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 348                         }\r
 349 \r
 350                 // Insert J=GapA+LetterB\r
 351                         {\r
 352                         SCORE scoreMJ = DPM(uPrefixLengthA, uPrefixLengthB-1) +\r
 353                           PB[uPrefixLengthB-1].m_scoreGapOpen2;\r
 354                         SCORE scoreJJ = DPJ(uPrefixLengthA, uPrefixLengthB-1) +\r
 355                           g_scoreGapExtend2;\r
 356 \r
 357                         SCORE scoreBest;\r
 358                         if (scoreMJ > scoreJJ)\r
 359                                 {\r
 360                                 scoreBest = scoreMJ;\r
 361                                 TBJ(uPrefixLengthA, uPrefixLengthB) = 'M';\r
 362                                 }\r
 363                         else \r
 364                                 {\r
 365                                 assert(scoreJJ >= scoreMJ);\r
 366                                 scoreBest = scoreJJ;\r
 367                                 TBJ(uPrefixLengthA, uPrefixLengthB) = 'J';\r
 368                                 }\r
 369                         DPJ(uPrefixLengthA, uPrefixLengthB) = scoreBest;\r
 370                         }\r
 371                         }\r
 372                 }\r
 373 \r
 374 // Special case: close gaps at end of alignment\r
 375         DPD(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose;\r
 376         DPE(uLengthA, uLengthB) += PA[uLengthA-1].m_scoreGapClose2;\r
 377 \r
 378         DPI(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose;\r
 379         DPJ(uLengthA, uLengthB) += PB[uLengthB-1].m_scoreGapClose2;\r
 380 \r
 381 #if TRACE\r
 382         Log("DPL:\n");\r
 383         ListDP(DPL_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 384         Log("DPM:\n");\r
 385         ListDP(DPM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 386         Log("DPD:\n");\r
 387         ListDP(DPD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 388         Log("DPE:\n");\r
 389         ListDP(DPE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 390         Log("DPI:\n");\r
 391         ListDP(DPI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 392         Log("DPJ:\n");\r
 393         ListDP(DPJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 394         Log("TBM:\n");\r
 395         ListTB(TBM_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 396         Log("TBD:\n");\r
 397         ListTB(TBD_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 398         Log("TBE:\n");\r
 399         ListTB(TBE_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 400         Log("TBI:\n");\r
 401         ListTB(TBI_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 402         Log("TBJ:\n");\r
 403         ListTB(TBJ_, PA, PB, uPrefixCountA, uPrefixCountB);\r
 404 #endif\r
 405 \r
 406 // ==========\r
 407 // Trace-back\r
 408 // ==========\r
 409 \r
 410         Path.Clear();\r
 411 \r
 412 // Find last edge\r
 413         char cEdgeType = '?';\r
 414         SCORE BestScore = MINUS_INFINITY;\r
 415         SCORE M = DPM(uLengthA, uLengthB);\r
 416         SCORE D = DPD(uLengthA, uLengthB);\r
 417         SCORE E = DPE(uLengthA, uLengthB);\r
 418         SCORE I = DPI(uLengthA, uLengthB);\r
 419         SCORE J = DPJ(uLengthA, uLengthB);\r
 420 \r
 421         if (M >= D && M >= E && M >= I && M >= J)\r
 422                 {\r
 423                 cEdgeType = 'M';\r
 424                 BestScore = M;\r
 425                 }\r
 426         else if (D >= M && D >= E && D >= I && D >= J)\r
 427                 {\r
 428                 cEdgeType = 'D';\r
 429                 BestScore = D;\r
 430                 }\r
 431         else if (E >= M && E >= D && E >= I && E >= J)\r
 432                 {\r
 433                 cEdgeType = 'E';\r
 434                 BestScore = E;\r
 435                 }\r
 436         else if (I >= M && I >= D && I >= E && I >= J)\r
 437                 {\r
 438                 cEdgeType = 'I';\r
 439                 BestScore = I;\r
 440                 }\r
 441         else if (J >= M && J >= D && J >= E && J >= I)\r
 442                 {\r
 443                 cEdgeType = 'J';\r
 444                 BestScore = J;\r
 445                 }\r
 446         else\r
 447                 Quit("Bad max");\r
 448 \r
 449         unsigned PLA = uLengthA;\r
 450         unsigned PLB = uLengthB;\r
 451         unsigned ECount = 0;\r
 452         unsigned JCount = 0;\r
 453         for (;;)\r
 454                 {\r
 455 #if     TRACE\r
 456                 Log("TraceBack: %c%u.%u\n", cEdgeType, PLA, PLB);\r
 457 #endif\r
 458                 PWEdge Edge;\r
 459                 Edge.cType = XlatEdgeType(cEdgeType);\r
 460                 Edge.uPrefixLengthA = PLA;\r
 461                 Edge.uPrefixLengthB = PLB;\r
 462                 Path.PrependEdge(Edge);\r
 463 \r
 464                 switch (cEdgeType)\r
 465                         {\r
 466                 case 'M':\r
 467                         assert(PLA > 0);\r
 468                         assert(PLB > 0);\r
 469                         cEdgeType = TBM(PLA, PLB);\r
 470                         --PLA;\r
 471                         --PLB;\r
 472                         break;\r
 473 \r
 474                 case 'D':\r
 475                         assert(PLA > 0);\r
 476                         cEdgeType = TBD(PLA, PLB);\r
 477                         --PLA;\r
 478                         break;\r
 479 \r
 480                 case 'E':\r
 481                         ++ECount;\r
 482                         assert(PLA > 0);\r
 483                         cEdgeType = TBE(PLA, PLB);\r
 484                         --PLA;\r
 485                         break;\r
 486 \r
 487                 case 'I':\r
 488                         assert(PLB > 0);\r
 489                         cEdgeType = TBI(PLA, PLB);\r
 490                         --PLB;\r
 491                         break;\r
 492                 \r
 493                 case 'J':\r
 494                         ++JCount;\r
 495                         assert(PLB > 0);\r
 496                         cEdgeType = TBJ(PLA, PLB);\r
 497                         --PLB;\r
 498                         break;\r
 499 \r
 500                 default:\r
 501                         Quit("Invalid edge %c", cEdgeType);\r
 502                         }\r
 503                 if (0 == PLA && 0 == PLB)\r
 504                         break;\r
 505                 }\r
 506         //if (ECount > 0 || JCount > 0)\r
 507         //      fprintf(stderr, "E=%d J=%d\n", ECount, JCount);\r
 508         Path.Validate();\r
 509         if (Path.GetMatchCount() + Path.GetDeleteCount() != uLengthA)\r
 510                 Quit("Path count A");\r
 511         if (Path.GetMatchCount() + Path.GetInsertCount() != uLengthB)\r
 512                 Quit("Path count B");\r
 513 \r
 514         if (g_bKeepSimpleDP)\r
 515                 {\r
 516                 g_DPM = DPM_;\r
 517                 g_DPD = DPD_;\r
 518                 g_DPE = DPE_;\r
 519                 g_DPI = DPI_;\r
 520                 g_DPJ = DPJ_;\r
 521 \r
 522                 g_TBM = TBM_;\r
 523                 g_TBD = TBD_;\r
 524                 g_TBE = TBE_;\r
 525                 g_TBI = TBI_;\r
 526                 g_TBJ = TBJ_;\r
 527                 }\r
 528         else\r
 529                 {\r
 530                 delete[] DPM_;\r
 531                 delete[] DPD_;\r
 532                 delete[] DPE_;\r
 533                 delete[] DPI_;\r
 534                 delete[] DPJ_;\r
 535 \r
 536                 delete[] TBM_;\r
 537                 delete[] TBD_;\r
 538                 delete[] TBE_;\r
 539                 delete[] TBI_;\r
 540                 delete[] TBJ_;\r
 541                 }\r
 542 \r
 543 #if     TRACE\r
 544         Log("BestScore=%.6g\n", BestScore);\r
 545 #endif\r
 546         return BestScore;\r
 547         }\r
 548 \r
 549 #endif  // DOUBLE_AFFINE\r