src/secdp.c

   1 #include <stdio.h>
   2 #include "blast.h"
   3
   4 /*      $Id: secdp.c,v 1.2 2006/02/20 08:35:37 akinjo Exp $      */
   5
   6 #ifndef lint
   7 static char vcid[] = "$Id: secdp.c,v 1.2 2006/02/20 08:35:37 akinjo Exp $";
   8 #endif /* lint */
   9
  10 static const double scagap = 1.0;
  11 static  char *sec = "HEC";
  12
  13 static double gap[3][3] =
  14 /* { */
  15 /*   {       1.0113036 ,      -5.8465444 ,      -1.6370586 , }, */
  16 /*   {      -4.2604538 ,       1.3119079 ,      -0.9302283 , }, */
  17 /*   {      -1.6658889 ,      -0.9093533 ,       0.6393491 , }, */
  18 /* }; */
  19 {
  20   {      -1.1949517 ,      -8.4615910 ,      -3.5438586 , },
  21   {      -6.8755005 ,      -1.7119301 ,      -3.2458195 , },
  22   {      -3.5726889 ,      -3.2249445 ,      -0.9679954 , },
  23 };
  24
  25
  26
  27 double mymax3(double st[], int *ind) {
  28   if(st[2] >= st[1] && st[2] >= st[0]) {
  29     *ind = 2;
  30     return st[2];
  31   }
  32   else if(st[1] >= st[2] && st[1] >= st[0]) {
  33     *ind = 1;
  34     return st[1];
  35   }
  36
  37   *ind = 0;
  38   return st[0];
  39 }
  40
  41 void secdp(int len, double y[][3], char psec[])
  42 {
  43   double t[MAXSEQ][3];
  44   int bt[MAXSEQ][3];
  45   double st[3];
  46
  47
  48   int i,j,ind;
  49   double tmax;
  50
  51   for(i=0; i < len; i++) {
  52     for(j=0; j<3; j++) {
  53       t[i][j] = 0.0;
  54       bt[i][j] = 0;
  55     }
  56   }
  57
  58   for(j=0;j<3;j++) t[0][j] = y[0][j];
  59
  60   for(i = 1; i<len; i++) {
  61     for(j=0; j<3; j++) {
  62       st[0] = t[i-1][0] + gap[0][j]*scagap;
  63       st[1] = t[i-1][1] + gap[1][j]*scagap;
  64       st[2] = t[i-1][2] + gap[2][j]*scagap;
  65       tmax = mymax3(st, &ind);
  66       t[i][j] = tmax + y[i][j];
  67       bt[i][j] = ind;
  68     }
  69   }
  70
  71 /*   for(i=0; i< len; i++) { */
  72 /*     printf ("%5d : ", i+1); */
  73 /*     for(j=0; j<3; j++) { */
  74 /*       printf("%8.3f %5d %8.3f: ", t[i][j],bt[i][j], y[i][j]); */
  75 /*     } */
  76 /*     printf("\n"); */
  77 /*   } */
  78
  79   tmax = mymax3(t[len-1], &ind);
  80   psec[len-1] = sec[ind];
  81   for(i=len-1; i> 0; i--) {
  82     psec[i-1] = sec[bt[i][ind]];
  83     ind = bt[i][ind];
  84   }
  85
  86   return;
  87 }
  88
  89 void secmax(int len, double y[][3], char psec[])
  90 {
  91   int i,ind;
  92   double tmax;
  93
  94   for(i=0; i< len; i++) {
  95     tmax = mymax3(y[i],&ind);
  96     psec[i] = sec[ind];
  97   }
  98   return;
  99 }