1 //uchime by Robert C. Edgar http://drive5.com/uchime This code is donated to the public domain.
10 static Mx<float> g_MxDPM;
11 static Mx<float> g_MxDPD;
12 static Mx<float> g_MxDPI;
14 static Mx<char> g_MxTBM;
15 static Mx<char> g_MxTBD;
16 static Mx<char> g_MxTBI;
27 static Mx<float> *g_DPMSimpleMx;
28 static Mx<float> *g_DPDSimpleMx;
29 static Mx<float> *g_DPISimpleMx;
30 static float **g_DPMSimple;
31 static float **g_DPDSimple;
32 static float **g_DPISimple;
34 #define cmpm(i, j, x) { if (!feq(x, g_DPMSimple[i][j])) \
36 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
37 __FILE__, __LINE__, x, i, j, g_DPMSimple[i][j]); \
41 #define cmpd(i, j, x) { if (!feq(x, g_DPDSimple[i][j])) \
43 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
44 __FILE__, __LINE__, x, i, j, g_DPDSimple[i][j]); \
48 #define cmpi(i, j, x) { if (!feq(x, g_DPISimple[i][j])) \
50 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
51 __FILE__, __LINE__, x, i, j, g_DPISimple[i][j]); \
57 #define cmpm(i, j, x) /* empty */
58 #define cmpd(i, j, x) /* empty */
59 #define cmpi(i, j, x) /* empty */
63 static void AllocSave(unsigned LA, unsigned LB)
66 GetSimpleDPMxs(&g_DPMSimpleMx, &g_DPDSimpleMx, &g_DPISimpleMx);
67 g_DPMSimple = g_DPMSimpleMx->GetData();
68 g_DPDSimple = g_DPDSimpleMx->GetData();
69 g_DPISimple = g_DPISimpleMx->GetData();
71 g_MxDPM.Alloc("FastM", LA+1, LB+1);
72 g_MxDPD.Alloc("FastD", LA+1, LB+1);
73 g_MxDPI.Alloc("FastI", LA+1, LB+1);
75 g_MxTBM.Alloc("FastTBM", LA+1, LB+1);
76 g_MxTBD.Alloc("FastTBD", LA+1, LB+1);
77 g_MxTBI.Alloc("FastTBI", LA+1, LB+1);
79 g_DPM = g_MxDPM.GetData();
80 g_DPD = g_MxDPD.GetData();
81 g_DPI = g_MxDPI.GetData();
83 g_TBM = g_MxTBM.GetData();
84 g_TBD = g_MxTBD.GetData();
85 g_TBI = g_MxTBI.GetData();
88 static void SAVE_DPM(unsigned i, unsigned j, float x)
93 asserta(feq(x, g_DPMSimple[i][j]));
97 static void SAVE_DPD(unsigned i, unsigned j, float x)
102 asserta(feq(x, g_DPDSimple[i][j]));
106 static void SAVE_DPI(unsigned i, unsigned j, float x)
111 asserta(feq(x, g_DPISimple[i][j]));
115 static void SAVE_TBM(unsigned i, unsigned j, char x)
120 static void SAVE_TBD(unsigned i, unsigned j, char x)
125 static void SAVE_TBI(unsigned i, unsigned j, char x)
130 void GetFastMxs(Mx<float> **M, Mx<float> **D, Mx<float> **I)
139 #define SAVE_DPM(i, j, x) /* empty */
140 #define SAVE_DPD(i, j, x) /* empty */
141 #define SAVE_DPI(i, j, x) /* empty */
143 #define SAVE_TBM(i, j, x) /* empty */
144 #define SAVE_TBD(i, j, x) /* empty */
145 #define SAVE_TBI(i, j, x) /* empty */
147 #define AllocSave(LA, LB) /* empty */
149 #define cmpm(i, j, x) /* empty */
150 #define cmpd(i, j, x) /* empty */
151 #define cmpi(i, j, x) /* empty */
155 float ViterbiFast(const byte *A, unsigned LA, const byte *B, unsigned LB,
156 const AlnParams &AP, PathData &PD)
158 if (LA*LB > 100*1000*1000)
159 Die("ViterbiFast, too long LA=%u, LB=%u", LA, LB);
164 StartTimer(ViterbiFast);
166 const float * const *Mx = AP.SubstMx;
167 float OpenA = AP.LOpenA;
168 float ExtA = AP.LExtA;
171 float *Mrow = g_DPRow1;
172 float *Drow = g_DPRow2;
174 // Use Mrow[-1], so...
175 Mrow[-1] = MINUS_INFINITY;
176 for (unsigned j = 0; j <= LB; ++j)
178 Mrow[j] = MINUS_INFINITY;
179 SAVE_DPM(0, j, MINUS_INFINITY);
182 Drow[j] = MINUS_INFINITY;
183 SAVE_DPD(0, j, MINUS_INFINITY);
188 float M0 = float (0);
190 for (unsigned i = 0; i < LA; ++i)
193 const float *MxRow = Mx[a];
194 float OpenB = AP.LOpenB;
195 float ExtB = AP.LExtB;
196 float I0 = MINUS_INFINITY;
200 SAVE_DPI(i, 0, MINUS_INFINITY);
201 SAVE_DPI(i, 1, MINUS_INFINITY);
207 for (unsigned j = 0; j < LB; ++j)
217 // Drow[j] = DPD[i][j]
223 SAVE_TBM(i+1, j+1, 'M');
227 TraceBits = TRACEBITS_DM;
228 SAVE_TBM(i+1, j+1, 'D');
233 TraceBits = TRACEBITS_IM;
234 SAVE_TBM(i+1, j+1, 'I');
239 Mrow[j] = xM + MxRow[b];
240 // Mrow[j] = DPM[i+1][j+1])
241 SAVE_DPM(i+1, j+1, Mrow[j]);
246 // SavedM0 = DPM[i][j]
247 // Drow[j] = DPD[i][j]
251 float md = SavedM0 + OpenB;
253 SAVE_TBD(i+1, j, 'D');
257 TraceBits |= TRACEBITS_MD;
258 SAVE_TBD(i+1, j, 'M');
260 // Drow[j] = DPD[i+1][j]
261 SAVE_DPD(i+1, j, Drow[j]);
266 // SavedM0 = DPM[i][j]
271 float mi = SavedM0 + OpenA;
273 SAVE_TBI(i, j+1, 'I');
277 TraceBits |= TRACEBITS_MI;
278 SAVE_TBI(i, j+1, 'M');
281 SAVE_DPI(i, j+1, I0);
287 TBrow[j] = TraceBits;
290 // Special case for end of Drow[]
293 // Drow[LB] = DPD[i][LB]
296 float md = M0 + AP.ROpenB;
297 Drow[LB] += AP.RExtB;
298 SAVE_TBD(i+1, LB, 'D');
302 TBrow[LB] = TRACEBITS_MD;
303 SAVE_TBD(i+1, LB, 'M');
305 // Drow[LB] = DPD[i+1][LB]
306 SAVE_DPD(i+1, LB, Drow[LB]);
309 SAVE_DPM(i+1, 0, MINUS_INFINITY);
316 SAVE_TBM(LA, 0, '?');
318 // Special case for last row of DPI
319 byte *TBrow = TB[LA];
320 float I1 = MINUS_INFINITY;
322 SAVE_DPI(LA, 0, MINUS_INFINITY);
323 SAVE_TBI(LA, 0, '?');
325 SAVE_DPI(LA, 1, MINUS_INFINITY);
326 SAVE_TBI(LA, 1, '?');
328 for (unsigned j = 1; j < LB; ++j)
330 // Mrow[j-1] = DPM[LA][j]
334 float mi = Mrow[int(j)-1] + AP.ROpenA;
336 SAVE_TBI(LA, j+1, 'I');
340 TBrow[j] = TRACEBITS_MI;
341 SAVE_TBI(LA, j+1, 'M');
343 SAVE_DPI(LA, j+1, I1);
346 float FinalM = Mrow[LB-1];
347 float FinalD = Drow[LB];
349 // FinalM = DPM[LA][LB]
350 // FinalD = DPD[LA][LB]
351 // FinalI = DPI[LA][LB]
353 float Score = FinalM;
366 EndTimer(ViterbiFast);
367 TraceBackBit(LA, LB, State, PD);