std_calculation_funcs.c

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

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>




#include "structs.h"
#include "funcs.h"

//#define DEBUG_LABELING_UPDATE
//#define DEBUG_DEALLOCATE_LABELINGS


#define REST_LETTER_INDEX 0.5

#define V_LIST_END  -99
#define V_LIST_NEXT  -9


double get_single_gaussian_statescore(double mu, double sigma_square, double letter)
{

  double res;

  if(sigma_square <= 0.0) {
    return 0.0;
  }
  else {
    res = exp(0.0 - (pow((letter-mu),2) / (2.0 * sigma_square))) / sqrt(sigma_square * 2.0 * 3.141592655);
    //printf("mu = %f, letter = %f, sigma_square = %f,  res = %f\n", mu, letter, sigma_square, res);
    return res;
  }
}

double get_dp_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 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 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;
      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);
      }
    }
    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 /
	(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);
      }
    }
    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;
      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);
      }
    }
  }
  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_normalizer);
    if(t_res_3 != 0.0) {
      t_res_3 = t_res_3 / (seq_normalizer * state_normalizer);
    }
  }

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


double get_dp_picasso_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 = 0.0;
  /* scoring using 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 /  *(aa_freqs + a_index);
      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);
      }
    }
    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 /
	((1.0 -(msa_seq + get_mtx_index(p,a_size, a_size+1))->share) *  *(aa_freqs + a_index));
      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);
      }
    }
    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 /  *(aa_freqs + a_index);
      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);
      }
    }
  }
  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_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);
    }
  }

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


double get_sjolander_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 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 sjolander score method */

  for(a_index = 0; a_index < a_size; a_index++) {
    if(use_prior_shares == YES) {
      t_res_3 *= pow(*(emissions + (vertex * a_size + a_index)), 
		    (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);
      }
     
    }
    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 *= pow(*(emissions + get_mtx_index(vertex, a_index, a_size)),
		     (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));
      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);
      }
    }
    else {
      t_res_3 *= pow(*(emissions + get_mtx_index(vertex, a_index, a_size)),
		     (msa_seq + get_mtx_index(p,a_index, a_size+1))->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);
      }
    }
  }

  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_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);
    }
  }
  if(t_res_3 < 0.0) {
    printf("Error: got strange geometric mean state result value\n");
    printf("t_res_3 = %f\n", t_res_3);
    exit(0);
  }
  return t_res_3;
}



double get_sjolander_reversed_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 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 sjolander score 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) {
      	t_res_3 *= pow((msa_seq + get_mtx_index(p,a_index, a_size+1))->prior_share,
      		       *(emissions + get_mtx_index(vertex, 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);
      }
    }
    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) {
	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),
		       *(emissions + get_mtx_index(vertex, a_index, a_size)));
      }
      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);
      }
    }
    else {
      if((msa_seq + get_mtx_index(p,a_index, a_size+1))->share != 0.0) {
	t_res_3 *= pow((msa_seq + get_mtx_index(p,a_index, a_size+1))->share,
		       *(emissions + get_mtx_index(vertex, 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);
      }
    }
  }
  if(normalize == YES) {
    seq_normalizer = sqrt(seq_normalizer);
    state_normalizer = sqrt(state_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);
    }
  }
  if(t_res_3 < 0.0) {
    printf("Error: got strange geometric mean state result value\n");
    printf("t_res_3 = %f\n", t_res_3);
    exit(0);
  }
  
  return t_res_3;
}


double get_picasso_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)));
      }
    }
    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)));
      }
    }
    else {
      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)));