std_funcs.c

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

	if(l == 2) {
	  *mtxpp_3 = (double*)(malloc_or_die(a_size * (a_size + 1) * sizeof(double)));
	  mtxp = *mtxpp_3;
	}
	if(l == 3) {
	  *mtxpp_4 = (double*)(malloc_or_die(a_size * (a_size + 1) * sizeof(double)));
	  mtxp = *mtxpp_4;
	}
	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 if(strncmp(s, "END", 3) == 0) {
	break;
      }
      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_prior_file(struct emission_dirichlet_s *em_di, struct hmm_s *hmmp, FILE *priorfile)
{
  int j,k;
  double q_value, alpha_value, alpha_sum, logbeta;
  char s[2048];
  char ps[2048];
  char *file_name;
  char *pri;
  
  rewind(priorfile);

  /* put default name */
  strcpy(em_di->name, "default");

  /* 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]);
  
  /* 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 *
				     hmmp->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 < hmmp->a_size; k++) {
      alpha_value = strtod(pri, &pri);
      alpha_sum += alpha_value;
      *((em_di->prior_values) +
	get_mtx_index(j, k, hmmp->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 < hmmp->a_size; k++) {
      logbeta += lgamma(*(em_di->prior_values +
			  get_mtx_index(j, k, hmmp->a_size)));
      
#ifdef DEBUG_PRI
      printf("prior_value = %f\n", *((em_di->prior_values) +
				     get_mtx_index(j, k, hmmp->a_size)));
      printf("lgamma_value = %f\n", lgamma(*((em_di->prior_values) +
					     get_mtx_index(j, k, hmmp->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
}

int read_frequencies(FILE *freqfile, double **aa_freqsp)
{
  char ps[2048];
  int cur;
  int a_size;
  double *aa_freqs;
  
  /* read frequencies */
  if(fgets(ps, 2048, freqfile) != NULL) {
  }
  else {
    return -1;
  }
  while(*ps == '#' || *ps == '\n') {
    if(fgets(ps, 2048, freqfile) != NULL) {
    }
    else {
      return -1;
    }
  }
  
  

  a_size = atoi(ps);

  aa_freqs = (double*)malloc_or_die(a_size * sizeof(double));
  *aa_freqsp = aa_freqs;
  if(freqfile == NULL) {
    printf("Could not read prior frequencies\n");
    exit(0);
  }

  /* read frequencies */
  for(cur = 0; cur < a_size; cur++) {
    if(fgets(ps, 2048, freqfile) != NULL) {
    }
    else {
      return -1;
    }
    while(*ps == '#' || *ps == '\n') {
      if(fgets(ps, 2048, freqfile) != NULL) {
      }
      else {
	return -1;
      }
    }
    *(aa_freqs + cur) = atof(ps);
    //printf("aa_freqs[%d] = %f\n", cur, *(aa_freqs + cur));
  }
  
  return 1;
  
}


int read_frequencies_multi(FILE *freqfile, double **aa_freqsp, double **aa_freqsp_2, double **aa_freqsp_3, double **aa_freqsp_4)
{
  char ps[2048];
  int cur;
  int a_size;
  double *aa_freqs, *aa_freqs_2, *aa_freqs_3, *aa_freqs_4, *aa_freqs_temp;
  int nr_alphabets;
  int i;

  /* read frequencies */
  if(fgets(ps, 2048, freqfile) != NULL) {
  }
  else {
    return -1;
  }
  while(*ps == '#' || *ps == '\n') {
    if(fgets(ps, 2048, freqfile) != NULL) {
    }
    else {
      return -1;
    }
  }
  
  nr_alphabets = atoi(ps);

  for(i = 0; i <= nr_alphabets; i++) {
    if(fgets(ps, 2048, freqfile) != NULL) {
    }
    else {
      return -1;
    }
    while(*ps == '#' || *ps == '\n') {
      if(fgets(ps, 2048, freqfile) != NULL) {
      }
      else {
	return -1;
      }
    }
    a_size = atoi(ps);
    if(i == 0) {
      aa_freqs = (double*)malloc_or_die(a_size * sizeof(double));
      *aa_freqsp = aa_freqs;
      aa_freqs_temp = aa_freqs;
    }
    if(i == 1) {
      aa_freqs_2 = (double*)malloc_or_die(a_size * sizeof(double));
      *aa_freqsp_2 = aa_freqs_2;
      aa_freqs_temp == aa_freqs_2;
    }
    if(i == 2) {
      aa_freqs_3 = (double*)malloc_or_die(a_size * sizeof(double));
      *aa_freqsp_3 = aa_freqs_3;
      aa_freqs_temp == aa_freqs_3;
    }
    if(i == 3) {
      aa_freqs_4 = (double*)malloc_or_die(a_size * sizeof(double));
      *aa_freqsp_4 = aa_freqs_4;
      aa_freqs_temp == aa_freqs_4;
    }
    
    
    /* read frequencies */
    for(cur = 0; cur < a_size; cur++) {
      if(fgets(ps, 2048, freqfile) != NULL) {
      }
      else {
	return -1;
      }
      while(*ps == '#' || *ps == '\n') {
	if(fgets(ps, 2048, freqfile) != NULL) {
	}
	else {
	  return -1;
	}
      }
      *(aa_freqs_temp + cur) = atof(ps);
      //printf("aa_freqs_temp[%d] = %f\n", cur, *(aa_freqs + cur));
    }
  }
  
  return 1;
  
}




int read_prior_file_multi(struct emission_dirichlet_s *em_di, struct hmm_multi_s *hmmp, FILE *priorfile, int alphabet)
{
  int j,k;
  double q_value, alpha_value, alpha_sum, logbeta;
  char s[2048];
  char ps[2048];
  char *file_name;
  char *pri;
  int a_size;
  
  if(alphabet == 1) {
    a_size = hmmp->a_size;
  }
  if(alphabet == 2) {
    a_size = hmmp->a_size_2;
    if(hmmp->nr_alphabets < 2) {
      printf("Trying to read priorfile for alphabet 2, but hmm only has one alphabet\n");
      exit(0);
    }
  }
  if(alphabet == 3) {
    a_size = hmmp->a_size_3;
    if(hmmp->nr_alphabets < 3) {
      printf("Trying to read priorfile for alphabet 3, but hmm only has two alphabets\n");
      exit(0);
    }
  }
  if(alphabet == 4) {
    a_size = hmmp->a_size_4;
    if(hmmp->nr_alphabets < 4) {
      printf("Trying to read priorfile for alphabet 4, but hmm only has three alphabets\n");
      exit(0);
    }
  }

  if(priorfile == NULL) {
    em_di->nr_components = 0;
    return 0;
  }
  rewind(priorfile);
  
  /* put default name */
  strcpy(em_di->name, "default");

  /* put nr of components in struct */
  if(fgets(ps, 2048, priorfile) != NULL) {
  }
  else {
    return -1;
  }
  while(*ps == '#') {
    if(fgets(ps, 2048, priorfile) != NULL) {
    }
    else {
      return -1;
    }
  }
  em_di->nr_components = atoi(&ps[0]);
  
  /* 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
 

}

int read_multi_prior_file_multi(struct emission_dirichlet_s *em_di, struct hmm_multi_s *hmmp, FILE *priorfile, int alphabet)
{

  /* returns negative value if an error in the priorfile is detected, 0 if there is no prior information for the alphabet
   * and a postive value if prior components were read successfully */

  int j,k;
  double q_value, alpha_value, alpha_sum, logbeta;
  char s[2048];
  char ps[2048];
  char *file_name;
  char *pri;
  int a_size;
  int cur_alphabet;

  cur_alphabet = -1;

  if(alphabet == 1) {
    a_size = hmmp->a_size;
  }
  if(alphabet == 2) {
    a_size = hmmp->a_size_2;
    if(hmmp->nr_alphabets < 2) {
      printf("Trying to read priorfile for alphabet 2, but hmm only has one alphabet\n");
      exit(0);
    }
  }
  if(alphabet == 3) {
    a_size = hmmp->a_size_3;
    if(hmmp->nr_alphabets < 3) {
      printf("Trying to read priorfile for alphabet 3, but hmm only has two alphabets\n");
      exit(0);
    }
  }
  if(alphabet == 4) {
    a_size = hmmp->a_size_4;
    if(hmmp->nr_alphabets < 4) {
      printf("Trying to read priorfile for alphabet 4, but hmm only has three alphabets\n");
      exit(0);
    }
  }
  
  if(priorfile == NULL) {
    em_di->nr_components = 0;
    return 0;
  }

  rewind(priorfile);

  /* put default name */
  strcpy(em_di->name, "default");