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