foba.c

来自「General Hidden Markov Model Library 一个通用」· C语言 代码 · 共 1,035 行 · 第 1/3 页

/*******************************************************************************
*
*       This file is part of the General Hidden Markov Model Library,
*       GHMM version 0.8_beta1, see http://ghmm.org
*
*       Filename: ghmm/ghmm/foba.c
*       Authors:  Utz J. Pape, Benjamin Georgi, Janne Grunau
*
*       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: 1932 $
*                       from $Date: 2007-11-01 19:24:29 +0100 (Thu, 01 Nov 2007) $
*             last change by $Author: grunau $.
*
*******************************************************************************/

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

#include <assert.h>

#include "ghmm.h"
#include "model.h"
#include "math.h"
#include "matrix.h"
#include "mes.h"
#include "mprintf.h"
#include "foba.h"
#include "ghmm_internals.h"

int ghmm_dmodel_forward_init (ghmm_dmodel * mo, double *alpha_1, int symb, double *scale)
{
# define CUR_PROC "ghmm_dmodel_forward_init"
  int res = -1;
  int i, j, id, in_id;
  double c_0;
  scale[0] = 0.0;

  /*printf(" *** foba_initforward\n");*/

  /*iterate over non-silent states*/
  /*printf(" *** iterate over non-silent states \n");*/
  for (i = 0; i < mo->N; i++) {
    if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {
      /*no starting in states with order > 0 !!!*/
      if (!(mo->model_type & GHMM_kHigherOrderEmissions) || mo->order[i] == 0) {
        alpha_1[i] = mo->s[i].pi * mo->s[i].b[symb];
        /*printf("\nalpha1[%i]=%f\n",i,alpha_1[i]);*/

        scale[0] += alpha_1[i];
      }
      else {
        alpha_1[i] = 0;
      }
    }
  }
  /*iterate over silent states*/
  /*printf(" *** iterate over silent states \n");*/
  if (mo->model_type & GHMM_kSilentStates) {
    for (i = 0; i < mo->topo_order_length; i++) {
      id = mo->topo_order[i];
      alpha_1[id] = mo->s[id].pi;

      /*printf("\nsilent_start alpha1[%i]=%f\n",id,alpha_1[id]);*/

      for (j = 0; j < mo->s[id].in_states; j++) {
        in_id = mo->s[id].in_id[j];
        alpha_1[id] += mo->s[id].in_a[j] * alpha_1[in_id];

        /*printf("\n\tsilent_run alpha1[%i]=%f\n",id,alpha_1[id]);*/

      }
      scale[0] += alpha_1[id];
    }
  }

  /* printf("\nwo label: scale[0] = %f\n",scale[0]); */

  if (scale[0] >= GHMM_EPS_PREC) {
    c_0 = 1 / scale[0];
    for (i = 0; i < mo->N; i++)
      alpha_1[i] *= c_0;
  }
  res = 0;
  return (0);                   /* attention: scale[0] might be 0 */
# undef CUR_PROC
}                               /* ghmm_dmodel_forward_init */

/*----------------------------------------------------------------------------*/

/** modified by Casillux to improve performance */
double ghmm_dmodel_forward_step (ghmm_dstate * s, double *alpha_t, const double b_symb)
{
  int i, id;
  double value = 0.0;

  if (b_symb < GHMM_EPS_PREC)
    return 0.;

  /*printf(" *** foba_stepforward\n");*/

  for (i = 0; i < s->in_states; i++) {
    id = s->in_id[i];
    value += s->in_a[i] * alpha_t[id];
    /*printf("    state %d, value %f, p_symb %f\n",id, value, b_symb); */
  }
  value *= b_symb;
  return (value);

}                               /* ghmm_dmodel_forward_step */

/*============================================================================*/

int ghmm_dmodel_forward (ghmm_dmodel * mo, const int *O, int len, double **alpha,
                  double *scale, double *log_p)
{
# define CUR_PROC "ghmm_dmodel_forward"
  char * str;
  int res = -1;
  int i, t, id;
  int e_index;
  double c_t;
  double log_scale_sum = 0.0;
  double non_silent_salpha_sum = 0.0;
  double salpha_log = 0.0;


  if (mo->model_type & GHMM_kSilentStates)
    ghmm_dmodel_order_topological(mo);

  ghmm_dmodel_forward_init (mo, alpha[0], O[0], scale);

  if (scale[0] < GHMM_EPS_PREC) {
    /* means: first symbol can't be generated by hmm */
    /*printf("\nscale kleiner als eps (line_no: 123)\n");*/
    *log_p = +1;
  }
  else {
    *log_p = -log (1 / scale[0]);
    for (t = 1; t < len; t++) {

      scale[t] = 0.0;
      update_emission_history (mo, O[t - 1]);

      /* printf("\n\nStep t=%i mit len=%i, O[i]=%i\n",t,len,O[t]);
         printf("iterate over non-silent state\n"); */
      /* iterate over non-silent states */
      for (i = 0; i < mo->N; i++) {
        if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {
          e_index = get_emission_index (mo, i, O[t], t);
          if (e_index != -1) {
            alpha[t][i] =
              ghmm_dmodel_forward_step (&mo->s[i], alpha[t - 1], mo->s[i].b[e_index]);
            scale[t] += alpha[t][i];
          }
          else {
            alpha[t][i] = 0;
          }
        }
      }
      /* iterate over silent states */
      /* printf("iterate over silent state\n"); */
      if (mo->model_type & GHMM_kSilentStates) {
        for (i = 0; i < mo->topo_order_length; i++) {
          /*printf("\nget id\n");*/
          id = mo->topo_order[i];
          /*printf("  akt_ state %d\n",id);*/
          /*printf("\nin stepforward\n");*/
          alpha[t][id] = ghmm_dmodel_forward_step (&mo->s[id], alpha[t], 1);
          /*printf("\nnach stepforward\n");*/
          scale[t] += alpha[t][id];
        }
      }

      if (scale[t] < GHMM_EPS_PREC) {
        /* O-string  can't be generated by hmm */
        str = ighmm_mprintf(NULL, 0, "scale smaller than epsilon (%g < %g) in position %d. Can't generate symbol %d\n", scale[t], GHMM_EPS_PREC, t, O[t]);
	GHMM_LOG(LCONVERTED, str);
	m_free (str);
        *log_p = +1.0;
        break;
      }
      c_t = 1 / scale[t];
      for (i = 0; i < mo->N; i++) {
        alpha[t][i] *= c_t;
      }

      if (!(mo->model_type & GHMM_kSilentStates) && *log_p != +1) {
        /* sum log(c[t]) scaling values to get  log( P(O|lambda) ) */

        /*printf("log_p %f -= log(%f) = ",*log_p,c_t);*/
        *log_p -= log (c_t);
        /*printf(" %f\n",*log_p); */
      }
    }
    if (mo->model_type & GHMM_kSilentStates && *log_p != +1) {
      /*printf("silent model\n");*/
      for (i = 0; i < len; i++) {
        log_scale_sum += log (scale[i]);
      }
      for (i = 0; i < mo->N; i++) {
        if (!(mo->silent[i])) {
          non_silent_salpha_sum += alpha[len - 1][i];
        }
      }
      salpha_log = log (non_silent_salpha_sum);
      *log_p = log_scale_sum + salpha_log;
    }
  }

  /*printf("\nin forward: log_p = %f\n",*log_p);*/

  if (*log_p == 1.0) {
    res = -1;
  }
  else {
    res = 0;
  }

  return res;
# undef CUR_PROC
}                               /* ghmm_dmodel_forward */

/*============================================================================*/

int ghmm_dmodel_forward_descale (double **alpha, double *scale, int t, int n,
                  double **newalpha)
{
# define CUR_PROC "ghmm_dmodel_forward_descale"
  int i, j, k;
  /*printf("\nAngekommen, t=%i, n=%i\n",t,n);*/
  for (i = 0; i < t; i++) {
    /*printf("i=%i\n",i);*/
    for (j = 0; j < n; j++) {
      /*printf("\tj=%i\n",j);*/
      newalpha[i][j] = alpha[i][j];
      /*newalpha[i][j] *= scale[j];*/
      /*printf(".");*/
      for (k = 0; k <= i; k++) {
        /*printf(",");*/
        newalpha[i][j] *= scale[k];
      }
    }
  }
  /*printf("\ndescale geschafft\n");*/
  return 0;
# undef CUR_PROC
}                               /* ghmm_dmodel_forward_descale */


/***************************** Backward Algorithm ******************************/
int ghmm_dmodel_backward (ghmm_dmodel * mo, const int *O, int len, double **beta,
                   const double *scale)
{
# define CUR_PROC "ghmm_dmodel_backward"
  /* beta_tmp holds beta-variables for silent states */
  double *beta_tmp=NULL;
  double sum, emission;
  int i, j, j_id, t, k, id;
  int res = -1;
  int e_index;


  for (t = 0; t < len; t++)
    mes_check_0 (scale[t], goto STOP);

  /* topological ordering for models with silent states and allocating
     temporary array needed for silent states */
  if (mo->model_type & GHMM_kSilentStates) {
    ARRAY_CALLOC (beta_tmp, mo->N);
    ghmm_dmodel_order_topological(mo);
  }

  /* initialize all states */
  for (i = 0; i < mo->N; i++) {
    beta[len - 1][i] = 1.0;
  }

  if (!(mo->model_type & GHMM_kHigherOrderEmissions)) {
    mo->maxorder = 0;
  }
  /* initialize emission history */
  for (t = len - (mo->maxorder); t < len; t++) {
    update_emission_history (mo, O[t]);
  }
  
  /* Backward Step for t = T-1, ..., 0 */
  /* loop over reverse topological ordering of silent states, non-silent states  */
  for (t = len - 2; t >= 0; t--) {
    /* printf(" ----------- *** t = %d ***  ---------- \n",t); */
    /* printf("\n*** O(%d) = %d\n",t+1,O[t+1]); */

    /* updating of emission_history with O[t] such that emission_history memorizes
       O[t - maxorder ... t] */
    if (0 <= t - mo->maxorder + 1)
      update_emission_history_front (mo, O[t - mo->maxorder + 1]);

    /* iterating over the the silent states and filling beta_tmp */
    if (mo->model_type & GHMM_kSilentStates) {
      for (k = mo->topo_order_length - 1; k >= 0; k--) {
        id = mo->topo_order[k];
        /* printf("  silent[%d] = %d\n",id,mo->silent[id]); */
        assert (mo->silent[id] == 1);

        sum = 0.0;
        for (j = 0; j < mo->s[id].out_states; j++) {
          j_id = mo->s[id].out_id[j];

	  /* out_state is not silent */
	  if (!mo->silent[j_id]) {
	    e_index = get_emission_index (mo, j_id, O[t + 1], t + 1);
	    if (e_index != -1) {
	      sum += mo->s[id].out_a[j] * mo->s[j_id].b[e_index] * beta[t+1][j_id];
	    }
	  }
	  /* out_state is silent, beta_tmp[j_id] is useful since we go through
	     the silent states in reversed topological order */
	  else {
	    sum += mo->s[id].out_a[j] * beta_tmp[j_id];
	  }
	}
	/* setting beta_tmp for the silent state
	   don't scale the betas for silent states now 
	   wait until the betas for non-silent states are complete to avoid