7 #include "likelihood.h"
8 #include <values.h> /* Has size of max_double, min_double. */
10 #include "debug_printout_flag.h"
12 /******************************************************************/
13 /*** mean[table] = 1/number_res * \Sum_res score_table(res) ****/
14 /*** covar[table][table2] = [1/number_res *
15 \Sum_res score_table(res)*score_table2(res)]
16 - mean_table * mean_table2 ****/
17 /** which equals.... ****/
18 /*** covar_def[table][table2] = 1/number_res *
19 \Sum_res [score_table(res)-mean_table]*
20 [score_table2(res)-mean_table2] ****/
21 /******************************************************************/
26 void multivariate_gauss_from_scfile2(char *sc_filename, int number_score_dim,
27 double mean[NUM_DIM_IN_ORIG_MATRIX],
28 double covar[NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX])
33 long n_below_mean[NUM_DIM_IN_ORIG_MATRIX];
34 long _counts[NUM_DIM_IN_ORIG_MATRIX][-4*MINSCORE +1];
35 long *counts[NUM_DIM_IN_ORIG_MATRIX];
37 int all_places_good, place[NUM_DIM_IN_ORIG_MATRIX];
38 int number_high=0, number_low=0;
39 int maxplaces[NUM_DIM_IN_ORIG_MATRIX], minplaces[NUM_DIM_IN_ORIG_MATRIX];
41 double score[NUM_DIM_IN_ORIG_MATRIX];
48 /********** Initialization. ******************************/
49 for (dim1=0; dim1<number_score_dim; dim1++) {
50 counts[dim1] = &(_counts[dim1][-2*MINSCORE]);
51 for (place[dim1]= 2*MINSCORE; place[dim1]<=-2*MINSCORE;
53 counts[dim1][place[dim1]]= 0;
54 n_below_mean[dim1] =0;
55 maxplaces[dim1] = 2*MINSCORE;
56 minplaces[dim1] = -2*MINSCORE;
59 for (dim2=0; dim2<number_score_dim; dim2++)
64 /***** Get a tuple of scores and the count on that tuple. ***/
66 sc_file = sopen(sc_filename,"r");
67 fprintf(stderr,"\nOpened file %s",sc_filename);
68 while (fgets(buf,500,sc_file)) {
70 /** fprintf(stderr,"\n%d",line_number); **/
73 for (dim1=0; dim1<number_score_dim; dim1++) {
74 sscanf(current_buf,"%lf", &score[dim1]);
75 current_buf= strchr(current_buf,' '); /* Get pointer to first space */
76 /* in current buf, and skip til */
77 while(current_buf[0] == ' ') /* not a space (so that next */
78 current_buf++; /* time 1st space comes after number */
81 sscanf(current_buf,"%d\n",&repetitions);
85 fprintf (stderr,"\nRead scores: ");
86 for (table=0; table < number_score_dim; table++)
87 fprintf (stderr," %lf",score[table]);
88 fprintf(stderr," %d ",repetitions);
91 all_places_good=1; /** If one of scores is out of range, don't include */
92 /** any of the scores. */
94 for (dim1=0; dim1<number_score_dim; dim1++) {
95 place[dim1] = ( (int)((score[dim1]-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
96 if (place[dim1] > -2*MINSCORE) {
97 fprintf(stderr, "Out of bounds: %d.\n", place[dim1]);
98 fprintf(stderr, "Score: %lf\n",score[dim1]);
99 /* counts[dim1][-2*MINSCORE]++;
100 maxplaces[dim1]=-2*MINSCORE; */
101 number_high+= repetitions;
104 else if (place[dim1] < 2*MINSCORE) {
105 fprintf(stderr, "Out of bounds: %d.\n", place[dim1]);
106 fprintf(stderr, "Score: %lf\n",score[dim1]);
107 /* counts[dim1][2*MINSCORE]++;
108 minplaces[dim1]=2*MINSCORE; */
109 number_low+= repetitions;
115 if (all_places_good) { /** Score the point in means and covariances. */
117 for (dim1=0; dim1<number_score_dim; dim1++) {
118 counts[dim1][place[dim1]]++;
119 if (place[dim1]> maxplaces[dim1])
120 maxplaces[dim1]=place[dim1];
121 if (place[dim1]< minplaces[dim1])
122 minplaces[dim1]= place[dim1];
123 mean[dim1] +=(double)repetitions*score[dim1];
125 for (dim2=0; dim2<number_score_dim; dim2++)
127 (double)repetitions*score[dim1]*score[dim2];
130 } /* Ends count on table after read in each scores. */
132 } /* Ends while scores to read in. */
134 fprintf(stderr,"\n%d out of high scores > %d were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE));
135 fprintf(stderr,"\n%d low scores < %d were scored as %d\n",number_low,MINSCORE, MINSCORE);
141 /** Now divide through by number of points did the expectations over. **/
143 for (dim1=0; dim1<number_score_dim; dim1++)
146 for (dim1=0; dim1<number_score_dim; dim1++)
147 for (dim2=0; dim2<number_score_dim; dim2++)
148 covar[dim1][dim2] = covar[dim1][dim2]/n -
149 mean[dim1]*mean[dim2];
155 int old_dimension(int old_num_dim_table0, int old_num_dim_table1,
156 int new_num_dim_table0, int new_num_dim_table1,
157 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
161 if (new_dim_number< new_num_dim_table0)
162 if (old_num_dim_table0 == 1)
165 return((int)multi_lib[0][new_dim_number]);
167 else /** Second (trimer) table. **/
168 if (old_num_dim_table1 == 1)
169 return(old_num_dim_table0);
171 return(old_num_dim_table0 +
172 (int)multi_lib[1][new_dim_number-new_num_dim_table0]);
176 void compute_submatrices2(
177 double means[MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX],
178 double means_submatrix[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
179 double covars[MAX_CLASS_NUMBER]
180 [NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX],
181 double inv[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
182 double det[MAX_CLASS_NUMBER],
183 char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM],
184 int old_num_dim_table0, int old_num_dim_table1,
185 int new_num_dim_table0, int new_num_dim_table1
188 int class, dim1, dim2, new_dim1, new_dim2;
189 double covars_sub[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
191 for (class=0; class<NUMBER_CLASSES; class++)
192 for(dim1=0; dim1< new_num_dim_table0 + new_num_dim_table1; dim1++) {
193 means_submatrix[class][dim1] = means[class]
194 [old_dimension(old_num_dim_table0, old_num_dim_table1,
195 new_num_dim_table0, new_num_dim_table1,
198 for(dim2=0; dim2< new_num_dim_table0 + new_num_dim_table1; dim2++)
199 covars_sub[class][dim1][dim2] = covars[class]
200 [old_dimension(old_num_dim_table0, old_num_dim_table1,
201 new_num_dim_table0, new_num_dim_table1,
203 [old_dimension(old_num_dim_table0, old_num_dim_table1,
204 new_num_dim_table0, new_num_dim_table1,
211 #ifdef DEBUG_PRINTOUT
212 for (class=0; class<NUMBER_CLASSES; class++) {
213 fprintf(stderr,"\n\nMean Submatrix %d:\n",class);
214 print_matrix(means_submatrix[class],1,
215 new_num_dim_table0 + new_num_dim_table1,stderr);
216 fprintf(stderr,"\n\nCovar Submatrix %d:\n", class);
217 print_matrix(covars_sub[class],
218 new_num_dim_table0 + new_num_dim_table1,
219 new_num_dim_table0 + new_num_dim_table1);
223 /*********************/
225 for (class=0; class<NUMBER_CLASSES; class++)
226 matrix_inverse_and_det(covars_sub[class], inv[class], &det[class],
227 new_num_dim_table0 + new_num_dim_table1);
234 /*** Init_prob[class] is the apriori prob. that a random residue */
235 /*** is in that class (i.e. our guess that a random residue should */
236 /*** score dimeric, trimeric, non-coiled. ********************************/
237 /*** Note that in estimating this, we want tetrameric coils to score **/
238 /*** somewhere between 2 and 3, so should put half their prob. in dimers **/
239 /*** and half in trimers. Want total of initial probs to be 1.... **/
240 /** THE first index in means and covars is the class number that that **/
241 /** gaussian is fit to. ****/
243 void calc_gauss_based_prob_like2(
244 double means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
245 double inverse_covars[MAX_CLASS_NUMBER]
246 [MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
247 double det[MAX_CLASS_NUMBER], int number_classes,
248 double init_prob[MAX_CLASS_NUMBER],
249 int number_score_dim,
250 double like_in_class[MAX_CLASS_NUMBER],
251 double scores[MAX_NUM_SCORE_DIM])
256 double score_diff_from_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
257 /* The vector is stored as a 1 x N matrix. **/
258 double total_like=0, max_like=0;
260 double inv[MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
261 double normalized_value;
262 double normalized_matrix[MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
263 double all_greater_than_minus_inf=1;
267 /*******Initialization. ****/
269 for (class =0 ; class< number_classes; class++)
270 for (dim =0; dim< number_score_dim; dim++) {
271 score_diff_from_means[class][dim][0] = scores[dim] -means[class][dim];
272 if (score_diff_from_means[class][dim][0] <= -HUGE_VAL)
273 all_greater_than_minus_inf =0;
278 /* if (!all_greater_than_minus_inf) return; */
280 /**************************/
283 for (class=0; class<number_classes; class++) {
284 /*** If x has some -HUGE_VAL scores, than all gaussian values are 0 (since */
285 /* gauss is 0 at infinity in any direction. So deal with that. ***/
287 if (normalize(number_score_dim, 1, number_score_dim, normalized_matrix,
288 inverse_covars[class], score_diff_from_means[class]) ) {
290 normalized_value= normalized_matrix[0][0];
291 like_in_class[class] = exp(-.5 * normalized_value) /
293 (double)number_score_dim/2)*
296 else { /* Score vector has some -HUGE_VAL scores. **/
297 /* So as a hack for now just give it 0 like. **/
298 like_in_class[class] =0;
301 /* fprintf(stderr,"\nnormalized_value = %lf",normalized_value); */
303 total_like += init_prob[class] * like_in_class[class];
304 if (like_in_class[class] > max_like) max_like = like_in_class[class];
307 for (class=0; all_one && (class<number_classes); class++)
308 if (init_prob[class] !=1) all_one=0;
310 if (!all_one) /** all_one is a flag to just output tuple of gauss_scores**/
311 /** Not the likelihoods/probs. **/
312 for (class=0; class<number_classes; class++) {
313 if (total_like >0) /* Otherwise all values are 0. **/
314 like_in_class[class] = init_prob[class] *like_in_class[class]/
324 /**** sc_filename(i) is the output of score_tuples for the known sequences in **/
325 /**** the ith class. **/
326 void get_gauss_param_all_classes2(char sc_filenames[MAX_CLASS_NUMBER][MAXLINE],
327 int number_score_dim,
328 double means[MAX_CLASS_NUMBER][NUM_DIM_IN_ORIG_MATRIX],
329 double covars[MAX_CLASS_NUMBER]
330 [NUM_DIM_IN_ORIG_MATRIX][NUM_DIM_IN_ORIG_MATRIX],
331 char gauss_param[MAXLINE])
339 int read_params_from_file =0;
341 /************************************************************************/
342 /**********************To scan the param from a file. *******************/
343 if (gauss_param[0] != ',') /* Either read or output to this file */
344 if (param_file=sopen(gauss_param,"r")) {
345 while (fgets(buf,1000,param_file)) {
346 /** Read in mean and covar. matrices. **/
347 if (!strncmp(buf,"Means for class",15)) {
348 read_params_from_file =1;
349 sscanf(buf,"Means for class %d", &class);
351 fgets(buf,1000,param_file);
353 for(dim1=0; dim1<number_score_dim; dim1++) {
354 current_buf= strchr(current_buf,' ');
355 /* Get pointer to first space */
356 /* in current buf, and skip til */
357 while(current_buf[0] == ' ') /* not a space (so that next */
358 current_buf++; /* time 1st space comes after number */
360 sscanf(current_buf,"%lf", &means[class][dim1]);
363 else if (!strncmp(buf,"Covariances for class",21)) {
364 sscanf(buf,"Covariances for class %d", &class);
365 for(dim1=0; dim1<number_score_dim; dim1++) {
366 fgets(buf,1000,param_file);
368 for(dim2=0; dim2<number_score_dim; dim2++) {
369 current_buf= strchr(current_buf,' ');
370 /* Get pointer to first space */
371 /* in current buf, and skip til */
372 while(current_buf[0] == ' ') /* not a space (so next */
373 current_buf++; /* time 1st space comes after number */
375 sscanf(current_buf,"%lf", &covars[class][dim1][dim2]);
383 /*************************************************************************/
385 if (read_params_from_file) {
386 #ifdef DEBUG_PRINTOUT
387 for (class=0; class<NUMBER_CLASSES; class++) {
388 fprintf(stderr,"\nRead in gauss parameters from %s\n",gauss_param);
389 fprintf(stderr,"\nMeans for class %d",class);
390 print_matrix(means[class],1,number_score_dim,stderr);
391 fprintf(stderr,"\nCovariances for class %d",class);
392 print_matrix(covars[class],number_score_dim,number_score_dim,stderr);
395 return; } /** All done. **/
397 if (gauss_param[0] != ',') /* output to this file */
398 (param_file=sopen(gauss_param,"w"));
400 for (class =0; class<NUMBER_CLASSES; class++) {
401 fprintf(stderr,"\n Class =%d",class);
402 multivariate_gauss_from_scfile2(sc_filenames[class], number_score_dim,
403 means[class], covars[class]);
405 fprintf(stderr,"\nMeans for class %d",class);
406 print_matrix2(means[class],1,number_score_dim,stderr);
407 fprintf(stderr,"\nCovariances for class %d",class);
408 print_matrix2(covars[class],number_score_dim,number_score_dim,stderr);
410 if (gauss_param[0] != ',') {
411 fprintf(param_file,"\nMeans for class %d",class);
412 print_matrix2(means[class],1,number_score_dim,param_file);
413 fprintf(param_file,"\nCovariances for class %d",class);
414 print_matrix2(covars[class],number_score_dim,number_score_dim,param_file);
417 if (gauss_param[0] != ',') /* output to this file */
423 double calc_total_coil_likelihood(double like_in_class[MAX_CLASS_NUMBER],
424 int number_classes, int reg,
426 double *max_class_like,
427 double init_class_prob[MAX_CLASS_NUMBER])
430 double total_coil_like=0;
431 double coil_init_prob=0;
432 /*** Note that in fact, the sum of */
433 /*** the like_in_class[positive_class] will be */
434 /*** 1-like_in_class[number_classes-1] = 1- like_in_pdb-. ****/
435 /*** We can sum the probabilities since we assume the classes are disjoint. */
437 for (positive_class = 0; positive_class<number_classes-1; positive_class++)
438 coil_init_prob += init_class_prob[positive_class];
440 if ( (Offset_to_Use != -1) || (reg != 7) ) {
441 for (positive_class = 0; positive_class<number_classes-1; positive_class++)
442 total_coil_like += like_in_class[positive_class];
444 if (total_coil_like > *max_class_like)
445 *max_class_like = total_coil_like;
447 else /* Find max scoring total coil reg for reg=7, max offset */
448 total_coil_like = *max_class_like;
451 return(total_coil_like);
454 /************************************************************************/
455 /*** Converts raw scores to probabilities. Stores dimer likelihoods in */
456 /*** score_dim_numbers 0 and 3, trimer likes in score_dim_numbers 1 and 4 */
457 /*** and non-coiled probs in score dim 2 and 5. This way, if we need the */
458 /*** residue scores later (say to change which coil scores are computed) */
459 /*** they will not need to be recomputed, since will still be in */
460 /*** tablenum + NUMBER_CLASSES. */
462 void convert_raw_scores_to_gauss_prob_like2(
463 double all_raw_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],int seq_len,
464 int number_tables, double class_means[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM],
465 double inverse_covars[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM],
466 double determinants[MAX_CLASS_NUMBER],
467 int number_classes, double init_class_prob[MAX_CLASS_NUMBER],
468 int number_score_dim, int consider_only_positive_classes, double bound,
469 int Offset_to_Use, int complement_last_class)
472 double like_in_class[MAX_CLASS_NUMBER];
473 double max_class_like[MAX_CLASS_NUMBER]; /** Used for computing max_reg **/
474 /** for output to ftotal_like. **/
475 double one_tuple_of_scores[MAX_NUM_SCORE_DIM];
480 int all_one =1, class;
482 int none_huge_val =1; /* Signal for when complementing last class. */
484 extern FILE *ftotal_like; /* If this is on, then only compute for the */
485 /* max register to save computation time. */
487 if ( (ftotal_like) && (Offset_to_Use != -1)) /* Saves time if computing */
488 first_reg =7; /* just ftotal_like, since then */
489 else first_reg=0; /* just need reg 7, unless in */
490 /* max Offset_to_Use. */
492 for (class=0; all_one && (class<number_classes); class++)
493 if (init_class_prob[class] !=1) all_one=0;
494 /** all_one is a flag to just output tuple of gauss_scores**/
495 /** Not the likelihoods. **/
498 for (i=0; i< seq_len; i++) {
500 if (Offset_to_Use == -1) /** Want to initialize for finding max reg. **/
501 for (tablenum=0; tablenum<number_classes; tablenum++)
502 max_class_like[tablenum] = 0;
505 for (reg=first_reg; reg<=7; reg++) { /** If outputting to ftotal_like **/
506 /** in max_reg mode, need this. **/
507 if ( (Offset_to_Use != -1) || (reg != 7) ) {
509 for (dim =0; dim<number_score_dim; dim++) {
510 one_tuple_of_scores[dim] = all_raw_scores[dim][i][reg];
511 if (one_tuple_of_scores[dim] <= -HUGE_VAL)
515 calc_gauss_based_prob_like2(class_means, inverse_covars,
516 determinants,number_classes,init_class_prob,
517 number_score_dim, like_in_class,
518 one_tuple_of_scores);
521 if (Offset_to_Use == -1) /* Want reg 7 to have max of other reg */
522 for (tablenum =0; tablenum<number_classes; tablenum++)
523 if (like_in_class[tablenum] > max_class_like[tablenum])
524 max_class_like[tablenum] = like_in_class[tablenum];
526 else /** Compute reg 7 as max of previous registers for offset -1 **/
527 for (tablenum=0; tablenum < number_classes; tablenum++)
528 like_in_class[tablenum] = max_class_like[tablenum];
532 /** Temporary hack to set the class we are interested to the table. **/
533 /** So make class 1 =cctb, class 2 = trimer-, class 3 = pdb- in .paircoil*/
534 /** Lets ASSUME that class= number_classes-1 is always pdb- */
535 /** so for table>=number_classes-1 output sum of other classes prob. **/
536 /** (which is the positive prob (prob of being coiled)). **/
537 for (tablenum=0; tablenum < number_classes; tablenum++) {
538 all_raw_scores[tablenum][i][reg] = like_in_class[tablenum];
539 } /** Count on tablenum **/
541 /*** Make last class have 1-* prob. (since sum of probs is 1 ****/
542 /*** Not all_raw_scores now contains probabilities. ****/
543 if (complement_last_class)
545 all_raw_scores[number_classes-1][i][reg]=
546 1 - all_raw_scores[number_classes -1][i][reg];
548 /************************************************************************/
549 /* Now make it so save a copy of the reg scores in tablenum+NUMBER_CLASSES */
550 for (tablenum=0; tablenum<number_classes; tablenum++)
551 all_raw_scores[tablenum+number_classes][i][reg]=
552 all_raw_scores[tablenum][i][reg];
555 /************************************************************************/
556 /*** If only concerned with positive classes (so the prob it is that oligim */
557 /*** given it is coiled. ****/
558 if ( (!all_one) && (consider_only_positive_classes))
559 /** Ignore class pdb-= num_class -1*/
560 for (tablenum=0; tablenum<number_classes-1; tablenum++)
561 if (all_raw_scores[number_classes-1][i][reg] >bound)
563 /*** Divide through by the prob are coiled to assume it is coiled. */
564 all_raw_scores[tablenum][i][reg] /=
565 all_raw_scores[number_classes -1][i][reg];
567 /************************************************************************/
571 } /** Count on reg.**/
578 int compute_dimension2(int tablenum, int libnumber,
579 int num_dist[MAX_TABLE_NUMBER],
580 int combine_dist[MAX_TABLE_NUMBER])
585 for (prev_table =0; prev_table<tablenum; prev_table++)
586 dimension+= num_table_dimensions(num_dist[prev_table],
587 combine_dist[prev_table]);
589 return (dimension + libnumber);
593 int num_table_dimensions(int num_dist_table,
594 int combine_dist_table)
597 return ((1-combine_dist_table)*num_dist_table +
605 /*** Note: See above in convert_raw_scores_to_gauss_like for what **/
606 /*** everything looks like. This takes the gauss_like and converts it **/
607 /*** from the old setting of only_coiled_classes to the new setting. **/
608 type_class_score_convert(Sequence sequence,
609 double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
610 double bound, int new_only_coiled_classes)
613 int i, tablenum, reg;
615 /** Note that the scores in tables above number_classes-2 is **/
616 /** the total positive likelihood always, and doesn't need changing. **/
618 for (i=0; i<sequence.seqlen; i++)
619 for (reg=0; reg<= POSNUM; reg++)
620 for (tablenum=0; tablenum<NUMBER_CLASSES-1; tablenum++)
621 if (all_like[NUMBER_CLASSES-1][i][reg]>bound)
622 if (new_only_coiled_classes ==0) { /* Make like for that class*/
623 all_like[tablenum][i][reg]= all_like[tablenum][i][reg] *
624 all_like[NUMBER_CLASSES-1][i][reg];
626 else{ /* Make into likelihood of being in that class given a pos. */
627 all_like[tablenum][i][reg]=
628 all_like[tablenum][i][reg]/all_like[NUMBER_CLASSES-1][i][reg];
634 void init_totals(double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
636 int dimer_box, trimer_box;
638 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
639 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
640 total_gauss_like[dimer_box][trimer_box] =0;
644 /** Note that the following assumes all_scores[0] is dimer gauss value and */
645 /** all_score[1] is trimer gauss value and all_scores[2] is pdb- gauss value */
646 /** so need to make sure of that. **/
647 void add_to_total_likes(double all_scores[MAX_NUM_SCORE_DIM]
648 [MAXSEQLEN][POSNUM+1],
649 double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER],
650 int seqlen, FILE *ftotal_like)
653 int dimer_box, trimer_box;
654 double init_prob_dimer, init_prob_trimer;
656 int need_to_adjust=0;
658 for (i=0; i<seqlen; i++) {
660 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++) {
661 init_prob_dimer = MIN_DIMER_PROB + (double)dimer_box*DIMER_GRID_STEP;
662 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++) {
663 init_prob_trimer = MIN_TRIMER_PROB +
664 (double)trimer_box*TRIMER_GRID_STEP;
665 total_gauss_like[dimer_box][trimer_box] +=
666 getdlog(init_prob_dimer * all_scores[0][i][7] +
667 init_prob_trimer * all_scores[1][i][7] +
668 (1-init_prob_dimer- init_prob_trimer) * all_scores[2][i][7]);
670 if (total_gauss_like[dimer_box][trimer_box] < -(MAXDOUBLE/2))
672 if (total_gauss_like[dimer_box][trimer_box] > (MAXDOUBLE/2))
677 /******* Make sure we don't over flow the double size by adding constants **/
678 /**** to each sum. ****/
680 if (need_to_adjust ==1) {
681 (*number_adjust_cuz_too_small)++;
682 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
683 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
684 total_gauss_like[dimer_box][trimer_box] += MAXDOUBLE/4;
686 if (need_to_adjust ==-1) {
687 (*number_adjust_cuz_too_big)++;
688 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
689 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
690 total_gauss_like[dimer_box][trimer_box] -= MAXDOUBLE/4;
696 fprintf(stderr,"\nError: Overflowed size of double in computeing total like. ");
703 void output_total_like(FILE *ftotal_like,
704 double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
707 int dimer_box, trimer_box;
708 double init_prob_dimer, init_prob_trimer;
709 double max_like=-HUGE_VAL;
712 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++)
713 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++)
714 if (total_gauss_like[dimer_box][trimer_box] > max_like)
715 max_like = total_gauss_like[dimer_box][trimer_box];
717 fprintf(ftotal_like, "The maximum sum of log likelihoods was %lf\n\n",max_like);
719 /** Normalize since values are near 0 ***/
721 for (dimer_box=0; dimer_box<GRID_SIZE_DIMER; dimer_box++) {
722 init_prob_dimer = MIN_DIMER_PROB + (double)dimer_box*DIMER_GRID_STEP;
723 for (trimer_box=0; trimer_box<GRID_SIZE_TRIMER; trimer_box++) {
724 init_prob_trimer = MIN_TRIMER_PROB +
725 (double)trimer_box*TRIMER_GRID_STEP;
726 fprintf(ftotal_like, "%lf %.4lf %.4lf\n",
727 total_gauss_like[dimer_box][trimer_box] - max_like,
728 init_prob_dimer,init_prob_trimer);