+++ /dev/null
-#include "muscle.h"\r
-#include "profile.h"\r
-#include "pwpath.h"\r
-\r
-#define OCC 1\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
-#if OCC\r
- FCOUNT *OccA;\r
- FCOUNT *OccB;\r
-#endif\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
-#if OCC\r
- delete[] DPM.OccA;\r
- delete[] DPM.OccB;\r
-#endif\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
-#if OCC\r
- DPM.OccA = new FCOUNT[uLength];\r
- DPM.OccB = new FCOUNT[uLength];\r
-#endif\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 GlobalAlignLE(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,\r
- unsigned uLengthB, PWPath &Path)\r
- {\r
- SetTermGaps(PA, uLengthA);\r
- SetTermGaps(PB, uLengthB);\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
-\r
-#if OCC\r
- FCOUNT *OccA = DPM.OccA;\r
- FCOUNT *OccB = DPM.OccB;\r
-#endif\r
-\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
-#if OCC\r
- OccA[i] = PA[i].m_fOcc;\r
-#endif\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
-#if OCC\r
- OccB[j] = PB[j].m_fOcc;\r
-#endif\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
- if (0 == scoreSum)\r
- MPrev[0] = -2.5;\r
- else\r
- {\r
-#if OCC\r
- MPrev[0] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[0];\r
-#else\r
- MPrev[0] = (logf(scoreSum) - g_scoreCenter);\r
-#endif\r
- }\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
- if (0 == scoreSum)\r
- MPrev[j] = -2.5;\r
- else\r
- {\r
-#if OCC\r
- MPrev[j] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[j] +\r
- GapOpenB[0] + GapCloseB[j-1];\r
-#else\r
- MPrev[j] = (logf(scoreSum) - g_scoreCenter) +\r
- GapOpenB[0] + GapCloseB[j-1];\r
-#endif\r
- }\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
-#if OCC\r
- const FCOUNT OccAi = OccA[i];\r
-#endif\r
- for (unsigned j = 0; j < uLengthB; ++j)\r
- {\r
- if (MCurr[j] == 0)\r
- MCurr[j] = -2.5;\r
- else\r
-#if OCC\r
- MCurr[j] = (logf(MCurr[j]) - g_scoreCenter)*OccAi*OccB[j];\r
-#else\r
- MCurr[j] = (logf(MCurr[j]) - g_scoreCenter);\r
-#endif\r
- }\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 = 0;\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