foba.c

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

    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
}                               /* foba_label_initforward */




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

int ghmm_dmodel_label_forward (ghmm_dmodel * mo, const int *O, const int *label, int len,
                        double **alpha, double *scale, double *log_p)
{
# define CUR_PROC "ghmm_dl_forward"
  int res = -1;
  int i, t;
  int e_index;
  double c_t;
  char *str;

  foba_label_initforward (mo, alpha[0], O[0], label[0], scale);
  if (scale[0] < GHMM_EPS_PREC) {
    /* means: first symbol can't be generated by hmm */
    *log_p = +1;
  }
  else {
    *log_p = -log (1 / scale[0]);

    for (t = 1; t < len; t++) {

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

      /* printf("\n\nStep t=%i mit len=%i, O[i]=%i\n",t,len,O[t]); */

      for (i = 0; i < mo->N; i++) {
        if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {
          if (mo->label[i] == label[t]) {
	    /*printf("%d: akt_ state %d, label: %d \t current Label: %d\n",
	      t, i, mo->label[i], label[t]);*/
            e_index = get_emission_index (mo, i, O[t], t);
            if (-1 != e_index) {
              alpha[t][i] = ghmm_dmodel_forward_step(&mo->s[i], alpha[t-1], mo->s[i].b[e_index]);
              /*if (alpha[t][i] < GHMM_EPS_PREC) {
                printf("alpha[%d][%d] = %g \t ", t, i, alpha[t][i]);
		printf("mo->s[%d].b[%d] = %g\n", i, e_index, mo->s[i].b[e_index]);
		}
	      else printf("alpha[%d][%d] = %g\n", t, i, alpha[t][i]);*/
            }
            else {
              alpha[t][i] = 0;
            }
          }
          else {
            alpha[t][i] = 0;
          }
          scale[t] += alpha[t][i];
        }
        else {
          GHMM_LOG(LCONVERTED, "ERROR: Silent state in foba_label_forward.\n");
        }
      }

      if (scale[t] < GHMM_EPS_PREC) {
        if (t > 4) {
          str = ighmm_mprintf (NULL, 0, "%g\t%g\t%g\t%g\t%g\n", scale[t-5],
                         scale[t-4], scale[t-3], scale[t-2], scale[t-1]);
          GHMM_LOG(LCONVERTED, str);
          m_free (str);
        }
        str = ighmm_mprintf (NULL, 0, "scale = %g smaller than eps = EPS_PREC "
			     "in the %d-th char.\ncannot generate emission: %d "
			     "with label: %d in sequence of length %d\n",
			     scale[t], t, O[t], label[t], len);
        GHMM_LOG(LCONVERTED, str);
        m_free (str);
        /* O-string  can't be generated by hmm */
        *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);*/
      }
    }
  }

  /*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_label_logp (ghmm_dmodel * mo, const int *O, const int *label, int len,
                     double *log_p)
{
# define CUR_PROC "ghmm_dl_logp"
  int res = -1;
  double **alpha, *scale = NULL;

  alpha = ighmm_cmatrix_stat_alloc (len, mo->N);
  if (!alpha) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }
  ARRAY_CALLOC (scale, len);
  /* run ghmm_dmodel_forward */
  if (ghmm_dmodel_label_forward (mo, O, label, len, alpha, scale, log_p) == -1) {
    GHMM_LOG_QUEUED(LCONVERTED);
    goto STOP;
  }

  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  ighmm_cmatrix_stat_free (&alpha);
  m_free (scale);
  return (res);
# undef CUR_PROC
}                               /* ghmm_dl_logp */


/*============================================================================*/
int ghmm_dmodel_label_backward (ghmm_dmodel * mo, const int *O, const int *label, int len,
                         double **beta, double *scale, double *log_p)
{
# define CUR_PROC "ghmm_dl_backward"

  double *beta_tmp, sum;
  int i, j, j_id, t;
  int res = -1;
  int e_index;
  /* int beta_out=0; */
  double emission;

  ARRAY_CALLOC (beta_tmp, mo->N);
  for (t = 0; t < len; t++)
    mes_check_0 (scale[t], goto STOP);

  /* check for silent states */
  if (mo->model_type & GHMM_kSilentStates) {
    GHMM_LOG(LCONVERTED, "ERROR: No silent states allowed in labelled HMM!\n");
    goto STOP;
  }

  /* initialize */
  for (i = 0; i < mo->N; i++) {
    /* start only in states with the correct label */
    if (label[len - 1] == mo->label[i])
      beta[len - 1][i] = 1.0;
    else
      beta[len - 1][i] = 0.0;

    beta_tmp[i] = beta[len - 1][i] / scale[len - 1];
  }

  /* initialize emission history */
  if (!(mo->model_type & GHMM_kHigherOrderEmissions))
    mo->maxorder = 0;
  for (t = len - (mo->maxorder); t < len; t++) {
    update_emission_history (mo, O[t]);
  }

  /* Backward Step for t = T-1, ..., 0
     beta_tmp: Vector for storage of scaled beta in one time step
     loop over reverse topological ordering of silent states, non-silent states */
  for (t = len - 2; t >= 0; t--) {

    /* 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]);

    for (i = 0; i < mo->N; i++) {
      sum = 0.0;
      for (j = 0; j < mo->s[i].out_states; j++) {
        j_id = mo->s[i].out_id[j];
        /* The state has only a emission with probability > 0, if the label matches */
        if (label[t] == mo->label[i]) {
          e_index = get_emission_index (mo, j_id, O[t + 1], t + 1);
          if (e_index != -1)
            emission = mo->s[j_id].b[e_index];
          else
            emission = 0.0;
        }
        else
          emission = 0.0;

        sum += mo->s[i].out_a[j] * emission * beta_tmp[j_id];
      }
      beta[t][i] = sum;
      /* if ((beta[t][i] > 0) && ((beta[t][i] < .01) || (beta[t][i] > 100))) beta_out++; */
    }

    for (i = 0; i < mo->N; i++)
      beta_tmp[i] = beta[t][i] / scale[t];
  }

  /* printf("labeled betas out of [.01 100]: %d (%f %)\n", beta_out, (float)beta_out/(mo->N*len)); */

  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  m_free (beta_tmp);
  return (res);
# undef CUR_PROC
}


/*===========  Lean forward algorithm for labeled states  ====================*/
int ghmm_dl_forward_lean (ghmm_dmodel * mo, const int *O, const int *label,
                             int len, double *log_p)
{
# define CUR_PROC "ghmm_dl_forward_lean"
  int res = -1;
  int i, t, id, e_index;
  double c_t;
  double log_scale_sum = 0.0;
  double non_silent_salpha_sum = 0.0;
  double salpha_log = 0.0;

  double *alpha_last_col=NULL;
  double *alpha_curr_col=NULL;
  double *switching_tmp=NULL;
  double *scale=NULL;

  /* Allocating */
  ARRAY_CALLOC (alpha_last_col, mo->N);
  ARRAY_CALLOC (alpha_curr_col, mo->N);
  ARRAY_CALLOC (scale, len);

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

  foba_label_initforward (mo, alpha_last_col, O[0], label[0], scale);
  if (scale[0] < GHMM_EPS_PREC) {
    /* means: first symbol can't be generated by hmm */
    *log_p = +1;
    goto STOP;
  }

  *log_p = -log (1 / scale[0]);

  for (t = 1; t < len; t++) {
    scale[t] = 0.0;

    /* iterate over non-silent states */
    for (i = 0; i < mo->N; i++) {
      if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {

        /* printf("  akt_ state %d\n",i);*/
        if (mo->label[i] == label[t]) {
          e_index = get_emission_index (mo, i, O[t], t);
          if (e_index != -1) {
            alpha_curr_col[i] = ghmm_dmodel_forward_step (&mo->s[i], alpha_last_col,
                                                  mo->s[i].b[e_index]);
            scale[t] += alpha_curr_col[i];
          }
          else
            alpha_curr_col[i] = 0;
        }
        else
          alpha_curr_col[i] = 0;
      }
    }
    /* iterate over silent states  */
    if (mo->model_type & GHMM_kSilentStates) {
      for (i = 0; i < mo->topo_order_length; i++) {
        id = mo->topo_order[i];
        alpha_curr_col[id] = ghmm_dmodel_forward_step (&mo->s[id], alpha_last_col, 1);
        scale[t] += alpha_curr_col[id];
      }
    }

    if (scale[t] < GHMM_EPS_PREC) {
      GHMM_LOG(LCONVERTED, "scale smaller than epsilon\n");
      /* O-string  can't be generated by hmm */
      *log_p = +1.0;
      break;
    }
    c_t = 1 / scale[t];
    for (i = 0; i < mo->N; i++)
      alpha_curr_col[i] *= c_t;

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

    /* switching pointers of alpha_curr_col and alpha_last_col
       don't set alpha_curr_col[i] to zero since its overwritten */
    switching_tmp = alpha_last_col;
    alpha_last_col = alpha_curr_col;
    alpha_curr_col = switching_tmp;
  }

  /* Termination step: compute log likelihood */
  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_curr_col[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;

STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  /* Deallocation */
  m_free (alpha_last_col);
  m_free (alpha_curr_col);
  m_free (scale);
  return res;
#undef CUR_PROC
}                               /* foba_forward_label_lean */