1 /*
2 * partiton function for RNA secondary structures
3 *
4 * Ivo L Hofacker + Ronny Lorenz
5 * Vienna RNA package
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h"
10 #endif
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <string.h>
15 #include <float.h>
16 #include <math.h>
17
18 #include "ViennaRNA/utils/basic.h"
19 #include "ViennaRNA/params/default.h"
20 #include "ViennaRNA/fold_vars.h"
21 #include "ViennaRNA/loops/all.h"
22 #include "ViennaRNA/gquad.h"
23 #include "ViennaRNA/constraints/hard.h"
24 #include "ViennaRNA/constraints/soft.h"
25 #include "ViennaRNA/alphabet.h"
26 #include "ViennaRNA/boltzmann_sampling.h"
27
28 #include "ViennaRNA/loops/external_sc_pf.inc"
29 #include "ViennaRNA/loops/internal_sc_pf.inc"
30 #include "ViennaRNA/loops/multibranch_sc_pf.inc"
31
32 #include "ViennaRNA/data_structures_nonred.inc"
33
34 /*
35 #################################
36 # PREPROCESSOR DEFININTIONS #
37 #################################
38 */
39
40 #ifdef VRNA_NR_SAMPLING_HASH
41 # define NR_NODE tr_node
42 # define NR_TOTAL_WEIGHT(a) total_weight_par(a)
43 # define NR_TOTAL_WEIGHT_TYPE(a, b) total_weight_par_type(a, b)
44 # define NR_GET_WEIGHT(a, b, c, d, e) tr_node_weight(a, c, d, e)
45 #else
46 # define NR_NODE tllr_node
47 # define NR_TOTAL_WEIGHT(a) get_weight_all(a)
48 # define NR_TOTAL_WEIGHT_TYPE(a, b) get_weight_type_spec(a, b)
49 # define NR_GET_WEIGHT(a, b, c, d, e) get_weight(b, c, d, e)
50 #endif
51
52
53 /* combination of soft constraint wrappers */
54 struct sc_wrappers {
55 struct sc_ext_exp_dat sc_wrapper_ext;
56 struct sc_int_exp_dat sc_wrapper_int;
57 struct sc_mb_exp_dat sc_wrapper_ml;
58 };
59
60 /*
61 * In the following:
62 * - q_remain is a pointer to value of sum of Boltzmann factors of still accessible solutions at that point
63 * - current_node is a pointer to current node in datastructure memorizing the solutions and paths taken
64 */
65 struct vrna_pbacktrack_memory_s {
66 double q_remain;
67 NR_NODE *root_node;
68 NR_NODE *current_node;
69 struct nr_memory *memory_dat;
70 };
71
72 /*
73 #################################
74 # GLOBAL VARIABLES #
75 #################################
76 */
77
78 /*
79 #################################
80 # PRIVATE VARIABLES #
81 #################################
82 */
83 PRIVATE char *info_set_uniq_ml =
84 "Unique multiloop decomposition is unset!\n"
85 "Activate unique multiloop decomposition by setting the"
86 " uniq_ML field of the model details structure to a non-zero"
87 " value before running vrna_pf()!";
88
89 PRIVATE char *info_call_pf =
90 "DP matrices are missing! Call vrna_pf() first!";
91
92
93 PRIVATE char *info_nr_duplicates =
94 "Duplicate structures detected, presumably due to numerical instabilities";
95
96
97 PRIVATE char *info_nr_overflow =
98 "Partition function overflow detected for forbidden structures,"
99 " presumably due to numerical instabilities.";
100
101
102 PRIVATE char *info_no_circ =
103 "No implementation for circular RNAs available.";
104
105
106 /*
107 #################################
108 # PRIVATE FUNCTION DECLARATIONS #
109 #################################
110 */
111
112 PRIVATE struct vrna_pbacktrack_memory_s *
113 nr_init(vrna_fold_compound_t *fc);
114
115
116 PRIVATE struct sc_wrappers *
117 sc_init(vrna_fold_compound_t *fc);
118
119
120 PRIVATE void
121 sc_free(struct sc_wrappers *sc_wrap);
122
123
124 PRIVATE unsigned int
125 wrap_pbacktrack(vrna_fold_compound_t *vc,
126 unsigned int length,
127 unsigned int num_samples,
128 vrna_boltzmann_sampling_callback *bs_cb,
129 void *data,
130 struct vrna_pbacktrack_memory_s *nr_mem);
131
132
133 PRIVATE int
134 backtrack(int i,
135 int j,
136 char *pstruc,
137 vrna_fold_compound_t *vc,
138 struct sc_wrappers *sc_wrap,
139 struct vrna_pbacktrack_memory_s *nr_mem);
140
141
142 PRIVATE int
143 backtrack_ext_loop(int init_val,
144 char *pstruc,
145 vrna_fold_compound_t *vc,
146 int length,
147 struct sc_wrappers *sc_wrap,
148 struct vrna_pbacktrack_memory_s *nr_mem);
149
150
151 PRIVATE int
152 backtrack_qm(int i,
153 int j,
154 char *pstruc,
155 vrna_fold_compound_t *vc,
156 struct sc_wrappers *sc_wrap,
157 struct vrna_pbacktrack_memory_s *nr_mem);
158
159
160 PRIVATE int
161 backtrack_qm1(int i,
162 int j,
163 char *pstruc,
164 vrna_fold_compound_t *vc,
165 struct sc_wrappers *sc_wrap,
166 struct vrna_pbacktrack_memory_s *nr_mem);
167
168
169 PRIVATE void
170 backtrack_qm2(int u,
171 int n,
172 char *pstruc,
173 vrna_fold_compound_t *vc,
174 struct sc_wrappers *sc_wrap);
175
176
177 PRIVATE unsigned int
178 pbacktrack_circ(vrna_fold_compound_t *fc,
179 unsigned int num_samples,
180 vrna_boltzmann_sampling_callback *bs_cb,
181 void *data);
182
183
184 /*
185 #################################
186 # BEGIN OF FUNCTION DEFINITIONS #
187 #################################
188 */
189 PUBLIC unsigned int
vrna_pbacktrack5_resume_cb(vrna_fold_compound_t * fc,unsigned int num_samples,unsigned int length,vrna_boltzmann_sampling_callback * bs_cb,void * data,vrna_pbacktrack_mem_t * nr_mem,unsigned int options)190 vrna_pbacktrack5_resume_cb(vrna_fold_compound_t *fc,
191 unsigned int num_samples,
192 unsigned int length,
193 vrna_boltzmann_sampling_callback *bs_cb,
194 void *data,
195 vrna_pbacktrack_mem_t *nr_mem,
196 unsigned int options)
197 {
198 unsigned int i = 0;
199
200 if (fc) {
201 vrna_mx_pf_t *matrices = fc->exp_matrices;
202
203 if (length > fc->length) {
204 vrna_message_warning("vrna_pbacktrack5*(): length exceeds sequence length");
205 } else if (length == 0) {
206 vrna_message_warning("vrna_pbacktrack5*(): length too small");
207 } else if ((!matrices) || (!matrices->q) || (!matrices->qb) || (!matrices->qm) ||
208 (!fc->exp_params)) {
209 vrna_message_warning("vrna_pbacktrack*(): %s", info_call_pf);
210 } else if ((!fc->exp_params->model_details.uniq_ML) || (!matrices->qm1)) {
211 vrna_message_warning("vrna_pbacktrack*(): %s", info_set_uniq_ml);
212 } else if ((fc->exp_params->model_details.circ) && (length < fc->length)) {
213 vrna_message_warning("vrna_pbacktrack5*(): %s", info_no_circ);
214 } else if (options & VRNA_PBACKTRACK_NON_REDUNDANT) {
215 if (fc->exp_params->model_details.circ) {
216 vrna_message_warning("vrna_pbacktrack5*(): %s", info_no_circ);
217 } else if (!nr_mem) {
218 vrna_message_warning("vrna_pbacktrack5*(): Pointer to nr_mem must not be NULL!");
219 } else {
220 if (*nr_mem == NULL)
221 *nr_mem = nr_init(fc);
222
223 i = wrap_pbacktrack(fc, length, num_samples, bs_cb, data, *nr_mem);
224
225 /* print warning if we've aborted backtracking too early */
226 if ((i > 0) && (i < num_samples)) {
227 vrna_message_warning("vrna_pbacktrack5*(): "
228 "Stopped non-redundant backtracking after %d samples"
229 " due to numeric instabilities!\n"
230 "Coverage of partition function so far: %.6f%%",
231 i,
232 100. *
233 return_node_weight((*nr_mem)->root_node) /
234 fc->exp_matrices->q[fc->iindx[1] - length]);
235 }
236 }
237 } else if (fc->exp_params->model_details.circ) {
238 i = pbacktrack_circ(fc, num_samples, bs_cb, data);
239 } else {
240 i = wrap_pbacktrack(fc, length, num_samples, bs_cb, data, NULL);
241 }
242 }
243
244 return i; /* actual number of structures backtraced */
245 }
246
247
248 PUBLIC void
vrna_pbacktrack_mem_free(struct vrna_pbacktrack_memory_s * s)249 vrna_pbacktrack_mem_free(struct vrna_pbacktrack_memory_s *s)
250 {
251 if (s) {
252 #ifdef VRNA_NR_SAMPLING_HASH
253 free_all_nr(s->current_node);
254 #else
255 free_all_nrll(&(s->memory_dat));
256 #endif
257 free(s);
258 }
259 }
260
261
262 /*
263 #####################################
264 # BEGIN OF STATIC HELPER FUNCTIONS #
265 #####################################
266 */
267 PRIVATE struct sc_wrappers *
sc_init(vrna_fold_compound_t * fc)268 sc_init(vrna_fold_compound_t *fc)
269 {
270 struct sc_wrappers *sc_wrap = (struct sc_wrappers *)vrna_alloc(sizeof(struct sc_wrappers));
271
272 init_sc_ext_exp(fc, &(sc_wrap->sc_wrapper_ext));
273 init_sc_int_exp(fc, &(sc_wrap->sc_wrapper_int));
274 init_sc_mb_exp(fc, &(sc_wrap->sc_wrapper_ml));
275
276 return sc_wrap;
277 }
278
279
280 PRIVATE void
sc_free(struct sc_wrappers * sc_wrap)281 sc_free(struct sc_wrappers *sc_wrap)
282 {
283 free_sc_ext_exp(&(sc_wrap->sc_wrapper_ext));
284 free_sc_int_exp(&(sc_wrap->sc_wrapper_int));
285 free_sc_mb_exp(&(sc_wrap->sc_wrapper_ml));
286
287 free(sc_wrap);
288 }
289
290
291 PRIVATE struct vrna_pbacktrack_memory_s *
nr_init(vrna_fold_compound_t * fc)292 nr_init(vrna_fold_compound_t *fc)
293 {
294 size_t block_size;
295 double pf;
296 struct vrna_pbacktrack_memory_s *s;
297
298 s = (struct vrna_pbacktrack_memory_s *)vrna_alloc(
299 sizeof(struct vrna_pbacktrack_memory_s));
300 pf = fc->exp_matrices->q[fc->iindx[1] - fc->length];
301 block_size = 5000 * sizeof(NR_NODE);
302
303 s->memory_dat = NULL;
304 s->q_remain = 0;
305
306 #ifdef VRNA_NR_SAMPLING_HASH
307 s->root_node = create_root(fc->length, pf);
308 #else
309 s->memory_dat = create_nr_memory(sizeof(NR_NODE), block_size, NULL); /* memory pre-allocation */
310 s->root_node = create_ll_root(&(s->memory_dat), pf);
311 #endif
312
313 s->current_node = s->root_node;
314
315 return s;
316 }
317
318
319 /* general expr of vrna5_pbacktrack with possibility of non-redundant sampling */
320 PRIVATE unsigned int
wrap_pbacktrack(vrna_fold_compound_t * vc,unsigned int length,unsigned int num_samples,vrna_boltzmann_sampling_callback * bs_cb,void * data,struct vrna_pbacktrack_memory_s * nr_mem)321 wrap_pbacktrack(vrna_fold_compound_t *vc,
322 unsigned int length,
323 unsigned int num_samples,
324 vrna_boltzmann_sampling_callback *bs_cb,
325 void *data,
326 struct vrna_pbacktrack_memory_s *nr_mem)
327 {
328 char *pstruc;
329 unsigned int i, n;
330 int ret, pf_overflow, is_dup, *my_iindx;
331 FLT_OR_DBL *q1k, *qln, *q;
332 vrna_mx_pf_t *matrices;
333 struct sc_wrappers *sc_wrap;
334
335 i = 0;
336 pf_overflow = 0;
337 sc_wrap = sc_init(vc);
338
339 n = vc->length;
340 my_iindx = vc->iindx;
341 matrices = vc->exp_matrices;
342 q = matrices->q;
343 q1k = matrices->q1k;
344 qln = matrices->qln;
345
346 if (!(q1k && qln)) {
347 matrices->q1k = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 1));
348 matrices->qln = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
349 q1k = matrices->q1k;
350 qln = matrices->qln;
351 for (i = 1; i <= n; i++) {
352 q1k[i] = q[my_iindx[1] - i];
353 qln[i] = q[my_iindx[i] - n];
354 }
355 q1k[0] = 1.0;
356 qln[n + 1] = 1.0;
357 }
358
359 for (i = 0; i < num_samples; i++) {
360 is_dup = 1;
361 pstruc = vrna_alloc((length + 1) * sizeof(char));
362
363 memset(pstruc, '.', sizeof(char) * length);
364
365 if (nr_mem)
366 nr_mem->q_remain = vc->exp_matrices->q[vc->iindx[1] - length];
367
368 #ifdef VRNA_WITH_BOUSTROPHEDON
369 ret = backtrack_ext_loop(length, pstruc, vc, length, sc_wrap, nr_mem);
370 #else
371 ret = backtrack_ext_loop(1, pstruc, vc, length, sc_wrap, nr_mem);
372 #endif
373
374 if (nr_mem) {
375 #ifdef VRNA_NR_SAMPLING_HASH
376 nr_mem->current_node = traceback_to_root(nr_mem->current_node,
377 nr_mem->q_remain,
378 &is_dup,
379 &pf_overflow);
380 #else
381 nr_mem->current_node = traceback_to_ll_root(nr_mem->current_node,
382 nr_mem->q_remain,
383 &is_dup,
384 &pf_overflow);
385 #endif
386
387 if (pf_overflow) {
388 vrna_message_warning("vrna_pbacktrack_nr*(): %s", info_nr_overflow);
389 free(pstruc);
390 break;
391 }
392
393 if (is_dup) {
394 vrna_message_warning("vrna_pbacktrack_nr*(): %s", info_nr_duplicates);
395 free(pstruc);
396 break;
397 }
398 }
399
400 if ((ret > 0) && (bs_cb))
401 bs_cb(pstruc, data);
402
403 free(pstruc);
404
405 if (ret == 0)
406 break;
407 }
408
409 sc_free(sc_wrap);
410
411 return i;
412 }
413
414
415 /* backtrack one external */
416 PRIVATE int
backtrack_ext_loop(int init_val,char * pstruc,vrna_fold_compound_t * vc,int length,struct sc_wrappers * sc_wrap,struct vrna_pbacktrack_memory_s * nr_mem)417 backtrack_ext_loop(int init_val,
418 char *pstruc,
419 vrna_fold_compound_t *vc,
420 int length,
421 struct sc_wrappers *sc_wrap,
422 struct vrna_pbacktrack_memory_s *nr_mem)
423 {
424 unsigned char *hard_constraints;
425 short *S1, *S2, **S, **S5, **S3;
426 unsigned int **a2s, s, n_seq;
427 int ret, i, j, ij, n, k, u, type, *my_iindx, hc_decompose, *hc_up_ext;
428 FLT_OR_DBL r, fbd, fbds, qt, q_temp, qkl, *q, *qb, *q1k, *qln, *scale;
429 double *q_remain;
430 vrna_mx_pf_t *matrices;
431 vrna_md_t *md;
432 vrna_hc_t *hc;
433 vrna_exp_param_t *pf_params;
434
435 struct nr_memory **memory_dat;
436 struct sc_ext_exp_dat *sc_wrapper_ext;
437
438 NR_NODE **current_node;
439
440 if (nr_mem) {
441 q_remain = &(nr_mem->q_remain);
442 current_node = &(nr_mem->current_node);
443 memory_dat = &(nr_mem->memory_dat);
444 } else {
445 q_remain = NULL;
446 current_node = NULL;
447 memory_dat = NULL;
448 }
449
450 #ifndef VRNA_NR_SAMPLING_HASH
451 /* non-redundant data-structure memorization nodes */
452 NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */
453 NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */
454 #endif
455
456 fbd = 0.; /* stores weight of forbidden terms for given q[ij]*/
457 fbds = 0.; /* stores weight of forbidden term for given motif */
458
459 n = vc->length;
460
461 pf_params = vc->exp_params;
462 md = &(vc->exp_params->model_details);
463 my_iindx = vc->iindx;
464 matrices = vc->exp_matrices;
465
466 hc = vc->hc;
467 if (vc->type == VRNA_FC_TYPE_SINGLE) {
468 n_seq = 1;
469 S1 = vc->sequence_encoding;
470 S2 = vc->sequence_encoding2;
471 S = NULL;
472 S5 = NULL;
473 S3 = NULL;
474 a2s = NULL;
475 } else {
476 n_seq = vc->n_seq;
477 S1 = NULL;
478 S2 = NULL;
479 S = vc->S;
480 S5 = vc->S5;
481 S3 = vc->S3;
482 a2s = vc->a2s;
483 }
484
485 hard_constraints = hc->mx;
486 hc_up_ext = hc->up_ext;
487 sc_wrapper_ext = &(sc_wrap->sc_wrapper_ext);
488
489 /* assume successful backtracing by default */
490 ret = 1;
491
492 q = matrices->q;
493 qb = matrices->qb;
494 q1k = matrices->q1k;
495 qln = matrices->qln;
496 scale = matrices->scale;
497
498 #ifndef VRNA_NR_SAMPLING_HASH
499 if (current_node) {
500 memorized_node_prev = NULL;
501 memorized_node_cur = (*current_node)->head;
502 }
503
504 #endif
505
506 q_temp = 0.;
507
508 #ifdef VRNA_WITH_BOUSTROPHEDON
509 j = init_val;
510 if (j > 1) {
511 /* find j position of first pair */
512 for (; j > 1; j--) {
513 if (hc_up_ext[j]) {
514 if (current_node) {
515 fbd = NR_TOTAL_WEIGHT(*current_node) *
516 q1k[j] /
517 (*q_remain);
518
519 #ifdef USE_FLOAT_PF
520 if (fabsf(NR_TOTAL_WEIGHT(*current_node) - (*q_remain)) / (*q_remain) <= FLT_EPSILON)
521 #else
522 if (fabs(NR_TOTAL_WEIGHT(*current_node) - (*q_remain)) / (*q_remain) <= DBL_EPSILON)
523 #endif
524 /* exhausted ensemble */
525 return 0;
526 }
527
528 r = vrna_urn() * (q1k[j] - fbd);
529 q_temp = q1k[j - 1] * scale[1];
530
531 if (sc_wrapper_ext->red_ext)
532 q_temp *= sc_wrapper_ext->red_ext(1, j, 1, j - 1, sc_wrapper_ext);
533
534 if (current_node) {
535 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_UNPAIRED_SG, j - 1, j) *
536 q1k[j] /
537 (*q_remain);
538 }
539
540 if (r > (q_temp - fbds)) {
541 break; /* j is paired */
542 } else if (current_node) {
543 /* j is unpaired */
544 *q_remain *= q_temp / q1k[j];
545 #ifdef VRNA_NR_SAMPLING_HASH
546 *current_node = add_if_nexists(NRT_UNPAIRED_SG, j - 1, j, *current_node, *q_remain);
547 #else
548 *current_node = add_if_nexists_ll(memory_dat,
549 NRT_UNPAIRED_SG,
550 j - 1,
551 j,
552 memorized_node_prev,
553 memorized_node_cur,
554 *current_node,
555 *q_remain);
556 reset_cursor(&memorized_node_prev, &memorized_node_cur, *current_node); /* resets cursor */
557 #endif
558 }
559 }
560 }
561 if (j <= md->min_loop_size + 1)
562 return ret; /* no more pairs, but still successful */
563
564 #ifndef VRNA_NR_SAMPLING_HASH
565 if (current_node)
566 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_UNPAIRED_SG, j - 1, j);
567
568 #endif
569 /* now find the pairing partner i */
570 if (current_node) {
571 fbd = NR_TOTAL_WEIGHT_TYPE(NRT_EXT_LOOP, *current_node) *
572 q1k[j] /
573 (*q_remain);
574 }
575
576 r = vrna_urn() * (q1k[j] - q_temp - fbd);
577 u = j - 1;
578 i = 2;
579
580 for (qt = 0, k = 1; k < j; k++) {
581 /* apply alternating boustrophedon scheme to variable i */
582 i = (int)(1 + (u - 1) * ((k - 1) % 2)) +
583 (int)((1 - (2 * ((k - 1) % 2))) * ((k - 1) / 2));
584 ij = my_iindx[i] - j;
585 hc_decompose = hard_constraints[n * j + i];
586 if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
587 qkl = qb[ij] *
588 q1k[i - 1];
589
590 if (vc->type == VRNA_FC_TYPE_SINGLE) {
591 type = vrna_get_ptype_md(S2[i], S2[j], md);
592 qkl *= vrna_exp_E_ext_stem(type,
593 (i > 1) ? S1[i - 1] : -1,
594 (j < n) ? S1[j + 1] : -1,
595 pf_params);
596 } else {
597 for (s = 0; s < n_seq; s++) {
598 type = vrna_get_ptype_md(S[s][i], S[s][j], md);
599 qkl *= vrna_exp_E_ext_stem(type,
600 (a2s[s][i] > 1) ? S5[s][i] : -1,
601 (a2s[s][j] < a2s[s][n]) ? S3[s][j] : -1,
602 pf_params);
603 }
604 }
605
606 if ((sc_wrapper_ext->red_stem) && (i == 1))
607 q_temp *= sc_wrapper_ext->red_stem(i, j, i, j, sc_wrapper_ext);
608 else if ((sc_wrapper_ext->split) && (i > 1))
609 q_temp *= sc_wrapper_ext->split(1, j, i, sc_wrapper_ext) *
610 sc_wrapper_ext->red_stem(i, j, i, j, sc_wrapper_ext);
611
612 if (current_node) {
613 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_EXT_LOOP, i, j) *
614 q1k[j] /
615 (*q_remain);
616 qt += qkl - fbds;
617 } else {
618 qt += qkl;
619 }
620
621 if (qt > r) {
622 if (current_node) {
623 *q_remain *= qkl / q1k[j];
624 #ifdef VRNA_NR_SAMPLING_HASH
625 *current_node = add_if_nexists(NRT_EXT_LOOP, i, j, *current_node, *q_remain);
626 #else
627 *current_node = add_if_nexists_ll(memory_dat,
628 NRT_EXT_LOOP,
629 i,
630 j,
631 memorized_node_prev,
632 memorized_node_cur,
633 *current_node,
634 *q_remain);
635 #endif
636 }
637
638 break; /* j is paired */
639 }
640
641 #ifndef VRNA_NR_SAMPLING_HASH
642 if (current_node)
643 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_EXT_LOOP, i, j);
644
645 #endif
646 }
647 }
648 if (k == j) {
649 if (current_node) {
650 /* exhausted ensemble */
651 return 0;
652 } else {
653 vrna_message_warning("backtracking failed in ext loop");
654 /* error */
655 return -1;
656 }
657 }
658
659 backtrack(i, j, pstruc, vc, sc_wrap, nr_mem);
660 j = i - 1;
661 ret = backtrack_ext_loop(j, pstruc, vc, length, sc_wrap, nr_mem);
662 }
663
664 #else
665 int start = init_val;
666 if (start < length) {
667 /* find i position of first pair */
668 for (i = start; i < length; i++) {
669 if (hc_up_ext[i]) {
670 if (current_node) {
671 fbd = NR_TOTAL_WEIGHT(*current_node) *
672 qln[i] /
673 (*q_remain);
674 }
675
676 r = vrna_urn() * (qln[i] - fbd);
677 q_temp = qln[i + 1] * scale[1];
678
679 if (sc_wrapper_ext->red_ext)
680 q_temp *= sc_wrapper_ext->red_ext(i, length, i + 1, length, sc_wrapper_ext);
681
682 if (current_node) {
683 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_UNPAIRED_SG, i, i + 1) *
684 qln[i] /
685 (*q_remain);
686 }
687
688 if (r > (q_temp - fbds)) {
689 break; /* i is paired */
690 } else if (current_node) {
691 *q_remain *= q_temp / qln[i];
692 #ifdef VRNA_NR_SAMPLING_HASH
693 *current_node = add_if_nexists(NRT_UNPAIRED_SG, i, i + 1, *current_node, *q_remain);
694 #else
695 *current_node = add_if_nexists_ll(memory_dat,
696 NRT_UNPAIRED_SG,
697 i,
698 i + 1,
699 memorized_node_prev,
700 memorized_node_cur,
701 *current_node,
702 *q_remain);
703 reset_cursor(&memorized_node_prev, &memorized_node_cur, *current_node); /* resets cursor */
704 #endif
705 }
706 }
707 }
708 if (i >= length)
709 return ret; /* no more pairs, but still successful */
710
711 #ifndef VRNA_NR_SAMPLING_HASH
712 if (current_node)
713 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_UNPAIRED_SG, i, i + 1);
714
715 #endif
716
717 /* now find the pairing partner j */
718 if (current_node) {
719 fbd = NR_TOTAL_WEIGHT_TYPE(NRT_EXT_LOOP, *current_node) *
720 qln[i] /
721 (*q_remain);
722 }
723
724 r = vrna_urn() * (qln[i] - q_temp - fbd);
725 for (qt = 0, j = i + 1; j <= length; j++) {
726 ij = my_iindx[i] - j;
727 hc_decompose = hard_constraints[n * i + j];
728 if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
729 qkl = qb[ij];
730 if (vc->type == VRNA_FC_TYPE_SINGLE) {
731 type = vrna_get_ptype_md(S2[i], S2[j], md);
732 qkl *= vrna_exp_E_ext_stem(type,
733 (i > 1) ? S1[i - 1] : -1,
734 (j < n) ? S1[j + 1] : -1,
735 pf_params);
736 } else {
737 for (s = 0; s < n_seq; s++) {
738 type = vrna_get_ptype_md(S[s][i], S[s][j], md);
739 qkl *= vrna_exp_E_ext_stem(type,
740 (a2s[s][i] > 1) ? S5[s][i] : -1,
741 (a2s[s][j] < a2s[s][n]) ? S3[s][j] : -1,
742 pf_params);
743 }
744 }
745
746 if (j < length) {
747 qkl *= qln[j + 1];
748 if (sc_wrapper_ext->split)
749 qkl *= sc_wrapper_ext->split(i, length, j + 1, sc_wrapper_ext) *
750 sc_wrapper_ext->red_stem(i, j, i, j, sc_wrapper_ext);
751 } else if (sc_wrapper_ext->red_stem) {
752 qkl *= sc_wrapper_ext->red_stem(i, j, i, j, sc_wrapper_ext);
753 }
754
755 if (current_node) {
756 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_EXT_LOOP, i, j) *
757 qln[i] /
758 (*q_remain);
759 qt += qkl - fbds;
760 } else {
761 qt += qkl;
762 }
763
764 if (qt > r) {
765 if (current_node) {
766 *q_remain *= qkl / qln[i];
767 #ifdef VRNA_NR_SAMPLING_HASH
768 *current_node = add_if_nexists(NRT_EXT_LOOP, i, j, *current_node, *q_remain);
769 #else
770 *current_node = add_if_nexists_ll(memory_dat,
771 NRT_EXT_LOOP,
772 i,
773 j,
774 memorized_node_prev,
775 memorized_node_cur,
776 *current_node,
777 *q_remain);
778 #endif
779 }
780
781 break; /* j is paired */
782 }
783
784 #ifndef VRNA_NR_SAMPLING_HASH
785 if (current_node)
786 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_EXT_LOOP, i, j);
787
788 #endif
789 }
790 }
791 if (j == length + 1) {
792 if (current_node) {
793 /* exhausted ensemble */
794 return 0;
795 } else {
796 vrna_message_warning("backtracking failed in ext loop");
797 /* error */
798 return -1;
799 }
800 }
801
802 start = j + 1;
803 ret = backtrack(i, j, pstruc, vc, sc_wrap, nr_mem);
804 if (!ret)
805 return ret;
806
807 ret = backtrack_ext_loop(start, pstruc, vc, length, sc_wrap, nr_mem);
808 }
809
810 #endif
811
812 return ret;
813 }
814
815
816 /* non redundant version of function bactrack_qm */
817 PRIVATE int
backtrack_qm(int i,int j,char * pstruc,vrna_fold_compound_t * vc,struct sc_wrappers * sc_wrap,struct vrna_pbacktrack_memory_s * nr_mem)818 backtrack_qm(int i,
819 int j,
820 char *pstruc,
821 vrna_fold_compound_t *vc,
822 struct sc_wrappers *sc_wrap,
823 struct vrna_pbacktrack_memory_s *nr_mem)
824 {
825 /* divide multiloop into qm and qm1 */
826 int k, u, cnt, span, turn, is_unpaired, *my_iindx, *jindx, *hc_up_ml, ret;
827 FLT_OR_DBL qmt, fbd, fbds, r, q_temp, *qm, *qm1, *expMLbase;
828 double *q_remain;
829 vrna_hc_t *hc;
830 vrna_mx_pf_t *matrices;
831
832 struct sc_mb_exp_dat *sc_wrapper_ml;
833 struct nr_memory **memory_dat;
834
835 NR_NODE **current_node;
836
837 if (nr_mem) {
838 q_remain = &(nr_mem->q_remain);
839 current_node = &(nr_mem->current_node);
840 memory_dat = &(nr_mem->memory_dat);
841 } else {
842 q_remain = NULL;
843 current_node = NULL;
844 memory_dat = NULL;
845 }
846
847 #ifndef VRNA_NR_SAMPLING_HASH
848 /* non-redundant data-structure memorization nodes */
849 NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */
850 NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */
851 #endif
852
853 ret = 1;
854 fbd = 0.; /* stores weight of forbidden terms for given q[ij]*/
855 fbds = 0.; /* stores weight of forbidden term for given motif */
856
857 is_unpaired = 0;
858
859 matrices = vc->exp_matrices;
860 my_iindx = vc->iindx;
861 jindx = vc->jindx;
862
863 hc = vc->hc;
864 hc_up_ml = hc->up_ml;
865 sc_wrapper_ml = &(sc_wrap->sc_wrapper_ml);
866
867 qm = matrices->qm;
868 qm1 = matrices->qm1;
869 expMLbase = matrices->expMLbase;
870
871 turn = vc->exp_params->model_details.min_loop_size;
872
873 #ifndef VRNA_NR_SAMPLING_HASH
874 if (current_node) {
875 memorized_node_prev = NULL;
876 memorized_node_cur = (*current_node)->head;
877 }
878
879 #endif
880
881 if (j > i) {
882 /* now backtrack [i ... j] in qm[] */
883
884 if (current_node) {
885 fbd = NR_TOTAL_WEIGHT(*current_node) *
886 qm[my_iindx[i] - j] /
887 (*q_remain);
888 }
889
890 r = vrna_urn() * (qm[my_iindx[i] - j] - fbd);
891 if (current_node) {
892 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_UNPAIR, i, 0) *
893 qm[my_iindx[i] - j] /
894 (*q_remain);
895 qmt = qm1[jindx[j] + i] - fbds;
896 } else {
897 qmt = qm1[jindx[j] + i];
898 }
899
900 k = cnt = i;
901 q_temp = qm1[jindx[j] + i];
902
903 if (qmt < r) {
904 #ifndef VRNA_NR_SAMPLING_HASH
905 if (current_node)
906 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_UNPAIR, i, 0);
907
908 #endif
909
910 for (span = j - i, cnt = i + 1; cnt <= j; cnt++) {
911 #ifdef VRNA_WITH_BOUSTROPHEDON
912 k = (int)(i + 1 + span * ((cnt - i - 1) % 2)) +
913 (int)((1 - (2 * ((cnt - i - 1) % 2))) * ((cnt - i) / 2));
914 #else
915 k = cnt;
916 #endif
917 q_temp = 0.;
918 u = k - i;
919 /* [i...k] is unpaired */
920 if (hc_up_ml[i] >= u) {
921 q_temp += expMLbase[u] * qm1[jindx[j] + k];
922
923 if (sc_wrapper_ml->red_ml)
924 q_temp *= sc_wrapper_ml->red_ml(i, j, k, j, sc_wrapper_ml);
925
926 if (current_node) {
927 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_UNPAIR, k, 0) *
928 qm[my_iindx[i] - j] /
929 (*q_remain);
930 qmt += q_temp - fbds;
931 } else {
932 qmt += q_temp;
933 }
934 }
935
936 if (qmt >= r) {
937 /* we have chosen unpaired version */
938 is_unpaired = 1;
939 break;
940 }
941
942 #ifndef VRNA_NR_SAMPLING_HASH
943 if (current_node)
944 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_UNPAIR, k, 0);
945
946 #endif
947
948 /* split between k-1, k */
949 q_temp = qm[my_iindx[i] - (k - 1)] *
950 qm1[jindx[j] + k];
951
952 if (sc_wrapper_ml->decomp_ml)
953 q_temp *= sc_wrapper_ml->decomp_ml(i, j, k - 1, k, sc_wrapper_ml);
954
955 if (current_node) {
956 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM_PAIR, k, 0) *
957 qm[my_iindx[i] - j] /
958 (*q_remain);
959 qmt += q_temp - fbds;
960 } else {
961 qmt += q_temp;
962 }
963
964 if (qmt >= r)
965 break;
966
967 #ifndef VRNA_NR_SAMPLING_HASH
968 if (current_node)
969 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM_PAIR, k, 0);
970
971 #endif
972 }
973 } else {
974 is_unpaired = 1;
975 }
976
977 if (current_node) {
978 *q_remain *= q_temp / qm[my_iindx[i] - j];
979 #ifdef VRNA_NR_SAMPLING_HASH
980 if (is_unpaired)
981 *current_node = add_if_nexists(NRT_QM_UNPAIR, k, 0, *current_node, *q_remain);
982 else
983 *current_node = add_if_nexists(NRT_QM_PAIR, k, 0, *current_node, *q_remain);
984
985 #else
986 if (is_unpaired)
987 *current_node = add_if_nexists_ll(memory_dat,
988 NRT_QM_UNPAIR,
989 k,
990 0,
991 memorized_node_prev,
992 memorized_node_cur,
993 *current_node,
994 *q_remain);
995 else
996 *current_node = add_if_nexists_ll(memory_dat,
997 NRT_QM_PAIR,
998 k,
999 0,
1000 memorized_node_prev,
1001 memorized_node_cur,
1002 *current_node,
1003 *q_remain);
1004
1005 #endif
1006 }
1007
1008 if (cnt > j)
1009 return 0;
1010
1011 ret = backtrack_qm1(k, j, pstruc, vc, sc_wrap, nr_mem);
1012
1013 if (ret == 0)
1014 return ret;
1015
1016 if (k < i + turn)
1017 return ret; /* no more pairs */
1018
1019 if (!is_unpaired) {
1020 /* if we've chosen creating a branch in [i..k-1] */
1021 ret = backtrack_qm(i, k - 1, pstruc, vc, sc_wrap, nr_mem);
1022
1023 if (ret == 0)
1024 return ret;
1025 }
1026 }
1027
1028 return ret;
1029 }
1030
1031
1032 PRIVATE int
backtrack_qm1(int i,int j,char * pstruc,vrna_fold_compound_t * vc,struct sc_wrappers * sc_wrap,struct vrna_pbacktrack_memory_s * nr_mem)1033 backtrack_qm1(int i,
1034 int j,
1035 char *pstruc,
1036 vrna_fold_compound_t *vc,
1037 struct sc_wrappers *sc_wrap,
1038 struct vrna_pbacktrack_memory_s *nr_mem)
1039 {
1040 /* i is paired to l, i<l<j; backtrack in qm1 to find l */
1041 unsigned char *hard_constraints;
1042 char *ptype;
1043 short *S1, **S, **S5, **S3;
1044 unsigned int n, s, n_seq;
1045 int ii, l, il, type, turn, u, *my_iindx, *jindx, *hc_up_ml;
1046 FLT_OR_DBL qt, fbd, fbds, r, q_temp, *qm1, *qb, *expMLbase;
1047 double *q_remain;
1048 vrna_exp_param_t *pf_params;
1049 vrna_md_t *md;
1050 vrna_hc_t *hc;
1051 vrna_mx_pf_t *matrices;
1052
1053 struct nr_memory **memory_dat;
1054 struct sc_mb_exp_dat *sc_wrapper_ml;
1055
1056 NR_NODE **current_node;
1057
1058 if (nr_mem) {
1059 q_remain = &(nr_mem->q_remain);
1060 current_node = &(nr_mem->current_node);
1061 memory_dat = &(nr_mem->memory_dat);
1062 } else {
1063 q_remain = NULL;
1064 current_node = NULL;
1065 memory_dat = NULL;
1066 }
1067
1068 #ifndef VRNA_NR_SAMPLING_HASH
1069 /* non-redundant data-structure memorization nodes */
1070 NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */
1071 NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */
1072 #endif
1073
1074 n = vc->length;
1075 fbd = 0.;
1076 fbds = 0.;
1077 pf_params = vc->exp_params;
1078 md = &(pf_params->model_details);
1079 my_iindx = vc->iindx;
1080 jindx = vc->jindx;
1081 hc = vc->hc;
1082 hc_up_ml = hc->up_ml;
1083 hard_constraints = hc->mx;
1084 sc_wrapper_ml = &(sc_wrap->sc_wrapper_ml);
1085
1086 matrices = vc->exp_matrices;
1087 qb = matrices->qb;
1088 qm1 = matrices->qm1;
1089 expMLbase = matrices->expMLbase;
1090 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1091 n_seq = 1;
1092 ptype = vc->ptype;
1093 S1 = vc->sequence_encoding;
1094 S = NULL;
1095 S5 = NULL;
1096 S3 = NULL;
1097 } else {
1098 n_seq = vc->n_seq;
1099 ptype = NULL;
1100 S1 = NULL;
1101 S = vc->S;
1102 S5 = vc->S5;
1103 S3 = vc->S3;
1104 }
1105
1106 turn = pf_params->model_details.min_loop_size;
1107
1108 #ifndef VRNA_NR_SAMPLING_HASH
1109 if (current_node) {
1110 memorized_node_prev = NULL;
1111 memorized_node_cur = (*current_node)->head;
1112 }
1113
1114 #endif
1115
1116 if (current_node) {
1117 fbd = NR_TOTAL_WEIGHT(*current_node) *
1118 qm1[jindx[j] + i] /
1119 (*q_remain);
1120 }
1121
1122 r = vrna_urn() * (qm1[jindx[j] + i] - fbd);
1123 ii = my_iindx[i];
1124 for (qt = 0., l = j; l > i + turn; l--) {
1125 il = jindx[l] + i;
1126 if (hard_constraints[n * i + l] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) {
1127 u = j - l;
1128 if (hc_up_ml[l + 1] >= u) {
1129 q_temp = qb[ii - l] *
1130 expMLbase[j - l];
1131
1132 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1133 type = vrna_get_ptype(il, ptype);
1134 q_temp *= exp_E_MLstem(type, S1[i - 1], S1[l + 1], pf_params);
1135 } else {
1136 for (s = 0; s < n_seq; s++) {
1137 type = vrna_get_ptype_md(S[s][i], S[s][l], md);
1138 q_temp *= exp_E_MLstem(type, S5[s][i], S3[s][l], pf_params);
1139 }
1140 }
1141
1142 if (sc_wrapper_ml->red_stem)
1143 q_temp *= sc_wrapper_ml->red_stem(i, j, i, l, sc_wrapper_ml);
1144
1145 if (current_node) {
1146 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_QM1_BRANCH, i, l) *
1147 qm1[jindx[j] + i] /
1148 (*q_remain);
1149 qt += q_temp - fbds;
1150 } else {
1151 qt += q_temp;
1152 }
1153
1154 if (qt >= r) {
1155 if (current_node) {
1156 *q_remain *= q_temp / qm1[jindx[j] + i];
1157 #ifdef VRNA_NR_SAMPLING_HASH
1158 *current_node = add_if_nexists(NRT_QM1_BRANCH, i, l, *current_node, *q_remain);
1159 #else
1160 *current_node = add_if_nexists_ll(memory_dat,
1161 NRT_QM1_BRANCH,
1162 i,
1163 l,
1164 memorized_node_prev,
1165 memorized_node_cur,
1166 *current_node,
1167 *q_remain);
1168 #endif
1169 }
1170
1171 break;
1172 }
1173
1174 #ifndef VRNA_NR_SAMPLING_HASH
1175 if (current_node)
1176 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_QM1_BRANCH, i, l);
1177
1178 #endif
1179 } else {
1180 l = i + turn;
1181 break;
1182 }
1183 }
1184 }
1185 if (l < i + turn + 1) {
1186 if (current_node) {
1187 return 0;
1188 } else {
1189 vrna_message_error("backtrack failed in qm1");
1190 return 0;
1191 }
1192 }
1193
1194 return backtrack(i, l, pstruc, vc, sc_wrap, nr_mem);
1195 }
1196
1197
1198 PRIVATE void
backtrack_qm2(int k,int n,char * pstruc,vrna_fold_compound_t * vc,struct sc_wrappers * sc_wrap)1199 backtrack_qm2(int k,
1200 int n,
1201 char *pstruc,
1202 vrna_fold_compound_t *vc,
1203 struct sc_wrappers *sc_wrap)
1204 {
1205 int u, turn, *jindx;
1206 FLT_OR_DBL qom2t, r, *qm1, *qm2;
1207 struct sc_mb_exp_dat *sc_wrapper_ml;
1208
1209 jindx = vc->jindx;
1210 qm1 = vc->exp_matrices->qm1;
1211 qm2 = vc->exp_matrices->qm2;
1212 turn = vc->exp_params->model_details.min_loop_size;
1213 sc_wrapper_ml = &(sc_wrap->sc_wrapper_ml);
1214
1215 r = vrna_urn() * qm2[k];
1216 /* we have to search for our barrier u between qm1 and qm1 */
1217 if (sc_wrapper_ml->decomp_ml) {
1218 for (qom2t = 0., u = k + turn + 1; u < n - turn - 1; u++) {
1219 qom2t += qm1[jindx[u] + k] *
1220 qm1[jindx[n] + (u + 1)] *
1221 sc_wrapper_ml->decomp_ml(k, n, u, u + 1, sc_wrapper_ml);
1222
1223 if (qom2t > r)
1224 break;
1225 }
1226 } else {
1227 for (qom2t = 0., u = k + turn + 1; u < n - turn - 1; u++) {
1228 qom2t += qm1[jindx[u] + k] * qm1[jindx[n] + (u + 1)];
1229 if (qom2t > r)
1230 break;
1231 }
1232 }
1233
1234 if (u == n - turn)
1235 vrna_message_error("backtrack failed in qm2");
1236
1237 backtrack_qm1(k, u, pstruc, vc, sc_wrap, NULL);
1238 backtrack_qm1(u + 1, n, pstruc, vc, sc_wrap, NULL);
1239 }
1240
1241
1242 PRIVATE int
backtrack(int i,int j,char * pstruc,vrna_fold_compound_t * vc,struct sc_wrappers * sc_wrap,struct vrna_pbacktrack_memory_s * nr_mem)1243 backtrack(int i,
1244 int j,
1245 char *pstruc,
1246 vrna_fold_compound_t *vc,
1247 struct sc_wrappers *sc_wrap,
1248 struct vrna_pbacktrack_memory_s *nr_mem)
1249 {
1250 unsigned char *hard_constraints, hc_decompose;
1251 char *ptype;
1252 short *S1, **S, **S5, **S3;
1253 unsigned int **a2s, s, n_seq, n, type, type_2, *types, u1_local, u2_local;
1254 int *my_iindx, *jindx, *hc_up_int, ret, *pscore, turn, *rtype,
1255 k, l, kl, u1, u2, max_k, min_l, ii, jj;
1256 FLT_OR_DBL *qb, *qm, *qm1, *scale, r, fbd, fbds, qbt1, qbr, qt, q_temp,
1257 kTn, closingPair, expMLclosing;
1258 double *q_remain;
1259 vrna_mx_pf_t *matrices;
1260 vrna_exp_param_t *pf_params;
1261 vrna_md_t *md;
1262 vrna_hc_t *hc;
1263
1264 struct nr_memory **memory_dat;
1265 struct sc_int_exp_dat *sc_wrapper_int;
1266 struct sc_mb_exp_dat *sc_wrapper_ml;
1267
1268 NR_NODE **current_node;
1269
1270 if (nr_mem) {
1271 q_remain = &(nr_mem->q_remain);
1272 current_node = &(nr_mem->current_node);
1273 memory_dat = &(nr_mem->memory_dat);
1274 } else {
1275 q_remain = NULL;
1276 current_node = NULL;
1277 memory_dat = NULL;
1278 }
1279
1280 #ifndef VRNA_NR_SAMPLING_HASH
1281 /* non-redundant data-structure memorization nodes */
1282 NR_NODE *memorized_node_prev = NULL; /* remembers previous-to-current node in linked list */
1283 NR_NODE *memorized_node_cur = NULL; /* remembers actual node in linked list */
1284 #endif
1285
1286 ret = 1; /* default is success */
1287 fbd = 0.; /* stores weight of forbidden terms for given q[ij] */
1288 fbds = 0.; /* stores weight of forbidden term for given motif */
1289 qbt1 = 0.;
1290 q_temp = 0.;
1291
1292 n = vc->length;
1293 pf_params = vc->exp_params;
1294 kTn = pf_params->kT / 10.;
1295 md = &(pf_params->model_details);
1296 my_iindx = vc->iindx;
1297 jindx = vc->jindx;
1298 turn = pf_params->model_details.min_loop_size;
1299 rtype = &(pf_params->model_details.rtype[0]);
1300
1301 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1302 n_seq = 1;
1303 ptype = vc->ptype;
1304 types = NULL;
1305 pscore = NULL;
1306 S1 = vc->sequence_encoding;
1307 S = NULL;
1308 S5 = NULL;
1309 S3 = NULL;
1310 a2s = NULL;
1311 expMLclosing = pf_params->expMLclosing;
1312 } else {
1313 n_seq = vc->n_seq;
1314 ptype = NULL;
1315 types = (unsigned int *)vrna_alloc(sizeof(unsigned int) * n_seq);
1316 pscore = vc->pscore;
1317 S1 = NULL;
1318 S = vc->S;
1319 S5 = vc->S5;
1320 S3 = vc->S3;
1321 a2s = vc->a2s;
1322 expMLclosing = pow(pf_params->expMLclosing, (double)n_seq);
1323 }
1324
1325 hc = vc->hc;
1326 hc_up_int = hc->up_int;
1327 hard_constraints = hc->mx;
1328 sc_wrapper_int = &(sc_wrap->sc_wrapper_int);
1329 sc_wrapper_ml = &(sc_wrap->sc_wrapper_ml);
1330
1331 matrices = vc->exp_matrices;
1332 qb = matrices->qb;
1333 qm = matrices->qm;
1334 qm1 = matrices->qm1;
1335 scale = matrices->scale;
1336
1337 #ifndef VRNA_NR_SAMPLING_HASH
1338 if (current_node) {
1339 memorized_node_prev = NULL;
1340 memorized_node_cur = (*current_node)->head;
1341 }
1342
1343 #endif
1344
1345 hc_decompose = hard_constraints[n * j + i];
1346
1347 do {
1348 k = i;
1349 l = j;
1350
1351 qbr = qb[my_iindx[i] - j];
1352
1353 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1354 type = vrna_get_ptype(jindx[j] + i, ptype);
1355 } else {
1356 qbr /= exp(pscore[jindx[j] + i] / kTn);
1357 for (s = 0; s < n_seq; s++)
1358 types[s] = vrna_get_ptype_md(S[s][i], S[s][j], md);
1359 }
1360
1361 if (current_node)
1362 fbd = NR_TOTAL_WEIGHT(*current_node) * qbr / (*q_remain);
1363
1364 pstruc[i - 1] = '(';
1365 pstruc[j - 1] = ')';
1366
1367 r = vrna_urn() * (qbr - fbd);
1368 qbt1 = 0.;
1369
1370 hc_decompose = hard_constraints[n * i + j];
1371
1372 /* hairpin contribution */
1373 q_temp = vrna_exp_E_hp_loop(vc, i, j);
1374
1375 if (current_node) {
1376 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_HAIRPIN, 0, 0) *
1377 qbr /
1378 (*q_remain);
1379 qbt1 += (q_temp - fbds);
1380 } else {
1381 qbt1 += q_temp;
1382 }
1383
1384 if (qbt1 >= r) {
1385 /* found the hairpin we're done */
1386 if (current_node) {
1387 *q_remain *= q_temp / qbr;
1388 #ifdef VRNA_NR_SAMPLING_HASH
1389 *current_node = add_if_nexists(NRT_HAIRPIN, 0, 0, *current_node, *q_remain);
1390 #else
1391 *current_node = add_if_nexists_ll(memory_dat,
1392 NRT_HAIRPIN,
1393 0,
1394 0,
1395 memorized_node_prev,
1396 memorized_node_cur,
1397 *current_node,
1398 *q_remain);
1399 #endif
1400 }
1401
1402 free(types);
1403
1404 return ret;
1405 }
1406
1407 #ifndef VRNA_NR_SAMPLING_HASH
1408 if (current_node)
1409 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_HAIRPIN, 0, 0);
1410
1411 #endif
1412
1413 if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
1414 /* interior loop contributions */
1415 max_k = i + MAXLOOP + 1;
1416 max_k = MIN2(max_k, j - turn - 2);
1417 max_k = MIN2(max_k, i + 1 + hc_up_int[i + 1]);
1418 for (k = i + 1; k <= max_k; k++) {
1419 u1 = k - i - 1;
1420 min_l = MAX2(k + turn + 1, j - 1 - MAXLOOP + u1);
1421 kl = my_iindx[k] - j + 1;
1422 for (u2 = 0, l = j - 1; l >= min_l; l--, kl++, u2++) {
1423 if (hc_up_int[l + 1] < u2)
1424 break;
1425
1426 if (hard_constraints[n * k + l] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC) {
1427 q_temp = qb[kl]
1428 * scale[u1 + u2 + 2];
1429
1430 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1431 type_2 = rtype[vrna_get_ptype(jindx[l] + k, ptype)];
1432
1433 /* add *scale[u1+u2+2] */
1434 q_temp *= exp_E_IntLoop(u1,
1435 u2,
1436 type,
1437 type_2,
1438 S1[i + 1],
1439 S1[j - 1],
1440 S1[k - 1],
1441 S1[l + 1],
1442 pf_params);
1443 } else {
1444 for (s = 0; s < n_seq; s++) {
1445 u1_local = a2s[s][k - 1] - a2s[s][i] /*??*/;
1446 u2_local = a2s[s][j - 1] - a2s[s][l];
1447 type_2 = vrna_get_ptype_md(S[s][l], S[s][k], md);
1448 q_temp *= exp_E_IntLoop(u1_local,
1449 u2_local,
1450 types[s],
1451 type_2,
1452 S3[s][i],
1453 S5[s][j],
1454 S5[s][k],
1455 S3[s][l],
1456 pf_params);
1457 }
1458 }
1459
1460 if (sc_wrapper_int->pair)
1461 q_temp *= sc_wrapper_int->pair(i, j, k, l, sc_wrapper_int);
1462
1463 if (current_node) {
1464 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_IT_LOOP, k, l) *
1465 qbr /
1466 (*q_remain);
1467 qbt1 += q_temp - fbds;
1468 } else {
1469 qbt1 += q_temp;
1470 }
1471
1472 if (qbt1 >= r)
1473 break;
1474
1475 #ifndef VRNA_NR_SAMPLING_HASH
1476 if (current_node)
1477 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_IT_LOOP, k, l);
1478
1479 #endif
1480 }
1481 }
1482 if (qbt1 >= r)
1483 break;
1484 }
1485 if (k <= max_k) {
1486 if (current_node) {
1487 *q_remain *= q_temp / qbr;
1488 #ifdef VRNA_NR_SAMPLING_HASH
1489 *current_node = add_if_nexists(NRT_IT_LOOP, k, l, *current_node, *q_remain);
1490 #else
1491 *current_node = add_if_nexists_ll(memory_dat,
1492 NRT_IT_LOOP,
1493 k,
1494 l,
1495 memorized_node_prev,
1496 memorized_node_cur,
1497 *current_node,
1498 *q_remain);
1499 #endif
1500 }
1501
1502 free(types);
1503
1504 return backtrack(k, l, pstruc, vc, sc_wrap, nr_mem); /* found the interior loop, repeat for inside */
1505 } else {
1506 /* interior loop contributions did not exceed threshold, so we break */
1507 break;
1508 }
1509 } else {
1510 /* must not be interior loop, so we break out */
1511 break;
1512 }
1513 } while (1);
1514
1515 /* backtrack in multi-loop */
1516 if (hard_constraints[n * j + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1517 closingPair = expMLclosing *
1518 scale[2];
1519
1520 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1521 type = rtype[vrna_get_ptype(jindx[j] + i, ptype)];
1522 closingPair *= exp_E_MLstem(type, S1[j - 1], S1[i + 1], pf_params);
1523 } else {
1524 for (s = 0; s < n_seq; s++) {
1525 type = vrna_get_ptype_md(S[s][j], S[s][i], md);
1526 closingPair *= exp_E_MLstem(type, S5[s][j], S3[s][i], pf_params);
1527 }
1528 }
1529
1530 if (sc_wrapper_ml->pair)
1531 closingPair *= sc_wrapper_ml->pair(i, j, sc_wrapper_ml);
1532
1533 i++;
1534 j--;
1535 /* find the first split index */
1536 ii = my_iindx[i]; /* ii-j=[i,j] */
1537 jj = jindx[j]; /* jj+i=[j,i] */
1538
1539 if (sc_wrapper_ml->decomp_ml) {
1540 for (qt = qbt1, k = i + 1; k < j; k++) {
1541 q_temp = qm[ii - (k - 1)] *
1542 qm1[jj + k] *
1543 closingPair *
1544 sc_wrapper_ml->decomp_ml(i,
1545 j,
1546 k - 1,
1547 k,
1548 sc_wrapper_ml);
1549
1550
1551 if (current_node) {
1552 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_MT_LOOP, k, 0) *
1553 qbr /
1554 (*q_remain);
1555 qt += q_temp - fbds;
1556 qbt1 += q_temp - fbds;
1557 } else {
1558 qt += q_temp;
1559 qbt1 += q_temp;
1560 }
1561
1562 if (qt >= r)
1563 break;
1564
1565 #ifndef VRNA_NR_SAMPLING_HASH
1566 if (current_node)
1567 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_MT_LOOP, k, 0);
1568
1569 #endif
1570 }
1571 } else {
1572 for (qt = qbt1, k = i + 1; k < j; k++) {
1573 q_temp = qm[ii - (k - 1)] *
1574 qm1[jj + k] *
1575 closingPair;
1576
1577 if (current_node) {
1578 fbds = NR_GET_WEIGHT(*current_node, memorized_node_cur, NRT_MT_LOOP, k, 0) *
1579 qbr /
1580 (*q_remain);
1581 qt += q_temp - fbds;
1582 qbt1 += q_temp - fbds;
1583 } else {
1584 qt += q_temp;
1585 qbt1 += q_temp;
1586 }
1587
1588 if (qt >= r)
1589 break;
1590
1591 #ifndef VRNA_NR_SAMPLING_HASH
1592 if (current_node)
1593 advance_cursor(&memorized_node_prev, &memorized_node_cur, NRT_MT_LOOP, k, 0);
1594
1595 #endif
1596 }
1597 }
1598
1599 if (k >= j) {
1600 if (current_node) {
1601 free(types);
1602 return 0; /* backtrack failed for non-redundant mode most likely due to numerical instabilities */
1603 } else {
1604 vrna_message_error("backtrack failed, can't find split index ");
1605 }
1606 }
1607
1608 if (current_node) {
1609 *q_remain *= q_temp / qbr;
1610 #ifdef VRNA_NR_SAMPLING_HASH
1611 *current_node = add_if_nexists(NRT_MT_LOOP, k, 0, *current_node, *q_remain);
1612 #else
1613 *current_node = add_if_nexists_ll(memory_dat,
1614 NRT_MT_LOOP,
1615 k,
1616 0,
1617 memorized_node_prev,
1618 memorized_node_cur,
1619 *current_node,
1620 *q_remain);
1621 #endif
1622 }
1623
1624 ret = backtrack_qm1(k, j, pstruc, vc, sc_wrap, nr_mem);
1625
1626 if (ret == 0) {
1627 free(types);
1628 return ret;
1629 }
1630
1631 j = k - 1;
1632
1633 ret = backtrack_qm(i, j, pstruc, vc, sc_wrap, nr_mem);
1634 }
1635
1636 free(types);
1637
1638 return ret;
1639 }
1640
1641
1642 PRIVATE unsigned int
pbacktrack_circ(vrna_fold_compound_t * vc,unsigned int num_samples,vrna_boltzmann_sampling_callback * bs_cb,void * data)1643 pbacktrack_circ(vrna_fold_compound_t *vc,
1644 unsigned int num_samples,
1645 vrna_boltzmann_sampling_callback *bs_cb,
1646 void *data)
1647 {
1648 unsigned char *hc_mx, eval_loop;
1649 char *pstruc;
1650 short *S1, *S2, **S, **S5, **S3;
1651 unsigned int type, type2, *tt, s, n_seq, **a2s, u1_local,
1652 u2_local, u3_local, count;
1653 int i, j, k, l, n, u, *hc_up, *my_iindx, turn,
1654 ln1, ln2, ln3, lstart;
1655 FLT_OR_DBL r, qt, q_temp, qo, qmo, *scale, *qb, *qm, *qm2,
1656 qb_ij, expMLclosing;
1657 vrna_exp_param_t *pf_params;
1658 vrna_md_t *md;
1659 vrna_mx_pf_t *matrices;
1660 struct sc_wrappers *sc_wrap;
1661 struct sc_ext_exp_dat *sc_wrapper_ext;
1662 struct sc_int_exp_dat *sc_wrapper_int;
1663 struct sc_mb_exp_dat *sc_wrapper_ml;
1664
1665 n = vc->length;
1666 pf_params = vc->exp_params;
1667 md = &(pf_params->model_details);
1668 matrices = vc->exp_matrices;
1669 my_iindx = vc->iindx;
1670 expMLclosing = pf_params->expMLclosing;
1671 turn = pf_params->model_details.min_loop_size;
1672
1673 qo = matrices->qo;
1674 qmo = matrices->qmo;
1675 qb = matrices->qb;
1676 qm = matrices->qm;
1677 qm2 = matrices->qm2;
1678 scale = matrices->scale;
1679
1680 hc_mx = vc->hc->mx;
1681 hc_up = vc->hc->up_int;
1682
1683 sc_wrap = sc_init(vc);
1684 sc_wrapper_ext = &(sc_wrap->sc_wrapper_ext);
1685 sc_wrapper_int = &(sc_wrap->sc_wrapper_int);
1686 sc_wrapper_ml = &(sc_wrap->sc_wrapper_ml);
1687
1688 if (vc->type == VRNA_FC_TYPE_SINGLE) {
1689 n_seq = 1;
1690 tt = NULL;
1691 S1 = vc->sequence_encoding;
1692 S2 = vc->sequence_encoding2;
1693 S = NULL;
1694 S5 = NULL;
1695 S3 = NULL;
1696 a2s = NULL;
1697 expMLclosing = pf_params->expMLclosing;
1698 } else {
1699 n_seq = vc->n_seq;
1700 tt = (unsigned int *)vrna_alloc(sizeof(unsigned int) * n_seq);
1701 S1 = NULL;
1702 S2 = NULL;
1703 S = vc->S;
1704 S5 = vc->S5;
1705 S3 = vc->S3;
1706 a2s = vc->a2s;
1707 expMLclosing = pow(pf_params->expMLclosing, (double)n_seq);
1708 }
1709
1710 for (count = 0; count < num_samples; count++) {
1711 pstruc = vrna_alloc((n + 1) * sizeof(char));
1712
1713 /* initialize pstruct with single bases */
1714 memset(pstruc, '.', sizeof(char) * n);
1715
1716 qt = 1.0 * scale[n];
1717
1718 /* add soft constraints for open chain configuration */
1719 if (sc_wrapper_ext->red_up)
1720 qt *= sc_wrapper_ext->red_up(1, n, sc_wrapper_ext);
1721
1722 r = vrna_urn() * qo;
1723
1724 /* open chain? */
1725 if (qt > r)
1726 goto pbacktrack_circ_loop_end;
1727
1728 for (i = 1; (i < n); i++) {
1729 for (j = i + turn + 1; (j <= n); j++) {
1730 u = n - j + i - 1;
1731
1732 if (u < turn)
1733 continue;
1734
1735 qb_ij = qb[my_iindx[i] - j];
1736
1737 qt += qb_ij *
1738 vrna_exp_E_hp_loop(vc, j, i);
1739
1740 /* found a hairpin? so backtrack in the enclosed part and we're done */
1741 if (qt > r) {
1742 backtrack(i, j, pstruc, vc, sc_wrap, NULL);
1743 goto pbacktrack_circ_loop_end;
1744 }
1745
1746 /* 2. search for (k,l) with which we can close an interior loop */
1747 if (hc_mx[n * i + j] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
1748 if (vc->type == VRNA_FC_TYPE_SINGLE)
1749 type = vrna_get_ptype_md(S2[j], S2[i], md);
1750 else
1751 for (s = 0; s < n_seq; s++)
1752 tt[s] = vrna_get_ptype_md(S[s][j], S[s][i], md);
1753
1754 for (k = j + 1; (k < n); k++) {
1755 ln1 = k - j - 1;
1756 if (ln1 + i - 1 > MAXLOOP)
1757 break;
1758
1759 if (hc_up[j + 1] < ln1)
1760 break;
1761
1762 lstart = ln1 + i - 1 + n - MAXLOOP;
1763 if (lstart < k + turn + 1)
1764 lstart = k + turn + 1;
1765
1766 for (l = lstart; (l <= n); l++) {
1767 ln2 = (i - 1);
1768 ln3 = (n - l);
1769
1770 if (hc_up[l + 1] < (ln2 + ln3))
1771 continue;
1772
1773 if ((ln1 + ln2 + ln3) > MAXLOOP)
1774 continue;
1775
1776 eval_loop = hc_mx[n * k + l] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP;
1777
1778 if (eval_loop) {
1779 q_temp = qb_ij *
1780 qb[my_iindx[k] - l] *
1781 scale[ln1 + ln2 + ln3];
1782
1783 switch (vc->type) {
1784 case VRNA_FC_TYPE_SINGLE:
1785 type2 = vrna_get_ptype_md(S2[l], S2[k], md);
1786 q_temp *= exp_E_IntLoop(ln2 + ln3,
1787 ln1,
1788 type2,
1789 type,
1790 S1[l + 1],
1791 S1[k - 1],
1792 S1[i - 1],
1793 S1[j + 1],
1794 pf_params);
1795 break;
1796 case VRNA_FC_TYPE_COMPARATIVE:
1797 for (s = 0; s < n_seq; s++) {
1798 type2 = vrna_get_ptype_md(S[s][l], S[s][k], md);
1799 u1_local = a2s[s][i - 1];
1800 u2_local = a2s[s][k - 1] - a2s[s][j];
1801 u3_local = a2s[s][n] - a2s[s][l];
1802 q_temp *= exp_E_IntLoop(u1_local + u3_local,
1803 u2_local,
1804 type2,
1805 tt[s],
1806 S3[s][l],
1807 S5[s][k],
1808 S5[s][i],
1809 S3[s][j],
1810 pf_params);
1811 }
1812 break;
1813 }
1814
1815 if (sc_wrapper_int->pair_ext)
1816 q_temp *= sc_wrapper_int->pair_ext(i, j, k, l, sc_wrapper_int);
1817
1818 qt += q_temp;
1819 /*
1820 * found an exterior interior loop? also this time, we can go straight
1821 * forward and backtracking the both enclosed parts and we're done
1822 */
1823 if (qt > r) {
1824 backtrack(i, j, pstruc, vc, sc_wrap, NULL);
1825 backtrack(k, l, pstruc, vc, sc_wrap, NULL);
1826 goto pbacktrack_circ_loop_end;
1827 }
1828 }
1829 }
1830 } /* end of kl double loop */
1831 }
1832 }
1833 } /* end of ij double loop */
1834 {
1835 /* as we reach this part, we have to search for our barrier between qm and qm2 */
1836 qt = 0.;
1837 r = vrna_urn() * qmo;
1838 if (sc_wrapper_ml->decomp_ml) {
1839 for (k = turn + 2; k < n - 2 * turn - 3; k++) {
1840 qt += qm[my_iindx[1] - k] *
1841 qm2[k + 1] *
1842 expMLclosing *
1843 sc_wrapper_ml->decomp_ml(1,
1844 n,
1845 k,
1846 k + 1,
1847 sc_wrapper_ml);
1848
1849
1850 /* backtrack in qm and qm2 if we've found a valid barrier k */
1851 if (qt > r) {
1852 backtrack_qm(1, k, pstruc, vc, sc_wrap, NULL);
1853 backtrack_qm2(k + 1, n, pstruc, vc, sc_wrap);
1854 goto pbacktrack_circ_loop_end;
1855 }
1856 }
1857 } else {
1858 for (k = turn + 2; k < n - 2 * turn - 3; k++) {
1859 qt += qm[my_iindx[1] - k] *
1860 qm2[k + 1] *
1861 expMLclosing;
1862 /* backtrack in qm and qm2 if we've found a valid barrier k */
1863 if (qt > r) {
1864 backtrack_qm(1, k, pstruc, vc, sc_wrap, NULL);
1865 backtrack_qm2(k + 1, n, pstruc, vc, sc_wrap);
1866 goto pbacktrack_circ_loop_end;
1867 }
1868 }
1869 }
1870 }
1871 /*
1872 * if we reach the actual end of this function, an error has occured
1873 * cause we HAVE TO find an exterior loop or an open chain!!!
1874 */
1875 vrna_message_error("backtracking failed in exterior loop");
1876
1877 pbacktrack_circ_loop_end:
1878
1879 if (bs_cb)
1880 bs_cb(pstruc, data);
1881
1882 free(pstruc);
1883 }
1884
1885 sc_free(sc_wrap);
1886
1887 return count;
1888 }
1889