binaries/src/ViennaRNA/lib/findpath.c

   1 /* Last changed Time-stamp: <2008-10-09 16:23:36 ivo> */
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <string.h>
   5 #include <limits.h>
   6
   7 #include "findpath.h"
   8 #include "fold.h"
   9 #include "fold_vars.h"
  10 #include "utils.h"
  11 #include "pair_mat.h"
  12
  13 #ifdef _OPENMP
  14 #include <omp.h>
  15 #endif
  16
  17 #define LOOP_EN
  18
  19 static char rcsid[] = "$Id: findpath.c,v 1.2 2008/10/09 15:42:45 ivo Exp $";
  20
  21 /*
  22 #################################
  23 # GLOBAL VARIABLES              #
  24 #################################
  25 */
  26
  27 /*
  28 #################################
  29 # PRIVATE VARIABLES             #
  30 #################################
  31 */
  32 PRIVATE const char  *seq=NULL;
  33 PRIVATE short       *S=NULL, *S1=NULL;
  34 PRIVATE int         BP_dist;
  35 PRIVATE move_t      *path=NULL;
  36 PRIVATE int         path_fwd; /* 1: struc1->struc2, else struc2 -> struc1 */
  37
  38 #ifdef _OPENMP
  39
  40 /* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
  41 */
  42 #pragma omp threadprivate(seq, S, S1, BP_dist, path, path_fwd)
  43
  44 #endif
  45
  46 /*
  47 #################################
  48 # PRIVATE FUNCTION DECLARATIONS #
  49 #################################
  50 */
  51 PRIVATE int     *pair_table_to_loop_index (short *pt);
  52 PRIVATE move_t  *copy_moves(move_t *mvs);
  53 PRIVATE int     compare_ptable(const void *A, const void *B);
  54 PRIVATE int     compare_energy(const void *A, const void *B);
  55 PRIVATE int     compare_moves_when(const void *A, const void *B);
  56 PRIVATE void    free_intermediate(intermediate_t *i);
  57 PRIVATE int     find_path_once(const char *struc1, const char *struc2, int maxE, int maxl);
  58 PRIVATE int     try_moves(intermediate_t c, int maxE, intermediate_t *next, int dist);
  59
  60 /*
  61 #################################
  62 # BEGIN OF FUNCTION DEFINITIONS #
  63 #################################
  64 */
  65 PUBLIC void free_path(path_t *path){
  66   path_t *tmp = path;
  67   if(tmp){
  68     while(tmp->s){ free(tmp->s); tmp++;}
  69     free(path);
  70   }
  71 }
  72
  73 PUBLIC int find_saddle(const char *sequence, const char *struc1, const char *struc2, int max) {
  74   int maxl, maxE, i;
  75   const char *tmp;
  76   move_t *bestpath=NULL;
  77   int dir;
  78
  79   path_fwd = 0;
  80   maxE = INT_MAX - 1;
  81   seq = sequence;
  82
  83   update_fold_params();
  84   make_pair_matrix();
  85
  86   /* nummerically encode sequence */
  87   S   = encode_sequence(seq, 0);
  88   S1  = encode_sequence(seq, 1);
  89
  90   maxl=1;
  91   do {
  92     int saddleE;
  93     path_fwd = !path_fwd;
  94     if (maxl>max) maxl=max;
  95     if(path) free(path);
  96     saddleE  = find_path_once(struc1, struc2, maxE, maxl);
  97     if (saddleE<maxE) {
  98       maxE = saddleE;
  99       if (bestpath) free(bestpath);
 100       bestpath = path;
 101       path = NULL;
 102       dir = path_fwd;
 103     } else{
 104       free(path);path=NULL;
 105     }
 106     tmp=struc1;
 107     struc1=struc2;
 108     struc2=tmp;
 109     maxl *=2;
 110   } while (maxl<2*max);
 111
 112   free(S); free(S1);
 113   /* (re)set some globals */
 114   path=bestpath;
 115   path_fwd = dir;
 116   return maxE;
 117 }
 118
 119 PUBLIC void print_path(const char *seq, const char *struc) {
 120   int d;
 121   char *s;
 122   s = strdup(struc);
 123   printf("%s\n%s %6.2f\n", seq, s, energy_of_structure(seq,s, 0));
 124   qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
 125   for (d=0; d<BP_dist; d++) {
 126     int i,j;
 127     i = path[d].i; j=path[d].j;
 128     if (i<0) { /* delete */
 129       s[(-i)-1] = s[(-j)-1] = '.';
 130     } else {
 131       s[i-1] = '('; s[j-1] = ')';
 132     }
 133     printf("%s %6.2f - %6.2f\n", s, energy_of_structure(seq,s, 0), path[d].E/100.0);
 134   }
 135   free(s);
 136 }
 137
 138 PUBLIC path_t *get_path(const char *seq, const char *s1, const char* s2, int maxkeep) {
 139   int E, d;
 140   path_t *route=NULL;
 141
 142   E = find_saddle(seq, s1, s2, maxkeep);
 143
 144   route = (path_t *)space((BP_dist+2)*sizeof(path_t));
 145
 146   qsort(path, BP_dist, sizeof(move_t), compare_moves_when);
 147
 148   if (path_fwd) {
 149     /* memorize start of path */
 150     route[0].s  = strdup(s1);
 151     route[0].en = energy_of_structure(seq, s1, 0);
 152
 153     for (d=0; d<BP_dist; d++) {
 154       int i,j;
 155       route[d+1].s = strdup(route[d].s);
 156       i = path[d].i; j=path[d].j;
 157       if (i<0) { /* delete */
 158         route[d+1].s[(-i)-1] = route[d+1].s[(-j)-1] = '.';
 159       } else {
 160         route[d+1].s[i-1] = '('; route[d+1].s[j-1] = ')';
 161       }
 162       route[d+1].en = path[d].E/100.0;
 163     }
 164   }
 165   else {
 166     /* memorize start of path */
 167
 168     route[BP_dist].s  = strdup(s2);
 169     route[BP_dist].en = energy_of_structure(seq, s2, 0);
 170
 171     for (d=0; d<BP_dist; d++) {
 172       int i,j;
 173       route[BP_dist-d-1].s = strdup(route[BP_dist-d].s);
 174       i = path[d].i;
 175       j = path[d].j;
 176       if (i<0) { /* delete */
 177         route[BP_dist-d-1].s[(-i)-1] = route[BP_dist-d-1].s[(-j)-1] = '.';
 178       } else {
 179         route[BP_dist-d-1].s[i-1] = '('; route[BP_dist-d-1].s[j-1] = ')';
 180       }
 181       route[BP_dist-d-1].en = path[d].E/100.0;
 182     }
 183   }
 184
 185 #if _DEBUG_FINDPATH_
 186   fprintf(stderr, "\n%s\n%s\n%s\n\n", seq, s1, s2);
 187   for (d=0; d<=BP_dist; d++)
 188     fprintf(stderr, "%s %6.2f\n", route[d].s, route[d].en);
 189   fprintf(stderr, "%d\n", *num_entry);
 190 #endif
 191
 192   free(path);path=NULL;
 193   return (route);
 194 }
 195
 196 PRIVATE int try_moves(intermediate_t c, int maxE, intermediate_t *next, int dist) {
 197   int *loopidx, len, num_next=0, en, oldE;
 198   move_t *mv;
 199   short *pt;
 200
 201   len = c.pt[0];
 202   loopidx = pair_table_to_loop_index(c.pt);
 203   oldE = c.Sen;
 204   for (mv=c.moves; mv->i!=0; mv++) {
 205     int i,j;
 206     if (mv->when>0) continue;
 207     i = mv->i; j = mv->j;
 208     pt = (short *) space(sizeof(short)*(len+1));
 209     memcpy(pt, c.pt,(len+1)*sizeof(short));
 210     if (j<0) { /*it's a delete move */
 211       pt[-i]=0;
 212       pt[-j]=0;
 213     } else { /* insert move */
 214       if ((loopidx[i] == loopidx[j]) && /* i and j belong to same loop */
 215           (pt[i] == 0) && (pt[j]==0)     /* ... and are unpaired */
 216           ) {
 217         pt[i]=j;
 218         pt[j]=i;
 219       } else {
 220         free(pt);
 221         continue; /* llegal move, try next; */
 222       }
 223     }
 224 #ifdef LOOP_EN
 225     en = c.curr_en + energy_of_move_pt(c.pt, S, S1, i, j);
 226 #else
 227     en = energy_of_structure_pt(seq, pt, S, S1, 0);
 228 #endif
 229     if (en<maxE) {
 230       next[num_next].Sen = (en>oldE)?en:oldE;
 231       next[num_next].curr_en = en;
 232       next[num_next].pt = pt;
 233       mv->when=dist;
 234       mv->E = en;
 235       next[num_next++].moves = copy_moves(c.moves);
 236       mv->when=0;
 237     }
 238     else free(pt);
 239   }
 240   free(loopidx);
 241   return num_next;
 242 }
 243
 244 PRIVATE int find_path_once(const char *struc1, const char *struc2, int maxE, int maxl) {
 245   short *pt1, *pt2;
 246   move_t *mlist;
 247   int i, len, d, dist=0, result;
 248   intermediate_t *current, *next;
 249
 250   pt1 = make_pair_table(struc1);
 251   pt2 = make_pair_table(struc2);
 252   len = (int) strlen(struc1);
 253
 254   mlist = (move_t *) space(sizeof(move_t)*len); /* bp_dist < n */
 255
 256   for (i=1; i<=len; i++) {
 257     if (pt1[i] != pt2[i]) {
 258       if (i<pt1[i]) { /* need to delete this pair */
 259         mlist[dist].i = -i;
 260         mlist[dist].j = -pt1[i];
 261         mlist[dist++].when = 0;
 262       }
 263       if (i<pt2[i]) { /* need to insert this pair */
 264         mlist[dist].i = i;
 265         mlist[dist].j = pt2[i];
 266         mlist[dist++].when = 0;
 267       }
 268     }
 269   }
 270   free(pt2);
 271   BP_dist = dist;
 272   current = (intermediate_t *) space(sizeof(intermediate_t)*(maxl+1));
 273   current[0].pt = pt1;
 274   current[0].Sen = current[0].curr_en = energy_of_structure_pt(seq, pt1, S, S1, 0);
 275   current[0].moves = mlist;
 276   next = (intermediate_t *) space(sizeof(intermediate_t)*(dist*maxl+1));
 277
 278   for (d=1; d<=dist; d++) { /* go through the distance classes */
 279     int c, u, num_next=0;
 280     intermediate_t *cc;
 281
 282     for (c=0; current[c].pt != NULL; c++) {
 283       num_next += try_moves(current[c], maxE, next+num_next, d);
 284     }
 285     if (num_next==0) {
 286       for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
 287       current[0].Sen=INT_MAX;
 288       break;
 289     }
 290     /* remove duplicates via sort|uniq
 291        if this becomes a bottleneck we can use a hash instead */
 292     qsort(next, num_next, sizeof(intermediate_t),compare_ptable);
 293     for (u=0,c=1; c<num_next; c++) {
 294       if (memcmp(next[u].pt,next[c].pt,sizeof(short)*len)!=0) {
 295         next[++u] = next[c];
 296       } else {
 297         free_intermediate(next+c);
 298       }
 299     }
 300     num_next = u+1;
 301     qsort(next, num_next, sizeof(intermediate_t),compare_energy);
 302     /* free the old stuff */
 303     for (cc=current; cc->pt != NULL; cc++) free_intermediate(cc);
 304     for (u=0; u<maxl && u<num_next; u++) {
 305       current[u] = next[u];
 306     }
 307     for (; u<num_next; u++)
 308       free_intermediate(next+u);
 309     num_next=0;
 310   }
 311   free(next);
 312   path = current[0].moves;
 313   result = current[0].Sen;
 314   free(current[0].pt); free(current);
 315   return(result);
 316 }
 317
 318
 319 PRIVATE int *pair_table_to_loop_index (short *pt){
 320   /* number each position by which loop it belongs to (positions start
 321      at 1) */
 322   int i,hx,l,nl;
 323   int length;
 324   int *stack = NULL;
 325   int *loop = NULL;
 326
 327   length = pt[0];
 328   stack  = (int *) space(sizeof(int)*(length+1));
 329   loop   = (int *) space(sizeof(int)*(length+2));
 330   hx=l=nl=0;
 331
 332   for (i=1; i<=length; i++) {
 333     if ((pt[i] != 0) && (i < pt[i])) { /* ( */
 334       nl++; l=nl;
 335       stack[hx++]=i;
 336     }
 337     loop[i]=l;
 338
 339     if ((pt[i] != 0) && (i > pt[i])) { /* ) */
 340       --hx;
 341       if (hx>0)
 342         l = loop[stack[hx-1]];  /* index of enclosing loop   */
 343       else l=0;                 /* external loop has index 0 */
 344       if (hx<0) {
 345         nrerror("unbalanced brackets in make_pair_table");
 346       }
 347     }
 348   }
 349   loop[0] = nl;
 350   free(stack);
 351
 352 #ifdef _DEBUG_LOOPIDX
 353   fprintf(stderr,"begin loop index\n");
 354   fprintf(stderr,
 355           "....,....1....,....2....,....3....,....4"
 356           "....,....5....,....6....,....7....,....8\n");
 357   print_structure(pt, loop[0]);
 358   for (i=1; i<=length; i++)
 359     fprintf(stderr,"%2d ", loop[i]);
 360   fprintf(stderr,"\n");
 361   fprintf(stderr, "end loop index\n");
 362   fflush(stderr);
 363 #endif
 364
 365   return (loop);
 366 }
 367
 368 PRIVATE void free_intermediate(intermediate_t *i) {
 369    free(i->pt);
 370    free(i->moves);
 371    i->pt = NULL;
 372    i->moves = NULL;
 373    i->Sen = INT_MAX;
 374  }
 375
 376 PRIVATE int compare_ptable(const void *A, const void *B) {
 377   intermediate_t *a, *b;
 378   int c;
 379   a = (intermediate_t *) A;
 380   b = (intermediate_t *) B;
 381
 382   c = memcmp(a->pt, b->pt, a->pt[0]*sizeof(short));
 383   if (c!=0) return c;
 384   if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
 385   return (a->curr_en - b->curr_en);
 386 }
 387
 388 PRIVATE int compare_energy(const void *A, const void *B) {
 389   intermediate_t *a, *b;
 390   a = (intermediate_t *) A;
 391   b = (intermediate_t *) B;
 392
 393   if ((a->Sen - b->Sen) != 0) return (a->Sen - b->Sen);
 394   return (a->curr_en - b->curr_en);
 395 }
 396
 397 PRIVATE int compare_moves_when(const void *A, const void *B) {
 398   move_t *a, *b;
 399   a = (move_t *) A;
 400   b = (move_t *) B;
 401
 402   return(a->when - b->when);
 403 }
 404
 405 PRIVATE move_t* copy_moves(move_t *mvs) {
 406   move_t *new;
 407   new = (move_t *) space(sizeof(move_t)*(BP_dist+1));
 408   memcpy(new,mvs,sizeof(move_t)*(BP_dist+1));
 409   return new;
 410 }
 411
 412 #ifdef TEST_FINDPATH
 413
 414 int main(int argc, char *argv[]) {
 415   char *seq, *s1, *s2;
 416   int E, maxkeep=10;
 417   int verbose=0, i;
 418   path_t *route, *r;
 419
 420   for (i=1; i<argc; i++) {
 421     switch ( argv[i][1] ) {
 422     case 'm': if (strcmp(argv[i],"-m")==0)
 423         sscanf(argv[++i], "%d", &maxkeep);
 424       break;
 425     case 'v':  verbose = !strcmp(argv[i],"-v");
 426       break;
 427     case 'd': if (strcmp(argv[i],"-d")==0)
 428         sscanf(argv[++i], "%d", &dangles);
 429       break;
 430     }
 431   }
 432   seq = get_line(stdin);
 433   s1 = get_line(stdin);
 434   s2 = get_line(stdin);
 435
 436   E = find_saddle(seq, s1, s2, maxkeep);
 437   printf("saddle_energy = %6.2f\n", E/100.);
 438   if (verbose) {
 439     if (path_fwd)
 440       print_path(seq,s1);
 441     else
 442       print_path(seq,s2);
 443     free(path);
 444     path = NULL;
 445     route = get_path(seq, s1, s2, maxkeep);
 446     for (r=route; r->s; r++) {
 447       printf("%s %6.2f - %6.2f\n", r->s, energy_of_struct(seq,r->s), r->en);
 448       free(r->s);
 449     }
 450   }
 451   free(seq); free(s1); free(s2);
 452   return(EXIT_SUCCESS);
 453 }
 454 #endif