--- /dev/null
+/*
+ Tree edit distances for RNA secondary structures
+ Walter Fontana, Ivo L Hofacker, Peter F Stadler
+ Vienna RNA Package
+*/
+/* Last changed Time-stamp: <97/10/27 15:23:48 ivo> */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+#include "edit_cost.h"
+#include "dist_vars.h"
+#include "utils.h"
+static char rcsid[] = "$Id: treedist.c,v 1.3 1997/11/03 10:39:43 ivo Rel $";
+
+#define PRIVATE static
+#define PUBLIC
+
+#define MNODES 4000 /* Maximal number of nodes for alignment */
+
+PUBLIC Tree *make_tree(char *struc);
+PUBLIC float tree_edit_distance(Tree *T1, Tree *T2);
+PUBLIC void print_tree(Tree *t);
+PUBLIC void free_tree(Tree *t);
+
+
+PRIVATE void tree_dist(int i, int j);
+PRIVATE int edit_cost(int i, int j);
+PRIVATE int *make_keyroots(Postorder_list *pl);
+PRIVATE void sort(int n, int *ra);
+PRIVATE Postorder_list *make_postorder_list(char *struc);
+PRIVATE int decode(char *id);
+PRIVATE int number_of_nodes(char *struc);
+PRIVATE void encode(int type, char *label);
+PRIVATE void print_keyroots(int *keyroots);
+PRIVATE void print_postorder_list(Postorder_list *pl);
+PRIVATE void backtracking(void);
+PRIVATE void sprint_aligned_trees(void);
+
+
+PRIVATE Tree *tree1, *tree2;
+PRIVATE int **tdist; /* contains distances between subtrees */
+PRIVATE int **fdist; /* contains distances between forests */
+PRIVATE int *alignment[2]; /* contains numeric information on the alignment:
+ alignment[0][p], aligment[1][p] are aligned postions.
+ INDELs have one 0.
+ alignment[0][0] contains the length of the alignment. */
+
+/*---------------------------------------------------------------------------*/
+
+PUBLIC float tree_edit_distance(Tree *T1, Tree *T2)
+{
+ int i1,j1,i,j, dist;
+ int n1, n2;
+
+ if (cost_matrix==0) EditCost = &UsualCost;
+ else EditCost = &ShapiroCost;
+
+ n1 = T1->postorder_list[0].sons;
+ n2 = T2->postorder_list[0].sons;
+
+ tdist = (int **) space(sizeof(int *) * (n1+1));
+ fdist = (int **) space(sizeof(int *) * (n1+1));
+ for (i=0; i<=n1; i++) {
+ tdist[i] = (int *) space(sizeof(int) * (n2+1));
+ fdist[i] = (int *) space(sizeof(int) * (n2+1));
+ }
+
+ tree1 = T1; tree2 = T2;
+
+ for (i1 = 1; i1 <= T1->keyroots[0]; i1++) {
+ i = T1->keyroots[i1];
+ for (j1 = 1; j1 <= T2->keyroots[0]; j1++) {
+ j = T2->keyroots[j1];
+
+ tree_dist(i,j);
+ }
+ }
+
+ if (edit_backtrack) {
+
+ if ((n1>MNODES)||(n2>MNODES)) nrerror("tree too large for alignment");
+
+ alignment[0] = (int *) space((n1+1)*sizeof(int));
+ alignment[1] = (int *) space((n2+1)*sizeof(int));
+
+ backtracking();
+ sprint_aligned_trees();
+ free(alignment[0]);
+ free(alignment[1]);
+ }
+ dist = tdist[n1][n2];
+ for (i=0; i<=n1; i++) {
+ free(tdist[i]);
+ free(fdist[i]);
+ }
+ free(tdist);
+ free(fdist);
+
+ return (float) dist;
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void tree_dist(int i, int j)
+{
+ int li,lj,i1,j1,i1_1,j1_1,li1_1,lj1_1,f1,f2,f3,f;
+ int cost, lleaf_i1, lleaf_j1;
+
+ fdist[0][0] = 0;
+
+ li = tree1->postorder_list[i].leftmostleaf;
+ lj = tree2->postorder_list[j].leftmostleaf;
+
+ for (i1 = li; i1 <= i; i1++) {
+ i1_1 = (li == i1 ? 0 : i1-1);
+ fdist[i1][0] = fdist[i1_1][0] + edit_cost(i1, 0);
+ }
+
+ for (j1 = lj; j1 <= j; j1++) {
+ j1_1 = (lj == j1 ? 0 : j1-1);
+ fdist[0][j1] = fdist[0][j1_1] + edit_cost(0, j1);
+ }
+
+ for (i1 = li; i1 <= i; i1++) {
+
+ lleaf_i1 = tree1->postorder_list[i1].leftmostleaf;
+ li1_1 = (li > lleaf_i1-1 ? 0 : lleaf_i1-1);
+ i1_1 = (i1 == li ? 0: i1-1);
+ cost = edit_cost(i1, 0);
+
+ for (j1 = lj; j1 <= j; j1++) {
+
+ lleaf_j1 = tree2->postorder_list[j1].leftmostleaf;
+ j1_1 = (j1 == lj ? 0: j1-1);
+
+ f1 = fdist[i1_1][j1] + cost;
+ f2 = fdist[i1][j1_1] + edit_cost(0, j1);
+
+ f = f1 < f2 ? f1 : f2;
+
+ if (lleaf_i1 == li && lleaf_j1 == lj) {
+
+ f3 = fdist[i1_1][j1_1] + edit_cost(i1, j1);
+
+ fdist[i1][j1] = f3 < f ? f3 : f;
+
+ tdist[i1][j1] = fdist[i1][j1]; /* store in array permanently */
+ }
+ else {
+ lj1_1 = (lj > lleaf_j1-1 ? 0 : lleaf_j1-1);
+
+ f3 = fdist[li1_1][lj1_1] + tdist[i1][j1];
+
+ fdist[i1][j1] = f3 < f ? f3 : f;
+ }
+ }
+ }
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE int edit_cost(int i, int j)
+{
+ int c, diff, cd, min, a, b;
+
+ c = (*EditCost)[tree1->postorder_list[i].type][tree2->postorder_list[j].type];
+
+ diff = abs((a=tree1->postorder_list[i].weight) - (b=tree2->postorder_list[j].weight));
+
+ min = (a < b ? a: b);
+ if (min == a) cd = (*EditCost)[0][tree2->postorder_list[j].type];
+ else cd = (*EditCost)[0][tree1->postorder_list[i].type];
+
+ return (c * min + cd * diff);
+
+}
+
+/*---------------------------------------------------------------------------*/
+
+PUBLIC Tree *make_tree(char *struc)
+{
+ Tree *tree;
+
+ tree = (Tree *) space(sizeof(Tree));
+
+ tree->postorder_list = make_postorder_list(struc);
+ tree->keyroots = make_keyroots(tree->postorder_list);
+
+ return (tree);
+}
+
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE int *make_keyroots(Postorder_list *pl)
+{
+ int i, k, keys;
+ int *keyroots;
+
+ keyroots = (int *) space(sizeof(int)*(pl[0].sons+1));
+ keys = 0;
+
+ for (i = 1; i <= pl[0].sons; i++) {
+ if (!pl[i].sons) {
+
+ /* leaf */
+
+ k = pl[0].sons;
+ while (pl[k].leftmostleaf != i) k--;
+ keyroots[++keys] = k;
+ }
+ }
+
+ sort(keys, keyroots);
+ keyroots[0] = keys;
+
+ return (keyroots);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void sort(int n, int *ra) /* heap sort, indices are 1..n !!! */
+{
+ int l,j,ir,i;
+ int rra;
+
+ if (n == 1) return;
+
+ l = (n >> 1)+1;
+ ir = n;
+ for (;;) {
+ if (l > 1)
+ rra = ra[--l];
+ else {
+ rra = ra[ir];
+ ra[ir] = ra[1];
+ if (--ir == 1) {
+ ra[1] = rra;
+ return;
+ }
+ }
+ i = l;
+ j = l << 1;
+ while (j <= ir) {
+ if (j < ir && ra[j] < ra[j+1]) ++j;
+ if (rra < ra[j]) {
+ ra[i] = ra[j];
+ j += (i = j);
+ }
+ else j = ir+1;
+ }
+ ra[i] = rra;
+ }
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE Postorder_list *make_postorder_list(char *struc)
+
+ /*
+ Convention for structure representation "struc":
+ Nodes are one pair of matching parentheses, with the type and possibly
+ a weight of the node immediately preceding the closing parentheses.
+
+ Types:
+
+ U....unpaired
+ P....paired
+ H....hairpin loop
+ B....bulge loop
+ I....internal loop
+ M....multiloop
+ S....stack
+ R....virtual root
+
+ Example:
+
+ .((..(((...)))..((..)))). in usual notation becomes:
+
+ full tree:
+ ((U)(((U)(U)((((U)(U)(U)P)P)P)(U)(U)(((U)(U)P)P)P)P)(U)R)
+ HIT:
+ ((U1)((U2)((U3)P3)(U2)((U2)P2)P2)(U1)R)
+ Shapiro:
+ (((((H)S)((H)S)M)S)R)
+
+ */
+{
+ int paren, i, l, order, local_order, w, sons, count;
+ int n_nodes, p;
+ char id[100];
+ Postorder_list *pl;
+ int match_pos[MNODES], match_order[MNODES];
+
+
+ n_nodes = number_of_nodes(struc);
+ if (n_nodes>MNODES) nrerror("structure too long in make_postorder_list");
+ pl = (Postorder_list *) space(sizeof(Postorder_list)*(n_nodes+1));
+ pl[0].sons = n_nodes;
+
+ paren = 1;
+ match_pos[paren] = 0;
+ match_order[paren] = 0;
+
+ i = 1;
+ l = 0;
+ order = 0;
+
+ while (paren) {
+ switch (struc[i]) {
+ case '(':
+ match_pos[++paren] = i;
+ match_order[paren] = order;
+ break;
+ case ')':
+ order++;
+ id[l] = '\0';
+ l = 0;
+ while (isalpha((int) id[l])) l++;
+ if (id[l]) sscanf(id+l, "%d", &w);
+ else w = 1;
+ id[l] = '\0';
+ pl[order].type = decode(id);
+ pl[order].weight = w;
+ pl[order].leftmostleaf = match_order[paren]+1;
+
+ sons = count = 0;
+ local_order = match_order[paren];
+ for (p = match_pos[paren]+1; p < i; p++) {
+ if (struc[p] == '(') count++;
+ else if (struc[p] == ')') {
+ local_order++;
+ if (count == 1) {
+ sons++;
+ pl[local_order].father = order;
+ }
+ count--;
+ }
+ }
+
+ pl[order].sons = sons;
+ paren--;
+ l = 0;
+ break;
+ default:
+ id[l++] = struc[i];
+ break;
+ }
+ i++;
+ }
+
+ return (pl);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE int decode(char *id)
+{
+ int n, quit, i;
+ char label[100], *code;
+
+ n = 0;
+
+ quit = 0;
+ code = coding;
+
+ while (!quit) {
+ for (i = 0; code[i] != sep; i++) {
+ if (code[i] == '\0') {
+ quit = 1;
+ break;
+ }
+ label[i] = code[i];
+ }
+ label[i] = '\0';
+ if (strcmp(id, label) == 0) return (n);
+ code += (i+1);
+ n++;
+ }
+
+ fprintf(stderr,"Syntax error: node identifier \"%s\" not found "
+ "in coding string \"%s\"\n", id, coding);
+ fprintf(stderr, "Exiting...");
+ exit(0);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void encode(int type, char label[])
+{
+ int i, l;
+
+ l = 0;
+ for (i = 0; i < type; i++) {
+ while (coding[l] != sep && coding[l]) l++;
+ l++;
+ }
+
+ for (i = 0; coding[l+i] != sep; i++) {
+ if (coding[l+i] == '\0') break;
+ label[i] = coding[l+i];
+ }
+ label[i] = '\0';
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE int number_of_nodes(char *struc)
+{
+ int l, c, i;
+
+ l = strlen(struc);
+ for (c = 0, i = 0; i < l; i++) if (struc[i] == ')') c++;
+ return (c);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void print_keyroots(int *keyroots)
+{
+ int i;
+
+ printf("---> key roots <---\n\n");
+
+ printf("entries: %d\n", keyroots[0]);
+ printf("{");
+ for (i = 1; i <= keyroots[0]; i++) printf(" %d", keyroots[i]);
+ printf(" }\n\n");
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void print_postorder_list(Postorder_list *pl)
+{
+ register i;
+ char label[100];
+
+ printf("---> postorder list <---\n\n");
+
+ for (i = 1; i <= pl[0].sons; i++) {
+ printf(" postorder: %3d\n", i);
+ *label = '\0';
+ encode(pl[i].type, label);
+ printf(" type: %3d (%s)\n", pl[i].type, label);
+ printf(" weight: %3d\n", pl[i].weight);
+ printf(" father: %3d\n", pl[i].father);
+ printf(" sons: %3d\n", pl[i].sons);
+ printf("leftmost leaf: %3d\n", pl[i].leftmostleaf);
+ printf("\n");
+ }
+}
+
+/*---------------------------------------------------------------------------*/
+
+PUBLIC void print_tree(Tree *t)
+{
+ print_postorder_list(t->postorder_list);
+ print_keyroots(t->keyroots);
+ fflush(stdout);
+}
+
+/*---------------------------------------------------------------------------*/
+
+PUBLIC void free_tree(Tree *t)
+{
+ free(t->postorder_list);
+ free(t->keyroots);
+ free(t);
+}
+
+/*---------------------------------------------------------------------------*/
+
+
+PRIVATE void backtracking(void)
+{
+ int li,lj,i1,j1,i1_1,j1_1,li1_1,lj1_1,f;
+ int cost, lleaf_i1, lleaf_j1, ss, i,j,k;
+ struct {int i,j;} sector[MNODES/2];
+
+ ss=0;
+
+ i=i1=tree1->postorder_list[0].sons;
+ j=j1=tree2->postorder_list[0].sons;
+
+ start:
+ li = tree1->postorder_list[i].leftmostleaf;
+ lj = tree2->postorder_list[j].leftmostleaf;
+
+
+ while ((i1>=li)&&(j1>=lj)) {
+
+ lleaf_i1 = tree1->postorder_list[i1].leftmostleaf;
+ li1_1 = (li > lleaf_i1-1 ? 0 : lleaf_i1-1);
+ i1_1 = (i1 == li ? 0: i1-1);
+ lleaf_j1 = tree2->postorder_list[j1].leftmostleaf;
+ lj1_1 = (lj > lleaf_j1-1 ? 0 : lleaf_j1-1);
+ j1_1 = (j1 == lj ? 0: j1-1);
+
+ f = fdist[i1][j1];
+
+ cost = edit_cost(i1, 0);
+ if (f == fdist[i1_1][j1] + cost) {
+ alignment[0][i1]=0;
+ i1=i1_1;
+ }
+ else {
+ if (f == fdist[i1][j1_1] + edit_cost(0, j1)) {
+ alignment[1][j1]=0;
+ j1=j1_1;
+ }
+ else if (lleaf_i1 == li && lleaf_j1 == lj) {
+ alignment[0][i1] = j1;
+ alignment[1][j1] = i1;
+ i1=i1_1; j1=j1_1;
+ } else {
+ sector[ss].i=i1;
+ sector[ss++].j=j1;
+ i1=li1_1;
+ j1=lj1_1;
+ }
+ }
+ }
+ for (; i1>=li; ) {
+ alignment[0][i1]=0;
+ i1 = (i1 == li ? 0: i1-1);
+ }
+ for (; j1>=lj; ) {
+ alignment[1][j1]=0;
+ j1 = (j1 == lj ? 0: j1-1);
+ }
+ while (ss>0) {
+ i1=sector[--ss].i;
+ j1=sector[ss].j;
+ for (k=1; 1; k++) {
+ i = tree1->keyroots[k];
+ if (tree1->postorder_list[i].leftmostleaf ==
+ tree1->postorder_list[i1].leftmostleaf) break;
+ }
+ for (k=1; 1; k++) {
+ j = tree2->keyroots[k];
+ if (tree2->postorder_list[j].leftmostleaf ==
+ tree2->postorder_list[j1].leftmostleaf) break;
+ }
+ tree_dist(i,j);
+ goto start;
+ }
+}
+
+/*---------------------------------------------------------------------------*/
+
+PRIVATE void sprint_aligned_trees(void)
+{
+ int i,j,n1,n2,k,l,p, ni, nj, weights;
+ char t1[2*MNODES+1], t2[2*MNODES+1], a1[8*MNODES], a2[8*MNODES], ll[20], ll1[20];
+
+ weights=0;
+ n1=tree1->postorder_list[0].sons;
+ n2=tree2->postorder_list[0].sons;
+ for (i=1; i<=n1; i++) weights |= (tree1->postorder_list[i].weight!=1);
+ for (i=1; i<=n2; i++) weights |= (tree2->postorder_list[i].weight!=1);
+
+ for (i=n1, l=2*n1-1; i>0; i--) {
+ if (alignment[0][i]!=0) t1[l--]=']';
+ else t1[l--]=')';
+ p=i;
+ while(i==tree1->postorder_list[p].leftmostleaf) {
+ if (alignment[0][p]!=0) t1[l--]='[';
+ else t1[l--]='(';
+ p=tree1->postorder_list[p].father;
+ }
+ }
+ t1[2*n1]='\0';
+ for (j=n2, l=2*n2-1; j>0; j--) {
+ if (alignment[1][j]!=0) t2[l--]=']';
+ else t2[l--]=')';
+ p=j;
+ while(j==tree2->postorder_list[p].leftmostleaf) {
+ if (alignment[1][p]!=0) t2[l--]='[';
+ else t2[l--]='(';
+ p=tree2->postorder_list[p].father;
+ }
+ }
+ t2[2*n2]='\0';
+
+ ni=nj=l=i=j=0;
+ while (t1[i]||t2[j]) {
+ while ((t1[i]=='(')||(t1[i]==')')) {
+ if (t1[i]==')') {
+ ni++;
+ encode(tree1->postorder_list[ni].type, ll);
+ if (weights)
+ sprintf(ll+strlen(ll), "%d", tree1->postorder_list[ni].weight);
+ for (k=0; k< strlen(ll); k++) {
+ a1[l]=ll[k];
+ a2[l++]='_';
+ }
+ a1[l]=')'; a2[l++]='_';
+ } else {
+ a1[l] = t1[i];
+ a2[l++] ='_';
+ }
+ i++;
+ }
+
+ while ((t2[j]=='(')||(t2[j]==')')) {
+ if (t2[j]==')') {
+ nj++;
+ encode(tree2->postorder_list[nj].type, ll);
+ if (weights)
+ sprintf(ll+strlen(ll), "%d", tree2->postorder_list[nj].weight);
+ for (k=0; k< strlen(ll); k++) {
+ a2[l]=ll[k];
+ a1[l++]='_';
+ }
+ a2[l]=')'; a1[l++]='_';
+ } else {
+ a2[l] = t2[j];
+ a1[l++] ='_';
+ }
+ j++;
+ }
+
+ if (t2[j]==']') {
+ ni++; nj++;
+ encode(tree2->postorder_list[nj].type, ll);
+ if (weights)
+ sprintf(ll+strlen(ll), "%d", tree2->postorder_list[nj].weight);
+ encode(tree1->postorder_list[ni].type, ll1);
+ if (weights)
+ sprintf(ll1+strlen(ll1), "%d", tree1->postorder_list[ni].weight);
+ if (strlen(ll)>strlen(ll1))
+ for (k=0; k<strlen(ll)-strlen(ll1); k++) strcat(ll1,"_");
+ if (strlen(ll)<strlen(ll1))
+ for (k=0; k<strlen(ll1)-strlen(ll); k++) strcat(ll,"_");
+ for (k=0; k< strlen(ll); k++) a2[l+k]=ll[k];
+ for (k=0; k< strlen(ll); k++) a1[l+k]=ll1[k];
+ l+=k;
+ a1[l]=a2[l]=')'; l++;
+ i++; j++;
+ } else if (t2[j]=='[') {
+ a1[l]=a2[l]='('; l++;
+ i++; j++;
+ }
+ }
+ a1[l]=a2[l]='\0';
+ if (l>8*MNODES) nrerror("structure too long in sprint_aligned_trees");
+ if (aligned_line[0]!= NULL) free(aligned_line[0]);
+ if (aligned_line[1]!= NULL) free(aligned_line[1]);
+ aligned_line[0] = (char *) space((l+1)*sizeof(char));
+ aligned_line[1] = (char *) space((l+1)*sizeof(char));
+ strcpy(aligned_line[0], a1);
+ strcpy(aligned_line[1], a2);
+}
+
+/*---------------------------------------------------------------------------*/
git://source.jalview.org / jabaws.git / blobdiff