pviterbi.c

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

  for (u = mo->max_offset_x + 1; u < X->length; u++) {
    for (v = -mo->max_offset_y; v < Y->length; v++) {
      for (j = 0; j < mo->N; j++) 
	{
	  /** initialization of phi (lookback matrix), psi (traceback) **/
	  set_phi(pv, u, v, j, +1);
	  set_psi(pv, u, v, j, -1);
	}
      
      for (i = 0; i < mo->N; i++) {
	/* Determine the maximum */
	/* max_phi = phi[i] + log_in_a[j][i] ... */
	if (!(mo->model_type & GHMM_kSilentStates) || !mo->silent[i] ) {
	  max_value = -DBL_MAX;
	  set_psi(pv, u, v, i, -1);
	  for (j = 0; j < mo->s[i].in_states; j++) {
	    /* look back in the phi matrix at the offsets */
	    previous_prob = get_phi(pv, u, v, mo->s[i].offset_x, mo->s[i].offset_y, mo->s[i].in_id[j]);
	    log_in_a_ij = (*log_in_a)(pv, i, j, X, Y, u, v);
	    if ( previous_prob != +1 && log_in_a_ij != +1) {
	      value = previous_prob + log_in_a_ij;
	      if (value > max_value) {
		max_value = value;
		set_psi(pv, u, v, i, mo->s[i].in_id[j]);
	      }
	    }
	    else
	      {;} /* fprintf(stderr, " %d --> %d = %f, \n", i,i,v->log_in_a[i][i]); */
	  }

	  emission = ghmm_dpmodel_pair(ghmm_dpseq_get_char(X, mo->s[i].alphabet, u), 
			      ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v),
			      mo->size_of_alphabet[mo->s[i].alphabet],
			      mo->s[i].offset_x, mo->s[i].offset_y);
#ifdef DEBUG
	  if (emission > ghmm_dpmodel_emission_table_size(mo, i)){
	    printf("State %i\n", i);
	    ghmm_dpmodel_state_print(&(mo->s[i]));
	    printf("charX: %i charY: %i alphabet size: %i emission table: %i emission index: %i\n", 
		   ghmm_dpseq_get_char(X, mo->s[i].alphabet, u),
		   ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v),
		   mo->size_of_alphabet[mo->s[i].alphabet],
		   ghmm_dpmodel_emission_table_size(mo, i), emission);
	  }
#endif
	  log_b_i = log_b(pv, i, ghmm_dpmodel_pair(ghmm_dpseq_get_char(X, mo->s[i].alphabet, u), 
				      ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v),
				      mo->size_of_alphabet[mo->s[i].alphabet],
				      mo->s[i].offset_x, mo->s[i].offset_y));

	  /* No maximum found (that is, state never reached)
	     or the output O[t] = 0.0: */
	  if (max_value == -DBL_MAX ||/* and then also: (v->psi[t][j] == -1) */
	      log_b_i == +1 ) {
	    set_phi(pv, u, v, i, +1);
	  }
	  else
	    set_phi(pv, u, v, i, max_value + log_b_i);
	}
      } /* complete time step for emitting states */
    
	/* last_osc = osc; */ 
        /* save last transition class */

      /*if (mo->model_type & kSilentStates) { 
	p__viterbi_silent( mo, t, v );
	}*/ /* complete time step for silent states */
      
      /**************
    for (j = 0; j < mo->N; j++) 
      {      
	printf("\npsi[%d],in:%d, phi=%f\n", t, v->psi[t][j], v->phi[j]);
       }
      
    for (i = 0; i < mo->N; i++){
      printf("%d\t", former_matchcount[i]);
    }

    for (i = 0; i < mo->N; i++){
      printf("%d\t", recent_matchcount[i]);
    }
      ****************/
    } /* End for v in Y */
    /* Next character in X */
    push_back_phi(pv, Y->length);
  } /* End for u in X */

  /* Termination */
  max_value = -DBL_MAX;
  ighmm_list_append(state_list, -1);
  /* if start_traceback_with is -1 (it is by default) search for the most 
     likely state at the end of both sequences */
  if (start_traceback_with == -1) {
    for (j = 0; j < mo->N; j++){
#ifdef DEBUG
      printf("phi(len_x)(len_y)(%i)=%f\n", j, get_phi(pv, u, Y->length-1, 0, 0, j));
#endif
      if ( get_phi(pv, u, Y->length-1, 0, 0, j) != +1 && 
	   get_phi(pv, u, Y->length-1, 0, 0, j) > max_value) { 
	max_value = get_phi(pv, X->length-1, Y->length-1, 0, 0, j);
	state_list->last->val = j;
      }
    }
  }
  /* this is the special traceback mode for the d & c algorithm that also 
     connects the traceback to the first state of the rest of the path */
  else {
#ifdef DEBUG
    printf("D & C traceback from state %i!\n", start_traceback_with);
    printf("Last characters emitted X: %i, Y: %i\n", 
	   ghmm_dpseq_get_char(X, mo->s[start_traceback_with].alphabet, 
			       X->length-1),
	   ghmm_dpseq_get_char(Y, mo->s[start_traceback_with].alphabet, 
			       Y->length-1));
    for (j = 0; j < mo->N; j++){
      printf("phi(len_x)(len_y)(%i)=%f\n", j, get_phi(pv, X->length-1, Y->length-1, 0, 0, j)); 
    }
#endif
    max_value = get_phi(pv, X->length-1, Y->length-1, 0, 0, start_traceback_with);
    if (max_value != 1 && max_value > -DBL_MAX)
      state_list->last->val = start_traceback_with;
  }
  if (max_value == -DBL_MAX) {
    /* Sequence can't be generated from the model! */
    *log_p = +1;
    /* Backtracing doesn't work, because state_seq[*] allocated with -1 */
    /* for (t = len - 2; t >= 0; t--)
       state_list->last->val = -1;    */
  }
  else {
    /* Backtracing, should put DEL path nicely */
    *log_p = max_value;
    /* removed the handling of silent states here */
    /* start trace back at the end of both sequences */
    u = X->length - 1;
    v = Y->length - 1;
    current_state_index = state_list->first->val;
    off_x = mo->s[current_state_index].offset_x;
    off_y = mo->s[current_state_index].offset_y;
    while (u - off_x >= -1 && v - off_y >= -1 && current_state_index != -1) { 
      /* while (u > 0 && v > 0) { */
      /* look up the preceding state and save it in the first position of the
	 state list */
      /* printf("Current state %i at (%i,%i) -> preceding state %i\n", 
	 current_state_index, u, v, get_psi(pv, u, v, current_state_index)); */
      /* update the current state */
      current_state_index = get_psi(pv, u, v, current_state_index);
      if (current_state_index != -1)
	ighmm_list_insert(state_list, current_state_index);
      /* move in the alignment matrix */
      u -= off_x;
      v -= off_y; 
      /* get the next offsets */
      off_x = mo->s[current_state_index].offset_x;
      off_y = mo->s[current_state_index].offset_y;
    }
  }
  
  /* Free the memory space */
  pviterbi_free(&pv, mo->N, X->length, Y->length, mo->max_offset_x , 
		mo->max_offset_y);
  /* printf("After traceback: last state = %i\n", state_list->last->val); */
  state_seq = ighmm_list_to_array(state_list);
  *path_length = state_list->length;
  /* PRINT PATH */
  
/*   fprintf(stderr, "Viterbi path: " ); */
/*   int t; */
/*   for(t=0; t < *path_length; t++) */
/*     if (state_seq[t] >= 0) fprintf(stderr, " %d ",  state_seq[t]); */
/*   fprintf(stderr, "\n Freeing ... \n");  */
  return (state_seq);
STOP:     /* Label STOP from ARRAY_[CM]ALLOC */
  /* Free the memory space */
  pviterbi_free(&pv, mo->N, X->length, Y->length, mo->max_offset_x, 
		mo->max_offset_y);
  m_free(state_seq);
  ighmm_list_free(state_list);
  return NULL;
#undef CUR_PROC
} /* viterbi */
  
  
/*============================================================================*/
double ghmm_dpmodel_viterbi_logp(ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y,
			    int *state_seq, int state_seq_len) {
#define CUR_PROC "ghmm_dpmodel_viterbi_logp"
  int s, t, i, j, u, v;
  double log_p = 0.0;
  double log_b_i = 1.0;
  double log_in_a = 1.0;
  plocal_store_t *pv;

  /* Allocate the matrices log_in_a, log_b,Vektor phi, phi_new, Matrix psi */
  pv = pviterbi_alloc(mo, 0, 0);
  pviterbi_precompute(mo, pv);
  /* initial state */
  t = 0; u = -1; v = -1;
  if (state_seq_len > t) {
    i = state_seq[t];
    /* initial state probability */
    /* log_p += log(mo->s[i].pi); */
    log_p += mo->s[i].log_pi;
    if (log_p == 1.0) {
      pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
      fprintf(stderr, "the initial probability of state %i is zero\n", i);
      return 1.0;/* the initial prob is zero */
    }
    /* consume the first characters of the sequences */
    u += mo->s[i].offset_x;
    v += mo->s[i].offset_y;
    /* get the emission probability */
    log_b_i = log_b(pv, i, ghmm_dpmodel_pair(ghmm_dpseq_get_char(X, mo->s[i].alphabet, u), 
				ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v),
				mo->size_of_alphabet[mo->s[i].alphabet],
				mo->s[i].offset_x, mo->s[i].offset_y));
    if (log_b_i == 1.0) { /* chars cant be emitted */
      pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
      fprintf(stderr, "characters (%i, %i) at position (%i, %i) cannot be emitted by state %i (t=%i)\n",
	      ghmm_dpseq_get_char(X, mo->s[i].alphabet, u),
	      ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v), u, v, i, t);
      return 1.0;
    }
    log_p += log_b_i;
  }
  else { /* there is no path.. */
    pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
    fprintf(stderr, "No path given!\n");
    return 1.0;
  }
  /* rest of the sequence */
  for (t=1; t<state_seq_len; t++) {
    j = i; /* predecessor state */
    i = state_seq[t]; /* current state */
    /* consume the next characters */
    u += mo->s[i].offset_x;
    v += mo->s[i].offset_y;
    if (u >= X->length || v >= Y->length) { /* path consumes too many chars */
      pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
      fprintf(stderr, "path consumes too many chars\n");
      return 1.0;
    }
    /* transition prob state j -> i*/
    log_in_a = 1.0;
    for (s=0; s<mo->s[i].in_states; s++) {
      if (mo->s[i].in_id[s] == j) {
	log_in_a = sget_log_in_a(pv, i, s, X, Y, u, v);
	break;
      }
    }
    if (log_in_a == 1.0) {
      pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
      fprintf(stderr, "transition (%i -> %i) at t=%i not possible\n", j, i,t); 
      return 1.0; /* transition not possible */
    }
    /* emission probability */
    log_b_i = log_b(pv, i, ghmm_dpmodel_pair(ghmm_dpseq_get_char(X, mo->s[i].alphabet, u), 
				ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v),
				mo->size_of_alphabet[mo->s[i].alphabet],
				mo->s[i].offset_x, mo->s[i].offset_y));
    if (log_b_i == 1.0) {
      pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, mo->max_offset_y);
      fprintf(stderr, "characters (%i, %i) at position (%i, %i) cannot be emitted by state %i (t=%i)\n",
	      ghmm_dpseq_get_char(X, mo->s[i].alphabet, u),
	      ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v), u, v, i, t);
      return 1.0; /* characters cant be emitted */
    }
    log_p += log_in_a + log_b_i;
  }
  pviterbi_free(&pv, mo->N, 0, 0, mo->max_offset_x, 
		mo->max_offset_y);
  /* check if all of the sequences has been consumed */
  if (u != X->length - 1 && v != Y->length - 1) {
    fprintf(stderr, "path consumes not all characters (%i of %i, %i of %i)\n",
	    u + 1, X->length, v + 1, Y->length);
    return 1.0;
  }
  return log_p;

#undef CUR_PROC
} /* ghmm_dpmodel_viterbi_logp */