readhmm_multialpha.c

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

      hmmp->vertex_emiss_prior_scalers_3 = (double*)(malloc_or_die(hmmp->nr_v * sizeof(double)));
    }
    if(hmmp->nr_alphabets > 3) {
      hmmp->vertex_emiss_prior_scalers_4 = (double*)(malloc_or_die(hmmp->nr_v * sizeof(double)));
    }
  }
  /* nr of transitions */
  if(fgets(s, MAX_LINE, file) != NULL) {
    hmmp->nr_t = atoi(&s[19]);
  }
  /* nr of distribution groups */
  if(fgets(s, MAX_LINE, file) != NULL) {
    hmmp->nr_d = atoi(&s[27]);
    hmmp->distrib_groups = (int*)(malloc_or_die((hmmp->nr_d + hmmp->nr_v) * sizeof(int)));
  }

  /* nr of trans tie groups */
  if(fgets(s, MAX_LINE, file) != NULL) {
    hmmp->nr_ttg = atoi(&s[29]);
    hmmp->trans_tie_groups = (int*)(malloc_or_die((hmmp->nr_t + hmmp->nr_ttg) * sizeof(struct transition_s)));
  }

  /* nr of emission priorfiles */
  if(fgets(s, MAX_LINE, file) != NULL) {
    nr_priorfiles = atoi(&s[29]);
    hmmp->nr_ed = nr_priorfiles;
    emission_priorsp = malloc_or_die(nr_priorfiles * sizeof(struct emission_dirichlet_s));
    hmmp->emission_dirichlets = emission_priorsp;
    hmmp->ed_ps = malloc_or_die(hmmp->nr_v * sizeof(struct emission_dirichlet_s*));
  }
  /* read the emission priorfiles */
  if(read_prior_files_multi(nr_priorfiles, emission_priorsp, hmmp->a_size, file) < 0) {
    printf("Could not read emission priorfiles\n");
    exit(-1);
  }

  if(hmmp->nr_alphabets > 1) {
    /* nr of emission priorfiles */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_priorfiles = atoi(&s[29]);
      hmmp->nr_ed_2 = nr_priorfiles;
      emission_priorsp_2 = malloc_or_die(nr_priorfiles * sizeof(struct emission_dirichlet_s));
      hmmp->emission_dirichlets_2 = emission_priorsp_2;
      hmmp->ed_ps_2 = malloc_or_die(hmmp->nr_v * sizeof(struct emission_dirichlet_s*));
    }
    /* read the emission priorfiles */
    if(read_prior_files_multi(nr_priorfiles, emission_priorsp_2, hmmp->a_size_2, file) < 0) {
      printf("Could not read emission priorfiles\n");
      exit(-1);
    }
  }

  if(hmmp->nr_alphabets > 2) {
    /* nr of emission priorfiles */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_priorfiles = atoi(&s[29]);
      hmmp->nr_ed_3 = nr_priorfiles;
      emission_priorsp_3 = malloc_or_die(nr_priorfiles * sizeof(struct emission_dirichlet_s));
      hmmp->emission_dirichlets_3 = emission_priorsp_3;
      hmmp->ed_ps_3 = malloc_or_die(hmmp->nr_v * sizeof(struct emission_dirichlet_s*));
    }
    /* read the emission priorfiles */
    if(read_prior_files_multi(nr_priorfiles, emission_priorsp_3, hmmp->a_size_3, file) < 0) {
      printf("Could not read emission priorfiles\n");
      exit(-1);
    }
  }

  if(hmmp->nr_alphabets > 3) {
    /* nr of emission priorfiles */
    if(fgets(s, MAX_LINE, file) != NULL) {
      nr_priorfiles = atoi(&s[29]);
      hmmp->nr_ed_4 = nr_priorfiles;
      emission_priorsp_4 = malloc_or_die(nr_priorfiles * sizeof(struct emission_dirichlet_s));
      hmmp->emission_dirichlets_4 = emission_priorsp_4;
      hmmp->ed_ps_4 = malloc_or_die(hmmp->nr_v * sizeof(struct emission_dirichlet_s*));
    }
    /* read the emission priorfiles */
    if(read_prior_files_multi(nr_priorfiles, emission_priorsp_4, hmmp->a_size_4, file) < 0) {
      printf("Could not read emission priorfiles\n");
      exit(-1);
    }
  }

   /* nr of transition priorfiles */
  if(fgets(s, MAX_LINE, file) != NULL) {
    nr_trans_priorfiles = atoi(&s[29]);
    transition_priorsp = NULL;
    /* not implemented yet */
  }
  /* read the transition priorfiles */
  if(read_trans_prior_files_multi(nr_trans_priorfiles, transition_priorsp, hmmp, file) < 0) {
    printf("Could not read transition priorfiles\n");
    exit(-1);
  }

  /* empty row */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  /* empty row */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  /* reads ****************Modules*****************/
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  /* read the modules */
  silent_counter = 0;
  locked_counter = 0;
  for(i = 0; i < hmmp->nr_m; i++) {
    *(hmmp->modules + i) = module;
    if((res = read_module_multi(s, file, hmmp, module, &silent_counter, &locked_counter)) < 0) {
      printf("Could not read modules\n");
      exit(-1);
    }
    module++;
  }

  *(hmmp->silent_vertices + silent_counter) = END;
  *(hmmp->locked_vertices + hmmp->nr_v) = END;
#ifdef DEBUG_RD
  //dump_locked_vertices(hmmp);
  //dump_silent_vertices(hmmp);
  //dump_multi_modules(hmmp);
#endif
  
  /* empty row */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  
  /* reads ****************Emission distribution groups*****************/
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  

  /* read the distribution groups */
  cur = hmmp->distrib_groups;
  for(i = 0; i < hmmp->nr_d; i++) {
    j = 0;
    if(fgets(s, MAX_LINE, file) != NULL) {
      while(1) {
	if(s[j] == ':') {
	  break;
	}
	j++;
      }
      j++;
      j++;
      while(1) {
	*cur = atoi(&s[j]);
	cur++;
	while(s[j] != ' ' && s[j] != '\n') {
	  j++;
	}
	while(s[j] == ' ') {
	  j++;
	}
	if(s[j] == '\n') {
	  break;
	}
      }
      *cur = END;
      cur++;
    }
  }
  
  /* empty row */
  if(fgets(s, MAX_LINE, file) != NULL) {
  }

  /* reads ****************Transition tie groups*****************/
  if(fgets(s, MAX_LINE, file) != NULL) {
  }
  /* read the trans tie groups */
  trans_ties = hmmp->trans_tie_groups;
  for(i = 0; i < hmmp->nr_ttg; i++) {
    if(fgets(s, MAX_LINE, file) != NULL && s[0] != '\n') {
      j = 0;
      while(1) {
	if(s[j] == ':') {
	  break;
	}
	j++;
      }
      j++;
      j++;
      while(1) {
	trans.from_v = atoi(&s[j]);
	while(s[j] != '>') {
	  j++;
	}
	j++;
	trans.to_v = atoi(&s[j]);
	memcpy(trans_ties, &trans, sizeof(struct transition_s));
	trans_ties++;
	while(s[j] != ' ' && s[j] != '\n') {
	  j++;
	}
	while(s[j] == ' ') {
	  j++;
	}
	if(s[j] == '\n') {
	  break;
	}
      }
      trans.to_v = END;
      trans.from_v = END;
      memcpy(trans_ties, &trans, sizeof(struct transition_s));
      trans_ties++;
    }
    else {
      hmmp->nr_ttg = i;
      break;
    }
  }
#ifdef DEBUG_RD
  dump_distrib_groups(hmmp->distrib_groups, hmmp->nr_d);
  dump_trans_tie_groups(hmmp->trans_tie_groups, hmmp->nr_ttg);
#endif
  /* create to_silent_trans_array */
  create_to_silent_trans_array_multi(hmmp);
  
  /* create from_trans_array */
  create_from_trans_array_multi(hmmp);
 
  /* create to_trans_array */
  create_to_trans_array_multi(hmmp);
 
  /* create tot_transitions */
  create_tot_transitions_multi(hmmp);
 
  /* create tot_to_trans_array and tot_from_trans_arrays*/
  create_tot_trans_arrays_multi(hmmp);

  /* get the set of labels and the number of labels */
  get_set_of_labels_multi(hmmp);

#ifdef DEBUG_RD 
  dump_emiss_matrix(hmmp->nr_v, hmmp->a_size, hmmp->emissions);
  dump_trans_matrix(hmmp->nr_v, hmmp->nr_v, hmmp->transitions);
  printf("hmmp->emission_dirichlets = %x\n", hmmp->emission_dirichlets);
  for(i = 0; i < hmmp->nr_v; i++) {
    printf("hmmp->ed_ps for vertex %d = %x\n", i, *(hmmp->ed_ps + i));
  }
#endif

  /* make sure all probabilities are legal*/
  //check_probs_multi(hmmp);

  if(verbose == YES) {
    printf("done\n");
  }
}


/************************read_module_multi *************************************/
int read_module_multi(char *s, FILE *file, struct hmm_multi_s *hmmp, struct module_multi_s *modulep,
		int *silent_counter, int *locked_counter)
{
  int nr_v, nr_t, nr_e, nr_e2, nr_e3, nr_e4, nr_et;
  int i,j,k;
  int from_v, to_v;
  double prob, log_prob;
  char type[50];
  char prifile_name[500], prifile_name_2[500], prifile_name_3[500], prifile_name_4[500];
  char *p, *probp;
  struct emission_dirichlet_s *priorsp, *priorsp_2, *priorsp_3, *priorsp_4;
  int silent_vertex;

  //i = 0;

  /* module name */
  if(fgets(s, MAX_LINE, file) != NULL) {
    strcpy(modulep->name, &s[8]);
#ifdef DEBUG_RD
    printf("module name %s", s);
#endif 
  }
  /* module type */
  if(fgets(s, MAX_LINE, file) != NULL) {
    strcpy(type, &s[6]);
    if(strncmp(type, "Singlenode", 10) == 0) {
      modulep->type = SINGLENODE;
    }
    else if(strncmp(type, "Cluster", 7) == 0) {
      modulep->type = CLUSTER;
    }
    else if(strncmp(type, "Forward_std", 11) == 0) {
      modulep->type = FORWARD_STD;
    }
    else if(strncmp(type, "Forward_alt", 11) == 0) {
      modulep->type = FORWARD_ALT;
    }
    else if(strncmp(type, "Singleloop", 10) == 0) {
      modulep->type = SINGLELOOP;
    }
    else if(strncmp(type, "Profile7", 8) == 0) {
      modulep->type = PROFILE7;
    }
    else if(strncmp(type, "Profile9", 8) == 0) {
      modulep->type = PROFILE9;
    }
    else {
      printf("Error: module is of unknown type\n");
      exit(-1);
    }
   
  }

  /* nr vertices */
  if(fgets(s, MAX_LINE, file) != NULL) {
    modulep->nr_v = atoi(&s[12]);
    modulep->vertices = (int*)(malloc_or_die(modulep->nr_v * sizeof(int)));
  }
  /* emission prior file */
  if(fgets(s, MAX_LINE, file) != NULL) {
    strcpy(prifile_name, (&s[23]));
    if((p = strstr(prifile_name, "\n")) != NULL) {
      *p = '\0';
    }
    if(strncmp(prifile_name, "null", 4) == 0) {
      strcpy(modulep->priorfile_name, "null");
      priorsp = NULL;
    }
    else { 
      strcpy(modulep->priorfile_name, prifile_name);
      strcat(modulep->priorfile_name, "\0");
      for(i = 0; i < hmmp->nr_ed; i++) {
	priorsp = (hmmp->emission_dirichlets + i);
	if((strncmp(prifile_name, priorsp->name, 200)) == 0) {
	  /* keep this priorsp */
	  break;
	}
	else {
#ifdef DEBUG_RD
	  printf("prifile_name = %s\n", prifile_name);
	  printf("priorsp->name = %s\n", priorsp->name);
#endif
	}
	if(i == hmmp->nr_ed-1) /* no name equals priorfile_name */{
	  printf("Couldn't find emission priorfile '%s'\n", prifile_name);
	  exit(-1);
	}
      }
    }
  }

  if(hmmp->nr_alphabets > 1) {
    /* emission prior file */
    if(fgets(s, MAX_LINE, file) != NULL) {
      strcpy(prifile_name_2, (&s[23]));
      if((p = strstr(prifile_name_2, "\n")) != NULL) {
	*p = '\0';
      }
      if(strncmp(prifile_name_2, "null", 4) == 0) {
	strcpy(modulep->priorfile_name_2, "null");
	priorsp_2 = NULL;
      }
      else { 
	strcpy(modulep->priorfile_name_2, prifile_name_2);
	strcat(modulep->priorfile_name_2, "\0");
	for(i = 0; i < hmmp->nr_ed_2; i++) {
	  priorsp_2 = (hmmp->emission_dirichlets_2 + i);
	  if((strncmp(prifile_name_2, priorsp_2->name, 200)) == 0) {
	    /* keep this priorsp */
	    break;
	  }
	  else {
#ifdef DEBUG_RD
	    printf("prifile_name_2 = %s\n", prifile_name_2);
	    printf("priorsp_2->name = %s\n", priorsp_2->name);
#endif
	  }
	  if(i == hmmp->nr_ed_2 - 1) /* no name equals priorfile_name */{
	    printf("Couldn't find emission priorfile '%s'\n", prifile_name_2);
	    exit(-1);
	  }
	}
      }
    }
  }
  
  if(hmmp->nr_alphabets > 2) {
    /* emission prior file */
    if(fgets(s, MAX_LINE, file) != NULL) {
      strcpy(prifile_name_3, (&s[23]));
      if((p = strstr(prifile_name_3, "\n")) != NULL) {
	*p = '\0';
      }
      if(strncmp(prifile_name_3, "null", 4) == 0) {
	strcpy(modulep->priorfile_name_3, "null");
	priorsp_3 = NULL;
      }
      else { 
	strcpy(modulep->priorfile_name_3, prifile_name_3);
	strcat(modulep->priorfile_name_3, "\0");
	for(i = 0; i < hmmp->nr_ed_3; i++) {
	  priorsp_3 = (hmmp->emission_dirichlets_3 + i);
	  if((strncmp(prifile_name_3, priorsp_3->name, 200)) == 0) {
	    /* keep this priorsp */
	    break;
	  }
	  else {
#ifdef DEBUG_RD
	    printf("prifile_name_3 = %s\n", prifile_name_3);
	    printf("priorsp_3->name = %s\n", priorsp_3->name);
#endif
	  }
	  if(i == hmmp->nr_ed_3 - 1) /* no name equals priorfile_name */{
	    printf("Couldn't find emission priorfile '%s'\n", prifile_name_3);
	    exit(-1);
	  }
	}
      }
    }
  }

  if(hmmp->nr_alphabets > 3) {
    /* emission prior file */
    if(fgets(s, MAX_LINE, file) != NULL) {
      strcpy(prifile_name_4, (&s[23]));
      if((p = strstr(prifile_name_4, "\n")) != NULL) {
	*p = '\0';
      }
      if(strncmp(prifile_name_4, "null", 4) == 0) {
	strcpy(modulep->priorfile_name_4, "null");
	priorsp_4 = NULL;
      }
      else { 
	strcpy(modulep->priorfile_name_4, prifile_name_4);
	strcat(modulep->priorfile_name_4, "\0");
	for(i = 0; i < hmmp->nr_ed_4; i++) {
	  priorsp_4 = (hmmp->emission_dirichlets_4 + i);
	  if((strncmp(prifile_name_4, priorsp_4->name, 200)) == 0) {
	    /* keep this priorsp */
	    break;
	  }
	  else {
#ifdef DEBUG_RD
	    printf("prifile_name_4 = %s\n", prifile_name_4);
	    printf("priorsp_4->name = %s\n", priorsp_4->name);
#endif
	  }
	  if(i == hmmp->nr_ed_4 - 1) /* no name equals priorfile_name */{
	    printf("Couldn't find emission priorfile '%s'\n", prifile_name_4);
	    exit(-1);
	  }
	}
      }
    }