+++ /dev/null
-#include "muscle.h"\r
-#include "profile.h"\r
-#include "pwpath.h"\r
-#include "seq.h"\r
-\r
-extern SCOREMATRIX VTML_SP;\r
-\r
-// #define SUBST(i, j) Subst(seqA, seqB, i, j)\r
-#define SUBST(i, j) MxRowA[i][seqB.GetLetter(j)]\r
-\r
-static SCORE Subst(const Seq &seqA, const Seq &seqB, unsigned i, unsigned j)\r
- {\r
- assert(i < seqA.Length());\r
- assert(j < seqB.Length());\r
-\r
- unsigned uLetterA = seqA.GetLetter(i);\r
- unsigned uLetterB = seqB.GetLetter(j);\r
- return VTML_SP[uLetterA][uLetterB] + g_scoreCenter;\r
- }\r
-\r
-struct DP_MEMORY\r
- {\r
- unsigned uLength;\r
- SCORE *MPrev;\r
- SCORE *MCurr;\r
- SCORE *MWork;\r
- SCORE *DPrev;\r
- SCORE *DCurr;\r
- SCORE *DWork;\r
- SCORE **MxRowA;\r
- unsigned *LettersB;\r
- unsigned *uDeletePos;\r
- int **TraceBack;\r
- };\r
-\r
-static struct DP_MEMORY DPM;\r
-\r
-static void AllocDPMem(unsigned uLengthA, unsigned uLengthB)\r
- {\r
-// Max prefix length\r
- unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1;\r
- if (uLength < DPM.uLength)\r
- return;\r
-\r
-// Add 256 to allow for future expansion and\r
-// round up to next multiple of 32.\r
- uLength += 256;\r
- uLength += 32 - uLength%32;\r
-\r
- const unsigned uOldLength = DPM.uLength;\r
- if (uOldLength > 0)\r
- {\r
- for (unsigned i = 0; i < uOldLength; ++i)\r
- delete[] DPM.TraceBack[i];\r
-\r
- delete[] DPM.MPrev;\r
- delete[] DPM.MCurr;\r
- delete[] DPM.MWork;\r
- delete[] DPM.DPrev;\r
- delete[] DPM.DCurr;\r
- delete[] DPM.DWork;\r
- delete[] DPM.MxRowA;\r
- delete[] DPM.LettersB;\r
- delete[] DPM.uDeletePos;\r
- delete[] DPM.TraceBack;\r
- }\r
-\r
- DPM.uLength = uLength;\r
-\r
- DPM.MPrev = new SCORE[uLength];\r
- DPM.MCurr = new SCORE[uLength];\r
- DPM.MWork = new SCORE[uLength];\r
-\r
- DPM.DPrev = new SCORE[uLength];\r
- DPM.DCurr = new SCORE[uLength];\r
- DPM.DWork = new SCORE[uLength];\r
- DPM.MxRowA = new SCORE *[uLength];\r
- DPM.LettersB = new unsigned[uLength];\r
- DPM.uDeletePos = new unsigned[uLength];\r
-\r
- DPM.TraceBack = new int*[uLength];\r
-\r
- for (unsigned i = 0; i < uLength; ++i)\r
- DPM.TraceBack[i] = new int[uLength];\r
- }\r
-\r
-static void RowFromSeq(const Seq &s, SCORE *Row[])\r
- {\r
- const unsigned uLength = s.Length();\r
- for (unsigned i = 0; i < uLength; ++i)\r
- {\r
- char c = s.GetChar(i);\r
- unsigned uLetter = CharToLetter(c);\r
- if (uLetter < 20)\r
- Row[i] = VTML_SP[uLetter];\r
- else\r
- Row[i] = VTML_SP[AX_X];\r
- }\r
- }\r
-\r
-static void LettersFromSeq(const Seq &s, unsigned Letters[])\r
- {\r
- const unsigned uLength = s.Length();\r
- for (unsigned i = 0; i < uLength; ++i)\r
- {\r
- char c = s.GetChar(i);\r
- unsigned uLetter = CharToLetter(c);\r
- if (uLetter < 20)\r
- Letters[i] = uLetter;\r
- else\r
- Letters[i] = AX_X;\r
- }\r
- }\r
-\r
-SCORE GlobalAlignSS(const Seq &seqA, const Seq &seqB, PWPath &Path)\r
- {\r
- const unsigned uLengthA = seqA.Length();\r
- const unsigned uLengthB = seqB.Length();\r
- const unsigned uPrefixCountA = uLengthA + 1;\r
- const unsigned uPrefixCountB = uLengthB + 1;\r
-\r
- AllocDPMem(uLengthA, uLengthB);\r
-\r
- SCORE *MPrev = DPM.MPrev;\r
- SCORE *MCurr = DPM.MCurr;\r
- SCORE *MWork = DPM.MWork;\r
-\r
- SCORE *DPrev = DPM.DPrev;\r
- SCORE *DCurr = DPM.DCurr;\r
- SCORE *DWork = DPM.DWork;\r
- SCORE **MxRowA = DPM.MxRowA;\r
- unsigned *LettersB = DPM.LettersB;\r
-\r
- RowFromSeq(seqA, MxRowA);\r
- LettersFromSeq(seqB, LettersB);\r
-\r
- unsigned *uDeletePos = DPM.uDeletePos;\r
-\r
- int **TraceBack = DPM.TraceBack;\r
-\r
-#if DEBUG\r
- for (unsigned i = 0; i < uPrefixCountA; ++i)\r
- memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));\r
-#endif\r
-\r
-// Special case for i=0\r
- TraceBack[0][0] = 0;\r
- MPrev[0] = MxRowA[0][LettersB[0]];\r
-\r
-// D(0,0) is -infinity (requires I->D).\r
- DPrev[0] = MINUS_INFINITY;\r
-\r
- for (unsigned j = 1; j < uLengthB; ++j)\r
- {\r
- unsigned uLetterB = LettersB[j];\r
-\r
- // Only way to get M(0, j) looks like this:\r
- // A ----X\r
- // B XXXXX\r
- // 0 j\r
- // So gap-open at j=0, gap-close at j-1.\r
- MPrev[j] = MxRowA[0][uLetterB] + g_scoreGapOpen/2; // term gaps half\r
- TraceBack[0][j] = -(int) j;\r
-\r
- // Assume no D->I transitions, then can't be a delete if only\r
- // one letter from A.\r
- DPrev[j] = MINUS_INFINITY;\r
- }\r
-\r
- SCORE IPrev_j_1;\r
- for (unsigned i = 1; i < uLengthA; ++i)\r
- {\r
- SCORE *ptrMCurr_j = MCurr;\r
- memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));\r
-\r
- const SCORE *RowA = MxRowA[i];\r
- const SCORE *ptrRowA = MxRowA[i];\r
- const SCORE *ptrMCurrEnd = ptrMCurr_j + uLengthB;\r
- unsigned *ptrLettersB = LettersB;\r
- for (; ptrMCurr_j != ptrMCurrEnd; ++ptrMCurr_j)\r
- {\r
- *ptrMCurr_j = RowA[*ptrLettersB];\r
- ++ptrLettersB;\r
- }\r
-\r
- unsigned *ptrDeletePos = uDeletePos;\r
-\r
- // Special case for j=0\r
- // Only way to get M(i, 0) looks like this:\r
- // 0 i\r
- // A XXXXX\r
- // B ----X\r
- // So gap-open at i=0, gap-close at i-1.\r
- ptrMCurr_j = MCurr;\r
- assert(ptrMCurr_j == &(MCurr[0]));\r
- *ptrMCurr_j += g_scoreGapOpen/2; // term gaps half\r
-\r
- ++ptrMCurr_j;\r
-\r
- int *ptrTraceBack_ij = TraceBack[i];\r
- *ptrTraceBack_ij++ = (int) i;\r
-\r
- SCORE *ptrMPrev_j = MPrev;\r
- SCORE *ptrDPrev = DPrev;\r
- SCORE d = *ptrDPrev;\r
- SCORE DNew = *ptrMPrev_j + g_scoreGapOpen;\r
- if (DNew > d)\r
- {\r
- d = DNew;\r
- *ptrDeletePos = i;\r
- }\r
-\r
- SCORE *ptrDCurr = DCurr;\r
-\r
- assert(ptrDCurr == &(DCurr[0]));\r
- *ptrDCurr = d;\r
-\r
- // Can't have an insert if no letters from B\r
- IPrev_j_1 = MINUS_INFINITY;\r
-\r
- unsigned uInsertPos;\r
-\r
- for (unsigned j = 1; j < uLengthB; ++j)\r
- {\r
- // Here, MPrev_j is preserved from previous\r
- // iteration so with current i,j is M[i-1][j-1]\r
- SCORE MPrev_j = *ptrMPrev_j;\r
- SCORE INew = MPrev_j + g_scoreGapOpen;\r
- if (INew > IPrev_j_1)\r
- {\r
- IPrev_j_1 = INew;\r
- uInsertPos = j;\r
- }\r
-\r
- SCORE scoreMax = MPrev_j;\r
-\r
- assert(ptrDPrev == &(DPrev[j-1]));\r
- SCORE scoreD = *ptrDPrev++;\r
- if (scoreD > scoreMax)\r
- {\r
- scoreMax = scoreD;\r
- assert(ptrDeletePos == &(uDeletePos[j-1]));\r
- *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos;\r
- assert(*ptrTraceBack_ij > 0);\r
- }\r
- ++ptrDeletePos;\r
-\r
- SCORE scoreI = IPrev_j_1;\r
- if (scoreI > scoreMax)\r
- {\r
- scoreMax = scoreI;\r
- *ptrTraceBack_ij = (int) uInsertPos - (int) j;\r
- assert(*ptrTraceBack_ij < 0);\r
- }\r
-\r
- *ptrMCurr_j += scoreMax;\r
- assert(ptrMCurr_j == &(MCurr[j]));\r
- ++ptrMCurr_j;\r
-\r
- MPrev_j = *(++ptrMPrev_j);\r
- assert(ptrDPrev == &(DPrev[j]));\r
- SCORE d = *ptrDPrev;\r
- SCORE DNew = MPrev_j + g_scoreGapOpen;\r
- if (DNew > d)\r
- {\r
- d = DNew;\r
- assert(ptrDeletePos == &uDeletePos[j]);\r
- *ptrDeletePos = i;\r
- }\r
- assert(ptrDCurr + 1 == &(DCurr[j]));\r
- *(++ptrDCurr) = d;\r
-\r
- ++ptrTraceBack_ij;\r
- }\r
-\r
- Rotate(MPrev, MCurr, MWork);\r
- Rotate(DPrev, DCurr, DWork);\r
- }\r
-\r
-// Special case for i=uLengthA\r
- SCORE IPrev = MINUS_INFINITY;\r
-\r
- unsigned uInsertPos;\r
-\r
- for (unsigned j = 1; j < uLengthB; ++j)\r
- {\r
- SCORE INew = MPrev[j-1];\r
- if (INew > IPrev)\r
- {\r
- uInsertPos = j;\r
- IPrev = INew;\r
- }\r
- }\r
-\r
-// Special case for i=uLengthA, j=uLengthB\r
- SCORE scoreMax = MPrev[uLengthB-1];\r
- int iTraceBack = 0;\r
-\r
- SCORE scoreD = DPrev[uLengthB-1] - g_scoreGapOpen/2; // term gaps half\r
- if (scoreD > scoreMax)\r
- {\r
- scoreMax = scoreD;\r
- iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];\r
- }\r
-\r
- SCORE scoreI = IPrev - g_scoreGapOpen/2;\r
- if (scoreI > scoreMax)\r
- {\r
- scoreMax = scoreI;\r
- iTraceBack = (int) uInsertPos - (int) uLengthB;\r
- }\r
-\r
- TraceBack[uLengthA][uLengthB] = iTraceBack;\r
-\r
- TraceBackToPath(TraceBack, uLengthA, uLengthB, Path);\r
-\r
- return scoreMax;\r
- }\r
git://source.jalview.org / jabaws.git / blobdiff