std_funcs.c

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

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
//#include <double.h>
#include <limits.h>

#include "structs.h"
#include "funcs.h"

//#define DEBUG_ALPHA
//#define DEBUG_REVERSE_SEQ
//#define DEBUG_PRI

#define POS 0


int verbose = NO;



void *malloc_or_die(int mem_size)
{
  void *mem;
  int i;
  
  if(mem_size > 6000 * 6000 * 8) {
    printf("Trying to allocate to much memory: %d bytes\n", mem_size);
    exit(0);
  }

  if((mem = malloc(mem_size)) != NULL) {
    memset(mem, 0, mem_size);
    return mem;
  }
  else {
    perror("Memory trouble");
    exit(-1);
  }
}

void init_float_mtx(double *mtx, double init_nr, int mtx_size)
{
  int i;
  for(i = 0; i < mtx_size; i++) {
    *mtx = init_nr;
    mtx++;
  }
}

void init_viterbi_s_mtx(struct viterbi_s *mtx, double init_nr, int mtx_size)
{
  int i;
  for(i = 0; i < mtx_size; i++) {
    mtx->prob = init_nr;
    mtx++;
  }
}

int get_mtx_index(int row, int col, int row_size)
{
  return row * row_size + col;
}

int get_alphabet_index(struct letter_s *c, char *alphabet, int a_size)
{
  int a_index;
  int found_index;
  int i;

  for(a_index = 0; a_index < a_size; a_index++) {
    
    found_index = NO;

#ifdef DEBUG_ALPHA
    printf("a_index = %d\n", a_index);
    printf("c = %s\n", c->letter);
    printf("Alphabet = %s\n", alphabet);   
#endif

    for(i = 0; c->letter[i] != '\0'; i++) {
      if(*alphabet == ';') {
	found_index = NO;
	break;
      }
      else if(c->letter[i] == *alphabet) {
	found_index = YES;
      }
      else {
	found_index = NO;
	break;
      }
      alphabet++;
    }
    if(found_index == YES) {
      break;
    }
    while(*alphabet != ';') {
      alphabet++;
    }
    alphabet++;
  }
  
  if(a_index >= a_size) {
    return -1;
  }
  return a_index;
}

int get_alphabet_index_msa_query(char *c, char *alphabet, int a_size)
{
  int a_index;
  int found_index;
  int i;

  for(a_index = 0; a_index < a_size; a_index++) {
    
    found_index = NO;

#ifdef DEBUG_ALPHA
    printf("a_index = %d\n", a_index);
    printf("c = %s\n", c);
    printf("Alphabet = %s\n", alphabet);   
#endif

    for(i = 0; c[i] != '\0'; i++) {
      if(*alphabet == ';') {
	found_index = NO;
	break;
      }
      else if(c[i] == *alphabet) {
	found_index = YES;
      }
      else {
	found_index = NO;
	break;
      }
      alphabet++;
    }
    if(found_index == YES) {
      break;
    }
    while(*alphabet != ';') {
      alphabet++;
    }
    alphabet++;
  }
  
  if(a_index >= a_size) {
    return -1;
  }
  return a_index;
}

/* if letter is known to be only one character */
int get_alphabet_index_single(char *alphabet, char letter, int a_size)
{
  int a_index;
  //printf("alphabet = %s\n", alphabet);
  //printf("letter = %c\n", letter);
  for(a_index = 0; a_index < a_size; a_index++) {
    if(letter == *alphabet) {
      return a_index;
    }
    alphabet++;
    alphabet++;
  }
  
  return -1;

}

int get_seq_length(struct letter_s *s)
{
  int i, len;
  
  len = 0;
  while((s->letter)[0] != '\0') {
    len++;
    s++;
  }
  return len;
}

/* checks if there is a path in the hmm from vertex 'v' to vertex 'w'
 * either directly or via silent states */
int path_length(int v, int w, struct hmm_s *hmmp, int length)
{
  int *xp;
  int tot_length;
  int temp_length;
  
  if(length > MAX_GAP_SIZE) {
    return 0;
  }
  
  tot_length = 0;
  if(*(hmmp->transitions + get_mtx_index(v, w, hmmp->nr_v)) != 0.0) {
    tot_length = length + 1; /* direct path */
  }
  for(xp = hmmp->silent_vertices; *xp != END; xp++) {
    if(*(hmmp->transitions + get_mtx_index(v, *xp, hmmp->nr_v)) != 0.0) {
      if((temp_length = path_length(*xp, w, hmmp, length + 1)) > 0) {
	tot_length += length + temp_length; /* path via a silent state */
      }
    }
  }
  return tot_length; /* if 0, then there is no path */
}

/* checks if there is a path in the hmm from vertex 'v' to vertex 'w'
 * either directly or via silent states */
int path_length_multi(int v, int w, struct hmm_multi_s *hmmp, int length)
{
  int *xp;
  int tot_length;
  int temp_length;
  

  if(length > MAX_GAP_SIZE) {
    return 0;
  }
  
  tot_length = 0;
  if(*(hmmp->transitions + get_mtx_index(v, w, hmmp->nr_v)) != 0.0) {
    tot_length = length + 1; /* direct path */
  }
  for(xp = hmmp->silent_vertices; *xp != END; xp++) {
    if(*(hmmp->transitions + get_mtx_index(v, *xp, hmmp->nr_v)) != 0.0) {
      if((temp_length = path_length_multi(*xp, w, hmmp, length + 1)) > 0) {
	tot_length += length + temp_length; /* path via a silent state */
      }
    }
  }
  return tot_length; /* if 0, then there is no path */
}


struct path_element* get_end_path_start(int l, struct hmm_s *hmmp)
{
  struct path_element *pep, *pep_const;
  
  pep_const = *(hmmp->to_trans_array + l);
  pep = pep_const;
  while(pep->vertex != END) {
    while(pep->next != NULL) {
      pep++;
    }
    if(pep->vertex == hmmp->nr_v-1) {
      return pep_const;
    }
    else {
      pep++;
      pep_const = pep;
    }
  }
  return NULL;
}

struct path_element* get_end_path_start_multi(int l, struct hmm_multi_s *hmmp)
{
  struct path_element *pep, *pep_const;
  
  pep_const = *(hmmp->to_trans_array + l);
  pep = pep_const;
  while(pep->vertex != END) {
    while(pep->next != NULL) {
      pep++;
    }
    if(pep->vertex == hmmp->nr_v-1) {
      return pep_const;
    }
    else {
      pep++;
      pep_const = pep;
    }
  }
  return NULL;
}

void print_seq(struct letter_s *seq, FILE *outfile, int seq_nr, char *name, int seq_length)
{
  int i,j;

  fprintf(outfile, "NR %d: %s\n", seq_nr, name);
  fprintf(outfile, ">");
  for(i = 0; i < seq_length; i++) {
    j = 0;
    while(*((seq+i)->letter + j) != '\0') {
      fputc((int)*((seq+i)->letter + j), outfile);
      fputc((int)';', outfile);
      j++;
    }
  }
  fputc('\n', outfile);
}


char* get_profile_vertex_type(int v, int *silent_vertices)
{
  while(*silent_vertices != END) {
    if(v == *silent_vertices) {
      return "silent\n";
    }
    silent_vertices++;
  }

  return "standard\n";
}

void get_aa_distrib_mtx(FILE *distribmtxfile, struct aa_distrib_mtx_s *aa_distrib_mtxp)
{
  int MAX_LINE = 500;
  char s[500];
  char *temp;
  int i;

  if(distribmtxfile == NULL) {
    aa_distrib_mtxp->a_size = -1;
    return;
  }

  i = 0;
  while(fgets(s, MAX_LINE, distribmtxfile) != NULL) {
    if(s[0] == '\n' || s[0] == '#') {
    }
    else {
      aa_distrib_mtxp->a_size = atoi(s);
      aa_distrib_mtxp->inside_values = (double*)malloc_or_die(aa_distrib_mtxp->a_size * sizeof(double));
      aa_distrib_mtxp->outside_values = (double*)malloc_or_die(aa_distrib_mtxp->a_size * sizeof(double));
      aa_distrib_mtxp->membrane_values = (double*)malloc_or_die(aa_distrib_mtxp->a_size * sizeof(double));
      break;
    }
  }
  while(fgets(s, MAX_LINE, distribmtxfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
      
    }
    else {
      temp = s+2;
      *(aa_distrib_mtxp->inside_values + i) = strtod(temp, &temp);
      *(aa_distrib_mtxp->outside_values + i) = strtod(temp, &temp);
      *(aa_distrib_mtxp->membrane_values + i) = strtod(temp, &temp);
      *(aa_distrib_mtxp->membrane_values + i) += strtod(temp, &temp);
      i++;
    }
  }
  //dump_aa_distrib_mtx(aa_distrib_mtxp);
}


void get_replacement_letters(FILE *replfile, struct replacement_letter_s *replacement_lettersp)
{
  int MAX_LINE = 1000;
  char s[1000];
  int i,j,k;
  int nr_repl_letters;
  int cur_letter;
  struct letter_s le;
  double prob;
  int a_index;
  int done;
  int a_size;
  char *alphabet;
 

  if(replfile == NULL) {
    replacement_lettersp->nr_rl = 0;
    return;
  }
  
  while(fgets(s, MAX_LINE, replfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
    }
    else {
      a_size = atoi(s);
      alphabet = (char*)malloc_or_die(a_size * sizeof(char) * 5);
      break;
    }
  }
  i = 0;
  while(fgets(s, MAX_LINE, replfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
      
    }
    else {
      strcpy(alphabet, s);
      break;
    }
  }

  
  while(fgets(s, MAX_LINE, replfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
      
    }
    else {
      nr_repl_letters = atoi(s);
      replacement_lettersp->nr_rl = nr_repl_letters;
      replacement_lettersp->letters = (struct letter_s*)(malloc_or_die(nr_repl_letters * sizeof(struct letter_s)));
      replacement_lettersp->probs = (double*)(malloc_or_die(nr_repl_letters * a_size * sizeof(double)));
      break;
    }
  }

  cur_letter = 0;
  while (fgets(s, MAX_LINE, replfile) != NULL) {
    if(s[0] == '#' || s[0] == '\n') {
      
    }
    else {
      i = 0;
      while(s[i] != ' ') {
	*((replacement_lettersp->letters + cur_letter)->letter + i) = s[i];
	i++;
      }
      *((replacement_lettersp->letters + 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]);
	a_index = get_alphabet_index(&le, alphabet, a_size);
	*(replacement_lettersp->probs + 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++;
    }
  }
  free(alphabet);
  //dump_replacement_letters(replacement_lettersp, a_size);
}

void get_replacement_letters_multi(FILE *replfile, struct replacement_letter_multi_s *replacement_lettersp)
{
  int MAX_LINE = 1000;
  char s[1000];
  int i,j,k,l,m;
  int nr_repl_letters;
  int cur_letter;
  struct letter_s le;
  double prob;
  int a_index;
  int done;
  int a_size;
  char *alphabet;
  int nr_alphabets;
  int letterindex;
  

  if(replfile == NULL) {
    replacement_lettersp->nr_alphabets = 0;
    return;
  }    
  replacement_lettersp->nr_rl_1 = 0;
  replacement_lettersp->nr_rl_2 = 0;
  replacement_lettersp->nr_rl_3 = 0;
  replacement_lettersp->nr_rl_4 = 0;
 
  while(fgets(s, MAX_LINE, replfile) != NULL){
    if(s[0] == '\n' || s[0] == '#') {
    }
    else {
      nr_alphabets = atoi(s);
      replacement_lettersp->nr_alphabets = nr_alphabets;
      break;
    }
  }  

  for(l = 0; l < nr_alphabets; l++) {
    while(fgets(s, MAX_LINE, replfile) != NULL){
      if(s[0] == '\n' || s[0] == '#') {
      }
      else {
	a_size = atoi(s);
	alphabet = (char*)malloc_or_die(a_size * sizeof(char) * 5);
	break;
      }
    }
    

    i = 0;
    while(fgets(s, MAX_LINE, replfile) != NULL){
      if(s[0] == '\n' || s[0] == '#') {
	
      }
      else {
	strcpy(alphabet, s);
	break;
      }
    }
    
    while(fgets(s, MAX_LINE, replfile) != NULL){
      if(s[0] == '\n' || s[0] == '#') {
      }
      else {
	if(l == 0) {
	  nr_repl_letters = atoi(s);
	  replacement_lettersp->nr_rl_1 = nr_repl_letters;