--- /dev/null
+/* 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);
+ }
+}
+
git://source.jalview.org / jpred.git / blobdiff