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