binaries/src/ViennaRNA/lib/svm_utils.c

   1 #include <config.h>
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <math.h>
   5 #include <ctype.h>
   6 #include <string.h>
   7
   8 #include "utils.h"
   9 #include "fold_vars.h"
  10 #include "pair_mat.h"
  11 #include "svm.h"
  12 #include "svm_utils.h"
  13
  14 #include "model_avg.inc"  /* defines avg_model_string */
  15 #include "model_sd.inc"   /* defines sd_model_string */
  16
  17
  18 PRIVATE svm_model *avg_model;
  19 PRIVATE svm_model *sd_model;
  20
  21 PRIVATE void    freeFields(char** fields);
  22 PRIVATE char**  splitFields(char* string);
  23 PRIVATE char**  splitLines(char* string);
  24
  25 PUBLIC float get_z(char *sequence, double energy) {
  26   double average_free_energy;
  27   double sd_free_energy;
  28   float my_z;
  29   int info_avg;
  30   make_pair_matrix();
  31   short *S      = encode_sequence(sequence, 0);
  32   int   length  = S[0];
  33   int   *AUGC   = get_seq_composition(S, 1, length);
  34   avg_model     = svm_load_model_string(avg_model_string);
  35   sd_model      = svm_load_model_string(sd_model_string);
  36   average_free_energy = avg_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4], avg_model, &info_avg);
  37
  38   if(info_avg == 0){
  39     double difference = (energy/* /100*/) - average_free_energy;
  40     sd_free_energy    = sd_regression(AUGC[0], AUGC[1], AUGC[2], AUGC[3], AUGC[4], sd_model);
  41     my_z              = difference / sd_free_energy;
  42   }
  43   else{
  44     fprintf(stderr,"warning: sequence out of bounds\n");
  45 #if 0
  46     my_z = shuffle_score(sequence, energy);
  47 #endif
  48   }
  49   free(AUGC);
  50   free(S);
  51   svm_destroy_model(avg_model);
  52   svm_destroy_model(sd_model);
  53   return my_z;
  54 }
  55
  56 PUBLIC int *get_seq_composition(short *S, unsigned int start, unsigned int stop){
  57   unsigned int i;
  58   int *ret = (int *)space(sizeof(int) * 6);
  59
  60   for (i=MAX2(start, 1); i <= MIN2(stop, S[0]); i++){
  61     if(S[i] > 4)  ret[0]++;
  62     else          ret[S[i]]++;
  63   }
  64   ret[5] = -1; /* indicate last entry */
  65   return ret;
  66 }
  67
  68 PUBLIC double sd_regression(int N, int A, int C, int G, int T,  svm_model *sd_model){
  69   double sd_free_energy = 0.0;
  70   int length = A + C + G + T + N;
  71   double GC_content  = (double) (G + C)/length;
  72   double AT_ratio    = (double) A/(A+T);
  73   double CG_ratio    = (double) C/(C+G);
  74   double norm_length = (double) (length-50)/350.0;
  75   struct svm_node node_mono[5];
  76
  77   node_mono[0].index = 1; node_mono[0].value = GC_content;
  78   node_mono[1].index = 2; node_mono[1].value = AT_ratio;
  79   node_mono[2].index = 3; node_mono[2].value = CG_ratio;
  80   node_mono[3].index = 4; node_mono[3].value = norm_length;
  81   node_mono[4].index =-1;
  82
  83   sd_free_energy = svm_predict(sd_model,node_mono);
  84
  85   sd_free_energy = (double) sd_free_energy * sqrt(length);
  86
  87   return sd_free_energy;
  88 }
  89
  90 PUBLIC double avg_regression(int N, int A, int C, int G, int T, struct svm_model *avg_model, int *info ){
  91   double average_free_energy = 0.0;
  92
  93   int length = A + C + G + T + N;
  94   double N_fraction = (double) N/length;
  95   double GC_content = (double) (G + C)/length;
  96   double AT_ratio   = (double) A/(A+T);
  97   double CG_ratio   = (double) C/(C+G);
  98
  99   double norm_length = (double) (length-50)/350.0;
 100
 101   struct svm_node node_mono[5];
 102   *info = 0;
 103   if ( length < 50 || length > 400 ) {
 104     *info = 1;
 105     return 0.0;
 106   }
 107   if ( N_fraction > 0.05 ) {
 108     *info = 2;
 109     return 0.0;
 110   }
 111   if ( GC_content < 0.20 || GC_content > 0.80 ) {
 112     *info = 3;
 113     return 0.0;
 114   }
 115   if ( AT_ratio < 0.20 || AT_ratio > 0.80 ) {
 116     *info = 4;
 117     return 0.0;
 118   }
 119   if ( CG_ratio < 0.20 || CG_ratio > 0.80 ) {
 120     *info = 5;
 121     return 0.0;
 122   }
 123
 124   node_mono[0].index = 1; node_mono[0].value = GC_content;
 125   node_mono[1].index = 2; node_mono[1].value = AT_ratio;
 126   node_mono[2].index = 3; node_mono[2].value = CG_ratio;
 127   node_mono[3].index = 4; node_mono[3].value = norm_length;
 128   node_mono[4].index =-1;
 129
 130   average_free_energy = svm_predict(avg_model,node_mono);
 131
 132   average_free_energy = (double) average_free_energy * length;
 133
 134   return average_free_energy;
 135 }
 136
 137 PUBLIC double minimal_sd(int N, int A, int C, int G, int T ){
 138   int length = A + C + G + T + N;
 139   if ( length <  60 ) return 0.450324;
 140   if ( length <  70 ) return 0.749771;
 141   if ( length <  80 ) return 1.029421;
 142   if ( length <  90 ) return 1.027517;
 143   if ( length <  100 ) return 1.347283;
 144   if ( length <  120 ) return 1.112086;
 145   if ( length <  150 ) return 1.574339;
 146   if ( length <  170 ) return 1.779043;
 147   if ( length <  200 ) return 1.922908;
 148   if ( length <  250 ) return 2.226856;
 149   if ( length <  300 ) return 2.349300;
 150   if ( length <  350 ) return 2.589703;
 151   if ( length <  400 ) return 2.791215;
 152
 153   return 0.450324;
 154 }
 155
 156 PUBLIC svm_model  *svm_load_model_string(char *modelString){
 157
 158   /* redefinition from svm.cpp */
 159   char *svm_type_table[]={"c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL};
 160   char *kernel_type_table[]={"linear","polynomial","rbf","sigmoid",NULL};
 161
 162   struct svm_model *model;
 163   char **lines, **fields;
 164   int i,j,k,l,m;
 165   char *key, *value, *field;
 166   char c;
 167   int dataStart, elements;
 168   int isColon;
 169   struct svm_node *x_space=NULL;
 170
 171   model = (struct svm_model*)space(sizeof(struct svm_model));
 172
 173   model->rho = NULL;
 174   model->probA = NULL;
 175   model->probB = NULL;
 176   model->label = NULL;
 177   model->nSV = NULL;
 178
 179
 180   /* Read header until support vectors start */
 181   lines=splitLines(modelString);
 182   i=0;
 183   while (strcmp(lines[i],"SV")!=0){
 184         fields=splitFields(lines[i]);
 185
 186         key=fields[0];
 187
 188         if(strcmp(key,"svm_type")==0){
 189           value=fields[1];
 190           for(j=0;svm_type_table[j];j++){
 191                 if(strcmp(svm_type_table[j],value)==0){
 192                   model->param.svm_type=j;
 193                   break;
 194                 }
 195           }
 196           if(svm_type_table[i] == NULL){
 197                 fprintf(stderr,"unknown svm type.\n");
 198                 free(model->rho);
 199                 free(model->label);
 200                 free(model->nSV);
 201                 free(model);
 202                 return NULL;
 203           }
 204         } else
 205
 206         if(strcmp(key,"kernel_type")==0){
 207           value=fields[1];
 208           for(j=0;kernel_type_table[j];j++){
 209                 if(strcmp(kernel_type_table[j],value)==0){
 210                   model->param.kernel_type=j;
 211                   break;
 212                 }
 213           }
 214           if(kernel_type_table[i] == NULL){
 215                 fprintf(stderr,"unknown kernel type.\n");
 216                 free(model->rho);
 217                 free(model->label);
 218                 free(model->nSV);
 219                 free(model);
 220                 return NULL;
 221           }
 222         } else
 223
 224         if (strcmp(key,"gamma")==0){
 225           value=fields[1];
 226           sscanf(value,"%lf",&model->param.gamma);
 227         }
 228
 229         if (strcmp(key,"degree")==0){
 230           value=fields[1];
 231           sscanf(value,"%d",&model->param.degree);
 232         } else
 233
 234         if (strcmp(key,"coef0")==0){
 235           value=fields[1];
 236           sscanf(value,"%lf",&model->param.coef0);
 237         } else
 238         if (strcmp(key,"nr_class")==0){
 239           value=fields[1];
 240           sscanf(value,"%d",&model->nr_class);
 241         } else
 242         if (strcmp(key,"total_sv")==0){
 243           value=fields[1];
 244           sscanf(value,"%d",&model->l);
 245         } else
 246
 247         if (strcmp(key,"rho")==0){
 248           int n = model->nr_class * (model->nr_class-1)/2;
 249           model->rho = (double*)space(sizeof(double)*n);
 250           for(j=0;j<n;j++){
 251                 sscanf(fields[j+1],"%lf",&model->rho[j]);
 252           }
 253         } else
 254
 255         if (strcmp(key,"nr_sv")==0){
 256           int n = model->nr_class;
 257           model->nSV = (int*)space(sizeof(int)*n);
 258           for(j=0;j<n;j++){
 259                 sscanf(fields[j+1],"%d",&model->nSV[j]);
 260           }
 261         } else
 262
 263         if (strcmp(key,"label")==0){
 264           int n = model->nr_class;
 265           model->label = (int*)space(sizeof(int)*n);
 266           for(j=0;j<n;j++){
 267                 sscanf(fields[j+1],"%d",&model->label[j]);
 268           }
 269         } else
 270
 271         if (strcmp(key,"probA")==0){
 272           int n = model->nr_class * (model->nr_class-1)/2;
 273           model->probA = (double*)space(sizeof(double)*n);
 274           for(j=0;j<n;j++){
 275                 sscanf(fields[j+1],"%lf",&model->probA[j]);
 276           }
 277         } else
 278
 279         if (strcmp(key,"probB")==0){
 280           int n = model->nr_class * (model->nr_class-1)/2;
 281           model->probB = (double*)space(sizeof(double)*n);
 282           for(j=0;j<n;j++){
 283                 sscanf(fields[j+1],"%lf",&model->probB[j]);
 284           }
 285         }
 286         i++;
 287         freeFields(fields);
 288   }
 289
 290   dataStart=i+1;
 291   elements=0;
 292
 293   /* Count number of nodes (by counting colons) in advance to allocate
 294          memory in one block */
 295   while (lines[i]!=NULL){
 296         j=0;
 297         while ((c=lines[i][j])!='\0'){
 298           if (c==':'){
 299                 elements++;
 300           }
 301           j++;
 302         }
 303         elements++;
 304         i++;
 305   }
 306
 307   /* allocate memory for SVs and coefficients */
 308   m = model->nr_class - 1;
 309   l = model->l;
 310   model->sv_coef = (double**)space(sizeof(double*)*m);
 311   for(i=0;i<m;i++){
 312         model->sv_coef[i] = (double*)space(sizeof(double)*l);
 313   }
 314   model->SV = (struct svm_node**)space(sizeof(struct svm_node*)*l);
 315
 316   if(l>0){
 317     x_space = (struct svm_node*)space(sizeof(struct svm_node)*(elements));
 318   }
 319
 320
 321   /* parse support vector data */
 322   j=0;
 323   for(i=0;i<l;i++){
 324         fields=splitFields(lines[dataStart+i]);
 325         model->SV[i] = &x_space[j];
 326         k=0;
 327         while ((field=fields[k])!=NULL){
 328           if (k<m){
 329             sscanf(fields[k],"%lf",&model->sv_coef[k][i]);
 330           } else {
 331             sscanf(fields[k],"%d:%lf",&(x_space[j].index),&(x_space[j].value));
 332             j++;
 333           }
 334           k++;
 335         }
 336         x_space[j++].index = -1;
 337         freeFields(fields);
 338   }
 339   freeFields(lines);
 340
 341   model->free_sv = 1;
 342
 343   return(model);
 344 }
 345
 346 PRIVATE char **splitFields(char* string){
 347
 348   char c;
 349   char* currField;
 350   char** output=NULL;
 351   int* seps;
 352   int nSep;
 353   int nField=0;
 354   int i=0;
 355
 356   if (strlen(string)==0 || string==NULL){
 357         return NULL;
 358   }
 359
 360   /* First find all characters which are whitespaces and store the
 361          positions in the array seps */
 362
 363   seps=(int *)space(sizeof(int));
 364   seps[0]=-1;
 365   nSep=1;
 366
 367   while ((c=string[i])!='\0' && (c!='\n')){
 368         if (isspace(c)){
 369           seps=(int*)xrealloc(seps,sizeof(int)*(nSep+1));
 370           seps[nSep++]=i;
 371         }
 372         i++;
 373   }
 374
 375   seps=(int*)xrealloc(seps,sizeof(int)*(nSep+1));
 376   seps[nSep]=strlen(string);
 377
 378
 379   /* Then go through all intervals in between of two whitespaces (or
 380          end or start of string) and store the fields in the array
 381          "output"; if there are two adjacent whitespaces this is ignored
 382          resulting in a behaviour like "split /\s+/" in perl */
 383
 384   for (i=0;i<nSep;i++){
 385
 386         int start=seps[i];
 387         int stop=seps[i+1];
 388         int length=(stop-start);
 389         int notSpace,j;
 390
 391
 392         currField=(char *)space(sizeof(char)*(length+1));
 393         strncpy(currField,string+start+1,length-1);
 394         currField[length]='\0';
 395
 396         /* check if field is not only whitespace */
 397         notSpace=0;
 398         j=0;
 399         while ((c=currField[j])!='\0'){
 400           if (!isspace(c)){
 401                 notSpace=1;
 402                 break;
 403           }
 404         }
 405
 406         if (notSpace){
 407           output=(char**)xrealloc(output,sizeof(char**)*(nField+1));
 408           output[nField++]=currField;
 409           currField=NULL;
 410         } else {
 411           free(currField);
 412           currField=NULL;
 413         }
 414
 415         /* printf("%s|\n",output[nField-1]); */
 416   }
 417
 418   if (nField==0){
 419         return NULL;
 420   }
 421
 422
 423   output=(char**)xrealloc(output,sizeof(char**)*(nField+1));
 424   output[nField]=NULL;
 425
 426   free(seps);
 427   return output;
 428
 429 }
 430
 431 PRIVATE char **splitLines(char* string){
 432
 433   char c;
 434   char* currLine=NULL;
 435   char** output=NULL;
 436   int i=0;
 437   int currLength=0;
 438   int lineN=0;
 439
 440   while ((c=string[i])!='\0'){
 441
 442         if (c=='\n'){
 443           output=(char**)xrealloc(output,sizeof(char**)*(lineN+1));
 444           currLine=(char*)xrealloc(currLine,sizeof(char)*(currLength+1));
 445           currLine[currLength]='\0';
 446           output[lineN]=currLine;
 447           currLength=0;
 448           currLine=NULL;
 449           lineN++;
 450         } else {
 451
 452           currLine=(char*)xrealloc(currLine,sizeof(char)*(currLength+1));
 453           currLine[currLength]=c;
 454           currLength++;
 455         }
 456         i++;
 457   }
 458
 459   output=(char**)xrealloc(output,sizeof(char**)*(lineN+1));
 460   output[lineN]=NULL;
 461
 462   return output;
 463
 464 }
 465
 466 /*  for both splitLines and splitFields */
 467 void freeFields(char** fields){
 468
 469   int i=0;
 470   while (fields[i]!=NULL){
 471         free(fields[i++]);
 472   }
 473   free(fields);
 474 }