/*  Ethan Wolf 1995 */

#include <stdio.h>
#include "scconst.h"
#include "math.h"

/* Compute scores for residues to be the max window score of windows containing
 * that residue. Windows are defined by where they start.  
 * windowscore[0..seqlen-window] are the window scores so window[i] is
 *     the score for the window starting in position i.
 *  The answer is placed in maxima[0..seqlen-1]
 */
/* Note that the last windowscore is windowscore[seqlen-1-(window-1)]. */



/*  A rough idea of the dynamic programming algorithm presented here:
*
*    Suppose have seq of scores (a_1,a_2,....,a_n)
*   and you want to get values  (m_1,m_2,....,m_n) from the window scores 
*   m_i= max{a_{i-(window-1)}, a_{i-window+1},...., a_i}.
*
*   Then divide the seq into blocks:
*   ( [a_1,....,a_window], [a_{window+1},...,a_{2window}],...,[....,a_n] )
*   and compute block_maxima[i]= max{a_start, a_{start+1},...., a_i},
*   which is the max of the first i scores in the block.

*   Also for k in a block ending at end define
*   prev_block_maxima[k]= max{a_end, a_{end-1},..., a_k},
*   which computes the max of the last k residues in the block.

*    Then m_j = max{block_maxima[j], prev_block_maxima[j-(window-1)]}
*/ 

/*  Note that for us a_i = windowscore[i] do not exist for                */
/*  i=seqlen - window+1 to seqlen-1 (i.e. number_windows=seqlen-window+1. */
/*  
/*  For these values of i, m_i=prev_block_maxima[i-window].         */




void max_over_windows (double maxima[],
		  double windowscore[],
		  int seqlen,
		  int window)
{
 int number_blocks= (int)seqlen/window;   /* Truncates */
 int number_windows=seqlen-window+1;
 int i,j;
 double prev_block_maxima[MAXSEQLEN], block_maxima[MAXSEQLEN];
 int start; 

  
  /*  First compute prev_block_max[k] defined above.  */ 

/*  Dealing with position j in block i from right to left. In the         */
/*  last block (last complete block) need to check if windowscore has any */
/*  meaning (i.e. if are further than window from sequence end).          */
/*  So start at position seqlen-window which is always in the last        */
/*  complete block in window position (seqlen mod window). */


 for (i=number_blocks-1, start=seqlen%window; i>=0; i--) {
                                /* First block, start at seqlen-window. */
   prev_block_maxima[i*window + start]=       /* Initialize block for   */
            windowscore[i*window+ start];     /* the window-1 position. */
   for (j=start-1; j>=0; j--)   /* Iterate through block passing values left.*/
     prev_block_maxima[i*window +j]= 
       MAX(windowscore[i*window +j],prev_block_maxima[i*window+j+1]);

   start= window-1;          /* For all but last block start at end of block*/
 }
 

 /* Now compute block_maxima[i] defined above with a_i=windowscore[i]. */
 for (i=0; i< number_blocks; i++)  {
   
   block_maxima[i*window] = windowscore[i*window];  /* Initialize block for */
                                                    /* the first position.  */
   for (j=1; j<window; j++){   /*Iterate through block passing values right.*/
     if ( (i*window+j +window-1) < seqlen) /* If there is a complete window */
                                           /* starting at the current pos.  */
       block_maxima[i*window +j] =         /* so winscore is defined...     */
            MAX(windowscore[i*window +j],block_maxima[i*window+j-1]) ;
     else
       block_maxima[i*window +j]= block_maxima[i*window+j-1];

     if (i>0)  /* so there is a previous window */
       maxima[i*window +j]= MAX(block_maxima[i*window+j],
			      prev_block_maxima[(i-1)*window+j+1]);
     else      /* In the first window. */
       maxima[i*window +j]= block_maxima[i*window+j];
   }    
 }

   
 
  /* Now calculate the scores for the last partial window.  */
   for (i= window*number_blocks; i<seqlen; i++)  
     maxima[i]=prev_block_maxima[i-window+1];

  /* Set maxima[window*i] which didn't get set before.  */
   maxima[0]=windowscore[0];
   for (i=1; i<number_blocks; i++)
     maxima[i*window]=MAX(prev_block_maxima[(i-1)*window +1],
			  windowscore[i*window]);

}

/*********************************************************************/
void min_over_windows (double minima[],
                  double windowscore[],
                  int seqlen,
                  int window)
{
 int number_blocks= (int)seqlen/window;   /* Truncates */
 int number_windows=seqlen-window+1;
 int i,j;
 double prev_block_minima[MAXSEQLEN], block_minima[MAXSEQLEN];
 int start; 

  
  /*  First compute prev_block_min[k] defined above.  */ 

/*  Dealing with position j in block i from right to left. In the         */
/*  last block (last complete block) need to check if windowscore has any */
/*  meaning (i.e. if are further than window from sequence end).          */
/*  So start at position seqlen-window which is always in the last        */
/*  complete block in window position (seqlen mod window). */


 for (i=number_blocks-1, start=seqlen%window; i>=0; i--) {
                                /* First block, start at seqlen-window. */
   prev_block_minima[i*window + start]=       /* Initialize block for   */
            windowscore[i*window+ start];     /* the window-1 position. */
   for (j=start-1; j>=0; j--)   /* Iterate through block passing values left.*/
     prev_block_minima[i*window +j]= 
       MIN(windowscore[i*window +j],prev_block_minima[i*window+j+1]);

   start= window-1;          /* For all but last block start at end of block*/
 }
 

 /* Now compute block_minima[i] defined above with a_i=windowscore[i]. */
 for (i=0; i< number_blocks; i++)  {
   
   block_minima[i*window] = windowscore[i*window];  /* Initialize block for */
                                                    /* the first position.  */
   for (j=1; j<window; j++){   /*Iterate through block passing values right.*/
     if ( (i*window+j +window-1) < seqlen) /* If there is a complete window */
                                           /* starting at the current pos.  */
       block_minima[i*window +j] =         /* so winscore is defined...     */
            MIN(windowscore[i*window +j],block_minima[i*window+j-1]) ;
     else
       block_minima[i*window +j]= block_minima[i*window+j-1];

     if (i>0)  /* so there is a previous window */
       minima[i*window +j]= MIN(block_minima[i*window+j],
                              prev_block_minima[(i-1)*window+j+1]);
     else      /* In the first window. */
       minima[i*window +j]= block_minima[i*window+j];
   }    
 }

   
  
  /* Now calculate the scores for the last partial window.  */
   for (i= window*number_blocks; i<seqlen; i++)  
     minima[i]=prev_block_minima[i-window+1];

  /* Set minima[window*i] which didn't get set before.  */
   minima[0]=windowscore[0];
   for (i=1; i<number_blocks; i++)
     minima[i*window]=MIN(prev_block_minima[(i-1)*window +1],
                          windowscore[i*window]);

}






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

void avg_over_windows (double average[],
		  double windowscore[],
		  int seqlen,
		  int window)
{


  int i;
  int number_windows=0;
  double total=0;


  for(i=0; i<seqlen; i++) {

    if ( i-window >= 0)  {
      total -= exp(windowscore[i-window]);
      number_windows--;
    }

    if ( i< seqlen - window+1) {    /** There are no window scores within **/
      total += exp(windowscore[i]); /* window residues of the end. **/
      number_windows++;
    }

    average[i] = getdlog(total/number_windows);
  }
}