std_calculation_funcs.c

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

      }
    }
  }
  
  if(t_res_3 < 0.0 || t_res_3 > 1000000000000.0) {
    printf("Error: got strange picasso product state result value\n");
    exit(0);
  }
  
  return t_res_3;
}


double get_picasso_sym_statescore(int a_size, int use_gap_shares, int use_prior_shares, struct msa_letter_s *msa_seq,
				  int p, double *emissions,  int vertex, int normalize, double *gap_shares, double *aa_freqs)
{
  
  double seq_normalizer;
  double state_normalizer;
  double subst_mtx_normalizer;
  int a_index;
  double t_res_3;
  
  seq_normalizer = 0.0;
  state_normalizer = 0.0;
  subst_mtx_normalizer = 0.0;
  

  t_res_3 = 1.0;
  /* scoring using picasso-product method */

  for(a_index = 0; a_index < a_size; a_index++) {
    if(use_prior_shares == YES) {
      if((msa_seq + get_mtx_index(p,a_index, a_size+1))->prior_share != 0.0  &&
	 *(emissions + get_mtx_index(vertex, a_index, a_size)) != SILENT) {
	t_res_3 *= pow((msa_seq + get_mtx_index(p,a_index, a_size+1))->prior_share / *(aa_freqs + a_index),
		       *(emissions + (vertex * a_size + a_index))) *
	  pow(*(emissions + (vertex * a_size + a_index)) / *(aa_freqs + a_index),
	      (msa_seq + get_mtx_index(p,a_index, a_size+1))->prior_share);
      }
    }
    else if(use_gap_shares == YES) {
      if((msa_seq + get_mtx_index(p,a_size, a_size+1))->share == 1.0) {
	printf("Error: all gap column in sequence\n");
	exit(0);
      }
      if((msa_seq + get_mtx_index(p,a_index, a_size+1))->share != 0.0  &&
	 *(emissions + get_mtx_index(vertex, a_index, a_size)) != SILENT) {
	t_res_3 *= pow(((msa_seq + get_mtx_index(p,a_index, a_size+1))->share /
			(1.0 -(msa_seq + get_mtx_index(p,a_size, a_size+1))->share)) / *(aa_freqs + a_index),
		       *(emissions + get_mtx_index(vertex, a_index, a_size))) *
	  pow(*(emissions + (vertex * a_size + a_index)) / *(aa_freqs + a_index),
	      (msa_seq + get_mtx_index(p,a_index, a_size+1))->share /
	      (1.0 - (msa_seq + get_mtx_index(p,a_size, a_size+1))->share));
      }
    }
    else {
      //printf(" *(aa_freqs + a_index) = %f\n",  *(aa_freqs + a_index));
      //printf("(msa_seq + (p * (a_size+1) + a_index))->share = %f\n", (msa_seq + get_mtx_index(p,a_index, a_size+1))->share);
      //printf("*(emissions + get_mtx_index(vertex, a_index, a_size)) = %f\n",*(emissions + get_mtx_index(vertex, a_index, a_size))); 
      if((msa_seq + get_mtx_index(p,a_index, a_size+1))->share != 0.0 &&
	 *(emissions + get_mtx_index(vertex, a_index, a_size)) != SILENT) {
	t_res_3 *= pow((msa_seq + get_mtx_index(p,a_index, a_size+1))->share / *(aa_freqs + a_index),
		       *(emissions + get_mtx_index(vertex, a_index, a_size))) *
	  pow(*(emissions + (vertex * a_size + a_index)) / *(aa_freqs + a_index),
	      (msa_seq + get_mtx_index(p,a_index, a_size+1))->share);
      }
    }
  }

  if(t_res_3 < 0.0 || t_res_3 > 1000000000000.0) {
    printf("Error: got strange picasso product state result value\n");
    exit(0);
  }
  
  return t_res_3;
}


double get_subst_mtx_product_statescore(int a_size, int use_gap_shares, int use_prior_shares, struct msa_letter_s *msa_seq,
					int p, double *emissions, int vertex, double *subst_mtx)
{
  int a_index, a_index2;
  double t_res_3;
  
  t_res_3 = 0.0;
  for(a_index = 0; a_index < a_size; a_index++) {
    for(a_index2 = 0; a_index2 < a_size; a_index2++) {
      if(use_gap_shares == YES) {
	t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	  (msa_seq + get_mtx_index(p,a_index2, a_size+1))->share /
	  (1.0 -(msa_seq + get_mtx_index(p,a_size, a_size+1))->share) *
	  *(subst_mtx + (a_index * a_size + a_index2));
      }
      else {
	t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	  (msa_seq + get_mtx_index(p,a_index2, a_size+1))->share /
	  *(subst_mtx + (a_index * a_size + a_index2));
      }
    }
  }
  
  if(t_res_3 < 0.0) {
    printf("Error: got strange subst mtx product state result value\n");
    exit(0);
  }
  
  return t_res_3;
}


double get_subst_mtx_dot_product_statescore(int a_size, int use_gap_shares, int use_prior_shares, struct msa_letter_s *msa_seq,
					    int p, double *emissions,  int vertex, int normalize, double *gap_shares,
					    int query_index, double *subst_mtx)
{
  double seq_normalizer;
  double state_normalizer;
  double subst_mtx_normalizer;
  int a_index;
  double t_res_3;
  
  seq_normalizer = 0.0;
  state_normalizer = 0.0;
  subst_mtx_normalizer = 0.0;
  
  t_res_3 = 0.0;
  /* scoring using subst_mtx_dot-product method */
  for(a_index = 0; a_index < a_size; a_index++) {
    if(use_prior_shares == YES) {
      t_res_3 += *(emissions + (vertex * a_size + a_index)) *
	(msa_seq + (p * (a_size+1) + a_index))->prior_share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));
      if(normalize == YES) {
	seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->prior_share, 2);
	state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
      }
    }
    else if(use_gap_shares == YES) {
      if((msa_seq + get_mtx_index(p,a_size, a_size+1))->share == 1.0) {
	printf("Error: all gap column in sequence\n");
	  exit(0);
      }
      t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	(msa_seq + get_mtx_index(p,a_index, a_size+1))->share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size)) /
	(1.0 -(msa_seq + get_mtx_index(p,a_size, a_size+1))->share);
      if(normalize == YES) {
	seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->share /
			      (1.0 - *(gap_shares + p)), 2);
	state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
      }
    }
    else {
      t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	(msa_seq + get_mtx_index(p,a_index, a_size+1))->share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));
      if(normalize == YES) {
	seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->share, 2);
	state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
      }
    }
  }
  
  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_normalizer);
    subst_mtx_normalizer = sqrt(subst_mtx_normalizer);
#ifdef DEBUG
    printf("state_normalizer = %f\n", state_normalizer);
      printf("seq_normalizer = %f\n", seq_normalizer);
#endif      
      if(t_res_3 != 0.0) {
	t_res_3 = t_res_3 / (seq_normalizer * state_normalizer * subst_mtx_normalizer);
      }
  }

  if(t_res_3 < 0.0) {
    printf("Error: got strange subst mtx dot product state result value\n");
    exit(0);
  }
  
  return t_res_3;
}


double get_subst_mtx_dot_product_prior_statescore(int a_size, int use_gap_shares, int use_prior_shares, struct msa_letter_s *msa_seq,
						  int p, double *emissions,  int vertex, int normalize, double *gap_shares,
						  int query_index, double *subst_mtx)
{
  double seq_normalizer;
  double state_normalizer;
  double subst_mtx_normalizer;
  int a_index;
  double t_res_3;
  double rest_share, default_share;
  
  seq_normalizer = 0.0;
  state_normalizer = 0.0;
  subst_mtx_normalizer = 0.0;
  default_share = 1.0 / (double)(a_size);

  t_res_3 = 0.0;
  
  /* scoring using dot-product method */
  rest_share = 1.0;
  for(a_index = 0; a_index < a_size; a_index++) {
    if(use_prior_shares == YES) {
      t_res_3 += *(emissions + (vertex * a_size + a_index)) *
	(msa_seq + (p * (a_size+1) + a_index))->prior_share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));
      if(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)) != 0.0) {
	rest_share = rest_share - (msa_seq + (p * (a_size+1) + a_index))->prior_share;
	if(normalize == YES) {
	  seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->prior_share, 2);
	  state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	  subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
	}
      }
    }
    else if(use_gap_shares == YES) {
      if((msa_seq + get_mtx_index(p,a_size, a_size+1))->share == 1.0) {
	printf("Error: all gap column in sequence\n");
	exit(0);
      }
      t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	(msa_seq + get_mtx_index(p,a_index, a_size+1))->share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size)) /
	(1.0 -(msa_seq + get_mtx_index(p,a_size, a_size+1))->share);
      if(*(subst_mtx + get_mtx_index(a_index, a_index, a_size)) != 0.0) {
	rest_share = rest_share - (msa_seq + (p * (a_size+1) + a_index))->share / 
	  (1.0 - *(gap_shares + p));
	if(normalize == YES) {
	  seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->share /
				(1.0 - *(gap_shares + p)), 2);
	  state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	  subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
	}
      }
      
    }
    else {
      t_res_3 += *(emissions + get_mtx_index(vertex, a_index, a_size)) *
	(msa_seq + get_mtx_index(p,a_index, a_size+1))->share *
	*(subst_mtx + get_mtx_index(query_index, a_index, a_size));
      if(*(subst_mtx + get_mtx_index(a_index, a_index, a_size)) != 0.0) {
	rest_share = rest_share - (msa_seq + (p * (a_size+1) + a_index))->share;
	if(normalize == YES) {
	  seq_normalizer += pow((msa_seq + (p * (a_size+1) + a_index))->share, 2);
	  state_normalizer += pow(*(emissions + (vertex * a_size + a_index)), 2);
	  subst_mtx_normalizer +=  pow(*(subst_mtx + get_mtx_index(query_index, a_index, a_size)), 2);
	}
      }
    }
  }
  if(rest_share < 0.0) {
    rest_share = 0.0;
  }
  t_res_3 += default_share * rest_share;
  seq_normalizer += pow(rest_share, 2);
  state_normalizer += pow(default_share, 2);
  
  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_normalizer);
    subst_mtx_normalizer = sqrt(subst_mtx_normalizer);
#ifdef DEBUG_BW
    printf("state_normalizer = %f\n", state_normalizer);
    printf("seq_normalizer = %f\n", seq_normalizer);
#endif
    if(t_res_3 != 0.0) {
      t_res_3 = t_res_3 / (seq_normalizer * state_normalizer * subst_mtx_normalizer);
    }
  }

  if(t_res_3 < 0.0) {
    printf("Error: got strange subst mtx dot product prior state result value\n");
    exit(0);
  }
  
  return t_res_3;
}



/************************************* add to E methods *********************************************/
void add_to_E_continuous(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
			 int k, int a_size, double *emissions)
{
  double mean_value, varians;
  int i,j;
  double continuous_score_all, continuous_score_j, gamma_p_j;

  mean_value = (msa_seq + get_mtx_index(p,0,a_size+1))->share;
  
  
  continuous_score_all = 0.0;
  for(j = 0; j < a_size / 3; j++) {
    continuous_score_all += 
      get_single_gaussian_statescore(*(emissions + get_mtx_index(k, (j * 3), a_size)),
				     *(emissions + get_mtx_index(k, (j * 3 + 1), a_size)),
				     mean_value) *
      *((emissions) + (k * (a_size)) + (j * 3 + 2));
  }
  
  for(j = 0; j < a_size / 3; j++) {
    continuous_score_j = 
      get_single_gaussian_statescore(*(emissions + get_mtx_index(k, (j * 3), a_size)),
				     *(emissions + get_mtx_index(k, (j * 3 + 1), a_size)),
				     mean_value) *
      *((emissions) + (k * (a_size)) + (j * 3 + 2));
    varians = pow((msa_seq + get_mtx_index(p,0,a_size+1))->share - *(emissions + get_mtx_index(k, j * 3, a_size)), 2);
    if(continuous_score_all > 0.0) {
      gamma_p_j = Eka_base * continuous_score_j / continuous_score_all;
    }
    else {
      gamma_p_j = 0.0;
    }
    *(E + get_mtx_index(k, j * 3, a_size + 1)) += mean_value * gamma_p_j;
    *(E + get_mtx_index(k, j * 3 + 1, a_size + 1)) += varians * gamma_p_j;
    *(E + get_mtx_index(k, j * 3 + 2, a_size + 1)) += gamma_p_j;
  }
  
  *(E + get_mtx_index(k, j * 3, a_size + 1)) += Eka_base;
}



void add_to_E_dot_product(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
			  int k, int a_size, int normalize)
{
  /* k is a state index, p is a sequence position index */
  
  int a_index;
  double prf_column_length;

  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);
      
    }
  }
  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);
    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) / prf_column_length;
    }
  }
}

void add_to_E_dot_product_picasso(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p,
				  int k, int a_size, int normalize)
{
  /* k is a state index, p is a sequence position index */
  
  int a_index;
  double prf_column_length;