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