sources/multicoil/dyn.c

   1 /*  Ethan Wolf 1995 */
   2
   3 #include <stdio.h>
   4 #include "scconst.h"
   5 #include "math.h"
   6
   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]
  12  */
  13 /* Note that the last windowscore is windowscore[seqlen-1-(window-1)]. */
  14
  15
  16
  17 /*  A rough idea of the dynamic programming algorithm presented here:
  18 *
  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}.
  22 *
  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.
  27
  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.
  31
  32 *    Then m_j = max{block_maxima[j], prev_block_maxima[j-(window-1)]}
  33 */
  34
  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. */
  37 /*
  38 /*  For these values of i, m_i=prev_block_maxima[i-window].         */
  39
  40
  41
  42
  43 void max_over_windows (double maxima[],
  44                   double windowscore[],
  45                   int seqlen,
  46                   int window)
  47 {
  48  int number_blocks= (int)seqlen/window;   /* Truncates */
  49  int number_windows=seqlen-window+1;
  50  int i,j;
  51  double prev_block_maxima[MAXSEQLEN], block_maxima[MAXSEQLEN];
  52  int start;
  53
  54
  55   /*  First compute prev_block_max[k] defined above.  */
  56
  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). */
  62
  63
  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]);
  71
  72    start= window-1;          /* For all but last block start at end of block*/
  73  }
  74
  75
  76  /* Now compute block_maxima[i] defined above with a_i=windowscore[i]. */
  77  for (i=0; i< number_blocks; i++)  {
  78
  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]) ;
  86      else
  87        block_maxima[i*window +j]= block_maxima[i*window+j-1];
  88
  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];
  94    }
  95  }
  96
  97
  98
  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];
 102
 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]);
 108
 109 }
 110
 111 /*********************************************************************/
 112 void min_over_windows (double minima[],
 113                   double windowscore[],
 114                   int seqlen,
 115                   int window)
 116 {
 117  int number_blocks= (int)seqlen/window;   /* Truncates */
 118  int number_windows=seqlen-window+1;
 119  int i,j;
 120  double prev_block_minima[MAXSEQLEN], block_minima[MAXSEQLEN];
 121  int start;
 122
 123
 124   /*  First compute prev_block_min[k] defined above.  */
 125
 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). */
 131
 132
 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]);
 140
 141    start= window-1;          /* For all but last block start at end of block*/
 142  }
 143
 144
 145  /* Now compute block_minima[i] defined above with a_i=windowscore[i]. */
 146  for (i=0; i< number_blocks; i++)  {
 147
 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]) ;
 155      else
 156        block_minima[i*window +j]= block_minima[i*window+j-1];
 157
 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];
 163    }
 164  }
 165
 166
 167
 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];
 171
 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]);
 177
 178 }
 179
 180
 181
 182
 183
 184
 185 /*********************************************************************/
 186
 187 void avg_over_windows (double average[],
 188                   double windowscore[],
 189                   int seqlen,
 190                   int window)
 191 {
 192
 193
 194   int i;
 195   int number_windows=0;
 196   double total=0;
 197
 198
 199   for(i=0; i<seqlen; i++) {
 200
 201     if ( i-window >= 0)  {
 202       total -= exp(windowscore[i-window]);
 203       number_windows--;
 204     }
 205
 206     if ( i< seqlen - window+1) {    /** There are no window scores within **/
 207       total += exp(windowscore[i]); /* window residues of the end. **/
 208       number_windows++;
 209     }
 210
 211     average[i] = getdlog(total/number_windows);
 212   }
 213 }
 214