viterbi.c

General Hidden Markov Model Library 一个通用的隐马尔科夫模型的C代码库

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/viterbi.c
*       Authors:  Wasinee Rungsarityotin, Benjamin Georgi
*
*       Copyright (C) 1998-2004 Alexander Schliep
*       Copyright (C) 1998-2001 ZAIK/ZPR, Universitaet zu Koeln
*       Copyright (C) 2002-2004 Max-Planck-Institut fuer Molekulare Genetik,
*                               Berlin
*
*       Contact: [e-mail]
*
*       This library is free software; you can redistribute it and/or
*       modify it under the terms of the GNU Library General Public
*       License as published by the Free Software Foundation; either
*       version 2 of the License, or (at your option) any later version.
*
*       This library is distributed in the hope that it will be useful,
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
*       Library General Public License for more details.
*
*       You should have received a copy of the GNU Library General Public
*       License along with this library; if not, write to the Free
*       Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*
*       This file is version $Revision: 1925 $
*                       from $Date: 2007-10-22 21:40:09 +0200 (Mon, 22 Oct 2007) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

#ifdef HAVE_CONFIG_H
#  include "../config.h"
#endif

#include <float.h>
#include <math.h>
#include <assert.h>

#include "ghmm.h"
#include "mes.h"
#include "mprintf.h"
#include "matrix.h"
#include "sdmodel.h"
#include "ghmm_internals.h"

typedef struct local_store_t {
    double **log_in_a;
    double **log_b;
    double *phi;
    double *phi_new;
    int **psi;

    int *topo_order;
    int topo_order_length;
} local_store_t;


/*----------------------------------------------------------------------------*/
static int viterbi_free(local_store_t **v, int n, int len)
{
#define CUR_PROC "viterbi_free"
    int j;
    mes_check_ptr(v, return -1);

    if (!*v)
        return 0;

    for (j = 0; j < n; j++)
        m_free((*v)->log_in_a[j]);

    m_free((*v)->log_in_a);
    ighmm_cmatrix_free(&((*v)->log_b), n);
    m_free((*v)->phi);
    m_free((*v)->phi_new);
    /*ighmm_dmatrix_free( &((*v)->psi), len ); */
    ighmm_dmatrix_stat_free(&((*v)->psi));
    m_free((*v)->topo_order);
    m_free(*v);

    return 0;
#undef CUR_PROC
}                               /* viterbi_free */

/*----------------------------------------------------------------------------*/
static local_store_t *viterbi_alloc(ghmm_dmodel *mo, int len)
{
#define CUR_PROC "sdviterbi_alloc"
    local_store_t *v = NULL;
    int j;
    ARRAY_CALLOC(v, 1);

    /* Allocate the log_in_a's -> individal lenghts */
    ARRAY_CALLOC(v->log_in_a, mo->N);
    for (j = 0; j < mo->N; j++) {
        ARRAY_CALLOC(v->log_in_a[j], mo->s[j].in_states);
    }

    v->log_b = ighmm_cmatrix_alloc(mo->N, len);
    if (!(v->log_b)) {
        GHMM_LOG_QUEUED(LCONVERTED);
        goto STOP;
    }
    ARRAY_CALLOC(v->phi, mo->N);
    ARRAY_CALLOC(v->phi_new, mo->N);
    v->psi = ighmm_dmatrix_stat_alloc(len, mo->N);
    if (!(v->psi)) {
        GHMM_LOG_QUEUED(LCONVERTED);
        goto STOP;
    }

    v->topo_order_length = 0;
    ARRAY_CALLOC(v->topo_order, mo->N);

    return v;
  STOP:                        /* Label STOP from ARRAY_[CM]ALLOC */
    viterbi_free(&v, mo->N, len);
    return NULL;
#undef CUR_PROC
}                               /* viterbi_alloc */

/*----------------------------------------------------------------------------*/
static void Viterbi_precompute(ghmm_dmodel *mo, int *o, int len, local_store_t *v)
{
#define CUR_PROC "viterbi_precompute"
    int i, j, t;

    /* Precomputing the log(a_ij) */
    for (j = 0; j < mo->N; j++)
        for (i = 0; i < mo->s[j].in_states; i++)
            if (mo->s[j].in_a[i] == 0.0)        /* DBL_EPSILON ? */
                v->log_in_a[j][i] = +1; /* Not used any further in the calculations */
            else
                v->log_in_a[j][i] = log(mo->s[j].in_a[i]);


    /* Precomputing the log(bj(ot)) */
    for (j = 0; j < mo->N; j++)
        for (t = 0; t < len; t++) {
            if (o[t] != mo->M) {
                if (mo->s[j].b[o[t]] == 0.0)    /* DBL_EPSILON ? */
                    v->log_b[j][t] = +1;
                else
                    v->log_b[j][t] = log(mo->s[j].b[o[t]]);
            } else {
                v->log_b[j][t] = 0.0;
            }
        }
#undef CUR_PROC
}                               /* viterbi_precompute */

/*----------------------------------------------------------------------------*/
static void viterbi_silent(ghmm_dmodel *mo, int t, local_store_t *v)
{
#define CUR_PROC "viterbi_silent"
    int topocount;
    int i, k;
    double max_value, value;

    for (topocount = 0; topocount < mo->topo_order_length; topocount++) {
        k = mo->topo_order[topocount];
        if (mo->silent[k]) {    /* Silent states */
            /* Determine the maximum */
            /* max_phi = phi[i] + log_in_a[j][i] ... */
            max_value = -DBL_MAX;
            v->psi[t][k] = -1;
            for (i = 0; i < mo->s[k].in_states; i++) {

                if (v->phi[mo->s[k].in_id[i]] != +1 && v->log_in_a[k][i] != +1) {
                    value = v->phi[mo->s[k].in_id[i]] + v->log_in_a[k][i];
                    if (value > max_value) {
                        max_value = value;
                        v->psi[t][k] = mo->s[k].in_id[i];
                    }
                }
            }
            /* No maximum found (that is, state never reached)
               or the output O[t] = 0.0: */
            if (max_value == -DBL_MAX) {
                v->phi[k] = +1;
            } else {
                v->phi[k] = max_value;
            }
        }
    }
#undef CUR_PROC
}

/*============================================================================*/
/* auxilary function. Reallocates an integer array to a given length
   and initialises the new positions with -1 */
static int extend_int_array(int *array, int cur_len, int extend)
{
#define CUR_PROC "extend_int_array"
    int j;
    cur_len += extend;
    ARRAY_REALLOC(array, cur_len);
    /* initialising new memory with -1  */
    for (j = cur_len - 1; j >= (cur_len - extend); j--) {
        array[j] = -1;
    }
    return cur_len;
  STOP:
    return -1;
#undef CUR_PROC
}

/*============================================================================*/
/** Return the viterbi path of the sequence.  */
int *ghmm_dmodel_viterbi(ghmm_dmodel * mo, int *o, int len, double *log_p)
{
#define CUR_PROC "ghmm_dmodel_viterbi"

    int *state_seq = NULL;
    int t, j, i, k, St;
    double value, max_value;
    local_store_t *v;

    /* length_factor determines the size of the memory allocation for the state
       path array in case the model contains silent states.
       The larger length_factor is chosen, the less reallocs will be necessary
       but it increases the amount of allocated memory that is not used. */
    int length_factor = 2;