foba.c

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

	   multiple scaling with the same scalingfactor in one term */
        beta_tmp[id] = sum;
      }
    }

    /* iterating over the the non-silent states */
    for (i = 0; i < mo->N; i++) {
      if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {
        sum = 0.0;
	
        for (j = 0; j < mo->s[i].out_states; j++) {
	  j_id = mo->s[i].out_id[j];
        
	  /* out_state is not silent: get the emission probability
	     and use beta[t+1]*/
	  if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[j_id])) {
            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;
	    sum += mo->s[i].out_a[j] * emission * beta[t+1][j_id];
          }
	  /* out_state is silent: use beta_tmp */
          else {
            sum += mo->s[i].out_a[j] * beta_tmp[j_id];
          }
        }
	/* updating beta[t] for non-silent state */
        beta[t][i] = sum / scale[t+1];
      }
    }
    /* updating beta[t] for silent states, finally scale them
       and resetting beta_tmp */
    if (mo->model_type & GHMM_kSilentStates)
      for (i=0; i<mo->N; i++) {
	if (mo->silent[i]) {
	  beta[t][i] = beta_tmp[i] / scale[t+1];
	  beta_tmp[i] = 0.0;
	}
      }
  }

  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  if (mo->model_type & GHMM_kSilentStates) m_free (beta_tmp);
  return (res);
# undef CUR_PROC
}                               /* ghmm_dmodel_backward */


/*============================================================================*/
int ghmm_dmodel_backward_termination (ghmm_dmodel *mo, const int *O, int length,
			       double **beta, double *scale, double *log_p) {
#define CUR_PROC "ghmm_dmodel_backward_termination"
  int i, id, j, j_id, k, res=-1;
  double log_scale_sum, sum;
  double * beta_tmp=NULL;

  /* topological ordering for models with silent states and precomputing
     the beta_tmp for silent states */
  if (mo->model_type & GHMM_kSilentStates) {
    ghmm_dmodel_order_topological(mo);

    ARRAY_CALLOC (beta_tmp, mo->N);
    for (k = mo->topo_order_length - 1; k >= 0; k--) {
      id = mo->topo_order[k];
      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]) {
	  /* no emission history for the first symbol */
	  if (!(mo->model_type & GHMM_kHigherOrderEmissions) || mo->order[id]==0) {
	    sum += mo->s[id].out_a[j] * mo->s[j_id].b[O[0]] * beta[0][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 */
      beta_tmp[id] = sum;
    }
  }

  sum = 0.0;
  /* iterating over all states with pi > 0.0 */
  for (i = 0; i < mo->N; i++) {
    if (mo->s[i].pi > 0.0) {
      /* silent states */
      if ((mo->model_type & GHMM_kSilentStates) && mo->silent[i]) {
	sum += mo->s[i].pi * beta_tmp[i];
      }
      /* non-silent states */
      else {
	/* no emission history for the first symbol */
	if (!(mo->model_type & GHMM_kHigherOrderEmissions) || mo->order[i]==0) {
	  sum += mo->s[i].pi * mo->s[i].b[O[0]] * beta[0][i];
	}
      }
    }
  }
  
  *log_p = log (sum / scale[0]);

  log_scale_sum = 0.0;
  for (i = 0; i < length; i++) {
    log_scale_sum += log (scale[i]);
  }

  *log_p += log_scale_sum;
  
  res = 0;
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  if (mo->model_type & GHMM_kSilentStates) m_free (beta_tmp);
  return (res);
# undef CUR_PROC
}


/*============================================================================*/
int ghmm_dmodel_logp (ghmm_dmodel * mo, const int *O, int len, double *log_p)
{
# define CUR_PROC "ghmm_dmodel_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_forward (mo, O, 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_dmodel_logp */

/*============================================================================*/
  int ghmm_dmodel_logp_joint(ghmm_dmodel * mo, const int *O, int len,
                            const int *S, int slen, double *log_p)
{
# define CUR_PROC "ghmm_dmodel_logp_joint"
    int prevstate, state, spos=0, pos=0, j;

    prevstate = state = S[0];
    *log_p = log(mo->s[state].pi);
    if (!(mo->model_type & GHMM_kSilentStates) || !mo->silent[state])
        *log_p += log(mo->s[state].b[O[pos++]]);
        
    for (spos=1; spos < slen || pos < len; spos++) {
        state = S[spos];
        for (j=0; j < mo->s[state].in_states; ++j) {
            if (prevstate == mo->s[state].in_id[j])
                break;
        }

        if (j == mo->s[state].in_states ||
            fabs(mo->s[state].in_a[j]) < GHMM_EPS_PREC) {
            GHMM_LOG_PRINTF(LERROR, LOC, "Sequence can't be built. There is no "
                            "transition from state %d to %d.", prevstate, state);
            goto STOP;
        }

        *log_p += log(mo->s[state].in_a[j]);

        if (!(mo->model_type & GHMM_kSilentStates) || !mo->silent[state]) {
            *log_p += log(mo->s[state].b[O[pos++]]);
        }
        
        prevstate = state;
    }

    if (pos < len)
        GHMM_LOG_PRINTF(LINFO, LOC, "state sequence too short! processed only %d symbols", pos);
    if (spos < slen)
        GHMM_LOG_PRINTF(LINFO, LOC, "sequence too short! visited only %d states", spos);

    return 0;
  STOP:
    return -1;
# undef CUR_PROC
}

/*============================================================================*/
/*====================== Lean forward algorithm  ====================*/
int ghmm_dmodel_forward_lean (ghmm_dmodel * mo, const int *O, int len, double *log_p)
{
# define CUR_PROC "ghmm_dmodel_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;
  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);

  ghmm_dmodel_forward_init (mo, alpha_last_col, O[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++) {
      scale[t] = 0.0;
      update_emission_history (mo, O[t - 1]);

      /* 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_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;
        }
      }
      /* 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_curr_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++) {
	/* use alpha_last_col since the columms are also in the last step
	   switched */
        if (!(mo->silent[i]))
          non_silent_salpha_sum += alpha_last_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
}                               /* ghmm_dmodel_forward_lean */


/*======================= labeled HMMs =======================================*/
static int foba_label_initforward (ghmm_dmodel * mo, double *alpha_1, int symb,
                                   int label, double *scale)
{
# define CUR_PROC "foba_label_initforward"
  int res = -1;
  int i;
  double c_0;
  scale[0] = 0.0;

  /* iterate over non-silent states */
  for (i = 0; i < mo->N; i++) {
    if (!(mo->model_type & GHMM_kSilentStates) || !(mo->silent[i])) {
      if (mo->label[i] == label)
        alpha_1[i] = mo->s[i].pi * mo->s[i].b[symb];
      else
        alpha_1[i] = 0.0;
    }
    else
      alpha_1[i] = 0.0;

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

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

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

  if (scale[0] >= GHMM_EPS_PREC) {