7 /* Compute scores for residues to be the max window score of windows containing
8 * that residue. Windows are defined by where they start.
9 * windowscore[0..seqlen-window] are the window scores so window[i] is
10 * the score for the window starting in position i.
11 * The answer is placed in maxima[0..seqlen-1]
13 /* Note that the last windowscore is windowscore[seqlen-1-(window-1)]. */
17 /* A rough idea of the dynamic programming algorithm presented here:
19 * Suppose have seq of scores (a_1,a_2,....,a_n)
20 * and you want to get values (m_1,m_2,....,m_n) from the window scores
21 * m_i= max{a_{i-(window-1)}, a_{i-window+1},...., a_i}.
23 * Then divide the seq into blocks:
24 * ( [a_1,....,a_window], [a_{window+1},...,a_{2window}],...,[....,a_n] )
25 * and compute block_maxima[i]= max{a_start, a_{start+1},...., a_i},
26 * which is the max of the first i scores in the block.
28 * Also for k in a block ending at end define
29 * prev_block_maxima[k]= max{a_end, a_{end-1},..., a_k},
30 * which computes the max of the last k residues in the block.
32 * Then m_j = max{block_maxima[j], prev_block_maxima[j-(window-1)]}
35 /* Note that for us a_i = windowscore[i] do not exist for */
36 /* i=seqlen - window+1 to seqlen-1 (i.e. number_windows=seqlen-window+1. */
38 /* For these values of i, m_i=prev_block_maxima[i-window]. */
43 void max_over_windows (double maxima[],
48 int number_blocks= (int)seqlen/window; /* Truncates */
49 int number_windows=seqlen-window+1;
51 double prev_block_maxima[MAXSEQLEN], block_maxima[MAXSEQLEN];
55 /* First compute prev_block_max[k] defined above. */
57 /* Dealing with position j in block i from right to left. In the */
58 /* last block (last complete block) need to check if windowscore has any */
59 /* meaning (i.e. if are further than window from sequence end). */
60 /* So start at position seqlen-window which is always in the last */
61 /* complete block in window position (seqlen mod window). */
64 for (i=number_blocks-1, start=seqlen%window; i>=0; i--) {
65 /* First block, start at seqlen-window. */
66 prev_block_maxima[i*window + start]= /* Initialize block for */
67 windowscore[i*window+ start]; /* the window-1 position. */
68 for (j=start-1; j>=0; j--) /* Iterate through block passing values left.*/
69 prev_block_maxima[i*window +j]=
70 MAX(windowscore[i*window +j],prev_block_maxima[i*window+j+1]);
72 start= window-1; /* For all but last block start at end of block*/
76 /* Now compute block_maxima[i] defined above with a_i=windowscore[i]. */
77 for (i=0; i< number_blocks; i++) {
79 block_maxima[i*window] = windowscore[i*window]; /* Initialize block for */
80 /* the first position. */
81 for (j=1; j<window; j++){ /*Iterate through block passing values right.*/
82 if ( (i*window+j +window-1) < seqlen) /* If there is a complete window */
83 /* starting at the current pos. */
84 block_maxima[i*window +j] = /* so winscore is defined... */
85 MAX(windowscore[i*window +j],block_maxima[i*window+j-1]) ;
87 block_maxima[i*window +j]= block_maxima[i*window+j-1];
89 if (i>0) /* so there is a previous window */
90 maxima[i*window +j]= MAX(block_maxima[i*window+j],
91 prev_block_maxima[(i-1)*window+j+1]);
92 else /* In the first window. */
93 maxima[i*window +j]= block_maxima[i*window+j];
99 /* Now calculate the scores for the last partial window. */
100 for (i= window*number_blocks; i<seqlen; i++)
101 maxima[i]=prev_block_maxima[i-window+1];
103 /* Set maxima[window*i] which didn't get set before. */
104 maxima[0]=windowscore[0];
105 for (i=1; i<number_blocks; i++)
106 maxima[i*window]=MAX(prev_block_maxima[(i-1)*window +1],
107 windowscore[i*window]);
111 /*********************************************************************/
112 void min_over_windows (double minima[],
113 double windowscore[],
117 int number_blocks= (int)seqlen/window; /* Truncates */
118 int number_windows=seqlen-window+1;
120 double prev_block_minima[MAXSEQLEN], block_minima[MAXSEQLEN];
124 /* First compute prev_block_min[k] defined above. */
126 /* Dealing with position j in block i from right to left. In the */
127 /* last block (last complete block) need to check if windowscore has any */
128 /* meaning (i.e. if are further than window from sequence end). */
129 /* So start at position seqlen-window which is always in the last */
130 /* complete block in window position (seqlen mod window). */
133 for (i=number_blocks-1, start=seqlen%window; i>=0; i--) {
134 /* First block, start at seqlen-window. */
135 prev_block_minima[i*window + start]= /* Initialize block for */
136 windowscore[i*window+ start]; /* the window-1 position. */
137 for (j=start-1; j>=0; j--) /* Iterate through block passing values left.*/
138 prev_block_minima[i*window +j]=
139 MIN(windowscore[i*window +j],prev_block_minima[i*window+j+1]);
141 start= window-1; /* For all but last block start at end of block*/
145 /* Now compute block_minima[i] defined above with a_i=windowscore[i]. */
146 for (i=0; i< number_blocks; i++) {
148 block_minima[i*window] = windowscore[i*window]; /* Initialize block for */
149 /* the first position. */
150 for (j=1; j<window; j++){ /*Iterate through block passing values right.*/
151 if ( (i*window+j +window-1) < seqlen) /* If there is a complete window */
152 /* starting at the current pos. */
153 block_minima[i*window +j] = /* so winscore is defined... */
154 MIN(windowscore[i*window +j],block_minima[i*window+j-1]) ;
156 block_minima[i*window +j]= block_minima[i*window+j-1];
158 if (i>0) /* so there is a previous window */
159 minima[i*window +j]= MIN(block_minima[i*window+j],
160 prev_block_minima[(i-1)*window+j+1]);
161 else /* In the first window. */
162 minima[i*window +j]= block_minima[i*window+j];
168 /* Now calculate the scores for the last partial window. */
169 for (i= window*number_blocks; i<seqlen; i++)
170 minima[i]=prev_block_minima[i-window+1];
172 /* Set minima[window*i] which didn't get set before. */
173 minima[0]=windowscore[0];
174 for (i=1; i<number_blocks; i++)
175 minima[i*window]=MIN(prev_block_minima[(i-1)*window +1],
176 windowscore[i*window]);
185 /*********************************************************************/
187 void avg_over_windows (double average[],
188 double windowscore[],
195 int number_windows=0;
199 for(i=0; i<seqlen; i++) {
201 if ( i-window >= 0) {
202 total -= exp(windowscore[i-window]);
206 if ( i< seqlen - window+1) { /** There are no window scores within **/
207 total += exp(windowscore[i]); /* window residues of the end. **/
211 average[i] = getdlog(total/number_windows);