readhmm_multialpha.c

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

  }

  /* transition prior file */
  if(fgets(s, MAX_LINE, file) != NULL) {
    /* not implemented yet */
  }
  /* empty line */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }

#ifdef DEBUG_RD
  printf("modulep->nr_v = %d\n", modulep->nr_v);
#endif  
/* loop over the vertices */
   for(i = 0; i < modulep->nr_v; i++) {
    /* Vertex nr */
    if(fgets(s, MAX_LINE, file) != NULL) {
#ifdef DEBUG_RD      
      printf("Vertex nr: %s", s);
#endif
      from_v = atoi(&s[7]);
      *(modulep->vertices + i) = from_v;
      /* connect this vertex to its priorfile */
      *(hmmp->ed_ps + from_v) = priorsp;
      if(hmmp->nr_alphabets > 1) {
	*(hmmp->ed_ps_2 + from_v) = priorsp_2;
      }
      if(hmmp->nr_alphabets > 2) {
	*(hmmp->ed_ps_3 + from_v) = priorsp_3;
      }
      if(hmmp->nr_alphabets > 3) {
	*(hmmp->ed_ps_4 + from_v) = priorsp_4;
      }
    }
    else {
      printf("Error in hmm spec\n");
      exit(0);
    }
    /* Vertex type */
    if(fgets(s, MAX_LINE, file) != NULL) {
#ifdef DEBUG_RD
      printf("Vertex type: %s", s);
#endif
      strcpy(type, &s[13]);
      if(modulep->type == PROFILE7 || modulep->type == PROFILE9) {
	modulep->v_type = PROFILEV;
      }
      if(strncmp(type, "standard", 8) == 0) {
	if(modulep->type != PROFILE7 && modulep->type != PROFILE9) {
	  modulep->v_type = STANDARDV;
	  silent_vertex = NO;
	}
      }
      else if(strncmp(type, "silent", 6) == 0) {
	silent_vertex = YES;
	if(modulep->type != PROFILE7 && modulep->type != PROFILE9) {
	  modulep->v_type = SILENTV;
	}
	*(hmmp->silent_vertices + *silent_counter) = from_v;
	*silent_counter = *silent_counter + 1;
      }
      else if(strncmp(type, "locked", 5) == 0) {
	modulep->v_type = LOCKEDV;
	*(hmmp->locked_vertices + from_v) = YES;
	*locked_counter = *locked_counter + 1;
	silent_vertex = NO;
      }
      else if(strncmp(type, "start", 5) == 0) {
	modulep->v_type = STARTV;
	hmmp->startnode = from_v;
      }
      else if(strncmp(type, "end", 3) == 0) {
	modulep->v_type = ENDV;
	hmmp->endnode = from_v;
      }
      else {
	printf("Error: vertex type is undefined\n");
	printf("vertex type = %s\n", type);
	printf("from_v = %d\n", from_v);
	exit(-1);
      }
    }
    /* Vertex label */
    if(fgets(s, MAX_LINE, file) != NULL) {
      *(hmmp->vertex_labels + from_v) = s[14];
    }
    /* transition prior scaler */
    if(fgets(s, MAX_LINE, file) != NULL) {
      *(hmmp->vertex_trans_prior_scalers + from_v) = atof(&(s[25]));
    }
    /* emission prior scaler */
    if(fgets(s, MAX_LINE, file) != NULL) {
      *(hmmp->vertex_emiss_prior_scalers + from_v) = atof(&(s[25]));
    }
    if(hmmp->nr_alphabets > 1) {
      /* emission prior scaler */
      if(fgets(s, MAX_LINE, file) != NULL) {
	*(hmmp->vertex_emiss_prior_scalers_2 + from_v) = atof(&(s[25]));
      }
    }
    if(hmmp->nr_alphabets > 2) {
      /* emission prior scaler */
      if(fgets(s, MAX_LINE, file) != NULL) {
	*(hmmp->vertex_emiss_prior_scalers_3 + from_v) = atof(&(s[25]));
      }
    }
    if(hmmp->nr_alphabets > 3) {
      /* emission prior scaler */
      if(fgets(s, MAX_LINE, file) != NULL) {
	*(hmmp->vertex_emiss_prior_scalers_4 + from_v) = atof(&(s[25]));
      }
    }
    /* Nr transitions */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_t = atoi(&s[17]);
    }
    /* Nr end transitions */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_et = atoi(&s[21]);
    }
    /* Nr emissions */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_e = atoi(&s[17]);
    }
    if(hmmp->nr_alphabets > 1) {
      /* Nr emissions */
      if(fgets(s, MAX_LINE, file) != NULL) {
	nr_e2 = atoi(&s[17]);
      }
    }
    if(hmmp->nr_alphabets > 2) {
      /* Nr emissions */
      if(fgets(s, MAX_LINE, file) != NULL) {
	nr_e3 = atoi(&s[17]);
      }
    }
    if(hmmp->nr_alphabets > 3) {
      /* Nr emissions */
      if(fgets(s, MAX_LINE, file) != NULL) {
	nr_e4 = atoi(&s[17]);
      }
    }
    /* read transition probabilities */
    fgets(s, MAX_LINE, file);
    for(j = 0; j < nr_t; j++) {
      if(fgets(s, MAX_LINE, file) != NULL) {
	to_v = atoi(&s[8]);
	if(to_v < 10 ) {
	  prob = (double)(atof(&s[11]));
	}
	else if(to_v < 100) {
	  prob = (double)(atof(&s[12]));
	 }
	else if(to_v < 1000) {
	  prob = (double)(atof(&s[13]));
	}
	else if(to_v < 10000) {
	  prob = (double)(atof(&s[14]));
	 }
	else {
	  printf("Sorry, reader cannot handle HMMs with more than 10000 states\n");
	   exit(0);
	}
	if(prob != 0.0) {
	  log_prob = log10(prob);
	}
	else {
	  log_prob = DEFAULT;
	}
#ifdef DEBUG_RD
	printf("prob from %d to %d = %f\n", from_v, to_v, prob);
#endif
	*(hmmp->transitions + get_mtx_index(from_v, to_v, hmmp->nr_v)) = prob;
	*(hmmp->log_transitions + get_mtx_index(from_v, to_v, hmmp->nr_v)) = log_prob;
      }
    }
    
    /* read end transition probabilities */
    fgets(s, MAX_LINE, file);
    for(j = 0; j < nr_et; j++) {
      if(fgets(s, MAX_LINE, file) != NULL) {
	to_v = atoi(&s[8]);
	if(to_v < 10 ) {
	  prob = (double)(atof(&s[11]));
	}
	else if(to_v < 100) {
	  prob = (double)(atof(&s[12]));
	}
	else if(to_v < 1000) {
	  prob = (double)(atof(&s[13]));
	}
	else if(to_v < 10000) {
	  prob = (double)(atof(&s[14]));
	}
	else {
	  printf("Sorry, reader cannot handle HMMs with more than 10000 states\n");
	  exit(0);
	}
	if(prob != 0.0) {
	  log_prob = log10(prob);
	}
	else {
	  log_prob = DEFAULT;
	}
#ifdef DEBUG_RD
	printf("prob from %d to %d = %f\n", from_v, to_v, prob);
#endif
	*(hmmp->transitions + get_mtx_index(from_v, to_v, hmmp->nr_v)) = prob;
	*(hmmp->log_transitions + get_mtx_index(from_v, to_v, hmmp->nr_v)) = log_prob;
      }
    }
    /* read emission probabilities */
    fgets(s, MAX_LINE, file);
    for(j = 0; j < nr_e; j++) {
      if(fgets(s, MAX_LINE, file) != NULL) {
#ifdef DEBUG_RD
	printf("%s", s);
#endif
	k = 0;
	while(s[k] != ' ') {
	  k++;
	}
	if(k > 7) {
	  printf("Cannot read hmm file, please check hmm specification\n");
	}
	k++;
	prob = (double)(atof(&s[k]));
	if(prob != 0.0) {
	  log_prob = log10(prob);
	}
	else {
	  log_prob = DEFAULT;
	}
	if(silent_vertex == YES) {
	  prob = SILENT;
	  log_prob = SILENT;
	}
	*(hmmp->emissions + get_mtx_index(from_v, j, hmmp->a_size)) = prob;
	*(hmmp->log_emissions + get_mtx_index(from_v, j, hmmp->a_size)) = log_prob;
      }
    }
    fgets(s, MAX_LINE, file);
    
    if(hmmp->nr_alphabets > 1) {
      fgets(s, MAX_LINE, file);
      for(j = 0; j < nr_e2; j++) {
	if(fgets(s, MAX_LINE, file) != NULL) {
	  k = 0;
	  while(s[k] != ' ') {
	    k++;
	  }
	  if(k > 7) {
	    printf("Cannot read hmm file, please check hmm specification\n");
	  }
	  k++;
	  prob = (double)(atof(&s[k]));
	  if(prob != 0.0) {
	    log_prob = log10(prob);
	  }
	  else {
	    log_prob = DEFAULT;
	  }
	  if(silent_vertex == YES) {
	    prob = SILENT;
	    log_prob = SILENT;
	  }
	  *(hmmp->emissions_2 + get_mtx_index(from_v, j, hmmp->a_size_2)) = prob;
	  *(hmmp->log_emissions_2 + get_mtx_index(from_v, j, hmmp->a_size_2)) = log_prob;
	}
      }
      fgets(s, MAX_LINE, file);
    }
    if(hmmp->nr_alphabets > 2) {
      fgets(s, MAX_LINE, file);
      for(j = 0; j < nr_e3; j++) {
	if(fgets(s, MAX_LINE, file) != NULL) {
	  k = 0;
	  while(s[k] != ' ') {
	    k++;
	  }
	  if(k > 7) {
	    printf("Cannot read hmm file, please check hmm specification\n");
	  }
	  k++;
	  prob = (double)(atof(&s[k]));
	  if(prob != 0.0) {
	    log_prob = log10(prob);
	  }
	  else {
	    log_prob = DEFAULT;
	  }
	  if(silent_vertex == YES) {
	    prob = SILENT;
	    log_prob = SILENT;
	  }
	  *(hmmp->emissions_3 + get_mtx_index(from_v, j, hmmp->a_size_3)) = prob;
	  *(hmmp->log_emissions_3 + get_mtx_index(from_v, j, hmmp->a_size_3)) = log_prob;
	}
      }
      fgets(s, MAX_LINE, file);
    }

    if(hmmp->nr_alphabets > 3) {
      fgets(s, MAX_LINE, file);
      for(j = 0; j < nr_e4; j++) {
	if(fgets(s, MAX_LINE, file) != NULL) {
	  k = 0;
	  while(s[k] != ' ') {
	    k++;
	  }
	  if(k > 7) {
	    printf("Cannot read hmm file, please check hmm specification\n");
	  }
	  k++;
	  prob = (double)(atof(&s[k]));
	  if(prob != 0.0) {
	    log_prob = log10(prob);
	  }
	  else {
	    log_prob = DEFAULT;
	  }
	  if(silent_vertex == YES) {
	    prob = SILENT;
	    log_prob = SILENT;
	  }
	  
	  *(hmmp->emissions_4 + get_mtx_index(from_v, j, hmmp->a_size_4)) = prob;
	  *(hmmp->log_emissions_4 + get_mtx_index(from_v, j, hmmp->a_size_4)) = log_prob;
	}
      }
      fgets(s, MAX_LINE, file);
    }
    
    silent_vertex = NO;
  }
  
  /* read ---------------------------------------- */
  fgets(s, MAX_LINE, file);
  return 0;
  
}


void create_to_silent_trans_array_multi(struct hmm_multi_s *hmmp)
{
  int v,w;
  int malloc_size;
  int *values;

  malloc_size = 0;
  for(v = 0; v < hmmp->nr_v; v++) {
    for(w = 0; w < hmmp->nr_v; w++) {
      if(*(hmmp->transitions + get_mtx_index(v,w,hmmp->nr_v)) != 0 && silent_vertex_multi(w,hmmp) == YES) {
	malloc_size++;
      }
    }
    malloc_size++;
  }
  
  hmmp->to_silent_trans_array = (int**)malloc_or_die(hmmp->nr_v * sizeof(int*) + malloc_size * sizeof(int));
  values = (int*)(hmmp->to_silent_trans_array + hmmp->nr_v);

  for(v = 0; v < hmmp->nr_v; v++) {
    *(hmmp->to_silent_trans_array + v) = values;
    for(w = 0; w < hmmp->nr_v; w++) {
      if(*(hmmp->transitions + get_mtx_index(v,w,hmmp->nr_v)) != 0 && silent_vertex_multi(w,hmmp) == YES) {
	*values = w;
	values++;
      }
    }
    *values = END;
    values++;
  }

#ifdef DEBUG_RD
  dump_to_silent_trans_array(hmmp->nr_v, hmmp->to_silent_trans_array);
#endif
}

/* Go through transmission matrix and get all probabilities that are not 0
 * into from_trans_array */
void create_from_trans_array_multi(struct hmm_multi_s *hmmp)
{
  int v,w,*xp;
  int has_to_trans;
  int array_head_size, array_tail_size;
  struct path_element **from_trans_array, *from_trans, *temp_path;

  array_tail_size = 0;
  array_head_size = hmmp->nr_v;

  /* estimate how much space we need to store transitions */
  array_tail_size = (hmmp->nr_t/hmmp->nr_v + 3 + MAX_GAP_SIZE) * MAX_GAP_SIZE/2 * hmmp->nr_v;
  
#ifdef DEBUG_RD
  printf("array_head_size, array_tail_size = %d, %d\n", array_head_size, array_tail_size);
#endif  
  from_trans_array = (struct path_element**)
    (malloc_or_die(array_head_size * sizeof(struct path_element*) +
		   (array_tail_size + hmmp->nr_v) * sizeof(struct path_element)));
  from_trans = (struct path_element*)(from_trans_array + hmmp->nr_v);
  hmmp->from_trans_array = from_trans_array;
  
  /* find all paths and add them to from_trans_array */
  for(v = 0; v < hmmp->nr_v; v++) /* to-vertex */ {
    *from_trans_array = from_trans;
    if(silent_vertex_multi(v, hmmp) == YES) {
      from_trans->vertex = END;
      from_trans->next = NULL;
      from_trans++;
      from_trans_array++;
      continue;
    }
    for(w = 0; w < hmmp->nr_v; w++) /* from-vertex */ {
      if(silent_vertex_multi(w,hmmp) == YES) {
	continue;
      }
      temp_path = (struct path_element*)(malloc_or_die(1000 * sizeof(struct path_element)));
      add_all_from_paths_multi(w, v, hmmp, &from_trans, temp_path, 0);
      free(temp_path);
    }
    from_trans->vertex = END;
    from_trans->next = NULL;
    from_trans++;
    from_trans_array++;
  }
#ifdef DEBUG_RD
  dump_trans_matrix(hmmp->nr_v, hmmp->nr_v, hmmp->transitions);
  dump_from_trans_array(hmmp->nr_v, hmmp->from_trans_array);
#endif
}

void add_all_from_paths_multi(int v, int w, struct hmm_multi_s *hmmp,
			struct path_element **from_transp, struct path_element *temp_pathp,
			int length)
{
  int i,j;
  int *xp;
  struct path_element p_el, cur_p_el;
  
  if(length > MAX_GAP_SIZE) {
    return;
  }

  cur_p_el.vertex = v;
  cur_p_el.next = NULL;
  memcpy(temp_pathp + length, &cur_p_el, sizeof(struct path_element));

  if(*(hmmp->transitions + get_mtx_index(v, w, hmmp->nr_v)) != 0.0) {