binaries/src/ViennaRNA/lib/move_set.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <limits.h>
   6 #include <time.h>
   7
   8 #include "fold.h"
   9
  10 #include "move_set.h"
  11
  12 /* maximum degeneracy value - if degeneracy is greater than this, program segfaults*/
  13 #define MAX_DEGEN 100
  14 #define MINGAP 3
  15
  16
  17 #define space(a) malloc(a)
  18
  19 #define bool int
  20 #define true 1
  21 #define false 0
  22
  23 int compare(short *lhs, short *rhs)
  24 {
  25
  26   /* printf("%d ", (int)lhs[0]); */
  27
  28   int i=1;
  29   char l=0,r=0;
  30   while (i<=lhs[0]) {
  31     l = (lhs[i]==0?'.':(lhs[i]<lhs[lhs[i]]?'(':')'));
  32     r = (rhs[i]==0?'.':(rhs[i]<rhs[rhs[i]]?'(':')'));
  33     if (l != r) break;
  34     i++;
  35   }
  36
  37   return (i<=lhs[0] && l>r);
  38 }
  39
  40 int find_min(short *arr[MAX_DEGEN], int begin, int end) {
  41   short *min = arr[begin];
  42   short min_num = begin;
  43   int i;
  44
  45   for (i=begin+1; i<end; i++) {
  46     if (compare(arr[i], min)) {
  47       min = arr[i];
  48       min_num = i;
  49     }
  50   }
  51   return min_num;
  52 }
  53
  54 int equals(const short *first, const short *second)
  55 {
  56   int i=1;
  57   while (i<=first[0] && first[i]==second[i]) {
  58     i++;
  59   }
  60   if (i>first[0]) return 1;
  61   else return 0;
  62 }
  63
  64 void copy_arr(short *dest, short *src)
  65 {
  66   if (!src || !dest) {
  67     fprintf(stderr, "Empty pointer in copying\n");
  68     return;
  69   }
  70   memcpy(dest, src, sizeof(short)*(src[0]+1));
  71 }
  72
  73 short *allocopy(short *src)
  74 {
  75   short *res = (short*) space(sizeof(short)*(src[0]+1));
  76   copy_arr(res, src);
  77   return res;
  78 }
  79
  80 /* ############################## DECLARATION #####################################*/
  81 /* private functions & declarations*/
  82
  83 static int cnt_move = 0;
  84 int count_move() {return cnt_move;}
  85
  86 /*check if base is lone*/
  87 int lone_base(short *pt, int i);
  88
  89 /* can move be done on this structure? (move is in the Enc)*/
  90 int check_insert(struct_en *str, int i, int j);
  91
  92 /* internal struct with moves, sequence, degeneracy and options*/
  93 typedef struct _Encoded {
  94   /* sequence*/
  95   short *s0;
  96   short *s1;
  97
  98   char  *seq;
  99
 100   /* moves*/
 101   int   bp_left;
 102   int   bp_right;
 103   int   bp_left2;   /* if noLP is enabled (and for shift moves)*/
 104   int   bp_right2;
 105
 106   /* options*/
 107   int noLP;
 108   int verbose_lvl;
 109   int first;
 110   int shift;
 111
 112   /* degeneracy*/
 113   int begin_unpr;
 114   int begin_pr;
 115   int end_unpr;
 116   int end_pr;
 117   short *processed[MAX_DEGEN];
 118   short *unprocessed[MAX_DEGEN];
 119   int current_en;
 120
 121   /* moves in random (needs to be freed afterwards)*/
 122   int *moves_from;
 123   int *moves_to;
 124   int num_moves;
 125
 126   /* function for flooding */
 127   int (*funct) (struct_en*, struct_en*);
 128
 129
 130 } Encoded;
 131
 132 /* frees all things allocated by degeneracy...*/
 133 void free_degen(Encoded *Enc)
 134 {
 135   int i;
 136   for (i=Enc->begin_unpr; i<Enc->end_unpr; i++) {
 137     if (Enc->unprocessed[i]) {
 138       free(Enc->unprocessed[i]);
 139       Enc->unprocessed[i]=NULL;
 140     }
 141   }
 142   for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
 143     if (Enc->processed[i]) {
 144       free(Enc->processed[i]);
 145       Enc->processed[i]=NULL;
 146     }
 147   }
 148   Enc->begin_pr=0;
 149   Enc->begin_unpr=0;
 150   Enc->end_pr=0;
 151   Enc->end_unpr=0;
 152 }
 153
 154 /* ############################## IMPLEMENTATION #####################################*/
 155
 156 /* reads a line no matter how long*/
 157 char* my_getline(FILE *fp)
 158 {
 159   char s[512], *line, *cp;
 160   line = NULL;
 161   do {
 162     if(fgets(s, 512, fp) == NULL) break;
 163     cp = strchr(s, '\n');
 164     if(cp != NULL) *cp = '\0';
 165     if(line == NULL) line = (char *) calloc(strlen(s) + 1, sizeof(char));
 166     else line = (char *) realloc(line, strlen(s) + strlen(line) + 1);
 167     strcat (line, s);
 168   } while (cp == NULL);
 169   return (line);
 170 }
 171
 172 inline void do_move(short *pt, int bp_left, int bp_right)
 173 {
 174   /* delete*/
 175   if (bp_left<0) {
 176     pt[-bp_left]=0;
 177     pt[-bp_right]=0;
 178   } else { /* insert*/
 179     pt[bp_left]=bp_right;
 180     pt[bp_right]=bp_left;
 181   }
 182 }
 183
 184 /* done with all structures along the way to deepest*/
 185 int update_deepest(Encoded *Enc, struct_en *str, struct_en *min)
 186 {
 187   /* apply move + get its energy*/
 188   int tmp_en;
 189   tmp_en = str->energy + energy_of_move_pt(str->structure, Enc->s0, Enc->s1, Enc->bp_left, Enc->bp_right);
 190   do_move(str->structure, Enc->bp_left, Enc->bp_right);
 191   if (Enc->bp_left2 != 0) {
 192     tmp_en += energy_of_move_pt(str->structure, Enc->s0, Enc->s1, Enc->bp_left2, Enc->bp_right2);
 193     do_move(str->structure, Enc->bp_left2, Enc->bp_right2);
 194   }
 195   int last_en = str->energy;
 196   str->energy = tmp_en;
 197
 198
 199   /* use f_point if we have it */
 200   if (Enc->funct) {
 201     int end = Enc->funct(str, min);
 202
 203     /*  undo moves */
 204     if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
 205     do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
 206     str->energy = last_en;
 207     Enc->bp_left=0;
 208     Enc->bp_right=0;
 209     Enc->bp_left2=0;
 210     Enc->bp_right2=0;
 211
 212     return (end?1:0);
 213   }
 214
 215   if (Enc->verbose_lvl>1) { fprintf(stderr, "  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n", tmp_en); }
 216
 217   /* better deepest*/
 218   if (tmp_en < min->energy) {
 219     min->energy = tmp_en;
 220     copy_arr(min->structure, str->structure);
 221
 222     /* delete degeneracy*/
 223     free_degen(Enc);
 224
 225     /* undo moves*/
 226     if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
 227     do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
 228     str->energy = last_en;
 229     Enc->bp_left=0;
 230     Enc->bp_right=0;
 231     Enc->bp_left2=0;
 232     Enc->bp_right2=0;
 233     return 1;
 234   }
 235
 236   /* degeneracy*/
 237   if ((str->energy == min->energy) && (Enc->current_en == min->energy)) {
 238     int found = 0;
 239     int i;
 240     for (i=Enc->begin_pr; i<Enc->end_pr; i++) {
 241       if (equals(Enc->processed[i], str->structure)) {
 242         found = 1;
 243         break;
 244       }
 245     }
 246     for (i=Enc->begin_unpr; !found && i<Enc->end_pr; i++) {
 247       if (equals(Enc->unprocessed[i], str->structure)) {
 248         found = 1;
 249         break;
 250       }
 251     }
 252
 253     if (!found) {
 254       Enc->unprocessed[Enc->end_unpr]=allocopy(str->structure);
 255       Enc->end_unpr++;
 256     }
 257   }
 258
 259   /* undo moves*/
 260   if (Enc->bp_left2!=0) do_move(str->structure, -Enc->bp_left2, -Enc->bp_right2);
 261   do_move(str->structure, -Enc->bp_left, -Enc->bp_right);
 262   str->energy = last_en;
 263   Enc->bp_left=0;
 264   Enc->bp_right=0;
 265   Enc->bp_left2=0;
 266   Enc->bp_right2=0;
 267   return 0;
 268 }
 269
 270
 271 /* deletions move set*/
 272 int deletions(Encoded *Enc, struct_en *str, struct_en *minim)
 273 {
 274   int cnt = 0;
 275   short *pt = str->structure;
 276   int len = pt[0];
 277   int i;
 278
 279   for (i=1; i<=len; i++) {
 280     if (pt[i]>pt[pt[i]]) {  /* '('*/
 281       Enc->bp_left=-i;
 282       Enc->bp_right=-pt[i];
 283
 284       /*if nolp enabled, make (maybe) 2nd delete*/
 285       if (Enc->noLP) {
 286         int lone = -1;
 287         if (lone_base(pt, i-1)) lone=i-1;
 288         else if (lone_base(pt, i+1)) lone=i+1;
 289         else if (lone_base(pt, pt[i]-1)) lone=pt[i]-1;
 290         else if (lone_base(pt, pt[i]+1)) lone=pt[i]+1;
 291
 292         /* check*/
 293         if (lone != -1 && (pt[lone]==0 || pt[pt[lone]]==0)) {
 294           fprintf(stderr, "WARNING: pt[%d(or %d)]!=\'.\'", lone, pt[lone]);
 295         }
 296
 297         if (lone != -1) {
 298           Enc->bp_left2=-lone-1;
 299           Enc->bp_right2=-pt[lone]-1;
 300         }
 301         if (!lone_base(pt, pt[lone]-1) && !lone_base(pt, pt[lone]+1)) {
 302           cnt += update_deepest(Enc, str, minim);
 303           /* in case useFirst is on and structure is found, end*/
 304           if (Enc->first && cnt > 0) return cnt;
 305         }
 306       } else {  /* nolp not enabled*/
 307         cnt += update_deepest(Enc, str, minim);
 308         /* in case useFirst is on and structure is found, end*/
 309         if (Enc->first && cnt > 0) return cnt;
 310       }
 311     }
 312   }
 313   return cnt;
 314 }
 315
 316   /* compatible base pair?*/
 317 inline bool compat(char a, char b) {
 318   if (a=='A' && b=='U') return true;
 319   if (a=='C' && b=='G') return true;
 320   if (a=='G' && b=='U') return true;
 321   if (a=='U' && b=='A') return true;
 322   if (a=='G' && b=='C') return true;
 323   if (a=='U' && b=='G') return true;
 324   /* and with T's*/
 325   if (a=='A' && b=='T') return true;
 326   if (a=='T' && b=='A') return true;
 327   if (a=='G' && b=='T') return true;
 328   if (a=='T' && b=='G') return true;
 329   return false;
 330 }
 331
 332 /* try insert base pair (i,j)*/
 333 inline bool try_insert(const short *pt, const char *seq, int i, int j)
 334 {
 335   if (i<=0 || j<=0 || i>pt[0] || j>pt[0]) return false;
 336   return (j-i>MINGAP && pt[j]==0 && pt[i]==0 && compat(seq[i-1], seq[j-1]));
 337 }
 338
 339 /*  try insert base pair (i,j) */
 340 inline bool try_insert_seq(const char *seq, int i, int j)
 341 {
 342   if (i<=0 || j<=0) return false;
 343   return (j-i>MINGAP && compat(seq[i-1], seq[j-1]));
 344 }
 345
 346 /* insertions move set*/
 347 int insertions(Encoded *Enc, struct_en *str, struct_en *minim)
 348 {
 349   int cnt = 0;
 350   short *pt = str->structure;
 351   int len = pt[0];
 352   int i,j;
 353
 354   for (i=1; i<=len; i++) {
 355     if (pt[i]==0) {
 356       for (j=i+1; j<=len; j++) {
 357         /* end if found closing bracket*/
 358         if (pt[j]!=0 && pt[j]<j) break;  /*')'*/
 359         if (pt[j]!=0 && pt[j]>j) {       /*'('*/
 360           j = pt[j];
 361           continue;
 362         }
 363         /* if conditions are met, do insert*/
 364         if (try_insert(pt, Enc->seq, i, j)) {
 365           Enc->bp_left=i;
 366           Enc->bp_right=j;
 367
 368           if (Enc->noLP) {
 369             /* if lone bases occur, try inserting one another base*/
 370             if (lone_base(pt, i) || lone_base(pt, j)) {
 371               /* inside*/
 372               if (try_insert(pt, Enc->seq, i+1, j-1)) {
 373                 Enc->bp_left2=i+1;
 374                 Enc->bp_right2=j-1;
 375                 cnt += update_deepest(Enc, str, minim);
 376                 /* in case useFirst is on and structure is found, end*/
 377                 if (Enc->first && cnt > 0) return cnt;
 378               } else  /*outside*/
 379               if (try_insert(pt, Enc->seq, i-1, j+1)) {
 380                 Enc->bp_left2=i-1;
 381                 Enc->bp_right2=j+1;
 382                 cnt += update_deepest(Enc, str, minim);
 383                 /* in case useFirst is on and structure is found, end*/
 384                 if (Enc->first && cnt > 0) return cnt;
 385               }
 386             } else {
 387               cnt += update_deepest(Enc, str, minim);
 388               /* in case useFirst is on and structure is found, end*/
 389               if (Enc->first && cnt > 0) return cnt;
 390             }
 391           } else {
 392             cnt += update_deepest(Enc, str, minim);
 393             /* in case useFirst is on and structure is found, end*/
 394             if (Enc->first && cnt > 0) return cnt;
 395           }
 396         }
 397       }
 398     }
 399   }
 400   return cnt;
 401 }
 402
 403 /*shift move set*/
 404 int shifts(Encoded *Enc, struct_en *str, struct_en *minim)
 405 {
 406   int cnt = 0;
 407   int brack_num = 0;
 408   short *pt = str->structure;
 409   int len = pt[0];
 410   int i, k;
 411
 412   for (i=1; i<=len; i++) {
 413     if (pt[i]!=0 && pt[i]>i) {  /*'('*/
 414       int j=pt[i];
 415
 416       /* outer switch left*/
 417       if (Enc->verbose_lvl>1) fprintf(stderr, "%2d bracket %2d position, outer switch left\n", brack_num+1, i);
 418       for (k=i-1; k>0; k--) {
 419         if (pt[k]!=0 && pt[k]>k/*'('*/) break;
 420         if (pt[k]!=0 && pt[k]<k/*')'*/) {
 421           k = pt[k];
 422           continue;
 423         }
 424         /* checks*/
 425         if (pt[k]!=0) {
 426           fprintf(stderr, "WARNING: \'%c\'should be \'.\' at pos %d!\n", pt[k], k);
 427         }
 428
 429         /* switch (i,j) to (k,j)*/
 430         if (j-k>MINGAP && compat(Enc->seq[k-1], Enc->seq[j-1])) {
 431           Enc->bp_left=-i;
 432           Enc->bp_right=-j;
 433           Enc->bp_left2=k;
 434           Enc->bp_right2=j;
 435           cnt += update_deepest(Enc, str, minim);
 436           /* in case useFirst is on and structure is found, end*/
 437           if (Enc->first && cnt > 0) return cnt;
 438         }
 439
 440         /* switch (i,j) to (k,i)*/
 441         if (i-k>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
 442           Enc->bp_left=-i;
 443           Enc->bp_right=-j;
 444           Enc->bp_left2=k;
 445           Enc->bp_right2=i;
 446           cnt += update_deepest(Enc, str, minim);
 447           /* in case useFirst is on and structure is found, end*/
 448           if (Enc->first && cnt > 0) return cnt;
 449
 450         }
 451       }
 452
 453       /* outer switch right*/
 454       if (Enc->verbose_lvl>1) fprintf(stderr, "%2d bracket %2d position, outer switch right\n", brack_num+1, i);
 455       for (k=j+1; k<=len; k++) {
 456         if (pt[k]!=0 && pt[k]<k/*')'*/) break;
 457         if (pt[k]!=0 && pt[k]>k/*'('*/) {
 458           k = pt[k];
 459           continue;
 460         }
 461
 462         /* check*/
 463         if (pt[k]!=0) {
 464           fprintf(stderr, "WARNING: \'%c\'should be \'.\' at pos %d!\n", pt[k], k);
 465         }
 466         /* switch (i,j) to (i,k)*/
 467         if (k-i>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
 468           Enc->bp_left=-i;
 469           Enc->bp_right=-j;
 470           Enc->bp_left2=i;
 471           Enc->bp_right2=k;
 472           cnt += update_deepest(Enc, str, minim);
 473           /* in case useFirst is on and structure is found, end*/
 474           if (Enc->first && cnt > 0) return cnt;
 475         }
 476         /* switch (i,j) to (j,k)*/
 477         if (k-j>MINGAP && compat(Enc->seq[j-1], Enc->seq[k-1])) {
 478           Enc->bp_left=-i;
 479           Enc->bp_right=-j;
 480           Enc->bp_left2=j;
 481           Enc->bp_right2=k;
 482           cnt += update_deepest(Enc, str, minim);
 483           /* in case useFirst is on and structure is found, end*/
 484           if (Enc->first && cnt > 0) return cnt;
 485         }
 486       }
 487
 488       if (Enc->verbose_lvl>1) fprintf(stderr, "%2d bracket %2d position, inner switch\n", brack_num+1, i);
 489       /* inner switch*/
 490       for (k=i+1; k<j; k++) {
 491         /* jump to end of the sub-bracketing*/
 492         if (pt[k]!=0 && pt[k]>k/*'('*/) {
 493             k=pt[k];
 494             continue;
 495         }
 496
 497         /* left switch (i,j) to (k,j)*/
 498         if (j-k>MINGAP && compat(Enc->seq[k-1], Enc->seq[j-1])) {
 499           Enc->bp_left=-i;
 500           Enc->bp_right=-j;
 501           Enc->bp_left2=k;
 502           Enc->bp_right2=j;
 503           cnt += update_deepest(Enc, str, minim);
 504           /* in case useFirst is on and structure is found, end*/
 505           if (Enc->first && cnt > 0) return cnt;
 506         }
 507
 508         /* right switch (i,j) to (i,k)*/
 509         if (k-i>MINGAP && compat(Enc->seq[i-1], Enc->seq[k-1])) {
 510           Enc->bp_left=-i;
 511           Enc->bp_right=-j;
 512           Enc->bp_left2=i;
 513           Enc->bp_right2=k;
 514           cnt += update_deepest(Enc, str, minim);
 515           /* in case useFirst is on and structure is found, end*/
 516           if (Enc->first && cnt > 0) return cnt;
 517         }
 518       } /* end inner switch for*/
 519       brack_num++;
 520     } /* end if (pt[i]=='(')*/
 521   } /* end for in switches*/
 522   return cnt;
 523 }
 524
 525 /* move to deepest (or first) neighbour*/
 526 int move_set(Encoded *Enc, struct_en *str)
 527 {
 528   /* count how many times called*/
 529   cnt_move++;
 530
 531   /* count better neighbours*/
 532   int cnt = 0;
 533
 534   /* deepest descent*/
 535   struct_en min;
 536   min.structure = allocopy(str->structure);
 537   min.energy = str->energy;
 538   Enc->current_en = str->energy;
 539
 540   if (Enc->verbose_lvl>0) { fprintf(stderr, "  start of MS:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }
 541
 542   /* if using first dont do all of them*/
 543   bool end = false;
 544   /* insertions*/
 545   if (!end) cnt += insertions(Enc, str, &min);
 546   if (Enc->first && cnt>0) end = true;
 547   if (Enc->verbose_lvl>1) fprintf(stderr, "\n");
 548
 549   /* deletions*/
 550   if (!end) cnt += deletions(Enc, str, &min);
 551   if (Enc->first && cnt>0) end = true;
 552
 553   /* shifts (only if enabled + noLP disabled)*/
 554   if (!end && Enc->shift && !Enc->noLP) {
 555     cnt += shifts(Enc, str, &min);
 556     if (Enc->first && cnt>0) end = true;
 557   }
 558
 559   /* if degeneracy occurs, solve it!*/
 560   if (!end && (Enc->end_unpr - Enc->begin_unpr)>0) {
 561     Enc->processed[Enc->end_pr] = str->structure;
 562     Enc->end_pr++;
 563     str->structure = Enc->unprocessed[Enc->begin_unpr];
 564     Enc->unprocessed[Enc->begin_unpr]=NULL;
 565     Enc->begin_unpr++;
 566     cnt += move_set(Enc, str);
 567   } else {
 568     /* write output to str*/
 569     copy_arr(str->structure, min.structure);
 570     str->energy = min.energy;
 571   }
 572   /* release minimal*/
 573   free(min.structure);
 574
 575   /* resolve degeneracy in local minima*/
 576   if ((Enc->end_pr - Enc->begin_pr)>0) {
 577     Enc->processed[Enc->end_pr]=str->structure;
 578     Enc->end_pr++;
 579
 580     int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
 581     short *tmp = Enc->processed[min];
 582     Enc->processed[min] = Enc->processed[Enc->begin_pr];
 583     Enc->processed[Enc->begin_pr] = tmp;
 584     str->structure = Enc->processed[Enc->begin_pr];
 585     Enc->begin_pr++;
 586     free_degen(Enc);
 587   }
 588
 589   if (Enc->verbose_lvl>1 && !(Enc->first)) { fprintf(stderr, "\n  end of MS:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }
 590
 591   return cnt;
 592 }
 593
 594 void construct_moves(Encoded *Enc, short *structure)
 595 {
 596   /* generate all possible moves (less than n^2)*/
 597   Enc->num_moves = 0;
 598   int i;
 599   for (i=1; i<=structure[0]; i++) {
 600     if (structure[i]!=0) {
 601       if (structure[i]<i) continue;
 602       Enc->moves_from[Enc->num_moves]=-i;
 603       Enc->moves_to[Enc->num_moves]=-structure[i];
 604       Enc->num_moves++;
 605       /* fprintf(stderr, "add  d(%d, %d)\n", i, str.structure[i]); */
 606     } else {
 607       int j;
 608       for (j=i+1; j<=structure[0]; j++) {
 609         /* fprintf(stderr, "check (%d, %d)\n", i, j); */
 610         if (structure[j]==0) {
 611           if (try_insert_seq(Enc->seq,i,j)) {
 612             Enc->moves_from[Enc->num_moves]=i;
 613             Enc->moves_to[Enc->num_moves]=j;
 614             Enc->num_moves++;
 615             /* fprintf(stderr, "add  i(%d, %d)\n", i, j); */
 616             continue;
 617           }
 618         } else if (structure[j]>j) { /*  '(' */
 619           j = structure[j];
 620         } else break;
 621       }
 622     }
 623   }
 624
 625   /* permute them */
 626   for (i=0; i<Enc->num_moves-1; i++) {
 627     int rnd = rand();
 628     rnd = rnd % (Enc->num_moves-i) + i;
 629     int swp;
 630     swp = Enc->moves_from[i];
 631     Enc->moves_from[i]=Enc->moves_from[rnd];
 632     Enc->moves_from[rnd]=swp;
 633     swp = Enc->moves_to[i];
 634     Enc->moves_to[i]=Enc->moves_to[rnd];
 635     Enc->moves_to[rnd]=swp;
 636   }
 637 }
 638
 639 int move_rset(Encoded *Enc, struct_en *str)
 640 {
 641   /* count how many times called*/
 642   cnt_move++;
 643
 644   /* count better neighbours*/
 645   int cnt = 0;
 646
 647   /* deepest descent*/
 648   struct_en min;
 649   min.structure = allocopy(str->structure);
 650   min.energy = str->energy;
 651   Enc->current_en = str->energy;
 652
 653   if (Enc->verbose_lvl>0) { fprintf(stderr, "  start of MR:\n  "); print_str(stderr, str->structure); fprintf(stderr, " %d\n\n", str->energy); }
 654
 655   /*  construct and permute possible moves */
 656   construct_moves(Enc, str->structure);
 657
 658   /* find first lower one*/
 659   int i;
 660   for (i=0; i<Enc->num_moves; i++) {
 661     Enc->bp_left = Enc->moves_from[i];
 662     Enc->bp_right = Enc->moves_to[i];
 663     cnt = update_deepest(Enc, str, &min);
 664     if (cnt) break;
 665   }
 666
 667   /* if degeneracy occurs, solve it!*/
 668   if (!cnt && (Enc->end_unpr - Enc->begin_unpr)>0) {
 669     Enc->processed[Enc->end_pr] = str->structure;
 670     Enc->end_pr++;
 671     str->structure = Enc->unprocessed[Enc->begin_unpr];
 672     Enc->unprocessed[Enc->begin_unpr]=NULL;
 673     Enc->begin_unpr++;
 674     cnt += move_rset(Enc, str);
 675   } else {
 676     /* write output to str*/
 677     copy_arr(str->structure, min.structure);
 678     str->energy = min.energy;
 679   }
 680   /* release minimal*/
 681   free(min.structure);
 682
 683   /* resolve degeneracy in local minima*/
 684   if ((Enc->end_pr - Enc->begin_pr)>0) {
 685     Enc->processed[Enc->end_pr]=str->structure;
 686     Enc->end_pr++;
 687
 688     int min = find_min(Enc->processed, Enc->begin_pr, Enc->end_pr);
 689     short *tmp = Enc->processed[min];
 690     Enc->processed[min] = Enc->processed[Enc->begin_pr];
 691     Enc->processed[Enc->begin_pr] = tmp;
 692     str->structure = Enc->processed[Enc->begin_pr];
 693     Enc->begin_pr++;
 694     free_degen(Enc);
 695   }
 696
 697   return cnt;
 698 }
 699
 700 /*check if base is lone*/
 701 int lone_base(short *pt, int i)
 702 {
 703   if (i<=0 || i>pt[0]) return 0;
 704   /* is not a base pair*/
 705   if (pt[i]==0) return 0;
 706
 707   /* base is lone:*/
 708   if (i-1>0) {
 709     /* is base pair and is the same bracket*/
 710     if (pt[i-1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
 711   }
 712
 713   if (i+1<=pt[0]) {
 714     if (pt[i+1]!=0 && ((pt[i-1]<pt[pt[i-1]]) == (pt[i]<pt[pt[i]]))) return 0;
 715   }
 716
 717   return 1;
 718 }
 719
 720 /* if the structure has lone pairs*/
 721 int find_lone_pair(short* str)
 722 {
 723   int i;
 724   for(i=1; i<str[0]; i++) {
 725     if (str[i]==0) continue; /* '.'*/
 726
 727     if (str[i]>str[str[i]]) {  /* '('*/
 728       if (i+1==str[0] || str[i+1]==0 || str[i+1]<str[str[i+1]]) {
 729         return i;
 730       } else while (i+1!=str[0] && str[i+1]!=0 && str[i+1]>str[str[i+1]]) i++;
 731     }
 732
 733     if (str[i]<str[str[i]]) {  /* ')'*/
 734       if (i+1==str[0] || str[i+1]==0 || str[i+1]>str[str[i+1]]) {
 735         return i;
 736       } else while (i+1!=str[0] && str[i+1]!=0 && str[i+1]<str[str[i+1]]) i++;
 737     }
 738   }
 739
 740   return -1;
 741 }
 742
 743
 744 int move_deepest( char *string,
 745                   short *ptable,
 746                   short *s,
 747                   short *s1,
 748                   int verbosity_level,
 749                   int shifts,
 750                   int noLP)
 751 {
 752   cnt_move = 0;
 753
 754   Encoded enc;
 755   enc.seq = string;
 756   enc.s0 = s;
 757   enc.s1 = s1;
 758
 759   /* moves*/
 760   enc.bp_left=0;
 761   enc.bp_right=0;
 762   enc.bp_left2=0;
 763   enc.bp_right2=0;
 764
 765   /* options*/
 766   enc.noLP=noLP;
 767   enc.verbose_lvl=verbosity_level;
 768   enc.first=0;
 769   enc.shift=shifts;
 770
 771   /* degeneracy*/
 772   enc.begin_unpr=0;
 773   enc.begin_pr=0;
 774   enc.end_unpr=0;
 775   enc.end_pr=0;
 776   enc.current_en=0;
 777
 778   /*  function */
 779   enc.funct=NULL;
 780
 781   int i;
 782   for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
 783
 784   struct_en str;
 785   str.structure = allocopy(ptable);
 786   str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);
 787
 788   while (move_set(&enc, &str)!=0) {
 789     free_degen(&enc);
 790   }
 791   free_degen(&enc);
 792
 793   copy_arr(ptable, str.structure);
 794   free(str.structure);
 795
 796   return str.energy;
 797 }
 798
 799 int move_first(   char *string,
 800                   short *ptable,
 801                   short *s,
 802                   short *s1,
 803                   int verbosity_level,
 804                   int shifts,
 805                   int noLP)
 806 {
 807   cnt_move = 0;
 808
 809   Encoded enc;
 810   enc.seq = string;
 811   enc.s0 = s;
 812   enc.s1 = s1;
 813
 814   /* moves*/
 815   enc.bp_left=0;
 816   enc.bp_right=0;
 817   enc.bp_left2=0;
 818   enc.bp_right2=0;
 819
 820   /* options*/
 821   enc.noLP=noLP;
 822   enc.verbose_lvl=verbosity_level;
 823   enc.first=1;
 824   enc.shift=shifts;
 825
 826   /* degeneracy*/
 827   enc.begin_unpr=0;
 828   enc.begin_pr=0;
 829   enc.end_unpr=0;
 830   enc.end_pr=0;
 831   enc.current_en=0;
 832
 833   /*  function */
 834   enc.funct=NULL;
 835
 836   int i;
 837   for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
 838
 839   struct_en str;
 840   str.structure = allocopy(ptable);
 841   str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);
 842
 843   while (move_set(&enc, &str)!=0) {
 844     free_degen(&enc);
 845   }
 846   free_degen(&enc);
 847
 848   copy_arr(ptable, str.structure);
 849   free(str.structure);
 850
 851   return str.energy;
 852 }
 853
 854 int move_rand(   char *string,
 855                   short *ptable,
 856                   short *s,
 857                   short *s1,
 858                   int verbosity_level)
 859 {
 860   srand(time(NULL));
 861
 862   cnt_move = 0;
 863
 864   Encoded enc;
 865   enc.seq = string;
 866   enc.s0 = s;
 867   enc.s1 = s1;
 868
 869   /* moves*/
 870   enc.bp_left=0;
 871   enc.bp_right=0;
 872   enc.bp_left2=0;
 873   enc.bp_right2=0;
 874
 875   /* options*/
 876   enc.noLP=0;
 877   enc.verbose_lvl=verbosity_level;
 878   enc.first=1;
 879   enc.shift=0;
 880
 881   /* degeneracy*/
 882   enc.begin_unpr=0;
 883   enc.begin_pr=0;
 884   enc.end_unpr=0;
 885   enc.end_pr=0;
 886   enc.current_en=0;
 887
 888   /*  function */
 889   enc.funct=NULL;
 890
 891   /*  allocate memory for moves */
 892   enc.moves_from = (int*) space(ptable[0]*ptable[0]*sizeof(int));
 893   enc.moves_to = (int*) space(ptable[0]*ptable[0]*sizeof(int));
 894
 895   int i;
 896   for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
 897
 898   struct_en str;
 899   str.structure = allocopy(ptable);
 900   str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);
 901
 902   while (move_rset(&enc, &str)!=0) {
 903     free_degen(&enc);
 904   }
 905   free_degen(&enc);
 906
 907   copy_arr(ptable, str.structure);
 908   free(str.structure);
 909   free(enc.moves_from);
 910   free(enc.moves_to);
 911
 912   return str.energy;
 913 }
 914
 915 int browse_neighs(   char *string,
 916                   short *ptable,
 917                   short *s,
 918                   short *s1,
 919                   int verbosity_level,
 920                   int shifts,
 921                   int noLP,
 922                   int (*funct) (struct_en*, struct_en*))
 923 {
 924   cnt_move = 0;
 925
 926   Encoded enc;
 927   enc.seq = string;
 928   enc.s0 = s;
 929   enc.s1 = s1;
 930
 931   /* moves*/
 932   enc.bp_left=0;
 933   enc.bp_right=0;
 934   enc.bp_left2=0;
 935   enc.bp_right2=0;
 936
 937   /* options*/
 938   enc.noLP=noLP;
 939   enc.verbose_lvl=verbosity_level;
 940   enc.first=1;
 941   enc.shift=shifts;
 942
 943   /* degeneracy*/
 944   enc.begin_unpr=0;
 945   enc.begin_pr=0;
 946   enc.end_unpr=0;
 947   enc.end_pr=0;
 948   enc.current_en=0;
 949
 950   /*  function */
 951   enc.funct=funct;
 952
 953   int i;
 954   for (i=0; i<MAX_DEGEN; i++) enc.processed[i]=enc.unprocessed[i]=NULL;
 955
 956   struct_en str;
 957   str.structure = allocopy(ptable);
 958   str.energy = energy_of_structure_pt(enc.seq, str.structure, enc.s0, enc.s1, 0);
 959
 960   move_set(&enc, &str);
 961   free_degen(&enc);
 962
 963   copy_arr(ptable, str.structure);
 964   free(str.structure);
 965
 966   return str.energy;
 967 }
 968
 969 /* printf*/
 970 void print_stren(FILE *out, struct_en *str) {
 971   print_str(out, str->structure);
 972   fprintf(out, " %6.2f\n", str->energy/100.0);
 973 }
 974
 975 void print_str(FILE *out, short *str) {
 976   int i;
 977   for (i=1; i<=str[0]; i++) {
 978     if (str[i]==0) fprintf(out, ".");
 979     else if (str[i]<i) fprintf(out, ")");
 980     else fprintf(out, "(");
 981   }
 982 }
 983
 984