training_algorithms_multialpha.c

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

	      a_index_3 = get_alphabet_index(&seq_3[p], hmmp->alphabet_3, hmmp->a_size_3);
	    }
	    if(hmmp->nr_alphabets > 3) {
	      a_index_4 = get_alphabet_index(&seq_4[p], hmmp->alphabet_4, hmmp->a_size_4);
	    }
	    /* get result and add to matrix */
	    *(E_ulab + 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);
	    if(hmmp->nr_alphabets > 1) {
	      *(E_ulab_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);
	    }
	    if(hmmp->nr_alphabets > 2) {
	      *(E_ulab_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);
	    }
	    if(hmmp->nr_alphabets > 3) {
	      *(E_ulab_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);
	    }
	  }
	}
      }

      t_res_ulab = (forw_mtx + get_mtx_index(seq_len+1, hmmp->nr_v-1, hmmp->nr_v))->prob;
      /* some garbage collection */
      free(forw_mtx);
      free(backw_mtx);
      free(forw_scale);

     
      /********* calculations using labels *************/
      
      /* 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, YES, 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, YES, 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_lab += 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*/
	    a_index = get_alphabet_index(&seq[p], hmmp->alphabet, hmmp->a_size);
	    if(hmmp->nr_alphabets > 1) {
	      a_index_2 = get_alphabet_index(&seq_2[p], hmmp->alphabet_2, hmmp->a_size_2);
	    }
	    if(hmmp->nr_alphabets > 2) {
	      a_index_3 = get_alphabet_index(&seq_3[p], hmmp->alphabet_3, hmmp->a_size_3);
	    }
	    if(hmmp->nr_alphabets > 3) {
	      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_lab, YES,
				       multi_scoring_method);
	  }
	  /* 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++) {
	    a_index = get_alphabet_index(&seq[p], hmmp->alphabet, hmmp->a_size);
	    if(hmmp->nr_alphabets > 1) {
	      a_index_2 = get_alphabet_index(&seq_2[p], hmmp->alphabet_2, hmmp->a_size_2);
	    }
	    if(hmmp->nr_alphabets > 2) {
	      a_index_3 = get_alphabet_index(&seq_3[p], hmmp->alphabet_3, hmmp->a_size_3);
	    }
	    if(hmmp->nr_alphabets > 3) {
	      a_index_4 = get_alphabet_index(&seq_4[p], hmmp->alphabet_4, hmmp->a_size_4);
	    } 
	    /* get result and add to matrix */
	    *(E_lab + 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);
	    if(hmmp->nr_alphabets > 1) {
	      *(E_lab_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);
	    }
	    if(hmmp->nr_alphabets > 2) {
	      *(E_lab_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);
	    }
	    if(hmmp->nr_alphabets > 3) {
	      *(E_lab_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);
	    }
	  }
	}
      }

      /* 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 diff rd %d: %f\n", iteration, new_log_likelihood_ulab - new_log_likelihood_lab);
    }
    
#ifdef DEBUG_BW
    dump_T_matrix(hmmp->nr_v, hmmp->nr_v, T);
    dump_E_matrix(hmmp->nr_v, hmmp->a_size, E);
#endif
    
    /* 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_lab, 1);
    if(hmmp->nr_alphabets > 1) {
      recalculate_emiss_expectations_multi(hmmp, E_lab_2, 2);
    }
    if(hmmp->nr_alphabets > 2) {
      recalculate_emiss_expectations_multi(hmmp, E_lab_3, 3);
    }
    if(hmmp->nr_alphabets > 3) {
      recalculate_emiss_expectations_multi(hmmp, E_lab_4, 4);
    }

    recalculate_emiss_expectations_multi(hmmp, E_ulab, 1);
    if(hmmp->nr_alphabets > 1) {
      recalculate_emiss_expectations_multi(hmmp, E_ulab_2, 2);
    }
    if(hmmp->nr_alphabets > 2) {
      recalculate_emiss_expectations_multi(hmmp, E_ulab_3, 3);
    }
    if(hmmp->nr_alphabets > 3) {
      recalculate_emiss_expectations_multi(hmmp, E_ulab_4, 4);
    }
    
    /* recalculate transition expectations for tied transitions according
     * to the same scheme as for emission distribution groups */
    recalculate_trans_expectations_multi(hmmp, T_lab);
    recalculate_trans_expectations_multi(hmmp, T_ulab);
    
    
    /* update real T end E matrices */
    calculate_TE_contributions_multi(T, E, E_2, E_3, E_4, T_lab, E_lab, E_lab_2, E_lab_3, E_lab_4, T_ulab, E_ulab,
				     E_ulab_2, E_ulab_3, E_ulab_4, hmmp->emissions, hmmp->emissions_2, hmmp->emissions_3,
				     hmmp->emissions_4, hmmp->transitions, hmmp->nr_v, hmmp->a_size,
				     hmmp->a_size_2, hmmp->a_size_3, hmmp->a_size_4, hmmp->vertex_emiss_prior_scalers,
				     hmmp->vertex_emiss_prior_scalers_2, hmmp->vertex_emiss_prior_scalers_3,
				     hmmp->vertex_emiss_prior_scalers_4, iteration, hmmp->nr_alphabets);
    
    /* check if likelihood change is small enough, then we are done */
    if(fabs((new_log_likelihood_ulab - new_log_likelihood_lab) - (old_log_likelihood_ulab - old_log_likelihood_lab))
       < CML_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;
    }

    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);
      }
      
      
#ifdef DEBUG_PRIORS
      printf("Starting emission matrix update\n");
#endif
      
      /* update emission matrix using Dirichlet prior files if they exist*/
      priorp = *(hmmp->ed_ps + k);
      if(priorp != NULL && use_prior == YES) {
#ifdef DEBUG_PRIORS	
	printf("k = %d\n", k);
	printf("value = %x\n", priorp);
#endif
	update_emiss_mtx_prior_multi(hmmp, E, k, priorp, 1);
      }
      else if(use_emission_pseudo_counts == YES) /* update emissions matrix "normally" when dirichlet file is missing */ {
	update_emiss_mtx_pseudocount_multi(hmmp, E, k, 1);
      }
      else {
	update_emiss_mtx_std_multi(hmmp, E, k, 1);
      }
      

      if(hmmp->nr_alphabets > 1) {
	priorp = *(hmmp->ed_ps_2 + k);
	if(priorp != NULL && use_prior == YES) {
	  update_emiss_mtx_prior_multi(hmmp, E_2, k, priorp, 2);
	}
	else if(use_emission_pseudo_counts == YES) /* update emissions matrix "normally" when dirichlet file is missing */ {
	  update_emiss_mtx_pseudocount_multi(hmmp, E_2, k, 2);
	}
	else {
	  update_emiss_mtx_std_multi(hmmp, E_2, k, 2);
	}
      }

      if(hmmp->nr_alphabets > 2) {
	priorp = *(hmmp->ed_ps_3 + k);
	if(priorp != NULL && use_prior == YES) {
	  update_emiss_mtx_prior_multi(hmmp, E_3, k, priorp, 3);
	}
	else if(use_emission_pseudo_counts == YES) /* update emissions matrix "normally" when dirichlet file is missing */ {
	  update_emiss_mtx_pseudocount_multi(hmmp, E_3, k, 3);
	}
	else {
	  update_emiss_mtx_std_multi(hmmp, E_3, k, 3);
	}
      }

      if(hmmp->nr_alphabets > 3) {
	priorp = *(hmmp->ed_ps_4 + k);
	if(priorp != NULL && use_prior == YES) {
	  update_emiss_mtx_prior_multi(hmmp, E_4, k, priorp, 4);
	}
	else if(use_emission_pseudo_counts == YES) /* update emissions matrix "normally" when dirichlet file is missing */ {
	  update_emiss_mtx_pseudocount_multi(hmmp, E_4, k, 4);
	}
	else {
	  update_emiss_mtx_std_multi(hmmp, E_4, k, 4);
	}
      }
    }
    
#ifdef DEBUG_BW
    dump_trans_matrix(hmmp->nr_v, hmmp->nr_v, hmmp->transitions);
    dump_emiss_matrix(hmmp->nr_v, hmmp->a_size, hmmp->emissions);
#endif    
    
    /* some garbage collection */
    free(E);
    if(hmmp->nr_alphabets > 1) {
      free(E_2);
    }
    if(hmmp->nr_alphabets > 2) {
      free(E_3);
    }
    if(hmmp->nr_alphabets > 3) {
      free(E_4);
    }
    free(T);
    free(T_lab);
    free(E_lab);
    if(hmmp->nr_alphabets > 1) {
      free(E_lab_2);
    }
    if(hmmp->nr_alphabets > 2) {
      free(E_lab_3);
    }
    if(hmmp->nr_alphabets > 3) {
      free(E_lab_4);
    }
    free(T_ulab);
    free(E_ulab);
    if(hmmp->nr_alphabets > 1) {
      free(E_ulab_2);
    }
    if(hmmp->nr_alphabets > 2) {
      free(E_ulab_3);
    }
    if(hmmp->nr_alphabets > 3) {
      free(E_ulab_4);
    }
    max_nr_iterations--;
    iteration++;
  }
  while(max_nr_iterations > 0); /* break condition is also when log_likelihood_difference is
				 * smaller than THRESHOLD, checked inside the loop for
				 * better efficiency */
#ifdef DEBUG_BW2
  printf("exiting\n");
#endif
#ifdef DEBUG_BW
  dump_trans_matrix(hmmp->nr_v, hmmp->nr_v, hmmp->transitions);
  dump_emiss_matrix(hmmp->nr_v, hmmp->a_size, hmmp->emissions);
#endif

}



/* implementation of the baum-welch training algorithm using dirichlet prior mixture to
 * calculate update of emission (and transition) matrices and using a multiple sequence
 * alignment as the training sequence */
void extended_msa_baum_welch_dirichlet_multi(struct hmm_multi_s *hmmp, struct msa_sequences_multi_s *msa_seq_infop,
					     int nr_seqs, int annealing,
					     int use_gap_shares, int use_lead_columns, int use_labels, int use_transition_pseudo_counts,
					     int use_emission_pseudo_counts, int normalize, int scoring_method, int use_nr_occ,
					     int multi_scoring_method, double *aa_freqs,
					     double *aa_freqs_2, double *aa_freqs_3, double *aa_freqs_4, int use_prior)
{
  struct msa_sequences_multi_s *msa_seq_infop_start;
  double *T, *E, *E_2, *E_3, *E_4; /* matrices for the estimated number of times
		 * each transition (T) and emission (E) is used */
  double *T_lab, *E_lab, *E_lab_2, *E_lab_3, *E_lab_4, *T_ulab, *E_ulab, *E_ulab_2, *E_ulab_3, *E_ulab_4;
  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,i; /* loop counters, s loops over the sequences, p over the
		    * positions in the sequence, k and l over states, a over the alphabet,
		    * d over the distribution groups and i is a slush variable  */
  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; /* pointer to current sequence */
  double old_log_likelihood_lab, new_log_likelihood_lab;
  double old_log_likelihood_ulab, new_log_likelihood_ulab; /* to calculate when to stop */
  double likelihood; /* temporary variable for calculating likelihood of a sequence */
  int max_nr_iterations, iteration;
  double Eka_base;

  /* dirichlet prior variables */
  struct emission_dirichlet_s *priorp;