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

   1 #include "muscle.h"\r
   2 #include "profile.h"\r
   3 #include "pwpath.h"\r
   4 \r
   5 #define OCC     1\r
   6 \r
   7 struct DP_MEMORY\r
   8         {\r
   9         unsigned uLength;\r
  10         SCORE *GapOpenA;\r
  11         SCORE *GapOpenB;\r
  12         SCORE *GapCloseA;\r
  13         SCORE *GapCloseB;\r
  14         SCORE *MPrev;\r
  15         SCORE *MCurr;\r
  16         SCORE *MWork;\r
  17         SCORE *DPrev;\r
  18         SCORE *DCurr;\r
  19         SCORE *DWork;\r
  20         SCORE **ScoreMxB;\r
  21 #if     OCC\r
  22         FCOUNT *OccA;\r
  23         FCOUNT *OccB;\r
  24 #endif\r
  25         unsigned **SortOrderA;\r
  26         unsigned *uDeletePos;\r
  27         FCOUNT **FreqsA;\r
  28         int **TraceBack;\r
  29         };\r
  30 \r
  31 static struct DP_MEMORY DPM;\r
  32 \r
  33 static void AllocDPMem(unsigned uLengthA, unsigned uLengthB)\r
  34         {\r
  35 // Max prefix length\r
  36         unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1;\r
  37         if (uLength < DPM.uLength)\r
  38                 return;\r
  39 \r
  40 // Add 256 to allow for future expansion and\r
  41 // round up to next multiple of 32.\r
  42         uLength += 256;\r
  43         uLength += 32 - uLength%32;\r
  44 \r
  45         const unsigned uOldLength = DPM.uLength;\r
  46         if (uOldLength > 0)\r
  47                 {\r
  48                 for (unsigned i = 0; i < uOldLength; ++i)\r
  49                         {\r
  50                         delete[] DPM.TraceBack[i];\r
  51                         delete[] DPM.FreqsA[i];\r
  52                         delete[] DPM.SortOrderA[i];\r
  53                         }\r
  54                 for (unsigned n = 0; n < 20; ++n)\r
  55                         delete[] DPM.ScoreMxB[n];\r
  56 \r
  57                 delete[] DPM.MPrev;\r
  58                 delete[] DPM.MCurr;\r
  59                 delete[] DPM.MWork;\r
  60                 delete[] DPM.DPrev;\r
  61                 delete[] DPM.DCurr;\r
  62                 delete[] DPM.DWork;\r
  63                 delete[] DPM.uDeletePos;\r
  64                 delete[] DPM.GapOpenA;\r
  65                 delete[] DPM.GapOpenB;\r
  66                 delete[] DPM.GapCloseA;\r
  67                 delete[] DPM.GapCloseB;\r
  68                 delete[] DPM.SortOrderA;\r
  69                 delete[] DPM.FreqsA;\r
  70                 delete[] DPM.ScoreMxB;\r
  71                 delete[] DPM.TraceBack;\r
  72 #if     OCC\r
  73                 delete[] DPM.OccA;\r
  74                 delete[] DPM.OccB;\r
  75 #endif\r
  76                 }\r
  77 \r
  78         DPM.uLength = uLength;\r
  79 \r
  80         DPM.GapOpenA = new SCORE[uLength];\r
  81         DPM.GapOpenB = new SCORE[uLength];\r
  82         DPM.GapCloseA = new SCORE[uLength];\r
  83         DPM.GapCloseB = new SCORE[uLength];\r
  84 #if     OCC\r
  85         DPM.OccA = new FCOUNT[uLength];\r
  86         DPM.OccB = new FCOUNT[uLength];\r
  87 #endif\r
  88 \r
  89         DPM.SortOrderA = new unsigned*[uLength];\r
  90         DPM.FreqsA = new FCOUNT*[uLength];\r
  91         DPM.ScoreMxB = new SCORE*[20];\r
  92         DPM.MPrev = new SCORE[uLength];\r
  93         DPM.MCurr = new SCORE[uLength];\r
  94         DPM.MWork = new SCORE[uLength];\r
  95 \r
  96         DPM.DPrev = new SCORE[uLength];\r
  97         DPM.DCurr = new SCORE[uLength];\r
  98         DPM.DWork = new SCORE[uLength];\r
  99         DPM.uDeletePos = new unsigned[uLength];\r
 100 \r
 101         DPM.TraceBack = new int*[uLength];\r
 102 \r
 103         for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
 104                 DPM.ScoreMxB[uLetter] = new SCORE[uLength];\r
 105 \r
 106         for (unsigned i = 0; i < uLength; ++i)\r
 107                 {\r
 108                 DPM.SortOrderA[i] = new unsigned[20];\r
 109                 DPM.FreqsA[i] = new FCOUNT[20];\r
 110                 DPM.TraceBack[i] = new int[uLength];\r
 111                 }\r
 112         }\r
 113 \r
 114 SCORE GlobalAlignLE(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
 115   unsigned uLengthB, PWPath &Path)\r
 116         {\r
 117         SetTermGaps(PA, uLengthA);\r
 118         SetTermGaps(PB, uLengthB);\r
 119 \r
 120         const unsigned uPrefixCountA = uLengthA + 1;\r
 121         const unsigned uPrefixCountB = uLengthB + 1;\r
 122 \r
 123         AllocDPMem(uLengthA, uLengthB);\r
 124 \r
 125         SCORE *GapOpenA = DPM.GapOpenA;\r
 126         SCORE *GapOpenB = DPM.GapOpenB;\r
 127         SCORE *GapCloseA = DPM.GapCloseA;\r
 128         SCORE *GapCloseB = DPM.GapCloseB;\r
 129 \r
 130         unsigned **SortOrderA = DPM.SortOrderA;\r
 131         FCOUNT **FreqsA = DPM.FreqsA;\r
 132         SCORE **ScoreMxB = DPM.ScoreMxB;\r
 133         SCORE *MPrev = DPM.MPrev;\r
 134         SCORE *MCurr = DPM.MCurr;\r
 135         SCORE *MWork = DPM.MWork;\r
 136 \r
 137         SCORE *DPrev = DPM.DPrev;\r
 138         SCORE *DCurr = DPM.DCurr;\r
 139         SCORE *DWork = DPM.DWork;\r
 140 \r
 141 #if     OCC\r
 142         FCOUNT *OccA = DPM.OccA;\r
 143         FCOUNT *OccB = DPM.OccB;\r
 144 #endif\r
 145 \r
 146         unsigned *uDeletePos = DPM.uDeletePos;\r
 147 \r
 148         int **TraceBack = DPM.TraceBack;\r
 149 \r
 150         for (unsigned i = 0; i < uLengthA; ++i)\r
 151                 {\r
 152                 GapOpenA[i] = PA[i].m_scoreGapOpen;\r
 153                 GapCloseA[i] = PA[i].m_scoreGapClose;\r
 154 #if     OCC\r
 155                 OccA[i] = PA[i].m_fOcc;\r
 156 #endif\r
 157 \r
 158                 for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
 159                         {\r
 160                         SortOrderA[i][uLetter] = PA[i].m_uSortOrder[uLetter];\r
 161                         FreqsA[i][uLetter] = PA[i].m_fcCounts[uLetter];\r
 162                         }\r
 163                 }\r
 164 \r
 165         for (unsigned j = 0; j < uLengthB; ++j)\r
 166                 {\r
 167                 GapOpenB[j] = PB[j].m_scoreGapOpen;\r
 168                 GapCloseB[j] = PB[j].m_scoreGapClose;\r
 169 #if     OCC\r
 170                 OccB[j] = PB[j].m_fOcc;\r
 171 #endif\r
 172                 }\r
 173 \r
 174         for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
 175                 {\r
 176                 for (unsigned j = 0; j < uLengthB; ++j)\r
 177                         ScoreMxB[uLetter][j] = PB[j].m_AAScores[uLetter];\r
 178                 }\r
 179 \r
 180         for (unsigned i = 0; i < uPrefixCountA; ++i)\r
 181                 memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));\r
 182 \r
 183 // Special case for i=0\r
 184         unsigned **ptrSortOrderA = SortOrderA;\r
 185         FCOUNT **ptrFreqsA = FreqsA;\r
 186         assert(ptrSortOrderA == &(SortOrderA[0]));\r
 187         assert(ptrFreqsA == &(FreqsA[0]));\r
 188         TraceBack[0][0] = 0;\r
 189 \r
 190         SCORE scoreSum = 0;\r
 191         unsigned *ptrSortOrderAi = SortOrderA[0];\r
 192         const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 20;\r
 193         FCOUNT *ptrFreqsAi = FreqsA[0];\r
 194         for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi)\r
 195                 {\r
 196                 const unsigned uLetter = *ptrSortOrderAi;\r
 197                 const FCOUNT fcLetter = ptrFreqsAi[uLetter];\r
 198                 if (0 == fcLetter)\r
 199                         break;\r
 200                 scoreSum += fcLetter*ScoreMxB[uLetter][0];\r
 201                 }\r
 202         if (0 == scoreSum)\r
 203                 MPrev[0] = -2.5;\r
 204         else\r
 205                 {\r
 206 #if     OCC\r
 207                 MPrev[0] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[0];\r
 208 #else\r
 209                 MPrev[0] = (logf(scoreSum) - g_scoreCenter);\r
 210 #endif\r
 211                 }\r
 212 \r
 213 // D(0,0) is -infinity (requires I->D).\r
 214         DPrev[0] = MINUS_INFINITY;\r
 215 \r
 216         for (unsigned j = 1; j < uLengthB; ++j)\r
 217                 {\r
 218         // Only way to get M(0, j) looks like this:\r
 219         //              A       ----X\r
 220         //              B       XXXXX\r
 221         //                      0   j\r
 222         // So gap-open at j=0, gap-close at j-1.\r
 223                 SCORE scoreSum = 0;\r
 224                 unsigned *ptrSortOrderAi = SortOrderA[0];\r
 225                 const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 20;\r
 226                 FCOUNT *ptrFreqsAi = FreqsA[0];\r
 227                 for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi)\r
 228                         {\r
 229                         const unsigned uLetter = *ptrSortOrderAi;\r
 230                         const FCOUNT fcLetter = ptrFreqsAi[uLetter];\r
 231                         if (0 == fcLetter)\r
 232                                 break;\r
 233                         scoreSum += fcLetter*ScoreMxB[uLetter][j];\r
 234                         }\r
 235                 if (0 == scoreSum)\r
 236                         MPrev[j] = -2.5;\r
 237                 else\r
 238                         {\r
 239 #if     OCC\r
 240                         MPrev[j] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[j] +\r
 241                           GapOpenB[0] + GapCloseB[j-1];\r
 242 #else\r
 243                         MPrev[j] = (logf(scoreSum) - g_scoreCenter) +\r
 244                           GapOpenB[0] + GapCloseB[j-1];\r
 245 #endif\r
 246                         }\r
 247                 TraceBack[0][j] = -(int) j;\r
 248 \r
 249         // Assume no D->I transitions, then can't be a delete if only\r
 250         // one letter from A.\r
 251                 DPrev[j] = MINUS_INFINITY;\r
 252                 }\r
 253 \r
 254         SCORE IPrev_j_1;\r
 255         for (unsigned i = 1; i < uLengthA; ++i)\r
 256                 {\r
 257                 ++ptrSortOrderA;\r
 258                 ++ptrFreqsA;\r
 259                 assert(ptrSortOrderA == &(SortOrderA[i]));\r
 260                 assert(ptrFreqsA == &(FreqsA[i]));\r
 261 \r
 262                 SCORE *ptrMCurr_j = MCurr;\r
 263                 memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));\r
 264                 const FCOUNT *FreqsAi = *ptrFreqsA;\r
 265 \r
 266                 const unsigned *SortOrderAi = *ptrSortOrderA;\r
 267                 const unsigned *ptrSortOrderAiEnd = SortOrderAi + 20;\r
 268                 const SCORE *ptrMCurrMax = MCurr + uLengthB;\r
 269                 for (const unsigned *ptrSortOrderAi = SortOrderAi;\r
 270                   ptrSortOrderAi != ptrSortOrderAiEnd;\r
 271                   ++ptrSortOrderAi)\r
 272                         {\r
 273                         const unsigned uLetter = *ptrSortOrderAi;\r
 274                         SCORE *NSBR_Letter = ScoreMxB[uLetter];\r
 275                         const FCOUNT fcLetter = FreqsAi[uLetter];\r
 276                         if (0 == fcLetter)\r
 277                                 break;\r
 278                         SCORE *ptrNSBR = NSBR_Letter;\r
 279                         for (SCORE *ptrMCurr = MCurr; ptrMCurr != ptrMCurrMax; ++ptrMCurr)\r
 280                                 *ptrMCurr += fcLetter*(*ptrNSBR++);\r
 281                         }\r
 282 \r
 283 #if     OCC\r
 284                 const FCOUNT OccAi = OccA[i];\r
 285 #endif\r
 286                 for (unsigned j = 0; j < uLengthB; ++j)\r
 287                         {\r
 288                         if (MCurr[j] == 0)\r
 289                                 MCurr[j] = -2.5;\r
 290                         else\r
 291 #if     OCC\r
 292                                 MCurr[j] = (logf(MCurr[j]) - g_scoreCenter)*OccAi*OccB[j];\r
 293 #else\r
 294                                 MCurr[j] = (logf(MCurr[j]) - g_scoreCenter);\r
 295 #endif\r
 296                         }\r
 297 \r
 298                 ptrMCurr_j = MCurr;\r
 299                 unsigned *ptrDeletePos = uDeletePos;\r
 300 \r
 301         // Special case for j=0\r
 302         // Only way to get M(i, 0) looks like this:\r
 303         //                      0   i\r
 304         //              A       XXXXX\r
 305         //              B       ----X\r
 306         // So gap-open at i=0, gap-close at i-1.\r
 307                 assert(ptrMCurr_j == &(MCurr[0]));\r
 308                 *ptrMCurr_j += GapOpenA[0] + GapCloseA[i-1];\r
 309 \r
 310                 ++ptrMCurr_j;\r
 311 \r
 312                 int *ptrTraceBack_ij = TraceBack[i];\r
 313                 *ptrTraceBack_ij++ = (int) i;\r
 314 \r
 315                 SCORE *ptrMPrev_j = MPrev;\r
 316                 SCORE *ptrDPrev = DPrev;\r
 317                 SCORE d = *ptrDPrev;\r
 318                 SCORE DNew = *ptrMPrev_j + GapOpenA[i];\r
 319                 if (DNew > d)\r
 320                         {\r
 321                         d = DNew;\r
 322                         *ptrDeletePos = i;\r
 323                         }\r
 324 \r
 325                 SCORE *ptrDCurr = DCurr;\r
 326 \r
 327                 assert(ptrDCurr == &(DCurr[0]));\r
 328                 *ptrDCurr = d;\r
 329 \r
 330         // Can't have an insert if no letters from B\r
 331                 IPrev_j_1 = MINUS_INFINITY;\r
 332 \r
 333                 unsigned uInsertPos = 0;\r
 334                 const SCORE scoreGapOpenAi = GapOpenA[i];\r
 335                 const SCORE scoreGapCloseAi_1 = GapCloseA[i-1];\r
 336 \r
 337                 for (unsigned j = 1; j < uLengthB; ++j)\r
 338                         {\r
 339                 // Here, MPrev_j is preserved from previous\r
 340                 // iteration so with current i,j is M[i-1][j-1]\r
 341                         SCORE MPrev_j = *ptrMPrev_j;\r
 342                         SCORE INew = MPrev_j + GapOpenB[j];\r
 343                         if (INew > IPrev_j_1)\r
 344                                 {\r
 345                                 IPrev_j_1 = INew;\r
 346                                 uInsertPos = j;\r
 347                                 }\r
 348 \r
 349                         SCORE scoreMax = MPrev_j;\r
 350 \r
 351                         assert(ptrDPrev == &(DPrev[j-1]));\r
 352                         SCORE scoreD = *ptrDPrev++ + scoreGapCloseAi_1;\r
 353                         if (scoreD > scoreMax)\r
 354                                 {\r
 355                                 scoreMax = scoreD;\r
 356                                 assert(ptrDeletePos == &(uDeletePos[j-1]));\r
 357                                 *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos;\r
 358                                 assert(*ptrTraceBack_ij > 0);\r
 359                                 }\r
 360                         ++ptrDeletePos;\r
 361 \r
 362                         SCORE scoreI = IPrev_j_1 + GapCloseB[j-1];\r
 363                         if (scoreI > scoreMax)\r
 364                                 {\r
 365                                 scoreMax = scoreI;\r
 366                                 *ptrTraceBack_ij = (int) uInsertPos - (int) j;\r
 367                                 assert(*ptrTraceBack_ij < 0);\r
 368                                 }\r
 369 \r
 370                         assert(ptrSortOrderA == &(SortOrderA[i]));\r
 371                         assert(ptrFreqsA == &(FreqsA[i]));\r
 372 \r
 373                         *ptrMCurr_j += scoreMax;\r
 374                         assert(ptrMCurr_j == &(MCurr[j]));\r
 375                         ++ptrMCurr_j;\r
 376 \r
 377                         MPrev_j = *(++ptrMPrev_j);\r
 378                         assert(ptrDPrev == &(DPrev[j]));\r
 379                         SCORE d = *ptrDPrev;\r
 380                         SCORE DNew = MPrev_j + scoreGapOpenAi;\r
 381                         if (DNew > d)\r
 382                                 {\r
 383                                 d = DNew;\r
 384                                 assert(ptrDeletePos == &uDeletePos[j]);\r
 385                                 *ptrDeletePos = i;\r
 386                                 }\r
 387                         assert(ptrDCurr + 1 == &(DCurr[j]));\r
 388                         *(++ptrDCurr) = d;\r
 389 \r
 390                         ++ptrTraceBack_ij;\r
 391                         }\r
 392 \r
 393                 Rotate(MPrev, MCurr, MWork);\r
 394                 Rotate(DPrev, DCurr, DWork);\r
 395                 }\r
 396 \r
 397 // Special case for i=uLengthA\r
 398         SCORE IPrev = MINUS_INFINITY;\r
 399 \r
 400         unsigned uInsertPos;\r
 401 \r
 402         for (unsigned j = 1; j < uLengthB; ++j)\r
 403                 {\r
 404                 SCORE INew = MPrev[j-1] + GapOpenB[j];\r
 405                 if (INew > IPrev)\r
 406                         {\r
 407                         uInsertPos = j;\r
 408                         IPrev = INew;\r
 409                         }\r
 410                 }\r
 411 \r
 412 // Special case for i=uLengthA, j=uLengthB\r
 413         SCORE scoreMax = MPrev[uLengthB-1];\r
 414         int iTraceBack = 0;\r
 415 \r
 416         SCORE scoreD = DPrev[uLengthB-1] + GapCloseA[uLengthA-1];\r
 417         if (scoreD > scoreMax)\r
 418                 {\r
 419                 scoreMax = scoreD;\r
 420                 iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];\r
 421                 }\r
 422 \r
 423         SCORE scoreI = IPrev + GapCloseB[uLengthB-1];\r
 424         if (scoreI > scoreMax)\r
 425                 {\r
 426                 scoreMax = scoreI;\r
 427                 iTraceBack = (int) uInsertPos - (int) uLengthB;\r
 428                 }\r
 429 \r
 430         TraceBack[uLengthA][uLengthB] = iTraceBack;\r
 431 \r
 432         TraceBackToPath(TraceBack, uLengthA, uLengthB, Path);\r
 433 \r
 434         return scoreMax;\r
 435         }\r