8 #define EQ(a, b) (fabs(a-b) < 0.1)
\r
10 void TraceBackSW(const ProfPos *PA, unsigned uLengthA, const ProfPos *PB,
\r
11 unsigned uLengthB, const SCORE *DPM_, const SCORE *DPD_, const SCORE *DPI_,
\r
12 unsigned uPrefixLengthAMax, unsigned uPrefixLengthBMax, PWPath &Path)
\r
16 Log("TraceBackSW LengthA=%u LengthB=%u PLAMax=%u PLBMax=%u\n",
\r
17 uLengthA, uLengthB, uPrefixLengthAMax, uPrefixLengthBMax);
\r
19 assert(uLengthB > 0 && uLengthA > 0);
\r
21 const unsigned uPrefixCountA = uLengthA + 1;
\r
22 const unsigned uPrefixCountB = uLengthB + 1;
\r
26 unsigned uPrefixLengthA = uPrefixLengthAMax;
\r
27 unsigned uPrefixLengthB = uPrefixLengthBMax;
\r
29 SCORE scoreMax = DPM(uPrefixLengthA, uPrefixLengthB);
\r
30 char cEdgeType = 'M';
\r
34 if ('S' == cEdgeType)
\r
38 Edge.cType = cEdgeType;
\r
39 Edge.uPrefixLengthA = uPrefixLengthA;
\r
40 Edge.uPrefixLengthB = uPrefixLengthB;
\r
41 Path.PrependEdge(Edge);
\r
44 unsigned uPrevPrefixLengthA = uPrefixLengthA;
\r
45 unsigned uPrevPrefixLengthB = uPrefixLengthB;
\r
51 assert(uPrefixLengthA > 0);
\r
52 assert(uPrefixLengthB > 0);
\r
53 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
54 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
56 const SCORE Score = DPM(uPrefixLengthA, uPrefixLengthB);
\r
57 const SCORE scoreMatch = ScoreProfPos2(PPA, PPB);
\r
60 if (1 == uPrefixLengthA && 1 == uPrefixLengthB)
\r
61 scoreSM = scoreMatch;
\r
63 scoreSM = MINUS_INFINITY;
\r
65 SCORE scoreMM = MINUS_INFINITY;
\r
66 SCORE scoreDM = MINUS_INFINITY;
\r
67 SCORE scoreIM = MINUS_INFINITY;
\r
68 if (uPrefixLengthA > 1 && uPrefixLengthB > 1)
\r
70 SCORE scoreTrans = DPM(uPrefixLengthA-1, uPrefixLengthB-1);
\r
71 scoreMM = scoreTrans + scoreMatch;
\r
73 if (uPrefixLengthA > 1)
\r
75 SCORE scoreTransDM = PA[uPrefixLengthA-2].m_scoreGapClose;
\r
76 scoreDM = DPD(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransDM + scoreMatch;
\r
78 if (uPrefixLengthB > 1)
\r
80 SCORE scoreTransIM = PB[uPrefixLengthB-2].m_scoreGapClose;
\r
81 scoreIM = DPI(uPrefixLengthA-1, uPrefixLengthB-1) + scoreTransIM + scoreMatch;
\r
84 if (EQ(scoreMM, Score))
\r
85 cPrevEdgeType = 'M';
\r
86 else if (EQ(scoreDM, Score))
\r
87 cPrevEdgeType = 'D';
\r
88 else if (EQ(scoreIM, Score))
\r
89 cPrevEdgeType = 'I';
\r
90 else if (EQ(scoreSM, Score))
\r
91 cPrevEdgeType = 'S';
\r
92 else if (EQ(scoreMatch, Score))
\r
93 cPrevEdgeType = 'S';
\r
95 Quit("TraceBack2: failed to match M score=%g M=%g D=%g I=%g S=%g",
\r
96 Score, scoreMM, scoreDM, scoreIM, scoreSM);
\r
98 --uPrevPrefixLengthA;
\r
99 --uPrevPrefixLengthB;
\r
105 assert(uPrefixLengthA > 0);
\r
106 const SCORE Score = DPD(uPrefixLengthA, uPrefixLengthB);
\r
108 SCORE scoreMD = MINUS_INFINITY;
\r
109 SCORE scoreDD = MINUS_INFINITY;
\r
110 SCORE scoreSD = MINUS_INFINITY;
\r
111 if (uPrefixLengthB == 0)
\r
113 if (uPrefixLengthA == 1)
\r
114 scoreSD = PA[0].m_scoreGapOpen;
\r
116 scoreSD = DPD(uPrefixLengthA - 1, 0);
\r
118 if (uPrefixLengthA > 1)
\r
120 const ProfPos &PPA = PA[uPrefixLengthA - 1];
\r
121 SCORE scoreTransMD = PPA.m_scoreGapOpen;
\r
122 scoreMD = DPM(uPrefixLengthA-1, uPrefixLengthB) + scoreTransMD;
\r
123 scoreDD = DPD(uPrefixLengthA-1, uPrefixLengthB);
\r
126 if (EQ(Score, scoreMD))
\r
127 cPrevEdgeType = 'M';
\r
128 else if (EQ(Score, scoreDD))
\r
129 cPrevEdgeType = 'D';
\r
130 else if (EQ(Score, scoreSD))
\r
131 cPrevEdgeType = 'S';
\r
133 Quit("TraceBack2: failed to match D");
\r
135 --uPrevPrefixLengthA;
\r
141 assert(uPrefixLengthB > 0);
\r
142 const SCORE Score = DPI(uPrefixLengthA, uPrefixLengthB);
\r
144 SCORE scoreMI = MINUS_INFINITY;
\r
145 SCORE scoreII = MINUS_INFINITY;
\r
146 SCORE scoreSI = MINUS_INFINITY;
\r
147 if (uPrefixLengthA == 0)
\r
149 if (uPrefixLengthB == 1)
\r
150 scoreSI = PB[0].m_scoreGapOpen;
\r
152 scoreSI = DPI(0, uPrefixLengthB - 1);
\r
154 if (uPrefixLengthB > 1)
\r
156 const ProfPos &PPB = PB[uPrefixLengthB - 1];
\r
157 SCORE scoreTransMI = PPB.m_scoreGapOpen;
\r
158 scoreMI = DPM(uPrefixLengthA, uPrefixLengthB-1) + scoreTransMI;
\r
159 scoreII = DPI(uPrefixLengthA, uPrefixLengthB-1);
\r
162 if (EQ(Score, scoreMI))
\r
163 cPrevEdgeType = 'M';
\r
164 else if (EQ(Score, scoreII))
\r
165 cPrevEdgeType = 'I';
\r
166 else if (EQ(Score, scoreSI))
\r
167 cPrevEdgeType = 'S';
\r
169 Quit("TraceBack2: failed to match I");
\r
171 --uPrevPrefixLengthB;
\r
179 Log("Edge %c%c%u.%u", cPrevEdgeType, cEdgeType, uPrefixLengthA, uPrefixLengthB);
\r
182 cEdgeType = cPrevEdgeType;
\r
183 uPrefixLengthA = uPrevPrefixLengthA;
\r
184 uPrefixLengthB = uPrevPrefixLengthB;
\r