#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 PairCoilDifferDimension(int mode, Sequence sequence, 
			     char lib[MAXFUNCTNUM], double *maxscore,
			     char preproc_offsets[MAXSEQLEN],
			     double preproc_scores[MAXSEQLEN][POSNUM+1],
		  double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
		  char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN], int dimension,
			     int offset_to_use,
			     int number_differ_score_dim, double bound)
{
  int i,j, pass,pos;
  double sc2scores[MAXSEQLEN][POSNUM];
  char diff_offsets[MAXSEQLEN];
  extern int which_differentiator;
  extern int differ_window_length;
  int number_passes=1;
  int start_coil =-1;

  double coil_score[POSNUM];

  /** which_differentiator = 1 means do window score starting at residue**/
  /** If it is 2 means find both max and min window.  ******/
  /** Passing 1 into scseqadj makes it return residue scores so the **/
  /** residue gets the window score for the window starting at that **/
  /** residue.  2 returns the max window score containing the residue **/
  /** and 3 returns the min window score. **/

  WINDOW = differ_window_length;  /** TEMPORARY.  **/

  if (which_differentiator ==2) number_passes=2;
  else number_passes=1;

  for ( pass=0; pass<number_passes; pass++) {
    if (which_differentiator ==3 || which_differentiator==4) { 
                                /** get coiled regions from preproc. **/
      i=0;
      while(i<=sequence.seqlen) {
	if ( (start_coil >= 0) && ( (i==sequence.seqlen) || 
	   (preproc_scores[i][7]<=bound) ||  /* End coil */
	   ( (i>0) && 
	    (preproc_offsets[i] != preproc_offsets[i-1])) )) { /*regshift*/
/*	  WINDOW= i-start_coil; */
/*******
	  scseqadj(lib,
		   &(sequence.seq[start_coil]), i-start_coil,
		   &(sc2scores[start_coil][0]),
		   &(diff_offsets[start_coil]),
		   maxscore, mode, 1);  
*******/
/**          score_differentiator_window(lib,
	       &(sequence.seq[start_coil]), i-start_coil,
	       coil_score,mode, 1);
**/
          score_differentiator_window(lib,
	       &(sequence.seq[start_coil]), i-start_coil,
	       coil_score,mode, 1);

/*	 fprintf(stderr,"\n start_coil = %d, end_coil =%d",start_coil, i-1); */

	  for(j=start_coil; j<i; j++)
	    for(pos=0;pos<POSNUM;pos++)
	      if (which_differentiator==3) 
		sc2scores[j][pos]= coil_score[pos];
	      else if (which_differentiator==4) 
		sc2scores[j][pos] = coil_score[pos]/(i-start_coil);

	  start_coil = -1;
	  
	}
	
	if ( (preproc_scores[i][7] > bound) && (start_coil ==-1))
	  start_coil =i;
	i++;
      }
    }
    else
      scseqadj(lib,
	       sequence.seq, sequence.seqlen,
	       sc2scores,
	       diff_offsets,
	       maxscore, mode, pass + which_differentiator); 
                               /* Get differentiator score.*/
                              /* pass + which_diff = 2 is max_window_score **/
                              /*                     3 is min_window_score **/

    for (i=0; i<sequence.seqlen; ++i) {
      /* Use the offset given by the preproc if offset_to_use is -1 or 7 */
      if ((offset_to_use == -1) || (offset_to_use ==7))  {
	all_scores[dimension+ pass*number_differ_score_dim/2][i][7] = 
	  sc2scores[i][preproc_offsets[i]];
	all_offsets[dimension+ pass*number_differ_score_dim/2][i] = 
	  preproc_offsets[i];
      }
      else {    /*  Offset to use is 0-6. */
	all_scores[dimension+ pass*number_differ_score_dim/2][i][7] = 
	  sc2scores[i][offset_to_use];
	all_offsets[dimension+ pass*number_differ_score_dim/2][i]= 
	  offset_to_use;
      }
    

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




/******************************************************************/
void zero_out_non_coils( 
		     double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
	             double preproc_like[MAXSEQLEN][POSNUM+1],
		     int seqlen, double bound)
{

  int res, reg, olig;

  for(olig=0; olig<2; olig++) 
    for (res=0; res<seqlen; res++)
      for (reg=0;reg<=POSNUM;reg++)
	if (preproc_like[res][7] > bound)
	  all_scores[olig][res][reg] = all_scores[olig+3][res][reg];
	else
	  all_scores[olig][res][reg] =0;
}


void zero_out_bad_windows(int differ_window_length, 
		     double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
	             double preproc_like[MAXSEQLEN][POSNUM+1],
		     int seqlen, double bound)
{
  int i,j;
  int last_bad_window_start=-1;
  int olig_state, pos;
/*** For now just look at score of first and last residue in window to **/
/**  see if they are above bound.  To really do it right should look at */
/** min score in window (using min_over_windows on preproc_like), but   */
/** the format expected by min_over_windows would require the POSNUM+1  */
/** indice of preproc_like be in the first, not second position....     */

  for (i=0; i<seqlen; i++)
    if (preproc_like[i][7] < bound)  { /** Zero out windows. **/
      for (j=MAX(last_bad_window_start+1,i-differ_window_length); j<=i; j++)
	for(olig_state=0; olig_state<2; olig_state++)
	  for(pos=0; pos<=POSNUM; pos++) 
	    all_scores[olig_state][j][pos]=0;

      last_bad_window_start=i;
    }
    else   /* Do this in case bound has changed....using stored value in +3 */
      for(olig_state=0; olig_state<2; olig_state++)
	for(pos=0; pos<=POSNUM; pos++)
	  all_scores[olig_state][i][pos]=all_scores[olig_state+3][i][pos];
}