8 static void LogPP(const ProfPos &PP)
\r
10 Log("ResidueGroup %u\n", PP.m_uResidueGroup);
\r
11 Log("AllGaps %d\n", PP.m_bAllGaps);
\r
12 Log("Occ %.3g\n", PP.m_fOcc);
\r
13 Log("LL=%.3g LG=%.3g GL=%.3g GG=%.3g\n", PP.m_LL, PP.m_LG, PP.m_GL, PP.m_GG);
\r
15 for (unsigned i = 0; i < 20; ++i)
\r
16 if (PP.m_fcCounts[i] > 0)
\r
17 Log("%c=%.3g ", LetterToChar(i), PP.m_fcCounts[i]);
\r
21 static void AssertProfPosEq(const ProfPos *PA, const ProfPos *PB, unsigned i)
\r
23 const ProfPos &PPA = PA[i];
\r
24 const ProfPos &PPB = PB[i];
\r
25 #define eq(x) if (PPA.m_##x != PPB.m_##x) { LogPP(PPA); LogPP(PPB); Quit("AssertProfPosEq." #x); }
\r
26 #define be(x) if (!BTEq(PPA.m_##x, PPB.m_##x)) { LogPP(PPA); LogPP(PPB); Quit("AssertProfPosEq." #x); }
\r
38 for (unsigned j = 0; j < 20; ++j)
\r
40 #define eqj(x) if (PPA.m_##x != PPB.m_##x) Quit("AssertProfPosEq j=%u " #x, j);
\r
41 #define bej(x) if (!BTEq(PPA.m_##x, PPB.m_##x)) Quit("AssertProfPosEq j=%u " #x, j);
\r
43 // eqj(uSortOrder[j]) // may differ due to ties, don't check?
\r
52 void AssertProfsEq(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
55 if (uLengthA != uLengthB)
\r
56 Quit("AssertProfsEq: lengths differ %u %u", uLengthA, uLengthB);
\r
57 for (unsigned i = 0; i < uLengthB; ++i)
\r
58 AssertProfPosEq(PA, PB, i);
\r
62 static void ValidateProf(const ProfPos *Prof, unsigned uLength)
\r
64 for (unsigned i = 0; i < uLength; ++i)
\r
66 const ProfPos &PP = Prof[i];
\r
68 FCOUNT s1 = PP.m_LL + PP.m_LG + PP.m_GL + PP.m_GG;
\r
69 assert(BTEq(s1, 1.0));
\r
73 const ProfPos &PPPrev = Prof[i-1];
\r
74 FCOUNT s2 = PPPrev.m_LL + PPPrev.m_GL;
\r
75 FCOUNT s3 = PP.m_LL + PP.m_LG;
\r
76 assert(BTEq(s2, s3));
\r
78 if (i < uLength - 1)
\r
80 const ProfPos &PPNext = Prof[i+1];
\r
81 FCOUNT s4 = PP.m_LL + PP.m_GL;
\r
82 FCOUNT s5 = PPNext.m_LL + PPNext.m_LG;
\r
83 assert(BTEq(s4, s5));
\r
88 #define ValidateProf(Prof, Length) /* empty */
\r
91 static void ScoresFromFreqsPos(ProfPos *Prof, unsigned uLength, unsigned uPos)
\r
93 ProfPos &PP = Prof[uPos];
\r
94 SortCounts(PP.m_fcCounts, PP.m_uSortOrder);
\r
95 PP.m_uResidueGroup = ResidueGroupFromFCounts(PP.m_fcCounts);
\r
98 PP.m_fOcc = PP.m_LL + PP.m_GL;
\r
100 // Frequency of gap-opens in this position (i)
\r
101 // Gap open = letter in i-1 and gap in i
\r
103 FCOUNT fcOpen = PP.m_LG;
\r
105 // Frequency of gap-closes in this position
\r
106 // Gap close = gap in i and letter in i+1
\r
109 if (uPos + 1 < uLength)
\r
110 fcClose = Prof[uPos + 1].m_GL;
\r
112 fcClose = PP.m_GG + PP.m_LG;
\r
114 PP.m_scoreGapOpen = (SCORE) ((1.0 - fcOpen)*g_scoreGapOpen/2.0);
\r
115 PP.m_scoreGapClose = (SCORE) ((1.0 - fcClose)*g_scoreGapOpen/2.0);
\r
117 PP.m_scoreGapOpen2 = (SCORE) ((1.0 - fcOpen)*g_scoreGapOpen2/2.0);
\r
118 PP.m_scoreGapClose2 = (SCORE) ((1.0 - fcClose)*g_scoreGapOpen2/2.0);
\r
121 for (unsigned i = 0; i < g_AlphaSize; ++i)
\r
123 SCORE scoreSum = 0;
\r
124 for (unsigned j = 0; j < g_AlphaSize; ++j)
\r
125 scoreSum += PP.m_fcCounts[j]*(*g_ptrScoreMatrix)[i][j];
\r
126 PP.m_AAScores[i] = scoreSum;
\r
130 void ProfScoresFromFreqs(ProfPos *Prof, unsigned uLength)
\r
132 for (unsigned i = 0; i < uLength; ++i)
\r
133 ScoresFromFreqsPos(Prof, uLength, i);
\r
136 static void AppendDelete(const MSA &msaA, unsigned &uColIndexA,
\r
137 unsigned uSeqCountA, unsigned uSeqCountB, MSA &msaCombined,
\r
138 unsigned &uColIndexCombined)
\r
141 Log("AppendDelete ColIxA=%u ColIxCmb=%u\n",
\r
142 uColIndexA, uColIndexCombined);
\r
144 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
146 char c = msaA.GetChar(uSeqIndexA, uColIndexA);
\r
147 msaCombined.SetChar(uSeqIndexA, uColIndexCombined, c);
\r
150 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
151 msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, '-');
\r
153 ++uColIndexCombined;
\r
157 static void AppendInsert(const MSA &msaB, unsigned &uColIndexB,
\r
158 unsigned uSeqCountA, unsigned uSeqCountB, MSA &msaCombined,
\r
159 unsigned &uColIndexCombined)
\r
162 Log("AppendInsert ColIxB=%u ColIxCmb=%u\n",
\r
163 uColIndexB, uColIndexCombined);
\r
165 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
166 msaCombined.SetChar(uSeqIndexA, uColIndexCombined, '-');
\r
168 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
170 char c = msaB.GetChar(uSeqIndexB, uColIndexB);
\r
171 msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, c);
\r
174 ++uColIndexCombined;
\r
178 static void AppendTplInserts(const MSA &msaA, unsigned &uColIndexA, unsigned uColCountA,
\r
179 const MSA &msaB, unsigned &uColIndexB, unsigned uColCountB, unsigned uSeqCountA,
\r
180 unsigned uSeqCountB, MSA &msaCombined, unsigned &uColIndexCombined)
\r
183 Log("AppendTplInserts ColIxA=%u ColIxB=%u ColIxCmb=%u\n",
\r
184 uColIndexA, uColIndexB, uColIndexCombined);
\r
186 const unsigned uLengthA = msaA.GetColCount();
\r
187 const unsigned uLengthB = msaB.GetColCount();
\r
189 unsigned uNewColCount = uColCountA;
\r
190 if (uColCountB > uNewColCount)
\r
191 uNewColCount = uColCountB;
\r
193 for (unsigned n = 0; n < uColCountA; ++n)
\r
195 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
197 char c = msaA.GetChar(uSeqIndexA, uColIndexA + n);
\r
198 c = UnalignChar(c);
\r
199 msaCombined.SetChar(uSeqIndexA, uColIndexCombined + n, c);
\r
202 for (unsigned n = uColCountA; n < uNewColCount; ++n)
\r
204 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
205 msaCombined.SetChar(uSeqIndexA, uColIndexCombined + n, '.');
\r
208 for (unsigned n = 0; n < uColCountB; ++n)
\r
210 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
212 char c = msaB.GetChar(uSeqIndexB, uColIndexB + n);
\r
213 c = UnalignChar(c);
\r
214 msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined + n, c);
\r
217 for (unsigned n = uColCountB; n < uNewColCount; ++n)
\r
219 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
220 msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined + n, '.');
\r
223 uColIndexCombined += uNewColCount;
\r
224 uColIndexA += uColCountA;
\r
225 uColIndexB += uColCountB;
\r
228 static void AppendMatch(const MSA &msaA, unsigned &uColIndexA, const MSA &msaB,
\r
229 unsigned &uColIndexB, unsigned uSeqCountA, unsigned uSeqCountB,
\r
230 MSA &msaCombined, unsigned &uColIndexCombined)
\r
233 Log("AppendMatch ColIxA=%u ColIxB=%u ColIxCmb=%u\n",
\r
234 uColIndexA, uColIndexB, uColIndexCombined);
\r
237 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
239 char c = msaA.GetChar(uSeqIndexA, uColIndexA);
\r
240 msaCombined.SetChar(uSeqIndexA, uColIndexCombined, c);
\r
243 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
245 char c = msaB.GetChar(uSeqIndexB, uColIndexB);
\r
246 msaCombined.SetChar(uSeqCountA + uSeqIndexB, uColIndexCombined, c);
\r
251 ++uColIndexCombined;
\r
254 void AlignTwoMSAsGivenPath(const PWPath &Path, const MSA &msaA, const MSA &msaB,
\r
257 msaCombined.Clear();
\r
260 Log("FastAlignProfiles\n");
\r
261 Log("Template A:\n");
\r
263 Log("Template B:\n");
\r
267 const unsigned uColCountA = msaA.GetColCount();
\r
268 const unsigned uColCountB = msaB.GetColCount();
\r
270 const unsigned uSeqCountA = msaA.GetSeqCount();
\r
271 const unsigned uSeqCountB = msaB.GetSeqCount();
\r
273 msaCombined.SetSeqCount(uSeqCountA + uSeqCountB);
\r
275 // Copy sequence names into combined MSA
\r
276 for (unsigned uSeqIndexA = 0; uSeqIndexA < uSeqCountA; ++uSeqIndexA)
\r
278 msaCombined.SetSeqName(uSeqIndexA, msaA.GetSeqName(uSeqIndexA));
\r
279 msaCombined.SetSeqId(uSeqIndexA, msaA.GetSeqId(uSeqIndexA));
\r
282 for (unsigned uSeqIndexB = 0; uSeqIndexB < uSeqCountB; ++uSeqIndexB)
\r
284 msaCombined.SetSeqName(uSeqCountA + uSeqIndexB, msaB.GetSeqName(uSeqIndexB));
\r
285 msaCombined.SetSeqId(uSeqCountA + uSeqIndexB, msaB.GetSeqId(uSeqIndexB));
\r
288 unsigned uColIndexA = 0;
\r
289 unsigned uColIndexB = 0;
\r
290 unsigned uColIndexCombined = 0;
\r
291 const unsigned uEdgeCount = Path.GetEdgeCount();
\r
292 for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
\r
294 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
\r
296 Log("\nEdge %u %c%u.%u\n",
\r
299 Edge.uPrefixLengthA,
\r
300 Edge.uPrefixLengthB);
\r
302 const char cType = Edge.cType;
\r
303 const unsigned uPrefixLengthA = Edge.uPrefixLengthA;
\r
304 unsigned uColCountA = 0;
\r
305 if (uPrefixLengthA > 0)
\r
307 const unsigned uNodeIndexA = uPrefixLengthA - 1;
\r
308 const unsigned uTplColIndexA = uNodeIndexA;
\r
309 if (uTplColIndexA > uColIndexA)
\r
310 uColCountA = uTplColIndexA - uColIndexA;
\r
313 const unsigned uPrefixLengthB = Edge.uPrefixLengthB;
\r
314 unsigned uColCountB = 0;
\r
315 if (uPrefixLengthB > 0)
\r
317 const unsigned uNodeIndexB = uPrefixLengthB - 1;
\r
318 const unsigned uTplColIndexB = uNodeIndexB;
\r
319 if (uTplColIndexB > uColIndexB)
\r
320 uColCountB = uTplColIndexB - uColIndexB;
\r
323 // TODO: This code looks like a hangover from HMM estimation -- can we delete it?
\r
324 assert(uColCountA == 0);
\r
325 assert(uColCountB == 0);
\r
326 AppendTplInserts(msaA, uColIndexA, uColCountA, msaB, uColIndexB, uColCountB,
\r
327 uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
\r
333 assert(uPrefixLengthA > 0);
\r
334 assert(uPrefixLengthB > 0);
\r
335 const unsigned uColA = uPrefixLengthA - 1;
\r
336 const unsigned uColB = uPrefixLengthB - 1;
\r
337 assert(uColIndexA == uColA);
\r
338 assert(uColIndexB == uColB);
\r
339 AppendMatch(msaA, uColIndexA, msaB, uColIndexB, uSeqCountA, uSeqCountB,
\r
340 msaCombined, uColIndexCombined);
\r
345 assert(uPrefixLengthA > 0);
\r
346 const unsigned uColA = uPrefixLengthA - 1;
\r
347 assert(uColIndexA == uColA);
\r
348 AppendDelete(msaA, uColIndexA, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
\r
353 assert(uPrefixLengthB > 0);
\r
354 const unsigned uColB = uPrefixLengthB - 1;
\r
355 assert(uColIndexB == uColB);
\r
356 AppendInsert(msaB, uColIndexB, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
\r
363 unsigned uInsertColCountA = uColCountA - uColIndexA;
\r
364 unsigned uInsertColCountB = uColCountB - uColIndexB;
\r
366 // TODO: This code looks like a hangover from HMM estimation -- can we delete it?
\r
367 assert(uInsertColCountA == 0);
\r
368 assert(uInsertColCountB == 0);
\r
369 AppendTplInserts(msaA, uColIndexA, uInsertColCountA, msaB, uColIndexB,
\r
370 uInsertColCountB, uSeqCountA, uSeqCountB, msaCombined, uColIndexCombined);
\r
372 assert(msaCombined.GetColCount() == uEdgeCount);
\r
375 static const ProfPos PPStart =
\r
377 false, //m_bAllGaps;
\r
378 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_uSortOrder[21];
\r
379 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_fcCounts[20];
\r
384 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // m_ALScores
\r
385 0, // m_uResidueGroup;
\r
387 0.0, // m_fcStartOcc;
\r
388 0.0, // m_fcEndOcc;
\r
389 0.0, // m_scoreGapOpen;
\r
390 0.0, // m_scoreGapClose;
\r
405 static void SetGapsMM(
\r
406 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
407 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
408 ProfPos *POut, unsigned uColIndexOut)
\r
410 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
411 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
412 ProfPos &PPO = POut[uColIndexOut];
\r
414 PPO.m_LL = wA*PPA.m_LL + wB*PPB.m_LL;
\r
415 PPO.m_LG = wA*PPA.m_LG + wB*PPB.m_LG;
\r
416 PPO.m_GL = wA*PPA.m_GL + wB*PPB.m_GL;
\r
417 PPO.m_GG = wA*PPA.m_GG + wB*PPB.m_GG;
\r
430 static void SetGapsMD(
\r
431 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
432 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
433 ProfPos *POut, unsigned uColIndexOut)
\r
435 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
436 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
437 ProfPos &PPO = POut[uColIndexOut];
\r
439 PPO.m_LL = wA*PPA.m_LL;
\r
440 PPO.m_LG = wA*PPA.m_LG + wB*(PPB.m_LL + PPB.m_GL);
\r
441 PPO.m_GL = wA*PPA.m_GL;
\r
442 PPO.m_GG = wA*PPA.m_GG + wB*(PPB.m_LG + PPB.m_GG);
\r
454 static void SetGapsDD(
\r
455 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
456 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
457 ProfPos *POut, unsigned uColIndexOut)
\r
459 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
460 ProfPos &PPO = POut[uColIndexOut];
\r
462 PPO.m_LL = wA*PPA.m_LL;
\r
463 PPO.m_LG = wA*PPA.m_LG;
\r
464 PPO.m_GL = wA*PPA.m_GL;
\r
465 PPO.m_GG = wA*PPA.m_GG + wB;
\r
478 static void SetGapsMI(
\r
479 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
480 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
481 ProfPos *POut, unsigned uColIndexOut)
\r
483 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
484 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
485 ProfPos &PPO = POut[uColIndexOut];
\r
487 PPO.m_LL = wB*PPB.m_LL;
\r
488 PPO.m_LG = wB*PPB.m_LG + wA*(PPA.m_LL + PPA.m_GL);
\r
489 PPO.m_GL = wB*PPB.m_GL;
\r
490 PPO.m_GG = wB*PPB.m_GG + wA*(PPA.m_LG + PPA.m_GG);
\r
503 static void SetGapsDM(
\r
504 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
505 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
506 ProfPos *POut, unsigned uColIndexOut)
\r
508 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
509 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
510 ProfPos &PPO = POut[uColIndexOut];
\r
512 PPO.m_LL = wA*PPA.m_LL;
\r
513 PPO.m_LG = wA*PPA.m_LG;
\r
514 PPO.m_GL = wA*PPA.m_GL + wB*(PPB.m_LL + PPB.m_GL);
\r
515 PPO.m_GG = wA*PPA.m_GG + wB*(PPB.m_LG + PPB.m_GG);
\r
528 static void SetGapsIM(
\r
529 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
530 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
531 ProfPos *POut, unsigned uColIndexOut)
\r
533 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
534 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
535 ProfPos &PPO = POut[uColIndexOut];
\r
537 PPO.m_LL = wB*PPB.m_LL;
\r
538 PPO.m_LG = wB*PPB.m_LG;
\r
539 PPO.m_GL = wB*PPB.m_GL + wA*(PPA.m_LL + PPA.m_GL);
\r
540 PPO.m_GG = wB*PPB.m_GG + wA*(PPA.m_LG + PPA.m_GG);
\r
551 static void SetGapsID(
\r
552 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
553 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
554 ProfPos *POut, unsigned uColIndexOut)
\r
556 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
557 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
558 ProfPos &PPO = POut[uColIndexOut];
\r
561 PPO.m_LG = wB*PPB.m_GL + wB*PPB.m_LL;
\r
562 PPO.m_GL = wA*PPA.m_GL + wA*PPA.m_LL;
\r
563 PPO.m_GG = wA*(PPA.m_LG + PPA.m_GG) + wB*(PPB.m_LG + PPB.m_GG);
\r
574 static void SetGapsDI(
\r
575 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
576 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
577 ProfPos *POut, unsigned uColIndexOut)
\r
579 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
580 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
581 ProfPos &PPO = POut[uColIndexOut];
\r
584 PPO.m_LG = wA*PPA.m_GL + wA*PPA.m_LL;
\r
585 PPO.m_GL = wB*PPB.m_GL + wB*PPB.m_LL;
\r
586 PPO.m_GG = wA*(PPA.m_LG + PPA.m_GG) + wB*(PPB.m_LG + PPB.m_GG);
\r
598 static void SetGapsII(
\r
599 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
600 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
601 ProfPos *POut, unsigned uColIndexOut)
\r
603 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
604 ProfPos &PPO = POut[uColIndexOut];
\r
606 PPO.m_LL = wB*PPB.m_LL;
\r
607 PPO.m_LG = wB*PPB.m_LG;
\r
608 PPO.m_GL = wB*PPB.m_GL;
\r
609 PPO.m_GG = wB*PPB.m_GG + wA;
\r
612 static void SetFreqs(
\r
613 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
614 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
615 ProfPos *POut, unsigned uColIndexOut)
\r
617 const ProfPos &PPA = uPrefixLengthA > 0 ? PA[uPrefixLengthA-1] : PPStart;
\r
618 const ProfPos &PPB = uPrefixLengthB > 0 ? PB[uPrefixLengthB-1] : PPStart;
\r
619 ProfPos &PPO = POut[uColIndexOut];
\r
621 if (g_bNormalizeCounts)
\r
623 const FCOUNT fA = PPA.m_fOcc*wA/(wA + wB);
\r
624 const FCOUNT fB = PPB.m_fOcc*wB/(wA + wB);
\r
626 for (unsigned i = 0; i < 20; ++i)
\r
628 const FCOUNT f = fA*PPA.m_fcCounts[i] + fB*PPB.m_fcCounts[i];
\r
629 PPO.m_fcCounts[i] = f;
\r
633 for (unsigned i = 0; i < 20; ++i)
\r
634 PPO.m_fcCounts[i] /= fTotal;
\r
638 for (unsigned i = 0; i < 20; ++i)
\r
639 PPO.m_fcCounts[i] = wA*PPA.m_fcCounts[i] + wB*PPB.m_fcCounts[i];
\r
643 void AlignTwoProfsGivenPath(const PWPath &Path,
\r
644 const ProfPos *PA, unsigned uPrefixLengthA, WEIGHT wA,
\r
645 const ProfPos *PB, unsigned uPrefixLengthB, WEIGHT wB,
\r
646 ProfPos **ptrPOut, unsigned *ptruLengthOut)
\r
649 Log("AlignTwoProfsGivenPath wA=%.3g wB=%.3g Path=\n", wA, wB);
\r
652 assert(BTEq(wA + wB, 1.0));
\r
654 unsigned uColIndexA = 0;
\r
655 unsigned uColIndexB = 0;
\r
656 unsigned uColIndexOut = 0;
\r
657 const unsigned uEdgeCount = Path.GetEdgeCount();
\r
658 ProfPos *POut = new ProfPos[uEdgeCount];
\r
659 char cPrevType = 'M';
\r
660 for (unsigned uEdgeIndex = 0; uEdgeIndex < uEdgeCount; ++uEdgeIndex)
\r
662 const PWEdge &Edge = Path.GetEdge(uEdgeIndex);
\r
663 const char cType = Edge.cType;
\r
665 const unsigned uPrefixLengthA = Edge.uPrefixLengthA;
\r
666 const unsigned uPrefixLengthB = Edge.uPrefixLengthB;
\r
669 Log("\nEdge %u %c%u.%u ColA=%u ColB=%u\n",
\r
672 Edge.uPrefixLengthA,
\r
673 Edge.uPrefixLengthB,
\r
678 POut[uColIndexOut].m_bAllGaps = false;
\r
683 assert(uPrefixLengthA > 0);
\r
684 assert(uPrefixLengthB > 0);
\r
686 PA, uPrefixLengthA, wA,
\r
687 PB, uPrefixLengthB, wB,
\r
688 POut, uColIndexOut);
\r
693 PA, uPrefixLengthA, wA,
\r
694 PB, uPrefixLengthB, wB,
\r
695 POut, uColIndexOut);
\r
699 PA, uPrefixLengthA, wA,
\r
700 PB, uPrefixLengthB, wB,
\r
701 POut, uColIndexOut);
\r
705 PA, uPrefixLengthA, wA,
\r
706 PB, uPrefixLengthB, wB,
\r
707 POut, uColIndexOut);
\r
710 Quit("Bad cPrevType");
\r
719 assert(uPrefixLengthA > 0);
\r
721 PA, uPrefixLengthA, wA,
\r
722 PB, uPrefixLengthB, 0,
\r
723 POut, uColIndexOut);
\r
728 PA, uPrefixLengthA, wA,
\r
729 PB, uPrefixLengthB, wB,
\r
730 POut, uColIndexOut);
\r
734 PA, uPrefixLengthA, wA,
\r
735 PB, uPrefixLengthB, wB,
\r
736 POut, uColIndexOut);
\r
740 PA, uPrefixLengthA, wA,
\r
741 PB, uPrefixLengthB, wB,
\r
742 POut, uColIndexOut);
\r
745 Quit("Bad cPrevType");
\r
753 assert(uPrefixLengthB > 0);
\r
755 PA, uPrefixLengthA, 0,
\r
756 PB, uPrefixLengthB, wB,
\r
757 POut, uColIndexOut);
\r
762 PA, uPrefixLengthA, wA,
\r
763 PB, uPrefixLengthB, wB,
\r
764 POut, uColIndexOut);
\r
768 PA, uPrefixLengthA, wA,
\r
769 PB, uPrefixLengthB, wB,
\r
770 POut, uColIndexOut);
\r
774 PA, uPrefixLengthA, wA,
\r
775 PB, uPrefixLengthB, wB,
\r
776 POut, uColIndexOut);
\r
779 Quit("Bad cPrevType");
\r
790 assert(uColIndexOut == uEdgeCount);
\r
792 ProfScoresFromFreqs(POut, uEdgeCount);
\r
793 ValidateProf(POut, uEdgeCount);
\r
796 *ptruLengthOut = uEdgeCount;
\r
799 Log("AlignTwoProfsGivenPath:\n");
\r
800 ListProfile(POut, uEdgeCount, 0);
\r