#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"

void average_score_over_coils2(Sequence seq, 
                               double scores[MAXSEQLEN][POSNUM+1],
                              double PreprocLike[MAXSEQLEN][POSNUM+1],
                              int CoilDiffMethod, int offset_to_use,
                              char offsets[MAXSEQLEN],
                              int start_at_reg_shift)
{
  int i=0, j,k;
  int start_coil[POSNUM+1];
  double total_coil_score[POSNUM+1];
  double avg_coil_score[POSNUM+1];
  int regist;

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

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

  if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
                                              /* to max offset locally. */

  while (i<=seq.seqlen) {
    for (k=0; k<=POSNUM; k++) {
      if (k!=7) regist = (i+k)%POSNUM;
      else regist =7;

      /*  Consider < 1/10 percent to be zero */
      if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
          ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift && 
           (PreprocLike[i-1][7]>= .001) ) )  {       
      /** 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];
          for  (j=start_coil[k]; j<i; j++) {
            if (k!=7)
              scores[j][(j+k)%POSNUM]= avg_coil_score[k];
            else
              scores[j][7]= avg_coil_score[k];

          }
        }
        total_coil_score[k] =0;
void average_score_over_coils2(Sequence seq, 
                               double scores[MAXSEQLEN][POSNUM+1],
                              double PreprocLike[MAXSEQLEN][POSNUM+1],
                              int CoilDiffMethod, int offset_to_use,
                              char offsets[MAXSEQLEN],
                              int start_at_reg_shift)
{
  int i=0, j,k;
  int start_coil[POSNUM+1];
  double total_coil_score[POSNUM+1];
  double avg_coil_score[POSNUM+1];
  int regist;

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

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

  if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
                                              /* to max offset locally. */

  while (i<=seq.seqlen) {
    for (k=0; k<=POSNUM; k++) {
      if (k!=7) regist = (i+k)%POSNUM;
      else regist =7;

      /*  Consider < 1/10 percent to be zero */
      if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
          ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift && 
           (PreprocLike[i-1][7]>= .001) ) )  {       
      /** 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];
          for  (j=start_coil[k]; j<i; j++) {
            if (k!=7)
              scores[j][(j+k)%POSNUM]= avg_coil_score[k];
            else
              scores[j][7]= avg_coil_score[k];

          }
        }
        total_coil_score[k] =0;
        if (i != seq.seqlen)   /* Consider < 1/10 percent to be 0 */
          if  ( PreprocLike[i][regist] < .001)  {  /* Not in coil anymore. */
            scores[i][regist]=0;
            start_coil[k] = -1;
            total_coil_length[k] =0;
          }
          else {  /* Changing registers. */
            start_coil[k]= i;
            if (CoilDiffMethod==0) {
              total_coil_length[k] = 1;
              total_coil_score[k] = scores[i][regist];
            }
            else if (CoilDiffMethod ==1) {
              total_coil_length[k]= PreprocLike[i][regist];
              total_coil_score[k] = scores[i][regist]*
                PreprocLike[i][regist];
            }
          }
      }   /** Ends if ending a coil.  **/

      else {    /* If in a coil.  */
        if (CoilDiffMethod ==0) {   /* Average scores over coil residues. */
          total_coil_score[k] += scores[i][regist];
          
          total_coil_length[k] += 1;
        }
        else if (CoilDiffMethod == 1) {  /* Weighted average scores.  */
          total_coil_score[k] += scores[i][regist]*
                 PreprocLike[i][regist];
          total_coil_length[k] += PreprocLike[i][regist];
        }
        if (start_coil[k] == -1) 
          start_coil[k] =i;
      }
    }

    i++;
  }

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



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

void PairCoilDiffer(int mode, Sequence sequence, char lib[MAXFUNCTNUM], 
		      double *maxscore,char offsets[MAXSEQLEN], 
		      double preproc_like[MAXSEQLEN][POSNUM+1], 
		      int offset_to_use, int table, 
		      double scores[MAXSEQLEN][POSNUM+1],
		      int Avg_Coil_Score, int CoilDiffMethod,
		      int start_coil_at_reg_shift)

{  
  int i,j,k;
  double sc2scores[MAXSEQLEN][POSNUM];
  char diff_offsets[MAXSEQLEN];
  double total_coilscore[POSNUM+1];
  int coil_length[POSNUM +1];
  int like2or3;
  
  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];


  int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
  
  if (table ==0) like2or3= 2;
  else like2or3=3;

  for (i=0; i<=POSNUM; i++) {
    coil_length[i] = 0;
    total_coilscore[i] =0;
  }



  /* score the sequence using sc2seq */
  scseqadj(lib,
	   sequence.seq, sequence.seqlen,
	   sc2scores,
	   diff_offsets,
	   maxscore, mode, 1);  /* Get differentiator score.*/





/*  *maxscore= exp(*maxscore);   */

  *maxscore =0;

  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 (like2or3 ==2)  
	scores[i][7] = 1 - exp(sc2scores[i][offsets[i]]); 
      else 
	scores[i][7] =  exp(sc2scores[i][offsets[i]]); 
     
    else    /*  Offset to use is 0-6. */
      if (like2or3 ==  2) 
	scores[i][7] = 1 - exp(sc2scores[i][offset_to_use]);
      else
	scores[i][7] = exp(sc2scores[i][offset_to_use]);
                /*****Trimer like is 1 - dimerlike ***/

    if (preproc_like[i][7] <.001)   /* Consider < 1/10 percent to be 0 */
      scores[i][7]=0;  /* Don't score non-coiled regions. */



    for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
			     i+j%POSNUM  **/
      if (preproc_like[i][(i+j)%POSNUM]> 0) { 
	if (like2or3==2) 
	  scores[i][(i+j)%POSNUM] = 1 - exp(sc2scores[i][j]);
	else 
	  scores[i][(i+j)%POSNUM] = exp(sc2scores[i][j]);
      }
      else 
	scores[i][(i+j)%POSNUM]= 0;
    }
    
  }
	

  if (Avg_Coil_Score)
    average_score_over_coils2(sequence, scores, preproc_like,
			      CoilDiffMethod, offset_to_use,offsets,
			      start_coil_at_reg_shift);



  for (i=0; i<sequence.seqlen; i++)
    if (scores[i][7] > *maxscore) *maxscore = scores[i][7];


}