pviterbi.c

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

#endif
  return pv->log_b[state][emission];
}

/*============================================================================*/
static void init_phi(plocal_store_t * pv, ghmm_dpseq * X, ghmm_dpseq * Y) {
#ifdef DEBUG
  int emission;
#endif
  int u, v, j, i, off_x, y;
  double log_in_a_ij;
  double value, max_value, previous_prob, log_b_i;  
  /* printf("ghmm_dpmodel_viterbi init\n"); */
  ghmm_dpmodel * mo = pv->mo;
  double (*log_in_a)(plocal_store_t*, int, int, ghmm_dpseq*, ghmm_dpseq*, 
		     int, int);
  log_in_a = &sget_log_in_a;

  /* Initialize the lookback matrix (for positions [-offsetX,0], [0, len_y]*/
  for (off_x=0; off_x<mo->max_offset_x + 1; off_x++)
    for (y=0; y<Y->length + mo->max_offset_y + 1; y++)
      for (j=0; j<mo->N; j++) {
	pv->phi[off_x][y][j] = +1;
      }
    if ( mo->model_type & GHMM_kSilentStates ) { /* could go into silent state at t=0 */

    /*p__viterbi_silent( mo, t=0, v);*/
  }
  /*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]);
  }*/
  
  /* initialize for offsets > 1 (u < max_offset_x, v < max_offset_y) */
  /* this is in principle the same as the main recurrence but adds initial
     probabilities to states that cannot be inhabitated at u=0, v=0 because
     of greater offsets than one 
     iteration start is u=-1 v=-1 to allow states with offset_x == 0 
     which corresponds to a series of gap states before reading the first 
     character of x at position x=0, y=v or equally for offset_y == 0 */
  /* u, v <= max offsets */
    for (u = -1; u <= mo->max_offset_x; 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 each state i */
	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]); */
	    }
#ifdef DEBUG
	    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);
	    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));
	    
	    /* this is the difference from the main loop:
	       check whether this state could be an initial state and add the
	       initial probability */
	    if (log_b_i == +1 ) {
	      set_phi(pv, u, v, i, +1);
	    }
	    else {
	      if (max_value == -DBL_MAX)
		set_phi(pv, u, v, i, +1);
	      else
		set_phi(pv, u, v, i, max_value);
	      /* if (mo->s[i].pi != 0 && mo->s[i].offset_x - 1 == u && 
		 mo->s[i].offset_y - 1 == v) { */
	      if (mo->s[i].log_pi != 1 && mo->s[i].offset_x - 1 == u && 
		  mo->s[i].offset_y - 1 == v) {
		set_phi(pv, u, v, i, mo->s[i].log_pi);
#ifdef DEBUG
		printf("Initial log prob state %i at (%i, %i) = %f\n", i, u, v, get_phi(pv, u, v, 0, 0, i));
		printf("Characters emitted X: %i, Y: %i\n", 
		       ghmm_dpseq_get_char(X, mo->s[i].alphabet, u),
		       ghmm_dpseq_get_char(Y, mo->s[i].alphabet, v));
#endif
	      }
	      if (get_phi(pv, u, v, 0, 0, i) != 1)
		set_phi(pv, u, v, i, get_phi(pv, u, v, 0, 0, i) + log_b_i);
	    }
	  }
	  /* if (v == 0) {
	     printf"(%i, %i, %i) preceding %i\n", u, v, i, pv->psi[u][v][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 the old phi values */
      push_back_phi(pv, Y->length);
    } /* End for u in X */
}

/*============================================================================*/
/* call this function when you are at position (x, y) and want to get the
   probability of state 'state' which is offset_x and offset_y away from x,y */
static double get_phi(plocal_store_t * pv, int x, int y, int offset_x, int offset_y, int state){
#ifdef DEBUG
  if (y > pv->len_y || state > pv->mo->N || y < - pv->mo->max_offset_y)
    fprintf(stderr, "get_phi: out of bounds %i %i %i\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, state);
#endif

  if (y - offset_y + pv->mo->max_offset_y >= 0)
    return pv->phi[offset_x][y - offset_y + pv->mo->max_offset_y][state];
  else
    return 1;
}

/*============================================================================*/
/* set the value of this matrix cell */
static void set_phi(plocal_store_t * pv, int x, int y, int state, double prob){
#ifdef DEBUG
  if (y > pv->len_y || state > pv->mo->N || y < - pv->mo->max_offset_y)
    fprintf(stderr, "set_phi: out of bounds %i %i %i %f\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, 
	    state, prob);
#endif
  pv->phi[0][y + pv->mo->max_offset_y][state] = prob;
}

/*============================================================================*/
/* since we only keep the frontier we have to push back when ever a row is
   complete */
static void push_back_phi(plocal_store_t * pv, int length_y){
  int off_x, y, j;
  /* push back the old phi values */
  for (off_x=pv->mo->max_offset_x; off_x>0; off_x--)
    for (y=0; y<length_y + pv->mo->max_offset_y + 1; y++)
      for (j=0; j<pv->mo->N; j++)
	pv->phi[off_x][y][j] = pv->phi[off_x-1][y][j];
}

/*============================================================================*/
static void set_psi(plocal_store_t * pv, int x, int y, int state, int from_state){
  /* shift by max_offsets for negative indices */
#ifdef DEBUG
  if (x > pv->len_x || y > pv->len_y || state > pv->mo->N || 
      x < - pv->mo->max_offset_x || y < - pv->mo->max_offset_y)
    fprintf(stderr, "set_psi: out of bounds %i %i %i %i\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, 
	    state, from_state);
#endif
  pv->psi[x + pv->mo->max_offset_x][y + pv->mo->max_offset_y][state] = from_state;
}

/*============================================================================*/
static int get_psi(plocal_store_t * pv, int x, int y, int state) {
  /* shift by max_offsets for negative indices*/
#ifdef DEBUG
  if (x > pv->len_x || y > pv->len_y || state > pv->mo->N || 
      x < - pv->mo->max_offset_x || y < - pv->mo->max_offset_y)
    fprintf(stderr, "get_psi: out of bounds %i %i %i\n", 
	    x + pv->mo->max_offset_x, y + pv->mo->max_offset_y, 
	    state);
#endif
  return pv->psi[x + pv->mo->max_offset_x][y + pv->mo->max_offset_y][state];
}

/*============================================================================*/
int *ghmm_dpmodel_viterbi_test(ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y,
			  double *log_p, int *path_length) {
  plocal_store_t *pv;
  printf("---- viterbi test -----\n");
  /*ghmm_dpmodel_print(mo);*/
  pv = pviterbi_alloc(mo, X->length, Y->length);
  printf("try free within pviterbi_test\n");
  pviterbi_free(&pv, mo->N, X->length, Y->length, mo->max_offset_x , 
		mo->max_offset_y);
  printf("OK\n");
  return NULL;
}


/*============================================================================*/
int *ghmm_dpmodel_viterbi(ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y, double *log_p,
		     int *path_length) {
  return ghmm_dpmodel_viterbi_variable_tb(mo, X, Y, log_p, path_length, -1);
}

/*============================================================================*/
int *ghmm_dpmodel_viterbi_variable_tb(ghmm_dpmodel *mo, ghmm_dpseq * X, ghmm_dpseq * Y,
				 double *log_p, int *path_length,
				 int start_traceback_with) {
#define CUR_PROC "ghmm_dpmodel_viterbi"
  int u, v, j, i, off_x, off_y, current_state_index;
  double value, max_value, previous_prob;  
  plocal_store_t *pv;
  int *state_seq = NULL;
  int emission;
  double log_b_i, log_in_a_ij;
  double (*log_in_a)(plocal_store_t*, int, int, ghmm_dpseq*, ghmm_dpseq*, int, int);

  /* printf("---- viterbi -----\n"); */
  i_list * state_list;
  state_list = ighmm_list_init_list();
  log_in_a = &sget_log_in_a;
  /* int len_path  = mo->N*len; the length of the path is not known apriori */

/*   if (mo->model_type & kSilentStates &&  */
/*       mo->silent != NULL &&  */
/*       mo->topo_order == NULL) { */
/*     ghmm_dmodel_topo_order( mo );  */
/*   } */

  /* Allocate the matrices log_in_a, log_b,Vektor phi, phi_new, Matrix psi */
  pv = pviterbi_alloc(mo, X->length, Y->length);
  if (!pv)                        { GHMM_LOG_QUEUED(LCONVERTED); goto STOP; }

  /* Precomputing the log(a_ij) and log(bj(ot)) */
  pviterbi_precompute(mo, pv);
  /* Initialize the lookback matrix (for positions [-offsetX,0], [-1, len_y]*/
  init_phi(pv, X, Y);
  
  /* u > max_offset_x , v starts -1 to allow states with offset_x == 0 
     which corresponds to a series of gap states before reading the first 
     character of x at position x=0, y=v */
  /** THIS IS THE MAIN RECURRENCE **/