sources/multicoil/multivariate_like2.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include "scio.h"
   6 #include "scconst.h"
   7 #include "likelihood.h"
   8 #include  <values.h>  /* Has size of max_double, min_double. */
   9
  10 #include "debug_printout_flag.h"
  11
  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 /******************************************************************/
  22
  23
  24
  25
  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])
  29
  30 {
  31   FILE *sc_file;
  32   long  n=0;
  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];
  36   int table,table2;
  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];
  40
  41   double score[NUM_DIM_IN_ORIG_MATRIX];
  42   int repetitions;
  43   char buf[500];
  44   char *current_buf;
  45   int line_number=0;
  46   int dim1, dim2;
  47
  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;
  52                   ++place[dim1])
  53         counts[dim1][place[dim1]]= 0;
  54     n_below_mean[dim1] =0;
  55     maxplaces[dim1] = 2*MINSCORE;
  56     minplaces[dim1] = -2*MINSCORE;
  57     mean[dim1]=0;
  58
  59     for (dim2=0; dim2<number_score_dim; dim2++)
  60       covar[dim1][dim2]=0;
  61   }
  62
  63
  64 /***** Get a tuple of scores and the count on that tuple. ***/
  65
  66   sc_file = sopen(sc_filename,"r");
  67   fprintf(stderr,"\nOpened file %s",sc_filename);
  68   while (fgets(buf,500,sc_file)) {
  69     line_number++;
  70     /**    fprintf(stderr,"\n%d",line_number);  **/
  71     current_buf = buf;
  72
  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 */
  79
  80     }
  81     sscanf(current_buf,"%d\n",&repetitions);
  82
  83
  84 /****
  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);
  89 ****/
  90
  91     all_places_good=1;   /** If one of scores is out of range, don't include */
  92                          /** any of the scores.                              */
  93
  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;
 102         all_places_good =0;
 103       }
 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;
 110         all_places_good =0;
 111       }
 112     }
 113
 114
 115     if (all_places_good)  {  /** Score the point in  means and covariances. */
 116       n += repetitions;
 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];
 124
 125         for (dim2=0; dim2<number_score_dim; dim2++)
 126           covar[dim1][dim2] +=
 127                 (double)repetitions*score[dim1]*score[dim2];
 128       }
 129
 130     }  /*  Ends count on table after read in each scores. */
 131
 132   }   /* Ends while scores to read in. */
 133
 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);
 136
 137   sclose(sc_file);
 138
 139
 140
 141   /** Now divide through by  number of points  did the expectations over.  **/
 142
 143   for (dim1=0; dim1<number_score_dim; dim1++)
 144     mean[dim1] /= n;
 145
 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];
 150 }
 151
 152
 153
 154
 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],
 158                   int new_dim_number)
 159 {
 160
 161   if (new_dim_number< new_num_dim_table0)
 162     if (old_num_dim_table0 == 1)
 163       return(0);
 164     else
 165       return((int)multi_lib[0][new_dim_number]);
 166
 167   else    /** Second (trimer) table. **/
 168     if (old_num_dim_table1 == 1)
 169       return(old_num_dim_table0);
 170     else
 171       return(old_num_dim_table0 +
 172              (int)multi_lib[1][new_dim_number-new_num_dim_table0]);
 173 }
 174
 175
 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
 186         )
 187 {
 188   int class, dim1, dim2, new_dim1, new_dim2;
 189   double covars_sub[MAX_CLASS_NUMBER][MAX_NUM_SCORE_DIM][MAX_NUM_SCORE_DIM];
 190
 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,
 196                            multi_lib, dim1)];
 197
 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,
 202                            multi_lib,dim1)]
 203             [old_dimension(old_num_dim_table0, old_num_dim_table1,
 204                            new_num_dim_table0, new_num_dim_table1,
 205                            multi_lib,dim2)];
 206
 207     }
 208
 209
 210
 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);
 220
 221   }
 222 #endif
 223 /*********************/
 224
 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);
 228
 229
 230 }
 231
 232
 233
 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. ****/
 242
 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])
 252
 253 {
 254   int dim, class;
 255   int table;
 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;
 259
 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;
 264
 265   int dist;
 266   int all_one =1;
 267 /*******Initialization. ****/
 268
 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;
 274     }
 275
 276
 277
 278 /*  if (!all_greater_than_minus_inf) return;  */
 279
 280 /**************************/
 281
 282
 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.  ***/
 286
 287     if  (normalize(number_score_dim, 1, number_score_dim, normalized_matrix,
 288                    inverse_covars[class], score_diff_from_means[class]) )  {
 289
 290       normalized_value= normalized_matrix[0][0];
 291       like_in_class[class] =  exp(-.5 * normalized_value) /
 292                                ( pow(2*Pi,
 293                                 (double)number_score_dim/2)*
 294                                   sqrt(det[class]) );
 295     }
 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;
 299     }
 300
 301     /*    fprintf(stderr,"\nnormalized_value = %lf",normalized_value);  */
 302
 303     total_like +=  init_prob[class] * like_in_class[class];
 304     if (like_in_class[class] > max_like) max_like = like_in_class[class];
 305   }
 306
 307   for (class=0; all_one && (class<number_classes); class++)
 308     if (init_prob[class] !=1) all_one=0;
 309
 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]/
 315                                                     total_like;
 316     }
 317
 318 }
 319
 320
 321
 322
 323
 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])
 332
 333 {
 334   int class;
 335   FILE *param_file;
 336   char buf[1000];
 337   char *current_buf;
 338   int dim1, dim2;
 339   int read_params_from_file =0;
 340
 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);
 350
 351           fgets(buf,1000,param_file);
 352           current_buf =buf;
 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 */
 359
 360             sscanf(current_buf,"%lf", &means[class][dim1]);
 361           }
 362         }
 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);
 367             current_buf=buf;
 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 */
 374
 375               sscanf(current_buf,"%lf", &covars[class][dim1][dim2]);
 376             }
 377           }
 378         }
 379       }
 380       sclose(param_file);
 381     }
 382
 383 /*************************************************************************/
 384
 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);
 393     }
 394 #endif
 395     return; }   /** All done. **/
 396
 397   if (gauss_param[0] != ',')    /* output to this file */
 398     (param_file=sopen(gauss_param,"w"));
 399
 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]);
 404
 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);
 409
 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);
 415     }
 416   }
 417   if (gauss_param[0] != ',')    /* output to this file */
 418     sclose(param_file);
 419
 420 }
 421
 422
 423 double calc_total_coil_likelihood(double like_in_class[MAX_CLASS_NUMBER],
 424                                   int number_classes, int reg,
 425                                   int Offset_to_Use,
 426                                   double *max_class_like,
 427                                   double init_class_prob[MAX_CLASS_NUMBER])
 428 {
 429   int positive_class;
 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. */
 436
 437   for (positive_class = 0; positive_class<number_classes-1; positive_class++)
 438     coil_init_prob += init_class_prob[positive_class];
 439
 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];
 443
 444       if (total_coil_like > *max_class_like)
 445       *max_class_like = total_coil_like;
 446   }
 447   else   /* Find max scoring total coil reg for reg=7, max offset */
 448     total_coil_like = *max_class_like;
 449
 450
 451   return(total_coil_like);
 452 }
 453
 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. */
 461
 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)
 470
 471 {
 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];
 476   int i;
 477   int reg;
 478   int positive_class;
 479   int dim, tablenum;
 480   int all_one =1, class;
 481   int first_reg;
 482   int none_huge_val =1;  /* Signal for when complementing last class.  */
 483
 484   extern FILE *ftotal_like;   /*  If this is on, then only compute for the */
 485                               /*  max register to save computation time.   */
 486
 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.         */
 491
 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.                                  **/
 496
 497
 498   for (i=0; i< seq_len; i++) {
 499
 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;
 503
 504
 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) ) {
 508         none_huge_val=1;
 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)
 512             none_huge_val=0;
 513         }
 514
 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);
 519
 520
 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];
 525       }
 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];
 529
 530
 531
 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 **/
 540
 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)
 544         if (none_huge_val)
 545           all_raw_scores[number_classes-1][i][reg]=
 546             1 - all_raw_scores[number_classes -1][i][reg];
 547
 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];
 553
 554
 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)
 562                                                 /* No longer raw*/
 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];
 566
 567 /************************************************************************/
 568
 569
 570
 571     }  /** Count on reg.**/
 572   }    /** Count on i **/
 573 }
 574
 575
 576
 577
 578 int compute_dimension2(int tablenum, int libnumber,
 579                        int num_dist[MAX_TABLE_NUMBER],
 580                        int combine_dist[MAX_TABLE_NUMBER])
 581 {
 582   int prev_table;
 583   int dimension=0;
 584
 585   for (prev_table =0; prev_table<tablenum; prev_table++)
 586     dimension+= num_table_dimensions(num_dist[prev_table],
 587                                      combine_dist[prev_table]);
 588
 589   return (dimension + libnumber);
 590 }
 591
 592
 593 int num_table_dimensions(int num_dist_table,
 594                          int combine_dist_table)
 595 {
 596
 597   return ((1-combine_dist_table)*num_dist_table +
 598           combine_dist_table);
 599 }
 600
 601
 602
 603
 604
 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)
 611 {
 612
 613   int i, tablenum, reg;
 614
 615   /** Note that the scores in tables above number_classes-2 is **/
 616   /** the total positive likelihood always, and doesn't need changing. **/
 617
 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];
 625           }
 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];
 629           }
 630 }
 631
 632
 633
 634 void init_totals(double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
 635 {
 636   int dimer_box, trimer_box;
 637
 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;
 641 }
 642
 643
 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)
 651 {
 652
 653   int dimer_box, trimer_box;
 654   double init_prob_dimer, init_prob_trimer;
 655   int i;
 656   int need_to_adjust=0;
 657
 658   for (i=0; i<seqlen; i++) {
 659     need_to_adjust=0;
 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]);
 669
 670         if (total_gauss_like[dimer_box][trimer_box] < -(MAXDOUBLE/2))
 671           need_to_adjust=1;
 672         if (total_gauss_like[dimer_box][trimer_box] > (MAXDOUBLE/2))
 673           need_to_adjust=-1;
 674       }
 675     }
 676
 677 /******* Make sure we don't over flow the double size by adding constants **/
 678 /****    to each sum. ****/
 679 /*************
 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;
 685     }
 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;
 691     }
 692 ******************/
 693   }
 694
 695   if (need_to_adjust)
 696     fprintf(stderr,"\nError:  Overflowed size of double in computeing total like. ");
 697
 698 }
 699
 700
 701
 702
 703 void output_total_like(FILE *ftotal_like,
 704                     double total_gauss_like[GRID_SIZE_DIMER][GRID_SIZE_TRIMER])
 705 {
 706
 707   int dimer_box, trimer_box;
 708   double init_prob_dimer, init_prob_trimer;
 709   double max_like=-HUGE_VAL;
 710
 711
 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];
 716
 717   fprintf(ftotal_like, "The maximum sum of log likelihoods was %lf\n\n",max_like);
 718
 719 /** Normalize since values are near 0 ***/
 720   if (max_like != 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);
 729     }
 730   }
 731
 732 }
 733
 734
 735
 736