binaries/src/disembl/disembl.c

   1 /********************************************************/
   2 /*  ____  _     _____ __  __ ____  _       _   _____    */
   3 /* |  _ \(_)___| ____|  \/  | __ )| |     / | |___ /    */
   4 /* | | | | / __|  _| | |\/| |  _ \| |     | |   |_ \    */
   5 /* | |_| | \__ \ |___| |  | | |_) | |___  | |_ ___) |   */
   6 /* |____/|_|___/_____|_|  |_|____/|_____| |_(_)____/    */
   7 /* DisEMBL is Copyright (C) 2003                        */
   8 /* Lars Juhl Jensen & Rune Linding - EMBL               */
   9 /* Licensed under GPL                                   */
  10 /********************************************************/
  11
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <string.h>
  15
  16 /* Define size of the alphabet */
  17 #define na 21
  18
  19 /* Define max number of neurons */
  20 #define mw 41
  21 #define mh 30
  22
  23 #include "russel.h"
  24 #include "bfactor.h"
  25 #include "missing.h"
  26
  27 float sigmoid[256] = {
  28   0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000,
  29   0.000000, 0.000000, 0.000000, 0.000000, 0.000001, 0.000001, 0.000001, 0.000001,
  30   0.000001, 0.000001, 0.000001, 0.000001, 0.000001, 0.000002, 0.000002, 0.000002,
  31   0.000002, 0.000003, 0.000003, 0.000003, 0.000004, 0.000004, 0.000005, 0.000005,
  32   0.000006, 0.000007, 0.000008, 0.000009, 0.000010, 0.000011, 0.000013, 0.000015,
  33   0.000017, 0.000019, 0.000021, 0.000024, 0.000028, 0.000031, 0.000035, 0.000040,
  34   0.000045, 0.000051, 0.000058, 0.000066, 0.000075, 0.000085, 0.000096, 0.000109,
  35   0.000123, 0.000140, 0.000158, 0.000180, 0.000203, 0.000231, 0.000261, 0.000296,
  36   0.000335, 0.000380, 0.000431, 0.000488, 0.000553, 0.000626, 0.000710, 0.000804,
  37   0.000911, 0.001032, 0.001170, 0.001325, 0.001501, 0.001701, 0.001927, 0.002183,
  38   0.002473, 0.002801, 0.003173, 0.003594, 0.004070, 0.004610, 0.005220, 0.005911,
  39   0.006693, 0.007577, 0.008577, 0.009708, 0.010987, 0.012432, 0.014064, 0.015906,
  40   0.017986, 0.020332, 0.022977, 0.025957, 0.029312, 0.033086, 0.037327, 0.042088,
  41   0.047426, 0.053403, 0.060087, 0.067547, 0.075858, 0.085099, 0.095349, 0.106691,
  42   0.119203, 0.132964, 0.148047, 0.164516, 0.182426, 0.201813, 0.222700, 0.245085,
  43   0.268941, 0.294215, 0.320821, 0.348645, 0.377541, 0.407333, 0.437823, 0.468791,
  44   0.500000, 0.531209, 0.562177, 0.592667, 0.622459, 0.651355, 0.679179, 0.705785,
  45   0.731059, 0.754915, 0.777300, 0.798187, 0.817574, 0.835484, 0.851953, 0.867036,
  46   0.880797, 0.893309, 0.904651, 0.914901, 0.924142, 0.932453, 0.939913, 0.946597,
  47   0.952574, 0.957912, 0.962673, 0.966914, 0.970688, 0.974043, 0.977023, 0.979668,
  48   0.982014, 0.984094, 0.985936, 0.987568, 0.989013, 0.990292, 0.991423, 0.992423,
  49   0.993307, 0.994089, 0.994780, 0.995390, 0.995930, 0.996406, 0.996827, 0.997199,
  50   0.997527, 0.997817, 0.998073, 0.998299, 0.998499, 0.998675, 0.998830, 0.998968,
  51   0.999089, 0.999196, 0.999290, 0.999374, 0.999447, 0.999512, 0.999569, 0.999620,
  52   0.999665, 0.999704, 0.999739, 0.999769, 0.999797, 0.999820, 0.999842, 0.999860,
  53   0.999877, 0.999891, 0.999904, 0.999915, 0.999925, 0.999934, 0.999942, 0.999949,
  54   0.999955, 0.999960, 0.999965, 0.999969, 0.999972, 0.999976, 0.999979, 0.999981,
  55   0.999983, 0.999985, 0.999987, 0.999989, 0.999990, 0.999991, 0.999992, 0.999993,
  56   0.999994, 0.999995, 0.999995, 0.999996, 0.999996, 0.999997, 0.999997, 0.999997,
  57   0.999998, 0.999998, 0.999998, 0.999998, 0.999999, 0.999999, 0.999999, 0.999999,
  58   0.999999, 0.999999, 0.999999, 0.999999, 0.999999, 1.000000, 1.000000, 1.000000,
  59   1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000
  60 };
  61
  62
  63
  64
  65 /*
  66  *  Perform feed forward evaluation of a neural network on a sequence window.
  67  */
  68
  69 float feed_forward(int const *s, float const w[], int nw, int nh) {
  70
  71   float h[mh], o[2], x;
  72   int i, j;
  73
  74   /* Shift input window to match network window size */
  75   s += (mw-nw)/2;
  76
  77   /* Feed input values to hidden layer making use of sparse encoding */
  78   for (i = 0; i < nh; ++i) {
  79     x = w[(na*nw+1)*(i+1)-1];
  80     for (j = 0; j < nw; ++j) {
  81       x += w[(na*nw+1)*i+na*j+s[j]];
  82     }
  83     if (x <= -16) {
  84       h[i] = 0;
  85     } else if (x >= 16) {
  86       h[i] = 1;
  87     } else {
  88       h[i] = sigmoid[(int)(8*x+128)];
  89     }
  90   }
  91
  92   /* Feed hidden layer values to output layer */
  93   for (i = 0; i <= 1; ++i) {
  94     x = w[(na*nw+1)*nh+(nh+1)*(i+1)-1];
  95     for (j = 0; j < nh; ++j) {
  96       x += w[(na*nw+1)*nh+(nh+1)*i+j]*h[j];
  97     }
  98     if (x <= -16) {
  99       o[i] = 0;
 100     } else if (x >= 16) {
 101       o[i] = 1;
 102     } else {
 103       o[i] = sigmoid[(int)(8*x+128)];
 104     }
 105   }
 106
 107   /* Combine the scores from the two output neurons */
 108   return((o[0]+1-o[1])/2);
 109
 110 }
 111
 112
 113
 114
 115 /*
 116  *  Calculate and print scores for a sequence window.
 117  */
 118
 119 void predict(int const *s) {
 120
 121   float sm, sb, sr;
 122
 123   if (s[(mw-1)/2] == na-1) {
 124     return;
 125   }
 126
 127   sr = 0;
 128   sr += feed_forward(s, r19_1, 19, 30);
 129   sr += feed_forward(s, r19_2, 19, 30);
 130   sr += feed_forward(s, r19_3, 19, 30);
 131   sr += feed_forward(s, r19_4, 19, 30);
 132   sr += feed_forward(s, r19_5, 19, 30);
 133   sr /= 5;
 134   sr = 0.07387214+0.8020778*sr;
 135
 136   sb = 0;
 137   sb += feed_forward(s, b41_1, 41, 5);
 138   sb += feed_forward(s, b41_2, 41, 5);
 139   sb += feed_forward(s, b41_3, 41, 5);
 140   sb += feed_forward(s, b41_4, 41, 5);
 141   sb += feed_forward(s, b41_5, 41, 5);
 142   sb /= 5;
 143   sb = 0.08016882+0.6282424*sb;
 144
 145   sm = 0;
 146   sm += feed_forward(s, m9_1, 9, 30);
 147   sm += feed_forward(s, m9_2, 9, 30);
 148   sm += feed_forward(s, m9_3, 9, 30);
 149   sm += feed_forward(s, m9_4, 9, 30);
 150   sm += feed_forward(s, m9_5, 9, 30);
 151   sm += feed_forward(s, m21_1, 21, 30);
 152   sm += feed_forward(s, m21_2, 21, 30);
 153   sm += feed_forward(s, m21_3, 21, 30);
 154   sm += feed_forward(s, m21_4, 21, 30);
 155   sm += feed_forward(s, m21_5, 21, 30);
 156   sm /= 10;
 157
 158   printf("%f\t%f\t%f\n", sr, sr*sb, sm);
 159
 160 }
 161
 162 int main(int ARGC, char *ARGV[]) {
 163
 164   char *alphabet = strdup("FIVWMLCHYAGNRTPDEQSK");
 165
 166   char *p;
 167   int c, i, j, s[mw];
 168
 169   for (i = 0; i < mw; ++i) {
 170     s[i] = na-1;
 171   }
 172
 173   while (!feof(stdin)) {
 174
 175     c = fgetc(stdin);
 176     p = strchr(alphabet, c);
 177     if (p != NULL) {
 178       for (i = 1; i < mw; ++i) {
 179         s[mw-i] = s[mw-i-1];
 180       }
 181       s[0] = p-alphabet;
 182       predict(s);
 183     } else if ((char)c == '\n') {
 184       for (i = 0; i < (mw-1)/2; i++) {
 185         for (j = 1; j < mw; ++j) {
 186           s[mw-j] = s[mw-j-1];
 187         }
 188         s[0] = na-1;
 189         predict(s);
 190       }
 191     }
 192
 193   }
 194
 195   return(0);
 196
 197 }