binaries/src/ViennaRNA/lib/stringdist.c

   1 /*
   2                 String alignment for RNA secondary structures
   3                       Peter F Stadler, Ivo Hofacker
   4                            Vienna RNA Package
   5 */
   6 /* Last changed Time-stamp: <2005-07-23 10:19:40 ivo> */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include <ctype.h>
  12 #include <math.h>
  13 #include  "edit_cost.h"
  14 #include  "dist_vars.h"
  15 #include  "utils.h"
  16
  17 static char rcsid[] = "$Id: stringdist.c,v 1.3 2005/07/24 08:37:15 ivo Exp $";
  18
  19 #define PUBLIC
  20 #define PRIVATE        static
  21
  22 #define INFINITY_DIST  10000
  23
  24 PUBLIC  float      string_edit_distance(swString *T1, swString *T2);
  25 PUBLIC  swString  *Make_swString(char *string);
  26 PUBLIC  void       print_swString(swString *x);
  27 PUBLIC  void       print_alignment_list(void);
  28
  29 PRIVATE void       sprint_aligned_swStrings(swString *T1, swString *T2);
  30 PRIVATE float      StrEditCost(int i, int j, swString *T1, swString *T2);
  31 PRIVATE void       DeCode(char *string, int k, int *tp, float *w);
  32 PRIVATE int        decode(char *id);
  33 PRIVATE void       encode(int type, char label[]);
  34 PRIVATE int       *alignment[2]; /* contains information from backtracking
  35                                     alignment[0][n] is the node in tree2
  36                                     matching node n in tree1               */
  37
  38
  39 /*---------------------------------------------------------------------------*/
  40
  41 PUBLIC float string_edit_distance(swString *T1, swString *T2)
  42
  43 {
  44     float  **distance;
  45     short    **i_point, **j_point;
  46
  47     int           i, j, i1, j1, pos, length1,length2;
  48     float         minus, plus, change, temp;
  49
  50     if (cost_matrix==0) EditCost = &UsualCost;
  51     else EditCost = &ShapiroCost;
  52
  53     length1 = T1[0].sign;
  54     length2 = T2[0].sign;
  55
  56     distance = (float **)  space((length1 +1)*sizeof(float *));
  57     if(edit_backtrack){
  58        i_point  = (short **)  space((length1 +1)*sizeof(short *));
  59        j_point  = (short **)  space((length1 +1)*sizeof(short *));
  60     }
  61     for(i=0; i<= length1; i++){
  62        distance[i] = (float *) space( (length2+1)*sizeof(float));
  63        if(edit_backtrack){
  64           i_point[i]  = (short *) space( (length2+1)*sizeof(short));
  65           j_point[i]  = (short *) space( (length2+1)*sizeof(short));
  66        }
  67     }
  68
  69     for(i = 1; i <= length1; i++) {
  70        if (edit_backtrack){
  71           i_point[i][0] = i-1;
  72           j_point[i][0] = 0;
  73        }
  74        distance[i][0] = distance[i-1][0]+StrEditCost(i,0,T1,T2);
  75     }
  76     for(j = 1; j <= length2; j++) {
  77        if (edit_backtrack){
  78           j_point[0][j] = j-1;
  79           i_point[0][j] = 0;
  80        }
  81        distance[0][j] = distance[0][j-1]+StrEditCost(0,j,T1,T2);
  82     }
  83
  84     for (i = 1; i <= length1; i++) {
  85        for (j = 1; j <= length2 ; j++) {
  86           minus  = distance[i-1][j]  + StrEditCost(i,0,T1,T2);
  87           plus   = distance[i][j-1]  + StrEditCost(0,j,T1,T2);
  88           change = distance[i-1][j-1]+ StrEditCost(i,j,T1,T2);
  89
  90           distance[i][j] = MIN3(minus, plus, change);
  91           /* printf("%g ", distance[i][j]); */
  92
  93           if(edit_backtrack){
  94              if(distance[i][j] == change) {
  95                 i_point[i][j]=i-1; j_point[i][j]=j-1;  }
  96              else if(distance[i][j] == plus) {
  97                 i_point[i][j]=i  ; j_point[i][j]=j-1;  }
  98              else {
  99                 i_point[i][j]=i-1; j_point[i][j]=j  ;  }
 100           }
 101        }
 102        /* printf("\n"); */
 103     }
 104     /* printf("\n"); */
 105     temp = distance[length1][length2];
 106     for(i=0;i<=length1;i++)
 107        free(distance[i]);
 108     free(distance);
 109
 110     if(edit_backtrack){
 111        if(alignment[0]!= NULL) free(alignment[0]);
 112        if(alignment[1]!= NULL) free(alignment[1]);
 113        alignment[0] = (int *) space((length1+length2+1)*sizeof(int));
 114        alignment[1] = (int *) space((length1+length2+1)*sizeof(int));
 115
 116        pos = length1+length2;
 117        i   = length1;
 118        j   = length2;
 119        while( (i>0)||(j>0) ) {
 120           i1 = i_point[i][j];
 121           j1 = j_point[i][j];
 122           if( ((i-i1)==1)&&((j-j1)==1) )  {  /* substitution    */
 123               alignment[0][pos] = i;
 124               alignment[1][pos] = j;
 125           }
 126           if( ((i-i1)==1)&&(j==j1) )      {  /* Deletion in [1] */
 127               alignment[0][pos] = i;
 128               alignment[1][pos] = 0;
 129           }
 130           if( (i==i1)&&((j-j1)==1)  )      {  /* Deletion in [0] */
 131               alignment[0][pos] = 0;
 132               alignment[1][pos] = j;
 133           }
 134           pos--;
 135           i = i1;
 136           j = j1;
 137        }
 138        for(i=pos+1; i<=length1+length2; i++){
 139           alignment[0][i-pos] = alignment[0][i];
 140           alignment[1][i-pos] = alignment[1][i];
 141        }
 142        alignment[0][0] = length1+length2-pos;   /* length of alignment */
 143
 144        for(i=0; i<=length1; i++){
 145           free(i_point[i]); free(j_point[i]);
 146        }
 147        free(i_point); free(j_point);
 148        sprint_aligned_swStrings(T1,T2);
 149
 150     }
 151
 152     return temp;
 153 }
 154
 155
 156 /*---------------------------------------------------------------------------*/
 157
 158 PRIVATE float StrEditCost(int i, int j, swString *T1, swString *T2)
 159 {
 160     float  c, diff, cd, min, a, b, dist;
 161
 162     if(i==0) {
 163        cd   =  (float) (*EditCost)[0][T2[j].type];
 164        diff =  T2[j].weight;
 165        dist =  cd*diff;
 166     }
 167     else
 168     if(j==0) {
 169        cd   =  (float) (*EditCost)[T1[i].type][0];
 170        diff =  T1[i].weight;
 171        dist =  cd*diff;
 172     }
 173     else
 174     if( ((T1[i].sign)*(T2[j].sign)) > 0) {
 175        c = (float) (*EditCost)[T1[i].type][T2[j].type];
 176        diff = (float) fabs((a=T1[i].weight) - (b=T2[j].weight));
 177        min = MIN2(a,b);
 178        if (min == a) cd = (float) (*EditCost)[0][T2[j].type];
 179        else          cd = (float) (*EditCost)[T1[i].type][0];
 180        dist = c * min + cd * diff;
 181     }
 182     else dist = (float) INFINITY_DIST;
 183     return dist;
 184 }
 185
 186 /*---------------------------------------------------------------------------*/
 187
 188 PUBLIC swString *Make_swString(char *string)
 189 {
 190    int i=0, j=0, k=0;
 191    int tp, len, l, length;
 192    float w;
 193    swString  *S;
 194
 195    length = strlen(string);
 196
 197    for(i=0; i<length; i++) {
 198       if( (string[i]=='(') || (string[i]==')') ) j++;
 199       if(string[i]=='.') j+=2;
 200    }
 201
 202    len = j;
 203
 204    S= (swString *) space(sizeof(swString)*(len+1));
 205    S[0].sign = j; /* number of entries */
 206    S[0].weight= 0.0;
 207    S[0].type=     0;
 208
 209    i=0;
 210    j=1;
 211    while(i<length){
 212       switch(string[i]){
 213        case '(' :
 214           S[j].sign = 1;
 215           l=1;
 216           k=i;
 217           while (l>0) {
 218              k++;
 219              if(string[k] == '(' ) l++;
 220              if(string[k] == ')' ) l--;
 221           }
 222           DeCode(string,k,&tp,&w);
 223           S[j].type   = tp;
 224           S[j].weight = w/2.0;
 225           j++;
 226           break;
 227        case ')' :
 228           k=i;
 229           S[j].sign = -1;
 230           DeCode(string,k,&tp,&w);
 231           S[j].type = tp;
 232           S[j].weight = w/2.0;
 233           j++;
 234           break;
 235        case '.' :
 236           S[j].sign = 1;
 237           S[j].type = 1;
 238           S[j].weight = 0.5;
 239           j++;
 240           S[j].sign = -1;
 241           S[j].type =  1;
 242           S[j].weight = 0.5;
 243           j++;
 244           break;
 245       }
 246       i++;
 247    }
 248    return S;
 249 }
 250
 251 /*---------------------------------------------------------------------------*/
 252
 253 PRIVATE void DeCode(char *string, int k, int *tp, float *w)
 254    /* retrieves type and weigth for a node closed  by a bracket at position k */
 255 {
 256    int i,j,l,m;
 257    char  label[20], id[20] ;
 258    i=k;
 259    label[0] = '\0';
 260    while(i>=0){
 261       i--;
 262       if( (string[i]=='(')||(string[i]==')')||(string[i]=='.') ) break;
 263       else {
 264          label[k-i-1] = string[i]; label[k-i] = '\0';
 265       }
 266    }
 267    l=strlen(label);
 268    if (l==0) {           /* Dot-Bracket notation */
 269      *w  = 1.0;
 270      *tp =   2;
 271    }
 272    else{
 273      for (i=0; i<l; i++) {
 274        if (!isalpha(label[l-i-1])) break;
 275        id[i] = label[l-i-1];
 276       }
 277       id[i] = '\0';
 278       *tp=decode(id);
 279       l=l-i-1;
 280       if(l>=0){
 281          for(j=0; j<=l; j++)
 282             id[j] = label[l-j];
 283          label[l+1] ='\0';
 284          m=-1;
 285          sscanf(label,"%d",&m);
 286          *w= (float) m;
 287          if(m==-1) {
 288             fprintf(stderr, "Warning: Non-integer weight in DeCode ignored\n");
 289             *w=1.0;
 290          }
 291       }
 292       else
 293          *w=1.0;
 294    }
 295 }
 296
 297 /*---------------------------------------------------------------------------*/
 298
 299 PRIVATE int decode(char *id)
 300 {
 301     int   n, quit, i;
 302     char  label[100], *code;
 303
 304     n = 0;
 305
 306     quit = 0;
 307     code = coding;
 308
 309     while (!quit) {
 310         for (i = 0; code[i] != sep; i++) {
 311             if (code[i] == '\0') {
 312                 quit = 1;
 313                 break;
 314             }
 315             label[i] = code[i];
 316         }
 317         label[i] = '\0';
 318         if (strcmp(id, label) == 0) return (n);
 319         code += (i+1);
 320         n++;
 321     }
 322
 323     fprintf(stderr,"Syntax error: node identifier \"%s\" not found "
 324                    "in coding string \"%s\"\n", id, coding);
 325     fprintf(stderr,"Exiting...");
 326     exit(0);
 327 }
 328
 329
 330 /*---------------------------------------------------------------------------*/
 331
 332 PRIVATE void encode( int type, char label[])
 333
 334 {
 335     int   i, l;
 336
 337     l = 0;
 338     for (i = 0; i < type; i++) {
 339         while (coding[l] != sep && coding[l]) l++;
 340         l++;
 341     }
 342
 343     for (i = 0; coding[l+i] != sep; i++) {
 344         if (coding[l+i] == '\0') break;
 345         label[i] = coding[l+i];
 346     }
 347     label[i] = '\0';
 348 }
 349
 350 /*---------------------------------------------------------------------------*/
 351
 352 PRIVATE void sprint_aligned_swStrings(swString *T1, swString *T2)
 353 {
 354    int i, j, l0, l1, ltmp=0, weights;
 355    char label[10], *a0, *a1, tmp0[20], tmp1[20];
 356
 357    weights = 0;
 358    for (i=1; i<=T1[0].sign; i++) weights = (weights||(T1[i].weight!=0.5));
 359    for (i=1; i<=T2[0].sign; i++) weights = (weights||(T2[i].weight!=0.5));
 360
 361    a0 = (char *) space(alignment[0][0]*4+2);
 362    a1 = (char *) space(alignment[0][0]*4+2);
 363    for(i=1; i<= alignment[0][0]; i++){
 364       tmp0[0] = '\0'; l0=0;
 365       if(alignment[0][i] > 0) {
 366          encode(T1[alignment[0][i]].type, label);
 367          if(T1[alignment[0][i]].sign > 0) {
 368             tmp0[0] = '(';
 369             tmp0[1] = '\0';
 370          }
 371          strcat(tmp0,label);
 372          if (weights)
 373             sprintf(tmp0+strlen(tmp0), "%d",
 374                     (int)(2*T1[alignment[0][i]].weight));
 375
 376          if(T1[alignment[0][i]].sign < 0) strcat(tmp0, ")");
 377          l0 = strlen(tmp0);
 378       }
 379       tmp1[0]= '\0'; l1=0;
 380       if(alignment[1][i] > 0) {
 381          encode(T2[alignment[1][i]].type, label);
 382          if(T2[alignment[1][i]].sign > 0) {
 383             tmp1[0] = '(';
 384             tmp1[1] = '\0';
 385          }
 386          strcat(tmp1,label);
 387          if (weights)
 388             sprintf(tmp1+strlen(tmp1), "%d",
 389                     (int)(2*T2[alignment[1][i]].weight));
 390
 391          if(T2[alignment[1][i]].sign < 0) strcat(tmp1, ")");
 392          l1 = strlen(tmp1);
 393       }
 394       ltmp = MAX2(l0,l1);
 395       for (j=l0; j<ltmp; j++) tmp0[j] = '_';
 396       for (j=l1; j<ltmp; j++) tmp1[j] = '_';
 397       tmp0[ltmp] = '\0'; tmp1[ltmp] = '\0';
 398
 399       strcat(a0,tmp0); strcat(a1,tmp1);
 400       ltmp = strlen(a0);
 401    }
 402    if (aligned_line[0]!= NULL) { free(aligned_line[0]); aligned_line[0]= NULL;}
 403    if (aligned_line[1]!= NULL) { free(aligned_line[1]); aligned_line[1]= NULL;}
 404    aligned_line[0] = strdup(a0);
 405    free(a0);
 406    aligned_line[1] = strdup(a1);
 407    free(a1);
 408 }
 409
 410 /*---------------------------------------------------------------------------*/
 411
 412
 413 PUBLIC void print_swString(swString *x)
 414 {
 415    int i;
 416    for (i=0; i<=x[0].sign; i++)
 417       printf("(%d,%d,%f\n) ", x[i].type, x[i].sign, x[i].weight );
 418    printf("\n");
 419 }
 420
 421 /*---------------------------------------------------------------------------*/
 422
 423 PUBLIC void print_alignment_list(void)
 424 {
 425    int i;
 426    printf("\n");
 427    for (i=1; i<= alignment[0][0]; i++)
 428       printf("%3d ",alignment[0][i]);
 429       printf("\n");
 430    for (i=1; i<= alignment[0][0]; i++)
 431       printf("%3d ",alignment[1][i]);
 432    printf("\n");
 433 }