readhmm.c

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

	  printf("Couldn't find emission priorfile '%s'\n", prifile_name);
	  exit(-1);
	}
      }
    }
  }
  /* transition prior file */
  if(fgets(s, MAX_LINE, file) != NULL) {
    /* not implemented yet */
  }
  /* empty line */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  /* loop over the vertices */
  for(i = 0; i < modulep->nr_v; i++) {
    /* Vertex nr */
    if(fgets(s, MAX_LINE, file) != NULL) {
      from_v = atoi(&s[7]);
      *(modulep->vertices + i) = from_v;

      /* connect this vertex to its priorfile */
      *(hmmp->ed_ps + from_v) = priorsp;
    }
    /* Vertex type */
    if(fgets(s, MAX_LINE, file) != NULL) {
      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[23]));
    }
    /* 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[15]);
    }
    /* 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("end 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) {
      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;
	}
#ifdef DEBUG_RD
	printf("emissprob in %d of letter %d = %f\n", from_v, j, prob);
#endif	
	*(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);
    silent_vertex = NO;
  }
  /* read ---------------------------------------- */
  fgets(s, MAX_LINE, file);
  
#ifdef DEBUG_RD
  printf("exiting read_module\n");
#endif	
  return 0;
  
}


void create_to_silent_trans_array(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(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(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(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(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(w,hmmp) == YES) {
	continue;
      }
      temp_path = (struct path_element*)(malloc_or_die(1000 * sizeof(struct path_element)));
      add_all_from_paths(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(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) {
#ifdef DEBUG_RD
    printf("adding path: ");
#endif
    /* direct path to w, add total path */
    for(i = 0; i < length; i++) {
      p_el.vertex = (temp_pathp + i)->vertex;
#ifdef DEBUG_RD
      printf("%d ", p_el.vertex);
#endif
      p_el.next = (*from_transp) + 1;
      memcpy(*from_transp, &p_el, sizeof(struct path_element));
      (*from_transp)++;
    }
    memcpy(*from_transp, &cur_p_el, sizeof(struct path_element));
#ifdef DEBUG_RD
    printf("%d %d\n", cur_p_el.vertex, w);
#endif    
    (*from_transp)++;
  }
  xp = *(hmmp->to_silent_trans_array + v);
  while(*xp != END) {
    add_all_from_paths(*xp, w, hmmp, from_transp, temp_pathp, length + 1);
    xp++;
  }

}


/* Go through transmission matrix and get all probabilities that are not 0
 * into to_trans_array */
void create_to_trans_array(struct hmm_multi_s *hmmp)
{
  int v,w,*xp;
  int has_to_trans;
  int array_head_size, array_tail_size;
  struct path_element **to_trans_array, *to_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_tail_size = %d\n", array_tail_size);
#endif
  to_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)));
  to_trans = (struct path_element*)(to_trans_array + hmmp->nr_v);
  hmmp->to_trans_array = to_trans_array;
  
  /* find all paths and add them to to_trans_array */
  for(v = 0; v < hmmp->nr_v; v++) /* from-vertex */ {
    *to_trans_array = to_trans;
    if(silent_vertex(v, hmmp) == YES) {
      to_trans->vertex = END;
      to_trans->next = NULL;
      to_trans++;
      to_trans_array++;
      continue;
    }
    for(w = 0; w < hmmp->nr_v; w++) /* to-vertex */ {
      if(silent_vertex(w,hmmp) == YES) {
	continue;
      }
      temp_path = (struct path_element*)(malloc_or_die(1000 * sizeof(struct path_element)));
      add_all_to_paths(v, w, hmmp, &to_trans, temp_path, 0);
      free(temp_path);
    }
    to_trans->vertex = END;
    to_trans->next = NULL;