binaries/src/ViennaRNA/lib/treedist.c

   1 /*
   2                 Tree edit distances for RNA secondary structures
   3                 Walter Fontana, Ivo L Hofacker, Peter F Stadler
   4                              Vienna RNA Package
   5 */
   6 /* Last changed Time-stamp: <97/10/27 15:23:48 ivo> */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <string.h>
  11 #include <ctype.h>
  12 #include  "edit_cost.h"
  13 #include  "dist_vars.h"
  14 #include  "utils.h"
  15 static char rcsid[] = "$Id: treedist.c,v 1.3 1997/11/03 10:39:43 ivo Rel $";
  16
  17 #define PRIVATE  static
  18 #define PUBLIC
  19
  20 #define MNODES    4000    /* Maximal number of nodes for alignment    */
  21
  22 PUBLIC  Tree     *make_tree(char *struc);
  23 PUBLIC  float     tree_edit_distance(Tree *T1, Tree *T2);
  24 PUBLIC  void      print_tree(Tree *t);
  25 PUBLIC  void      free_tree(Tree *t);
  26
  27
  28 PRIVATE void            tree_dist(int i, int j);
  29 PRIVATE int             edit_cost(int i, int j);
  30 PRIVATE int            *make_keyroots(Postorder_list *pl);
  31 PRIVATE void            sort(int n, int *ra);
  32 PRIVATE Postorder_list *make_postorder_list(char *struc);
  33 PRIVATE int             decode(char *id);
  34 PRIVATE int             number_of_nodes(char *struc);
  35 PRIVATE void            encode(int type, char *label);
  36 PRIVATE void            print_keyroots(int *keyroots);
  37 PRIVATE void            print_postorder_list(Postorder_list *pl);
  38 PRIVATE void            backtracking(void);
  39 PRIVATE void            sprint_aligned_trees(void);
  40
  41
  42 PRIVATE  Tree  *tree1, *tree2;
  43 PRIVATE  int  **tdist;       /* contains distances between subtrees */
  44 PRIVATE  int  **fdist;       /* contains distances between forests */
  45 PRIVATE  int  *alignment[2]; /* contains numeric information on the alignment:
  46                                 alignment[0][p], aligment[1][p] are aligned postions.
  47                                 INDELs have one 0.
  48                                 alignment[0][0] contains the length of the alignment. */
  49
  50 /*---------------------------------------------------------------------------*/
  51
  52 PUBLIC float tree_edit_distance(Tree *T1, Tree *T2)
  53 {
  54    int i1,j1,i,j, dist;
  55    int n1, n2;
  56
  57    if (cost_matrix==0) EditCost = &UsualCost;
  58    else EditCost = &ShapiroCost;
  59
  60    n1 = T1->postorder_list[0].sons;
  61    n2 = T2->postorder_list[0].sons;
  62
  63    tdist = (int **) space(sizeof(int *) * (n1+1));
  64    fdist = (int **) space(sizeof(int *) * (n1+1));
  65    for (i=0; i<=n1; i++) {
  66       tdist[i] = (int *) space(sizeof(int) * (n2+1));
  67       fdist[i] = (int *) space(sizeof(int) * (n2+1));
  68    }
  69
  70    tree1 = T1;  tree2 = T2;
  71
  72    for (i1 = 1; i1 <= T1->keyroots[0]; i1++) {
  73       i = T1->keyroots[i1];
  74       for (j1 = 1; j1 <= T2->keyroots[0]; j1++) {
  75          j = T2->keyroots[j1];
  76
  77          tree_dist(i,j);
  78       }
  79    }
  80
  81    if (edit_backtrack) {
  82
  83       if ((n1>MNODES)||(n2>MNODES)) nrerror("tree too large for alignment");
  84
  85       alignment[0] = (int *) space((n1+1)*sizeof(int));
  86       alignment[1] = (int *) space((n2+1)*sizeof(int));
  87
  88       backtracking();
  89       sprint_aligned_trees();
  90       free(alignment[0]);
  91       free(alignment[1]);
  92    }
  93    dist = tdist[n1][n2];
  94    for (i=0; i<=n1; i++) {
  95       free(tdist[i]);
  96       free(fdist[i]);
  97    }
  98    free(tdist);
  99    free(fdist);
 100
 101    return (float) dist;
 102 }
 103
 104 /*---------------------------------------------------------------------------*/
 105
 106 PRIVATE void tree_dist(int i, int j)
 107 {
 108    int li,lj,i1,j1,i1_1,j1_1,li1_1,lj1_1,f1,f2,f3,f;
 109    int cost, lleaf_i1, lleaf_j1;
 110
 111    fdist[0][0] = 0;
 112
 113    li = tree1->postorder_list[i].leftmostleaf;
 114    lj = tree2->postorder_list[j].leftmostleaf;
 115
 116    for (i1 = li; i1 <= i; i1++) {
 117       i1_1 = (li == i1 ? 0 : i1-1);
 118       fdist[i1][0] = fdist[i1_1][0] + edit_cost(i1, 0);
 119    }
 120
 121    for (j1 = lj; j1 <= j; j1++) {
 122       j1_1 = (lj == j1 ? 0 : j1-1);
 123       fdist[0][j1] = fdist[0][j1_1] + edit_cost(0, j1);
 124    }
 125
 126    for (i1 = li; i1 <= i; i1++) {
 127
 128       lleaf_i1 = tree1->postorder_list[i1].leftmostleaf;
 129       li1_1 = (li > lleaf_i1-1 ? 0 : lleaf_i1-1);
 130       i1_1 = (i1 == li ? 0: i1-1);
 131       cost = edit_cost(i1, 0);
 132
 133       for (j1 = lj; j1 <= j; j1++) {
 134
 135          lleaf_j1 = tree2->postorder_list[j1].leftmostleaf;
 136          j1_1 = (j1 == lj ? 0: j1-1);
 137
 138          f1 = fdist[i1_1][j1] + cost;
 139          f2 = fdist[i1][j1_1] + edit_cost(0, j1);
 140
 141          f = f1 < f2 ? f1 : f2;
 142
 143          if (lleaf_i1 == li && lleaf_j1 == lj)   {
 144
 145             f3 = fdist[i1_1][j1_1] + edit_cost(i1, j1);
 146
 147             fdist[i1][j1] = f3 < f ? f3 : f;
 148
 149             tdist[i1][j1] = fdist[i1][j1];   /* store in array permanently */
 150          }
 151          else {
 152             lj1_1 = (lj > lleaf_j1-1 ? 0 : lleaf_j1-1);
 153
 154             f3 = fdist[li1_1][lj1_1] + tdist[i1][j1];
 155
 156             fdist[i1][j1] = f3 < f ? f3 : f;
 157          }
 158       }
 159    }
 160 }
 161
 162 /*---------------------------------------------------------------------------*/
 163
 164 PRIVATE int edit_cost(int i, int j)
 165 {
 166    int  c, diff, cd, min, a, b;
 167
 168    c = (*EditCost)[tree1->postorder_list[i].type][tree2->postorder_list[j].type];
 169
 170    diff = abs((a=tree1->postorder_list[i].weight) - (b=tree2->postorder_list[j].weight));
 171
 172    min = (a < b ? a: b);
 173    if (min == a) cd = (*EditCost)[0][tree2->postorder_list[j].type];
 174    else          cd = (*EditCost)[0][tree1->postorder_list[i].type];
 175
 176    return (c * min + cd * diff);
 177
 178 }
 179
 180 /*---------------------------------------------------------------------------*/
 181
 182 PUBLIC Tree *make_tree(char *struc)
 183 {
 184    Tree *tree;
 185
 186    tree = (Tree *) space(sizeof(Tree));
 187
 188    tree->postorder_list = make_postorder_list(struc);
 189    tree->keyroots       = make_keyroots(tree->postorder_list);
 190
 191    return (tree);
 192 }
 193
 194
 195 /*---------------------------------------------------------------------------*/
 196
 197 PRIVATE int  *make_keyroots(Postorder_list *pl)
 198 {
 199    int   i, k, keys;
 200    int  *keyroots;
 201
 202    keyroots = (int *) space(sizeof(int)*(pl[0].sons+1));
 203    keys      = 0;
 204
 205    for (i = 1; i <= pl[0].sons; i++) {
 206       if (!pl[i].sons) {
 207
 208          /* leaf */
 209
 210          k = pl[0].sons;
 211          while (pl[k].leftmostleaf != i) k--;
 212          keyroots[++keys] = k;
 213       }
 214    }
 215
 216    sort(keys, keyroots);
 217    keyroots[0] = keys;
 218
 219    return (keyroots);
 220 }
 221
 222 /*---------------------------------------------------------------------------*/
 223
 224 PRIVATE void sort(int n, int *ra)       /* heap sort,  indices are 1..n !!! */
 225 {
 226    int l,j,ir,i;
 227    int rra;
 228
 229    if (n == 1) return;
 230
 231    l = (n >> 1)+1;
 232    ir = n;
 233    for (;;) {
 234       if (l > 1)
 235          rra = ra[--l];
 236       else {
 237          rra = ra[ir];
 238          ra[ir] = ra[1];
 239          if (--ir == 1) {
 240             ra[1] = rra;
 241             return;
 242          }
 243       }
 244       i = l;
 245       j = l << 1;
 246       while (j <= ir) {
 247          if (j < ir && ra[j] < ra[j+1]) ++j;
 248          if (rra < ra[j]) {
 249             ra[i] = ra[j];
 250             j += (i = j);
 251          }
 252          else j = ir+1;
 253       }
 254       ra[i] = rra;
 255    }
 256 }
 257
 258 /*---------------------------------------------------------------------------*/
 259
 260 PRIVATE Postorder_list *make_postorder_list(char *struc)
 261
 262      /*
 263        Convention for structure representation "struc":
 264        Nodes are one pair of matching parentheses, with the type and possibly
 265        a weight of the node immediately preceding the closing parentheses.
 266
 267        Types:
 268
 269        U....unpaired
 270        P....paired
 271        H....hairpin loop
 272        B....bulge loop
 273        I....internal loop
 274        M....multiloop
 275        S....stack
 276        R....virtual root
 277
 278        Example:
 279
 280        .((..(((...)))..((..)))). in usual notation becomes:
 281
 282        full tree:
 283        ((U)(((U)(U)((((U)(U)(U)P)P)P)(U)(U)(((U)(U)P)P)P)P)(U)R)
 284        HIT:
 285        ((U1)((U2)((U3)P3)(U2)((U2)P2)P2)(U1)R)
 286        Shapiro:
 287        (((((H)S)((H)S)M)S)R)
 288
 289        */
 290 {
 291    int  paren, i, l, order, local_order, w, sons, count;
 292    int  n_nodes, p;
 293    char id[100];
 294    Postorder_list *pl;
 295    int  match_pos[MNODES], match_order[MNODES];
 296
 297
 298    n_nodes = number_of_nodes(struc);
 299    if (n_nodes>MNODES) nrerror("structure too long in make_postorder_list");
 300    pl = (Postorder_list *) space(sizeof(Postorder_list)*(n_nodes+1));
 301    pl[0].sons = n_nodes;
 302
 303    paren = 1;
 304    match_pos[paren]   = 0;
 305    match_order[paren] = 0;
 306
 307    i     = 1;
 308    l     = 0;
 309    order = 0;
 310
 311    while (paren) {
 312       switch (struc[i]) {
 313        case '(':
 314          match_pos[++paren] = i;
 315          match_order[paren] = order;
 316          break;
 317        case ')':
 318          order++;
 319          id[l] = '\0';
 320          l = 0;
 321          while (isalpha((int) id[l])) l++;
 322          if (id[l]) sscanf(id+l, "%d", &w);
 323          else  w = 1;
 324          id[l] = '\0';
 325          pl[order].type         = decode(id);
 326          pl[order].weight       = w;
 327          pl[order].leftmostleaf = match_order[paren]+1;
 328
 329          sons = count = 0;
 330          local_order  = match_order[paren];
 331          for (p = match_pos[paren]+1; p < i; p++) {
 332             if (struc[p] == '(') count++;
 333             else if (struc[p] == ')') {
 334                local_order++;
 335                if (count == 1) {
 336                   sons++;
 337                   pl[local_order].father = order;
 338                }
 339                count--;
 340             }
 341          }
 342
 343          pl[order].sons         = sons;
 344          paren--;
 345          l = 0;
 346          break;
 347        default:
 348          id[l++] = struc[i];
 349          break;
 350       }
 351       i++;
 352    }
 353
 354    return (pl);
 355 }
 356
 357 /*---------------------------------------------------------------------------*/
 358
 359 PRIVATE int decode(char *id)
 360 {
 361    int   n, quit, i;
 362    char  label[100], *code;
 363
 364    n = 0;
 365
 366    quit = 0;
 367    code = coding;
 368
 369    while (!quit) {
 370       for (i = 0; code[i] != sep; i++) {
 371          if (code[i] == '\0') {
 372             quit = 1;
 373             break;
 374          }
 375          label[i] = code[i];
 376       }
 377       label[i] = '\0';
 378       if (strcmp(id, label) == 0) return (n);
 379       code += (i+1);
 380       n++;
 381    }
 382
 383    fprintf(stderr,"Syntax error: node identifier \"%s\" not found "
 384                   "in coding string \"%s\"\n", id, coding);
 385    fprintf(stderr, "Exiting...");
 386    exit(0);
 387 }
 388
 389 /*---------------------------------------------------------------------------*/
 390
 391 PRIVATE void encode(int type, char label[])
 392 {
 393    int   i, l;
 394
 395    l = 0;
 396    for (i = 0; i < type; i++) {
 397       while (coding[l] != sep && coding[l]) l++;
 398       l++;
 399    }
 400
 401    for (i = 0; coding[l+i] != sep; i++) {
 402       if (coding[l+i] == '\0') break;
 403       label[i] = coding[l+i];
 404    }
 405    label[i] = '\0';
 406 }
 407
 408 /*---------------------------------------------------------------------------*/
 409
 410 PRIVATE int number_of_nodes(char *struc)
 411 {
 412    int  l, c, i;
 413
 414    l = strlen(struc);
 415    for (c = 0, i = 0; i < l; i++) if (struc[i] == ')') c++;
 416    return (c);
 417 }
 418
 419 /*---------------------------------------------------------------------------*/
 420
 421 PRIVATE void print_keyroots(int *keyroots)
 422 {
 423    int i;
 424
 425    printf("--->  key roots  <---\n\n");
 426
 427    printf("entries: %d\n", keyroots[0]);
 428    printf("{");
 429    for (i = 1; i <= keyroots[0]; i++) printf(" %d", keyroots[i]);
 430    printf(" }\n\n");
 431 }
 432
 433 /*---------------------------------------------------------------------------*/
 434
 435 PRIVATE void print_postorder_list(Postorder_list *pl)
 436 {
 437    register i;
 438    char     label[100];
 439
 440    printf("--->  postorder list  <---\n\n");
 441
 442    for (i = 1; i <= pl[0].sons; i++) {
 443       printf("    postorder: %3d\n", i);
 444       *label = '\0';
 445       encode(pl[i].type, label);
 446       printf("         type: %3d (%s)\n", pl[i].type, label);
 447       printf("       weight: %3d\n", pl[i].weight);
 448       printf("       father: %3d\n", pl[i].father);
 449       printf("         sons: %3d\n", pl[i].sons);
 450       printf("leftmost leaf: %3d\n", pl[i].leftmostleaf);
 451       printf("\n");
 452    }
 453 }
 454
 455 /*---------------------------------------------------------------------------*/
 456
 457 PUBLIC void print_tree(Tree *t)
 458 {
 459    print_postorder_list(t->postorder_list);
 460    print_keyroots(t->keyroots);
 461    fflush(stdout);
 462 }
 463
 464 /*---------------------------------------------------------------------------*/
 465
 466 PUBLIC void free_tree(Tree *t)
 467 {
 468    free(t->postorder_list);
 469    free(t->keyroots);
 470    free(t);
 471 }
 472
 473 /*---------------------------------------------------------------------------*/
 474
 475
 476 PRIVATE void backtracking(void)
 477 {
 478    int li,lj,i1,j1,i1_1,j1_1,li1_1,lj1_1,f;
 479    int cost, lleaf_i1, lleaf_j1, ss, i,j,k;
 480    struct {int i,j;} sector[MNODES/2];
 481
 482    ss=0;
 483
 484    i=i1=tree1->postorder_list[0].sons;
 485    j=j1=tree2->postorder_list[0].sons;
 486
 487  start:
 488    li = tree1->postorder_list[i].leftmostleaf;
 489    lj = tree2->postorder_list[j].leftmostleaf;
 490
 491
 492    while ((i1>=li)&&(j1>=lj)) {
 493
 494       lleaf_i1 = tree1->postorder_list[i1].leftmostleaf;
 495       li1_1 = (li > lleaf_i1-1 ? 0 : lleaf_i1-1);
 496       i1_1 = (i1 == li ? 0: i1-1);
 497       lleaf_j1 = tree2->postorder_list[j1].leftmostleaf;
 498       lj1_1 = (lj > lleaf_j1-1 ? 0 : lleaf_j1-1);
 499       j1_1 = (j1 == lj ? 0: j1-1);
 500
 501       f = fdist[i1][j1];
 502
 503       cost = edit_cost(i1, 0);
 504       if (f == fdist[i1_1][j1] + cost) {
 505          alignment[0][i1]=0;
 506          i1=i1_1;
 507       }
 508       else {
 509          if (f ==  fdist[i1][j1_1] + edit_cost(0, j1)) {
 510             alignment[1][j1]=0;
 511             j1=j1_1;
 512          }
 513          else if (lleaf_i1 == li && lleaf_j1 == lj) {
 514             alignment[0][i1] = j1;
 515             alignment[1][j1] = i1;
 516             i1=i1_1; j1=j1_1;
 517          } else  {
 518             sector[ss].i=i1;
 519             sector[ss++].j=j1;
 520             i1=li1_1;
 521             j1=lj1_1;
 522          }
 523       }
 524    }
 525    for (; i1>=li; ) {
 526       alignment[0][i1]=0;
 527       i1 = (i1 == li ? 0: i1-1);
 528    }
 529    for (; j1>=lj; ) {
 530       alignment[1][j1]=0;
 531       j1 = (j1 == lj ? 0: j1-1);
 532    }
 533    while (ss>0) {
 534       i1=sector[--ss].i;
 535       j1=sector[ss].j;
 536       for (k=1; 1; k++) {
 537          i = tree1->keyroots[k];
 538          if (tree1->postorder_list[i].leftmostleaf ==
 539              tree1->postorder_list[i1].leftmostleaf) break;
 540       }
 541       for (k=1; 1; k++) {
 542          j = tree2->keyroots[k];
 543          if (tree2->postorder_list[j].leftmostleaf ==
 544              tree2->postorder_list[j1].leftmostleaf) break;
 545       }
 546       tree_dist(i,j);
 547       goto start;
 548    }
 549 }
 550
 551 /*---------------------------------------------------------------------------*/
 552
 553 PRIVATE void sprint_aligned_trees(void)
 554 {
 555    int i,j,n1,n2,k,l,p, ni, nj, weights;
 556    char t1[2*MNODES+1], t2[2*MNODES+1], a1[8*MNODES], a2[8*MNODES], ll[20], ll1[20];
 557
 558    weights=0;
 559    n1=tree1->postorder_list[0].sons;
 560    n2=tree2->postorder_list[0].sons;
 561    for (i=1; i<=n1; i++) weights |= (tree1->postorder_list[i].weight!=1);
 562    for (i=1; i<=n2; i++) weights |= (tree2->postorder_list[i].weight!=1);
 563
 564    for (i=n1, l=2*n1-1; i>0; i--) {
 565       if (alignment[0][i]!=0) t1[l--]=']';
 566       else t1[l--]=')';
 567       p=i;
 568       while(i==tree1->postorder_list[p].leftmostleaf) {
 569          if (alignment[0][p]!=0) t1[l--]='[';
 570          else t1[l--]='(';
 571          p=tree1->postorder_list[p].father;
 572       }
 573    }
 574    t1[2*n1]='\0';
 575    for (j=n2, l=2*n2-1; j>0; j--) {
 576       if (alignment[1][j]!=0) t2[l--]=']';
 577       else t2[l--]=')';
 578       p=j;
 579       while(j==tree2->postorder_list[p].leftmostleaf) {
 580          if (alignment[1][p]!=0) t2[l--]='[';
 581          else t2[l--]='(';
 582          p=tree2->postorder_list[p].father;
 583       }
 584    }
 585    t2[2*n2]='\0';
 586
 587    ni=nj=l=i=j=0;
 588    while (t1[i]||t2[j]) {
 589       while ((t1[i]=='(')||(t1[i]==')')) {
 590          if (t1[i]==')') {
 591             ni++;
 592             encode(tree1->postorder_list[ni].type, ll);
 593             if (weights)
 594                sprintf(ll+strlen(ll), "%d", tree1->postorder_list[ni].weight);
 595             for (k=0; k< strlen(ll); k++) {
 596                a1[l]=ll[k];
 597                a2[l++]='_';
 598             }
 599             a1[l]=')'; a2[l++]='_';
 600          } else {
 601             a1[l] = t1[i];
 602             a2[l++] ='_';
 603          }
 604          i++;
 605       }
 606
 607       while ((t2[j]=='(')||(t2[j]==')')) {
 608          if (t2[j]==')') {
 609             nj++;
 610             encode(tree2->postorder_list[nj].type, ll);
 611             if (weights)
 612                sprintf(ll+strlen(ll), "%d", tree2->postorder_list[nj].weight);
 613             for (k=0; k< strlen(ll); k++) {
 614                a2[l]=ll[k];
 615                a1[l++]='_';
 616             }
 617             a2[l]=')'; a1[l++]='_';
 618          } else {
 619             a2[l] = t2[j];
 620             a1[l++] ='_';
 621          }
 622          j++;
 623       }
 624
 625       if (t2[j]==']') {
 626          ni++; nj++;
 627          encode(tree2->postorder_list[nj].type, ll);
 628          if (weights)
 629             sprintf(ll+strlen(ll), "%d", tree2->postorder_list[nj].weight);
 630          encode(tree1->postorder_list[ni].type, ll1);
 631          if (weights)
 632             sprintf(ll1+strlen(ll1), "%d", tree1->postorder_list[ni].weight);
 633          if (strlen(ll)>strlen(ll1))
 634             for (k=0; k<strlen(ll)-strlen(ll1); k++) strcat(ll1,"_");
 635          if (strlen(ll)<strlen(ll1))
 636             for (k=0; k<strlen(ll1)-strlen(ll); k++) strcat(ll,"_");
 637          for (k=0; k< strlen(ll); k++) a2[l+k]=ll[k];
 638          for (k=0; k< strlen(ll); k++) a1[l+k]=ll1[k];
 639          l+=k;
 640          a1[l]=a2[l]=')'; l++;
 641          i++; j++;
 642       } else if (t2[j]=='[') {
 643          a1[l]=a2[l]='('; l++;
 644          i++; j++;
 645       }
 646    }
 647    a1[l]=a2[l]='\0';
 648    if (l>8*MNODES) nrerror("structure too long in sprint_aligned_trees");
 649    if (aligned_line[0]!= NULL)  free(aligned_line[0]);
 650    if (aligned_line[1]!= NULL)  free(aligned_line[1]);
 651    aligned_line[0] = (char *) space((l+1)*sizeof(char));
 652    aligned_line[1] = (char *) space((l+1)*sizeof(char));
 653    strcpy(aligned_line[0], a1);
 654    strcpy(aligned_line[1], a2);
 655 }
 656
 657 /*---------------------------------------------------------------------------*/