std_funcs.c

马尔科夫模型的java版本实现

	  replacement_lettersp->letters_1 = (struct letter_s*)(malloc_or_die(nr_repl_letters * sizeof(struct letter_s)));
	  replacement_lettersp->probs_1 = (double*)(malloc_or_die(nr_repl_letters * a_size * sizeof(double)));
	  break;
	}
	if(l == 1) {
	  nr_repl_letters = atoi(s);
	  replacement_lettersp->nr_rl_2 = nr_repl_letters;
	  replacement_lettersp->letters_2 = (struct letter_s*)(malloc_or_die(nr_repl_letters * sizeof(struct letter_s)));
	  replacement_lettersp->probs_2 = (double*)(malloc_or_die(nr_repl_letters * a_size * sizeof(double)));
	  break;
	}
	if(l == 2) {
	  nr_repl_letters = atoi(s);
	  replacement_lettersp->nr_rl_3 = nr_repl_letters;
	  replacement_lettersp->letters_3 = (struct letter_s*)(malloc_or_die(nr_repl_letters * sizeof(struct letter_s)));
	  replacement_lettersp->probs_3 = (double*)(malloc_or_die(nr_repl_letters * a_size * sizeof(double)));
	  break;
	}
	if(l == 3) {
	  nr_repl_letters = atoi(s);
	  replacement_lettersp->nr_rl_4 = nr_repl_letters;
	  replacement_lettersp->letters_4 = (struct letter_s*)(malloc_or_die(nr_repl_letters * sizeof(struct letter_s)));
	  replacement_lettersp->probs_4 = (double*)(malloc_or_die(nr_repl_letters * a_size * sizeof(double)));
	  break;
	}
      }
    }

    cur_letter = 0;
    letterindex = 0;
    while (fgets(s, MAX_LINE, replfile) != NULL) {
      if(letterindex >= nr_repl_letters) {
	break;
      }
      if(s[0] == '#' || s[0] == '\n') {
	
      }
      else {
	i = 0;
	while(s[i] != ' ') {
	  if(l == 0) {
	    *((replacement_lettersp->letters_1 + cur_letter)->letter + i) = s[i]; 
	  }
	  if(l == 1) {
	    *((replacement_lettersp->letters_2 + cur_letter)->letter + i) = s[i];
	  }
	  if(l == 2) {
	    *((replacement_lettersp->letters_3 + cur_letter)->letter + i) = s[i];
	    }
	  if(l == 3) {
	    *((replacement_lettersp->letters_4 + cur_letter)->letter + i) = s[i];
	  }
	  i++;
	}
	
	if(l == 0) {
	  *((replacement_lettersp->letters_1 + cur_letter)->letter + i) = '\0';
	}
	if(l == 1) {
	  *((replacement_lettersp->letters_2 + cur_letter)->letter + i) = '\0';
	}
	if(l == 2) {
	  *((replacement_lettersp->letters_3 + cur_letter)->letter + i) = '\0';
	}
	if(l == 3) {
	  *((replacement_lettersp->letters_4 + cur_letter)->letter + i) = '\0';
	}
	while(s[i] == ' ' || s[i] == '=') {
	  i++;
	}
	done = NO;
	while(done == NO) {
	  j = 0;
	  while(s[i] != ':') {
	    *(le.letter + j) = s[i];
	    i++;
	    j++;
	  }
	  *(le.letter + j) = '\0';
	  i++;
	  
	  prob = atof(&s[i]);
	  /* get alphabet index for correct alphabet (must hardcode this to the same order as in the hmm */
	  a_index = get_alphabet_index(&le, alphabet, a_size);

	  if(l == 0) {
	    *(replacement_lettersp->probs_1 + get_mtx_index(cur_letter, a_index, a_size)) = prob;
	  }
	  if(l == 1) {
	    *(replacement_lettersp->probs_2 + get_mtx_index(cur_letter, a_index, a_size)) = prob;
	  }
	  if(l == 2) {
	    *(replacement_lettersp->probs_3 + get_mtx_index(cur_letter, a_index, a_size)) = prob;
	  }
	  if(l == 3) {
	    *(replacement_lettersp->probs_4 + get_mtx_index(cur_letter, a_index, a_size)) = prob;
	  }
	  while(s[i] != ' ' && s[i] != '\n') {
	    i++;
	  }
	  if(s[i] == '\n') {
	    done = YES;
	  }
	  i++;
	}
	cur_letter++;
	letterindex++;
      }
      
    }
    free(alphabet);
    //dump_replacement_letters_multi(replacement_lettersp, l+1, a_size);
  }
}

int get_replacement_letter_index(struct letter_s *c, struct replacement_letter_s *replacement_letters)
{
  int a_index;
  for(a_index = 0; a_index < replacement_letters->nr_rl; a_index++) {
    if(strcmp((replacement_letters->letters + a_index)->letter, c->letter) == 0) {
      return a_index;
    }  
  }
  
  return -1;
}

int get_replacement_letter_index_multi(struct letter_s *c, struct replacement_letter_multi_s *replacement_letters, int alphabet)
{
  int a_index;
  if(alphabet == 1) {
    for(a_index = 0; a_index < replacement_letters->nr_rl_1; a_index++) {
      if(strcmp((replacement_letters->letters_1 + a_index)->letter, c->letter) == 0) {
	return a_index;
      }  
    }
  }
  if(alphabet == 2) {
    for(a_index = 0; a_index < replacement_letters->nr_rl_2; a_index++) {
      if(strcmp((replacement_letters->letters_2 + a_index)->letter, c->letter) == 0) {
	return a_index;
      }  
    }
  }
  if(alphabet == 3) {
    for(a_index = 0; a_index < replacement_letters->nr_rl_3; a_index++) {
      if(strcmp((replacement_letters->letters_3 + a_index)->letter, c->letter) == 0) {
	return a_index;
      }  
    }
  }
  if(alphabet == 4) {
    for(a_index = 0; a_index < replacement_letters->nr_rl_4; a_index++) {
      if(strcmp((replacement_letters->letters_4 + a_index)->letter, c->letter) == 0) {
	return a_index;
      }  
    }
  }
  
  return -1;
}

int get_replacement_letter_index_single(char *c, struct replacement_letter_s *replacement_letters)
{
  int a_index;
  for(a_index = 0; a_index < replacement_letters->nr_rl; a_index++) {
    if(((replacement_letters->letters + a_index)->letter)[0] == c) {
      return a_index;
    }  
  }
  
  return -1;
}

char* sequence_as_string(struct letter_s *sequence)
{
  struct letter_s *c;
  char *s;
  int s_length;

  s_length = 0;
  c = sequence;
  while(c->letter[0] != '\0') {
    s_length++;
    c++;
  }
  
  s = (char*)(malloc_or_die(s_length * 5 * sizeof(char)));

  c = sequence;
  while(c->letter[0] != '\0') {
    strcat(s,c->letter);
    c++;
  }
  
  printf("%s\n", s);
  free(s);
  return NULL;
}

void get_viterbi_label_path(struct viterbi_s *cur, struct hmm_s *hmmp,
			    struct viterbi_s *viterbi_mtxp, int row, int row_size, char *labels, int *ip)
{
  struct path_element *p_el;
  
  if(cur->prev == 0) {
    p_el = cur->prevp;
    while(p_el->next != NULL) {
      labels[*ip] = *(hmmp->vertex_labels + p_el->vertex);
      *ip = (*ip) + 1;
      p_el++;
    }
  }
  else {
    get_viterbi_label_path(viterbi_mtxp + get_mtx_index(row-1, cur->prev, row_size), hmmp,
		      viterbi_mtxp, row-1, row_size, labels, ip);
    p_el = cur->prevp;
    labels[*ip] = *(hmmp->vertex_labels + (int)(cur->prev));
    *ip = (*ip) + 1;
    while(p_el->next != NULL) {
      p_el++;
      labels[*ip] = *(hmmp->vertex_labels + p_el->vertex);
      *ip = (*ip) + 1;
    }
  }
}

void get_viterbi_label_path_multi(struct viterbi_s *cur, struct hmm_multi_s *hmmp,
			    struct viterbi_s *viterbi_mtxp, int row, int row_size, char *labels, int *ip)
{
  struct path_element *p_el;
  
  if(cur->prev == 0) {
    p_el = cur->prevp;
    while(p_el->next != NULL) {
      labels[*ip] = *(hmmp->vertex_labels + p_el->vertex);
      *ip = (*ip) + 1;
      p_el++;
    }
  }
  else {
    get_viterbi_label_path_multi(viterbi_mtxp + get_mtx_index(row-1, cur->prev, row_size), hmmp,
				 viterbi_mtxp, row-1, row_size, labels, ip);
    p_el = cur->prevp;
    labels[*ip] = *(hmmp->vertex_labels + (int)(cur->prev));
    *ip = (*ip) + 1;
    while(p_el->next != NULL) {
      p_el++;
      labels[*ip] = *(hmmp->vertex_labels + p_el->vertex);
      *ip = (*ip) + 1;
    }
  }
}


void get_viterbi_path(struct viterbi_s *cur, struct hmm_s *hmmp,
		      struct viterbi_s *viterbi_mtxp, int row, int row_size, int *path, int *ip)
{
  struct path_element *p_el;
  
  //printf("cur->prev = %d\n", cur->prev);
  if(cur->prev == 0) {
    p_el = cur->prevp;
    while(p_el->next != NULL) {
      path[*ip] = p_el->vertex;
      *ip = (*ip) + 1;
      p_el++;
    }
  }
  else {
    get_viterbi_path(viterbi_mtxp + get_mtx_index(row-1, cur->prev, row_size), hmmp,
			   viterbi_mtxp, row-1, row_size, path, ip);
    p_el = cur->prevp;
    path[*ip] = p_el->vertex;
    //printf("%d ", path[*ip]);
    *ip = (*ip) + 1;
    while(p_el->next != NULL) {
      p_el++;
      path[*ip] = p_el->vertex;
      *ip = (*ip) + 1;
    }
  }
}

void get_viterbi_path_multi(struct viterbi_s *cur, struct hmm_multi_s *hmmp,
			    struct viterbi_s *viterbi_mtxp, int row, int row_size, int *path, int *ip)
{
  struct path_element *p_el;
  
  //printf("cur->prev = %d\n", cur->prev);
  if(cur->prev == 0) {
    p_el = cur->prevp;
    while(p_el->next != NULL) {
      path[*ip] = p_el->vertex;
      *ip = (*ip) + 1;
      p_el++;
    }
  }
  else {
    get_viterbi_path_multi(viterbi_mtxp + get_mtx_index(row-1, cur->prev, row_size), hmmp,
			   viterbi_mtxp, row-1, row_size, path, ip);
    p_el = cur->prevp;
    path[*ip] = p_el->vertex;
    //printf("%d ", path[*ip]);
    *ip = (*ip) + 1;
    while(p_el->next != NULL) {
      p_el++;
      path[*ip] = p_el->vertex;
      *ip = (*ip) + 1;
    }
  }
}


void loosen_labels(char *labels, char *loose_labels, int label_looseness, int seq_len)
{
  int *locked_labels;
  int reg_len;
  int i,j;
  char cur;
  
  /* initial memory and copy stuff */
  locked_labels = (int*)(malloc_or_die(seq_len * sizeof(int)));
  memcpy(loose_labels, labels, seq_len * sizeof(char));
  

  /* lock middle labels */
  reg_len = 1;
  cur = labels[0];
  locked_labels[0] = 1;
  for(i = 1; i < seq_len; i++) {
    if(labels[i] != cur) {
      cur = labels[i];
      if(reg_len == i) {
	/* first reg shift, do nothing */
      }
      else {
	if(reg_len % 2 == 1) {
	  locked_labels[i - (reg_len/2 + 1)] = 1;
	}
	else {
	  locked_labels[i - reg_len/2] = 1;
	  locked_labels[i - (reg_len/2 + 1)] = 1;
	}
      }
      reg_len = 1;
    }
    else {
      reg_len++;
    }
  }

  /* loosen labels */
  reg_len = 1;
  cur = labels[0];
  for(i = 1; i < seq_len; i++) {
    if(labels[i] != cur) {
      for(j = 0; j < label_looseness; j++) {
	if(locked_labels[i + j] == 0 && locked_labels[i-1-j] == 0) {
	  //printf("i = %d, j = %d\n", i , j);
	  loose_labels[i + j] = '.';
	  loose_labels[i-1-j] = '.';
	}
	else {
	  break;
	}
      }

      cur = labels[i];
      reg_len = 1;
    }
    else {
      reg_len++;
    }

  }
  free(locked_labels);

  //dump_labeling(labels, seq_len);
  //dump_labeling(loose_labels, seq_len);
}

int read_subst_matrix(double **mtxpp, FILE *substmtxfile)
{
  int MAX_LINE = 1000;
  char s[1000];
  int i,j,k;
  int nr_rows;
  int row_le_index, col_le_index;
  struct letter_s row_le, col_le;
  double prob;
  int row_a_index, col_a_index;
  int done;
  int a_size;
  char *alphabet;
  double *mtxp;

  if(substmtxfile == NULL) {
    return NO;
  }
  
  while(fgets(s, MAX_LINE, substmtxfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
    }
    else {
      a_size = atoi(s);
      alphabet = (char*)malloc_or_die(a_size * sizeof(char) * 5);
      *mtxpp = (double*)(malloc_or_die(a_size * (a_size + 1) * sizeof(double)));
      mtxp = *mtxpp;
      break;
    }
  }
  i = 0;
  while(fgets(s, MAX_LINE, substmtxfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
      
    }
    else {
      strcpy(alphabet, s);
      break;
    }
  }

  while (fgets(s, MAX_LINE, substmtxfile) != NULL) {
    if(s[0] == '#' || s[0] == '\n') {
      
    }
    else {
      i = 0;
      while(s[i] != ' ') {
	*(row_le.letter + i) = s[i];
	i++;
      }
      *(row_le.letter + i) = '\0';
      while(s[i] == ' ' || s[i] == '=') {
	i++;
      }
      done = NO;
      while(done == NO) {
	j = 0;
	while(s[i] != ':') {
	  *(col_le.letter + j) = s[i];
	  i++;
	  j++;
	}
	*(col_le.letter + j) = '\0';
	i++;
	
	prob = atof(&s[i]);
	row_a_index = get_alphabet_index(&row_le, alphabet, a_size);
	if(row_a_index < 0) {
	  row_a_index = a_size;
	}
	col_a_index = get_alphabet_index(&col_le, alphabet, a_size);
	*(mtxp + get_mtx_index(row_a_index, col_a_index, a_size)) = prob;
	while(s[i] != ' ' && s[i] != '\n') {
	  i++;
	}
	if(s[i] == '\n') {
	  done = YES;
	}
	i++;
      }
    }
  }
  free(alphabet);
  //dump_subst_mtx(mtxp, a_size);
  //exit(0);  
  return YES;
}


int read_subst_matrix_multi(double **mtxpp, double **mtxpp_2, double **mtxpp_3, double **mtxpp_4, FILE *substmtxfile)
{
  int MAX_LINE = 1000;
  char s[1000];
  int i,j,k, l;
  int nr_rows;
  int row_le_index, col_le_index;
  struct letter_s row_le, col_le;
  double prob;
  int row_a_index, col_a_index;
  int done;
  int a_size;
  char *alphabet;
  double *mtxp;
  int nr_alphabets;
  
  if(substmtxfile == NULL) {
    return NO;
  }
  
  while(fgets(s, MAX_LINE, substmtxfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
    }
    else {
      nr_alphabets = atoi(s);
      break;
    }
  }
  
  for(l = 0; l < nr_alphabets; l++) {
    while(fgets(s, MAX_LINE, substmtxfile) != NULL){
      if(s[0] == '\n' || s[0] == '#') {
      }
      else {
	a_size = atoi(s);
	alphabet = (char*)malloc_or_die(a_size * sizeof(char) * 5);
	if(l == 0) {
	  *mtxpp = (double*)(malloc_or_die(a_size * (a_size + 1) * sizeof(double)));
	  mtxp = *mtxpp;
	}
	if(l == 1) {
	  *mtxpp_2 = (double*)(malloc_or_die(a_size * (a_size + 1) * sizeof(double)));
	  mtxp = *mtxpp_2;
	}