PocketSphinx  0.6
src/libpocketsphinx/hmm.c
00001 /* -*- c-basic-offset: 4; indent-tabs-mode: nil -*- */
00002 /* ====================================================================
00003  * Copyright (c) 1999-2004 Carnegie Mellon University.  All rights
00004  * reserved.
00005  *
00006  * Redistribution and use in source and binary forms, with or without
00007  * modification, are permitted provided that the following conditions
00008  * are met:
00009  *
00010  * 1. Redistributions of source code must retain the above copyright
00011  *    notice, this list of conditions and the following disclaimer. 
00012  *
00013  * 2. Redistributions in binary form must reproduce the above copyright
00014  *    notice, this list of conditions and the following disclaimer in
00015  *    the documentation and/or other materials provided with the
00016  *    distribution.
00017  *
00018  * This work was supported in part by funding from the Defense Advanced 
00019  * Research Projects Agency and the National Science Foundation of the 
00020  * United States of America, and the CMU Sphinx Speech Consortium.
00021  *
00022  * THIS SOFTWARE IS PROVIDED BY CARNEGIE MELLON UNIVERSITY ``AS IS'' AND 
00023  * ANY EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
00024  * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00025  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL CARNEGIE MELLON UNIVERSITY
00026  * NOR ITS EMPLOYEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00027  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
00028  * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
00029  * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
00030  * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
00031  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
00032  * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00033  *
00034  * ====================================================================
00035  *
00036  */
00037 
00042 /* System headers. */
00043 #include <assert.h>
00044 #include <stdlib.h>
00045 #include <string.h>
00046 #include <limits.h>
00047 
00048 /* SphinxBase headers. */
00049 #include <sphinxbase/ckd_alloc.h>
00050 #include <sphinxbase/err.h>
00051 
00052 /* Local headers. */
00053 #include "hmm.h"
00054 
00055 hmm_context_t *
00056 hmm_context_init(int32 n_emit_state,
00057                  uint8 ** const *tp,
00058                  int16 const *senscore,
00059                  uint16 * const *sseq)
00060 {
00061     hmm_context_t *ctx;
00062 
00063     assert(n_emit_state > 0);
00064     if (n_emit_state > HMM_MAX_NSTATE) {
00065         E_ERROR("Number of emitting states must be <= %d\n", HMM_MAX_NSTATE);
00066         return NULL;
00067     }
00068 
00069     ctx = ckd_calloc(1, sizeof(*ctx));
00070     ctx->n_emit_state = n_emit_state;
00071     ctx->tp = tp;
00072     ctx->senscore = senscore;
00073     ctx->sseq = sseq;
00074     ctx->st_sen_scr = ckd_calloc(n_emit_state, sizeof(*ctx->st_sen_scr));
00075 
00076     return ctx;
00077 }
00078 
00079 void
00080 hmm_context_free(hmm_context_t *ctx)
00081 {
00082     if (ctx == NULL)
00083         return;
00084     ckd_free(ctx->st_sen_scr);
00085     ckd_free(ctx);
00086 }
00087 
00088 void
00089 hmm_init(hmm_context_t *ctx, hmm_t *hmm, int mpx, int ssid, int tmatid)
00090 {
00091     hmm->ctx = ctx;
00092     hmm->mpx = mpx;
00093     hmm->n_emit_state = ctx->n_emit_state;
00094     if (mpx) {
00095         int i;
00096         hmm->ssid = BAD_SSID;
00097         hmm->senid[0] = ssid;
00098         for (i = 1; i < hmm_n_emit_state(hmm); ++i) {
00099             hmm->senid[i] = BAD_SSID;
00100         }
00101     }
00102     else {
00103         hmm->ssid = ssid;
00104         memcpy(hmm->senid, ctx->sseq[ssid], hmm->n_emit_state * sizeof(*hmm->senid));
00105     }
00106     hmm->tmatid = tmatid;
00107     hmm_clear(hmm);
00108 }
00109 
00110 void
00111 hmm_deinit(hmm_t *hmm)
00112 {
00113 }
00114 
00115 void
00116 hmm_dump(hmm_t * hmm,
00117          FILE * fp)
00118 {
00119     int32 i;
00120 
00121     if (hmm_is_mpx(hmm)) {
00122         fprintf(fp, "MPX   ");
00123         for (i = 0; i < hmm_n_emit_state(hmm); i++)
00124             fprintf(fp, " %11d", hmm_senid(hmm, i));
00125         fprintf(fp, " ( ");
00126         for (i = 0; i < hmm_n_emit_state(hmm); i++)
00127             fprintf(fp, "%d ", hmm_ssid(hmm, i));
00128         fprintf(fp, ")\n");
00129     }
00130     else {
00131         fprintf(fp, "SSID  ");
00132         for (i = 0; i < hmm_n_emit_state(hmm); i++)
00133             fprintf(fp, " %11d", hmm_senid(hmm, i));
00134         fprintf(fp, " (%d)\n", hmm_ssid(hmm, 0));
00135     }
00136 
00137     if (hmm->ctx->senscore) {
00138         fprintf(fp, "SENSCR");
00139         for (i = 0; i < hmm_n_emit_state(hmm); i++)
00140             fprintf(fp, " %11d", hmm_senscr(hmm, i));
00141         fprintf(fp, "\n");
00142     }
00143 
00144     fprintf(fp, "SCORES %11d", hmm_in_score(hmm));
00145     for (i = 1; i < hmm_n_emit_state(hmm); i++)
00146         fprintf(fp, " %11d", hmm_score(hmm, i));
00147     fprintf(fp, " %11d", hmm_out_score(hmm));
00148     fprintf(fp, "\n");
00149 
00150     fprintf(fp, "HISTID %11d", hmm_in_history(hmm));
00151     for (i = 1; i < hmm_n_emit_state(hmm); i++)
00152         fprintf(fp, " %11d", hmm_history(hmm, i));
00153     fprintf(fp, " %11d", hmm_out_history(hmm));
00154     fprintf(fp, "\n");
00155 
00156     if (hmm_in_score(hmm) > 0)
00157         fprintf(fp,
00158                 "ALERT!! The input score %d is large than 0. Probably wrap around.\n",
00159                 hmm_in_score(hmm));
00160     if (hmm_out_score(hmm) > 0)
00161         fprintf(fp,
00162                 "ALERT!! The output score %d is large than 0. Probably wrap around\n.",
00163                 hmm_out_score(hmm));
00164 
00165     fflush(fp);
00166 }
00167 
00168 
00169 void
00170 hmm_clear_scores(hmm_t * h)
00171 {
00172     int32 i;
00173 
00174     hmm_in_score(h) = WORST_SCORE;
00175     for (i = 1; i < hmm_n_emit_state(h); i++)
00176         hmm_score(h, i) = WORST_SCORE;
00177     hmm_out_score(h) = WORST_SCORE;
00178 
00179     h->bestscore = WORST_SCORE;
00180 }
00181 
00182 void
00183 hmm_clear(hmm_t * h)
00184 {
00185     int32 i;
00186 
00187     hmm_in_score(h) = WORST_SCORE;
00188     hmm_in_history(h) = -1;
00189     for (i = 1; i < hmm_n_emit_state(h); i++) {
00190         hmm_score(h, i) = WORST_SCORE;
00191         hmm_history(h, i) = -1;
00192     }
00193     hmm_out_score(h) = WORST_SCORE;
00194     hmm_out_history(h) = -1;
00195 
00196     h->bestscore = WORST_SCORE;
00197     h->frame = -1;
00198 }
00199 
00200 void
00201 hmm_enter(hmm_t *h, int32 score, int32 histid, int frame)
00202 {
00203     hmm_in_score(h) = score;
00204     hmm_in_history(h) = histid;
00205     hmm_frame(h) = frame;
00206 }
00207 
00208 void
00209 hmm_normalize(hmm_t *h, int32 bestscr)
00210 {
00211     int32 i;
00212 
00213     for (i = 0; i < hmm_n_emit_state(h); i++) {
00214         if (hmm_score(h, i) BETTER_THAN WORST_SCORE)
00215             hmm_score(h, i) -= bestscr;
00216     }
00217     if (hmm_out_score(h) BETTER_THAN WORST_SCORE)
00218         hmm_out_score(h) -= bestscr;
00219 }
00220 
00221 #define hmm_tprob_5st(i, j) (-tp[(i)*6+(j)])
00222 #define nonmpx_senscr(i) (-senscore[sseq[i]])
00223 
00224 static int32
00225 hmm_vit_eval_5st_lr(hmm_t * hmm)
00226 {
00227     int16 const *senscore = hmm->ctx->senscore;
00228     uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
00229     uint16 const *sseq = hmm->senid;
00230     int32 s5, s4, s3, s2, s1, s0, t2, t1, t0, bestScore;
00231 
00232     /* It was the best of scores, it was the worst of scores. */
00233     bestScore = WORST_SCORE;
00234 
00235     /* Cache problem here! */
00236     s4 = hmm_score(hmm, 4) + nonmpx_senscr(4);
00237     s3 = hmm_score(hmm, 3) + nonmpx_senscr(3);
00238     /* Transitions into non-emitting state 5 */
00239     if (s3 BETTER_THAN WORST_SCORE) {
00240         t1 = s4 + hmm_tprob_5st(4, 5);
00241         t2 = s3 + hmm_tprob_5st(3, 5);
00242         if (t1 BETTER_THAN t2) {
00243             s5 = t1;
00244             hmm_out_history(hmm)  = hmm_history(hmm, 4);
00245         } else {
00246             s5 = t2;
00247             hmm_out_history(hmm)  = hmm_history(hmm, 3);
00248         }
00249         if (s5 WORSE_THAN WORST_SCORE) s5 = WORST_SCORE;
00250         hmm_out_score(hmm) = s5;
00251         bestScore = s5;
00252     }
00253 
00254     s2 = hmm_score(hmm, 2) + nonmpx_senscr(2);
00255     /* All transitions into state 4 */
00256     if (s2 BETTER_THAN WORST_SCORE) {
00257         t0 = s4 + hmm_tprob_5st(4, 4);
00258         t1 = s3 + hmm_tprob_5st(3, 4);
00259         t2 = s2 + hmm_tprob_5st(2, 4);
00260         if (t0 BETTER_THAN t1) {
00261             if (t2 BETTER_THAN t0) {
00262                 s4 = t2;
00263                 hmm_history(hmm, 4)  = hmm_history(hmm, 2);
00264             } else
00265                 s4 = t0;
00266         } else {
00267             if (t2 BETTER_THAN t1) {
00268                 s4 = t2;
00269                 hmm_history(hmm, 4)  = hmm_history(hmm, 2);
00270             } else {
00271                 s4 = t1;
00272                 hmm_history(hmm, 4)  = hmm_history(hmm, 3);
00273             }
00274         }
00275         if (s4 WORSE_THAN WORST_SCORE) s4 = WORST_SCORE;
00276         if (s4 BETTER_THAN bestScore) bestScore = s4;
00277         hmm_score(hmm, 4) = s4;
00278     }
00279 
00280     s1 = hmm_score(hmm, 1) + nonmpx_senscr(1);
00281     /* All transitions into state 3 */
00282     if (s1 BETTER_THAN WORST_SCORE) {
00283         t0 = s3 + hmm_tprob_5st(3, 3);
00284         t1 = s2 + hmm_tprob_5st(2, 3);
00285         t2 = s1 + hmm_tprob_5st(1, 3);
00286         if (t0 BETTER_THAN t1) {
00287             if (t2 BETTER_THAN t0) {
00288                 s3 = t2;
00289                 hmm_history(hmm, 3)  = hmm_history(hmm, 1);
00290             } else
00291                 s3 = t0;
00292         } else {
00293             if (t2 BETTER_THAN t1) {
00294                 s3 = t2;
00295                 hmm_history(hmm, 3)  = hmm_history(hmm, 1);
00296             } else {
00297                 s3 = t1;
00298                 hmm_history(hmm, 3)  = hmm_history(hmm, 2);
00299             }
00300         }
00301         if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
00302         if (s3 BETTER_THAN bestScore) bestScore = s3;
00303         hmm_score(hmm, 3) = s3;
00304     }
00305 
00306     s0 = hmm_in_score(hmm) + nonmpx_senscr(0);
00307     /* All transitions into state 2 (state 0 is always active) */
00308     t0 = s2 + hmm_tprob_5st(2, 2);
00309     t1 = s1 + hmm_tprob_5st(1, 2);
00310     t2 = s0 + hmm_tprob_5st(0, 2);
00311     if (t0 BETTER_THAN t1) {
00312         if (t2 BETTER_THAN t0) {
00313             s2 = t2;
00314             hmm_history(hmm, 2)  = hmm_in_history(hmm);
00315         } else
00316             s2 = t0;
00317     } else {
00318         if (t2 BETTER_THAN t1) {
00319             s2 = t2;
00320             hmm_history(hmm, 2)  = hmm_in_history(hmm);
00321         } else {
00322             s2 = t1;
00323             hmm_history(hmm, 2)  = hmm_history(hmm, 1);
00324         }
00325     }
00326     if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
00327     if (s2 BETTER_THAN bestScore) bestScore = s2;
00328     hmm_score(hmm, 2) = s2;
00329 
00330 
00331     /* All transitions into state 1 */
00332     t0 = s1 + hmm_tprob_5st(1, 1);
00333     t1 = s0 + hmm_tprob_5st(0, 1);
00334     if (t0 BETTER_THAN t1) {
00335         s1 = t0;
00336     } else {
00337         s1 = t1;
00338         hmm_history(hmm, 1)  = hmm_in_history(hmm);
00339     }
00340     if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
00341     if (s1 BETTER_THAN bestScore) bestScore = s1;
00342     hmm_score(hmm, 1) = s1;
00343 
00344     /* All transitions into state 0 */
00345     s0 = s0 + hmm_tprob_5st(0, 0);
00346     if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
00347     if (s0 BETTER_THAN bestScore) bestScore = s0;
00348     hmm_in_score(hmm) = s0;
00349 
00350     hmm_bestscore(hmm) = bestScore;
00351     return bestScore;
00352 }
00353 
00354 #define mpx_senid(st) sseq[ssid[st]][st]
00355 #define mpx_senscr(st) (-senscore[mpx_senid(st)])
00356 
00357 static int32
00358 hmm_vit_eval_5st_lr_mpx(hmm_t * hmm)
00359 {
00360     uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
00361     int16 const *senscore = hmm->ctx->senscore;
00362     uint16 * const *sseq = hmm->ctx->sseq;
00363     uint16 *ssid = hmm->senid;
00364     int32 bestScore;
00365     int32 s5, s4, s3, s2, s1, s0, t2, t1, t0;
00366 
00367     /* Don't propagate WORST_SCORE */
00368     if (ssid[4] == BAD_SSID)
00369         s4 = t1 = WORST_SCORE;
00370     else {
00371         s4 = hmm_score(hmm, 4) + mpx_senscr(4);
00372         t1 = s4 + hmm_tprob_5st(4, 5);
00373     }
00374     if (ssid[3] == BAD_SSID)
00375         s3 = t2 = WORST_SCORE;
00376     else {
00377         s3 = hmm_score(hmm, 3) + mpx_senscr(3);
00378         t2 = s3 + hmm_tprob_5st(3, 5);
00379     }
00380     if (t1 BETTER_THAN t2) {
00381         s5 = t1;
00382         hmm_out_history(hmm) = hmm_history(hmm, 4);
00383     }
00384     else {
00385         s5 = t2;
00386         hmm_out_history(hmm) = hmm_history(hmm, 3);
00387     }
00388     if (s5 WORSE_THAN WORST_SCORE) s5 = WORST_SCORE;
00389     hmm_out_score(hmm) = s5;
00390     bestScore = s5;
00391 
00392     /* Don't propagate WORST_SCORE */
00393     if (ssid[2] == BAD_SSID)
00394         s2 = t2 = WORST_SCORE;
00395     else {
00396         s2 = hmm_score(hmm, 2) + mpx_senscr(2);
00397         t2 = s2 + hmm_tprob_5st(2, 4);
00398     }
00399 
00400     t0 = t1 = WORST_SCORE;
00401     if (s4 != WORST_SCORE)
00402         t0 = s4 + hmm_tprob_5st(4, 4);
00403     if (s3 != WORST_SCORE)
00404         t1 = s3 + hmm_tprob_5st(3, 4);
00405     if (t0 BETTER_THAN t1) {
00406         if (t2 BETTER_THAN t0) {
00407             s4 = t2;
00408             hmm_history(hmm, 4) = hmm_history(hmm, 2);
00409             ssid[4] = ssid[2];
00410         }
00411         else
00412             s4 = t0;
00413     }
00414     else {
00415         if (t2 BETTER_THAN t1) {
00416             s4 = t2;
00417             hmm_history(hmm, 4) = hmm_history(hmm, 2);
00418             ssid[4] = ssid[2];
00419         }
00420         else {
00421             s4 = t1;
00422             hmm_history(hmm, 4) = hmm_history(hmm, 3);
00423             ssid[4] = ssid[3];
00424         }
00425     }
00426     if (s4 WORSE_THAN WORST_SCORE) s4 = WORST_SCORE;
00427     if (s4 BETTER_THAN bestScore)
00428         bestScore = s4;
00429     hmm_score(hmm, 4) = s4;
00430 
00431     /* Don't propagate WORST_SCORE */
00432     if (ssid[1] == BAD_SSID)
00433         s1 = t2 = WORST_SCORE;
00434     else {
00435         s1 = hmm_score(hmm, 1) + mpx_senscr(1);
00436         t2 = s1 + hmm_tprob_5st(1, 3);
00437     }
00438     t0 = t1 = WORST_SCORE;
00439     if (s3 != WORST_SCORE)
00440         t0 = s3 + hmm_tprob_5st(3, 3);
00441     if (s2 != WORST_SCORE)
00442         t1 = s2 + hmm_tprob_5st(2, 3);
00443     if (t0 BETTER_THAN t1) {
00444         if (t2 BETTER_THAN t0) {
00445             s3 = t2;
00446             hmm_history(hmm, 3) = hmm_history(hmm, 1);
00447             ssid[3] = ssid[1];
00448         }
00449         else
00450             s3 = t0;
00451     }
00452     else {
00453         if (t2 BETTER_THAN t1) {
00454             s3 = t2;
00455             hmm_history(hmm, 3) = hmm_history(hmm, 1);
00456             ssid[3] = ssid[1];
00457         }
00458         else {
00459             s3 = t1;
00460             hmm_history(hmm, 3) = hmm_history(hmm, 2);
00461             ssid[3] = ssid[2];
00462         }
00463     }
00464     if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
00465     if (s3 BETTER_THAN bestScore) bestScore = s3;
00466     hmm_score(hmm, 3) = s3;
00467 
00468     /* State 0 is always active */
00469     s0 = hmm_in_score(hmm) + mpx_senscr(0);
00470 
00471     /* Don't propagate WORST_SCORE */
00472     t0 = t1 = WORST_SCORE;
00473     if (s2 != WORST_SCORE)
00474         t0 = s2 + hmm_tprob_5st(2, 2);
00475     if (s1 != WORST_SCORE)
00476         t1 = s1 + hmm_tprob_5st(1, 2);
00477     t2 = s0 + hmm_tprob_5st(0, 2);
00478     if (t0 BETTER_THAN t1) {
00479         if (t2 BETTER_THAN t0) {
00480             s2 = t2;
00481             hmm_history(hmm, 2) = hmm_in_history(hmm);
00482             ssid[2] = ssid[0];
00483         }
00484         else
00485             s2 = t0;
00486     }
00487     else {
00488         if (t2 BETTER_THAN t1) {
00489             s2 = t2;
00490             hmm_history(hmm, 2) = hmm_in_history(hmm);
00491             ssid[2] = ssid[0];
00492         }
00493         else {
00494             s2 = t1;
00495             hmm_history(hmm, 2) = hmm_history(hmm, 1);
00496             ssid[2] = ssid[1];
00497         }
00498     }
00499     if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
00500     if (s2 BETTER_THAN bestScore) bestScore = s2;
00501     hmm_score(hmm, 2) = s2;
00502 
00503     /* Don't propagate WORST_SCORE */
00504     t0 = WORST_SCORE;
00505     if (s1 != WORST_SCORE)
00506         t0 = s1 + hmm_tprob_5st(1, 1);
00507     t1 = s0 + hmm_tprob_5st(0, 1);
00508     if (t0 BETTER_THAN t1) {
00509         s1 = t0;
00510     }
00511     else {
00512         s1 = t1;
00513         hmm_history(hmm, 1) = hmm_in_history(hmm);
00514         ssid[1] = ssid[0];
00515     }
00516     if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
00517     if (s1 BETTER_THAN bestScore) bestScore = s1;
00518     hmm_score(hmm, 1) = s1;
00519 
00520     s0 += hmm_tprob_5st(0, 0);
00521     if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
00522     if (s0 BETTER_THAN bestScore) bestScore = s0;
00523     hmm_in_score(hmm) = s0;
00524 
00525     hmm_bestscore(hmm) = bestScore;
00526     return bestScore;
00527 }
00528 
00529 #define hmm_tprob_3st(i, j) (-tp[(i)*4+(j)])
00530 
00531 static int32
00532 hmm_vit_eval_3st_lr(hmm_t * hmm)
00533 {
00534     int16 const *senscore = hmm->ctx->senscore;
00535     uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
00536     uint16 const *sseq = hmm->senid;
00537     int32 s3, s2, s1, s0, t2, t1, t0, bestScore;
00538 
00539     s2 = hmm_score(hmm, 2) + nonmpx_senscr(2);
00540     s1 = hmm_score(hmm, 1) + nonmpx_senscr(1);
00541     s0 = hmm_in_score(hmm) + nonmpx_senscr(0);
00542 
00543     /* It was the best of scores, it was the worst of scores. */
00544     bestScore = WORST_SCORE;
00545     t2 = INT_MIN; /* Not used unless skipstate is true */
00546 
00547     /* Transitions into non-emitting state 3 */
00548     if (s1 BETTER_THAN WORST_SCORE) {
00549         t1 = s2 + hmm_tprob_3st(2, 3);
00550         if (hmm_tprob_3st(1,3) BETTER_THAN TMAT_WORST_SCORE)
00551             t2 = s1 + hmm_tprob_3st(1, 3);
00552         if (t1 BETTER_THAN t2) {
00553             s3 = t1;
00554             hmm_out_history(hmm)  = hmm_history(hmm, 2);
00555         } else {
00556             s3 = t2;
00557             hmm_out_history(hmm)  = hmm_history(hmm, 1);
00558         }
00559         if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
00560         hmm_out_score(hmm) = s3;
00561         bestScore = s3;
00562     }
00563 
00564     /* All transitions into state 2 (state 0 is always active) */
00565     t0 = s2 + hmm_tprob_3st(2, 2);
00566     t1 = s1 + hmm_tprob_3st(1, 2);
00567     if (hmm_tprob_3st(0, 2) BETTER_THAN TMAT_WORST_SCORE)
00568         t2 = s0 + hmm_tprob_3st(0, 2);
00569     if (t0 BETTER_THAN t1) {
00570         if (t2 BETTER_THAN t0) {
00571             s2 = t2;
00572             hmm_history(hmm, 2)  = hmm_in_history(hmm);
00573         } else
00574             s2 = t0;
00575     } else {
00576         if (t2 BETTER_THAN t1) {
00577             s2 = t2;
00578             hmm_history(hmm, 2)  = hmm_in_history(hmm);
00579         } else {
00580             s2 = t1;
00581             hmm_history(hmm, 2)  = hmm_history(hmm, 1);
00582         }
00583     }
00584     if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
00585     if (s2 BETTER_THAN bestScore) bestScore = s2;
00586     hmm_score(hmm, 2) = s2;
00587 
00588     /* All transitions into state 1 */
00589     t0 = s1 + hmm_tprob_3st(1, 1);
00590     t1 = s0 + hmm_tprob_3st(0, 1);
00591     if (t0 BETTER_THAN t1) {
00592         s1 = t0;
00593     } else {
00594         s1 = t1;
00595         hmm_history(hmm, 1)  = hmm_in_history(hmm);
00596     }
00597     if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
00598     if (s1 BETTER_THAN bestScore) bestScore = s1;
00599     hmm_score(hmm, 1) = s1;
00600 
00601     /* All transitions into state 0 */
00602     s0 = s0 + hmm_tprob_3st(0, 0);
00603     if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
00604     if (s0 BETTER_THAN bestScore) bestScore = s0;
00605     hmm_in_score(hmm) = s0;
00606 
00607     hmm_bestscore(hmm) = bestScore;
00608     return bestScore;
00609 }
00610 
00611 static int32
00612 hmm_vit_eval_3st_lr_mpx(hmm_t * hmm)
00613 {
00614     uint8 const *tp = hmm->ctx->tp[hmm->tmatid][0];
00615     int16 const *senscore = hmm->ctx->senscore;
00616     uint16 * const *sseq = hmm->ctx->sseq;
00617     uint16 *ssid = hmm->senid;
00618     int32 bestScore;
00619     int32 s3, s2, s1, s0, t2, t1, t0;
00620 
00621     /* Don't propagate WORST_SCORE */
00622     t2 = INT_MIN; /* Not used unless skipstate is true */
00623     if (ssid[2] == BAD_SSID)
00624         s2 = t1 = WORST_SCORE;
00625     else {
00626         s2 = hmm_score(hmm, 2) + mpx_senscr(2);
00627         t1 = s2 + hmm_tprob_3st(2, 3);
00628     }
00629     if (ssid[1] == BAD_SSID)
00630         s1 = t2 = WORST_SCORE;
00631     else {
00632         s1 = hmm_score(hmm, 1) + mpx_senscr(1);
00633         if (hmm_tprob_3st(1,3) BETTER_THAN TMAT_WORST_SCORE)
00634             t2 = s1 + hmm_tprob_3st(1, 3);
00635     }
00636     if (t1 BETTER_THAN t2) {
00637         s3 = t1;
00638         hmm_out_history(hmm) = hmm_history(hmm, 2);
00639     }
00640     else {
00641         s3 = t2;
00642         hmm_out_history(hmm) = hmm_history(hmm, 1);
00643     }
00644     if (s3 WORSE_THAN WORST_SCORE) s3 = WORST_SCORE;
00645     hmm_out_score(hmm) = s3;
00646     bestScore = s3;
00647 
00648     /* State 0 is always active */
00649     s0 = hmm_in_score(hmm) + mpx_senscr(0);
00650 
00651     /* Don't propagate WORST_SCORE */
00652     t0 = t1 = WORST_SCORE;
00653     if (s2 != WORST_SCORE)
00654         t0 = s2 + hmm_tprob_3st(2, 2);
00655     if (s1 != WORST_SCORE)
00656         t1 = s1 + hmm_tprob_3st(1, 2);
00657     if (hmm_tprob_3st(0,2) BETTER_THAN TMAT_WORST_SCORE)
00658         t2 = s0 + hmm_tprob_3st(0, 2);
00659     if (t0 BETTER_THAN t1) {
00660         if (t2 BETTER_THAN t0) {
00661             s2 = t2;
00662             hmm_history(hmm, 2) = hmm_in_history(hmm);
00663             ssid[2] = ssid[0];
00664         }
00665         else
00666             s2 = t0;
00667     }
00668     else {
00669         if (t2 BETTER_THAN t1) {
00670             s2 = t2;
00671             hmm_history(hmm, 2) = hmm_in_history(hmm);
00672             ssid[2] = ssid[0];
00673         }
00674         else {
00675             s2 = t1;
00676             hmm_history(hmm, 2) = hmm_history(hmm, 1);
00677             ssid[2] = ssid[1];
00678         }
00679     }
00680     if (s2 WORSE_THAN WORST_SCORE) s2 = WORST_SCORE;
00681     if (s2 BETTER_THAN bestScore) bestScore = s2;
00682     hmm_score(hmm, 2) = s2;
00683 
00684     /* Don't propagate WORST_SCORE */
00685     t0 = WORST_SCORE;
00686     if (s1 != WORST_SCORE)
00687         t0 = s1 + hmm_tprob_3st(1, 1);
00688     t1 = s0 + hmm_tprob_3st(0, 1);
00689     if (t0 BETTER_THAN t1) {
00690         s1 = t0;
00691     }
00692     else {
00693         s1 = t1;
00694         hmm_history(hmm, 1) = hmm_in_history(hmm);
00695         ssid[1] = ssid[0];
00696     }
00697     if (s1 WORSE_THAN WORST_SCORE) s1 = WORST_SCORE;
00698     if (s1 BETTER_THAN bestScore) bestScore = s1;
00699     hmm_score(hmm, 1) = s1;
00700 
00701     /* State 0 is always active */
00702     s0 += hmm_tprob_3st(0, 0);
00703     if (s0 WORSE_THAN WORST_SCORE) s0 = WORST_SCORE;
00704     if (s0 BETTER_THAN bestScore) bestScore = s0;
00705     hmm_in_score(hmm) = s0;
00706 
00707     hmm_bestscore(hmm) = bestScore;
00708     return bestScore;
00709 }
00710 
00711 static int32
00712 hmm_vit_eval_anytopo(hmm_t * hmm)
00713 {
00714     hmm_context_t *ctx = hmm->ctx;
00715     int32 to, from, bestfrom;
00716     int32 newscr, scr, bestscr;
00717     int final_state;
00718 
00719     /* Compute previous state-score + observation output prob for each emitting state */
00720     ctx->st_sen_scr[0] = hmm_in_score(hmm) + hmm_senscr(hmm, 0);
00721     for (from = 1; from < hmm_n_emit_state(hmm); ++from) {
00722         if ((ctx->st_sen_scr[from] =
00723              hmm_score(hmm, from) + hmm_senscr(hmm, from)) WORSE_THAN WORST_SCORE)
00724             ctx->st_sen_scr[from] = WORST_SCORE;
00725     }
00726 
00727     /* FIXME/TODO: Use the BLAS for all this. */
00728     /* Evaluate final-state first, which does not have a self-transition */
00729     final_state = hmm_n_emit_state(hmm);
00730     to = final_state;
00731     scr = WORST_SCORE;
00732     bestfrom = -1;
00733     for (from = to - 1; from >= 0; --from) {
00734         if ((hmm_tprob(hmm, from, to) BETTER_THAN TMAT_WORST_SCORE) &&
00735             ((newscr = ctx->st_sen_scr[from]
00736               + hmm_tprob(hmm, from, to)) BETTER_THAN scr)) {
00737             scr = newscr;
00738             bestfrom = from;
00739         }
00740     }
00741     hmm_out_score(hmm) = scr;
00742     if (bestfrom >= 0)
00743         hmm_out_history(hmm) = hmm_history(hmm, bestfrom);
00744     bestscr = scr;
00745 
00746     /* Evaluate all other states, which might have self-transitions */
00747     for (to = final_state - 1; to >= 0; --to) {
00748         /* Score from self-transition, if any */
00749         scr =
00750             (hmm_tprob(hmm, to, to) BETTER_THAN TMAT_WORST_SCORE)
00751             ? ctx->st_sen_scr[to] + hmm_tprob(hmm, to, to)
00752             : WORST_SCORE;
00753 
00754         /* Scores from transitions from other states */
00755         bestfrom = -1;
00756         for (from = to - 1; from >= 0; --from) {
00757             if ((hmm_tprob(hmm, from, to) BETTER_THAN TMAT_WORST_SCORE) &&
00758                 ((newscr = ctx->st_sen_scr[from]
00759                   + hmm_tprob(hmm, from, to)) BETTER_THAN scr)) {
00760                 scr = newscr;
00761                 bestfrom = from;
00762             }
00763         }
00764 
00765         /* Update new result for state to */
00766         if (to == 0) {
00767             hmm_in_score(hmm) = scr;
00768             if (bestfrom >= 0)
00769                 hmm_in_history(hmm) = hmm_history(hmm, bestfrom);
00770         }
00771         else {
00772             hmm_score(hmm, to) = scr;
00773             if (bestfrom >= 0)
00774                 hmm_history(hmm, to) = hmm_history(hmm, bestfrom);
00775         }
00776         /* Propagate ssid for multiplex HMMs */
00777         if (bestfrom >= 0 && hmm_is_mpx(hmm))
00778             hmm->senid[to] = hmm->senid[bestfrom];
00779 
00780         if (bestscr WORSE_THAN scr)
00781             bestscr = scr;
00782     }
00783 
00784     hmm_bestscore(hmm) = bestscr;
00785     return bestscr;
00786 }
00787 
00788 int32
00789 hmm_vit_eval(hmm_t * hmm)
00790 {
00791     if (hmm_is_mpx(hmm)) {
00792         if (hmm_n_emit_state(hmm) == 5)
00793             return hmm_vit_eval_5st_lr_mpx(hmm);
00794         else if (hmm_n_emit_state(hmm) == 3)
00795             return hmm_vit_eval_3st_lr_mpx(hmm);
00796         else
00797             return hmm_vit_eval_anytopo(hmm);
00798     }
00799     else {
00800         if (hmm_n_emit_state(hmm) == 5)
00801             return hmm_vit_eval_5st_lr(hmm);
00802         else if (hmm_n_emit_state(hmm) == 3)
00803             return hmm_vit_eval_3st_lr(hmm);
00804         else
00805             return hmm_vit_eval_anytopo(hmm);
00806     }
00807 }
00808 
00809 int32
00810 hmm_dump_vit_eval(hmm_t * hmm, FILE * fp)
00811 {
00812     int32 bs = 0;
00813 
00814     if (fp) {
00815         fprintf(fp, "BEFORE:\n");
00816         hmm_dump(hmm, fp);
00817     }
00818     bs = hmm_vit_eval(hmm);
00819     if (fp) {
00820         fprintf(fp, "AFTER:\n");
00821         hmm_dump(hmm, fp);
00822     }
00823 
00824     return bs;
00825 }