training_algorithms_multialpha.c

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

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


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

#define INNER_BW_THRESHOLD 0.1
#define OUTER_BW_THRESHOLD 0.1
#define CML_THRESHOLD 0.0002
#define TRUE 1
#define FALSE -1
#define STARTRANDOM 0

#define REST_LETTER_INDEX 0.5

/* for simulated annealing */
#define INIT_TEMP 1.0 
#define INIT_COOL 0.8
#define ANNEAL_THRESHOLD 0.1
#define DONE -10
#define ACTIVE 25

/* for transition matrix pseudo count */
#define TRANSITION_PSEUDO_VALUE 0.1
#define EMISSION_PSEUDO_VALUE 1.0

//#define DEBUG_BW
//#define DEBUG_BW_TRANS
//#define DEBUG_BW2
//#define DEBUG_EXTBW
//#define DEBUG_PRIORS
//#define DEBUG_Tkl

extern int verbose;

void update_emiss_mtx_std_multi(struct hmm_multi_s*, double*, int, int);
void update_emiss_mtx_std_continuous_multi(struct hmm_multi_s*, double*, int, int);
void update_emiss_mtx_pseudocount_multi(struct hmm_multi_s*, double*, int, int);
void update_emiss_mtx_prior_multi(struct hmm_multi_s*, double*, int, struct emission_dirichlet_s*, int);
void update_trans_mtx_std_multi(struct hmm_multi_s*, double*, int);
void update_trans_mtx_pseudocount_multi(struct hmm_multi_s*, double*, int);
void update_tot_trans_mtx_multi(struct hmm_multi_s*);
void recalculate_emiss_expectations_multi(struct hmm_multi_s*, double*, int);
void recalculate_trans_expectations_multi(struct hmm_multi_s *hmmp, double *T);
double add_Eka_contribution_multi(struct hmm_multi_s*, struct letter_s*, struct forward_s*,
				  struct backward_s*, int, int, int);
double add_Eka_contribution_continuous_multi(struct hmm_multi_s*, struct letter_s*, struct forward_s*,
					     struct backward_s*, int, int, int, double*, int);
double add_Eka_contribution_msa_multi(struct hmm_multi_s*, struct msa_sequences_multi_s*, struct forward_s*,
				      struct backward_s*, int, int, int, int);
void add_Tkl_contribution_multi(struct hmm_multi_s*, struct letter_s*, struct letter_s*, struct letter_s*,
				struct letter_s*, struct forward_s*,
				struct backward_s*, double*, int, int,
				struct path_element*, int, int, int, int, double*, int use_labels, int multi_scoring_method);
void add_Tkl_contribution_msa_multi(struct hmm_multi_s*, struct msa_sequences_multi_s*, struct forward_s*,
				    struct backward_s*, double*,
				    int, int, struct path_element*, double*, int use_gap_shares, int use_lead_columns, int i,
				    int use_labels, int scoring_method, int normalize, int multi_scoring_method,
				    double *aa_freqs, double *aa_freqs_2, double *aa_freqs_3, double *aa_freqs_4);
void random_walk_multi(struct hmm_multi_s*, double*, double*, char**, int, int, int);
int silent_state_multi(int, struct hmm_multi_s*);
void anneal_E_matrix_multi(double temperature, double *E, struct hmm_multi_s *hmmp, int alphabet);
void anneal_T_matrix_multi(double temperature, double *T, struct hmm_multi_s *hmmp);
void calculate_TE_contributions_multi(double *T, double *E, double *E_2, double *E_3, double *E_4,
				      double *T_lab, double *E_lab, double *E_lab_2, double *E_lab_3, double *E_lab_4,
				      double *T_ulab, double *E_ulab, double *E_ulab_2, double *E_ulab_3, double *E_ulab_4,
				      double *emissions, double *emissions_2, double *emissions_3, double *emissions_4,
				      double *transitions, int nr_v, int a_size, int a_size_2, int a_size_3, int a_size_4,
				      double *emiss_prior_scalers, double *emiss_prior_scalers_2, double *emiss_prior_scalers_3,
				      double *emiss_prior_scalers_4, int rd, int nr_alphabets);
void add_to_E_multi(double *E, double Eka_base, struct msa_letter_s *msa_seq, int p, int k, int a_size, int normalize,
		    double *subst_mtx, int alphabet, int scoring_method, int use_nr_occ, int alphabet_type, double *emissions);



/************** baum-welch training algorithm ************************/




/* implementation of the baum-welch training algorithm using dirichlet prior mixture to
 * calculate update of emission (and transition) matrices */
void baum_welch_dirichlet_multi(struct hmm_multi_s *hmmp, struct sequence_multi_s *seqsp, int nr_seqs, int annealing, int use_labels,
				int use_transition_pseudo_counts, int use_emission_pseudo_counts,
				int multi_scoring_method, int use_prior)
{
  double *T, *E, *E_2, *E_3, *E_4; /* matrices for the estimated number of times
		 * each transition (T) and emission (E) is used */
  struct forward_s *forw_mtx; /* forward matrix */
  struct backward_s *backw_mtx; /* backward matrix */
  double *forw_scale; /* scaling array */
  int s,p,k,l,a,d; /* loop counters, s loops over the sequences, p over the
			* positions in the sequence, k and l over states, a over the alphabet
			* and d over the distribution groups */
  struct path_element *lp;
  double t_res, t_res_1, t_res_2, t_res_3; /* for temporary results */
  double t_res_4, t_res_5, t_res_6; /* for temporary results */
  double e_res, e_res_1, e_res_2, e_res_3; /* for temporary results */

  int seq_len; /* length of the seqences */
  int a_index, a_index_2, a_index_3, a_index_4; /* holds current letters index in the alphabet */
  struct letter_s *seq, *seq_2, *seq_3, *seq_4; /* pointer to current sequence */
  double old_log_likelihood, new_log_likelihood; /* to calculate when to stop */
  double likelihood; /* temporary variable for calculating likelihood of a sequence */
  int max_nr_iterations, iteration;

  /* dirichlet prior variables */
  struct emission_dirichlet_s *priorp;
  struct emission_dirichlet_s *priorp_2;
  struct emission_dirichlet_s *priorp_3;
  struct emission_dirichlet_s *priorp_4;

  /* simulated annealing variables */
  double temperature;
  double cooling_factor;
  int annealing_status;
 

  /* some initialization */
  old_log_likelihood = 9999.0;
  new_log_likelihood = 9999.0;
  max_nr_iterations = 20;
  iteration = 1;
  if(annealing == YES) {
    temperature = INIT_TEMP;
    cooling_factor = INIT_COOL;
    annealing_status = ACTIVE;
  }
  else {
    annealing_status = DONE;
  }
  

  do {  
    T = (double*)(malloc_or_die(hmmp->nr_v * hmmp->nr_v * 
			       sizeof(double)));
    if(hmmp->alphabet_type == DISCRETE) {
      E = (double*)(malloc_or_die(hmmp->nr_v * hmmp->a_size * 
				  sizeof(double)));
    }
    else {
      E = (double*)(malloc_or_die(hmmp->nr_v * (hmmp->a_size + 1) * 
				  sizeof(double)));
    }
    if(hmmp->nr_alphabets > 1) {
      if(hmmp->alphabet_type_2 == DISCRETE) {
	E_2 = (double*)(malloc_or_die(hmmp->nr_v * hmmp->a_size_2 * 
				      sizeof(double)));
      }
      else {
	E_2 = (double*)(malloc_or_die(hmmp->nr_v * (hmmp->a_size_2 + 1) * 
				      sizeof(double)));
      }
    }
    if(hmmp->nr_alphabets > 2) {
      if(hmmp->alphabet_type_3 == DISCRETE) {
	E_3 = (double*)(malloc_or_die(hmmp->nr_v * hmmp->a_size_3 * 
				      sizeof(double)));
      }
      else {
	E_3 = (double*)(malloc_or_die(hmmp->nr_v * (hmmp->a_size_3 + 1) * 
				      sizeof(double)));
      }
    }
    if(hmmp->nr_alphabets > 3) {
      if(hmmp->alphabet_type_4 == DISCRETE) {
	E_4 = (double*)(malloc_or_die(hmmp->nr_v * hmmp->a_size_4 * 
				      sizeof(double)));
      }
      else {
	E_4 = (double*)(malloc_or_die(hmmp->nr_v * (hmmp->a_size_4 + 1) * 
				      sizeof(double)));
      }
    }

    old_log_likelihood = new_log_likelihood;
    new_log_likelihood = 0.0;
    for(s = 0; s < nr_seqs; s++) {
      /* Convert sequence to 1...L for easier indexing */
      seq_len = (seqsp + s)->length;
      seq = (struct letter_s*) (malloc_or_die((seq_len + 2) * sizeof(struct letter_s)));
      memcpy(seq+1, (seqsp + s)->seq_1, seq_len * sizeof(struct letter_s));
      if(hmmp->nr_alphabets > 1) {
	seq_2 = (struct letter_s*) (malloc_or_die((seq_len + 2) * sizeof(struct letter_s)));
	memcpy(seq_2+1, (seqsp + s)->seq_2, seq_len * sizeof(struct letter_s));
      }
      if(hmmp->nr_alphabets > 2) {
	seq_3 = (struct letter_s*) malloc_or_die(((seq_len + 2) * sizeof(struct letter_s)));
	memcpy(seq_3+1, (seqsp + s)->seq_3, seq_len * sizeof(struct letter_s));
      }
      if(hmmp->nr_alphabets > 3) {
	seq_4 = (struct letter_s*) malloc_or_die(((seq_len + 2) * sizeof(struct letter_s)));
	memcpy(seq_4+1, (seqsp + s)->seq_4, seq_len * sizeof(struct letter_s));
      }

      
      /* calculate forward and backward matrices */
      forward_multi(hmmp, (seqsp + s)->seq_1,(seqsp + s)->seq_2, (seqsp + s)->seq_3, (seqsp + s)->seq_4,
		    &forw_mtx, &forw_scale, use_labels, multi_scoring_method);
      backward_multi(hmmp, (seqsp + s)->seq_1, (seqsp + s)->seq_2, (seqsp + s)->seq_3, (seqsp + s)->seq_4,
		     &backw_mtx, forw_scale, use_labels, multi_scoring_method);
      /* memory for forw_mtx, scale_mtx and
       * backw_mtx is allocated in the functions */
      
      /* update new_log_likelihood */
      likelihood = log10((forw_mtx +
			  get_mtx_index(seq_len+1, hmmp->nr_v-1, hmmp->nr_v))->prob);
      for(k = 0; k <= seq_len; k++) {
	likelihood = likelihood + log10(*(forw_scale + k));
      }
#ifdef DEBUG_BW
      dump_scaling_array(k-1,forw_scale);
      printf("likelihood = %f\n", likelihood);
#endif      
      new_log_likelihood += likelihood;
      
      for(k = 0; k < hmmp->nr_v-1; k++) /* k = from vertex */ {
      lp = *(hmmp->to_trans_array + k);
	while(lp->vertex != END) /* l = to-vertex */ {
	  for(p = 1; p <= seq_len; p++) {
	    
	     /* get alphabet index for c (add replacement letter stuff here) */
	    if(hmmp->alphabet_type == DISCRETE) {
	      a_index = get_alphabet_index(&seq[p], hmmp->alphabet, hmmp->a_size);
	    }	    
	    if(hmmp->nr_alphabets > 1 && hmmp->alphabet_type_2 == DISCRETE) {
	      a_index_2 = get_alphabet_index(&seq_2[p], hmmp->alphabet_2, hmmp->a_size_2);
	    }
	    if(hmmp->nr_alphabets > 2 && hmmp->alphabet_type_3 == DISCRETE) {
	      a_index_3 = get_alphabet_index(&seq_3[p], hmmp->alphabet_3, hmmp->a_size_3);
	    }
	    if(hmmp->nr_alphabets > 3 && hmmp->alphabet_type_4 == DISCRETE) {
	      a_index_4 = get_alphabet_index(&seq_4[p], hmmp->alphabet_4, hmmp->a_size_4);
	    } 

	    /* add T[k][l] contribution for this sequence */
	    add_Tkl_contribution_multi(hmmp, seq+1, seq_2+1, seq_3+1, seq_4+1, forw_mtx, backw_mtx,
				       forw_scale, p, k, lp, a_index, a_index_2, a_index_3, a_index_4, T, use_labels,
				       multi_scoring_method);
	    
	    /* continuous? */
	    
	  }
	  /* move on to next path */
	  while(lp->next != NULL) {
	    lp++;
	  }
	  lp++;
	}
	/* calculate E[k][a] contribution from this sequence */
	if(silent_state_multi(k, hmmp) != 0) {
	  for(p = 1; p <= seq_len; p++) {
	    if(hmmp->alphabet_type == DISCRETE) {
	      a_index = get_alphabet_index(&seq[p], hmmp->alphabet, hmmp->a_size);
	      *(E + get_mtx_index(k, a_index, hmmp->a_size)) +=
		add_Eka_contribution_multi(hmmp, seq+1, forw_mtx, backw_mtx, p, k, multi_scoring_method);
	    }
	    else {
	      add_Eka_contribution_continuous_multi(hmmp, seq+1, forw_mtx, backw_mtx, p, k, multi_scoring_method, E, 1);
	    }
	    if(hmmp->nr_alphabets > 1 && hmmp->alphabet_type_2 == DISCRETE) {
	      a_index_2 = get_alphabet_index(&seq_2[p], hmmp->alphabet_2, hmmp->a_size_2);
	      *(E_2 + get_mtx_index(k, a_index_2, hmmp->a_size_2)) +=
		add_Eka_contribution_multi(hmmp, seq_2+1, forw_mtx, backw_mtx, p, k, multi_scoring_method);
	    }
	    else if(hmmp->nr_alphabets > 1) {
	      add_Eka_contribution_continuous_multi(hmmp, seq_2+1, forw_mtx, backw_mtx, p, k, multi_scoring_method, E_2, 2);
	    }
	    if(hmmp->nr_alphabets > 2 && hmmp->alphabet_type_3 == DISCRETE) {
	      a_index_3 = get_alphabet_index(&seq_3[p], hmmp->alphabet_3, hmmp->a_size_3);
	      *(E_3 + get_mtx_index(k, a_index_3, hmmp->a_size_3)) +=
		add_Eka_contribution_multi(hmmp, seq_3+1, forw_mtx, backw_mtx, p, k, multi_scoring_method);
	    }
	    else  if(hmmp->nr_alphabets > 2) {
	      add_Eka_contribution_continuous_multi(hmmp, seq_3+1, forw_mtx, backw_mtx, p, k, multi_scoring_method, E_3, 3);
	    }
	    if(hmmp->nr_alphabets > 3 && hmmp->alphabet_type_4 == DISCRETE) {
	      a_index_4 = get_alphabet_index(&seq_4[p], hmmp->alphabet_4, hmmp->a_size_4);
	      *(E_4 + get_mtx_index(k, a_index_4, hmmp->a_size_4)) +=
		add_Eka_contribution_multi(hmmp, seq_4+1, forw_mtx, backw_mtx, p, k, multi_scoring_method);
	    }
	    else if(hmmp->nr_alphabets > 3) {
	      add_Eka_contribution_continuous_multi(hmmp, seq_4+1, forw_mtx, backw_mtx, p, k, multi_scoring_method, E_4, 4);
	    }
	  }
	}
      }
      /* some garbage collection */
      free(seq);
      if(hmmp->nr_alphabets > 1) {
	free(seq_2);
      }
      if(hmmp->nr_alphabets > 2) {
	free(seq_3);
      }
      if(hmmp->nr_alphabets > 3) {
	free(seq_4);
      }
      free(forw_mtx);
      free(backw_mtx);
      free(forw_scale); 
    }
    if(verbose == YES) {
      printf("log likelihood rd %d: %f\n", iteration, new_log_likelihood);
    }
    
#ifdef DEBUG_BW2
    dump_T_matrix(hmmp->nr_v, hmmp->nr_v, T);
    dump_E_matrix(hmmp->nr_v, hmmp->a_size, E);
    //dump_E_matrix(hmmp->nr_v, hmmp->a_size_2 + 1, E_2);
#endif
    
    /* check if likelihood change is small enough, then we are done */
    if(fabs(new_log_likelihood - old_log_likelihood) < INNER_BW_THRESHOLD && annealing_status == DONE) {
	break;
    }
    
    /* if simulated annealing is used, scramble results in E and T matrices */
    if(annealing == YES && temperature > ANNEAL_THRESHOLD) {
      anneal_E_matrix_multi(temperature, E, hmmp, 1);
      if(hmmp->nr_alphabets > 1) {
	anneal_E_matrix_multi(temperature, E_2, hmmp, 2);
      }
      if(hmmp->nr_alphabets > 2) {
	anneal_E_matrix_multi(temperature, E_3, hmmp, 3);
      }
      if(hmmp->nr_alphabets > 3) {
	anneal_E_matrix_multi(temperature, E_4, hmmp, 4);
      }
      anneal_T_matrix_multi(temperature, T, hmmp);
      temperature = temperature * cooling_factor;
    }
    
    if(temperature < ANNEAL_THRESHOLD) {
      annealing_status = DONE;
    }
    
    /* recalculate emission expectations according to distribution groups 
     * by simply taking the mean of the expected emissions within this group
     * for each letter in the alphabet and replacing each expectation for the
     * letter with this value for every member of the distribution group */
    recalculate_emiss_expectations_multi(hmmp, E, 1);
    if(hmmp->nr_alphabets > 1) {
      recalculate_emiss_expectations_multi(hmmp, E_2, 2);
    }
    if(hmmp->nr_alphabets > 2) {
      recalculate_emiss_expectations_multi(hmmp, E_3, 3);
    }
    if(hmmp->nr_alphabets > 3) {
      recalculate_emiss_expectations_multi(hmmp, E_4, 4);
    }
    
    /* recalculate transition expectations for tied transitions according
     * to the same scheme as for emission distribution groups */
    recalculate_trans_expectations_multi(hmmp, T);
    
    for(k = 0; k < hmmp->nr_v-1; k++) /* k = from-vertex */ {
      /* update transition matrix */
      if(use_transition_pseudo_counts == YES) {
	update_trans_mtx_pseudocount_multi(hmmp, T, k);
      }
      else {
	update_trans_mtx_std_multi(hmmp, T, k);
      }