--- /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