25 unsigned **SortOrderA;
\r
26 unsigned *uDeletePos;
\r
31 static struct DP_MEMORY DPM;
\r
33 static void AllocDPMem(unsigned uLengthA, unsigned uLengthB)
\r
35 // Max prefix length
\r
36 unsigned uLength = (uLengthA > uLengthB ? uLengthA : uLengthB) + 1;
\r
37 if (uLength < DPM.uLength)
\r
40 // Add 256 to allow for future expansion and
\r
41 // round up to next multiple of 32.
\r
43 uLength += 32 - uLength%32;
\r
45 const unsigned uOldLength = DPM.uLength;
\r
48 for (unsigned i = 0; i < uOldLength; ++i)
\r
50 delete[] DPM.TraceBack[i];
\r
51 delete[] DPM.FreqsA[i];
\r
52 delete[] DPM.SortOrderA[i];
\r
54 for (unsigned n = 0; n < 20; ++n)
\r
55 delete[] DPM.ScoreMxB[n];
\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
78 DPM.uLength = uLength;
\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
85 DPM.OccA = new FCOUNT[uLength];
\r
86 DPM.OccB = new FCOUNT[uLength];
\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
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
101 DPM.TraceBack = new int*[uLength];
\r
103 for (unsigned uLetter = 0; uLetter < 20; ++uLetter)
\r
104 DPM.ScoreMxB[uLetter] = new SCORE[uLength];
\r
106 for (unsigned i = 0; i < uLength; ++i)
\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
114 SCORE GlobalAlignLE(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
115 unsigned uLengthB, PWPath &Path)
\r
117 SetTermGaps(PA, uLengthA);
\r
118 SetTermGaps(PB, uLengthB);
\r
120 const unsigned uPrefixCountA = uLengthA + 1;
\r
121 const unsigned uPrefixCountB = uLengthB + 1;
\r
123 AllocDPMem(uLengthA, uLengthB);
\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
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
137 SCORE *DPrev = DPM.DPrev;
\r
138 SCORE *DCurr = DPM.DCurr;
\r
139 SCORE *DWork = DPM.DWork;
\r
142 FCOUNT *OccA = DPM.OccA;
\r
143 FCOUNT *OccB = DPM.OccB;
\r
146 unsigned *uDeletePos = DPM.uDeletePos;
\r
148 int **TraceBack = DPM.TraceBack;
\r
150 for (unsigned i = 0; i < uLengthA; ++i)
\r
152 GapOpenA[i] = PA[i].m_scoreGapOpen;
\r
153 GapCloseA[i] = PA[i].m_scoreGapClose;
\r
155 OccA[i] = PA[i].m_fOcc;
\r
158 for (unsigned uLetter = 0; uLetter < 20; ++uLetter)
\r
160 SortOrderA[i][uLetter] = PA[i].m_uSortOrder[uLetter];
\r
161 FreqsA[i][uLetter] = PA[i].m_fcCounts[uLetter];
\r
165 for (unsigned j = 0; j < uLengthB; ++j)
\r
167 GapOpenB[j] = PB[j].m_scoreGapOpen;
\r
168 GapCloseB[j] = PB[j].m_scoreGapClose;
\r
170 OccB[j] = PB[j].m_fOcc;
\r
174 for (unsigned uLetter = 0; uLetter < 20; ++uLetter)
\r
176 for (unsigned j = 0; j < uLengthB; ++j)
\r
177 ScoreMxB[uLetter][j] = PB[j].m_AAScores[uLetter];
\r
180 for (unsigned i = 0; i < uPrefixCountA; ++i)
\r
181 memset(TraceBack[i], 0, uPrefixCountB*sizeof(int));
\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
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
196 const unsigned uLetter = *ptrSortOrderAi;
\r
197 const FCOUNT fcLetter = ptrFreqsAi[uLetter];
\r
200 scoreSum += fcLetter*ScoreMxB[uLetter][0];
\r
207 MPrev[0] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[0];
\r
209 MPrev[0] = (logf(scoreSum) - g_scoreCenter);
\r
213 // D(0,0) is -infinity (requires I->D).
\r
214 DPrev[0] = MINUS_INFINITY;
\r
216 for (unsigned j = 1; j < uLengthB; ++j)
\r
218 // Only way to get M(0, j) looks like this:
\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
229 const unsigned uLetter = *ptrSortOrderAi;
\r
230 const FCOUNT fcLetter = ptrFreqsAi[uLetter];
\r
233 scoreSum += fcLetter*ScoreMxB[uLetter][j];
\r
240 MPrev[j] = (logf(scoreSum) - g_scoreCenter)*OccA[0]*OccB[j] +
\r
241 GapOpenB[0] + GapCloseB[j-1];
\r
243 MPrev[j] = (logf(scoreSum) - g_scoreCenter) +
\r
244 GapOpenB[0] + GapCloseB[j-1];
\r
247 TraceBack[0][j] = -(int) j;
\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
255 for (unsigned i = 1; i < uLengthA; ++i)
\r
259 assert(ptrSortOrderA == &(SortOrderA[i]));
\r
260 assert(ptrFreqsA == &(FreqsA[i]));
\r
262 SCORE *ptrMCurr_j = MCurr;
\r
263 memset(ptrMCurr_j, 0, uLengthB*sizeof(SCORE));
\r
264 const FCOUNT *FreqsAi = *ptrFreqsA;
\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
273 const unsigned uLetter = *ptrSortOrderAi;
\r
274 SCORE *NSBR_Letter = ScoreMxB[uLetter];
\r
275 const FCOUNT fcLetter = FreqsAi[uLetter];
\r
278 SCORE *ptrNSBR = NSBR_Letter;
\r
279 for (SCORE *ptrMCurr = MCurr; ptrMCurr != ptrMCurrMax; ++ptrMCurr)
\r
280 *ptrMCurr += fcLetter*(*ptrNSBR++);
\r
284 const FCOUNT OccAi = OccA[i];
\r
286 for (unsigned j = 0; j < uLengthB; ++j)
\r
292 MCurr[j] = (logf(MCurr[j]) - g_scoreCenter)*OccAi*OccB[j];
\r
294 MCurr[j] = (logf(MCurr[j]) - g_scoreCenter);
\r
298 ptrMCurr_j = MCurr;
\r
299 unsigned *ptrDeletePos = uDeletePos;
\r
301 // Special case for j=0
\r
302 // Only way to get M(i, 0) looks like this:
\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
312 int *ptrTraceBack_ij = TraceBack[i];
\r
313 *ptrTraceBack_ij++ = (int) i;
\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
325 SCORE *ptrDCurr = DCurr;
\r
327 assert(ptrDCurr == &(DCurr[0]));
\r
330 // Can't have an insert if no letters from B
\r
331 IPrev_j_1 = MINUS_INFINITY;
\r
333 unsigned uInsertPos = 0;
\r
334 const SCORE scoreGapOpenAi = GapOpenA[i];
\r
335 const SCORE scoreGapCloseAi_1 = GapCloseA[i-1];
\r
337 for (unsigned j = 1; j < uLengthB; ++j)
\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
349 SCORE scoreMax = MPrev_j;
\r
351 assert(ptrDPrev == &(DPrev[j-1]));
\r
352 SCORE scoreD = *ptrDPrev++ + scoreGapCloseAi_1;
\r
353 if (scoreD > scoreMax)
\r
356 assert(ptrDeletePos == &(uDeletePos[j-1]));
\r
357 *ptrTraceBack_ij = (int) i - (int) *ptrDeletePos;
\r
358 assert(*ptrTraceBack_ij > 0);
\r
362 SCORE scoreI = IPrev_j_1 + GapCloseB[j-1];
\r
363 if (scoreI > scoreMax)
\r
366 *ptrTraceBack_ij = (int) uInsertPos - (int) j;
\r
367 assert(*ptrTraceBack_ij < 0);
\r
370 assert(ptrSortOrderA == &(SortOrderA[i]));
\r
371 assert(ptrFreqsA == &(FreqsA[i]));
\r
373 *ptrMCurr_j += scoreMax;
\r
374 assert(ptrMCurr_j == &(MCurr[j]));
\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
384 assert(ptrDeletePos == &uDeletePos[j]);
\r
387 assert(ptrDCurr + 1 == &(DCurr[j]));
\r
393 Rotate(MPrev, MCurr, MWork);
\r
394 Rotate(DPrev, DCurr, DWork);
\r
397 // Special case for i=uLengthA
\r
398 SCORE IPrev = MINUS_INFINITY;
\r
400 unsigned uInsertPos;
\r
402 for (unsigned j = 1; j < uLengthB; ++j)
\r
404 SCORE INew = MPrev[j-1] + GapOpenB[j];
\r
412 // Special case for i=uLengthA, j=uLengthB
\r
413 SCORE scoreMax = MPrev[uLengthB-1];
\r
414 int iTraceBack = 0;
\r
416 SCORE scoreD = DPrev[uLengthB-1] + GapCloseA[uLengthA-1];
\r
417 if (scoreD > scoreMax)
\r
420 iTraceBack = (int) uLengthA - (int) uDeletePos[uLengthB-1];
\r
423 SCORE scoreI = IPrev + GapCloseB[uLengthB-1];
\r
424 if (scoreI > scoreMax)
\r
427 iTraceBack = (int) uInsertPos - (int) uLengthB;
\r
430 TraceBack[uLengthA][uLengthB] = iTraceBack;
\r
432 TraceBackToPath(TraceBack, uLengthA, uLengthB, Path);
\r