std_calculation_funcs.c

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

    }
    prf_column_length = sqrt(prf_column_length);
    subst_mtx_row_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      subst_mtx_row_length += pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)),2);
    }
    subst_mtx_row_length = sqrt(subst_mtx_row_length);
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->share) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size))/ 
	(prf_column_length * subst_mtx_row_length);
    }
  }
}

void add_to_E_subst_mtx_dot_product_nr_occ(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
					   int k, int a_size, int normalize, double *subst_mtx, char *alphabet)
{
  /* k is a state index, p is a sequence position index */
  
  int a_index;
  double prf_column_length, subst_mtx_row_length;
  int query_index;
  
  query_index = get_alphabet_index_msa_query((msa_seq + (p * (a_size+1)))->query_letter, alphabet, a_size);
  if(normalize == NO) {
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));
    }
  }
  else {
    prf_column_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      prf_column_length += pow(((msa_seq + get_mtx_index(p, a_index, a_size+1))->share),2);
    }
    prf_column_length = sqrt(prf_column_length);
    subst_mtx_row_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      subst_mtx_row_length += pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)),2);
    }
    subst_mtx_row_length = sqrt(subst_mtx_row_length);
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size))/ 
	(prf_column_length * subst_mtx_row_length);
    }
  }
}

void add_to_E_subst_mtx_dot_product_prior(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
					  int k, int a_size, int normalize, double *subst_mtx, char *alphabet)
{
  /* k is a state index, p is a sequence position index */
  
  int a_index;
  double prf_column_length, subst_mtx_row_length;
  int query_index;
  double rli;
  
  rli = REST_LETTER_INDEX;
  
  query_index = get_alphabet_index_msa_query((msa_seq + (p * (a_size+1)))->query_letter, alphabet, a_size);
  if(normalize == NO) {
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->share) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));

      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->share) * rli;
    }
    
  }
  else {
    prf_column_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      prf_column_length += pow(((msa_seq + get_mtx_index(p, a_index, a_size+1))->share),2);
    }
    prf_column_length = sqrt(prf_column_length);
    subst_mtx_row_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      subst_mtx_row_length += pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)),2);
    }
    subst_mtx_row_length = sqrt(subst_mtx_row_length);
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->share) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size))/ 
	(prf_column_length * subst_mtx_row_length);

      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->share) * rli / 
	(prf_column_length * subst_mtx_row_length);
    }
  }
}


void add_to_E_subst_mtx_dot_product_prior_nr_occ(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
						 int k, int a_size, int normalize, double *subst_mtx, char *alphabet)
{
  /* k is a state index, p is a sequence position index */
  
  int a_index;
  double prf_column_length, subst_mtx_row_length;
  int query_index;
  double rli;
  
  rli = REST_LETTER_INDEX;
  
  query_index = get_alphabet_index_msa_query((msa_seq + (p * (a_size+1)))->query_letter, alphabet, a_size);
  if(normalize == NO) {
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));

      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) * rli;
    }
    
  }
  else {
    prf_column_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      prf_column_length += pow(((msa_seq + get_mtx_index(p, a_index, a_size+1))->share),2);
    }
    prf_column_length = sqrt(prf_column_length);
    subst_mtx_row_length = 0.0;
    for(a_index = 0; a_index < a_size; a_index++) {
      subst_mtx_row_length += pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)),2);
    }
    subst_mtx_row_length = sqrt(subst_mtx_row_length);
    for(a_index = 0; a_index < a_size; a_index++) {
      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size))/ 
	(prf_column_length * subst_mtx_row_length);

      *(E + get_mtx_index(k, a_index, a_size)) += Eka_base *
	(double)((msa_seq + get_mtx_index(p, a_index, a_size+1))->nr_occurences) * rli / 
	(prf_column_length * subst_mtx_row_length);
    }
  }
}



/* General versions of the functions needed for keeping track of the labeleings in the one-best algorithm */

void update_labelings(struct one_best_s *cur_rowp, char *vertex_labels, 
		      int *sorted_v_list, int seq_len, int c, char *labels, int nr_of_labels, int nr_v)
{
  int v,w;
  int v_list_index;
  int cur_address;
  int first;
  char *tmp_labeling;
  char cur_label;
  int **same_labeling_lists;
  int *same_labeling_list_indices;
  int i;
  

#ifdef DEBUG_LABELING_UPDATE
  dump_v_list(sorted_v_list);
  printf("nr of labels = %d\n", nr_of_labels);
  printf("dump of nr three\n");
#endif

  v_list_index = 0;

  same_labeling_list_indices = (int*)(malloc_or_die(nr_of_labels * sizeof(int)));
  
  same_labeling_lists = (int**)(malloc_or_die(nr_of_labels * sizeof(int*)));


  for(i = 0; i < nr_of_labels; i++) {
    same_labeling_lists[i] = (int*)(malloc_or_die((nr_v * 2 + 1) * sizeof(int)));
  }
 

  for(v = 0; v < nr_v; v++) {
    (cur_rowp + v)->is_updated = NO;
  }
  
  for(i = 0; i < nr_of_labels; i++) {
    same_labeling_list_indices[i] = 0;
  }


  for(v = 1; v < nr_v-1; v++) {
    /* find all states with same labeling as this state up to previous row */
    
    if((cur_rowp + v)->is_updated == NO && (cur_rowp + v)->labeling != NULL) {
      cur_address = (int)((cur_rowp+v)->labeling);
#ifdef DEBUG_LABELING_UPDATE
      printf("searching vertex %d\n", v);
#endif
      for(i = 0; i < nr_of_labels; i++) {
	if(*(vertex_labels + v) == *(labels + i)) {
	  *(*(same_labeling_lists + i) + same_labeling_list_indices[i]) = v;
	  same_labeling_list_indices[i] += 1;
	  break;
	}
      }
      
      (cur_rowp+v)->is_updated = YES;
      for(w = v+1; w < nr_v-1; w++) {
	if((int)((cur_rowp+w)->labeling) == cur_address && (cur_rowp + w)->is_updated == NO) {
#ifdef DEBUG_LABELING_UPDATE	  
	  printf("found same address, vertex nr = %d\n", w);
#endif
	  for(i = 0; i < nr_of_labels; i++) {
	    if(*(vertex_labels + w) == *(labels + i)) {
	      *(*(same_labeling_lists + i) + same_labeling_list_indices[i]) = w;
	      same_labeling_list_indices[i] += 1;
	      break;
	    }
	  }

	  (cur_rowp+w)->is_updated = YES;
	}
      }

      for(i = 0; i < nr_of_labels; i++) {
	*(*(same_labeling_lists + i) + same_labeling_list_indices[i]) = END;
	same_labeling_list_indices[i] += 1;
      }
    }
  }
  for(i = 0; i < nr_of_labels; i++) {
    *(*(same_labeling_lists + i) + same_labeling_list_indices[i]) = TOT_END;
  }


#ifdef DEBUG_LABELING_UPDATE
  for(i = 0; i < nr_of_labels; i++) {
    printf("same_labeling_lists, label: %c\n", labels[i]);
    dump_label_tmp_list(*(same_labeling_lists + i));
  }
  //exit(0);
#endif
  
  for(i = 0; i < nr_of_labels; i++) {
    same_labeling_list_indices[i] = 0;
    while(*(*(same_labeling_lists + i) + same_labeling_list_indices[i]) != TOT_END) {
      first = YES;
      while(*(*(same_labeling_lists + i) + same_labeling_list_indices[i]) != END) {
	/* update sorted_v_list */
	*(sorted_v_list + v_list_index) = *(*(same_labeling_lists + i) + same_labeling_list_indices[i]);
	v_list_index++;

	/* update pointers and label paths */
	if(first == YES) {
	  tmp_labeling = (cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling;
	  (cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling =
	    (char*)malloc_or_die((c+1) * sizeof(char));
	  memcpy((cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling,
		 tmp_labeling, (c) * sizeof(char));
	  ((cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling)[c] = labels[i];
#ifdef DEBUG_LABELING_UPDATE	 
	  printf("added label;labels[%d] = %c\n", i, labels[i]);
#endif
	  first = NO;
	  tmp_labeling = (cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling;
	}
	else {
	  (cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling = tmp_labeling;
	}
#ifdef DEBUG_LABELING_UPDATE
      printf("label length c = %d\n", c);
      dump_labeling((cur_rowp + *(*(same_labeling_lists + i) + same_labeling_list_indices[i]))->labeling, c);
#endif


	same_labeling_list_indices[i] += 1;
      }
      if(first == NO) {
	*(sorted_v_list + v_list_index) = V_LIST_NEXT;
	v_list_index++;
      }
      same_labeling_list_indices[i] += 1;
    }
  }
  
  *(sorted_v_list + v_list_index) = V_LIST_END;
  
  for(v = 1; v < nr_v; v++) {
    (cur_rowp + v)->is_updated = NO;
  }
  

  /* garbage collection */
  for(i = 0; i < nr_of_labels; i++) {
    free(same_labeling_lists[i]);
  }
  free(same_labeling_list_indices);
  free(same_labeling_lists);
}

void deallocate_row_labelings(struct one_best_s *prev_rowp, int nr_v)
{
  int dealloc_index;
  int cur_address;
  int v,w;
  int *dealloc_list;
  
#ifdef DEBUG_DEALLOCATE_LABELINGS
  printf("starting dealloc\n");
  printf("nr_v = %d\n", nr_v);
#endif

  dealloc_list = (int*)(malloc_or_die((nr_v+1) * sizeof(int)));

  for(v = 0; v < nr_v; v++) {
    (prev_rowp + v)->is_updated = NO;
  }
  
  
  /* deallocate last row's labelings */
  dealloc_index = 0;
  for(v = 0; v < nr_v; v++) {
    /* find all states with same labeling as this state up to previous row */
    if((prev_rowp + v)->is_updated == NO && (prev_rowp + v)->labeling != NULL) {
      cur_address = (int)((prev_rowp+v)->labeling);
      dealloc_list[dealloc_index] = v;
      dealloc_index++;
      (prev_rowp+v)->is_updated = YES;
      for(w = v+1; w < nr_v; w++) {
	if((int)((prev_rowp+w)->labeling) == cur_address) {
#ifdef DEBUG_DEALLOCATE_LABELINGS
	  printf("found same address, vertices %d and %d: %x\n", v, w, (prev_rowp + w)->labeling);
#endif	
	  (prev_rowp+w)->is_updated = YES;
	  (prev_rowp+w)->labeling = NULL;
	}
      }
    }
  }
  dealloc_list[dealloc_index] = END;
  
  for(dealloc_index = 0; dealloc_list[dealloc_index] != END; dealloc_index++) {
#ifdef DEBUG_DEALLOCATE_LABELINGS    
    printf("dealloc_index = %d\n", dealloc_index);
    printf("freeing labeling of vertex %d\n", dealloc_list[dealloc_index]);
#endif
    free((prev_rowp + dealloc_list[dealloc_index])->labeling);
#ifdef DEBUG_DEALLOCATE_LABELINGS    
    printf("done\n");
#endif
  } 
  
  free(dealloc_list);
#ifdef DEBUG_DEALLOCATE_LABELINGS  
  printf("exiting dealloc\n");
#endif
}