#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);
    }
  }

}