binaries/src/ViennaRNA/Kinfold/baum.c

   1 /*
   2   Last changed Time-stamp: <2010-06-24 15:12:09 ivo>
   3   c  Christoph Flamm and Ivo L Hofacker
   4   {xtof,ivo}@tbi.univie.ac.at
   5   Kinfold: $Name:  $
   6   $Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $
   7 */
   8
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <assert.h>
  13 #include "fold_vars.h"
  14 #include "fold.h"
  15 #include "utils.h"
  16 #include "pair_mat.h"
  17 #include "nachbar.h"
  18 #include "globals.h"
  19
  20 #define MYTURN 1
  21 #define SAME_STRAND(I,J) (((I)>=cut_point)||((J)<cut_point))
  22 #define ORDER(x,y) if ((x)->nummer>(y)->nummer) {tempb=x; x=y; y=tempb;}
  23
  24 /* item of structure ringlist */
  25 typedef struct _baum {
  26   int nummer; /* number of base in sequence */
  27   char typ;   /* 'r' virtualroot, 'p' or 'q' paired, 'u' unpaired */
  28   unsigned short base; /* 0<->unknown, 1<->A, 2<->C, 3<->G, 4<->U */
  29   int loop_energy;
  30   struct _baum *up;
  31   struct _baum *next;
  32   struct _baum *prev;
  33   struct _baum *down;
  34 } baum;
  35
  36 static char UNUSED rcsid[]="$Id: baum.c,v 1.9 2008/05/21 10:15:45 ivo Exp $";
  37 static short *pairList = NULL;
  38 static short *typeList = NULL;
  39 static short *aliasList = NULL;
  40 static baum *rl = NULL;         /* ringlist */
  41 static baum *wurzl = NULL;      /* virtualroot of ringlist-tree */
  42 static char **ptype = NULL;
  43
  44 static int comp_struc(const void *A, const void *B);
  45 /* PUBLIC FUNCTIONES */
  46 void ini_or_reset_rl (void);
  47 void move_it (void);
  48 void update_tree (int i, int j);
  49 void clean_up_rl (void);
  50
  51 /* PRIVATE FUNCTIONES */
  52 static void ini_ringlist(void);
  53 static void reset_ringlist(void);
  54 static void struc2tree (char *struc);
  55 static void close_bp_en (baum *i, baum *j);
  56 static void close_bp (baum *i, baum *j);
  57 static void open_bp (baum *i);
  58 static void open_bp_en (baum *i);
  59 static void inb (baum *root);
  60 static void inb_nolp (baum *root);
  61 static void dnb (baum *rli);
  62 static void dnb_nolp (baum *rli);
  63 static void fnb (baum *rli);
  64 static void make_ptypes(const short *S);
  65 /* debugging tool(s) */
  66 #if 0
  67 static void rl_status(void);
  68 #endif
  69
  70 /* convert structure in bracked-dot-notation to a ringlist-tree */
  71 static void struc2tree(char *struc) {
  72   char* struc_copy;
  73   int ipos, jpos, balance = 0;
  74   baum *rli, *rlj;
  75
  76   struc_copy = (char *)calloc(GSV.len+1, sizeof(char));
  77   assert(struc_copy);
  78   strcpy(struc_copy,struc);
  79
  80   for (ipos = 0; ipos < GSV.len; ipos++) {
  81     if (struc_copy[ipos] == ')') {
  82       jpos = ipos;
  83       struc_copy[ipos] = '.';
  84       balance++;
  85       while (struc_copy[--ipos] != '(');
  86       struc_copy[ipos] = '.';
  87       balance--;
  88       rli = &rl[ipos];
  89       rlj = &rl[jpos];
  90       close_bp(rli, rlj);
  91     }
  92   }
  93
  94   if (balance) {
  95     fprintf(stderr,
  96             "struc2tree(): start structure is not balanced !\n%s\n%s\n",
  97             GAV.farbe, struc);
  98     exit(1);
  99   }
 100
 101   GSV.currE = GSV.startE =
 102     (float )energy_of_structure_pt(GAV.farbe,
 103                                 pairList, typeList, aliasList, 0) / 100.0;
 104   {
 105     int i;
 106     for(i = 0; i < GSV.len; i++) {
 107       if (pairList[i+1]>i+1)
 108         rl[i].loop_energy = loop_energy(pairList, typeList, aliasList,i+1);
 109     }
 110     wurzl->loop_energy = loop_energy(pairList, typeList, aliasList,0);
 111   }
 112
 113   free(struc_copy);
 114 }
 115
 116 /**/
 117 static void ini_ringlist(void) {
 118   int i;
 119
 120   /* needed by function energy_of_struct_pt() from Vienna-RNA-1.4 */
 121   pairList = (short *)calloc(GSV.len + 2, sizeof(short));
 122   assert(pairList != NULL);
 123   typeList = (short *)calloc(GSV.len + 2, sizeof(short));
 124   assert(typeList != NULL);
 125   aliasList = (short *)calloc(GSV.len + 2, sizeof(short));
 126   assert(aliasList != NULL);
 127   pairList[0] = typeList[0] = aliasList[0] = GSV.len;
 128   ptype =  (char **)calloc(GSV.len + 2, sizeof(char *));
 129   assert(ptype != NULL);
 130   for (i=0; i<=GSV.len; i++) {
 131     ptype[i] =   (char*)calloc(GSV.len + 2, sizeof(char));
 132     assert(ptype[i] != NULL);
 133   }
 134
 135   /* allocate virtual root */
 136   wurzl = (baum *)calloc(1, sizeof(baum));
 137   assert(wurzl != NULL);
 138   /* allocate ringList */
 139   rl = (baum *)calloc(GSV.len+1, sizeof(baum));
 140   assert(rl != NULL);
 141   /* allocate PostOrderList */
 142
 143   /* initialize virtualroot */
 144   wurzl->typ = 'r';
 145   wurzl->nummer = -1;
 146   /* connect virtualroot to ringlist-tree in down direction */
 147   wurzl->down = &rl[GSV.len];
 148   /* initialize post-order list */
 149
 150   make_pair_matrix();
 151
 152   /* initialize rest of ringlist-tree */
 153   for(i = 0; i < GSV.len; i++) {
 154     int c;
 155     GAV.currform[i] = '.';
 156     GAV.prevform[i] = 'x';
 157     pairList[i+1] = 0;
 158     rl[i].typ = 'u';
 159     /* decode base to numeric value */
 160     c = encode_char(GAV.farbe[i]);
 161     rl[i].base = typeList[i+1] = c;
 162     aliasList[i+1] = alias[typeList[i+1]];
 163     /* astablish links for node of the ringlist-tree */
 164     rl[i].nummer = i;
 165     rl[i].next = &rl[i+1];
 166     rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i-1]);
 167     rl[i].up = rl[i].down = NULL;
 168   }
 169   GAV.currform[GSV.len] =   GAV.prevform[GSV.len] = '\0';
 170   make_ptypes(aliasList);
 171
 172   rl[i].nummer = i;
 173   rl[i].base = 0;
 174   /* make ringlist circular in next, prev direction */
 175   rl[i].next = &rl[0];
 176   rl[i].prev = &rl[i-1];
 177   /* make virtual basepair for virtualroot */
 178   rl[i].up = wurzl;
 179   rl[i].typ = 'x';
 180
 181 }
 182
 183 /**/
 184 void ini_or_reset_rl(void) {
 185
 186   /* if there is no ringList-tree make a new one */
 187   if (wurzl == NULL) {
 188     ini_ringlist();
 189
 190     /* start structure */
 191     struc2tree(GAV.startform);
 192     GSV.currE = GSV.startE = energy_of_structure(GAV.farbe, GAV.startform, 0);
 193
 194
 195     /* stop structure(s) */
 196     if ( GTV.stop )  {
 197       int i;
 198
 199       qsort(GAV.stopform, GSV.maxS, sizeof(char *), comp_struc);
 200       for (i = 0; i< GSV.maxS; i++)
 201         GAV.sE[i] = energy_of_structure(GAV.farbe_full, GAV.stopform[i], 0);
 202     }
 203     else {
 204       if(GTV.noLP)
 205         noLonelyPairs=1;
 206       initialize_cofold(GSV.len);
 207       /* fold sequence to get Minimum free energy structure (Mfe) */
 208       GAV.sE[0] = cofold(GAV.farbe_full, GAV.stopform[0]);
 209       free_arrays();
 210       /* revaluate energy of Mfe (maye differ if --logML=logarthmic */
 211       GAV.sE[0] = energy_of_structure(GAV.farbe_full, GAV.stopform[0], 0);
 212     }
 213     GSV.stopE = GAV.sE[0];
 214     ini_nbList(strlen(GAV.farbe_full)*strlen(GAV.farbe_full));
 215   }
 216   else {
 217     /* reset ringlist-tree to start conditions */
 218     reset_ringlist();
 219     if(GTV.start) struc2tree(GAV.startform);
 220     else {
 221       GSV.currE = GSV.startE;
 222     }
 223   }
 224 }
 225
 226 /**/
 227 static void reset_ringlist(void) {
 228   int i;
 229
 230   for(i = 0; i < GSV.len; i++) {
 231     GAV.currform[i] = '.';
 232     GAV.prevform[i] = 'x';
 233     pairList[i+1] = 0;
 234     rl[i].typ = 'u';
 235     rl[i].next = &rl[i + 1];
 236     rl[i].prev = ((i == 0) ? &rl[GSV.len] : &rl[i - 1]);
 237     rl[i].up = rl[i].down = NULL;
 238   }
 239   rl[i].next = &rl[0];
 240   rl[i].prev = &rl[i-1];
 241   rl[i].up = wurzl;
 242 }
 243
 244 /* update ringlist-tree */
 245 void update_tree(int i, int j) {
 246
 247   baum *rli, *rlj, *tempb;
 248
 249   if ( abs(i) < GSV.len) { /* >> single basepair move */
 250     if ((i > 0) && (j > 0)) { /* insert */
 251       rli = &rl[i-1];
 252       rlj = &rl[j-1];
 253       close_bp_en(rli, rlj);
 254     }
 255     else if ((i < 0)&&(j < 0)) { /* delete */
 256       i = -i;
 257       rli = &rl[i-1];
 258       open_bp_en(rli);
 259     }
 260     else { /* shift */
 261       if (i > 0) { /* i remains the same, j shifts */
 262         j=-j;
 263         rli=&rl[i-1];
 264         rlj=&rl[j-1];
 265         open_bp_en(rli);
 266         ORDER(rli, rlj);
 267         close_bp_en(rli, rlj);
 268       }
 269       else { /* j remains the same, i shifts */
 270         baum *old_rli;
 271         i = -i;
 272         rli = &rl[i-1];
 273         rlj = &rl[j-1];
 274         old_rli = rlj->up;
 275         open_bp_en(old_rli);
 276         ORDER(rli, rlj);
 277         close_bp_en(rli, rlj);
 278       }
 279     }
 280   } /* << single basepair move */
 281   else { /* >> double basepair move */
 282     if ((i > 0) && (j > 0)) { /* insert */
 283       rli = &rl[i-GSV.len-2];
 284       rlj = &rl[j-GSV.len-2];
 285       close_bp_en(rli->next, rlj->prev);
 286       close_bp_en(rli, rlj);
 287     }
 288     else if ((i < 0)&&(j < 0)) { /* delete */
 289       i = -i;
 290       rli = &rl[i-GSV.len-2];
 291       open_bp_en(rli);
 292       open_bp_en(rli->next);
 293     }
 294   } /* << double basepair move */
 295
 296 }
 297
 298 /* open a particular base pair */
 299 void open_bp(baum *i) {
 300
 301   baum *in; /* points to i->next */
 302
 303   /* change string representation */
 304   GAV.currform[i->nummer] = '.';
 305   GAV.currform[i->down->nummer] = '.';
 306
 307   /* change pairtable representation */
 308   pairList[1 + i->nummer] = 0;
 309   pairList[1 + i->down->nummer] = 0;
 310
 311   /* change tree representation */
 312   in = i->next;
 313   i->typ = 'u';
 314   i->down->typ = 'u';
 315   i->next = i->down->next;
 316   i->next->prev = i;
 317   in->prev = i->down;
 318   i->down->next = in;
 319   i->down = in->prev->up = NULL;
 320 }
 321
 322 /* close a particular base pair */
 323 void close_bp (baum *i, baum *j) {
 324
 325   baum *jn; /* points to j->next */
 326
 327   /* change string representation */
 328   GAV.currform[i->nummer] = '(';
 329   GAV.currform[j->nummer] = ')';
 330
 331   /* change pairtable representation */
 332   pairList[1 + i->nummer] = 1+ j->nummer;
 333   pairList[1 + j->nummer] = 1 + i->nummer;
 334
 335   /* change tree representation */
 336   jn = j->next;
 337   i->typ = 'p';
 338   j->typ = 'q';
 339   i->down = j;
 340   j->up = i;
 341   i->next->prev = j;
 342   j->next->prev = i;
 343   j->next = i->next;
 344   i->next = jn;
 345 }
 346
 347 # if 0
 348 /* for a given tree, generate postorder-list */
 349 static void make_poList (baum *root) {
 350
 351   baum *stop, *rli;
 352
 353   if (!root) root = wurzl;
 354   stop = root->down;
 355
 356   /* foreach base in ringlist ... */
 357   for (rli = stop->next; rli != stop; rli = rli->next) {
 358     /* ... explore subtee if bp found */
 359     if (rli->typ == 'p') {
 360       /*  fprintf(stderr, "%d >%d<\n", poListop, rli->nummer); */
 361       poList[poListop++] = rli;
 362       if ( poListop > GSV.len+1 ) {
 363         fprintf(stderr, "Something went wrong in make_poList()\n");
 364         exit(1);
 365       }
 366       make_poList(rli);
 367     }
 368   }
 369   return;
 370 }
 371 #endif
 372
 373 /* for a given ringlist, generate all structures
 374    with one additional basepair */
 375 static void inb(baum *root) {
 376
 377   int EoT;
 378   int E_old, E_new_in, E_new_out;
 379   baum *stop,*rli,*rlj;
 380
 381   E_old = root->loop_energy;
 382   stop=root->down;
 383   /* loop ringlist over all possible i positions */
 384   for(rli=stop->next;rli!=stop;rli=rli->next){
 385     /* potential i-position is already paired */
 386     if(rli->typ=='p') continue;
 387     /* loop ringlist over all possible j positions */
 388     for(rlj=rli->next;rlj!=stop;rlj=rlj->next){
 389       /* base pair must enclose at least 3 bases */
 390       if(rlj->nummer - rli->nummer < MYTURN) continue;
 391       /* potential j-position is already paired */
 392       if(rlj->typ=='p') continue;
 393       /* if i-j can form a base pair ... */
 394       if(ptype[rli->nummer][rlj->nummer]){
 395         /* close the base bair and ... */
 396         close_bp(rli,rlj);
 397         E_new_in  = loop_energy(pairList, typeList, aliasList,rli->nummer+1);
 398         E_new_out = loop_energy(pairList, typeList, aliasList,root->nummer+1);
 399         /* ... evaluate energy of the structure */
 400         EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) +  E_new_in + E_new_out - E_old ;
 401         /* assert(EoT ==  energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList)); */
 402         /* open the base pair again... */
 403         open_bp(rli);
 404         /* ... and put the move and the enegy
 405            of the structure into the neighbour list */
 406         update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
 407       }
 408     }
 409   }
 410 }
 411
 412 /* for a given ringlist, generate all structures (canonical)
 413    with one additional base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
 414 static void inb_nolp(baum *root) {
 415
 416   int EoT = 0;
 417   baum *stop, *rli, *rlj;
 418
 419   stop = root->down;
 420     /* loop ringlist over all possible i positions */
 421   for (rli=stop->next;rli!=stop;rli=rli->next) {
 422     /* potential i-position is already paired */
 423     if (rli->typ=='p') continue;
 424     /* loop ringlist over all possible j positions */
 425     for (rlj=rli->next;rlj!=stop;rlj=rlj->next) {
 426       /* base pair must enclose at least 3 bases */
 427       if (rlj->nummer - rli->nummer < MYTURN) continue;
 428       /* potential j-position is already paired */
 429       if (rlj->typ=='p') continue;
 430       /* if i-j can form a base pair ... */
 431       if (ptype[rli->nummer][rlj->nummer]) {
 432         /* ... and extends a helix ... */
 433         if (((rli->prev==stop && rlj->next==stop) && stop->typ != 'x') ||
 434             (rli->next == rlj->prev)) {
 435           /* ... close the base bair and ... */
 436           close_bp(rli,rlj);
 437           /* ... evaluate energy of the structure */
 438           EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 439           /* open the base pair again... */
 440           open_bp(rli);
 441           /* ... and put the move and the enegy
 442              of the structure into the neighbour list */
 443           update_nbList(1 + rli->nummer, 1 + rlj->nummer, EoT);
 444         }
 445         /* if double insertion is possible ... */
 446         else if ((rlj->nummer - rli->nummer >= MYTURN+2)&&
 447                  (rli->next->typ != 'p' && rlj->prev->typ != 'p') &&
 448                  (rli->next->next != rlj->prev->prev) &&
 449                  (ptype[rli->next->nummer][rlj->prev->nummer])) {
 450           /* close the two base bair and ... */
 451           close_bp(rli->next, rlj->prev);
 452           close_bp(rli, rlj);
 453           /* ... evaluate energy of the structure */
 454           EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 455           /* open the two base pair again ... */
 456           open_bp(rli);
 457           open_bp(rli->next);
 458           /* ... and put the move and the enegy
 459              of the structure into the neighbour list */
 460           update_nbList(1+rli->nummer+GSV.len+1, 1+rlj->nummer+GSV.len+1, EoT);
 461         }
 462       }
 463     }
 464   }
 465 }
 466
 467 /* for a given ringlist, generate all structures
 468  with one less base pair */
 469 static void dnb(baum *rli){
 470
 471   int EoT, E_old_in, E_old_out, E_new;
 472
 473   baum *rlj, *r;
 474
 475   rlj=rli->down;
 476   open_bp(rli);
 477   /* ... evaluate energy of the structure */
 478
 479   for (r=rli->next; r->up==NULL; r=r->next);
 480   E_old_in = rli->loop_energy;
 481   E_old_out = r->up->loop_energy;
 482   E_new = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
 483   EoT = (int) (GSV.currE*100 + ((GSV.currE<0)?-0.4:0.4)) -
 484     E_old_in - E_old_out + E_new;
 485
 486   /* assert(EoT== energy_of_struct_pt(GAV.farbe, pairList, typeList, aliasList));*/
 487   close_bp(rli,rlj);
 488   update_nbList(-(1 + rli->nummer), -(1 + rlj->nummer), EoT);
 489 }
 490
 491 /* for a given ringlist, generate all structures (canonical)
 492  with one less base pair (BUT WITHOUT ISOLATED BASE PAIRS) */
 493 static void dnb_nolp(baum *rli) {
 494
 495   int EoT = 0;
 496   baum *rlj;
 497   baum *rlin = NULL; /* pointers to following pair in helix, if any */
 498   baum *rljn = NULL;
 499   baum *rlip = NULL; /* pointers to preceding pair in helix, if any */
 500   baum *rljp = NULL;
 501
 502   rlj = rli->down;
 503
 504   /* immediate interior base pair ? */
 505   if (rlj->next == rlj->prev) {
 506     rlin = rlj->next;
 507     rljn = rlin->down;
 508   }
 509
 510   /* immediate exterior base pair and not virtualroot ? */
 511   if (rli->prev == rli->next && rli->next->typ != 'x') {
 512     rlip = rli->next->up;
 513     rljp = rli->next;
 514   }
 515
 516   /* double delete ? */
 517   if (rlip==NULL && rlin && rljn->next != rljn->prev ) {
 518     /* open the two base pairs ... */
 519     open_bp(rli);
 520     open_bp(rlin);
 521     /* ... evaluate energy of the structure ... */
 522     EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 523     /* ... and put the move and the enegy
 524        of the structure into the neighbour list ... */
 525     update_nbList(-(1+rli->nummer+GSV.len+1),-(1+rlj->nummer+GSV.len+1), EoT);
 526     /* ... and close the two base pairs again */
 527     close_bp(rlin, rljn);
 528     close_bp(rli, rlj);
 529   } else { /* single delete */
 530     /* the following will work only if boolean expr are shortcicuited */
 531     if (rlip==NULL || (rlip->prev == rlip->next && rlip->prev->typ != 'x'))
 532       if (rlin ==NULL || (rljn->next == rljn->prev)) {
 533         /* open the base pair ... */
 534         open_bp(rli);
 535         /* ... evaluate energy of the structure ... */
 536         EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 537         /* ... and put the move and the enegy
 538            of the structure into the neighbour list ... */
 539         update_nbList(-(1 + rli->nummer),-(1 + rlj->nummer), EoT);
 540         /* and close the base pair again */
 541         close_bp(rli, rlj);
 542       }
 543   }
 544 }
 545
 546 /* for a given ringlist, generate all structures
 547  with one shifted base pair */
 548 static void fnb(baum *rli) {
 549
 550   int EoT = 0, x;
 551   baum *rlj, *stop, *help_rli, *help_rlj;
 552
 553   stop = rli->down;
 554
 555   /* examin interior loop of bp(ij); (.......)
 556      i of j move                      ->   <- */
 557   for (rlj = stop->next; rlj != stop; rlj = rlj->next) {
 558     /* prevent shifting to paired position */
 559     if ((rlj->typ=='p')||(rlj->typ=='q')) continue;
 560     /* j-position of base pair shifts to k position (ij)->(ik) i<k<j */
 561     if ( (rlj->nummer-rli->nummer >= MYTURN)
 562          && (ptype[rli->nummer][rlj->nummer]) ) {
 563       /* open original basepair */
 564       open_bp(rli);
 565       /* close shifted version of original basepair */
 566       close_bp(rli, rlj);
 567       /* evaluate energy of the structure */
 568       EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 569       /* put the move and the enegy of the structure into the neighbour list */
 570       update_nbList(1+rli->nummer, -(1+rlj->nummer), EoT);
 571       /* open shifted basepair */
 572       open_bp(rli);
 573       /* restore original basepair */
 574       close_bp(rli, stop);
 575     }
 576     /* i-position of base pair shifts to position k (ij)->(kj) i<k<j */
 577     if ( (stop->nummer-rlj->nummer >= MYTURN)
 578          && (ptype[stop->nummer][rlj->nummer]) ) {
 579       /* open original basepair */
 580       open_bp(rli);
 581       /* close shifted version of original basepair */
 582       close_bp(rlj, stop);
 583       /* evaluate energy of the structure */
 584       EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 585       /* put the move and the enegy of the structure into the neighbour list */
 586       update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
 587       /* open shifted basepair */
 588       open_bp(rlj);
 589       /* restore original basepair */
 590       close_bp(rli, stop);
 591     }
 592   }
 593   /* examin exterior loop of bp(ij);   (.......)
 594      i or j moves                    <-         -> */
 595   for (rlj=rli->next;rlj!=rli;rlj=rlj->next) {
 596     if ((rlj->typ=='p') || (rlj->typ=='q' ) || (rlj->typ=='x')) continue;
 597     x=rlj->nummer-rli->nummer;
 598     if (x<0) x=-x;
 599     /* j-position of base pair shifts to position k */
 600     if ((x >= MYTURN) && (ptype[rli->nummer][rlj->nummer])) {
 601       if (rli->nummer<rlj->nummer) {
 602         help_rli=rli;
 603         help_rlj=rlj;
 604       }
 605       else {
 606         help_rli=rlj;
 607         help_rlj=rli;
 608       }
 609       /* open original basepair */
 610       open_bp(rli);
 611       /* close shifted version of original basepair */
 612       close_bp(help_rli,help_rlj);
 613       /* evaluate energy of the structure */
 614       EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 615       /* put the move and the enegy of the structure into the neighbour list */
 616       update_nbList(1 + rli->nummer, -(1 + rlj->nummer), EoT);
 617       /* open shifted base pair */
 618       open_bp(help_rli);
 619       /* restore original basepair */
 620       close_bp(rli,stop);
 621     }
 622     x = rlj->nummer-stop->nummer;
 623     if (x < 0) x = -x;
 624     /* i-position of base pair shifts to position k */
 625     if ((x >= MYTURN) && (ptype[stop->nummer][rlj->nummer])) {
 626       if (stop->nummer < rlj->nummer) {
 627         help_rli = stop;
 628         help_rlj = rlj;
 629       }
 630       else {
 631         help_rli = rlj;
 632         help_rlj = stop;
 633       }
 634       /* open original basepair */
 635       open_bp(rli);
 636        /* close shifted version of original basepair */
 637       close_bp(help_rli, help_rlj);
 638       /* evaluate energy of the structure */
 639       EoT = energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0);
 640       /* put the move and the enegy of the structure into the neighbour list */
 641       update_nbList(-(1 + rlj->nummer), 1 + stop->nummer, EoT);
 642       /* open shifted basepair */
 643       open_bp(help_rli);
 644       /* restore original basepair */
 645       close_bp(rli,stop);
 646     }
 647   }
 648 }
 649
 650 /* for a given tree (structure),
 651    generate all neighbours according to moveset */
 652 void move_it (void) {
 653   int i;
 654
 655   GSV.currE =
 656     energy_of_structure_pt(GAV.farbe, pairList, typeList, aliasList, 0)/100.;
 657
 658   if ( GTV.noLP ) { /* canonical neighbours only */
 659     inb_nolp(wurzl);
 660     for (i = 0; i < GSV.len; i++) {
 661
 662       if (pairList[i+1]>i+1) {
 663         inb_nolp(rl+i);      /* insert pair neighbours */
 664         dnb_nolp(rl+i);  /* delete pair neighbour */
 665       }
 666     }
 667   }
 668   else { /* all neighbours */
 669     inb(wurzl);
 670     for (i = 0; i < GSV.len; i++) {
 671
 672       if (pairList[i+1]>i+1) {
 673         inb(rl+i);       /* insert pair neighbours */
 674         dnb(rl+i);  /* delete pair neighbour */
 675         if ( GTV.noShift == 0 ) fnb(rl+i);
 676       }
 677     }
 678   }
 679 }
 680
 681
 682 /**/
 683 void clean_up_rl(void) {
 684   int i;
 685   free(pairList); pairList=NULL;
 686   free(typeList); typeList = NULL;
 687   free(aliasList); aliasList = NULL;
 688   free(rl); rl=NULL;
 689   free(wurzl);  wurzl=NULL;
 690   for (i=0; i<=GSV.len; i++)
 691     free(ptype[i]);
 692   free(ptype);
 693   ptype=NULL;
 694 }
 695
 696 /**/
 697 static int comp_struc(const void *A, const void *B) {
 698   int aE, bE;
 699   aE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)A)[0], 0));
 700   bE = (int)(100 * energy_of_structure(GAV.farbe_full, ((char **)B)[0], 0));
 701   return (aE-bE);
 702 }
 703
 704 #if 0
 705 /**/
 706 static void rl_status(void) {
 707
 708   int i;
 709
 710   printf("\n%s\n%s\n", GAV.farbe, GAV.currform);
 711   for (i=0; i <= GSV.len; i++) {
 712     printf("%2d %c %c %2d %2d %2d %2d\n",
 713            rl[i].nummer,
 714            i == GSV.len ? 'X': GAV.farbe[i],
 715            rl[i].typ,
 716            rl[i].up==NULL?0:(rl[i].up)->nummer,
 717            rl[i].down==NULL?0:(rl[i].down)->nummer,
 718            (rl[i].prev)->nummer,
 719            (rl[i].next)->nummer);
 720   }
 721   printf("---\n");
 722 }
 723 #endif
 724
 725 #define TURN 3
 726 static void make_ptypes(const short *S) {
 727   int n,i,j,k,l;
 728   n=S[0];
 729   for (k=1; k<n; k++)
 730     for (l=1; l<=2; l++) {
 731       int type,ntype=0,otype=0;
 732       i=k; j = i+l+TURN;
 733       if (!SAME_STRAND(i,j)) j=cut_point;
 734       if (j>n) continue;
 735       type = pair[S[i]][S[j]];
 736       while ((i>=1)&&(j<=n)) {
 737         if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
 738         if (noLonelyPairs && (!otype) && (!ntype))
 739           type = 0; /* i.j can only form isolated pairs */
 740         ptype[i-1][j-1] = ptype[j-1][i-1] = (char) type;
 741         otype =  type;
 742         type  = ntype;
 743         i--; j++;
 744       }
 745     }
 746 }
 747
 748 static void close_bp_en (baum *i, baum *j) {
 749   /* close bp and update energy */
 750   baum *r;
 751   close_bp(i,j);
 752   i->loop_energy = loop_energy(pairList,typeList,aliasList,i->nummer+1);
 753   for (r=i->next; r->up==NULL; r=r->next);
 754   r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
 755 };
 756
 757 static void open_bp_en (baum *i) {
 758   /* open bp and update energy */
 759   baum *r;
 760   i->loop_energy=0;
 761   open_bp(i);
 762   for (r=i->next; r->up==NULL; r=r->next);
 763   r->up->loop_energy = loop_energy(pairList,typeList,aliasList,r->up->nummer+1);
 764 };