--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "scio.h"
+#include "scconst.h"
+#include "likelihood.h"
+#include <values.h> /* Has size of max_double, min_double. */
+
+#include "debug_printout_flag.h"
+
+/******************************************************************/
+/*** mean[table] = 1/number_res * \Sum_res score_table(res) ****/
+/*** covar[table][table2] = [1/number_res *
+ \Sum_res score_table(res)*score_table2(res)]
+ - mean_table * mean_table2 ****/
+/** which equals.... ****/
+/*** covar_def[table][table2] = 1/number_res *
+ \Sum_res [score_table(res)-mean_table]*
+ [score_table2(res)-mean_table2] ****/
+/******************************************************************/
+
+
+
+
+void multivariate_gauss_from_scfile2(char *sc_filename, int number_score_dim,
+ double mean[NUM_DIM_IN_ORIG_MATRIX],
+ double covar[NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX])
+
+{
+ FILE *sc_file;
+ long n=0;
+ long n_below_mean[NUM_DIM_IN_ORIG_MATRIX];
+ long _counts[NUM_DIM_IN_ORIG_MATRIX][-4*MINSCORE +1];
+ long *counts[NUM_DIM_IN_ORIG_MATRIX];
+ int table,table2;
+ int all_places_good, place[NUM_DIM_IN_ORIG_MATRIX];
+ int number_high=0, number_low=0;
+ int maxplaces[NUM_DIM_IN_ORIG_MATRIX], minplaces[NUM_DIM_IN_ORIG_MATRIX];
+
+ double score[NUM_DIM_IN_ORIG_MATRIX];
+ int repetitions;
+ char buf[500];
+ char *current_buf;
+ int line_number=0;
+ int dim1, dim2;
+
+/********** Initialization. ******************************/
+ for (dim1=0; dim1<number_score_dim; dim1++) {
+ counts[dim1] = &(_counts[dim1][-2*MINSCORE]);
+ for (place[dim1]= 2*MINSCORE; place[dim1]<=-2*MINSCORE;
+ ++place[dim1])
+ counts[dim1][place[dim1]]= 0;
+ n_below_mean[dim1] =0;
+ maxplaces[dim1] = 2*MINSCORE;
+ minplaces[dim1] = -2*MINSCORE;
+ mean[dim1]=0;
+
+ for (dim2=0; dim2<number_score_dim; dim2++)
+ covar[dim1][dim2]=0;
+ }
+
+
+/***** Get a tuple of scores and the count on that tuple. ***/
+
+ sc_file = sopen(sc_filename,"r");
+ fprintf(stderr,"\nOpened file %s",sc_filename);
+ while (fgets(buf,500,sc_file)) {
+ line_number++;
+ /** fprintf(stderr,"\n%d",line_number); **/
+ current_buf = buf;
+
+ for (dim1=0; dim1<number_score_dim; dim1++) {
+ sscanf(current_buf,"%lf", &score[dim1]);
+ current_buf= strchr(current_buf,' '); /* Get pointer to first space */
+ /* in current buf, and skip til */
+ while(current_buf[0] == ' ') /* not a space (so that next */
+ current_buf++; /* time 1st space comes after number */
+
+ }
+ sscanf(current_buf,"%d\n",&repetitions);
+
+
+/****
+ fprintf (stderr,"\nRead scores: ");
+ for (table=0; table < number_score_dim; table++)
+ fprintf (stderr," %lf",score[table]);
+ fprintf(stderr," %d ",repetitions);
+****/
+
+ all_places_good=1; /** If one of scores is out of range, don't include */
+ /** any of the scores. */
+
+ for (dim1=0; dim1<number_score_dim; dim1++) {
+ place[dim1] = ( (int)((score[dim1]-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
+ if (place[dim1] > -2*MINSCORE) {
+ fprintf(stderr, "Out of bounds: %d.\n", place[dim1]);
+ fprintf(stderr, "Score: %lf\n",score[dim1]);
+ /* counts[dim1][-2*MINSCORE]++;
+ maxplaces[dim1]=-2*MINSCORE; */
+ number_high+= repetitions;
+ all_places_good =0;
+ }
+ else if (place[dim1] < 2*MINSCORE) {
+ fprintf(stderr, "Out of bounds: %d.\n", place[dim1]);
+ fprintf(stderr, "Score: %lf\n",score[dim1]);
+ /* counts[dim1][2*MINSCORE]++;
+ minplaces[dim1]=2*MINSCORE; */
+ number_low+= repetitions;
+ all_places_good =0;
+ }
+ }
+
+
+ if (all_places_good) { /** Score the point in means and covariances. */
+ n += repetitions;
+ for (dim1=0; dim1<number_score_dim; dim1++) {
+ counts[dim1][place[dim1]]++;
+ if (place[dim1]> maxplaces[dim1])
+ maxplaces[dim1]=place[dim1];
+ if (place[dim1]< minplaces[dim1])
+ minplaces[dim1]= place[dim1];
+ mean[dim1] +=(double)repetitions*score[dim1];
+
+ for (dim2=0; dim2<number_score_dim; dim2++)
+ covar[dim1][dim2] +=
+ (double)repetitions*score[dim1]*score[dim2];
+ }
+
+ } /* Ends count on table after read in each scores. */
+
+ } /* Ends while scores to read in. */
+
+ fprintf(stderr,"\n%d out of high scores > %d were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE));
+ fprintf(stderr,"\n%d low scores < %d were scored as %d\n",number_low,MINSCORE, MINSCORE);
+
+ sclose(sc_file);
+
+
+
+ /** Now divide through by number of points did the expectations over. **/
+
+ for (dim1=0; dim1<number_score_dim; dim1++)
+ mean[dim1] /= n;
+
+ for (dim1=0; dim1<number_score_dim; dim1++)
+ for (dim2=0; dim2<number_score_dim; dim2++)
+ covar[dim1][dim2] = covar[dim1][dim2]/n -
+ mean[dim1]*mean[dim2];
+}
+
+
+
+
+int old_dimension(int old_num_dim_table0, int old_num_dim_table1,
+ int new_num_dim_table0, int new_num_dim_table1,
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ int new_dim_number)
+{
+
+ if (new_dim_number< new_num_dim_table0)
+ if (old_num_dim_table0 == 1)
+ return(0);
+ else
+ return((int)multi_lib[0][new_dim_number]);
+
+ else /** Second (trimer) table. **/
+ if (old_num_dim_table1 == 1)
+ return(old_num_dim_table0);
+ else
+ return(old_num_dim_table0 +
+ (int)multi_lib[1][new_dim_number-new_num_dim_table0]);
+}
+
+
+void compute_submatrices2(
+ double means[MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX],
+ double means_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
+ double covars[MAX_CLASS_NUMBER]
+ [NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX],
+ double inv[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
+ double det[MAX_CLASS_NUMBER],
+ char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
+ int old_num_dim_table0, int old_num_dim_table1,
+ int new_num_dim_table0, int new_num_dim_table1
+ )
+{
+ int class, dim1, dim2, new_dim1, new_dim2;
+ double covars_sub[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+
+ for (class=0; class<NUMBER_CLASSES; class++)
+ for(dim1=0; dim1< new_num_dim_table0 + new_num_dim_table1; dim1++) {
+ means_submatrix[class][dim1] = means[class]
+ [old_dimension(old_num_dim_table0, old_num_dim_table1,
+ new_num_dim_table0, new_num_dim_table1,
+ multi_lib, dim1)];
+
+ for(dim2=0; dim2< new_num_dim_table0 + new_num_dim_table1; dim2++)
+ covars_sub[class][dim1][dim2] = covars[class]
+ [old_dimension(old_num_dim_table0, old_num_dim_table1,
+ new_num_dim_table0, new_num_dim_table1,
+ multi_lib,dim1)]
+ [old_dimension(old_num_dim_table0, old_num_dim_table1,
+ new_num_dim_table0, new_num_dim_table1,
+ multi_lib,dim2)];
+
+ }
+
+
+
+#ifdef DEBUG_PRINTOUT
+ for (class=0; class<NUMBER_CLASSES; class++) {
+ fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
+ print_matrix(means_submatrix[class],1,
+ new_num_dim_table0 + new_num_dim_table1,stderr);
+ fprintf(stderr,"\n\nCovar Submatrix %d:\n", class);
+ print_matrix(covars_sub[class],
+ new_num_dim_table0 + new_num_dim_table1,
+ new_num_dim_table0 + new_num_dim_table1);
+
+ }
+#endif
+/*********************/
+
+ for (class=0; class<NUMBER_CLASSES; class++)
+ matrix_inverse_and_det(covars_sub[class], inv[class], &det[class],
+ new_num_dim_table0 + new_num_dim_table1);
+
+
+}
+
+
+
+/*** Init_prob[class] is the apriori prob. that a random residue */
+/*** is in that class (i.e. our guess that a random residue should */
+/*** score dimeric, trimeric, non-coiled. ********************************/
+/*** Note that in estimating this, we want tetrameric coils to score **/
+/*** somewhere between 2 and 3, so should put half their prob. in dimers **/
+/*** and half in trimers. Want total of initial probs to be 1.... **/
+/** THE first index in means and covars is the class number that that **/
+/** gaussian is fit to. ****/
+
+void calc_gauss_based_prob_like2(
+ double means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
+ double inverse_covars[MAX_CLASS_NUMBER]
+ [MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
+ double det[MAX_CLASS_NUMBER], int number_classes,
+ double init_prob[MAX_CLASS_NUMBER],
+ int number_score_dim,
+ double like_in_class[MAX_CLASS_NUMBER],
+ double scores[MAX_NUM_SCORE_DIM])
+
+{
+ int dim, class;
+ int table;
+ double score_diff_from_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ /* The vector is stored as a 1 x N matrix. **/
+ double total_like=0, max_like=0;
+
+ double inv[MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ double normalized_value;
+ double normalized_matrix[MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
+ double all_greater_than_minus_inf=1;
+
+ int dist;
+ int all_one =1;
+/*******Initialization. ****/
+
+ for (class =0 ; class< number_classes; class++)
+ for (dim =0; dim< number_score_dim; dim++) {
+ score_diff_from_means[class][dim][0] = scores[dim] -means[class][dim];
+ if (score_diff_from_means[class][dim][0] <= -HUGE_VAL)
+ all_greater_than_minus_inf =0;
+ }
+
+
+
+/* if (!all_greater_than_minus_inf) return; */
+
+/**************************/
+
+
+ for (class=0; class<number_classes; class++) {
+/*** If x has some -HUGE_VAL scores, than all gaussian values are 0 (since */
+/* gauss is 0 at infinity in any direction. So deal with that. ***/
+
+ if (normalize(number_score_dim, 1, number_score_dim, normalized_matrix,
+ inverse_covars[class], score_diff_from_means[class]) ) {
+
+ normalized_value= normalized_matrix[0][0];
+ like_in_class[class] = exp(-.5 * normalized_value) /
+ ( pow(2*Pi,
+ (double)number_score_dim/2)*
+ sqrt(det[class]) );
+ }
+ else { /* Score vector has some -HUGE_VAL scores. **/
+ /* So as a hack for now just give it 0 like. **/
+ like_in_class[class] =0;
+ }
+
+ /* fprintf(stderr,"\nnormalized_value = %lf",normalized_value); */
+
+ total_like += init_prob[class] * like_in_class[class];
+ if (like_in_class[class] > max_like) max_like = like_in_class[class];
+ }
+
+ for (class=0; all_one && (class<number_classes); class++)
+ if (init_prob[class] !=1) all_one=0;
+
+ if (!all_one) /** all_one is a flag to just output tuple of gauss_scores**/
+ /** Not the likelihoods/probs. **/
+ for (class=0; class<number_classes; class++) {
+ if (total_like >0) /* Otherwise all values are 0. **/
+ like_in_class[class] = init_prob[class] *like_in_class[class]/
+ total_like;
+ }
+
+}
+
+
+
+
+
+/**** sc_filename(i) is the output of score_tuples for the known sequences in **/
+/**** the ith class. **/
+void get_gauss_param_all_classes2(char sc_filenames[MAX_CLASS_NUMBER][MAXLINE],
+ int number_score_dim,
+ double means[MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX],
+ double covars[MAX_CLASS_NUMBER]
+ [NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX],
+ char gauss_param[MAXLINE])
+
+{
+ int class;
+ FILE *param_file;
+ char buf[1000];
+ char *current_buf;
+ int dim1, dim2;
+ int read_params_from_file =0;
+
+/************************************************************************/
+/**********************To scan the param from a file. *******************/
+ if (gauss_param[0] != ',') /* Either read or output to this file */
+ if (param_file=sopen(gauss_param,"r")) {
+ while (fgets(buf,1000,param_file)) {
+ /** Read in mean and covar. matrices. **/
+ if (!strncmp(buf,"Means for class",15)) {
+ read_params_from_file =1;
+ sscanf(buf,"Means for class %d", &class);
+
+ fgets(buf,1000,param_file);
+ current_buf =buf;
+ for(dim1=0; dim1<number_score_dim; dim1++) {
+ current_buf= strchr(current_buf,' ');
+ /* Get pointer to first space */
+ /* in current buf, and skip til */
+ while(current_buf[0] == ' ') /* not a space (so that next */
+ current_buf++; /* time 1st space comes after number */
+
+ sscanf(current_buf,"%lf", &means[class][dim1]);
+ }
+ }
+ else if (!strncmp(buf,"Covariances for class",21)) {
+ sscanf(buf,"Covariances for class %d", &class);
+ for(dim1=0; dim1<number_score_dim; dim1++) {
+ fgets(buf,1000,param_file);
+ current_buf=buf;
+ for(dim2=0; dim2<number_score_dim; dim2++) {
+ current_buf= strchr(current_buf,' ');
+ /* Get pointer to first space */
+ /* in current buf, and skip til */
+ while(current_buf[0] == ' ') /* not a space (so next */
+ current_buf++; /* time 1st space comes after number */
+
+ sscanf(current_buf,"%lf", &covars[class][dim1][dim2]);
+ }
+ }
+ }
+ }
+ sclose(param_file);
+ }
+
+/*************************************************************************/
+
+ if (read_params_from_file) {
+#ifdef DEBUG_PRINTOUT
+ for (class=0; class<NUMBER_CLASSES; class++) {
+ fprintf(stderr,"\nRead in gauss parameters from %s\n",gauss_param);
+ fprintf(stderr,"\nMeans for class %d",class);
+ print_matrix(means[class],1,number_score_dim,stderr);
+ fprintf(stderr,"\nCovariances for class %d",class);
+ print_matrix(covars[class],number_score_dim,number_score_dim,stderr);
+ }
+#endif
+ return; } /** All done. **/
+
+ if (gauss_param[0] != ',') /* output to this file */
+ (param_file=sopen(gauss_param,"w"));
+
+ for (class =0; class<NUMBER_CLASSES; class++) {
+ fprintf(stderr,"\n Class =%d",class);
+ multivariate_gauss_from_scfile2(sc_filenames[class], number_score_dim,
+ means[class], covars[class]);
+
+ fprintf(stderr,"\nMeans for class %d",class);
+ print_matrix2(means[class],1,number_score_dim,stderr);
+ fprintf(stderr,"\nCovariances for class %d",class);
+ print_matrix2(covars[class],number_score_dim,number_score_dim,stderr);
+
+ if (gauss_param[0] != ',') {
+ fprintf(param_file,"\nMeans for class %d",class);
+ print_matrix2(means[class],1,number_score_dim,param_file);
+ fprintf(param_file,"\nCovariances for class %d",class);
+ print_matrix2(covars[class],number_score_dim,number_score_dim,param_file);
+ }
+ }
+ if (gauss_param[0] != ',') /* output to this file */
+ sclose(param_file);
+
+}
+
+
+double calc_total_coil_likelihood(double like_in_class[MAX_CLASS_NUMBER],
+ int number_classes, int reg,
+ int Offset_to_Use,
+ double *max_class_like,
+ double init_class_prob[MAX_CLASS_NUMBER])
+{
+ int positive_class;
+ double total_coil_like=0;
+ double coil_init_prob=0;
+/*** Note that in fact, the sum of */
+/*** the like_in_class[positive_class] will be */
+/*** 1-like_in_class[number_classes-1] = 1- like_in_pdb-. ****/
+/*** We can sum the probabilities since we assume the classes are disjoint. */
+
+ for (positive_class = 0; positive_class<number_classes-1; positive_class++)
+ coil_init_prob += init_class_prob[positive_class];
+
+ if ( (Offset_to_Use != -1) || (reg != 7) ) {
+ for (positive_class = 0; positive_class<number_classes-1; positive_class++)
+ total_coil_like += like_in_class[positive_class];
+
+ if (total_coil_like > *max_class_like)
+ *max_class_like = total_coil_like;
+ }
+ else /* Find max scoring total coil reg for reg=7, max offset */
+ total_coil_like = *max_class_like;
+
+
+ return(total_coil_like);
+}
+
+/************************************************************************/
+/*** Converts raw scores to probabilities. Stores dimer likelihoods in */
+/*** score_dim_numbers 0 and 3, trimer likes in score_dim_numbers 1 and 4 */
+/*** and non-coiled probs in score dim 2 and 5. This way, if we need the */
+/*** residue scores later (say to change which coil scores are computed) */
+/*** they will not need to be recomputed, since will still be in */
+/*** tablenum + NUMBER_CLASSES. */
+
+void convert_raw_scores_to_gauss_prob_like2(
+ double all_raw_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],int seq_len,
+ int number_tables, double class_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
+ double inverse_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
+ double determinants[MAX_CLASS_NUMBER],
+ int number_classes, double init_class_prob[MAX_CLASS_NUMBER],
+ int number_score_dim, int consider_only_positive_classes, double bound,
+ int Offset_to_Use, int complement_last_class)
+
+{
+ double like_in_class[MAX_CLASS_NUMBER];
+ double max_class_like[MAX_CLASS_NUMBER]; /** Used for computing max_reg **/
+ /** for output to ftotal_like. **/
+ double one_tuple_of_scores[MAX_NUM_SCORE_DIM];
+ int i;
+ int reg;
+ int positive_class;
+ int dim, tablenum;
+ int all_one =1, class;
+ int first_reg;
+ int none_huge_val =1; /* Signal for when complementing last class. */
+
+ extern FILE *ftotal_like; /* If this is on, then only compute for the */
+ /* max register to save computation time. */
+
+ if ( (ftotal_like) && (Offset_to_Use != -1)) /* Saves time if computing */
+ first_reg =7; /* just ftotal_like, since then */
+ else first_reg=0; /* just need reg 7, unless in */
+ /* max Offset_to_Use. */
+
+ for (class=0; all_one && (class<number_classes); class++)
+ if (init_class_prob[class] !=1) all_one=0;
+ /** all_one is a flag to just output tuple of gauss_scores**/
+ /** Not the likelihoods. **/
+
+
+ for (i=0; i< seq_len; i++) {
+
+ if (Offset_to_Use == -1) /** Want to initialize for finding max reg. **/
+ for (tablenum=0; tablenum<number_classes; tablenum++)
+ max_class_like[tablenum] = 0;
+
+
+ for (reg=first_reg; reg<=7; reg++) { /** If outputting to ftotal_like **/
+ /** in max_reg mode, need this. **/
+ if ( (Offset_to_Use != -1) || (reg != 7) ) {
+ none_huge_val=1;
+ for (dim =0; dim<number_score_dim; dim++) {
+ one_tuple_of_scores[dim] = all_raw_scores[dim][i][reg];
+ if (one_tuple_of_scores[dim] <= -HUGE_VAL)
+ none_huge_val=0;
+ }
+
+ calc_gauss_based_prob_like2(class_means, inverse_covars,
+ determinants,number_classes,init_class_prob,
+ number_score_dim, like_in_class,
+ one_tuple_of_scores);
+
+
+ if (Offset_to_Use == -1) /* Want reg 7 to have max of other reg */
+ for (tablenum =0; tablenum<number_classes; tablenum++)
+ if (like_in_class[tablenum] > max_class_like[tablenum])
+ max_class_like[tablenum] = like_in_class[tablenum];
+ }
+ else /** Compute reg 7 as max of previous registers for offset -1 **/
+ for (tablenum=0; tablenum < number_classes; tablenum++)
+ like_in_class[tablenum] = max_class_like[tablenum];
+
+
+
+ /** Temporary hack to set the class we are interested to the table. **/
+ /** So make class 1 =cctb, class 2 = trimer-, class 3 = pdb- in .paircoil*/
+ /** Lets ASSUME that class= number_classes-1 is always pdb- */
+ /** so for table>=number_classes-1 output sum of other classes prob. **/
+ /** (which is the positive prob (prob of being coiled)). **/
+ for (tablenum=0; tablenum < number_classes; tablenum++) {
+ all_raw_scores[tablenum][i][reg] = like_in_class[tablenum];
+ } /** Count on tablenum **/
+
+/*** Make last class have 1-* prob. (since sum of probs is 1 ****/
+/*** Not all_raw_scores now contains probabilities. ****/
+ if (complement_last_class)
+ if (none_huge_val)
+ all_raw_scores[number_classes-1][i][reg]=
+ 1 - all_raw_scores[number_classes -1][i][reg];
+
+/************************************************************************/
+/* Now make it so save a copy of the reg scores in tablenum+NUMBER_CLASSES */
+ for (tablenum=0; tablenum<number_classes; tablenum++)
+ all_raw_scores[tablenum+number_classes][i][reg]=
+ all_raw_scores[tablenum][i][reg];
+
+
+/************************************************************************/
+/*** If only concerned with positive classes (so the prob it is that oligim */
+/*** given it is coiled. ****/
+ if ( (!all_one) && (consider_only_positive_classes))
+ /** Ignore class pdb-= num_class -1*/
+ for (tablenum=0; tablenum<number_classes-1; tablenum++)
+ if (all_raw_scores[number_classes-1][i][reg] >bound)
+ /* No longer raw*/
+ /*** Divide through by the prob are coiled to assume it is coiled. */
+ all_raw_scores[tablenum][i][reg] /=
+ all_raw_scores[number_classes -1][i][reg];
+
+/************************************************************************/
+
+
+
+ } /** Count on reg.**/
+ } /** Count on i **/
+}
+
+
+
+
+int compute_dimension2(int tablenum, int libnumber,
+ int num_dist[MAX_TABLE_NUMBER],
+ int combine_dist[MAX_TABLE_NUMBER])
+{
+ int prev_table;
+ int dimension=0;
+
+ for (prev_table =0; prev_table<tablenum; prev_table++)
+ dimension+= num_table_dimensions(num_dist[prev_table],
+ combine_dist[prev_table]);
+
+ return (dimension + libnumber);
+}
+
+
+int num_table_dimensions(int num_dist_table,
+ int combine_dist_table)
+{
+
+ return ((1-combine_dist_table)*num_dist_table +
+ combine_dist_table);
+}
+
+
+
+
+
+/*** Note: See above in convert_raw_scores_to_gauss_like for what **/
+/*** everything looks like. This takes the gauss_like and converts it **/
+/*** from the old setting of only_coiled_classes to the new setting. **/
+type_class_score_convert(Sequence sequence,
+ double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
+ double bound, int new_only_coiled_classes)
+{
+
+ int i, tablenum, reg;
+
+ /** Note that the scores in tables above number_classes-2 is **/
+ /** the total positive likelihood always, and doesn't need changing. **/
+
+ for (i=0; i<sequence.seqlen; i++)
+ for (reg=0; reg<= POSNUM; reg++)
+ for (tablenum=0; tablenum<NUMBER_CLASSES-1; tablenum++)
+ if (all_like[NUMBER_CLASSES-1][i][reg]>bound)
+ if (new_only_coiled_classes ==0) { /* Make like for that class*/
+ all_like[tablenum][i][reg]= all_like[tablenum][i][reg] *
+ all_like[NUMBER_CLASSES-1][i][reg];
+ }
+ else{ /* Make into likelihood of being in that class given a pos. */
+ all_like[tablenum][i][reg]=
+ all_like[tablenum][i][reg]/all_like[NUMBER_CLASSES-1][i][reg];
+ }
+}
+
+
+
+void init_totals(double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
+{
+ int dimer_box, trimer_box;
+
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
+ total_gauss_like[dimer_box][trimer_box] =0;
+}
+
+
+/** Note that the following assumes all_scores[0] is dimer gauss value and */
+/** all_score[1] is trimer gauss value and all_scores[2] is pdb- gauss value */
+/** so need to make sure of that. **/
+void add_to_total_likes(double all_scores[MAX_NUM_SCORE_DIM]
+ [MAXSEQLEN][POSNUM+1],
+ double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER],
+ int seqlen, FILE *ftotal_like)
+{
+
+ int dimer_box, trimer_box;
+ double init_prob_dimer, init_prob_trimer;
+ int i;
+ int need_to_adjust=0;
+
+ for (i=0; i<seqlen; i++) {
+ need_to_adjust=0;
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++) {
+ init_prob_dimer = MIN_DIMER_PROB + (double)dimer_box*DIMER_GRID_STEP;
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++) {
+ init_prob_trimer = MIN_TRIMER_PROB +
+ (double)trimer_box*TRIMER_GRID_STEP;
+ total_gauss_like[dimer_box][trimer_box] +=
+ getdlog(init_prob_dimer * all_scores[0][i][7] +
+ init_prob_trimer * all_scores[1][i][7] +
+ (1-init_prob_dimer- init_prob_trimer) * all_scores[2][i][7]);
+
+ if (total_gauss_like[dimer_box][trimer_box] < -(MAXDOUBLE/2))
+ need_to_adjust=1;
+ if (total_gauss_like[dimer_box][trimer_box] > (MAXDOUBLE/2))
+ need_to_adjust=-1;
+ }
+ }
+
+/******* Make sure we don't over flow the double size by adding constants **/
+/**** to each sum. ****/
+/*************
+ if (need_to_adjust ==1) {
+ (*number_adjust_cuz_too_small)++;
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
+ total_gauss_like[dimer_box][trimer_box] += MAXDOUBLE/4;
+ }
+ if (need_to_adjust ==-1) {
+ (*number_adjust_cuz_too_big)++;
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
+ total_gauss_like[dimer_box][trimer_box] -= MAXDOUBLE/4;
+ }
+******************/
+ }
+
+ if (need_to_adjust)
+ fprintf(stderr,"\nError: Overflowed size of double in computeing total like. ");
+
+}
+
+
+
+
+void output_total_like(FILE *ftotal_like,
+ double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
+{
+
+ int dimer_box, trimer_box;
+ double init_prob_dimer, init_prob_trimer;
+ double max_like=-HUGE_VAL;
+
+
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
+ if (total_gauss_like[dimer_box][trimer_box] > max_like)
+ max_like = total_gauss_like[dimer_box][trimer_box];
+
+ fprintf(ftotal_like, "The maximum sum of log likelihoods was %lf\n\n",max_like);
+
+/** Normalize since values are near 0 ***/
+ if (max_like != 0)
+ for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++) {
+ init_prob_dimer = MIN_DIMER_PROB + (double)dimer_box*DIMER_GRID_STEP;
+ for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++) {
+ init_prob_trimer = MIN_TRIMER_PROB +
+ (double)trimer_box*TRIMER_GRID_STEP;
+ fprintf(ftotal_like, "%lf %.4lf %.4lf\n",
+ total_gauss_like[dimer_box][trimer_box] - max_like,
+ init_prob_dimer,init_prob_trimer);
+ }
+ }
+
+}
+
+
+
+
git://source.jalview.org / jpred.git / blobdiff