std_funcs.c

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

  /* find correct alphabet */
  
  while(fgets(ps, 2048, priorfile) != NULL) {
    if(strncmp(ps, "ALPHABET:", 9) == 0) {
      cur_alphabet = atoi(&ps[9]);
      if(cur_alphabet == alphabet) {
	break;
      }
    }
  }
  if(cur_alphabet != alphabet) {
    return 0;
  }
  while(fgets(ps, 2048, priorfile) != NULL) {
    if(strncmp(ps, "START", 5) == 0) {
      break;
    }
  }

  /* put nr of components in struct */
  if(fgets(ps, 2048, priorfile) != NULL) {
  }
  else {
    return -1;
  }
  while(*ps == '#' || *ps == '\n') {
    if(fgets(ps, 2048, priorfile) != NULL) {
    }
    else {
      return -1;
    }
  }
  em_di->nr_components = atoi(&ps[0]);
  em_di->alphabet_size = a_size;
  
  /* allocate memory for arrays and matrix to this prior struct */
  em_di->q_values = malloc_or_die(em_di->nr_components *
				  sizeof(double));
  em_di->alpha_sums = malloc_or_die(em_di->nr_components *
				    sizeof(double));
  em_di->logbeta_values =
    malloc_or_die(em_di->nr_components * sizeof(double));
  em_di->prior_values = malloc_or_die(em_di->nr_components *
				      a_size * sizeof(double));
  
  for(j = 0; j < em_di->nr_components; j++) {
    /* put q-value in array  */
    if(fgets(ps, 2048, priorfile) != NULL) {
    }
    else {
      return -1;
    }
    while(*ps == '#' || *ps == '\n') {
      if(fgets(ps, 2048, priorfile) != NULL) {
      }
      else {
	return -1;
      }
    }
    q_value = atof(&ps[0]);
    *(em_di->q_values + j) = q_value; 
#ifdef DEBUG_PRI
    printf("q_value = %f\n", *(em_di->q_values + j));
#endif
    
    /* put alpha-values of this component in matrix */
    alpha_sum = 0.0;
    k = 0;
    if(fgets(ps, 2048, priorfile) != NULL) {
    }
    else {
      return -1;
    }
    while(*ps == '#' || *ps == '\n') {
      if(fgets(ps, 2048, priorfile) != NULL) {
      }
      else {
	return -1;
      }
    }
    pri = &ps[0];
    for(k = 0; k < a_size; k++) {
      alpha_value = strtod(pri, &pri);
      alpha_sum += alpha_value;
      *((em_di->prior_values) +
	get_mtx_index(j, k, a_size)) = alpha_value;
    }
    
    /* put sum of alphavalues in array */
    *((em_di->alpha_sums) + j) = alpha_sum; 
    
    /* calculate logB(alpha) for this compoment, store in array*/
    logbeta = 0;
    for(k = 0; k < a_size; k++) {
      logbeta += lgamma(*(em_di->prior_values +
			  get_mtx_index(j, k, a_size)));
      
#ifdef DEBUG_PRI
      printf("prior_value = %f\n", *((em_di->prior_values) +
				     get_mtx_index(j, k, a_size)));
      printf("lgamma_value = %f\n", lgamma(*((em_di->prior_values) +
					     get_mtx_index(j, k, a_size))));
#endif
    }
    logbeta = logbeta - lgamma(*(em_di->alpha_sums + j));
    *(em_di->logbeta_values + j) = logbeta;
  }
  
#ifdef DEBUG_PRI
  dump_prior_struct(em_di);
  exit(0);
#endif
  
  return 1;

}




int locked_state(struct hmm_s *hmmp, int v)
{
  if(*(hmmp->locked_vertices + v) == YES) {
    return YES;
  }
  else {
    return NO;
  }
}

int locked_state_multi(struct hmm_multi_s *hmmp, int v)
{
  if(*(hmmp->locked_vertices + v) == YES) {
    return YES;
  }
  else {
    return NO;
  }
}

int get_best_reliability_score(double reliability_score_1, double reliability_score_2, double reliability_score_3)
{
  int max = 2;
  if(reliability_score_1 > reliability_score_2 && reliability_score_1 > reliability_score_3) {
    max = 1;
  }
  else if(reliability_score_3 > reliability_score_2) {
    max = 3;
  }
  return max;
}





void itoa(char* s, int nr) {
  int dec, sign;
  char temp[30];
  
  strcpy(temp, fcvt(nr, 0, &dec, &sign));
  if(sign == POS) {
    strncpy(s, temp, dec);
    s[dec] = '\0';
  }
  else {
    s[0] = '-';
    strncpy(s+1, temp, dec);
    s[dec+1] = '\0';
  }


} 

void ftoa(char* s, double nr, int prec) {
  int dec, sign;
  char temp[30];
  int i, pos;

  strcpy(temp, fcvt(nr, prec, &dec, &sign));
  if(sign == POS) {
    if(dec <= 0) {
      s[0] = '0';
      s[1] = '.';
      pos = 2;
      for(i = 0; i > dec; i--) {
	s[pos] = '0';
	pos++;
      } 
      strncpy(s+pos, temp, prec);
      s[pos+prec] = '\0';
    }
    else {
      strncpy(s, temp, dec);
      s[dec] = '.';
      strncpy(s+dec+1, temp+dec, prec);
      s[dec+1+prec] = '\0';
    }
  }
  else {
     if(dec <= 0) {
      s[0] = '-';
      s[1] = '0';
      s[2] = '.';
      pos = 3;
      for(i = 0; i > dec; i--) {
	s[pos] = '0';
	pos++;
      }
      strncpy(s+pos, temp, prec);
      s[pos+prec] = '\0';
    }
    else {
      s[0] = '-';
      strncpy(s+1, temp, dec);
      s[dec+1] = '.';
      strncpy(s+dec+2, temp+dec, prec);
      s[dec+2+prec] = '\0';
    }
  }
}


void hmm_garbage_collection(FILE *hmmfile, struct hmm_s *hmmp)
{
  int i;

  free(hmmp->transitions);
  free(hmmp->log_transitions);
  free(hmmp->emissions);
  free(hmmp->log_emissions);
  for(i = 0; i < hmmp->nr_m; i++) {
    free((*(hmmp->modules + i))->vertices);
  }
  free(hmmp->silent_vertices);
  free(hmmp->vertex_labels);
  free(hmmp->labels);
  free(hmmp->vertex_trans_prior_scalers);
  free(hmmp->vertex_emiss_prior_scalers);
  free(hmmp->modules);
  free(hmmp->to_trans_array);
  free(hmmp->from_trans_array);
  free(hmmp->to_silent_trans_array);
  free(hmmp->tot_transitions);
  free(hmmp->max_log_transitions);
  free(hmmp->tot_from_trans_array);
  free(hmmp->tot_to_trans_array);
  free(hmmp->distrib_groups);
  free(hmmp->trans_tie_groups);
  for(i = 0; i < hmmp->nr_ed; i++) {
    free(hmmp->emission_dirichlets->q_values);
    free(hmmp->emission_dirichlets->alpha_sums);
    free(hmmp->emission_dirichlets->logbeta_values);
    free(hmmp->emission_dirichlets->prior_values);
  }
  free(hmmp->emission_dirichlets);
  free(hmmp->ed_ps);
  fclose(hmmfile);
}

void hmm_garbage_collection_multi(FILE *hmmfile, struct hmm_multi_s *hmmp)
{
  int i;

  free(hmmp->transitions);
  free(hmmp->log_transitions);
  free(hmmp->emissions);
  free(hmmp->log_emissions);
  if(hmmp->nr_alphabets > 1) {
    free(hmmp->emissions_2);
    free(hmmp->log_emissions_2);
  }
  if(hmmp->nr_alphabets > 2) {
    free(hmmp->emissions_3);
    free(hmmp->log_emissions_3);
  }
  if(hmmp->nr_alphabets > 3) {
    free(hmmp->emissions_4);
    free(hmmp->log_emissions_4);
  }
  for(i = 0; i < hmmp->nr_m; i++) {
    free((*(hmmp->modules + i))->vertices);
  }
  free(hmmp->silent_vertices);
  free(hmmp->vertex_labels);
  free(hmmp->labels);
  free(hmmp->vertex_trans_prior_scalers);
  free(hmmp->vertex_emiss_prior_scalers);
  if(hmmp->nr_alphabets > 1) {
    free(hmmp->vertex_emiss_prior_scalers_2);
  }
  if(hmmp->nr_alphabets > 2) {
    free(hmmp->vertex_emiss_prior_scalers_3);
  }
  if(hmmp->nr_alphabets > 3) {
    free(hmmp->vertex_emiss_prior_scalers_4);
  }
  free(hmmp->modules);
  free(hmmp->to_trans_array);
  free(hmmp->from_trans_array);
  free(hmmp->to_silent_trans_array);
  free(hmmp->tot_transitions);
  free(hmmp->max_log_transitions);
  free(hmmp->tot_from_trans_array);
  free(hmmp->tot_to_trans_array);
  free(hmmp->distrib_groups);
  free(hmmp->trans_tie_groups);
  for(i = 0; i < hmmp->nr_ed; i++) {
    free(hmmp->emission_dirichlets->q_values);
    free(hmmp->emission_dirichlets->alpha_sums);
    free(hmmp->emission_dirichlets->logbeta_values);
    free(hmmp->emission_dirichlets->prior_values);
  }
  free(hmmp->emission_dirichlets);
  free(hmmp->ed_ps);
  if(hmmp->nr_alphabets > 1) {
    for(i = 0; i < hmmp->nr_ed_2; i++) {
      free(hmmp->emission_dirichlets_2->q_values);
      free(hmmp->emission_dirichlets_2->alpha_sums);
      free(hmmp->emission_dirichlets_2->logbeta_values);
      free(hmmp->emission_dirichlets_2->prior_values);
    }
    free(hmmp->emission_dirichlets_2);
    free(hmmp->ed_ps_2);
  }
  if(hmmp->nr_alphabets > 2) {
    for(i = 0; i < hmmp->nr_ed_3; i++) {
      free(hmmp->emission_dirichlets_3->q_values);
      free(hmmp->emission_dirichlets_3->alpha_sums);
      free(hmmp->emission_dirichlets_3->logbeta_values);
      free(hmmp->emission_dirichlets_3->prior_values);
    }
    free(hmmp->emission_dirichlets_3);
    free(hmmp->ed_ps_3);
  }
  if(hmmp->nr_alphabets > 3) {
    for(i = 0; i < hmmp->nr_ed_4; i++) {
      free(hmmp->emission_dirichlets_4->q_values);
      free(hmmp->emission_dirichlets_4->alpha_sums);
      free(hmmp->emission_dirichlets_4->logbeta_values);
      free(hmmp->emission_dirichlets_4->prior_values);
    }
    free(hmmp->emission_dirichlets_4);
    free(hmmp->ed_ps_4);
  }
  if(hmmfile != NULL) {
    fclose(hmmfile);
  }
}

void hmm_garbage_collection_multi_no_dirichlet(FILE *hmmfile, struct hmm_multi_s *hmmp)
{
  int i;

  //printf("transitions\n");
  free(hmmp->transitions);
  free(hmmp->log_transitions);
  free(hmmp->emissions);
  //printf("log_emissions\n");
  free(hmmp->log_emissions);
  if(hmmp->nr_alphabets > 1) {
    free(hmmp->emissions_2);
    free(hmmp->log_emissions_2);
  }
  if(hmmp->nr_alphabets > 2) {
    free(hmmp->emissions_3);
    free(hmmp->log_emissions_3);
  }
  if(hmmp->nr_alphabets > 3) {
    free(hmmp->emissions_4);
    free(hmmp->log_emissions_4);
  }
  for(i = 0; i < hmmp->nr_m; i++) {
    //free((*(hmmp->modules + i))->vertices);
  }
  //printf("silent_vertices\n");
  free(hmmp->silent_vertices);
  //printf("vertex_labels\n");
  free(hmmp->vertex_labels);
  //printf("labels\n");
  free(hmmp->labels);
  //printf("trans_prior_scalers\n");
  free(hmmp->vertex_trans_prior_scalers);
  //printf("emiss_prior_scalers\n");
  free(hmmp->vertex_emiss_prior_scalers);
  if(hmmp->nr_alphabets > 1) {
    free(hmmp->vertex_emiss_prior_scalers_2);
  }
  if(hmmp->nr_alphabets > 2) {
    free(hmmp->vertex_emiss_prior_scalers_3);
  }
  if(hmmp->nr_alphabets > 3) {
    free(hmmp->vertex_emiss_prior_scalers_4);
  }
  //printf("modules\n");
  free(hmmp->modules);
  free(hmmp->to_trans_array);
  free(hmmp->from_trans_array);
  free(hmmp->to_silent_trans_array);
  //printf("tot_transitionss\n");
  free(hmmp->tot_transitions);
  free(hmmp->max_log_transitions);
  free(hmmp->tot_from_trans_array);
  free(hmmp->tot_to_trans_array);
  free(hmmp->distrib_groups);
  //printf("trans_tie_groups\n");
  free(hmmp->trans_tie_groups);
  
  if(hmmfile != NULL) {
    fclose(hmmfile);
  }
}

void msa_seq_garbage_collection_multi(struct msa_sequences_multi_s *msa_seq_info, int nr_alphabets)
{
  
}

void seq_garbage_collection_multi(struct sequences_multi_s *seq_info, int nr_alphabets)
{
  
}

struct letter_s* get_reverse_seq(struct letter_s *seq_const, int seq_len)
{
  struct letter_s *reverse_seq, *seq;
  int i;

  seq = seq_const;
  reverse_seq = (struct letter_s*)(malloc_or_die((seq_len+1)* sizeof(struct letter_s)));
  for(i = seq_len - 1; i >= 0; i--) {
    memcpy(reverse_seq + i, seq, sizeof(struct letter_s));
    seq++;
  }
  
  memcpy(reverse_seq + seq_len, seq, sizeof(struct letter_s));
  

#ifdef DEBUG_REVERSE_SEQ
  dump_seq(reverse_seq);
  dump_seq(seq_const);
  printf("reverse_seq = %x\n", reverse_seq);
#endif
  
  return reverse_seq;
}

void get_reverse_seq_multi(struct sequence_multi_s *seqs, struct letter_s **reverse_seq_1,
			   struct letter_s **reverse_seq_2, struct letter_s **reverse_seq_3,
			   struct letter_s **reverse_seq_4, struct hmm_multi_s *hmmp, int seq_len)
{
  /* letter_s pointers allocated here must be freed by caller */

  struct letter_s *reverse_seq, *seq, *seq_const;
  int i,j;
  *reverse_seq_1 = NULL;
  *reverse_seq_2 = NULL;
  *reverse_seq_3 = NULL;
  *reverse_seq_4 = NULL;

  for(j = 0; j < hmmp->nr_alphabets; j++) {
    reverse_seq = (struct letter_s*)(malloc_or_die((seq_len+1)* sizeof(struct letter_s)));
    if(j == 0) {
      seq = seqs->seq_1;
      *reverse_seq_1 = reverse_seq;
    }
    if(j == 1) {
      seq = seqs->seq_2;
      *reverse_seq_2 = reverse_seq;
    }
    if(j == 1) {
      seq = seqs->seq_3;
      *reverse_seq_3 = reverse_seq;
    }
    if(j == 1) {
      seq = seqs->seq_4;
      *reverse_seq_4 = reverse_seq;
    }
    seq_const = seq;
    for(i = seq_len - 1; i >= 0; i--) {
      memcpy(reverse_seq + i, seq, sizeof(struct letter_s));
      seq++;
    }
    
    memcpy(reverse_seq + seq_len, seq, sizeof(struct letter_s));
    
#ifdef DEBUG_REVERSE_SEQ
    printf("seq_len = %d\n", seq_len);
    dump_seq(reverse_seq);
    dump_seq(*reverse_seq_1);
    dump_seq(seq_const);
    printf("reverse_seq = %x\n", reverse_seq);
#endif
  }
}


void get_reverse_msa_seq(struct msa_sequences_s *msa_seq_infop, struct msa_sequences_s *reverse_msa_seq_infop, int a_size)
{

  /* note that this function does not implement the gaps-pointing, everything points to END */
  
  int i,j;

  reverse_msa_seq_infop->nr_seqs = msa_seq_infop->nr_seqs;
  reverse_msa_seq_infop->msa_seq_length = msa_seq_infop->msa_seq_length;
  reverse_msa_seq_infop->nr_lead_columns = msa_seq_infop->nr_lead_columns;
  reverse_msa_seq_infop->msa_seq = (struct msa_letter_s*)(malloc_or_die(msa_seq_infop->msa_seq_length * (a_size + 1) *
									sizeof(struct msa_letter_s)));
  reverse_msa_seq_infop->gaps = (int**)(malloc_or_die(msa_seq_infop->msa_seq_length * sizeof(int*) +
						      msa_seq_infop->msa_seq_length * sizeof(int)));
  reverse_msa_seq_infop->lead_columns_start = (int*)(malloc_or_die((msa_seq_infop->nr_lead_columns +1)* sizeof(int)));
  reverse_msa_seq_infop->lead_columns_end = reverse_msa_seq_infop->lead_columns_start + (msa_seq_infop->nr_lead_columns);
  reverse_msa_seq_infop->gap_shares = (double*)(malloc_or_die(msa_seq_infop->msa_seq_length * sizeof(double)));