1 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
2 /* ====================================================================
3 * Copyright (c) 1999-2004 Carnegie Mellon University. All rights
4 * reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 *
13 * 2. Redistributions in binary form must reproduce the above copyright
14 * notice, this list of conditions and the following disclaimer in
15 * the documentation and/or other materials provided with the
16 * distribution.
17 *
18 * This work was supported in part by funding from the Defense Advanced
19 * Research Projects Agency and the National Science Foundation of the
20 * United States of America, and the CMU Sphinx Speech Consortium.
21 *
22 * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND
23 * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
24 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
25 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
26 * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
27 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
28 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
29 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
30 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
31 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
32 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 *
34 * ====================================================================
35 *
36 */
37
38 /**
39 * @file hmm.h Implementation of HMM base structure.
40 */
41
42 /* System headers. */
43 #include <assert.h>
44 #include <stdlib.h>
45 #include <string.h>
46 #include <limits.h>
47
48 /* SphinxBase headers. */
49 #include <sphinxbase/ckd_alloc.h>
50 #include <sphinxbase/err.h>
51
52 /* Local headers. */
53 #include "hmm.h"
54
55 hmm_context_t *
hmm_context_init(int32 n_emit_state,uint8 ** const * tp,int16 const * senscore,uint16 * const * sseq)56 hmm_context_init(int32 n_emit_state,
57 uint8 ** const *tp,
58 int16 const *senscore,
59 uint16 * const *sseq)
60 {
61 hmm_context_t *ctx;
62
63 assert(n_emit_state > 0);
64 if (n_emit_state > HMM_MAX_NSTATE) {
65 E_ERROR("Number of emitting states must be <= %d\n", HMM_MAX_NSTATE);
66 return NULL;
67 }
68
69 ctx = ckd_calloc(1, sizeof(*ctx));
70 ctx->n_emit_state = n_emit_state;
71 ctx->tp = tp;
72 ctx->senscore = senscore;
73 ctx->sseq = sseq;
74 ctx->st_sen_scr = ckd_calloc(n_emit_state, sizeof(*ctx->st_sen_scr));
75
76 return ctx;
77 }
78
79 void
hmm_context_free(hmm_context_t * ctx)80 hmm_context_free(hmm_context_t *ctx)
81 {
82 if (ctx == NULL)
83 return;
84 ckd_free(ctx->st_sen_scr);
85 ckd_free(ctx);
86 }
87
88 void
hmm_init(hmm_context_t * ctx,hmm_t * hmm,int mpx,int ssid,int tmatid)89 hmm_init(hmm_context_t *ctx, hmm_t *hmm, int mpx, int ssid, int tmatid)
90 {
91 hmm->ctx = ctx;
92 hmm->mpx = mpx;
93 hmm->n_emit_state = ctx->n_emit_state;
94 if (mpx) {
95 int i;
96 hmm->ssid = BAD_SSID;
97 hmm->senid[0] = ssid;
98 for (i = 1; i < hmm_n_emit_state(hmm); ++i) {
99 hmm->senid[i] = BAD_SSID;
100 }
101 }
102 else {
103 hmm->ssid = ssid;
104 memcpy(hmm->senid, ctx->sseq[ssid], hmm->n_emit_state * sizeof(*hmm->senid));
105 }
106 hmm->tmatid = tmatid;
107 hmm_clear(hmm);
108 }
109
110 void
hmm_deinit(hmm_t * hmm)111 hmm_deinit(hmm_t *hmm)
112 {
113 }
114
115 void
hmm_dump(hmm_t * hmm,FILE * fp)116 hmm_dump(hmm_t * hmm,
117 FILE * fp)
118 {
119 int32 i;
120
121 if (hmm_is_mpx(hmm)) {
122 fprintf(fp, "MPX ");
123 for (i = 0; i < hmm_n_emit_state(hmm); i++)
124 fprintf(fp, " %11d", hmm_senid(hmm, i));
125 fprintf(fp, " ( ");
126 for (i = 0; i < hmm_n_emit_state(hmm); i++)
127 fprintf(fp, "%d ", hmm_ssid(hmm, i));
128 fprintf(fp, ")\n");
129 }
130 else {
131 fprintf(fp, "SSID ");
132 for (i = 0; i < hmm_n_emit_state(hmm); i++)
133 fprintf(fp, " %11d", hmm_senid(hmm, i));
134 fprintf(fp, " (%d)\n", hmm_ssid(hmm, 0));
135 }
136
137 if (hmm->ctx->senscore) {
138 fprintf(fp, "SENSCR");
139 for (i = 0; i < hmm_n_emit_state(hmm); i++)
140 fprintf(fp, " %11d", hmm_senscr(hmm, i));
141 fprintf(fp, "\n");
142 }
143
144 fprintf(fp, "SCORES %11d", hmm_in_score(hmm));
145 for (i = 1; i < hmm_n_emit_state(hmm); i++)
146 fprintf(fp, " %11d", hmm_score(hmm, i));
147 fprintf(fp, " %11d", hmm_out_score(hmm));
148 fprintf(fp, "\n");
149
150 fprintf(fp, "HISTID %11d", hmm_in_history(hmm));
151 for (i = 1; i < hmm_n_emit_state(hmm); i++)
152 fprintf(fp, " %11d", hmm_history(hmm, i));
153 fprintf(fp, " %11d", hmm_out_history(hmm));
154 fprintf(fp, "\n");
155
156 if (hmm_in_score(hmm) > 0)
157 fprintf(fp,
158 "ALERT!! The input score %d is large than 0. Probably wrap around.\n",
159 hmm_in_score(hmm));
160 if (hmm_out_score(hmm) > 0)
161 fprintf(fp,
162 "ALERT!! The output score %d is large than 0. Probably wrap around\n.",
163 hmm_out_score(hmm));
164
165 fflush(fp);
166 }
167
168
169 void
hmm_clear_scores(hmm_t * h)170 hmm_clear_scores(hmm_t * h)
171 {
172 int32 i;
173
174 hmm_in_score(h) = WORST_SCORE;
175 for (i = 1; i < hmm_n_emit_state(h); i++)
176 hmm_score(h, i) = WORST_SCORE;
177 hmm_out_score(h) = WORST_SCORE;
178
179 h->bestscore = WORST_SCORE;
180 }
181
182 void
hmm_clear(hmm_t * h)183 hmm_clear(hmm_t * h)
184 {
185 int32 i;
186
187 hmm_in_score(h) = WORST_SCORE;
188 hmm_in_history(h) = -1;
189 for (i = 1; i < hmm_n_emit_state(h); i++) {
190 hmm_score(h, i) = WORST_SCORE;
191 hmm_history(h, i) = -1;
192 }
193 hmm_out_score(h) = WORST_SCORE;
194 hmm_out_history(h) = -1;
195
196 h->bestscore = WORST_SCORE;
197 h->frame = -1;
198 }
199
200 void
hmm_enter(hmm_t * h,int32 score,int32 histid,int frame)201 hmm_enter(hmm_t *h, int32 score, int32 histid, int frame)
202 {
203 hmm_in_score(h) = score;
204 hmm_in_history(h) = histid;
205 hmm_frame(h) = frame;
206 }
207
208 void
hmm_normalize(hmm_t * h,int32 bestscr)209 hmm_normalize(hmm_t *h, int32 bestscr)
210 {
211 int32 i;
212
213 for (i = 0; i < hmm_n_emit_state(h); i++) {
214 if (hmm_score(h, i) BETTER_THAN WORST_SCORE)
215 hmm_score(h, i) -= bestscr;
216 }
217 if (hmm_out_score(h) BETTER_THAN WORST_SCORE)
218 hmm_out_score(h) -= bestscr;
219 }
220
221 #define hmm_tprob_5st(i, j) (-tp[(i)*6+(j)])
222 #define nonmpx_senscr(i) (-senscore[sseq[i]])
223
224 static int32
hmm_vit_eval_5st_lr(hmm_t * hmm)225 hmm_vit_eval_5st_lr(hmm_t * hmm)
226 {
227 int16 const *senscore = hmm->ctx->senscore;
228 uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
229 uint16 const *sseq = hmm->senid;
230 int32 s5, s4, s3, s2, s1, s0, t2, t1, t0, bestScore;
231
232 /* It was the best of scores, it was the worst of scores. */
233 bestScore = WORST_SCORE;
234
235 /* Cache problem here! */
236 s4 = hmm_score(hmm, 4) + nonmpx_senscr(4);
237 s3 = hmm_score(hmm, 3) + nonmpx_senscr(3);
238 /* Transitions into non-emitting state 5 */
239 if (s3 BETTER_THAN WORST_SCORE) {
240 t1 = s4 + hmm_tprob_5st(4, 5);
241 t2 = s3 + hmm_tprob_5st(3, 5);
242 if (t1 BETTER_THAN t2) {
243 s5 = t1;
244 hmm_out_history(hmm) = hmm_history(hmm, 4);
245 } else {
246 s5 = t2;
247 hmm_out_history(hmm) = hmm_history(hmm, 3);
248 }
249 if (s5 WORSE_THAN WORST_SCORE) s5 = WORST_SCORE;
250 hmm_out_score(hmm) = s5;
251 bestScore = s5;
252 }
253
254 s2 = hmm_score(hmm, 2) + nonmpx_senscr(2);
255 /* All transitions into state 4 */
256 if (s2 BETTER_THAN WORST_SCORE) {
257 t0 = s4 + hmm_tprob_5st(4, 4);
258 t1 = s3 + hmm_tprob_5st(3, 4);
259 t2 = s2 + hmm_tprob_5st(2, 4);
260 if (t0 BETTER_THAN t1) {
261 if (t2 BETTER_THAN t0) {
262 s4 = t2;
263 hmm_history(hmm, 4) = hmm_history(hmm, 2);
264 } else
265 s4 = t0;
266 } else {
267 if (t2 BETTER_THAN t1) {
268 s4 = t2;
269 hmm_history(hmm, 4) = hmm_history(hmm, 2);
270 } else {
271 s4 = t1;
272 hmm_history(hmm, 4) = hmm_history(hmm, 3);
273 }
274 }
275 if (s4 WORSE_THAN WORST_SCORE) s4 = WORST_SCORE;
276 if (s4 BETTER_THAN bestScore) bestScore = s4;
277 hmm_score(hmm, 4) = s4;
278 }
279
280 s1 = hmm_score(hmm, 1) + nonmpx_senscr(1);
281 /* All transitions into state 3 */
282 if (s1 BETTER_THAN WORST_SCORE) {
283 t0 = s3 + hmm_tprob_5st(3, 3);
284 t1 = s2 + hmm_tprob_5st(2, 3);
285 t2 = s1 + hmm_tprob_5st(1, 3);
286 if (t0 BETTER_THAN t1) {
287 if (t2 BETTER_THAN t0) {
288 s3 = t2;
289 hmm_history(hmm, 3) = hmm_history(hmm, 1);
290 } else
291 s3 = t0;
292 } else {
293 if (t2 BETTER_THAN t1) {
294 s3 = t2;
295 hmm_history(hmm, 3) = hmm_history(hmm, 1);
296 } else {
297 s3 = t1;
298 hmm_history(hmm, 3) = hmm_history(hmm, 2);
299 }
300 }
301 if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
302 if (s3 BETTER_THAN bestScore) bestScore = s3;
303 hmm_score(hmm, 3) = s3;
304 }
305
306 s0 = hmm_in_score(hmm) + nonmpx_senscr(0);
307 /* All transitions into state 2 (state 0 is always active) */
308 t0 = s2 + hmm_tprob_5st(2, 2);
309 t1 = s1 + hmm_tprob_5st(1, 2);
310 t2 = s0 + hmm_tprob_5st(0, 2);
311 if (t0 BETTER_THAN t1) {
312 if (t2 BETTER_THAN t0) {
313 s2 = t2;
314 hmm_history(hmm, 2) = hmm_in_history(hmm);
315 } else
316 s2 = t0;
317 } else {
318 if (t2 BETTER_THAN t1) {
319 s2 = t2;
320 hmm_history(hmm, 2) = hmm_in_history(hmm);
321 } else {
322 s2 = t1;
323 hmm_history(hmm, 2) = hmm_history(hmm, 1);
324 }
325 }
326 if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
327 if (s2 BETTER_THAN bestScore) bestScore = s2;
328 hmm_score(hmm, 2) = s2;
329
330
331 /* All transitions into state 1 */
332 t0 = s1 + hmm_tprob_5st(1, 1);
333 t1 = s0 + hmm_tprob_5st(0, 1);
334 if (t0 BETTER_THAN t1) {
335 s1 = t0;
336 } else {
337 s1 = t1;
338 hmm_history(hmm, 1) = hmm_in_history(hmm);
339 }
340 if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
341 if (s1 BETTER_THAN bestScore) bestScore = s1;
342 hmm_score(hmm, 1) = s1;
343
344 /* All transitions into state 0 */
345 s0 = s0 + hmm_tprob_5st(0, 0);
346 if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
347 if (s0 BETTER_THAN bestScore) bestScore = s0;
348 hmm_in_score(hmm) = s0;
349
350 hmm_bestscore(hmm) = bestScore;
351 return bestScore;
352 }
353
354 #define mpx_senid(st) sseq[ssid[st]][st]
355 #define mpx_senscr(st) (-senscore[mpx_senid(st)])
356
357 static int32
hmm_vit_eval_5st_lr_mpx(hmm_t * hmm)358 hmm_vit_eval_5st_lr_mpx(hmm_t * hmm)
359 {
360 uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
361 int16 const *senscore = hmm->ctx->senscore;
362 uint16 * const *sseq = hmm->ctx->sseq;
363 uint16 *ssid = hmm->senid;
364 int32 bestScore;
365 int32 s5, s4, s3, s2, s1, s0, t2, t1, t0;
366
367 /* Don't propagate WORST_SCORE */
368 if (ssid[4] == BAD_SSID)
369 s4 = t1 = WORST_SCORE;
370 else {
371 s4 = hmm_score(hmm, 4) + mpx_senscr(4);
372 t1 = s4 + hmm_tprob_5st(4, 5);
373 }
374 if (ssid[3] == BAD_SSID)
375 s3 = t2 = WORST_SCORE;
376 else {
377 s3 = hmm_score(hmm, 3) + mpx_senscr(3);
378 t2 = s3 + hmm_tprob_5st(3, 5);
379 }
380 if (t1 BETTER_THAN t2) {
381 s5 = t1;
382 hmm_out_history(hmm) = hmm_history(hmm, 4);
383 }
384 else {
385 s5 = t2;
386 hmm_out_history(hmm) = hmm_history(hmm, 3);
387 }
388 if (s5 WORSE_THAN WORST_SCORE) s5 = WORST_SCORE;
389 hmm_out_score(hmm) = s5;
390 bestScore = s5;
391
392 /* Don't propagate WORST_SCORE */
393 if (ssid[2] == BAD_SSID)
394 s2 = t2 = WORST_SCORE;
395 else {
396 s2 = hmm_score(hmm, 2) + mpx_senscr(2);
397 t2 = s2 + hmm_tprob_5st(2, 4);
398 }
399
400 t0 = t1 = WORST_SCORE;
401 if (s4 != WORST_SCORE)
402 t0 = s4 + hmm_tprob_5st(4, 4);
403 if (s3 != WORST_SCORE)
404 t1 = s3 + hmm_tprob_5st(3, 4);
405 if (t0 BETTER_THAN t1) {
406 if (t2 BETTER_THAN t0) {
407 s4 = t2;
408 hmm_history(hmm, 4) = hmm_history(hmm, 2);
409 ssid[4] = ssid[2];
410 }
411 else
412 s4 = t0;
413 }
414 else {
415 if (t2 BETTER_THAN t1) {
416 s4 = t2;
417 hmm_history(hmm, 4) = hmm_history(hmm, 2);
418 ssid[4] = ssid[2];
419 }
420 else {
421 s4 = t1;
422 hmm_history(hmm, 4) = hmm_history(hmm, 3);
423 ssid[4] = ssid[3];
424 }
425 }
426 if (s4 WORSE_THAN WORST_SCORE) s4 = WORST_SCORE;
427 if (s4 BETTER_THAN bestScore)
428 bestScore = s4;
429 hmm_score(hmm, 4) = s4;
430
431 /* Don't propagate WORST_SCORE */
432 if (ssid[1] == BAD_SSID)
433 s1 = t2 = WORST_SCORE;
434 else {
435 s1 = hmm_score(hmm, 1) + mpx_senscr(1);
436 t2 = s1 + hmm_tprob_5st(1, 3);
437 }
438 t0 = t1 = WORST_SCORE;
439 if (s3 != WORST_SCORE)
440 t0 = s3 + hmm_tprob_5st(3, 3);
441 if (s2 != WORST_SCORE)
442 t1 = s2 + hmm_tprob_5st(2, 3);
443 if (t0 BETTER_THAN t1) {
444 if (t2 BETTER_THAN t0) {
445 s3 = t2;
446 hmm_history(hmm, 3) = hmm_history(hmm, 1);
447 ssid[3] = ssid[1];
448 }
449 else
450 s3 = t0;
451 }
452 else {
453 if (t2 BETTER_THAN t1) {
454 s3 = t2;
455 hmm_history(hmm, 3) = hmm_history(hmm, 1);
456 ssid[3] = ssid[1];
457 }
458 else {
459 s3 = t1;
460 hmm_history(hmm, 3) = hmm_history(hmm, 2);
461 ssid[3] = ssid[2];
462 }
463 }
464 if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
465 if (s3 BETTER_THAN bestScore) bestScore = s3;
466 hmm_score(hmm, 3) = s3;
467
468 /* State 0 is always active */
469 s0 = hmm_in_score(hmm) + mpx_senscr(0);
470
471 /* Don't propagate WORST_SCORE */
472 t0 = t1 = WORST_SCORE;
473 if (s2 != WORST_SCORE)
474 t0 = s2 + hmm_tprob_5st(2, 2);
475 if (s1 != WORST_SCORE)
476 t1 = s1 + hmm_tprob_5st(1, 2);
477 t2 = s0 + hmm_tprob_5st(0, 2);
478 if (t0 BETTER_THAN t1) {
479 if (t2 BETTER_THAN t0) {
480 s2 = t2;
481 hmm_history(hmm, 2) = hmm_in_history(hmm);
482 ssid[2] = ssid[0];
483 }
484 else
485 s2 = t0;
486 }
487 else {
488 if (t2 BETTER_THAN t1) {
489 s2 = t2;
490 hmm_history(hmm, 2) = hmm_in_history(hmm);
491 ssid[2] = ssid[0];
492 }
493 else {
494 s2 = t1;
495 hmm_history(hmm, 2) = hmm_history(hmm, 1);
496 ssid[2] = ssid[1];
497 }
498 }
499 if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
500 if (s2 BETTER_THAN bestScore) bestScore = s2;
501 hmm_score(hmm, 2) = s2;
502
503 /* Don't propagate WORST_SCORE */
504 t0 = WORST_SCORE;
505 if (s1 != WORST_SCORE)
506 t0 = s1 + hmm_tprob_5st(1, 1);
507 t1 = s0 + hmm_tprob_5st(0, 1);
508 if (t0 BETTER_THAN t1) {
509 s1 = t0;
510 }
511 else {
512 s1 = t1;
513 hmm_history(hmm, 1) = hmm_in_history(hmm);
514 ssid[1] = ssid[0];
515 }
516 if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
517 if (s1 BETTER_THAN bestScore) bestScore = s1;
518 hmm_score(hmm, 1) = s1;
519
520 s0 += hmm_tprob_5st(0, 0);
521 if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
522 if (s0 BETTER_THAN bestScore) bestScore = s0;
523 hmm_in_score(hmm) = s0;
524
525 hmm_bestscore(hmm) = bestScore;
526 return bestScore;
527 }
528
529 #define hmm_tprob_3st(i, j) (-tp[(i)*4+(j)])
530
531 static int32
hmm_vit_eval_3st_lr(hmm_t * hmm)532 hmm_vit_eval_3st_lr(hmm_t * hmm)
533 {
534 int16 const *senscore = hmm->ctx->senscore;
535 uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
536 uint16 const *sseq = hmm->senid;
537 int32 s3, s2, s1, s0, t2, t1, t0, bestScore;
538
539 s2 = hmm_score(hmm, 2) + nonmpx_senscr(2);
540 s1 = hmm_score(hmm, 1) + nonmpx_senscr(1);
541 s0 = hmm_in_score(hmm) + nonmpx_senscr(0);
542
543 /* It was the best of scores, it was the worst of scores. */
544 bestScore = WORST_SCORE;
545 t2 = INT_MIN; /* Not used unless skipstate is true */
546
547 /* Transitions into non-emitting state 3 */
548 if (s1 BETTER_THAN WORST_SCORE) {
549 t1 = s2 + hmm_tprob_3st(2, 3);
550 if (hmm_tprob_3st(1,3) BETTER_THAN TMAT_WORST_SCORE)
551 t2 = s1 + hmm_tprob_3st(1, 3);
552 if (t1 BETTER_THAN t2) {
553 s3 = t1;
554 hmm_out_history(hmm) = hmm_history(hmm, 2);
555 } else {
556 s3 = t2;
557 hmm_out_history(hmm) = hmm_history(hmm, 1);
558 }
559 if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
560 hmm_out_score(hmm) = s3;
561 bestScore = s3;
562 }
563
564 /* All transitions into state 2 (state 0 is always active) */
565 t0 = s2 + hmm_tprob_3st(2, 2);
566 t1 = s1 + hmm_tprob_3st(1, 2);
567 if (hmm_tprob_3st(0, 2) BETTER_THAN TMAT_WORST_SCORE)
568 t2 = s0 + hmm_tprob_3st(0, 2);
569 if (t0 BETTER_THAN t1) {
570 if (t2 BETTER_THAN t0) {
571 s2 = t2;
572 hmm_history(hmm, 2) = hmm_in_history(hmm);
573 } else
574 s2 = t0;
575 } else {
576 if (t2 BETTER_THAN t1) {
577 s2 = t2;
578 hmm_history(hmm, 2) = hmm_in_history(hmm);
579 } else {
580 s2 = t1;
581 hmm_history(hmm, 2) = hmm_history(hmm, 1);
582 }
583 }
584 if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
585 if (s2 BETTER_THAN bestScore) bestScore = s2;
586 hmm_score(hmm, 2) = s2;
587
588 /* All transitions into state 1 */
589 t0 = s1 + hmm_tprob_3st(1, 1);
590 t1 = s0 + hmm_tprob_3st(0, 1);
591 if (t0 BETTER_THAN t1) {
592 s1 = t0;
593 } else {
594 s1 = t1;
595 hmm_history(hmm, 1) = hmm_in_history(hmm);
596 }
597 if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
598 if (s1 BETTER_THAN bestScore) bestScore = s1;
599 hmm_score(hmm, 1) = s1;
600
601 /* All transitions into state 0 */
602 s0 = s0 + hmm_tprob_3st(0, 0);
603 if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
604 if (s0 BETTER_THAN bestScore) bestScore = s0;
605 hmm_in_score(hmm) = s0;
606
607 hmm_bestscore(hmm) = bestScore;
608 return bestScore;
609 }
610
611 static int32
hmm_vit_eval_3st_lr_mpx(hmm_t * hmm)612 hmm_vit_eval_3st_lr_mpx(hmm_t * hmm)
613 {
614 uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
615 int16 const *senscore = hmm->ctx->senscore;
616 uint16 * const *sseq = hmm->ctx->sseq;
617 uint16 *ssid = hmm->senid;
618 int32 bestScore;
619 int32 s3, s2, s1, s0, t2, t1, t0;
620
621 /* Don't propagate WORST_SCORE */
622 t2 = INT_MIN; /* Not used unless skipstate is true */
623 if (ssid[2] == BAD_SSID)
624 s2 = t1 = WORST_SCORE;
625 else {
626 s2 = hmm_score(hmm, 2) + mpx_senscr(2);
627 t1 = s2 + hmm_tprob_3st(2, 3);
628 }
629 if (ssid[1] == BAD_SSID)
630 s1 = t2 = WORST_SCORE;
631 else {
632 s1 = hmm_score(hmm, 1) + mpx_senscr(1);
633 if (hmm_tprob_3st(1,3) BETTER_THAN TMAT_WORST_SCORE)
634 t2 = s1 + hmm_tprob_3st(1, 3);
635 }
636 if (t1 BETTER_THAN t2) {
637 s3 = t1;
638 hmm_out_history(hmm) = hmm_history(hmm, 2);
639 }
640 else {
641 s3 = t2;
642 hmm_out_history(hmm) = hmm_history(hmm, 1);
643 }
644 if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
645 hmm_out_score(hmm) = s3;
646 bestScore = s3;
647
648 /* State 0 is always active */
649 s0 = hmm_in_score(hmm) + mpx_senscr(0);
650
651 /* Don't propagate WORST_SCORE */
652 t0 = t1 = WORST_SCORE;
653 if (s2 != WORST_SCORE)
654 t0 = s2 + hmm_tprob_3st(2, 2);
655 if (s1 != WORST_SCORE)
656 t1 = s1 + hmm_tprob_3st(1, 2);
657 if (hmm_tprob_3st(0,2) BETTER_THAN TMAT_WORST_SCORE)
658 t2 = s0 + hmm_tprob_3st(0, 2);
659 if (t0 BETTER_THAN t1) {
660 if (t2 BETTER_THAN t0) {
661 s2 = t2;
662 hmm_history(hmm, 2) = hmm_in_history(hmm);
663 ssid[2] = ssid[0];
664 }
665 else
666 s2 = t0;
667 }
668 else {
669 if (t2 BETTER_THAN t1) {
670 s2 = t2;
671 hmm_history(hmm, 2) = hmm_in_history(hmm);
672 ssid[2] = ssid[0];
673 }
674 else {
675 s2 = t1;
676 hmm_history(hmm, 2) = hmm_history(hmm, 1);
677 ssid[2] = ssid[1];
678 }
679 }
680 if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
681 if (s2 BETTER_THAN bestScore) bestScore = s2;
682 hmm_score(hmm, 2) = s2;
683
684 /* Don't propagate WORST_SCORE */
685 t0 = WORST_SCORE;
686 if (s1 != WORST_SCORE)
687 t0 = s1 + hmm_tprob_3st(1, 1);
688 t1 = s0 + hmm_tprob_3st(0, 1);
689 if (t0 BETTER_THAN t1) {
690 s1 = t0;
691 }
692 else {
693 s1 = t1;
694 hmm_history(hmm, 1) = hmm_in_history(hmm);
695 ssid[1] = ssid[0];
696 }
697 if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
698 if (s1 BETTER_THAN bestScore) bestScore = s1;
699 hmm_score(hmm, 1) = s1;
700
701 /* State 0 is always active */
702 s0 += hmm_tprob_3st(0, 0);
703 if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
704 if (s0 BETTER_THAN bestScore) bestScore = s0;
705 hmm_in_score(hmm) = s0;
706
707 hmm_bestscore(hmm) = bestScore;
708 return bestScore;
709 }
710
711 static int32
hmm_vit_eval_anytopo(hmm_t * hmm)712 hmm_vit_eval_anytopo(hmm_t * hmm)
713 {
714 hmm_context_t *ctx = hmm->ctx;
715 int32 to, from, bestfrom;
716 int32 newscr, scr, bestscr;
717 int final_state;
718
719 /* Compute previous state-score + observation output prob for each emitting state */
720 ctx->st_sen_scr[0] = hmm_in_score(hmm) + hmm_senscr(hmm, 0);
721 for (from = 1; from < hmm_n_emit_state(hmm); ++from) {
722 if ((ctx->st_sen_scr[from] =
723 hmm_score(hmm, from) + hmm_senscr(hmm, from)) WORSE_THAN WORST_SCORE)
724 ctx->st_sen_scr[from] = WORST_SCORE;
725 }
726
727 /* FIXME/TODO: Use the BLAS for all this. */
728 /* Evaluate final-state first, which does not have a self-transition */
729 final_state = hmm_n_emit_state(hmm);
730 to = final_state;
731 scr = WORST_SCORE;
732 bestfrom = -1;
733 for (from = to - 1; from >= 0; --from) {
734 if ((hmm_tprob(hmm, from, to) BETTER_THAN TMAT_WORST_SCORE) &&
735 ((newscr = ctx->st_sen_scr[from]
736 + hmm_tprob(hmm, from, to)) BETTER_THAN scr)) {
737 scr = newscr;
738 bestfrom = from;
739 }
740 }
741 hmm_out_score(hmm) = scr;
742 if (bestfrom >= 0)
743 hmm_out_history(hmm) = hmm_history(hmm, bestfrom);
744 bestscr = scr;
745
746 /* Evaluate all other states, which might have self-transitions */
747 for (to = final_state - 1; to >= 0; --to) {
748 /* Score from self-transition, if any */
749 scr =
750 (hmm_tprob(hmm, to, to) BETTER_THAN TMAT_WORST_SCORE)
751 ? ctx->st_sen_scr[to] + hmm_tprob(hmm, to, to)
752 : WORST_SCORE;
753
754 /* Scores from transitions from other states */
755 bestfrom = -1;
756 for (from = to - 1; from >= 0; --from) {
757 if ((hmm_tprob(hmm, from, to) BETTER_THAN TMAT_WORST_SCORE) &&
758 ((newscr = ctx->st_sen_scr[from]
759 + hmm_tprob(hmm, from, to)) BETTER_THAN scr)) {
760 scr = newscr;
761 bestfrom = from;
762 }
763 }
764
765 /* Update new result for state to */
766 if (to == 0) {
767 hmm_in_score(hmm) = scr;
768 if (bestfrom >= 0)
769 hmm_in_history(hmm) = hmm_history(hmm, bestfrom);
770 }
771 else {
772 hmm_score(hmm, to) = scr;
773 if (bestfrom >= 0)
774 hmm_history(hmm, to) = hmm_history(hmm, bestfrom);
775 }
776 /* Propagate ssid for multiplex HMMs */
777 if (bestfrom >= 0 && hmm_is_mpx(hmm))
778 hmm->senid[to] = hmm->senid[bestfrom];
779
780 if (bestscr WORSE_THAN scr)
781 bestscr = scr;
782 }
783
784 hmm_bestscore(hmm) = bestscr;
785 return bestscr;
786 }
787
788 int32
hmm_vit_eval(hmm_t * hmm)789 hmm_vit_eval(hmm_t * hmm)
790 {
791 if (hmm_is_mpx(hmm)) {
792 if (hmm_n_emit_state(hmm) == 5)
793 return hmm_vit_eval_5st_lr_mpx(hmm);
794 else if (hmm_n_emit_state(hmm) == 3)
795 return hmm_vit_eval_3st_lr_mpx(hmm);
796 else
797 return hmm_vit_eval_anytopo(hmm);
798 }
799 else {
800 if (hmm_n_emit_state(hmm) == 5)
801 return hmm_vit_eval_5st_lr(hmm);
802 else if (hmm_n_emit_state(hmm) == 3)
803 return hmm_vit_eval_3st_lr(hmm);
804 else
805 return hmm_vit_eval_anytopo(hmm);
806 }
807 }
808
809 int32
hmm_dump_vit_eval(hmm_t * hmm,FILE * fp)810 hmm_dump_vit_eval(hmm_t * hmm, FILE * fp)
811 {
812 int32 bs = 0;
813
814 if (fp) {
815 fprintf(fp, "BEFORE:\n");
816 hmm_dump(hmm, fp);
817 }
818 bs = hmm_vit_eval(hmm);
819 if (fp) {
820 fprintf(fp, "AFTER:\n");
821 hmm_dump(hmm, fp);
822 }
823
824 return bs;
825 }
826