1 /*
2 * suboptimal folding - Stefan Wuchty, Walter Fontana & Ivo Hofacker
3 *
4 * Vienna RNA package
5 */
6
7 #ifdef HAVE_CONFIG_H
8 #include "config.h"
9 #endif
10
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <unistd.h>
14 #include <ctype.h>
15 #include <string.h>
16 #include <math.h>
17 #include "ViennaRNA/fold.h"
18 #include "ViennaRNA/constraints/hard.h"
19 #include "ViennaRNA/constraints/soft.h"
20 #include "ViennaRNA/utils/basic.h"
21 #include "ViennaRNA/utils/strings.h"
22 #include "ViennaRNA/params/default.h"
23 #include "ViennaRNA/fold_vars.h"
24 #include "ViennaRNA/datastructures/lists.h"
25 #include "ViennaRNA/eval.h"
26 #include "ViennaRNA/params/basic.h"
27 #include "ViennaRNA/loops/all.h"
28 #include "ViennaRNA/cofold.h"
29 #include "ViennaRNA/gquad.h"
30 #include "ViennaRNA/alphabet.h"
31 #include "ViennaRNA/subopt.h"
32
33 /* hack */
34 #include "ViennaRNA/color_output.inc"
35
36 #ifdef _OPENMP
37 #include <omp.h>
38 #endif
39
40 #define true 1
41 #define false 0
42
43 #ifndef ON_SAME_STRAND
44 #define ON_SAME_STRAND(I, J, C) (((I) >= (C)) || ((J) < (C)))
45 #endif
46
47 /**
48 * @brief Sequence interval stack element used in subopt.c
49 */
50 typedef struct INTERVAL {
51 int i;
52 int j;
53 int array_flag;
54 } INTERVAL;
55
56 typedef struct {
57 char *structure;
58 LIST *Intervals;
59 int partial_energy;
60 int is_duplex;
61 /* int best_energy; */ /* best attainable energy */
62 } STATE;
63
64 typedef struct {
65 LIST *Intervals;
66 LIST *Stack;
67 int nopush;
68 } subopt_env;
69
70
71 struct old_subopt_dat {
72 unsigned long max_sol;
73 unsigned long n_sol;
74 SOLUTION *SolutionList;
75 FILE *fp;
76 int cp;
77 };
78
79 /*
80 #################################
81 # GLOBAL VARIABLES #
82 #################################
83 */
84 PUBLIC int subopt_sorted = 0; /* output sorted by energy */
85 PUBLIC int density_of_states[MAXDOS + 1];
86 PUBLIC double print_energy = 9999; /* printing threshold for use with logML */
87
88 /*
89 #################################
90 # PRIVATE VARIABLES #
91 #################################
92 */
93
94 /* some backward compatibility stuff */
95 PRIVATE int backward_compat = 0;
96 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
97
98 #ifdef _OPENMP
99
100 #pragma omp threadprivate(backward_compat_compound, backward_compat)
101
102 #endif
103
104 /*
105 #################################
106 # PRIVATE FUNCTION DECLARATIONS #
107 #################################
108 */
109
110 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
111
112 PRIVATE SOLUTION *
113 wrap_subopt(char *seq,
114 char *structure,
115 vrna_param_t *parameters,
116 int delta,
117 int is_constrained,
118 int is_circular,
119 FILE *fp);
120
121
122 #endif
123
124 PRIVATE void
125 make_pair(int i,
126 int j,
127 STATE *state);
128
129
130 /* mark a gquadruplex in the resulting dot-bracket structure */
131 PRIVATE void
132 make_gquad(int i,
133 int L,
134 int l[3],
135 STATE *state);
136
137
138 PRIVATE INTERVAL *
139 make_interval(int i,
140 int j,
141 int ml);
142
143
144 PRIVATE STATE *
145 make_state(LIST *Intervals,
146 char *structure,
147 int partial_energy,
148 int is_duplex,
149 int length);
150
151
152 PRIVATE STATE *
153 copy_state(STATE *state);
154
155
156 PRIVATE void
157 print_state(STATE *state);
158
159
160 PRIVATE void
161 UNUSED print_stack(LIST *list);
162
163
164 PRIVATE LIST *
165 make_list(void);
166
167
168 PRIVATE void
169 push(LIST *list,
170 void *data);
171
172
173 PRIVATE void
174 *pop(LIST *list);
175
176
177 PRIVATE int
178 best_attainable_energy(vrna_fold_compound_t *vc,
179 STATE *state);
180
181
182 PRIVATE void
183 scan_interval(vrna_fold_compound_t *vc,
184 int i,
185 int j,
186 int array_flag,
187 int threshold,
188 STATE *state,
189 subopt_env *env);
190
191
192 PRIVATE void
193 free_interval_node(INTERVAL *node);
194
195
196 PRIVATE void
197 free_state_node(STATE *node);
198
199
200 PRIVATE void
201 push_back(LIST *Stack,
202 STATE *state);
203
204
205 PRIVATE char *
206 get_structure(STATE *state);
207
208
209 PRIVATE int
210 compare(const void *solution1,
211 const void *solution2);
212
213 PRIVATE int
214 compare_en(const void *solution1,
215 const void *solution2);
216
217
218 PRIVATE void
219 make_output(SOLUTION *SL,
220 int cp,
221 FILE *fp);
222
223
224 PRIVATE void
225 repeat(vrna_fold_compound_t *vc,
226 int i,
227 int j,
228 STATE *state,
229 int part_energy,
230 int temp_energy,
231 int best_energy,
232 int threshold,
233 subopt_env *env);
234
235
236 PRIVATE void
237 repeat_gquad(vrna_fold_compound_t *vc,
238 int i,
239 int j,
240 STATE *state,
241 int part_energy,
242 int temp_energy,
243 int best_energy,
244 int threshold,
245 subopt_env *env);
246
247
248 PRIVATE void
249 old_subopt_print(const char *structure,
250 float energy,
251 void *data);
252
253
254 PRIVATE void
255 old_subopt_store(const char *structure,
256 float energy,
257 void *data);
258
259
260 PRIVATE void
261 old_subopt_store_compressed(const char *structure,
262 float energy,
263 void *data);
264
265
266 /*
267 #################################
268 # BEGIN OF FUNCTION DEFINITIONS #
269 #################################
270 */
271
272
273 /*---------------------------------------------------------------------------*/
274 /*List routines--------------------------------------------------------------*/
275 /*---------------------------------------------------------------------------*/
276
277 PRIVATE void
make_pair(int i,int j,STATE * state)278 make_pair(int i,
279 int j,
280 STATE *state)
281 {
282 state->structure[i - 1] = '(';
283 state->structure[j - 1] = ')';
284 }
285
286
287 PRIVATE void
make_gquad(int i,int L,int l[3],STATE * state)288 make_gquad(int i,
289 int L,
290 int l[3],
291 STATE *state)
292 {
293 int x;
294
295 for (x = 0; x < L; x++) {
296 state->structure[i - 1 + x] = '+';
297 state->structure[i - 1 + x + L + l[0]] = '+';
298 state->structure[i - 1 + x + 2 * L + l[0] + l[1]] = '+';
299 state->structure[i - 1 + x + 3 * L + l[0] + l[1] + l[2]] = '+';
300 }
301 }
302
303
304 /*---------------------------------------------------------------------------*/
305
306 PRIVATE INTERVAL *
make_interval(int i,int j,int array_flag)307 make_interval(int i,
308 int j,
309 int array_flag)
310 {
311 INTERVAL *interval;
312
313 interval = lst_newnode(sizeof(INTERVAL));
314 interval->i = i;
315 interval->j = j;
316 interval->array_flag = array_flag;
317 return interval;
318 }
319
320
321 /*---------------------------------------------------------------------------*/
322
323 PRIVATE void
free_interval_node(INTERVAL * node)324 free_interval_node(INTERVAL *node)
325 {
326 lst_freenode(node);
327 }
328
329
330 /*---------------------------------------------------------------------------*/
331
332 PRIVATE void
free_state_node(STATE * node)333 free_state_node(STATE *node)
334 {
335 free(node->structure);
336 if (node->Intervals)
337 lst_kill(node->Intervals, lst_freenode);
338
339 lst_freenode(node);
340 }
341
342
343 /*---------------------------------------------------------------------------*/
344
345 PRIVATE STATE *
make_state(LIST * Intervals,char * structure,int partial_energy,int is_duplex,int length)346 make_state(LIST *Intervals,
347 char *structure,
348 int partial_energy,
349 int is_duplex,
350 int length)
351 {
352 STATE *state;
353
354 state = lst_newnode(sizeof(STATE));
355
356 if (Intervals)
357 state->Intervals = Intervals;
358 else
359 state->Intervals = lst_init();
360
361 if (structure) {
362 state->structure = structure;
363 } else {
364 int i;
365 state->structure = (char *)vrna_alloc(length + 1);
366 for (i = 0; i < length; i++)
367 state->structure[i] = '.';
368 }
369
370 state->partial_energy = partial_energy;
371
372 return state;
373 }
374
375
376 /*---------------------------------------------------------------------------*/
377
378 PRIVATE STATE *
copy_state(STATE * state)379 copy_state(STATE *state)
380 {
381 STATE *new_state;
382 void *after;
383 INTERVAL *new_interval, *next;
384
385 new_state = lst_newnode(sizeof(STATE));
386 new_state->Intervals = lst_init();
387 new_state->partial_energy = state->partial_energy;
388 /* new_state->best_energy = state->best_energy; */
389
390 if (state->Intervals->count) {
391 after = LST_HEAD(new_state->Intervals);
392 for (next = lst_first(state->Intervals); next; next = lst_next(next)) {
393 new_interval = lst_newnode(sizeof(INTERVAL));
394 *new_interval = *next;
395 lst_insertafter(new_state->Intervals, new_interval, after);
396 after = new_interval;
397 }
398 }
399
400 new_state->structure = strdup(state->structure);
401 if (!new_state->structure)
402 vrna_message_error("out of memory");
403
404 return new_state;
405 }
406
407
408 /*---------------------------------------------------------------------------*/
409
410 /*@unused @*/ PRIVATE void
print_state(STATE * state)411 print_state(STATE *state)
412 {
413 INTERVAL *next;
414
415 if (state->Intervals->count) {
416 printf("%d intervals:\n", state->Intervals->count);
417 for (next = lst_first(state->Intervals); next; next = lst_next(next))
418 printf("[%d,%d],%d ", next->i, next->j, next->array_flag);
419 printf("\n");
420 }
421
422 printf("partial structure: %s\n", state->structure);
423 printf("\n");
424 printf(" partial_energy: %d\n", state->partial_energy);
425 /* printf(" best_energy: %d\n", state->best_energy); */
426 (void)fflush(stdout);
427 }
428
429
430 /*---------------------------------------------------------------------------*/
431
432 /*@unused @*/ PRIVATE void
print_stack(LIST * list)433 print_stack(LIST *list)
434 {
435 void *rec;
436
437 printf("================\n");
438 printf("%d states\n", list->count);
439 for (rec = lst_first(list); rec; rec = lst_next(rec)) {
440 printf("state-----------\n");
441 print_state(rec);
442 }
443 printf("================\n");
444 }
445
446
447 /*---------------------------------------------------------------------------*/
448
449 PRIVATE LIST *
make_list(void)450 make_list(void)
451 {
452 return lst_init();
453 }
454
455
456 /*---------------------------------------------------------------------------*/
457
458 PRIVATE void
push(LIST * list,void * data)459 push(LIST *list,
460 void *data)
461 {
462 lst_insertafter(list, data, LST_HEAD(list));
463 }
464
465
466 /* PRIVATE void */
467 /* push_stack(STATE *state) { */ /* keep the stack sorted by energy */
468 /* STATE *after, *next; */
469 /* nopush = false; */
470 /* next = after = LST_HEAD(Stack); */
471 /* while ( next = lst_next(next)) { */
472 /* if ( next->best_energy >= state->best_energy ) break; */
473 /* after = next; */
474 /* } */
475 /* lst_insertafter(Stack, state, after); */
476 /* } */
477
478 /*---------------------------------------------------------------------------*/
479
480 PRIVATE void *
pop(LIST * list)481 pop(LIST *list)
482 {
483 void *data;
484
485 data = lst_deletenext(list, LST_HEAD(list));
486 return data;
487 }
488
489
490 /*---------------------------------------------------------------------------*/
491 /*auxiliary routines---------------------------------------------------------*/
492 /*---------------------------------------------------------------------------*/
493
494 PRIVATE int
best_attainable_energy(vrna_fold_compound_t * vc,STATE * state)495 best_attainable_energy(vrna_fold_compound_t *vc,
496 STATE *state)
497 {
498 /* evaluation of best possible energy attainable within remaining intervals */
499
500 register int sum;
501 INTERVAL *next;
502 vrna_md_t *md;
503 vrna_mx_mfe_t *matrices;
504 int *indx;
505
506 md = &(vc->params->model_details);
507 matrices = vc->matrices;
508 indx = vc->jindx;
509
510 sum = state->partial_energy; /* energy of already found elements */
511
512 for (next = lst_first(state->Intervals); next; next = lst_next(next)) {
513 if (next->array_flag == 0)
514 sum += (md->circ) ? matrices->Fc : matrices->f5[next->j];
515 else if (next->array_flag == 1)
516 sum += matrices->fML[indx[next->j] + next->i];
517 else if (next->array_flag == 2)
518 sum += matrices->c[indx[next->j] + next->i];
519 else if (next->array_flag == 3)
520 sum += matrices->fM1[indx[next->j] + next->i];
521 else if (next->array_flag == 4)
522 sum += matrices->fc[next->i];
523 else if (next->array_flag == 5)
524 sum += matrices->fc[next->j];
525 else if (next->array_flag == 6)
526 sum += matrices->ggg[indx[next->j] + next->i];
527 }
528
529 return sum;
530 }
531
532
533 /*---------------------------------------------------------------------------*/
534
535 PRIVATE void
push_back(LIST * Stack,STATE * state)536 push_back(LIST *Stack,
537 STATE *state)
538 {
539 push(Stack, copy_state(state));
540 return;
541 }
542
543
544 /*---------------------------------------------------------------------------*/
545
546 PRIVATE char *
get_structure(STATE * state)547 get_structure(STATE *state)
548 {
549 char *structure;
550
551 structure = strdup(state->structure);
552 return structure;
553 }
554
555
556 /*---------------------------------------------------------------------------*/
557 PRIVATE int
compare(const void * solution1,const void * solution2)558 compare(const void *solution1,
559 const void *solution2)
560 {
561 if (((SOLUTION *)solution1)->energy > ((SOLUTION *)solution2)->energy)
562 return 1;
563
564 if (((SOLUTION *)solution1)->energy < ((SOLUTION *)solution2)->energy)
565 return -1;
566
567 return strcmp(((SOLUTION *)solution1)->structure,
568 ((SOLUTION *)solution2)->structure);
569 }
570
571
572 PRIVATE int
compare_en(const void * solution1,const void * solution2)573 compare_en(const void *solution1,
574 const void *solution2)
575 {
576 if (((SOLUTION *)solution1)->energy > ((SOLUTION *)solution2)->energy)
577 return 1;
578
579 if (((SOLUTION *)solution1)->energy < ((SOLUTION *)solution2)->energy)
580 return -1;
581
582 return 0;
583 }
584
585
586 /*---------------------------------------------------------------------------*/
587
588 PRIVATE void
make_output(SOLUTION * SL,int cp,FILE * fp)589 make_output(SOLUTION *SL,
590 int cp,
591 FILE *fp) /* prints stuff */
592 {
593 SOLUTION *sol;
594
595 for (sol = SL; sol->structure != NULL; sol++) {
596 char *e_string = vrna_strdup_printf(" %6.2f", sol->energy);
597 char *ss = vrna_db_unpack(sol->structure);
598 char *s = vrna_cut_point_insert(ss, cp);
599 print_structure(fp, s, e_string);
600 free(s);
601 free(ss);
602 free(e_string);
603 }
604 }
605
606
607 PRIVATE STATE *
derive_new_state(int i,int j,STATE * s,int e,int flag)608 derive_new_state(int i,
609 int j,
610 STATE *s,
611 int e,
612 int flag)
613 {
614 STATE *s_new = copy_state(s);
615 INTERVAL *ival = make_interval(i, j, flag);
616
617 push(s_new->Intervals, ival);
618
619 s_new->partial_energy += e;
620
621 return s_new;
622 }
623
624
625 PRIVATE void
fork_state(int i,int j,STATE * s,int e,int flag,subopt_env * env)626 fork_state(int i,
627 int j,
628 STATE *s,
629 int e,
630 int flag,
631 subopt_env *env)
632 {
633 STATE *s_new = derive_new_state(i, j, s, e, flag);
634
635 push(env->Stack, s_new);
636 env->nopush = false;
637 }
638
639
640 PRIVATE void
fork_int_state(int i,int j,int p,int q,STATE * s,int e,subopt_env * env)641 fork_int_state(int i,
642 int j,
643 int p,
644 int q,
645 STATE *s,
646 int e,
647 subopt_env *env)
648 {
649 STATE *s_new = derive_new_state(p, q, s, e, 2);
650
651 make_pair(i, j, s_new);
652 make_pair(p, q, s_new);
653 push(env->Stack, s_new);
654 env->nopush = false;
655 }
656
657
658 PRIVATE void
fork_state_pair(int i,int j,STATE * s,int e,subopt_env * env)659 fork_state_pair(int i,
660 int j,
661 STATE *s,
662 int e,
663 subopt_env *env)
664 {
665 STATE *new_state;
666
667 new_state = copy_state(s);
668 make_pair(i, j, new_state);
669 new_state->partial_energy += e;
670 push(env->Stack, new_state);
671 env->nopush = false;
672 }
673
674
675 PRIVATE void
fork_two_states_pair(int i,int j,int k,STATE * s,int e,int flag1,int flag2,subopt_env * env)676 fork_two_states_pair(int i,
677 int j,
678 int k,
679 STATE *s,
680 int e,
681 int flag1,
682 int flag2,
683 subopt_env *env)
684 {
685 INTERVAL *interval1, *interval2;
686 STATE *new_state;
687
688 new_state = copy_state(s);
689 interval1 = make_interval(i + 1, k - 1, flag1);
690 interval2 = make_interval(k, j - 1, flag2);
691 if (k - i < j - k) {
692 /* push larger interval first */
693 push(new_state->Intervals, interval1);
694 push(new_state->Intervals, interval2);
695 } else {
696 push(new_state->Intervals, interval2);
697 push(new_state->Intervals, interval1);
698 }
699
700 make_pair(i, j, new_state);
701 new_state->partial_energy += e;
702
703 push(env->Stack, new_state);
704 env->nopush = false;
705 }
706
707
708 PRIVATE void
fork_two_states(int i,int j,int p,int q,STATE * s,int e,int flag1,int flag2,subopt_env * env)709 fork_two_states(int i,
710 int j,
711 int p,
712 int q,
713 STATE *s,
714 int e,
715 int flag1,
716 int flag2,
717 subopt_env *env)
718 {
719 INTERVAL *interval1, *interval2;
720 STATE *new_state;
721
722 new_state = copy_state(s);
723 interval1 = make_interval(i, j, flag1);
724 interval2 = make_interval(p, q, flag2);
725
726 if ((j - i) < (q - p)) {
727 push(new_state->Intervals, interval1);
728 push(new_state->Intervals, interval2);
729 } else {
730 push(new_state->Intervals, interval2);
731 push(new_state->Intervals, interval1);
732 }
733
734 new_state->partial_energy += e;
735
736 push(env->Stack, new_state);
737 env->nopush = false;
738 }
739
740
741 /*---------------------------------------------------------------------------*/
742 /* start of subopt backtracking ---------------------------------------------*/
743 /*---------------------------------------------------------------------------*/
744
745 PUBLIC SOLUTION *
vrna_subopt(vrna_fold_compound_t * vc,int delta,int sorted,FILE * fp)746 vrna_subopt(vrna_fold_compound_t *vc,
747 int delta,
748 int sorted,
749 FILE *fp)
750 {
751 struct old_subopt_dat data;
752 vrna_subopt_callback *cb;
753
754 data.SolutionList = NULL;
755 data.max_sol = 128;
756 data.n_sol = 0;
757 data.fp = fp;
758 data.cp = vc->cutpoint;
759
760 if (vc) {
761 /* SolutionList stores the suboptimal structures found */
762
763 data.SolutionList = (SOLUTION *)vrna_alloc(data.max_sol * sizeof(SOLUTION));
764
765 /* end initialize ------------------------------------------------------- */
766
767 if (fp) {
768 float min_en;
769 char *SeQ, *energies = NULL;
770 if (vc->strands > 1)
771 min_en = vrna_mfe_dimer(vc, NULL);
772 else
773 min_en = vrna_mfe(vc, NULL);
774
775 SeQ = vrna_cut_point_insert(vc->sequence, vc->cutpoint);
776 energies = vrna_strdup_printf(" %6.2f %6.2f", min_en, (float)delta / 100.);
777 print_structure(fp, SeQ, energies);
778 free(SeQ);
779 free(energies);
780
781 vrna_mx_mfe_free(vc);
782 }
783
784 cb = old_subopt_store;
785
786 if (fp)
787 cb = (sorted) ? old_subopt_store_compressed : old_subopt_print;
788
789 /* call subopt() */
790 vrna_subopt_cb(vc, delta, cb, (void *)&data);
791
792 if (sorted) {
793 /* sort structures by energy */
794 if (data.n_sol > 0) {
795 int (*compare_fun)(const void *a, const void *b);
796
797 switch (sorted) {
798 case VRNA_SORT_BY_ENERGY_ASC:
799 compare_fun = compare_en;
800 break;
801
802 default: /* a.k.a. VRNA_SORT_BY_ENERGY_LEXICOGRAPHIC_ASC */
803 compare_fun = compare;
804 break;
805 }
806
807 qsort(data.SolutionList, data.n_sol - 1, sizeof(SOLUTION), compare_fun);
808 }
809
810 if (fp)
811 make_output(data.SolutionList, vc->cutpoint, fp);
812 }
813
814 if (fp) {
815 /* we've printed everything -- free solutions */
816 SOLUTION *sol;
817 for (sol = data.SolutionList; sol->structure != NULL; sol++)
818 free(sol->structure);
819 free(data.SolutionList);
820 data.SolutionList = NULL;
821 }
822 }
823
824 return data.SolutionList;
825 }
826
827
828 PUBLIC void
vrna_subopt_cb(vrna_fold_compound_t * vc,int delta,vrna_subopt_callback * cb,void * data)829 vrna_subopt_cb(vrna_fold_compound_t *vc,
830 int delta,
831 vrna_subopt_callback *cb,
832 void *data)
833 {
834 subopt_env *env;
835 STATE *state;
836 INTERVAL *interval;
837 unsigned int *so, *ss, *se;
838 int maxlevel, count, partial_energy, old_dangles, logML, dangle_model, length, circular,
839 threshold;
840 double structure_energy, min_en, eprint;
841 char *struc, *structure;
842 float correction;
843 vrna_param_t *P;
844 vrna_md_t *md;
845 int minimal_energy;
846 int Fc;
847 int *f5;
848
849 vrna_fold_compound_prepare(vc, VRNA_OPTION_MFE | VRNA_OPTION_HYBRID);
850
851 length = vc->length;
852 so = vc->strand_order;
853 ss = vc->strand_start;
854 se = vc->strand_end;
855 P = vc->params;
856 md = &(P->model_details);
857
858 /* do mfe folding to get fill arrays and get ground state energy */
859 /* in case dangles is neither 0 or 2, set dangles=2 while folding */
860
861 circular = md->circ;
862 logML = md->logML;
863 old_dangles = dangle_model = md->dangles;
864
865 if (md->uniq_ML != 1) /* failsafe mechanism to enforce valid fM1 array */
866 md->uniq_ML = 1;
867
868 /* temporarily set dangles to 2 if necessary */
869 if ((md->dangles != 0) && (md->dangles != 2))
870 md->dangles = 2;
871
872 struc = (char *)vrna_alloc(sizeof(char) * (length + 1));
873
874 if (circular) {
875 min_en = vrna_mfe(vc, struc);
876 Fc = vc->matrices->Fc;
877 f5 = vc->matrices->f5;
878
879 /* restore dangle model */
880 md->dangles = old_dangles;
881 /* re-evaluate in case we're using logML etc */
882 min_en = vrna_eval_structure(vc, struc);
883 } else {
884 min_en = vrna_mfe_dimer(vc, struc);
885
886 f5 = vc->matrices->f5;
887
888 /* restore dangle model */
889 md->dangles = old_dangles;
890 /* re-evaluate in case we're using logML etc */
891 min_en = vrna_eval_structure(vc, struc);
892 }
893
894 free(struc);
895 eprint = print_energy + min_en;
896
897 correction = (min_en < 0) ? -0.1 : 0.1;
898
899 /* Initialize ------------------------------------------------------------ */
900
901 maxlevel = 0;
902 count = 0;
903 partial_energy = 0;
904
905 /* Initialize the stack ------------------------------------------------- */
906
907 minimal_energy = (circular) ? Fc : f5[length];
908 threshold = minimal_energy + delta;
909 if (threshold >= INF) {
910 vrna_message_warning("Energy range too high, limiting to reasonable value");
911 threshold = INF - EMAX;
912 }
913
914 /* init env data structure */
915 env = (subopt_env *)vrna_alloc(sizeof(subopt_env));
916 env->Stack = NULL;
917 env->nopush = true;
918 env->Stack = make_list(); /* anchor */
919 env->Intervals = make_list(); /* initial state: */
920 interval = make_interval(1, length, 0); /* interval [1,length,0] */
921 push(env->Intervals, interval);
922 env->nopush = false;
923 state = make_state(env->Intervals, NULL, partial_energy, 0, length);
924 /* state->best_energy = minimal_energy; */
925 push(env->Stack, state);
926 env->nopush = false;
927
928 /* end initialize ------------------------------------------------------- */
929
930
931 while (1) {
932 /* forever, til nothing remains on stack */
933
934 maxlevel = (env->Stack->count > maxlevel ? env->Stack->count : maxlevel);
935
936 if (LST_EMPTY(env->Stack)) {
937 /* we are done! clean up and quit */
938 /* fprintf(stderr, "maxlevel: %d\n", maxlevel); */
939
940 lst_kill(env->Stack, free_state_node);
941
942 cb(NULL, 0, data); /* NULL (last time to call callback function */
943
944 break;
945 }
946
947 /* pop the last element ---------------------------------------------- */
948
949 state = pop(env->Stack); /* current state to work with */
950
951 if (LST_EMPTY(state->Intervals)) {
952 int e;
953 /* state has no intervals left: we got a solution */
954
955 count++;
956 structure = get_structure(state);
957 structure_energy = state->partial_energy / 100.;
958
959 #ifdef CHECK_ENERGY
960 structure_energy = vrna_eval_structure(vc, structure);
961
962 if (!logML)
963 if ((double)(state->partial_energy / 100.) != structure_energy) {
964 vrna_message_error("%s %6.2f %6.2f",
965 structure,
966 state->partial_energy / 100.,
967 structure_energy);
968 exit(1);
969 }
970
971 #endif
972 if (logML || (dangle_model == 1) || (dangle_model == 3)) /* recalc energy */
973 structure_energy = vrna_eval_structure(vc, structure);
974
975 e = (int)((structure_energy - min_en) * 10. - correction); /* avoid rounding errors */
976 if (e > MAXDOS)
977 e = MAXDOS;
978
979 density_of_states[e]++;
980 if (structure_energy <= eprint) {
981 char *outstruct = vrna_cut_point_insert(structure, (vc->strands > 1) ? ss[so[1]] : -1);
982 cb((const char *)outstruct, structure_energy, data);
983 free(outstruct);
984 }
985
986 free(structure);
987 } else {
988 /* get (and remove) next interval of state to analyze */
989
990 interval = pop(state->Intervals);
991 scan_interval(vc, interval->i, interval->j, interval->array_flag, threshold, state, env);
992
993 free_interval_node(interval); /* free the current interval */
994 }
995
996 free_state_node(state); /* free the current state */
997 } /* end of while (1) */
998
999 /* cleanup memory */
1000 free(env);
1001 }
1002
1003
1004 PRIVATE void
scan_interval(vrna_fold_compound_t * vc,int i,int j,int array_flag,int threshold,STATE * state,subopt_env * env)1005 scan_interval(vrna_fold_compound_t *vc,
1006 int i,
1007 int j,
1008 int array_flag,
1009 int threshold,
1010 STATE *state,
1011 subopt_env *env)
1012 {
1013 /* real backtrack routine */
1014
1015 /* array_flag = 0: trace back in f5-array */
1016 /* array_flag = 1: trace back in fML-array */
1017 /* array_flag = 2: trace back in repeat() */
1018 /* array_flag = 3: trace back in fM1-array */
1019
1020 STATE *new_state, *temp_state;
1021 INTERVAL *new_interval;
1022 vrna_param_t *P;
1023 vrna_md_t *md;
1024 register int k, fi, cij, ij;
1025 register int type;
1026 register int dangle_model;
1027 register int noLP;
1028 unsigned int *sn, *so, *ss, *se;
1029 int element_energy, best_energy;
1030 int *fc, *f5, *c, *fML, *fM1, *ggg;
1031 int FcH, FcI, FcM, *fM2;
1032 int length, *indx, *rtype, circular, with_gquad, turn;
1033 char *ptype;
1034 short *S1;
1035 unsigned char *hard_constraints, hc_decompose;
1036 vrna_hc_t *hc;
1037 vrna_sc_t *sc;
1038
1039 length = vc->length;
1040 sn = vc->strand_number;
1041 so = vc->strand_order;
1042 ss = vc->strand_start;
1043 se = vc->strand_end;
1044 indx = vc->jindx;
1045 ptype = vc->ptype;
1046 S1 = vc->sequence_encoding;
1047 P = vc->params;
1048 md = &(P->model_details);
1049 rtype = &(md->rtype[0]);
1050
1051 dangle_model = md->dangles;
1052 noLP = md->noLP;
1053 circular = md->circ;
1054 with_gquad = md->gquad;
1055 turn = md->min_loop_size;
1056
1057 fc = vc->matrices->fc;
1058 f5 = vc->matrices->f5;
1059 c = vc->matrices->c;
1060 fML = vc->matrices->fML;
1061 fM1 = vc->matrices->fM1;
1062 ggg = vc->matrices->ggg;
1063 FcH = vc->matrices->FcH;
1064 FcI = vc->matrices->FcI;
1065 FcM = vc->matrices->FcM;
1066 fM2 = vc->matrices->fM2;
1067
1068 hc = vc->hc;
1069 hard_constraints = hc->mx;
1070
1071 sc = vc->sc;
1072
1073 best_energy = best_attainable_energy(vc, state); /* .. on remaining intervals */
1074 env->nopush = true;
1075
1076 if ((i > 1) && (!array_flag))
1077 vrna_message_error("Error while backtracking!");
1078
1079 if ((j < i + turn + 1) &&
1080 ((sn[i] == so[1]) || (sn[j] == so[0]))) {
1081 /* minimal structure element */
1082 if (array_flag == 0)
1083 /* do not forget to add f5[j], since it may contain pseudo energies from soft constraining */
1084 state->partial_energy += f5[j];
1085
1086 if (env->nopush) {
1087 push_back(env->Stack, state);
1088 env->nopush = false;
1089 }
1090
1091 return;
1092 }
1093
1094 ij = indx[j] + i;
1095
1096 /* 13131313131313131313131313131313131313131313131313131313131313131313131 */
1097 if (array_flag == 3 || array_flag == 1) {
1098 /* array_flag = 3: interval i,j was generated during */
1099 /* a multiloop decomposition using array fM1 in repeat() */
1100 /* or in this block */
1101
1102 /* array_flag = 1: interval i,j was generated from a */
1103 /* stack, bulge, or internal loop in repeat() */
1104 /* or in this block */
1105
1106 if ((hc->up_ml[j]) &&
1107 (((array_flag == 3) && (fM1[indx[j - 1] + i] != INF)) ||
1108 (fML[indx[j - 1] + i] != INF))) {
1109
1110 if (array_flag == 3)
1111 fi = fM1[indx[j - 1] + i] + P->MLbase;
1112 else
1113 fi = fML[indx[j - 1] + i] + P->MLbase;
1114
1115 if (sc) {
1116 if (sc->energy_up)
1117 fi += sc->energy_up[j][1];
1118
1119 if (sc->f)
1120 fi += sc->f(i, j, i, j - 1, VRNA_DECOMP_ML_ML, sc->data);
1121 }
1122
1123 if ((fi + best_energy <= threshold) && (sn[j - 1] == sn[j])) {
1124 /* no basepair, nibbling of 3'-end */
1125 if ((sc) &&
1126 (sc->energy_up))
1127 fork_state(i, j - 1, state, P->MLbase + sc->energy_up[j][1], array_flag, env);
1128 else
1129 fork_state(i, j - 1, state, P->MLbase, array_flag, env);
1130 }
1131 }
1132
1133 hc_decompose = hard_constraints[length * i + j];
1134
1135 if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) {
1136 /* i,j may pair */
1137 cij = c[ij];
1138 if (cij != INF) {
1139 type = vrna_get_ptype(ij, ptype);
1140
1141 switch (dangle_model) {
1142 case 0:
1143 element_energy = E_MLstem(type, -1, -1, P);
1144 break;
1145 default:
1146 element_energy = E_MLstem(type,
1147 (((i > 1) && (sn[i - 1] == sn[i])) || circular) ? S1[i - 1] : -1,
1148 (((j < length) && (sn[j] == sn[j + 1])) || circular) ? S1[j + 1] : -1,
1149 P);
1150 break;
1151 }
1152
1153 if (sc) {
1154 if (sc->f)
1155 element_energy += sc->f(i, j, i, j, VRNA_DECOMP_ML_STEM, sc->data);
1156 }
1157
1158 cij += element_energy;
1159
1160 if (cij + best_energy <= threshold)
1161 repeat(vc, i, j, state, element_energy, 0, best_energy, threshold, env);
1162 }
1163 } else if ((with_gquad) && (ggg[ij] != INF)) {
1164 element_energy = E_MLstem(0, -1, -1, P);
1165 cij = ggg[ij] + element_energy;
1166
1167 if (cij + best_energy <= threshold)
1168 repeat_gquad(vc, i, j, state, element_energy, 0, best_energy, threshold, env);
1169 }
1170 } /* array_flag == 3 || array_flag == 1 */
1171
1172 /* 11111111111111111111111111111111111111111111111111111111111111111111111 */
1173
1174 if (array_flag == 1) {
1175 /* array_flag = 1: interval i,j was generated from a */
1176 /* stack, bulge, or internal loop in repeat() */
1177 /* or in this block */
1178
1179 int stopp, k1j;
1180 if ((sn[i - 1] == sn[i]) && (sn[j] == sn[j + 1])) {
1181 /*backtrack in FML only if multiloop is possible*/
1182 for (k = i + turn + 1; k <= j - 1 - turn; k++) {
1183 /* Multiloop decomposition if i,j contains more than 1 stack */
1184
1185 if ((with_gquad) &&
1186 (sn[k] == sn[k + 1]) &&
1187 (fML[indx[k] + i] != INF) &&
1188 (ggg[indx[j] + k + 1] != INF)) {
1189 element_energy = E_MLstem(0, -1, -1, P);
1190
1191 if (fML[indx[k] + i] + ggg[indx[j] + k + 1] + element_energy + best_energy <= threshold) {
1192 temp_state = derive_new_state(i, k, state, 0, array_flag);
1193 env->nopush = false;
1194 repeat_gquad(vc,
1195 k + 1,
1196 j,
1197 temp_state,
1198 element_energy,
1199 fML[indx[k] + i],
1200 best_energy,
1201 threshold,
1202 env);
1203 free_state_node(temp_state);
1204 }
1205 }
1206
1207 k1j = indx[j] + k + 1;
1208
1209 if ((hard_constraints[length * j + k + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) &&
1210 (fML[indx[k] + i] != INF) &&
1211 (c[k1j] != INF)) {
1212 short s5, s3;
1213 type = vrna_get_ptype(k1j, ptype);
1214
1215 switch (dangle_model) {
1216 case 0:
1217 s5 = s3 = -1;
1218 break;
1219
1220 default:
1221 s5 = (sn[i - 1] == sn[i]) ? S1[k] : -1;
1222 s3 = (sn[j] == sn[j + 1]) ? S1[j + 1] : -1;
1223 break;
1224 }
1225
1226 element_energy = E_MLstem(type, s5, s3, P);
1227
1228 if (sc) {
1229 if (sc->f)
1230 element_energy += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_STEM, sc->data);
1231 }
1232
1233 if (sn[k] == sn[k + 1]) {
1234 if (fML[indx[k] + i] + c[k1j] + element_energy + best_energy <= threshold) {
1235 temp_state = derive_new_state(i, k, state, 0, array_flag);
1236 env->nopush = false;
1237 repeat(vc,
1238 k + 1,
1239 j,
1240 temp_state,
1241 element_energy,
1242 fML[indx[k] + i],
1243 best_energy,
1244 threshold,
1245 env);
1246 free_state_node(temp_state);
1247 }
1248 }
1249 }
1250 }
1251 }
1252
1253 if (vc->strands > 1) {
1254 stopp = se[so[0]] - 1; /*if cp -1: k on cut, => no ml*/
1255 stopp = MIN2(stopp, j - 1 - turn);
1256 if (i > ss[so[1]])
1257 stopp = j - 1 - turn;
1258 else if (i == ss[so[1]])
1259 stopp = 0; /*not a multi loop*/
1260 } else {
1261 stopp = j - 1 - turn;
1262 }
1263 int up = 1;
1264 for (k = i; k <= stopp; k++, up++) {
1265 if (hc->up_ml[i] >= up) {
1266 k1j = indx[j] + k + 1;
1267
1268 /* Multiloop decomposition if i,j contains only 1 stack */
1269 if ((with_gquad) && (ggg[k1j] != INF)) {
1270 element_energy = E_MLstem(0, -1, -1, P) + P->MLbase * up;
1271
1272 if (sc)
1273 if (sc->energy_up)
1274 element_energy += sc->energy_up[i][up];
1275
1276 if (ggg[k1j] + element_energy + best_energy <= threshold)
1277 repeat_gquad(vc, k + 1, j, state, element_energy, 0, best_energy, threshold, env);
1278 }
1279
1280 if ((hard_constraints[length * j + k + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) &&
1281 (c[k1j] != INF)) {
1282 int s5, s3;
1283
1284 type = vrna_get_ptype(k1j, ptype);
1285
1286 switch (dangle_model) {
1287 case 0:
1288 s5 = s3 = -1;
1289 break;
1290 default:
1291 s5 = (sn[k - 1] == sn[k]) ? S1[k] : -1;
1292 s3 = (sn[j] == sn[j + 1]) ? S1[j + 1] : -1;
1293 break;
1294 }
1295
1296 element_energy = E_MLstem(type, s5, s3, P);
1297
1298 element_energy += P->MLbase * up;
1299
1300 if (sc) {
1301 if (sc->energy_up)
1302 element_energy += sc->energy_up[i][up];
1303
1304 if (sc->f)
1305 element_energy += sc->f(i, j, k + 1, j, VRNA_DECOMP_ML_STEM, sc->data);
1306 }
1307
1308 if (c[k1j] + element_energy + best_energy <= threshold)
1309 repeat(vc, k + 1, j, state, element_energy, 0, best_energy, threshold, env);
1310 }
1311 }
1312 }
1313 } /* array_flag == 1 */
1314
1315 /* 22222222222222222222222222222222222222222222222222 */
1316 /* */
1317 /* array_flag = 2: interval i,j was generated from a */
1318 /* stack, bulge, or internal loop in repeat() */
1319 /* */
1320 /* 22222222222222222222222222222222222222222222222222 */
1321 if (array_flag == 2) {
1322 repeat(vc, i, j, state, 0, 0, best_energy, threshold, env);
1323
1324 if (env->nopush)
1325 if (!noLP)
1326 vrna_message_warning("%d,%d\nOops, no solution in repeat!", i, j);
1327
1328 return;
1329 }
1330
1331 /* 00000000000000000000000000000000000000000000000000 */
1332 /* */
1333 /* array_flag = 0: interval i,j was found while */
1334 /* tracing back through f5-array and c-array */
1335 /* or within this block */
1336 /* */
1337 /* 00000000000000000000000000000000000000000000000000 */
1338 if ((array_flag == 0) && !circular) {
1339 int s5, s3, kj, tmp_en;
1340
1341 if ((hc->up_ext[j]) &&
1342 (f5[j - 1] != INF)) {
1343 tmp_en = 0;
1344
1345 if (sc) {
1346 if (sc->energy_up)
1347 tmp_en += sc->energy_up[j][1];
1348
1349 if (sc->f)
1350 tmp_en += sc->f(1, j, 1, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
1351 }
1352
1353 if (f5[j - 1] + tmp_en + best_energy <= threshold)
1354 /* no basepair, nibbling of 3'-end */
1355 fork_state(i, j - 1, state, tmp_en, 0, env);
1356 }
1357
1358 for (k = j - turn - 1; k > 1; k--) {
1359 kj = indx[j] + k;
1360
1361 if ((with_gquad) &&
1362 (sn[k] == sn[j]) &&
1363 (f5[k - 1] != INF) &&
1364 (ggg[kj] != INF)) {
1365 element_energy = 0;
1366
1367 if (f5[k - 1] + ggg[kj] + element_energy + best_energy <= threshold) {
1368 temp_state = derive_new_state(1, k - 1, state, 0, 0);
1369 env->nopush = false;
1370 /* backtrace the quadruplex */
1371 repeat_gquad(vc,
1372 k,
1373 j,
1374 temp_state,
1375 element_energy,
1376 f5[k - 1],
1377 best_energy,
1378 threshold,
1379 env);
1380 free_state_node(temp_state);
1381 }
1382 }
1383
1384 if ((hard_constraints[length * j + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1385 (f5[k - 1] != INF) &&
1386 (c[kj] != INF)) {
1387 type = vrna_get_ptype(kj, ptype);
1388
1389 /* k and j pair */
1390 switch (dangle_model) {
1391 case 0:
1392 s5 = s3 = -1;
1393 break;
1394 default:
1395 s5 = (sn[k - 1] == sn[k]) ? S1[k - 1] : -1;
1396 s3 = ((j < length) && (sn[j] == sn[j + 1])) ? S1[j + 1] : -1;
1397 break;
1398 }
1399
1400 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1401
1402 if (sn[k] != sn[j]) /*&&(state->is_duplex==0))*/
1403 element_energy += P->DuplexInit;
1404
1405 /*state->is_duplex=1;*/
1406 if (sc) {
1407 if (sc->f)
1408 element_energy += sc->f(1, j, k - 1, k, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
1409 }
1410
1411 if (f5[k - 1] + c[kj] + element_energy + best_energy <= threshold) {
1412 temp_state = derive_new_state(1, k - 1, state, 0, 0);
1413 env->nopush = false;
1414 repeat(vc, k, j, temp_state, element_energy, f5[k - 1], best_energy, threshold, env);
1415 free_state_node(temp_state);
1416 }
1417 }
1418 }
1419
1420 kj = indx[j] + 1;
1421
1422 if ((with_gquad) &&
1423 (sn[k] == sn[j]) &&
1424 (ggg[kj] != INF)) {
1425 element_energy = 0;
1426
1427 if (ggg[kj] + element_energy + best_energy <= threshold)
1428 /* backtrace the quadruplex */
1429 repeat_gquad(vc, 1, j, state, element_energy, 0, best_energy, threshold, env);
1430 }
1431
1432 if ((hard_constraints[length + j] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1433 (c[kj] != INF)) {
1434 type = vrna_get_ptype(kj, ptype);
1435 s5 = -1;
1436
1437 switch (dangle_model) {
1438 case 0:
1439 s3 = -1;
1440 break;
1441 default:
1442 s3 = (j < length) && (sn[j] == sn[j + 1]) ? S1[j + 1] : -1;
1443 break;
1444 }
1445
1446 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1447
1448 if (sn[1] != sn[j])
1449 element_energy += P->DuplexInit;
1450
1451 if (sc) {
1452 if (sc->f)
1453 element_energy += sc->f(1, j, 1, j, VRNA_DECOMP_EXT_STEM, sc->data);
1454 }
1455
1456 if (c[kj] + element_energy + best_energy <= threshold)
1457 repeat(vc, 1, j, state, element_energy, 0, best_energy, threshold, env);
1458 }
1459 } /* end array_flag == 0 && !circular*/
1460 /* or do we subopt circular? */
1461 else if (array_flag == 0) {
1462 int k, l, p, q, tmp_en;
1463 /* if we've done everything right, we will never reach this case more than once */
1464 /* right after the initilization of the stack with ([1,n], empty, 0) */
1465 /* lets check, if we can have an open chain without breaking the threshold */
1466 /* this is an ugly work-arround cause in case of an open chain we do not have to */
1467 /* backtrack anything further... */
1468 if (hc->up_ext[1] >= length) {
1469 tmp_en = 0;
1470
1471 if (sc) {
1472 if (sc->energy_up)
1473 tmp_en += sc->energy_up[1][length];
1474
1475 if (sc->f)
1476 tmp_en += sc->f(1, j, 1, j, VRNA_DECOMP_EXT_UP, sc->data);
1477 }
1478
1479 if (tmp_en <= threshold) {
1480 new_state = derive_new_state(1, 2, state, 0, 0);
1481 new_state->partial_energy = 0;
1482 push(env->Stack, new_state);
1483 env->nopush = false;
1484 }
1485 }
1486
1487 /* ok, lets check if we can do an exterior hairpin without breaking the threshold */
1488 /* best energy should be 0 if we are here */
1489 if (FcH + best_energy <= threshold) {
1490 /* lets search for all exterior hairpin cases, that fit into our threshold barrier */
1491 /* we use index k,l to avoid confusion with i,j index of our state... */
1492 /* if we reach here, i should be 1 and j should be n respectively */
1493 for (k = i; k < j; k++) {
1494 if (hc->up_hp[1] < k)
1495 break;
1496
1497 for (l = j; l >= k + turn + 1; l--) {
1498 int kl, tmpE;
1499
1500 kl = indx[l] + k;
1501
1502 if (c[kl] != INF) {
1503 tmpE = vrna_E_hp_loop(vc, l, k);
1504
1505 if (c[kl] + tmpE + best_energy <= threshold) {
1506 /* what we really have to do is something like this, isn't it? */
1507 /* we have to create a new state, with interval [k,l], then we */
1508 /* add our loop energy as initial energy of this state and put */
1509 /* the state onto the stack R... for further refinement... */
1510 /* we also denote this new interval to be scanned in C */
1511 fork_state(k, l, state, tmpE, 2, env);
1512 }
1513 }
1514 }
1515 }
1516 }
1517
1518 /* now lets see, if we can do an exterior interior loop without breaking the threshold */
1519 if (FcI + best_energy <= threshold) {
1520 /* now we search for our exterior interior loop possibilities */
1521 for (k = i; k < j; k++) {
1522 for (l = j; l >= k + turn + 1; l--) {
1523 int kl, type, tmpE;
1524
1525 kl = indx[l] + k; /* just confusing these indices ;-) */
1526
1527 if ((hard_constraints[length * k + l] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) &&
1528 (c[kl] != INF)) {
1529 type = rtype[vrna_get_ptype(kl, ptype)];
1530
1531 for (p = l + 1; p < j; p++) {
1532 int u1, qmin;
1533 u1 = p - l - 1;
1534 if (u1 + k - 1 > MAXLOOP)
1535 break;
1536
1537 if (hc->up_int[l + 1] < u1)
1538 break;
1539
1540 qmin = u1 + k - 1 + j - MAXLOOP;
1541 if (qmin < p + turn + 1)
1542 qmin = p + turn + 1;
1543
1544 for (q = j; q >= qmin; q--) {
1545 int u2, type_2;
1546
1547 if (hc->up_int[q + 1] < (j - q + k - 1))
1548 break;
1549
1550 if ((hard_constraints[length * p + q] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) &&
1551 (c[indx[q] + p] != INF)) {
1552 type_2 = rtype[vrna_get_ptype(indx[q] + p, ptype)];
1553
1554 u2 = k - 1 + j - q;
1555 if (u1 + u2 > MAXLOOP)
1556 continue;
1557
1558 tmpE = E_IntLoop(u1,
1559 u2,
1560 type,
1561 type_2,
1562 S1[l + 1],
1563 S1[k - 1],
1564 S1[p - 1],
1565 S1[q + 1],
1566 P);
1567
1568 if (sc) {
1569 if (sc->energy_up)
1570 tmpE += sc->energy_up[l + 1][p - l - 1]
1571 + sc->energy_up[q + 1][j - q]
1572 + sc->energy_up[1][k - 1];
1573
1574 if (sc->energy_stack) {
1575 if (u1 + u2 == 0) {
1576 tmpE += sc->energy_stack[k]
1577 + sc->energy_stack[l]
1578 + sc->energy_stack[p]
1579 + sc->energy_stack[q];
1580 }
1581 }
1582 }
1583
1584 if (c[kl] + c[indx[q] + p] + tmpE + best_energy <= threshold) {
1585 /* ok, similar to the hairpin stuff, we add new states onto the stack R */
1586 /* but in contrast to the hairpin decomposition, we have to add two new */
1587 /* intervals, enclosed by k,l and p,q respectively and we also have to */
1588 /* add the partial energy, that comes from the exterior interior loop */
1589 fork_two_states(k, l, p, q, state, tmpE, 2, 2, env);
1590 }
1591 }
1592 }
1593 }
1594 }
1595 }
1596 }
1597 }
1598
1599 /* and last but not least, we have a look, if we can do an exterior multiloop within the energy threshold */
1600 if (FcM <= threshold) {
1601 /* this decomposition will be somehow more complicated...so lets see what we do here... */
1602 /* first we want to find out which split inidices we can use without exceeding the threshold */
1603 int tmpE2;
1604 for (k = turn + 1; k < j - 2 * turn; k++) {
1605 if ((fML[indx[k] + 1] != INF) &&
1606 (fM2[k + 1] != INF)) {
1607 tmpE2 = fML[indx[k] + 1] + fM2[k + 1] + P->MLclosing;
1608
1609 if (tmpE2 + best_energy <= threshold) {
1610 /* grmpfh, we have found a possible split index k so we have to split fM2 and fML now */
1611 /* lets do it first in fM2 anyway */
1612 for (l = k + turn + 2; l < j - turn - 1; l++) {
1613 tmpE2 = fM1[indx[l] + k + 1] + fM1[indx[j] + l + 1];
1614 if (tmpE2 + fML[indx[k] + 1] + P->MLclosing <= threshold) {
1615 /* we've (hopefully) found a valid decomposition of fM2 and therefor we have all */
1616 /* three intervals for our new state to be pushed on stack R */
1617 new_state = copy_state(state);
1618
1619 /* first interval leads for search in fML array */
1620 new_interval = make_interval(1, k, 1);
1621 push(new_state->Intervals, new_interval);
1622 env->nopush = false;
1623
1624 /* next, we have the first interval that has to be traced in fM1 */
1625 new_interval = make_interval(k + 1, l, 3);
1626 push(new_state->Intervals, new_interval);
1627 env->nopush = false;
1628
1629 /* and the last of our three intervals is also one to be traced within fM1 array... */
1630 new_interval = make_interval(l + 1, j, 3);
1631 push(new_state->Intervals, new_interval);
1632 env->nopush = false;
1633
1634 /* mmh, we add the energy for closing the multiloop now... */
1635 new_state->partial_energy += P->MLclosing;
1636 /* next we push our state onto the R stack */
1637 push(env->Stack, new_state);
1638 env->nopush = false;
1639 }
1640
1641 /* else we search further... */
1642 }
1643
1644 /* ok, we have to decompose fML now... */
1645 }
1646 }
1647 }
1648 }
1649 } /* thats all folks for the circular case... */
1650
1651 /* 44444444444444444444444444444444444444444444444444 */
1652 /* */
1653 /* array_flag = 4: interval i,j was found while */
1654 /* tracing back through fc-array smaller than than cp */
1655 /* or within this block */
1656 /* */
1657 /* 44444444444444444444444444444444444444444444444444 */
1658 if (array_flag == 4) {
1659 int ik, s5, s3, tmp_en;
1660
1661 if ((hc->up_ext[i]) &&
1662 (fc[i + 1] != INF)) {
1663 tmp_en = 0;
1664
1665 if (sc) {
1666 if (sc->energy_up)
1667 tmp_en += sc->energy_up[i][1];
1668
1669 if (sc->f)
1670 tmp_en += sc->f(i, j, i + 1, j, VRNA_DECOMP_EXT_EXT, sc->data);
1671 }
1672
1673 if (fc[i + 1] + tmp_en + best_energy <= threshold)
1674 /* no basepair, nibbling of 5'-end */
1675 fork_state(i + 1, j, state, tmp_en, 4, env);
1676 }
1677
1678 for (k = i + turn + 1; k < j; k++) {
1679 ik = indx[k] + i;
1680
1681 if ((with_gquad) &&
1682 (fc[k + 1] != INF) &&
1683 (ggg[ik] != INF)) {
1684 if (fc[k + 1] + ggg[ik] + best_energy <= threshold) {
1685 temp_state = derive_new_state(k + 1, j, state, 0, 4);
1686 env->nopush = false;
1687 repeat_gquad(vc, i, k, temp_state, 0, fc[k + 1], best_energy, threshold, env);
1688 free_state_node(temp_state);
1689 }
1690 }
1691
1692 if ((hard_constraints[length * i + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1693 (fc[k + 1] != INF) &&
1694 (c[ik] != INF)) {
1695 type = vrna_get_ptype(ik, ptype);
1696
1697 switch (dangle_model) {
1698 case 0:
1699 s5 = s3 = -1;
1700 break;
1701 default:
1702 s5 = (i > 1) ? S1[i - 1] : -1;
1703 s3 = S1[k + 1];
1704 break;
1705 }
1706
1707 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1708
1709 if (sc) {
1710 if (sc->f)
1711 element_energy += sc->f(i, j, k, k + 1, VRNA_DECOMP_EXT_STEM_EXT, sc->data);
1712 }
1713
1714 if (fc[k + 1] + c[ik] + element_energy + best_energy <= threshold) {
1715 temp_state = derive_new_state(k + 1, j, state, 0, 4);
1716 env->nopush = false;
1717 repeat(vc, i, k, temp_state, element_energy, fc[k + 1], best_energy, threshold, env);
1718 free_state_node(temp_state);
1719 }
1720 }
1721 }
1722
1723 ik = indx[se[so[0]]] + i; /* indx[j] + i; */
1724
1725 if ((with_gquad) && (ggg[ik] != INF))
1726 if (ggg[ik] + best_energy <= threshold)
1727 repeat_gquad(vc, i, se[so[0]], state, 0, 0, best_energy, threshold, env);
1728
1729 if ((hard_constraints[length * i + se[so[0]]] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1730 (c[ik] != INF)) {
1731 type = vrna_get_ptype(ik, ptype);
1732 s3 = -1;
1733
1734 switch (dangle_model) {
1735 case 0:
1736 s5 = -1;
1737 break;
1738 default:
1739 s5 = (i > 1) ? S1[i - 1] : -1;
1740 break;
1741 }
1742
1743 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1744
1745 if (sc) {
1746 if (sc->f)
1747 element_energy += sc->f(i, se[so[0]], i, se[so[0]], VRNA_DECOMP_EXT_STEM, sc->data);
1748 }
1749
1750 if (c[ik] + element_energy + best_energy <= threshold)
1751 repeat(vc, i, se[so[0]], state, element_energy, 0, best_energy, threshold, env);
1752 }
1753 } /* array_flag == 4 */
1754
1755 /* 55555555555555555555555555555555555555555555555555 */
1756 /* */
1757 /* array_flag = 5: interval cp=i,j was found while */
1758 /* tracing back through fc-array greater than cp */
1759 /* or within this block */
1760 /* */
1761 /* 55555555555555555555555555555555555555555555555555 */
1762 if (array_flag == 5) {
1763 int kj, s5, s3, tmp_en;
1764
1765 if ((hc->up_ext[j]) &&
1766 (fc[j - 1] != INF)) {
1767 tmp_en = 0;
1768
1769 if (sc) {
1770 if (sc->energy_up)
1771 tmp_en += sc->energy_up[j][1];
1772
1773 if (sc->f)
1774 tmp_en += sc->f(i, j, i, j - 1, VRNA_DECOMP_EXT_EXT, sc->data);
1775 }
1776
1777 if (fc[j - 1] + tmp_en + best_energy <= threshold)
1778 /* no basepair, nibbling of 3'-end */
1779 fork_state(i, j - 1, state, tmp_en, 5, env);
1780 }
1781
1782 for (k = j - turn - 1; k > i; k--) {
1783 kj = indx[j] + k;
1784
1785 if ((with_gquad) &&
1786 (fc[k - 1] != INF) &&
1787 (ggg[kj] != INF)) {
1788 if (fc[k - 1] + ggg[kj] + best_energy <= threshold) {
1789 temp_state = derive_new_state(i, k - 1, state, 0, 5);
1790 env->nopush = false;
1791 repeat_gquad(vc, k, j, temp_state, 0, fc[k - 1], best_energy, threshold, env);
1792 free_state_node(temp_state);
1793 }
1794 }
1795
1796 if ((hard_constraints[length * j + k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1797 (fc[k - 1] != INF) &&
1798 (c[kj] != INF)) {
1799 type = vrna_get_ptype(kj, ptype);
1800 element_energy = 0;
1801
1802 switch (dangle_model) {
1803 case 0:
1804 s3 = s5 = -1;
1805 break;
1806 default:
1807 s5 = S1[k - 1];
1808 s3 = (j < length) ? S1[j + 1] : -1;
1809 break;
1810 }
1811
1812 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1813
1814 if (sc) {
1815 if (sc->f)
1816 element_energy += sc->f(i, j, k - 1, k, VRNA_DECOMP_EXT_EXT_STEM, sc->data);
1817 }
1818
1819 if (fc[k - 1] + c[kj] + element_energy + best_energy <= threshold) {
1820 temp_state = derive_new_state(i, k - 1, state, 0, 5);
1821 env->nopush = false;
1822 repeat(vc, k, j, temp_state, element_energy, fc[k - 1], best_energy, threshold, env);
1823 free_state_node(temp_state);
1824 }
1825 }
1826 }
1827
1828 kj = indx[j] + ss[so[1]]; /* indx[j] + i; */
1829
1830 if ((with_gquad) && (ggg[kj] != INF))
1831 if (ggg[kj] + best_energy <= threshold)
1832 repeat_gquad(vc, ss[so[1]], j, state, 0, 0, best_energy, threshold, env);
1833
1834 if ((hard_constraints[length * ss[so[1]] + j] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
1835 (c[kj] != INF)) {
1836 type = vrna_get_ptype(kj, ptype);
1837 s5 = -1;
1838
1839 switch (dangle_model) {
1840 case 0:
1841 s3 = -1;
1842 break;
1843 default:
1844 s3 = (j < length) ? S1[j + 1] : -1;
1845 break;
1846 }
1847
1848 element_energy = vrna_E_ext_stem(type, s5, s3, P);
1849
1850 if (sc) {
1851 if (sc->f)
1852 element_energy += sc->f(ss[so[1]], j, ss[so[1]], j, VRNA_DECOMP_EXT_STEM, sc->data);
1853 }
1854
1855 if (c[kj] + element_energy + best_energy <= threshold)
1856 repeat(vc, ss[so[1]], j, state, element_energy, 0, best_energy, threshold, env);
1857 }
1858 } /* array_flag == 5 */
1859
1860 if (array_flag == 6) {
1861 /* we have a gquad */
1862 repeat_gquad(vc, i, j, state, 0, 0, best_energy, threshold, env);
1863 if (env->nopush)
1864 vrna_message_warning("%d,%d\nOops, no solution in gquad-repeat!", i, j);
1865
1866 return;
1867 }
1868
1869 if (env->nopush) {
1870 push_back(env->Stack, state);
1871 env->nopush = false;
1872 }
1873
1874 return;
1875 }
1876
1877
1878 /*---------------------------------------------------------------------------*/
1879 PRIVATE void
repeat_gquad(vrna_fold_compound_t * vc,int i,int j,STATE * state,int part_energy,int temp_energy,int best_energy,int threshold,subopt_env * env)1880 repeat_gquad(vrna_fold_compound_t *vc,
1881 int i,
1882 int j,
1883 STATE *state,
1884 int part_energy,
1885 int temp_energy,
1886 int best_energy,
1887 int threshold,
1888 subopt_env *env)
1889 {
1890 unsigned int *sn;
1891 int *ggg, *indx, element_energy;
1892 short *S1;
1893 vrna_param_t *P;
1894
1895 indx = vc->jindx;
1896 sn = vc->strand_number;
1897 ggg = vc->matrices->ggg;
1898 S1 = vc->sequence_encoding;
1899 P = vc->params;
1900
1901
1902 /* find all gquads that fit into the energy range and the interval [i,j] */
1903 STATE *new_state;
1904 best_energy += part_energy; /* energy of current structural element */
1905 best_energy += temp_energy; /* energy from unpushed interval */
1906
1907 if (sn[i] == sn[j]) {
1908 element_energy = ggg[indx[j] + i];
1909 if ((element_energy != INF) &&
1910 (element_energy + best_energy <= threshold)) {
1911 int cnt;
1912 int *L;
1913 int *l;
1914 /* find out how many gquads we might expect in the interval [i,j] */
1915 int num_gquads = get_gquad_count(S1, i, j);
1916 num_gquads++;
1917 L = (int *)vrna_alloc(sizeof(int) * num_gquads);
1918 l = (int *)vrna_alloc(sizeof(int) * num_gquads * 3);
1919 L[0] = -1;
1920
1921 get_gquad_pattern_exhaustive(S1, i, j, P, L, l, threshold - best_energy);
1922
1923 for (cnt = 0; L[cnt] != -1; cnt++) {
1924 new_state = copy_state(state);
1925
1926 make_gquad(i, L[cnt], &(l[3 * cnt]), new_state);
1927 new_state->partial_energy += part_energy;
1928 new_state->partial_energy += element_energy;
1929 /* new_state->best_energy =
1930 * hairpin[unpaired] + element_energy + best_energy; */
1931 push(env->Stack, new_state);
1932 env->nopush = false;
1933 }
1934 free(L);
1935 free(l);
1936 }
1937 }
1938
1939 best_energy -= part_energy;
1940 best_energy -= temp_energy;
1941 return;
1942 }
1943
1944
1945 PRIVATE void
repeat(vrna_fold_compound_t * vc,int i,int j,STATE * state,int part_energy,int temp_energy,int best_energy,int threshold,subopt_env * env)1946 repeat(vrna_fold_compound_t *vc,
1947 int i,
1948 int j,
1949 STATE *state,
1950 int part_energy,
1951 int temp_energy,
1952 int best_energy,
1953 int threshold,
1954 subopt_env *env)
1955 {
1956 /* routine to find stacks, bulges, internal loops and multiloops */
1957 /* within interval closed by basepair i,j */
1958
1959 STATE *new_state;
1960 vrna_param_t *P;
1961 vrna_md_t *md;
1962
1963 register int ij, k, p, q, energy, new;
1964 register int mm;
1965 register int no_close, type, type_2;
1966 char *ptype;
1967 unsigned int n, *sn, *so, *ss, *se;
1968 int element_energy;
1969 int *fc, *c, *fML, *fM1, *ggg;
1970 int rt, *indx, *rtype, noGUclosure, noLP, with_gquad, dangle_model, turn;
1971 short *S1;
1972 vrna_hc_t *hc;
1973 vrna_sc_t *sc;
1974
1975 n = vc->length;
1976 S1 = vc->sequence_encoding;
1977 ptype = vc->ptype;
1978 indx = vc->jindx;
1979 sn = vc->strand_number;
1980 so = vc->strand_order;
1981 ss = vc->strand_start;
1982 se = vc->strand_end;
1983 P = vc->params;
1984 md = &(P->model_details);
1985 rtype = &(md->rtype[0]);
1986
1987 noGUclosure = md->noGUclosure;
1988 noLP = md->noLP;
1989 with_gquad = md->gquad;
1990 dangle_model = md->dangles;
1991 turn = md->min_loop_size;
1992
1993 fc = vc->matrices->fc;
1994 c = vc->matrices->c;
1995 fML = vc->matrices->fML;
1996 fM1 = vc->matrices->fM1;
1997 ggg = vc->matrices->ggg;
1998
1999 hc = vc->hc;
2000 sc = vc->sc;
2001
2002 ij = indx[j] + i;
2003
2004 type = vrna_get_ptype(ij, ptype);
2005 /*
2006 * if (type==0) fprintf(stderr, "repeat: Warning: %d %d can't pair\n", i,j);
2007 */
2008
2009 no_close = (((type == 3) || (type == 4)) && noGUclosure);
2010
2011 if (hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
2012 if (noLP) {
2013 /* always consider the structure with additional stack */
2014 if (i + turn + 2 < j) {
2015 if (hc->mx[n * (i + 1) + j - 1] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC) {
2016 type_2 = rtype[vrna_get_ptype(indx[j - 1] + i + 1, ptype)];
2017 energy = 0;
2018
2019 if ((sn[i] == sn[i + 1]) && (sn[j - 1] == sn[j])) {
2020 energy = E_IntLoop(0, 0, type, type_2, S1[i + 1], S1[j - 1], S1[i + 1], S1[j - 1], P);
2021
2022 if (sc) {
2023 if (sc->energy_bp)
2024 energy += sc->energy_bp[ij];
2025
2026 if (sc->energy_stack) {
2027 energy += sc->energy_stack[i]
2028 + sc->energy_stack[i + 1]
2029 + sc->energy_stack[j - 1]
2030 + sc->energy_stack[j];
2031 }
2032
2033 if (sc->f)
2034 energy += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_IL, sc->data);
2035 }
2036
2037 new_state = derive_new_state(i + 1, j - 1, state, part_energy + energy, 2);
2038 make_pair(i, j, new_state);
2039 make_pair(i + 1, j - 1, new_state);
2040
2041 /* new_state->best_energy = new + best_energy; */
2042 push(env->Stack, new_state);
2043 env->nopush = false;
2044 if (i == 1 || state->structure[i - 2] != '(' || state->structure[j] != ')')
2045 /* adding a stack is the only possible structure */
2046 return;
2047 }
2048 }
2049 }
2050 }
2051 }
2052
2053 best_energy += part_energy; /* energy of current structural element */
2054 best_energy += temp_energy; /* energy from unpushed interval */
2055
2056 if (hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
2057 for (p = i + 1; p <= MIN2(j - 2 - turn, i + MAXLOOP + 1); p++) {
2058 int minq = j - i + p - MAXLOOP - 2;
2059 if (minq < p + 1 + turn)
2060 minq = p + 1 + turn;
2061
2062 if (hc->up_int[i + 1] < (p - i - 1))
2063 break;
2064
2065 for (q = j - 1; q >= minq; q--) {
2066 if (hc->up_int[q + 1] < (j - q - 1))
2067 break;
2068
2069 /* skip stack if noLP, since we've already processed it above */
2070 if ((noLP) && (p == i + 1) && (q == j - 1))
2071 continue;
2072
2073 if (!(hc->mx[n * p + q] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC))
2074 continue;
2075
2076 if (c[indx[q] + p] == INF)
2077 continue;
2078
2079 type_2 = vrna_get_ptype(indx[q] + p, ptype);
2080
2081 if (noGUclosure)
2082 if (no_close || (type_2 == 3) || (type_2 == 4))
2083 if ((p > i + 1) || (q < j - 1))
2084 continue;
2085
2086 /* continue unless stack */
2087
2088 if ((sn[i] == sn[p]) && (sn[q] == sn[j])) {
2089 energy = E_IntLoop(p - i - 1, j - q - 1, type, rtype[type_2],
2090 S1[i + 1], S1[j - 1], S1[p - 1], S1[q + 1], P);
2091
2092 new = energy + c[indx[q] + p];
2093
2094 if (sc) {
2095 if (sc->energy_up)
2096 energy += sc->energy_up[i + 1][p - i - 1]
2097 + sc->energy_up[q + 1][j - q - 1];
2098
2099 if (sc->energy_bp)
2100 energy += sc->energy_bp[ij];
2101
2102 if (sc->energy_stack) {
2103 if ((p == i + 1) && (q == j - 1)) {
2104 energy += sc->energy_stack[i]
2105 + sc->energy_stack[p]
2106 + sc->energy_stack[q]
2107 + sc->energy_stack[j];
2108 }
2109 }
2110
2111 if (sc->f)
2112 energy += sc->f(i, j, p, q, VRNA_DECOMP_PAIR_IL, sc->data);
2113 }
2114
2115 new = energy + c[indx[q] + p];
2116
2117 if (new + best_energy <= threshold)
2118 /* stack, bulge, or interior loop */
2119 fork_int_state(i, j, p, q, state, part_energy + energy, env);
2120 } /*end of if block */
2121 } /* end of q-loop */
2122 } /* end of p-loop */
2123 }
2124
2125 if (sn[i] != sn[j]) {
2126 /*look in fc*/
2127 if ((hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) &&
2128 (fc[i + 1] != INF) &&
2129 (fc[j - 1] != INF)) {
2130 rt = rtype[type];
2131
2132 element_energy = 0;
2133 switch (dangle_model) {
2134 case 0:
2135 element_energy = vrna_E_ext_stem(rt, -1, -1, P);
2136 break;
2137 default:
2138 element_energy =
2139 vrna_E_ext_stem(rt,
2140 (sn[j - 1] == sn[j]) ?
2141 S1[j - 1] :
2142 -1,
2143 (sn[i] == sn[i + 1]) ?
2144 S1[i + 1] :
2145 -1,
2146 P);
2147 break;
2148 }
2149
2150 if (fc[i + 1] + fc[j - 1] + element_energy + best_energy <= threshold)
2151 fork_two_states_pair(i, j, ss[so[1]], state, part_energy + element_energy, 4, 5, env);
2152 }
2153 }
2154
2155 mm = P->MLclosing;
2156 rt = rtype[type];
2157
2158 if ((hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) &&
2159 ((vc->strands < 2) || ((i != se[so[0]]) && (j != ss[so[1]])))) {
2160 element_energy = mm;
2161 switch (dangle_model) {
2162 case 0:
2163 element_energy = E_MLstem(rt, -1, -1, P) + mm;
2164 break;
2165 default:
2166 element_energy = E_MLstem(rt, S1[j - 1], S1[i + 1], P) + mm;
2167 break;
2168 }
2169
2170 if (sc) {
2171 if (sc->energy_bp)
2172 element_energy += sc->energy_bp[ij];
2173
2174 if (sc->f)
2175 element_energy += sc->f(i, j, i + 1, j - 1, VRNA_DECOMP_PAIR_ML, sc->data);
2176 }
2177
2178 /* multiloop decomposition */
2179 if ((sc) && (sc->f)) {
2180 for (k = i + turn + 2; k <= j - turn - 2; k++) {
2181 int eee = fML[indx[k - 1] + i + 1];
2182
2183 if ((eee != INF) && (fM1[indx[j - 1] + k] != INF)) {
2184 eee += fM1[indx[j - 1] + k] +
2185 best_energy;
2186 int aux_eee = element_energy +
2187 sc->f(i + 1, j - 1, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data);
2188 if ((eee + aux_eee) <= threshold)
2189 fork_two_states_pair(i, j, k, state, part_energy + aux_eee, 1, 3, env);
2190 }
2191 }
2192 } else {
2193 for (k = i + turn + 2; k <= j - turn - 2; k++) {
2194 int eee = fML[indx[k - 1] + i + 1];
2195
2196 if ((eee != INF) && (fM1[indx[j - 1] + k] != INF)) {
2197 /* multiloop decomposition */
2198 if ((eee + fM1[indx[j - 1] + k] +
2199 element_energy + best_energy) <= threshold)
2200 fork_two_states_pair(i, j, k, state, part_energy + element_energy, 1, 3, env);
2201 }
2202 }
2203 }
2204 }
2205
2206 if (sn[i] == sn[j]) {
2207 if ((hc->mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) &&
2208 (!no_close)) {
2209
2210 element_energy = vrna_E_hp_loop(vc, i, j);
2211
2212 if (element_energy != INF) {
2213 if (element_energy + best_energy <= threshold)
2214 /* hairpin structure */
2215 fork_state_pair(i, j, state, part_energy + element_energy, env);
2216 }
2217 }
2218
2219 if (with_gquad) {
2220 /* now we have to find all loops where (i,j) encloses a gquad in an interior loops style */
2221 int cnt, *p, *q, *en, tmp_en;
2222 p = q = en = NULL;
2223 en =
2224 E_GQuad_IntLoop_exhaustive(i, j, &p, &q, type, S1, ggg, threshold - best_energy, indx, P);
2225 for (cnt = 0; p[cnt] != -1; cnt++) {
2226 if ((hc->up_int[i + 1] >= p[cnt] - i - 1) && (hc->up_int[q[cnt] + 1] >= j - q[cnt] - 1)) {
2227 tmp_en = en[cnt];
2228
2229 if (sc) {
2230 if (sc->energy_bp)
2231 tmp_en += sc->energy_bp[ij];
2232
2233 if (sc->energy_up)
2234 tmp_en += sc->energy_up[i + 1][p[cnt] - i - 1]
2235 + sc->energy_up[q[cnt] + 1][j - q[cnt] - 1];
2236 }
2237
2238 new_state = derive_new_state(p[cnt], q[cnt], state, tmp_en + part_energy, 6);
2239
2240 make_pair(i, j, new_state);
2241
2242 /* new_state->best_energy = new + best_energy; */
2243 push(env->Stack, new_state);
2244 env->nopush = false;
2245 }
2246 }
2247 free(en);
2248 free(p);
2249 free(q);
2250 }
2251 }
2252
2253 best_energy -= part_energy;
2254 best_energy -= temp_energy;
2255 return;
2256 }
2257
2258
2259 PRIVATE void
old_subopt_print(const char * structure,float energy,void * data)2260 old_subopt_print(const char *structure,
2261 float energy,
2262 void *data)
2263 {
2264 struct old_subopt_dat *d = (struct old_subopt_dat *)data;
2265
2266 if (structure && d->fp) {
2267 char *e_string = vrna_strdup_printf(" %6.2f", energy);
2268 print_structure(d->fp, structure, e_string);
2269 free(e_string);
2270 }
2271 }
2272
2273
2274 PRIVATE void
old_subopt_store(const char * structure,float energy,void * data)2275 old_subopt_store(const char *structure,
2276 float energy,
2277 void *data)
2278 {
2279 struct old_subopt_dat *d = (struct old_subopt_dat *)data;
2280
2281 /* store solution */
2282 if (d->n_sol + 1 == d->max_sol) {
2283 d->max_sol *= 2;
2284 d->SolutionList = (SOLUTION *)vrna_realloc(d->SolutionList, d->max_sol * sizeof(SOLUTION));
2285 }
2286
2287 if (structure) {
2288 d->SolutionList[d->n_sol].energy = energy;
2289 d->SolutionList[d->n_sol++].structure = strdup(structure);
2290 } else {
2291 d->SolutionList[d->n_sol].energy = 0;
2292 d->SolutionList[d->n_sol++].structure = NULL;
2293 }
2294 }
2295
2296
2297 PRIVATE void
old_subopt_store_compressed(const char * structure,float energy,void * data)2298 old_subopt_store_compressed(const char *structure,
2299 float energy,
2300 void *data)
2301 {
2302 struct old_subopt_dat *d = (struct old_subopt_dat *)data;
2303
2304 /* store solution */
2305 if (d->n_sol + 1 == d->max_sol) {
2306 d->max_sol *= 2;
2307 d->SolutionList = (SOLUTION *)vrna_realloc(d->SolutionList, d->max_sol * sizeof(SOLUTION));
2308 }
2309
2310 if (structure) {
2311 d->SolutionList[d->n_sol].energy = energy;
2312 if (d->cp > 0) {
2313 int cp = d->cp;
2314 char *s = vrna_cut_point_remove(structure, &cp);
2315 d->SolutionList[d->n_sol++].structure = vrna_db_pack(s);
2316 free(s);
2317 } else {
2318 d->SolutionList[d->n_sol++].structure = vrna_db_pack(structure);
2319 }
2320 } else {
2321 d->SolutionList[d->n_sol].energy = 0;
2322 d->SolutionList[d->n_sol++].structure = NULL;
2323 }
2324 }
2325
2326
2327 /*###########################################*/
2328 /*# deprecated functions below #*/
2329 /*###########################################*/
2330 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
2331
2332 PUBLIC SOLUTION *
subopt(char * seq,char * structure,int delta,FILE * fp)2333 subopt(char *seq,
2334 char *structure,
2335 int delta,
2336 FILE *fp)
2337 {
2338 return wrap_subopt(seq, structure, NULL, delta, fold_constrained, 0, fp);
2339 }
2340
2341
2342 PUBLIC SOLUTION *
subopt_circ(char * seq,char * structure,int delta,FILE * fp)2343 subopt_circ(char *seq,
2344 char *structure,
2345 int delta,
2346 FILE *fp)
2347 {
2348 return wrap_subopt(seq, structure, NULL, delta, fold_constrained, 1, fp);
2349 }
2350
2351
2352 PUBLIC SOLUTION *
subopt_par(char * seq,char * structure,vrna_param_t * parameters,int delta,int is_constrained,int is_circular,FILE * fp)2353 subopt_par(char *seq,
2354 char *structure,
2355 vrna_param_t *parameters,
2356 int delta,
2357 int is_constrained,
2358 int is_circular,
2359 FILE *fp)
2360 {
2361 return wrap_subopt(seq, structure, parameters, delta, is_constrained, is_circular, fp);
2362 }
2363
2364
2365 PRIVATE SOLUTION *
wrap_subopt(char * string,char * structure,vrna_param_t * parameters,int delta,int is_constrained,int is_circular,FILE * fp)2366 wrap_subopt(char *string,
2367 char *structure,
2368 vrna_param_t *parameters,
2369 int delta,
2370 int is_constrained,
2371 int is_circular,
2372 FILE *fp)
2373 {
2374 vrna_fold_compound_t *vc;
2375 vrna_param_t *P;
2376 char *seq;
2377
2378 #ifdef _OPENMP
2379 /* Explicitly turn off dynamic threads */
2380 omp_set_dynamic(0);
2381 #endif
2382
2383 /* we need the parameter structure for hard constraints */
2384 if (parameters) {
2385 P = vrna_params_copy(parameters);
2386 } else {
2387 vrna_md_t md;
2388 set_model_details(&md);
2389 md.temperature = temperature;
2390 P = vrna_params(&md);
2391 }
2392
2393 P->model_details.circ = is_circular;
2394 P->model_details.uniq_ML = uniq_ML = 1;
2395
2396 /* what about cofold sequences here? Is it safe to call the below cut_point_insert() ? */
2397 /* dirty hack to reinsert the '&' according to the global variable 'cut_point' */
2398 seq = vrna_cut_point_insert(string, cut_point);
2399
2400 vc =
2401 vrna_fold_compound(seq,
2402 &(P->model_details),
2403 ((is_circular == 0) ? VRNA_OPTION_HYBRID : VRNA_OPTION_DEFAULT));
2404
2405 if (parameters) {
2406 /* replace params if necessary */
2407 free(vc->params);
2408 vc->params = P;
2409 } else {
2410 free(P);
2411 }
2412
2413 /* handle hard constraints in pseudo dot-bracket format if passed via simple interface */
2414 if (is_constrained && structure) {
2415 unsigned int constraint_options = 0;
2416 constraint_options |= VRNA_CONSTRAINT_DB
2417 | VRNA_CONSTRAINT_DB_PIPE
2418 | VRNA_CONSTRAINT_DB_DOT
2419 | VRNA_CONSTRAINT_DB_X
2420 | VRNA_CONSTRAINT_DB_ANG_BRACK
2421 | VRNA_CONSTRAINT_DB_RND_BRACK
2422 | VRNA_CONSTRAINT_DB_INTRAMOL
2423 | VRNA_CONSTRAINT_DB_INTERMOL;
2424
2425 vrna_constraints_add(vc, (const char *)structure, constraint_options);
2426 }
2427
2428 if (backward_compat_compound && backward_compat)
2429 vrna_fold_compound_free(backward_compat_compound);
2430
2431 backward_compat_compound = vc;
2432 backward_compat = 1;
2433
2434 /* cleanup */
2435 free(seq);
2436
2437 return vrna_subopt(vc, delta, subopt_sorted, fp);
2438 }
2439
2440
2441 #endif
2442
2443 /*---------------------------------------------------------------------------*/
2444 /* Well, that is the end!----------------------------------------------------*/
2445 /*---------------------------------------------------------------------------*/
2446