/* Interface to sc2seq created by David Wilson and modified by Ethan Wolf */
/* Other interfaces by Ethan Wolf 1995. */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "scio.h"
#include "switches.h"
#include "likelihood.h"
#include "sc.h"
#include "scconst.h"
#include "options.h"
#include "compute_like.h"
#include "interface.h"
#include "stats.h"

#include "debug_printout_flag.h"

#define PIR_LIKES "/tmp/pir.likes"
#define PS_OF_LIKELIHOOD "~/likelihood.ps"
#define PAIR_PIR_TXT "/tmp/pir1.txt" 
#define SINGLE_PIR_TXT "/tmp/pir2.txt"

double many_pprobs[MAX_TABLE_NUMBER][NUM_RES_TYPE][POSNUM], 
         many_pprobp[MAX_TABLE_NUMBER][NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM];


double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];

double (*pprobs)[POSNUM];
double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];


long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM]; 
long ptotals[MAX_TABLE_NUMBER][POSNUM];
long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM]; 
long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];

double getdlog();

#define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))



void *Malloc (n)
     int n;
{
  void *ptr;
  ptr = (void*)malloc(n);
  if (ptr) return ptr;
  fprintf(stderr,"Memory shortage!!\n");
  exit(-1);
}




void switch_tables(int table)
{
  extern int WINDOW, window_length[MAX_TABLE_NUMBER];

  pprobs= many_pprobs[table];
  pprobp= many_pprobp[table];
  WINDOW = window_length[table];
}



/******** Used to subtract off a sequence from the table.  ***/

void recalc_prob(Sequence sequence, int input_table, int mode, long pfreqp2[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
{
  int ProlineFreeWindow = mode & PROLINE_FREE_MODE;

  /* Remove the sequence from the positive table before computing on it. */
  add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table], 
             ptotalp[input_table], 
             sequence.reg, sequence.seq, sequence.seqlen, -1);

  if (!(mode & WEIGHTED_PROBS)) 
    estimate_database_probs2(input_table, pfreqs[input_table], 
			    ptotals[input_table], many_pprobs,
			    pfreqp[input_table], ptotalp[input_table], 
			    many_pprobp, ProlineFreeWindow);
/** Add seq back into the frequence counts. **/
  add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table], 
             ptotalp[input_table], 
             sequence.reg, sequence.seq, sequence.seqlen, 1);

  #ifdef AVG_WEIRD_PROBS 
    calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
  #else
    calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
  #endif  

}



void normalize_gen(double gprobs[NUM_RES_TYPE], 
		   double gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM])
{
  int res1, res2, dist;

  for (res1 =0; res1<Proline; res1++) {
    gprobs[res1] = getdlog(exp(gprobs[res1])/(1-exp(gprobs[Proline])));
    for (res2= 0; res2<Proline; res2++)
      for (dist=0;dist<7;dist++) {
	/** Subtract off for Proline with anything in other position once **/
	/** for each position, then add back on the extra bit subtracted. **/
	gprobp[res1][res2][dist] = getdlog(exp(gprobp[res1][res2][dist]) /
	   ((1-exp(gprobs[Proline])-exp(gprobs[Proline]) + 
	   exp(gprobp[Proline][Proline][dist]))) );
      }
  }

}




void PairCoilData (FILE *fgin, FILE **fpin,
		   int ProlineFree, int table_to_remove_from,
		   char *input_filename, int mode,
		   long pfreqp2[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM])
{
  long gfreqs[AANUM]; long gtotals;  
  long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM]; 

  extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM]; 
  extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
  //extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM]; 
  extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
  extern double SCALE0;


  int i=0;

  if (fgin) {
    readgenfreq(fgin, &gtotals, gfreqs, gtotalp, gfreqp, stderr);
    calcgprobs(gfreqs, gtotals, gprobs);
    calcgprobp(gfreqp, gtotalp, gprobp);
#ifdef NoGenProline
    normalize_gen(gprobs, gprobp);
#endif
  }
  else { setgprob_to_uniform(gprobs, gprobp); }

  #ifdef AVG_WEIRD_PROBS 
    calc_weird_gprob_avg(gprobp, gprobs);
  #else
    calc_weird_gprob_sum(gprobp, gprobs);
  #endif


  for (i=0; i<NUMBER_TABLES; i++)  { /*  Go through all input tables. */
    if (fpin[i] != NULL) {  /* Calculate probabilities for that table.  */
      readposfreq(fpin[i],ptotals[i],pfreqs[i], ptotalp[i], pfreqp[i], stderr);

      if (i == table_to_remove_from)
	    remove_all_input_seq_from_table(pfreqs[table_to_remove_from],
	      ptotals[table_to_remove_from],	          
	      pfreqp[table_to_remove_from],
	      ptotalp[table_to_remove_from], input_filename, -1);  

      estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
				pfreqp[i], ptotalp[i], many_pprobp,
				ProlineFree);
    }
    else {
      fprintf(stderr,"\nCould not make input table number %d.",i);
      exit(-1);
    }
    

#ifdef AVG_WEIRD_PROBS        /* make table.     */
    calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);   
#else
    calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);   
#endif

  }

  switch_tables(0);

}


NEWCOILSinit ()
{
  initialize_relative();
}

#define NEGCCRATIO 30
#define NEGMEAN 0.77
#define NEGSD 0.20
#define CCMEAN 1.63
#define CCSD 0.24
#define SQRT2PI 2.5066282746310002



/** To compute the paircoil likelihood line. */


void likelihood_from_txt(char *txt_filename, double *m, double *b1)
{
  FILE *txt_file;
  int mode2=0;
  double percent_coil=PERCENT_COILED;
  double score1;
  int filenum=1;
  long n[2]={0,0}, n_halfgauss[2]={0,0};
  double sum[2]={0,0}, sum2[2]={0,0};
  double a[2],b[2];  /* ax+b is the likelihood line approximation */
  double mean[2], sd[2], down[2];
  long _counts[2][-4*MINSCORE +1];
  double  _means[2][-4*MINSCORE+1];
  long *counts[2];
  double  *means[2];
  int i, maxplaces[2]= {2*MINSCORE, 2*MINSCORE}, minplaces[2]={-2*MINSCORE,-2*MINSCORE};
  int minplace, maxplace;
  int number_strange=0, number_high=0;
  FILE *devnull, *outfile, *outpsfile;
  int place;

  counts[0] = &(_counts[0][-2*MINSCORE]);
  counts[1] = &(_counts[1][-2*MINSCORE]);
  means[0] = &(_means[0][-2*MINSCORE]);
  means[1] = &(_means[1][-2*MINSCORE]);


  mode2 |= RATIO_TRAP_METHOD;
  mode2 |= PERCENT_COIL_MODE;


 
   for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
     counts[0][i] = 0;
     means[0][i]=0;
    }

  
  txt_file = sopen (txt_filename, "r");
  while (fread(&score1,sizeof(double),1,txt_file)) {
    if (score1<MINSCORE || score1>-(MINSCORE)) {
/*
  fprintf(stderr,"Strange score: %lf\n",score1);
  exit(-1);
*/
      number_strange++;
    }
    place = ( (int)((score1-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
    if (place > -2*MINSCORE) {
      fprintf(stderr, "Out of bounds: %d.\n", place);
      fprintf(stderr, "Score: %lf\n",score1);
      /* exit(-1); */
      counts[0][-2*MINSCORE]++;
      means[0][-2*MINSCORE]++;
      number_high++;
      maxplaces[0]=-2*MINSCORE;
    }
    else if (place < 2*MINSCORE) {
      fprintf(stderr, "Out of bounds: %d.\n", place);
      fprintf(stderr, "Score: %lf\n",score1);
      counts[0][2*MINSCORE]++;
      means[0][2*MINSCORE]++;
      minplaces[0]=2*MINSCORE;
	/* exit(-1); */
    }

    else {
      counts[0][place]++;
      means[0][place] += score1;
      if (place>maxplaces[0]) maxplaces[0]=place;
      if (place<minplaces[0]) minplaces[0]=place;
      ++n[0];
    }
  }
  
   fprintf(stderr,"\n%d out of bounds scores > %d  were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE)); 
   fprintf(stderr,"\n%d strange scores < %d were scored as %d\n",number_strange-number_high,MINSCORE, MINSCORE);
   sclose(txt_file);
  
   /** Compute the means of the scores in each box.          **/ 
   for (i= minplaces[0]; i<=maxplaces[0]; ++i) 
     if (counts[0][i] != 0) 
       means[0][i] /= counts[0][i];  /** Each bucket gets average score*/

   
   /*** We need a mean***/

   mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
			 maxplaces[0],n[0]);

  
  /*** Calculate gaussian fit to scores below the mean. **/
  txt_file=sopen(txt_filename, "r");
  while (fread(&score1,sizeof(double),1,txt_file)) 
    if ( (score1 > MINSCORE) && (score1 < mean[0])) {
      ++n_halfgauss[0];
      sum2[0] += (score1-mean[0])*(score1-mean[0]);
    }
  sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
  down[0]= (double) 2*n_halfgauss[0]/n[0]; 

  fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
  sclose(txt_file);

  minplace=minplaces[0];  maxplace=maxplaces[0];

  devnull=sopen("/dev/null","r+");
  outfile=sopen(PIR_LIKES,"w"); 
  outpsfile=sopen(PS_OF_LIKELIHOOD,"w"); 
/*  outfile=sopen("~/pir.out","w"); */
  likelihood(n,mode2,  filenum, 0.0 , minplace, maxplace, means, 
             counts, mean, sd, down, outfile,outpsfile,devnull,devnull,
	     a,b,.20,.80,.03);
  sclose(outfile);
  sclose(outpsfile);
  {
    char clean_command[200];
     
    strcpy(clean_command, "rm -f ");

/*** The following removes the pir.txt output file. */    
    system(strcat(clean_command,txt_filename));   
	
    }                            /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
                                 /** WRITE TO SAME FILESNAMES.             **/
  {
    char clean_command[200];
     
    strcpy(clean_command, "rm -f ");
    
    system(strcat(clean_command,PIR_LIKES));   
	
    }  

  *m=a[0]; *b1=b[0];
}





/*  Look in file .likelihood for likelihood line.   If not there, compute */
/*  it from pir.                                               */
void get_likelihood_line(char likelihoods[MAX_TABLE_NUMBER][MAXLINE], double *m, double *b, double *m_single, double *b_single, char lib[MAXFUNCTNUM], int functnum, char *pir_name, int mode, int number_tables)
{
 char *tmp;
 int libnumber=0;
 int  i, cnt=0;
 int current_lib;
 char buf[80];

 FILE *pir;
 FILE *linefile;
 FILE *txt_file1, *txt_file2;
 double score, scs[MAXSEQLEN][POSNUM];
 Sequence sequence;
 char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
 char title[500], code[500];
 int seqlen;
 int likefile;
 int found_like_lib=0;  int found_single=0;

 sequence.seq=seq;
 sequence.reg=reg;
 sequence.title=title;
 sequence.code=code;


 for (i=0; i<functnum; i++)
   libnumber |= (int) pow(2,lib[i]);
#ifdef DEBUG_PRINTOUT
 fprintf(stderr,"\nlibnumber= %d",libnumber);
#endif

  /*  What directory to look for paircoil defauts in.  Set environment  */
  /*  variable LIKELINE_DIR to be a path to this file.               */
 for (likefile=0; likefile <number_tables; likefile++) {
   current_lib = -1;
   found_like_lib =0;
   found_single =0;
#ifdef DEBUG_PRINTOUT
   fprintf(stderr,"\nLikelihood %d=%s",likefile,likelihoods[likefile]);
#endif
   if (likelihoods[likefile][0] == ',') { /* Don't bother loading likefile */
#ifdef DEBUG_PRINTOUT
     fprintf(stderr,
	     "\n Not doing likelihood conversion for paircoil with table %d\n",
	     likefile);
#endif
   }
   else {
     if (linefile=sopen(likelihoods[likefile],"r")) {
       while ((fgets(buf,80,linefile)) && !(found_single && found_like_lib))  {
	 if (buf[0] == '0') {
	   sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
	   found_single=1;
	 }
	 else if (!found_like_lib) {
	   sscanf(buf, "%d %lf %lf",&current_lib, &m[likefile], &b[likefile]);
	   if (current_lib == libnumber)  found_like_lib=1;;
	 }
       }
       sclose(linefile);

       if (!(found_like_lib && found_single)) {
	 if (!found_single)
	   fprintf(stderr,"\nSingles likelihood function %d unknown.",
		   likefile);

	 if (current_lib != libnumber )  /* Wasn't in file so should create it
				          by running on pir.  Add this later.*/
	   fprintf(stderr, "\nLikelihood function %d unknown...", likefile);


	 if (pir_name[0]==',') {
	   fprintf(stderr,"\nUnable to continue.  Please input the location of the pir in the .paircoil file.");
	   exit(-1);
	 }
	 else {
	   fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
	 
	   pir = sopen(pir_name, "r");

	   if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
	   if (!found_single)   txt_file2 = sopen(SINGLE_PIR_TXT,"w");
	   
	   for (cnt=1;getseq2(pir,&sequence); cnt++)
	     if (sequence.seqlen >= MINWINDOW) {
	       seqlen= sequence.seqlen;
	       if (!(cnt %10)) {
		 putc('.',stderr);
		 fflush(stderr);
	       }
	       
	       WINDOW=PAIRWINDOW;
	       if (!found_like_lib) {
		 switch_tables(likefile);
		 scseqadj(lib,seq,seqlen,scs,offsets, &score,
			  mode,0);
		 txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
	       }
	       if (!found_single) {
		 switch_tables(likefile);

		 scorestock(seq,seqlen,scs,offsets, &score, 
			    mode & PROLINE_FREE_MODE, 0);   
		 /* Score stock = singles method using our table. */
		 txt_output(sequence,0,txt_file2, scs, NULL,offsets,0);
	       }
       
	     }  /* Ends if sequence > MINWINDOW, and hence the for */
	         /* loop that reads from pir. */
	   sclose(pir);


	   if (!found_like_lib) {    /* Now compute the lib likelihood line. */
	     sclose(txt_file1);	  

	     likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
	     
	     linefile=sopen(likelihoods[likefile],"a");
	     fprintf(linefile,"%d %lf %lf\n",
		     libnumber,m[likefile], b[likefile]);
	 
	     fprintf(linefile,"\t (lib =");  
	     /* Now print human readable lib number. */
	     for (i=0; i<functnum; i++)
	       fprintf(linefile," %d",lib[i]);
	     fprintf(linefile," )\n");

	     sclose(linefile);
	   }

	   if (!found_single) { /*  Now compute the singles method like line */
	     sclose(txt_file2);	  

	     likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
			       &b_single[likefile]);
		       
	     linefile=sopen(likelihoods[likefile],"a");
	     fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile], 
		     b_single[likefile]);
	 
	     fprintf(linefile,"\t (lib =");  
	     /* Now print human readable lib number. */
	     fprintf(linefile," Singles Method");
	     fprintf(linefile," )\n");

	     sclose(linefile);
	   }
	 }      /* Ends else from if (pir_name == ',') */

       }       /*  Ends section for file not having lib or singles like line */
  
     }         /*  Ends section where search likelihoods[likefile] file. */
     else   {  
       fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
       exit(-1);
     }
   }  /** Ends else the likelihood file exists to load a likelihood line... **/

 }          /*  Ends count through tables=linefiles.  */
}



double stocksc2likelihood (double sc)
{
  double gg, gcc;

  sc = exp(sc/28);
  gg = (sc - NEGMEAN)/NEGSD;
  gcc = (sc - CCMEAN)/CCSD;
  gg = exp((-gg * gg) / 2) / NEGSD / SQRT2PI;
  gcc = exp((-gcc * gcc) / 2) / CCSD / SQRT2PI;
  return(gcc / (NEGCCRATIO * gg + gcc));
}



/*  Likelihood line is mx +b */
double sc2likelihood (double sc, double m, double b)
{
  double likelihood;

  likelihood = m*sc + b;
 
  if (likelihood>1) return(1);
  else if (likelihood<0) return(0);
  else return(likelihood);
}


static int current_seq_number=0;

void PairCoilScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double m, double b, double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use, int raw_score)
{
  int i,j;
  double sc2scores[MAXSEQLEN][POSNUM];

  int ProlineFreeWindow= mode & PROLINE_FREE_MODE;

  double  frac_coiled_this_table;
  extern double init_class_prob[MAX_CLASS_NUMBER];

  frac_coiled_this_table = init_class_prob[table];


  switch_tables(table);

  /* score the sequence using sc2seq */
  scseqadj(lib,
	    sequence.seq, sequence.seqlen,
	    sc2scores,
	    offsets,
	    maxscore, mode, 0);


  if (!raw_score)  
    if (mode & USE_LIKE_LINE)
      *maxscore= sc2likelihood(*maxscore,m,b);
    else {
      *maxscore= frac_coiled_this_table * exp(*maxscore);
    }

  for (i=0; i<sequence.seqlen; ++i) {
    /*  Do the most probabable register if offset_to_use is -1 or 7. */
    if ((offset_to_use == -1) || (offset_to_use ==7)) 
      if (raw_score)  /*** Return rawscore. **/
	scores[i][7] = sc2scores[i][offsets[i]];
      else
	if (mode & USE_LIKE_LINE)
	  scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
	else   {
	  scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
	  if (scores[i][7] >1) scores[i][7]=1;
	}

    else    /*  Offset to use is 0-6. */
      if (raw_score)
	scores[i][7] = sc2scores[i][offset_to_use];
      else 
	if (mode & USE_LIKE_LINE)
	  scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
	else {
	  scores[i][7] = frac_coiled_this_table *
	                  exp(sc2scores[i][offset_to_use]);
	  if (scores[i][7] >1) scores[i][7]=1;
	}
    
    /* do other registers */
    for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
			   i+j%POSNUM  **/
      if (raw_score)
	scores[i][(i+j)%POSNUM] = sc2scores[i][j];
      else
	if (mode & USE_LIKE_LINE)
	  scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
	else  {
	  scores[i][(i+j)%POSNUM] = frac_coiled_this_table * 
	                                exp(sc2scores[i][j]);
	  if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
	}
    }
  }

}





/***************************************************************************/

/** Score using the singles method on our tables.   **/
double *SingleCoilScore(int mode, Sequence sequence, double m, double b,
			double *maxscore, int table, 
			double scores[MAXSEQLEN][POSNUM+1], int offset_to_use,
		        int raw_score)
{
  int i,j;
  double singlescores[MAXSEQLEN][POSNUM];
  char offsets[MAXSEQLEN];

  int ProlineFreeWindow= mode & PROLINE_FREE_MODE; 

  double  frac_coiled_this_table;


  switch_tables(table);

  /* score the thing */
  scorestock(
	    sequence.seq, sequence.seqlen,
	    singlescores,
	    offsets,
	    maxscore,ProlineFreeWindow, 0);  /** Score using our table. */


  if (!raw_score)
    if (mode & USE_LIKE_LINE)
      *maxscore= sc2likelihood(*maxscore,m,b);
    else {
      /*    fprintf(stderr,"\nmaxscore = %lf", *maxscore);
       */
      *maxscore= frac_coiled_this_table * exp(*maxscore);
      /*    fprintf(stderr,"  maxlike  = %lf",  *maxscore);
       */
    }

  for (i=0; i<sequence.seqlen; ++i) {
    /*  Do the most likely register.  */
    if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
      if (!raw_score)
	if (mode & USE_LIKE_LINE)
	  scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b); 
	else {
	  scores[i][7] = frac_coiled_this_table* 
	              exp(singlescores[i][offsets[i]]);
	  if (scores[i][7] >1) scores[i][7]=1;
	}
      else scores[i][7] = singlescores[i][offsets[i]];

    else                  /*  Use offset_to_use. */
      if (!raw_score)
	if (mode & USE_LIKE_LINE)
	  scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
	else {
	  scores[i][7] = frac_coiled_this_table* 
	             exp(singlescores[i][offsets[i]]);
	  if (scores[i][7] >1) scores[i][7]=1;
	}
      else scores[i][7] = singlescores[i][offset_to_use];

    /* do the other registers */
    for (j=0; j<7; ++j)
      if (!raw_score)
	if (mode & USE_LIKE_LINE)
	  scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b); 
	else {
	  scores[i][(i+j)%POSNUM] = frac_coiled_this_table*
	                                 exp(singlescores[i][j]);
	  if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
	}
      else scores[i][(i+j)%POSNUM] = singlescores[i][j];  /*Raw_score*/
    
  }


}
/***************************************************************************/



void MultiCoilDimensionScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use)
{
  int i,j;
  double sc2scores[MAXSEQLEN][POSNUM];

  int ProlineFreeWindow= mode & PROLINE_FREE_MODE;


  switch_tables(table);

  /* score the sequence using sc2seq */
    scseqadj(lib,
            sequence.seq, sequence.seqlen,
            sc2scores,
            offsets,
            maxscore, mode, 0);

/*** Do a SingleCoil type score. ****/
/***
    scorestock(sequence.seq, sequence.seqlen,
	    sc2scores,offsets,
	    maxscore,ProlineFreeWindow, 0);  **/ /** Score using our table. */


  for (i=0; i<sequence.seqlen; ++i) {
    /*  Do the most probabable register if offset_to_use is -1 or 7. */
    if ((offset_to_use == -1) || (offset_to_use ==7)) 
      scores[i][7] = sc2scores[i][offsets[i]];
    else    /*  Offset to use is 0-6. */
      scores[i][7] = sc2scores[i][offset_to_use];

    /* do other registers */
    for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
                           i+j%POSNUM  **/
      scores[i][(i+j)%POSNUM] = sc2scores[i][j];
    }
  }

}





/***************************************************************************/

/*** Score using the singles method on stock table and likelihood.  **/
void NEWCOILSScore(int mode, Sequence sequence, double *maxscore,
		   double scores[MAXSEQLEN][POSNUM+1], 
		   int offset_to_use)  
{
  int i,j;
  double stockscores[MAXSEQLEN][POSNUM];
  char offsets[MAXSEQLEN];


  int ProlineFreeWindow= mode & PROLINE_FREE_MODE; 

  /* score the thing */
  scorestock(
	    sequence.seq, sequence.seqlen,
	    stockscores,
	    offsets,
	    maxscore,ProlineFreeWindow, 1);  /** Score using STOCK table. */


  *maxscore= stocksc2likelihood(*maxscore);


  for (i=0; i<sequence.seqlen; ++i) {
    /*  Do the most likely register.  */
    if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
      scores[i][7] = stocksc2likelihood(stockscores[i][offsets[i]]); 
    else                  /*  Use offset_to_use. */
      scores[i][7] = stocksc2likelihood(stockscores[i][offset_to_use]);
    /* do the other registers */
    for (j=0; j<7; ++j)
      scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]); 
  }


}


/***********************************************************/






void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1], 
		  int offset_to_use, int preprocessor)
{
  int i,j;

  if (!sequence.reg) 
    fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");


  for (i=0; i<sequence.seqlen; i++) {
    for (j=0; j<POSNUM; j++)
      if (!sequence.reg) scores[i][j] =0;  /* No coils in input file. */

      else if ( (offset_to_use ==7) ||  /* max offset is pos file reg */
	       ((offset_to_use == -1) && !preprocessor) ) 
	if (j== sequence.reg[i])  
	  scores[i][j] = 1;
	else scores[i][j]=0;

      else if ( (offset_to_use == -1) &&  preprocessor) /* Do all registers */
	if (sequence.reg[i] != '.')
	  scores[i][j] = 1;
	else scores[i][j]=0;

      else  /* just give likelihood 1 to offset_to_use. */
	if ( (sequence.reg[i] != '.') && (j ==  (i+offset_to_use)%POSNUM))
    	  scores[i][j] = 1;
	else scores[i][j]=0;

    if (!sequence.reg) scores[i][7]=0;
    else if (sequence.reg[i] != '.') 
      scores[i][7]=1;
    else scores[i][7]=0;
  }

}
          






/********************************************************/

/*** This averages score over coil, where a register shift does not start **/
/*** a new coil if start_at_reg_shift =0.                          ****/
/*** This defines whether in a coil by if likelihood is > bound.    ***/
/** Value of 0 for avg_max means do average, 1 means do max.        ***/
/*** For max score, find the maximum scoring residue in the CoilLike ***/
/*** and take the corresponding StructureScore.                      ***/
/*** To compute the max version, average over all residues in the coil **/
/*** that score the max coil likelihood.  ***/
/*** BE VERY  CAREFUL in making changes!!! i.e. need to count down on k **/
/*** since offset 7 needs to make use of previous offsets, so don't want **/
/*** them to have been changed to coil score yet.  Also, we sometimes   ***/
/*** pass in CoilLike as the StructureLike too, in which case have to be**/
/*** doubly careful.   There is also ONE place where make a reference   **/
/** to offset 7 when dealing with other registers, so had to save an    **/
/** old version of offset 7  in case it changes (when StructureLike     **/
/** is the same as CoilLike.                      **/
/*** ACTUALLY, all that should be moot now that we changed it so residue **/
/*** scores are always saved in tablenum +3, so pass those in and save   **/
/*** coil scores in tablenum. ***/
void average_score_over_coils3(Sequence seq, 
			       double StructureResLike[MAXSEQLEN][POSNUM+1],
			       double StructureCoilLike[MAXSEQLEN][POSNUM+1],
			       double CoilLike[MAXSEQLEN][POSNUM+1],
			       int offset_to_use,
			       char offsets[MAXSEQLEN],
			       int start_at_reg_shift, double bound,
			       int avg_max)
{
  int i=0, j,k;   /* k is the offset. */
  int start_coil[POSNUM+1];
  double total_coil_score[POSNUM+1];
  double avg_coil_score[POSNUM+1];
  double max_CoilLike[POSNUM+1];
  double total_max_StructLike[POSNUM+1];
  int number_max_residues[POSNUM+1];
  double max_StructLike[POSNUM+1];
  int regist;

  double total_coil_length[POSNUM +1];  /* A weighted length. */
  double CoilLike_MaxOffset[MAXSEQLEN];



  for (j=0; j<POSNUM+1; j++) {
    start_coil[j] = -1;
    total_coil_score[j]=0;
    avg_coil_score[j]=0;
    total_coil_length[j]=0;
    max_CoilLike[j]=0;
    total_max_StructLike[j]=0;
    number_max_residues[j]=0;
  }
  
  while (i<=seq.seqlen) {
    if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];

    for (k=POSNUM; k>=0; k--) {  /** Count down since k=7 needs to use **/
                                 /** scores from other offset (i.e see **/
                                 /** use  of  regist_to_use).          **/
      if (k!=7) regist = (i+k)%POSNUM;
      else regist =7;
      
      /*  Consider <=  bound  to be out of coil */
      if ( ( CoilLike[i][regist] <= bound) || (i==seq.seqlen) ||
	  ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift && 
           (CoilLike_MaxOffset[i-1]> bound) ) )  {       
	/** If not in a coil (value of -HUGE_VAL signals not in coil).  */
	if (start_coil[k] != -1)  {  /* Ended coil at i-1. */
	  avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
	  max_StructLike[k] = total_max_StructLike[k]/number_max_residues[k];
	  for  (j=start_coil[k]; j<i; j++) {
	    if (k!=7)
	      if (avg_max == 0)
		StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
	      else
		StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
	    else
	      if (avg_max == 0)
		StructureCoilLike[j][7]= avg_coil_score[k];
	      else
		StructureCoilLike[j][7]= max_StructLike[k];
	  }
	}
	total_coil_score[k] =0;      /***Resets for "not in coil"  ***/
        max_CoilLike[k]=0;

	if (i != seq.seqlen)   /* Consider <= bound to be out of coil */
	  if  ( CoilLike[i][regist] <= bound)  {  /* Not in coil anymore. */
	    StructureCoilLike[i][regist]=0;  /** Make coil score under bound */
                                             /** be 0**/

/***	    StructureCoilLike[i][regist]= StructureResLike[i][regist];  **/
	                              /** If not in coil above bound make   **/
	                              /** the coil like the same as the     **/
	                              /** regist like. **/
	    start_coil[k] = -1;
	    total_coil_length[k] =0;
            max_CoilLike[k]=0;
	  }
	  else {  /* Changing registers. */
	    start_coil[k]= i;
	    if (CoilLike[i][regist] > max_CoilLike[k]) {
	      max_CoilLike[k] = CoilLike[i][regist];
	      total_max_StructLike[k] = StructureResLike[i][regist];
	      number_max_residues[k] =1;
	    }
	    else if (CoilLike[i][regist] == max_CoilLike[k]) {
	      total_max_StructLike[k] += StructureResLike[i][regist];
	      number_max_residues[k]++;
	    }	

	    if (1) {   /** Do weighted_avg **/
	      int reg_to_use=regist;
	      if (regist ==7) 
		reg_to_use= (i+offsets[i])%POSNUM;

	      total_coil_length[k]= CoilLike[i][reg_to_use];
	      total_coil_score[k] = StructureResLike[i][regist]*
		CoilLike[i][reg_to_use];
/*********
	      if (total_coil_length[k] ==0)
		fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
			i, 'a'+reg_to_use);
*************/
	    }
	  }
      }   /** Ends if ending a coil.  **/

      else {    /* If in a coil.  */
	if (CoilLike[i][regist] > max_CoilLike[k]) {
	  max_CoilLike[k] = CoilLike[i][regist];
	  total_max_StructLike[k] = StructureResLike[i][regist];
	  number_max_residues[k] =1;
	}
	else  if (CoilLike[i][regist] == max_CoilLike[k]) {
	  total_max_StructLike[k] += StructureResLike[i][regist];
	  number_max_residues[k]++;
	}	
	
	if (1) {  /* Weighted average scores.  */
	  int reg_to_use=regist;
	  if (regist ==7) 
	    reg_to_use= (i+offsets[i])%POSNUM;

	  total_coil_score[k] += StructureResLike[i][regist]*
	                                   CoilLike[i][reg_to_use];
          total_coil_length[k] += CoilLike[i][reg_to_use];
	}
	if (start_coil[k] == -1) 
	  start_coil[k] =i;
      }
    }

    i++;
  }

 /*  Need the following hack to get the correct register in combo off method.*/
/************
  if ( (offset_to_use == 7) && (!start_at_reg_shift)) 
    for (i=0; i<seq.seqlen; i++) {
      for (j=0; j<POSNUM; j++)
	StructureCoilLike[i][j]= -HUGE_VAL;
      StructureCoilLike[i][(i+offsets[i]) %POSNUM] = StructureCoilLike[i][7];
    }
  else
*************/

/** Need the following hack to get scores for max offset method. **/
  if (offset_to_use == -1)
    for (i=0; i<seq.seqlen; i++) 
      StructureCoilLike[i][7]= StructureCoilLike[i][(i+offsets[i]) %POSNUM];

}