8 static Mx<float> g_MxDPM;
9 static Mx<float> g_MxDPD;
10 static Mx<float> g_MxDPI;
12 static Mx<char> g_MxTBM;
13 static Mx<char> g_MxTBD;
14 static Mx<char> g_MxTBI;
25 static Mx<float> *g_DPMSimpleMx;
26 static Mx<float> *g_DPDSimpleMx;
27 static Mx<float> *g_DPISimpleMx;
28 static float **g_DPMSimple;
29 static float **g_DPDSimple;
30 static float **g_DPISimple;
32 #define cmpm(i, j, x) { if (!feq(x, g_DPMSimple[i][j])) \
34 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
35 __FILE__, __LINE__, x, i, j, g_DPMSimple[i][j]); \
39 #define cmpd(i, j, x) { if (!feq(x, g_DPDSimple[i][j])) \
41 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
42 __FILE__, __LINE__, x, i, j, g_DPDSimple[i][j]); \
46 #define cmpi(i, j, x) { if (!feq(x, g_DPISimple[i][j])) \
48 Die("%s:%d %.1f != DPMSimple[%u][%u] = %.1f", \
49 __FILE__, __LINE__, x, i, j, g_DPISimple[i][j]); \
55 #define cmpm(i, j, x) /* empty */
56 #define cmpd(i, j, x) /* empty */
57 #define cmpi(i, j, x) /* empty */
61 static void AllocSave(unsigned LA, unsigned LB)
64 GetSimpleDPMxs(&g_DPMSimpleMx, &g_DPDSimpleMx, &g_DPISimpleMx);
65 g_DPMSimple = g_DPMSimpleMx->GetData();
66 g_DPDSimple = g_DPDSimpleMx->GetData();
67 g_DPISimple = g_DPISimpleMx->GetData();
69 g_MxDPM.Alloc("FastM", LA+1, LB+1);
\r
70 g_MxDPD.Alloc("FastD", LA+1, LB+1);
\r
71 g_MxDPI.Alloc("FastI", LA+1, LB+1);
\r
73 g_MxTBM.Alloc("FastTBM", LA+1, LB+1);
\r
74 g_MxTBD.Alloc("FastTBD", LA+1, LB+1);
\r
75 g_MxTBI.Alloc("FastTBI", LA+1, LB+1);
\r
77 g_DPM = g_MxDPM.GetData();
\r
78 g_DPD = g_MxDPD.GetData();
\r
79 g_DPI = g_MxDPI.GetData();
\r
81 g_TBM = g_MxTBM.GetData();
\r
82 g_TBD = g_MxTBD.GetData();
\r
83 g_TBI = g_MxTBI.GetData();
\r
86 static void SAVE_DPM(unsigned i, unsigned j, float x)
91 asserta(feq(x, g_DPMSimple[i][j]));
95 static void SAVE_DPD(unsigned i, unsigned j, float x)
100 asserta(feq(x, g_DPDSimple[i][j]));
104 static void SAVE_DPI(unsigned i, unsigned j, float x)
109 asserta(feq(x, g_DPISimple[i][j]));
113 static void SAVE_TBM(unsigned i, unsigned j, char x)
118 static void SAVE_TBD(unsigned i, unsigned j, char x)
123 static void SAVE_TBI(unsigned i, unsigned j, char x)
128 void GetFastMxs(Mx<float> **M, Mx<float> **D, Mx<float> **I)
137 #define SAVE_DPM(i, j, x) /* empty */
138 #define SAVE_DPD(i, j, x) /* empty */
139 #define SAVE_DPI(i, j, x) /* empty */
141 #define SAVE_TBM(i, j, x) /* empty */
142 #define SAVE_TBD(i, j, x) /* empty */
143 #define SAVE_TBI(i, j, x) /* empty */
145 #define AllocSave(LA, LB) /* empty */
147 #define cmpm(i, j, x) /* empty */
148 #define cmpd(i, j, x) /* empty */
149 #define cmpi(i, j, x) /* empty */
153 float ViterbiFast(const byte *A, unsigned LA, const byte *B, unsigned LB,
154 const AlnParams &AP, PathData &PD)
156 if (LA*LB > 100*1000*1000)
157 Die("ViterbiFast, too long LA=%u, LB=%u", LA, LB);
162 StartTimer(ViterbiFast);
164 const float * const *Mx = AP.SubstMx;
165 float OpenA = AP.LOpenA;
166 float ExtA = AP.LExtA;
169 float *Mrow = g_DPRow1;
170 float *Drow = g_DPRow2;
172 // Use Mrow[-1], so...
173 Mrow[-1] = MINUS_INFINITY;
174 for (unsigned j = 0; j <= LB; ++j)
176 Mrow[j] = MINUS_INFINITY;
177 SAVE_DPM(0, j, MINUS_INFINITY);
180 Drow[j] = MINUS_INFINITY;
181 SAVE_DPD(0, j, MINUS_INFINITY);
186 float M0 = float (0);
188 for (unsigned i = 0; i < LA; ++i)
191 const float *MxRow = Mx[a];
192 float OpenB = AP.LOpenB;
193 float ExtB = AP.LExtB;
194 float I0 = MINUS_INFINITY;
198 SAVE_DPI(i, 0, MINUS_INFINITY);
199 SAVE_DPI(i, 1, MINUS_INFINITY);
205 for (unsigned j = 0; j < LB; ++j)
215 // Drow[j] = DPD[i][j]
221 SAVE_TBM(i+1, j+1, 'M');
225 TraceBits = TRACEBITS_DM;
226 SAVE_TBM(i+1, j+1, 'D');
231 TraceBits = TRACEBITS_IM;
232 SAVE_TBM(i+1, j+1, 'I');
237 Mrow[j] = xM + MxRow[b];
238 // Mrow[j] = DPM[i+1][j+1])
239 SAVE_DPM(i+1, j+1, Mrow[j]);
244 // SavedM0 = DPM[i][j]
245 // Drow[j] = DPD[i][j]
249 float md = SavedM0 + OpenB;
251 SAVE_TBD(i+1, j, 'D');
255 TraceBits |= TRACEBITS_MD;
256 SAVE_TBD(i+1, j, 'M');
258 // Drow[j] = DPD[i+1][j]
259 SAVE_DPD(i+1, j, Drow[j]);
264 // SavedM0 = DPM[i][j]
269 float mi = SavedM0 + OpenA;
271 SAVE_TBI(i, j+1, 'I');
275 TraceBits |= TRACEBITS_MI;
276 SAVE_TBI(i, j+1, 'M');
279 SAVE_DPI(i, j+1, I0);
285 TBrow[j] = TraceBits;
288 // Special case for end of Drow[]
291 // Drow[LB] = DPD[i][LB]
294 float md = M0 + AP.ROpenB;
295 Drow[LB] += AP.RExtB;
296 SAVE_TBD(i+1, LB, 'D');
300 TBrow[LB] = TRACEBITS_MD;
301 SAVE_TBD(i+1, LB, 'M');
303 // Drow[LB] = DPD[i+1][LB]
304 SAVE_DPD(i+1, LB, Drow[LB]);
307 SAVE_DPM(i+1, 0, MINUS_INFINITY);
314 SAVE_TBM(LA, 0, '?');
316 // Special case for last row of DPI
317 byte *TBrow = TB[LA];
318 float I1 = MINUS_INFINITY;
320 SAVE_DPI(LA, 0, MINUS_INFINITY);
321 SAVE_TBI(LA, 0, '?');
323 SAVE_DPI(LA, 1, MINUS_INFINITY);
324 SAVE_TBI(LA, 1, '?');
326 for (unsigned j = 1; j < LB; ++j)
328 // Mrow[j-1] = DPM[LA][j]
332 float mi = Mrow[int(j)-1] + AP.ROpenA;
334 SAVE_TBI(LA, j+1, 'I');
338 TBrow[j] = TRACEBITS_MI;
339 SAVE_TBI(LA, j+1, 'M');
341 SAVE_DPI(LA, j+1, I1);
344 float FinalM = Mrow[LB-1];
345 float FinalD = Drow[LB];
347 // FinalM = DPM[LA][LB]
348 // FinalD = DPD[LA][LB]
349 // FinalI = DPI[LA][LB]
351 float Score = FinalM;
364 EndTimer(ViterbiFast);
365 TraceBackBit(LA, LB, State, PD);