--- /dev/null
+#include "muscle.h"\r
+#include "profile.h"\r
+#include "pwpath.h"\r
+\r
+struct DP_MEMORY\r
+ {\r
+ unsigned uLength;\r
+ SCORE *GapOpenA;\r
+ SCORE *GapOpenB;\r
+ SCORE *GapCloseA;\r
+ SCORE *GapCloseB;\r
+ SCORE *MPrev;\r
+ SCORE *MCurr;\r
+ SCORE *MWork;\r
+ SCORE *DPrev;\r
+ SCORE *DCurr;\r
+ SCORE *DWork;\r
+ SCORE **ScoreMxB;\r
+ unsigned **SortOrderA;\r
+ unsigned *uDeletePos;\r
+ FCOUNT **FreqsA;\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
+ {\r
+ delete[] DPM.TraceBack[i];\r
+ delete[] DPM.FreqsA[i];\r
+ delete[] DPM.SortOrderA[i];\r
+ }\r
+ for (unsigned n = 0; n < 20; ++n)\r
+ delete[] DPM.ScoreMxB[n];\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.uDeletePos;\r
+ delete[] DPM.GapOpenA;\r
+ delete[] DPM.GapOpenB;\r
+ delete[] DPM.GapCloseA;\r
+ delete[] DPM.GapCloseB;\r
+ delete[] DPM.SortOrderA;\r
+ delete[] DPM.FreqsA;\r
+ delete[] DPM.ScoreMxB;\r
+ delete[] DPM.TraceBack;\r
+ }\r
+\r
+ DPM.uLength = uLength;\r
+\r
+ DPM.GapOpenA = new SCORE[uLength];\r
+ DPM.GapOpenB = new SCORE[uLength];\r
+ DPM.GapCloseA = new SCORE[uLength];\r
+ DPM.GapCloseB = new SCORE[uLength];\r
+\r
+ DPM.SortOrderA = new unsigned*[uLength];\r
+ DPM.FreqsA = new FCOUNT*[uLength];\r
+ DPM.ScoreMxB = new SCORE*[20];\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.uDeletePos = new unsigned[uLength];\r
+\r
+ DPM.TraceBack = new int*[uLength];\r
+\r
+ for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
+ DPM.ScoreMxB[uLetter] = new SCORE[uLength];\r
+\r
+ for (unsigned i = 0; i < uLength; ++i)\r
+ {\r
+ DPM.SortOrderA[i] = new unsigned[20];\r
+ DPM.FreqsA[i] = new FCOUNT[20];\r
+ DPM.TraceBack[i] = new int[uLength];\r
+ }\r
+ }\r
+\r
+SCORE GlobalAlignNS(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
+ unsigned uLengthB, PWPath &Path)\r
+ {\r
+ const unsigned uPrefixCountA = uLengthA + 1;\r
+ const unsigned uPrefixCountB = uLengthB + 1;\r
+\r
+ AllocDPMem(uLengthA, uLengthB);\r
+\r
+ SCORE *GapOpenA = DPM.GapOpenA;\r
+ SCORE *GapOpenB = DPM.GapOpenB;\r
+ SCORE *GapCloseA = DPM.GapCloseA;\r
+ SCORE *GapCloseB = DPM.GapCloseB;\r
+\r
+ unsigned **SortOrderA = DPM.SortOrderA;\r
+ FCOUNT **FreqsA = DPM.FreqsA;\r
+ SCORE **ScoreMxB = DPM.ScoreMxB;\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
+ unsigned *uDeletePos = DPM.uDeletePos;\r
+\r
+ int **TraceBack = DPM.TraceBack;\r
+\r
+ for (unsigned i = 0; i < uLengthA; ++i)\r
+ {\r
+ GapOpenA[i] = PA[i].m_scoreGapOpen;\r
+ GapCloseA[i] = PA[i].m_scoreGapClose;\r
+\r
+ for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
+ {\r
+ SortOrderA[i][uLetter] = PA[i].m_uSortOrder[uLetter];\r
+ FreqsA[i][uLetter] = PA[i].m_fcCounts[uLetter];\r
+ }\r
+ }\r
+\r
+ for (unsigned j = 0; j < uLengthB; ++j)\r
+ {\r
+ GapOpenB[j] = PB[j].m_scoreGapOpen;\r
+ GapCloseB[j] = PB[j].m_scoreGapClose;\r
+ }\r
+\r
+ for (unsigned uLetter = 0; uLetter < 20; ++uLetter)\r
+ {\r
+ for (unsigned j = 0; j < uLengthB; ++j)\r
+ ScoreMxB[uLetter][j] = PB[j].m_AAScores[uLetter];\r
+ }\r
+\r
+ for (unsigned i = 0; i < uPrefixCountA; ++i)\r
+ memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));\r
+\r
+// Special case for i=0\r
+ unsigned **ptrSortOrderA = SortOrderA;\r
+ FCOUNT **ptrFreqsA = FreqsA;\r
+ assert(ptrSortOrderA == &(SortOrderA[0]));\r
+ assert(ptrFreqsA == &(FreqsA[0]));\r
+ TraceBack[0][0] = 0;\r
+\r
+ SCORE scoreSum = 0;\r
+ unsigned *ptrSortOrderAi = SortOrderA[0];\r
+ const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 20;\r
+ FCOUNT *ptrFreqsAi = FreqsA[0];\r
+ for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi)\r
+ {\r
+ const unsigned uLetter = *ptrSortOrderAi;\r
+ const FCOUNT fcLetter = ptrFreqsAi[uLetter];\r
+ if (0 == fcLetter)\r
+ break;\r
+ scoreSum += fcLetter*ScoreMxB[uLetter][0];\r
+ }\r
+ MPrev[0] = scoreSum - g_scoreCenter;\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
+ // 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
+ SCORE scoreSum = 0;\r
+ unsigned *ptrSortOrderAi = SortOrderA[0];\r
+ const unsigned *ptrSortOrderAEnd = ptrSortOrderAi + 20;\r
+ FCOUNT *ptrFreqsAi = FreqsA[0];\r
+ for (; ptrSortOrderAi != ptrSortOrderAEnd; ++ptrSortOrderAi)\r
+ {\r
+ const unsigned uLetter = *ptrSortOrderAi;\r
+ const FCOUNT fcLetter = ptrFreqsAi[uLetter];\r
+ if (0 == fcLetter)\r
+ break;\r
+ scoreSum += fcLetter*ScoreMxB[uLetter][j];\r
+ }\r
+ MPrev[j] = scoreSum - g_scoreCenter + GapOpenB[0] + GapCloseB[j-1];\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
+ ++ptrSortOrderA;\r
+ ++ptrFreqsA;\r
+ assert(ptrSortOrderA == &(SortOrderA[i]));\r
+ assert(ptrFreqsA == &(FreqsA[i]));\r
+\r
+ SCORE *ptrMCurr_j = MCurr;\r
+ memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));\r
+ const FCOUNT *FreqsAi = *ptrFreqsA;\r
+\r
+ const unsigned *SortOrderAi = *ptrSortOrderA;\r
+ const unsigned *ptrSortOrderAiEnd = SortOrderAi + 20;\r
+ const SCORE *ptrMCurrMax = MCurr + uLengthB;\r
+ for (const unsigned *ptrSortOrderAi = SortOrderAi;\r
+ ptrSortOrderAi != ptrSortOrderAiEnd;\r
+ ++ptrSortOrderAi)\r
+ {\r
+ const unsigned uLetter = *ptrSortOrderAi;\r
+ SCORE *NSBR_Letter = ScoreMxB[uLetter];\r
+ const FCOUNT fcLetter = FreqsAi[uLetter];\r
+ if (0 == fcLetter)\r
+ break;\r
+ SCORE *ptrNSBR = NSBR_Letter;\r
+ for (SCORE *ptrMCurr = MCurr; ptrMCurr != ptrMCurrMax; ++ptrMCurr)\r
+ *ptrMCurr += fcLetter*(*ptrNSBR++);\r
+ }\r
+\r
+ for (unsigned j = 0; j < uLengthB; ++j)\r
+ MCurr[j] -= g_scoreCenter;\r
+\r
+ ptrMCurr_j = MCurr;\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
+ assert(ptrMCurr_j == &(MCurr[0]));\r
+ *ptrMCurr_j += GapOpenA[0] + GapCloseA[i-1];\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 + GapOpenA[i];\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
+ const SCORE scoreGapOpenAi = GapOpenA[i];\r
+ const SCORE scoreGapCloseAi_1 = GapCloseA[i-1];\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 + GapOpenB[j];\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++ + scoreGapCloseAi_1;\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 + GapCloseB[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
+ assert(ptrSortOrderA == &(SortOrderA[i]));\r
+ assert(ptrFreqsA == &(FreqsA[i]));\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 + scoreGapOpenAi;\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] + GapOpenB[j];\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] + GapCloseA[uLengthA-1];\r
+ if (scoreD > scoreMax)\r
+ {\r
+ scoreMax = scoreD;\r
+ iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];\r
+ }\r
+\r
+ SCORE scoreI = IPrev + GapCloseB[uLengthB-1];\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