1 /*
2 * minimum free energy
3 * RNA secondary structure prediction
4 *
5 * c Ivo Hofacker, Chrisoph Flamm
6 * original implementation by
7 * Walter Fontana
8 * g-quadruplex support and threadsafety
9 * by Ronny Lorenz
10 *
11 * Vienna RNA package
12 */
13
14 #ifdef HAVE_CONFIG_H
15 #include "config.h"
16 #endif
17
18 #include <stdio.h>
19 #include <stdlib.h>
20 #include <math.h>
21 #include <ctype.h>
22 #include <string.h>
23 #include <limits.h>
24
25 #include "ViennaRNA/utils/basic.h"
26 #include "ViennaRNA/utils/structures.h"
27 #include "ViennaRNA/params/default.h"
28 #include "ViennaRNA/datastructures/basic.h"
29 #include "ViennaRNA/fold_vars.h"
30 #include "ViennaRNA/params/basic.h"
31 #include "ViennaRNA/constraints/hard.h"
32 #include "ViennaRNA/constraints/soft.h"
33 #include "ViennaRNA/gquad.h"
34 #include "ViennaRNA/structured_domains.h"
35 #include "ViennaRNA/unstructured_domains.h"
36 #include "ViennaRNA/loops/all.h"
37 #include "ViennaRNA/alphabet.h"
38 #include "ViennaRNA/mfe.h"
39
40 #ifdef __GNUC__
41 # define INLINE inline
42 #else
43 # define INLINE
44 #endif
45
46 #define MAXSECTORS 500 /* dimension for a backtrack array */
47
48 struct aux_arrays {
49 int *cc; /* auxilary arrays for canonical structures */
50 int *cc1; /* auxilary arrays for canonical structures */
51 int *Fmi; /* holds row i of fML (avoids jumps in memory) */
52 int *DMLi; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
53 int *DMLi1; /* MIN(fML[i+1,k]+fML[k+1,j]) */
54 int *DMLi2; /* MIN(fML[i+2,k]+fML[k+1,j]) */
55 };
56
57
58 /*
59 #################################
60 # GLOBAL VARIABLES #
61 #################################
62 */
63
64 /*
65 #################################
66 # PRIVATE VARIABLES #
67 #################################
68 */
69
70 /*
71 #################################
72 # PRIVATE FUNCTION DECLARATIONS #
73 #################################
74 */
75
76 PRIVATE int
77 fill_arrays(vrna_fold_compound_t *fc);
78
79
80 PRIVATE int
81 postprocess_circular(vrna_fold_compound_t *fc,
82 sect bt_stack[],
83 int *bt);
84
85
86 PRIVATE INLINE void
87 fill_fM_d5(vrna_fold_compound_t *fc,
88 int *fM_d5);
89
90
91 PRIVATE INLINE void
92 fill_fM_d3(vrna_fold_compound_t *fc,
93 int *fM_d3);
94
95
96 PRIVATE int
97 backtrack(vrna_fold_compound_t *fc,
98 vrna_bp_stack_t *bp_stack,
99 sect bt_stack[],
100 int s);
101
102
103 PRIVATE INLINE int
104 decompose_pair(vrna_fold_compound_t *fc,
105 int i,
106 int j,
107 struct aux_arrays *aux);
108
109
110 PRIVATE INLINE struct aux_arrays *
111 get_aux_arrays(unsigned int length);
112
113
114 PRIVATE INLINE void
115 rotate_aux_arrays(struct aux_arrays *aux,
116 unsigned int length);
117
118
119 PRIVATE INLINE void
120 free_aux_arrays(struct aux_arrays *aux);
121
122
123 /*
124 #################################
125 # BEGIN OF FUNCTION DEFINITIONS #
126 #################################
127 */
128 PUBLIC float
vrna_mfe(vrna_fold_compound_t * fc,char * structure)129 vrna_mfe(vrna_fold_compound_t *fc,
130 char *structure)
131 {
132 char *ss;
133 int length, energy, s;
134 float mfe;
135 sect bt_stack[MAXSECTORS]; /* stack of partial structures for backtracking */
136 vrna_bp_stack_t *bp;
137
138 s = 0;
139 mfe = (float)(INF / 100.);
140
141 if (fc) {
142 length = (int)fc->length;
143
144 if (!vrna_fold_compound_prepare(fc, VRNA_OPTION_MFE)) {
145 vrna_message_warning("vrna_mfe@mfe.c: Failed to prepare vrna_fold_compound");
146 return mfe;
147 }
148
149 /* call user-defined recursion status callback function */
150 if (fc->stat_cb)
151 fc->stat_cb(VRNA_STATUS_MFE_PRE, fc->auxdata);
152
153 /* call user-defined grammar pre-condition callback function */
154 if ((fc->aux_grammar) && (fc->aux_grammar->cb_proc))
155 fc->aux_grammar->cb_proc(fc, VRNA_STATUS_MFE_PRE, fc->aux_grammar->data);
156
157 energy = fill_arrays(fc);
158
159 if (fc->params->model_details.circ)
160 energy = postprocess_circular(fc, bt_stack, &s);
161
162 if (structure && fc->params->model_details.backtrack) {
163 /* add a guess of how many G's may be involved in a G quadruplex */
164 bp = (vrna_bp_stack_t *)vrna_alloc(sizeof(vrna_bp_stack_t) * (4 * (1 + length / 2)));
165
166 if (backtrack(fc, bp, bt_stack, s) != 0) {
167 ss = vrna_db_from_bp_stack(bp, length);
168 strncpy(structure, ss, length + 1);
169 free(ss);
170 } else {
171 memset(structure, '\0', sizeof(char) * (length + 1));
172 }
173
174 free(bp);
175 }
176
177 /* call user-defined recursion status callback function */
178 if (fc->stat_cb)
179 fc->stat_cb(VRNA_STATUS_MFE_POST, fc->auxdata);
180
181 /* call user-defined grammar post-condition callback function */
182 if ((fc->aux_grammar) && (fc->aux_grammar->cb_proc))
183 fc->aux_grammar->cb_proc(fc, VRNA_STATUS_MFE_POST, fc->aux_grammar->data);
184
185 switch (fc->params->model_details.backtrack_type) {
186 case 'C':
187 mfe = (float)fc->matrices->c[fc->jindx[length] + 1] / 100.;
188 break;
189
190 case 'M':
191 mfe = (float)fc->matrices->fML[fc->jindx[length] + 1] / 100.;
192 break;
193
194 default:
195 if (fc->type == VRNA_FC_TYPE_COMPARATIVE)
196 mfe = (float)energy / (100. * (float)fc->n_seq);
197 else
198 mfe = (float)energy / 100.;
199
200 break;
201 }
202 }
203
204 return mfe;
205 }
206
207
208 PUBLIC int
vrna_backtrack_from_intervals(vrna_fold_compound_t * fc,vrna_bp_stack_t * bp_stack,sect bt_stack[],int s)209 vrna_backtrack_from_intervals(vrna_fold_compound_t *fc,
210 vrna_bp_stack_t *bp_stack,
211 sect bt_stack[],
212 int s)
213 {
214 if (fc)
215 return backtrack(fc, bp_stack, bt_stack, s);
216
217 return 0;
218 }
219
220
221 PUBLIC float
vrna_backtrack5(vrna_fold_compound_t * fc,unsigned int length,char * structure)222 vrna_backtrack5(vrna_fold_compound_t *fc,
223 unsigned int length,
224 char *structure)
225 {
226 char *ss;
227 int s;
228 float mfe;
229 sect bt_stack[MAXSECTORS]; /* stack of partial structures for backtracking */
230 vrna_bp_stack_t *bp;
231
232 s = 0;
233 mfe = (float)(INF / 100.);
234
235 if ((fc) && (structure) && (fc->matrices) && (fc->matrices->f5) &&
236 (!fc->params->model_details.circ)) {
237 memset(structure, '\0', sizeof(char) * (length + 1));
238
239 if (length > fc->length)
240 return mfe;
241
242 /* add a guess of how many G's may be involved in a G quadruplex */
243 bp = (vrna_bp_stack_t *)vrna_alloc(sizeof(vrna_bp_stack_t) * (4 * (1 + length / 2)));
244
245 bt_stack[++s].i = 1;
246 bt_stack[s].j = length;
247 bt_stack[s].ml = 0;
248
249
250 if (backtrack(fc, bp, bt_stack, s) != 0) {
251 ss = vrna_db_from_bp_stack(bp, length);
252 strncpy(structure, ss, length + 1);
253 free(ss);
254
255 if (fc->type == VRNA_FC_TYPE_COMPARATIVE)
256 mfe = (float)fc->matrices->f5[length] / (100. * (float)fc->n_seq);
257 else
258 mfe = (float)fc->matrices->f5[length] / 100.;
259 }
260
261 free(bp);
262 }
263
264 return mfe;
265 }
266
267
268 /*
269 #####################################
270 # BEGIN OF STATIC HELPER FUNCTIONS #
271 #####################################
272 */
273
274 /* fill DP matrices */
275 PRIVATE int
fill_arrays(vrna_fold_compound_t * fc)276 fill_arrays(vrna_fold_compound_t *fc)
277 {
278 int i, j, ij, length, turn, uniq_ML, *indx, *f5, *c, *fML, *fM1;
279 vrna_param_t *P;
280 vrna_mx_mfe_t *matrices;
281 vrna_ud_t *domains_up;
282 struct aux_arrays *helper_arrays;
283
284 length = (int)fc->length;
285 indx = fc->jindx;
286 P = fc->params;
287 uniq_ML = P->model_details.uniq_ML;
288 turn = P->model_details.min_loop_size;
289 matrices = fc->matrices;
290 f5 = matrices->f5;
291 c = matrices->c;
292 fML = matrices->fML;
293 fM1 = matrices->fM1;
294 domains_up = fc->domains_up;
295
296 /* allocate memory for all helper arrays */
297 helper_arrays = get_aux_arrays(length);
298
299 if ((turn < 0) || (turn > length))
300 turn = length; /* does this make any sense? */
301
302 /* pre-processing ligand binding production rule(s) */
303 if (domains_up && domains_up->prod_cb)
304 domains_up->prod_cb(fc, domains_up->data);
305
306 /* prefill matrices with init contributions */
307 for (j = 1; j <= length; j++)
308 for (i = (j > turn ? (j - turn) : 1); i <= j; i++) {
309 c[indx[j] + i] = fML[indx[j] + i] = INF;
310 if (uniq_ML)
311 fM1[indx[j] + i] = INF;
312 }
313
314 /* start recursion */
315 if (length <= turn) {
316 /* clean up memory */
317 free_aux_arrays(helper_arrays);
318
319 /* return free energy of unfolded chain */
320 return 0;
321 }
322
323 for (i = length - turn - 1; i >= 1; i--) {
324 for (j = i + turn + 1; j <= length; j++) {
325 ij = indx[j] + i;
326
327 /* decompose subsegment [i, j] with pair (i, j) */
328 c[ij] = decompose_pair(fc, i, j, helper_arrays);
329
330 /* decompose subsegment [i, j] that is multibranch loop part with at least one branch */
331 fML[ij] = vrna_E_ml_stems_fast(fc, i, j, helper_arrays->Fmi, helper_arrays->DMLi);
332
333 /* decompose subsegment [i, j] that is multibranch loop part with exactly one branch */
334 if (uniq_ML)
335 fM1[ij] = E_ml_rightmost_stem(i, j, fc);
336
337 if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux))
338 fc->aux_grammar->cb_aux(fc, i, j, fc->aux_grammar->data);
339 } /* end of j-loop */
340
341 rotate_aux_arrays(helper_arrays, length);
342 } /*
343 * end of i-loop
344 * calculate energies of 5' fragments
345 */
346 (void)vrna_E_ext_loop_5(fc);
347
348 /* clean up memory */
349 free_aux_arrays(helper_arrays);
350
351 return f5[length];
352 }
353
354
355 /* post-processing step for circular RNAs */
356 PRIVATE int
postprocess_circular(vrna_fold_compound_t * fc,sect bt_stack[],int * bt)357 postprocess_circular(vrna_fold_compound_t *fc,
358 sect bt_stack[],
359 int *bt)
360 {
361 /*
362 * auxiliarry arrays:
363 * fM2 = multiloop region with exactly two stems, extending to 3' end
364 * for stupid dangles=1 case we also need:
365 * fM_d3 = multiloop region with >= 2 stems, starting at pos 2
366 * (a pair (k,n) will form 3' dangle with pos 1)
367 * fM_d5 = multiloop region with >= 2 stems, extending to pos n-1
368 * (a pair (1,k) will form a 5' dangle with pos n)
369 */
370 unsigned char *hard_constraints, eval;
371 char *ptype;
372 short *S1, **SS, **S5, **S3;
373 unsigned int **a2s;
374 int Hi, Hj, Ii, Ij, Ip, Iq, ip, iq, Mi, *fM_d3, *fM_d5, Md3i,
375 Md5i, FcMd3, FcMd5, FcH, FcI, FcM, Fc, *fM2, i, j, ij, u,
376 length, new_c, fm, type, *my_c, *my_fML, *indx, FcO, tmp,
377 dangle_model, turn, s, n_seq;
378 vrna_param_t *P;
379 vrna_md_t *md;
380 vrna_hc_t *hc;
381 vrna_sc_t *sc, **scs;
382
383 length = fc->length;
384 n_seq = (fc->type == VRNA_FC_TYPE_SINGLE) ? 1 : fc->n_seq;
385 P = fc->params;
386 md = &(P->model_details);
387 ptype = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->ptype : NULL;
388 indx = fc->jindx;
389 S1 = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->sequence_encoding : NULL;
390 SS = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->S;
391 S5 = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->S5;
392 S3 = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->S3;
393 a2s = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->a2s;
394 hc = fc->hc;
395 sc = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->sc : NULL;
396 scs = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->scs;
397 dangle_model = md->dangles;
398 turn = md->min_loop_size;
399 hard_constraints = hc->mx;
400 my_c = fc->matrices->c;
401 my_fML = fc->matrices->fML;
402 fM2 = fc->matrices->fM2;
403
404 Fc = FcO = FcH = FcI = FcM = FcMd3 = FcMd5 = INF;
405 Mi = Md5i = Md3i = Iq = Ip = Ij = Ii = Hj = Hi = 0;
406
407 /* unfolded state */
408 eval = (hc->up_ext[1] >= length) ? 1 : 0;
409 if (hc->f)
410 eval = (hc->f(1, length, 1, length, VRNA_DECOMP_EXT_UP, hc->data)) ? eval : 0;
411
412 if (eval) {
413 Fc = 0; /* base line for unfolded state */
414
415 switch (fc->type) {
416 case VRNA_FC_TYPE_SINGLE:
417 if (sc) {
418 if (sc->energy_up)
419 Fc += sc->energy_up[1][length];
420
421 if (sc->f)
422 Fc += sc->f(1, length, 1, length, VRNA_DECOMP_EXT_UP, sc->data);
423 }
424
425 break;
426
427 case VRNA_FC_TYPE_COMPARATIVE:
428 if (scs) {
429 for (s = 0; s < fc->n_seq; s++)
430 if (scs[s])
431 if (scs[s]->energy_up)
432 Fc += scs[s]->energy_up[1][a2s[s][length]];
433 }
434
435 break;
436 }
437 FcO = Fc;
438 } else {
439 Fc = INF;
440 }
441
442 for (i = 1; i < length; i++)
443 for (j = i + turn + 1; j <= length; j++) {
444 u = length - j + i - 1;
445 if (u < turn)
446 continue;
447
448 ij = indx[j] + i;
449
450 if (!hard_constraints[length * i + j])
451 continue;
452
453 /* exterior hairpin case */
454 new_c = vrna_E_hp_loop(fc, j, i);
455 if (new_c != INF)
456 new_c += my_c[ij];
457
458 if (new_c < FcH) {
459 FcH = new_c;
460 Hi = i;
461 Hj = j;
462 }
463
464 /* exterior interior loop case */
465 ip = iq = 0;
466 new_c = vrna_E_ext_int_loop(fc, i, j, &ip, &iq);
467 if (new_c != INF)
468 new_c += my_c[ij];
469
470 if (ip != 0) {
471 if (new_c < FcI) {
472 FcI = new_c;
473 Ii = i;
474 Ij = j;
475 Ip = ip;
476 Iq = iq;
477 }
478 }
479 } /* end of i,j loop */
480 Fc = MIN2(Fc, FcH);
481 Fc = MIN2(Fc, FcI);
482
483 /*
484 * compute the fM2 array (multi loops with exactly 2 helices)
485 * to get a unique ML decomposition, just use fM1 instead of fML
486 * below. However, that will not work with dangle_model==1
487 */
488 int *fml_tmp = my_fML;
489
490 /* some pre-processing to reduce redundant code */
491 if ((hc->f) ||
492 ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc)) ||
493 ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)))
494 fml_tmp = vrna_alloc(sizeof(int) * (length + 2));
495 else
496 fml_tmp += indx[length];
497
498 for (i = 1; i < length - turn; i++) {
499 fM2[i] = INF;
500 /* some preparations in case we have to comply/add certain constraints */
501 if ((fml_tmp - indx[length]) != my_fML) {
502 /* copy original data */
503 for (u = i + turn; u < length - turn; u++)
504 fml_tmp[u + 1] = my_fML[indx[length] + u + 1];
505
506 /* apply hard constraints */
507 if (hc->f) {
508 for (u = i + turn; u < length - turn; u++)
509 if (!hc->f(i, length, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
510 fml_tmp[u + 1] = INF;
511 }
512
513 switch (fc->type) {
514 case VRNA_FC_TYPE_SINGLE:
515 if ((sc) && (sc->f)) {
516 for (u = i + turn; u < length - turn; u++) {
517 fm = sc->f(i, length, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
518 if ((fm != INF) && (fml_tmp[u + 1] != INF))
519 fml_tmp[u + 1] += fm;
520 }
521 }
522
523 break;
524
525 case VRNA_FC_TYPE_COMPARATIVE:
526 if (scs) {
527 for (u = i + turn; u < length - turn; u++) {
528 fm = 0;
529 for (s = 0; s < n_seq; s++)
530 if ((scs[s]) && (scs[s]->f))
531 fm += scs[s]->f(i, length, u, u + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
532
533 if ((fm != INF) && (fml_tmp[u + 1] != INF))
534 fml_tmp[u + 1] += fm;
535 }
536 }
537
538 break;
539 }
540 }
541
542 /* actual decomposition */
543 for (u = i + turn; u < length - turn; u++) {
544 fm = my_fML[indx[u] + i];
545 if ((fm != INF) && (fml_tmp[u + 1] != INF)) {
546 fm += fml_tmp[u + 1];
547 fM2[i] = MIN2(fM2[i], fm);
548 }
549 }
550 }
551
552 if ((fml_tmp - indx[length]) != my_fML)
553 free(fml_tmp);
554
555 /*
556 * Now, process all exterior multibranch loop configurations
557 */
558 int *fm2_tmp = fM2;
559
560 if (hc->f) {
561 fm2_tmp = vrna_alloc(sizeof(int) * (length + 1));
562
563 /* copy data */
564 for (i = turn + 1; i < length - 2 * turn; i++)
565 fm2_tmp[i + 1] = fM2[i + 1];
566
567 /* apply hard constraints */
568 for (i = turn + 1; i < length - 2 * turn; i++)
569 if (!hc->f(1, length, i, i + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
570 fm2_tmp[i + 1] = INF;
571 }
572
573 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->f)) {
574 if (fm2_tmp == fM2) {
575 fm2_tmp = vrna_alloc(sizeof(int) * (length + 1));
576
577 /* copy data */
578 for (i = turn + 1; i < length - 2 * turn; i++)
579 fm2_tmp[i + 1] = fM2[i + 1];
580 }
581
582 /* add soft constraints */
583 for (i = turn + 1; i < length - 2 * turn; i++) {
584 fm = sc->f(1, length, i, i + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
585 if ((fm != INF) && (fm2_tmp[i + 1] != INF))
586 fm2_tmp[i + 1] += fm;
587 }
588 }
589
590 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
591 if (fm2_tmp == fM2) {
592 fm2_tmp = vrna_alloc(sizeof(int) * (length + 1));
593
594 /* copy data */
595 for (i = turn + 1; i < length - 2 * turn; i++)
596 fm2_tmp[i + 1] = fM2[i + 1];
597 }
598
599 /* add soft constraints */
600 for (i = turn + 1; i < length - 2 * turn; i++) {
601 fm = 0;
602 for (s = 0; s < n_seq; s++)
603 if ((scs[s]) && (scs[s]->f))
604 fm += scs[s]->f(1, length, i, i + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
605
606 if ((fm != INF) && (fm2_tmp[i + 1] != INF))
607 fm2_tmp[i + 1] += fm;
608 }
609 }
610
611 /* actual decomposition */
612 for (i = turn + 1; i < length - 2 * turn; i++) {
613 if ((my_fML[indx[i] + 1] != INF) && (fm2_tmp[i + 1] != INF)) {
614 fm = my_fML[indx[i] + 1] +
615 fm2_tmp[i + 1];
616
617 if (fm < FcM) {
618 FcM = fm;
619 Mi = i;
620 }
621 }
622 }
623
624 if (fm2_tmp != fM2)
625 free(fm2_tmp);
626
627 if (FcM != INF) {
628 switch (fc->type) {
629 case VRNA_FC_TYPE_SINGLE:
630 FcM += P->MLclosing;
631 break;
632
633 case VRNA_FC_TYPE_COMPARATIVE:
634 FcM += n_seq * P->MLclosing;
635 break;
636 }
637 }
638
639 Fc = MIN2(Fc, FcM);
640
641 /* add multibranch loop configurations for odd dangle models */
642 if ((dangle_model == 1) || (dangle_model == 3)) {
643 fM_d5 = (int *)vrna_alloc(sizeof(int) * (length + 2));
644 fM_d3 = (int *)vrna_alloc(sizeof(int) * (length + 2));
645
646 for (i = turn + 1; i < length - turn; i++)
647 fM_d5[i] = INF;
648
649 for (i = turn + 1; i < length - turn; i++)
650 fM_d3[i] = INF;
651
652 if (hc->up_ml[1]) {
653 fill_fM_d3(fc, fM_d3);
654
655 /*
656 **********************************************************************
657 * 1. test ML loops with closing pair (i + 1, length), 3' dangle pos 1
658 **********************************************************************
659 */
660 int *c_tmp, *c_tmp2;
661
662 c_tmp = vrna_alloc(sizeof(int) * (length + 2));
663 c_tmp2 = my_c + indx[length];
664
665 /* copy contributions for enclosed structure part */
666 for (i = 2 * turn + 1; i < length - turn; i++)
667 c_tmp[i + 1] = c_tmp2[i + 1];
668
669 /* obey hard constraints */
670 if (hc->f) {
671 for (i = 2 * turn + 1; i < length - turn; i++)
672 if (!hc->f(i + 1, length, 2, i, VRNA_DECOMP_PAIR_ML_EXT, hc->data))
673 c_tmp[i + 1] = INF;
674 }
675
676 /* add soft constraints */
677 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->f)) {
678 for (i = 2 * turn + 1; i < length - turn; i++) {
679 tmp = sc->f(i + 1, length, 2, i,
680 VRNA_DECOMP_PAIR_ML_EXT,
681 sc->data);
682 if (tmp != INF) {
683 if (c_tmp[i + 1] != INF)
684 c_tmp[i + 1] += tmp;
685 } else {
686 c_tmp[i + 1] = INF;
687 }
688 }
689 }
690
691 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
692 for (i = 2 * turn + 1; i < length - turn; i++) {
693 tmp = 0;
694 for (s = 0; s < n_seq; s++)
695 if ((scs[s]) && (scs[s]->f))
696 tmp += scs[s]->f(i + 1, length, 2, i,
697 VRNA_DECOMP_PAIR_ML_EXT,
698 scs[s]->data);
699
700 if (tmp != INF) {
701 if (c_tmp[i + 1] != INF)
702 c_tmp[i + 1] += tmp;
703 } else {
704 c_tmp[i + 1] = INF;
705 }
706 }
707 }
708
709 /* add contributions for enclosing pair */
710 for (i = 2 * turn + 1; i < length - turn; i++) {
711 if (c_tmp[i + 1] != INF) {
712 /* obey internal hard constraints */
713 if (hard_constraints[length * length + i + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
714 tmp = 0;
715 switch (fc->type) {
716 case VRNA_FC_TYPE_SINGLE:
717 type = vrna_get_ptype(indx[length] + i + 1, ptype);
718 tmp = E_MLstem(type, -1, S1[1], P) +
719 P->MLclosing;
720 break;
721
722 case VRNA_FC_TYPE_COMPARATIVE:
723 tmp = P->MLclosing * n_seq;
724 for (s = 0; s < n_seq; s++) {
725 type = vrna_get_ptype_md(SS[s][i + 1], SS[s][length], md);
726 tmp += E_MLstem(type, -1, S3[s][length], P);
727 }
728 break;
729 }
730
731 c_tmp[i + 1] += tmp;
732 } else {
733 c_tmp[i + 1] = INF;
734 }
735 }
736 }
737
738 /* actual decomposition */
739 for (i = 2 * turn + 1; i < length - turn; i++) {
740 fm = fM_d3[i];
741 if ((fm != INF) && (c_tmp[i + 1] != INF)) {
742 fm += c_tmp[i + 1];
743 if (fm < FcMd3) {
744 FcMd3 = fm;
745 Md3i = i;
746 }
747 }
748 }
749
750 /*
751 **********************************************************************
752 * 2. test ML loops with closing pair (i + 1, length), 5' dangle
753 * pos i, 3' dangle pos 1
754 **********************************************************************
755 */
756
757 /* copy contributions for enclosed structure part */
758 for (i = 2 * turn + 1; i < length - turn; i++)
759 c_tmp[i + 1] = c_tmp2[i + 1];
760
761
762 /* obey hard constraints */
763 if (hc->f) {
764 for (i = 2 * turn + 1; i < length - turn; i++)
765 if (!hc->f(i + 1, length, 2, i - 1, VRNA_DECOMP_PAIR_ML_EXT, hc->data))
766 c_tmp[i + 1] = INF;
767 }
768
769 /* add soft constraints (static and user-defined) */
770 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc)) {
771 if (sc->energy_up) {
772 for (i = 2 * turn + 1; i < length - turn; i++)
773 if (c_tmp[i + 1] != INF)
774 c_tmp[i + 1] += sc->energy_up[i][1];
775 }
776
777 if (sc->f) {
778 for (i = 2 * turn + 1; i < length - turn; i++) {
779 tmp = sc->f(i + 1, length, 2, i - 1,
780 VRNA_DECOMP_PAIR_ML_EXT,
781 sc->data);
782 if (tmp == INF)
783 c_tmp[i + 1] = INF;
784 else if (c_tmp[i + 1] != INF)
785 c_tmp[i + 1] += tmp;
786 }
787 }
788 }
789
790 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
791 for (i = 2 * turn + 1; i < length - turn; i++) {
792 if (c_tmp[i + 1] != INF) {
793 for (s = 0; s < n_seq; s++) {
794 if (scs[s]) {
795 if (scs[s]->energy_up)
796 c_tmp[i + 1] += scs[s]->energy_up[a2s[s][i]][1];
797
798 if (scs[s]->f) {
799 c_tmp[i + 1] += scs[s]->f(i + 1, length, 2, i - 1,
800 VRNA_DECOMP_PAIR_ML_EXT,
801 scs[s]->data);
802 }
803 }
804 }
805 }
806 }
807 }
808
809 /* add contributions for enclosing pair */
810 for (i = 2 * turn + 1; i < length - turn; i++) {
811 if (c_tmp[i + 1] != INF) {
812 /* obey internal hard constraints */
813 if ((hard_constraints[length * length + i + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) &&
814 (hc->up_ml[i])) {
815 tmp = 0;
816 switch (fc->type) {
817 case VRNA_FC_TYPE_SINGLE:
818 type = vrna_get_ptype(indx[length] + i + 1, ptype);
819 tmp = E_MLstem(type, S1[i], S1[1], P) +
820 P->MLclosing;
821 break;
822
823 case VRNA_FC_TYPE_COMPARATIVE:
824 tmp = P->MLclosing * n_seq;
825 for (s = 0; s < n_seq; s++) {
826 type = vrna_get_ptype_md(SS[s][i + 1], SS[s][length], md);
827 tmp += E_MLstem(type, S5[s][i + 1], S3[s][length], P);
828 }
829 break;
830 }
831
832 c_tmp[i + 1] += tmp;
833 } else {
834 c_tmp[i + 1] = INF;
835 }
836 }
837 }
838
839 /* actual decomposition */
840 for (i = 2 * turn + 1; i < length - turn; i++) {
841 fm = fM_d3[i - 1];
842 if ((fm != INF) && (c_tmp[i + 1] != INF)) {
843 fm += c_tmp[i + 1];
844 if (fm < FcMd3) {
845 FcMd3 = fm;
846 Md3i = -i;
847 }
848 }
849 }
850
851 free(c_tmp);
852 }
853
854 if (hc->up_ml[length]) {
855 fill_fM_d5(fc, fM_d5);
856
857 /*
858 **********************************************************************
859 * 1. test ML loops with closing pair (1, i), 5' dangle pos n
860 **********************************************************************
861 */
862 int *fmd5_tmp = vrna_alloc(sizeof(int) * (length + 2));
863
864 /* copy contributions */
865 for (i = turn + 1; i < length - turn; i++)
866 fmd5_tmp[i + 1] = fM_d5[i + 1];
867
868 /* obey hard constraints */
869 if (hc->f) {
870 for (i = turn + 1; i < length - turn; i++)
871 if (!hc->f(1, i, i + 1, length - 1, VRNA_DECOMP_PAIR_ML_EXT, hc->data))
872 fmd5_tmp[i + 1] = INF;
873 }
874
875 /* add soft constraints */
876 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->f)) {
877 for (i = turn + 1; i < length - turn; i++) {
878 fm = sc->f(1, i, i + 1, length - 1,
879 VRNA_DECOMP_PAIR_ML_EXT,
880 sc->data);
881 if ((fm != INF) && (fmd5_tmp[i + 1] != INF))
882 fmd5_tmp += fm;
883 }
884 }
885
886 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
887 for (i = turn + 1; i < length - turn; i++) {
888 fm = 0;
889 for (s = 0; s < n_seq; s++)
890 if ((scs[s]) && (scs[s]->f))
891 fm += scs[s]->f(1, i, i + 1, length - 1,
892 VRNA_DECOMP_PAIR_ML_EXT,
893 scs[s]->data);
894
895 if ((fm != INF) && (fmd5_tmp[i + 1] != INF))
896 fmd5_tmp += fm;
897 }
898 }
899
900 /* add contributions for enclosing pair */
901 for (i = turn + 1; i < length - turn; i++) {
902 if (fmd5_tmp[i + 1] != INF) {
903 if (hard_constraints[length + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
904 tmp = 0;
905 switch (fc->type) {
906 case VRNA_FC_TYPE_SINGLE:
907 type = vrna_get_ptype(indx[i] + 1, ptype);
908 tmp = E_MLstem(type, S1[length], -1, P) +
909 P->MLclosing;
910 break;
911
912 case VRNA_FC_TYPE_COMPARATIVE:
913 tmp = P->MLclosing * n_seq;
914 for (s = 0; s < n_seq; s++) {
915 type = vrna_get_ptype_md(SS[s][1], SS[s][i], md);
916 tmp += E_MLstem(type, S5[s][1], -1, P);
917 }
918 break;
919 }
920
921 fmd5_tmp[i + 1] += tmp;
922 } else {
923 fmd5_tmp[i + 1] = INF;
924 }
925 }
926 }
927 /* actual decomposition */
928 for (i = turn + 1; i < length - turn; i++) {
929 fm = fmd5_tmp[i + 1];
930 if ((fm != INF) && (my_c[indx[i] + 1] != INF)) {
931 fm += my_c[indx[i] + 1];
932 if (fm < FcMd5) {
933 FcMd5 = fm;
934 Md5i = i;
935 }
936 }
937 }
938
939 /*
940 **********************************************************************
941 * 2. test ML loops with closing pair (1, i), 5' dangle pos n, 3' dangle
942 * pos i + 1
943 **********************************************************************
944 */
945
946 /* copy contributions */
947 for (i = turn + 1; i < length - turn; i++)
948 fmd5_tmp[i + 2] = fM_d5[i + 2];
949
950 /* obey hard constraints */
951 if (hc->f) {
952 for (i = turn + 1; i < length - turn; i++)
953 if (!hc->f(1, i, i + 2, length - 1, VRNA_DECOMP_PAIR_ML_EXT, hc->data))
954 fmd5_tmp[i + 2] = INF;
955 }
956
957 /* add soft constraints (static and user-defined) */
958 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc)) {
959 if (sc->energy_up) {
960 for (i = turn + 1; i < length - turn; i++)
961 if (fmd5_tmp[i + 2] != INF)
962 fmd5_tmp[i + 2] += sc->energy_up[i + 1][1];
963 }
964
965 if (sc->f) {
966 for (i = turn + 1; i < length - turn; i++) {
967 tmp = sc->f(1, i, i + 2, length - 1,
968 VRNA_DECOMP_PAIR_ML_EXT,
969 sc->data);
970 if (tmp == INF)
971 fmd5_tmp[i + 2] = INF;
972 else if (fmd5_tmp[i + 2] != INF)
973 fmd5_tmp[i + 2] += tmp;
974 }
975 }
976 }
977
978 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
979 for (i = turn + 1; i < length - turn; i++) {
980 if (fmd5_tmp[i + 2] != INF) {
981 for (s = 0; s < n_seq; s++) {
982 if (scs[s]) {
983 if (scs[s]->energy_up)
984 fmd5_tmp[i + 2] += scs[s]->energy_up[a2s[s][i + 1]][1];
985
986 if (scs[s]->f) {
987 fmd5_tmp[i + 2] += scs[s]->f(1, i, i + 2, length - 1,
988 VRNA_DECOMP_PAIR_ML_EXT,
989 scs[s]->data);
990 }
991 }
992 }
993 }
994 }
995 }
996
997 /* add contributions for enclosing pair */
998 for (i = turn + 1; i < length - turn; i++) {
999 if (fmd5_tmp[i + 2] != INF) {
1000 /* obey internal hard constraints */
1001 if ((hard_constraints[length + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) &&
1002 (hc->up_ml[i + 1])) {
1003 tmp = 0;
1004 switch (fc->type) {
1005 case VRNA_FC_TYPE_SINGLE:
1006 type = vrna_get_ptype(indx[i] + 1, ptype);
1007 tmp = E_MLstem(type, S1[length], S1[i + 1], P) +
1008 P->MLclosing;
1009 break;
1010
1011 case VRNA_FC_TYPE_COMPARATIVE:
1012 tmp = P->MLclosing * n_seq;
1013 for (s = 0; s < n_seq; s++) {
1014 type = vrna_get_ptype_md(SS[s][1], SS[s][i], md);
1015 tmp += E_MLstem(type, S5[s][1], S3[s][i], P);
1016 }
1017 break;
1018 }
1019
1020 fmd5_tmp[i + 2] += tmp;
1021 } else {
1022 fmd5_tmp[i + 2] = INF;
1023 }
1024 }
1025 }
1026
1027 /* actual decomposition */
1028 for (i = turn + 1; i < length - turn; i++) {
1029 fm = fmd5_tmp[i + 2];
1030 if ((fm != INF) && (my_c[indx[i] + 1] != INF)) {
1031 fm += my_c[indx[i] + 1];
1032
1033 if (fm < FcMd5) {
1034 FcMd5 = fm;
1035 Md5i = -i;
1036 }
1037 }
1038 }
1039
1040 free(fmd5_tmp);
1041 }
1042
1043 if (FcMd5 < MIN2(Fc, FcMd3)) {
1044 int real_i, sc_en = 0;
1045
1046 /* looks like we have to do this ... */
1047 bt_stack[++(*bt)].i = 1;
1048 bt_stack[(*bt)].j = (Md5i > 0) ? Md5i : -Md5i;
1049 bt_stack[(*bt)].ml = 2;
1050 i = (Md5i > 0) ? Md5i + 1 : -Md5i + 2; /* let's backtrack fm_d5[Md5i+1] */
1051 real_i = (Md5i > 0) ? i : i - 1;
1052
1053 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->energy_up)) {
1054 sc_en += sc->energy_up[length][1];
1055 } else if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
1056 for (s = 0; s < n_seq; s++)
1057 if ((scs[s]) && (scs[s]->energy_up))
1058 sc_en += scs[s]->energy_up[a2s[s][length]][1];
1059 }
1060
1061 for (u = i + turn; u < length - turn; u++) {
1062 fm = my_fML[indx[u] + i] +
1063 my_fML[indx[length - 1] + u + 1] +
1064 sc_en;
1065
1066 switch (fc->type) {
1067 case VRNA_FC_TYPE_SINGLE:
1068 if (sc) {
1069 if (sc->energy_up)
1070 fm += sc->energy_up[real_i][i - real_i];
1071
1072 if (sc->f) {
1073 fm += sc->f(real_i, length, i, length - 1,
1074 VRNA_DECOMP_ML_ML,
1075 sc->data) +
1076 sc->f(i, length - 1, u, u + 1,
1077 VRNA_DECOMP_ML_ML_ML,
1078 sc->data);
1079 }
1080 }
1081
1082 break;
1083
1084 case VRNA_FC_TYPE_COMPARATIVE:
1085 if (scs) {
1086 for (s = 0; s < n_seq; s++) {
1087 if (scs[s]) {
1088 if (scs[s]->energy_up)
1089 fm += scs[s]->energy_up[a2s[s][real_i]][i - real_i];
1090
1091 if (scs[s]->f) {
1092 fm += scs[s]->f(real_i, length, i, length - 1,
1093 VRNA_DECOMP_ML_ML,
1094 scs[s]->data) +
1095 scs[s]->f(i, length - 1, u, u + 1,
1096 VRNA_DECOMP_ML_ML_ML,
1097 scs[s]->data);
1098 }
1099 }
1100 }
1101 }
1102
1103 break;
1104 }
1105
1106 if (fM_d5[i] == fm) {
1107 bt_stack[++(*bt)].i = i;
1108 bt_stack[(*bt)].j = u;
1109 bt_stack[(*bt)].ml = 1;
1110 bt_stack[++(*bt)].i = u + 1;
1111 bt_stack[(*bt)].j = length - 1;
1112 bt_stack[(*bt)].ml = 1;
1113 break;
1114 }
1115 }
1116 Fc = FcMd5;
1117 } else if (FcMd3 < Fc) {
1118 int real_i, sc_en = 0;
1119 /* here we go again... */
1120 bt_stack[++(*bt)].i = (Md3i > 0) ? Md3i + 1 : -Md3i + 1;
1121 bt_stack[(*bt)].j = length;
1122 bt_stack[(*bt)].ml = 2;
1123 i = (Md3i > 0) ? Md3i : -Md3i - 1; /* let's backtrack fm_d3[Md3i] */
1124 real_i = (Md3i > 0) ? i : i + 1;
1125
1126 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->energy_up)) {
1127 sc_en += sc->energy_up[1][1];
1128 } else if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
1129 for (s = 0; s < n_seq; s++)
1130 if ((scs[s]) && (scs[s]->energy_up))
1131 sc_en += scs[s]->energy_up[a2s[s][1]][1];
1132 }
1133
1134 for (u = 2 + turn; u < i - turn; u++) {
1135 fm = my_fML[indx[u] + 2] +
1136 my_fML[indx[i] + u + 1] +
1137 sc_en;
1138
1139 switch (fc->type) {
1140 case VRNA_FC_TYPE_SINGLE:
1141 if (sc) {
1142 if (sc->energy_up)
1143 fm += sc->energy_up[real_i][real_i - i];
1144
1145 if (sc->f) {
1146 fm += sc->f(1, real_i, 2, i,
1147 VRNA_DECOMP_ML_ML,
1148 sc->data) +
1149 sc->f(2, i, u, u + 1,
1150 VRNA_DECOMP_ML_ML_ML,
1151 sc->data);
1152 }
1153 }
1154
1155 break;
1156
1157 case VRNA_FC_TYPE_COMPARATIVE:
1158 if (scs) {
1159 for (s = 0; s < n_seq; s++) {
1160 if (scs[s]) {
1161 if (scs[s]->energy_up)
1162 fm += scs[s]->energy_up[a2s[s][real_i]][real_i - i];
1163
1164 if (scs[s]->f) {
1165 fm += scs[s]->f(1, real_i, 2, i,
1166 VRNA_DECOMP_ML_ML,
1167 scs[s]->data) +
1168 scs[s]->f(2, i, u, u + 1,
1169 VRNA_DECOMP_ML_ML_ML,
1170 scs[s]->data);
1171 }
1172 }
1173 }
1174 }
1175
1176 break;
1177 }
1178
1179 if (fM_d3[i] == fm) {
1180 bt_stack[++(*bt)].i = 2;
1181 bt_stack[(*bt)].j = u;
1182 bt_stack[(*bt)].ml = 1;
1183 bt_stack[++(*bt)].i = u + 1;
1184 bt_stack[(*bt)].j = i;
1185 bt_stack[(*bt)].ml = 1;
1186 break;
1187 }
1188 }
1189 Fc = FcMd3;
1190 }
1191
1192 free(fM_d3);
1193 free(fM_d5);
1194 }
1195
1196 if (Fc < INF) {
1197 if (FcH == Fc) {
1198 bt_stack[++(*bt)].i = Hi;
1199 bt_stack[(*bt)].j = Hj;
1200 bt_stack[(*bt)].ml = 2;
1201 } else if (FcI == Fc) {
1202 bt_stack[++(*bt)].i = Ii;
1203 bt_stack[(*bt)].j = Ij;
1204 bt_stack[(*bt)].ml = 2;
1205 bt_stack[++(*bt)].i = Ip;
1206 bt_stack[(*bt)].j = Iq;
1207 bt_stack[(*bt)].ml = 2;
1208 } else if (FcM == Fc) {
1209 /* grumpf we found a Multiloop */
1210 int eee;
1211 /* backtrack in fM2 */
1212 fm = fM2[Mi + 1];
1213 for (u = Mi + turn + 1; u < length - turn; u++) {
1214 eee = my_fML[indx[u] + Mi + 1] +
1215 my_fML[indx[length] + u + 1];
1216
1217 switch (fc->type) {
1218 case VRNA_FC_TYPE_SINGLE:
1219 if (sc) {
1220 if (sc->f)
1221 eee += sc->f(Mi + 1, length, u, u + 1,
1222 VRNA_DECOMP_ML_ML_ML,
1223 sc->data);
1224 }
1225
1226 break;
1227
1228 case VRNA_FC_TYPE_COMPARATIVE:
1229 if (scs) {
1230 for (s = 0; s < n_seq; s++)
1231 if ((scs[s]) && (scs[s]->f))
1232 eee += scs[s]->f(Mi + 1, length, u, u + 1,
1233 VRNA_DECOMP_ML_ML_ML,
1234 scs[s]->data);
1235 }
1236
1237 break;
1238 }
1239
1240 if (fm == eee) {
1241 bt_stack[++(*bt)].i = Mi + 1;
1242 bt_stack[(*bt)].j = u;
1243 bt_stack[(*bt)].ml = 1;
1244 bt_stack[++(*bt)].i = u + 1;
1245 bt_stack[(*bt)].j = length;
1246 bt_stack[(*bt)].ml = 1;
1247 break;
1248 }
1249 }
1250 bt_stack[++(*bt)].i = 1;
1251 bt_stack[(*bt)].j = Mi;
1252 bt_stack[(*bt)].ml = 1;
1253 } else if (Fc == FcO) {
1254 /* unstructured */
1255 bt_stack[++(*bt)].i = 1;
1256 bt_stack[(*bt)].j = 1;
1257 bt_stack[(*bt)].ml = 0;
1258 }
1259 } else {
1260 /* forbidden, i.e. no configuration fulfills constraints */
1261 }
1262
1263 fc->matrices->FcH = FcH;
1264 fc->matrices->FcI = FcI;
1265 fc->matrices->FcM = FcM;
1266 fc->matrices->Fc = Fc;
1267 return Fc;
1268 }
1269
1270
1271 /*
1272 **************************
1273 * Construct fM_d5 array
1274 **************************
1275 */
1276 PRIVATE INLINE void
fill_fM_d5(vrna_fold_compound_t * fc,int * fM_d5)1277 fill_fM_d5(vrna_fold_compound_t *fc,
1278 int *fM_d5)
1279 {
1280 unsigned int s, n_seq, **a2s;
1281 int *fm_tmp, *fm_tmp2, sc_base, fm, tmp, i, u,
1282 length, turn, *my_fML, *indx;
1283 vrna_md_t *md;
1284 vrna_param_t *P;
1285 vrna_hc_t *hc;
1286 vrna_sc_t *sc, **scs;
1287
1288 n_seq = (fc->type == VRNA_FC_TYPE_SINGLE) ? 1 : fc->n_seq;
1289 length = fc->length;
1290 a2s = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->a2s;
1291 P = fc->params;
1292 md = &(P->model_details);
1293 my_fML = fc->matrices->fML;
1294 hc = fc->hc;
1295 sc = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->sc : NULL;
1296 scs = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->scs;
1297 indx = fc->jindx;
1298 turn = md->min_loop_size;
1299
1300 fm_tmp2 = vrna_alloc(sizeof(int) * (length + 2));
1301 sc_base = 0;
1302
1303 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->energy_up))
1304 sc_base += sc->energy_up[length][1];
1305 else if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs))
1306 for (s = 0; s < n_seq; s++)
1307 if ((scs[s]) && (scs[s]->energy_up))
1308 sc_base += scs[s]->energy_up[a2s[s][length]][1];
1309
1310 for (i = turn + 1; i < length - turn; i++) {
1311 fm_tmp = my_fML + indx[length - 1];
1312
1313 if (sc_base != 0) {
1314 fm_tmp = fm_tmp2;
1315
1316 for (u = 2 + turn; u < i - turn; u++)
1317 fm_tmp[u + 1] = my_fML[indx[length - 1] + u + 1] +
1318 sc_base;
1319 }
1320
1321 if (hc->f) {
1322 if (!hc->f(i, length, i, length - 1, VRNA_DECOMP_ML_ML, hc->data))
1323 continue;
1324
1325 if (fm_tmp != fm_tmp2) {
1326 fm_tmp = fm_tmp2;
1327
1328 for (u = 2 + turn; u < i - turn; u++)
1329 fm_tmp[u + 1] = my_fML[indx[length - 1] + u + 1];
1330 }
1331
1332 for (u = 2 + turn; u < i - turn; u++)
1333 if (!hc->f(i, length - 1, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
1334 fm_tmp[u + 1] = INF;
1335 }
1336
1337 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->f)) {
1338 if (fm_tmp != fm_tmp2) {
1339 fm_tmp = fm_tmp2;
1340
1341 for (u = 2 + turn; u < i - turn; u++)
1342 fm_tmp[u + 1] = my_fML[indx[length - 1] + u + 1];
1343 }
1344
1345 fm = sc->f(i, length, i, length - 1,
1346 VRNA_DECOMP_ML_ML,
1347 sc->data);
1348
1349 if (fm != INF) {
1350 for (u = 2 + turn; u < i - turn; u++) {
1351 if (fm_tmp[u + 1] != INF) {
1352 tmp = sc->f(i, length - 1, u, u + 1,
1353 VRNA_DECOMP_ML_ML_ML,
1354 sc->data);
1355 if (tmp != INF)
1356 tmp += fm;
1357
1358 fm_tmp[u + 1] += tmp;
1359 }
1360 }
1361 } else {
1362 for (u = 2 + turn; u < i - turn; u++)
1363 fm_tmp[u + 1] = INF;
1364 }
1365 }
1366
1367 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
1368 if (fm_tmp != fm_tmp2) {
1369 fm_tmp = fm_tmp2;
1370
1371 for (u = 2 + turn; u < i - turn; u++)
1372 fm_tmp[u + 1] = my_fML[indx[length - 1] + u + 1];
1373 }
1374
1375 fm = 0;
1376 for (s = 0; s < n_seq; s++)
1377 if ((scs[s]) && (scs[s]->f))
1378 fm += scs[s]->f(i, length, i, length - 1,
1379 VRNA_DECOMP_ML_ML,
1380 scs[s]->data);
1381
1382 if (fm != INF) {
1383 for (u = 2 + turn; u < i - turn; u++) {
1384 if (fm_tmp[u + 1] != INF) {
1385 tmp = 0;
1386 for (s = 0; s < n_seq; s++)
1387 if ((scs[s]) && (scs[s]->f))
1388 tmp += scs[s]->f(i, length - 1, u, u + 1,
1389 VRNA_DECOMP_ML_ML_ML,
1390 scs[s]->data);
1391
1392 tmp += fm;
1393 fm_tmp[u + 1] += tmp;
1394 }
1395 }
1396 } else {
1397 for (u = 2 + turn; u < i - turn; u++)
1398 fm_tmp[u + 1] = INF;
1399 }
1400 }
1401
1402 /* actually compute the entries for fM_d3 array */
1403 for (u = i + turn; u < length - turn; u++) {
1404 fm = my_fML[indx[u] + i];
1405
1406 /* skip configurations that violate (hard) constraints */
1407 if ((fm == INF) || (fm_tmp[u + 1] == INF))
1408 continue;
1409
1410 fm += fm_tmp[u + 1];
1411
1412 fM_d5[i] = MIN2(fM_d5[i], fm);
1413 }
1414 }
1415
1416 free(fm_tmp2);
1417 }
1418
1419
1420 /*
1421 **************************
1422 * Construct fM_d3 array
1423 **************************
1424 */
1425 PRIVATE INLINE void
fill_fM_d3(vrna_fold_compound_t * fc,int * fM_d3)1426 fill_fM_d3(vrna_fold_compound_t *fc,
1427 int *fM_d3)
1428 {
1429 unsigned int s, n_seq, **a2s;
1430 int *fm_tmp, *fm_tmp2, sc_base, fm, tmp, i, u,
1431 length, turn, *my_fML, *indx;
1432 vrna_md_t *md;
1433 vrna_param_t *P;
1434 vrna_hc_t *hc;
1435 vrna_sc_t *sc, **scs;
1436
1437 n_seq = (fc->type == VRNA_FC_TYPE_SINGLE) ? 1 : fc->n_seq;
1438 length = fc->length;
1439 a2s = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->a2s;
1440 P = fc->params;
1441 md = &(P->model_details);
1442 my_fML = fc->matrices->fML;
1443 hc = fc->hc;
1444 sc = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->sc : NULL;
1445 scs = (fc->type == VRNA_FC_TYPE_SINGLE) ? NULL : fc->scs;
1446 indx = fc->jindx;
1447 turn = md->min_loop_size;
1448
1449 fm_tmp2 = vrna_alloc(sizeof(int) * (length + 2));
1450 sc_base = 0;
1451
1452 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->energy_up))
1453 sc_base += sc->energy_up[1][1];
1454 else if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs))
1455 for (s = 0; s < n_seq; s++)
1456 if ((scs[s]) && (scs[s]->energy_up))
1457 sc_base += scs[s]->energy_up[a2s[s][1]][1];
1458
1459 for (i = turn + 1; i < length - turn; i++) {
1460 fm_tmp = my_fML + indx[i];
1461
1462 if (sc_base != 0) {
1463 fm_tmp = fm_tmp2;
1464
1465 for (u = 2 + turn; u < i - turn; u++)
1466 fm_tmp[u + 1] = my_fML[indx[i] + u + 1] +
1467 sc_base;
1468 }
1469
1470 if (hc->f) {
1471 if (!hc->f(1, i, 2, i, VRNA_DECOMP_ML_ML, hc->data))
1472 continue;
1473
1474 if (fm_tmp != fm_tmp2) {
1475 fm_tmp = fm_tmp2;
1476
1477 for (u = 2 + turn; u < i - turn; u++)
1478 fm_tmp[u + 1] = my_fML[indx[i] + u + 1];
1479 }
1480
1481 for (u = 2 + turn; u < i - turn; u++)
1482 if (!hc->f(2, i, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
1483 fm_tmp[u + 1] = INF;
1484 }
1485
1486 if ((fc->type == VRNA_FC_TYPE_SINGLE) && (sc) && (sc->f)) {
1487 if (fm_tmp != fm_tmp2) {
1488 fm_tmp = fm_tmp2;
1489
1490 for (u = 2 + turn; u < i - turn; u++)
1491 fm_tmp[u + 1] = my_fML[indx[i] + u + 1];
1492 }
1493
1494 fm = sc->f(1, i, 2, i, VRNA_DECOMP_ML_ML, sc->data);
1495
1496 if (fm != INF) {
1497 for (u = 2 + turn; u < i - turn; u++) {
1498 if (fm_tmp[u + 1] != INF) {
1499 tmp = sc->f(2, i, u, u + 1,
1500 VRNA_DECOMP_ML_ML_ML,
1501 sc->data);
1502 if (tmp != INF) {
1503 tmp += fm;
1504 fm_tmp[u + 1] += tmp;
1505 } else {
1506 fm_tmp[u + 1] = INF;
1507 }
1508 }
1509 }
1510 } else {
1511 for (u = 2 + turn; u < i - turn; u++)
1512 fm_tmp[u + 1] = INF;
1513 }
1514 }
1515
1516 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) && (scs)) {
1517 if (fm_tmp != fm_tmp2) {
1518 fm_tmp = fm_tmp2;
1519
1520 for (u = 2 + turn; u < i - turn; u++)
1521 fm_tmp[u + 1] = my_fML[indx[i] + u + 1];
1522 }
1523
1524 fm = 0;
1525
1526 for (s = 0; s < n_seq; s++)
1527 if ((scs[s]) && (scs[s]->f))
1528 fm += scs[s]->f(1, i, 2, i,
1529 VRNA_DECOMP_ML_ML,
1530 scs[s]->data);
1531
1532 for (u = 2 + turn; u < i - turn; u++) {
1533 if (fm_tmp[u + 1] != INF) {
1534 tmp = fm;
1535 for (s = 0; s < n_seq; s++)
1536 if ((scs[s]) && (scs[s]->f))
1537 tmp += scs[s]->f(2, i, u, u + 1,
1538 VRNA_DECOMP_ML_ML_ML,
1539 scs[s]->data);
1540
1541 fm_tmp[u + 1] += tmp;
1542 }
1543 }
1544 }
1545
1546 /* actually compute the entries for fM_d3 array */
1547 for (u = 2 + turn; u < i - turn; u++) {
1548 fm = my_fML[indx[u] + 2];
1549
1550 /* skip configurations that violate (hard) constraints */
1551 if ((fm == INF) || (fm_tmp[u + 1] == INF))
1552 continue;
1553
1554 fm += fm_tmp[u + 1];
1555
1556 fM_d3[i] = MIN2(fM_d3[i], fm);
1557 }
1558 }
1559
1560 free(fm_tmp2);
1561 }
1562
1563
1564 /**
1565 *** trace back through the "c", "f5" and "fML" arrays to get the
1566 *** base pairing list. No search for equivalent structures is done.
1567 *** This is fast, since only few structure elements are recalculated.
1568 ***
1569 *** normally s=0.
1570 *** If s>0 then s items have been already pushed onto the bt_stack
1571 **/
1572 PRIVATE int
backtrack(vrna_fold_compound_t * fc,vrna_bp_stack_t * bp_stack,sect bt_stack[],int s)1573 backtrack(vrna_fold_compound_t *fc,
1574 vrna_bp_stack_t *bp_stack,
1575 sect bt_stack[],
1576 int s)
1577 {
1578 char backtrack_type;
1579 int i, j, ij, k, length, b, *my_c, *indx, noLP, *pscore, ret;
1580 vrna_param_t *P;
1581
1582 ret = 1;
1583 b = 0;
1584 length = fc->length;
1585 my_c = fc->matrices->c;
1586 indx = fc->jindx;
1587 P = fc->params;
1588 noLP = P->model_details.noLP;
1589 pscore = fc->pscore; /* covariance scores for comparative structure prediction */
1590 backtrack_type = P->model_details.backtrack_type;
1591
1592 if (s == 0) {
1593 bt_stack[++s].i = 1;
1594 bt_stack[s].j = length;
1595 bt_stack[s].ml = (backtrack_type == 'M') ? 1 : ((backtrack_type == 'C') ? 2 : 0);
1596 }
1597
1598 while (s > 0) {
1599 int ml, cij;
1600 int canonical = 1; /* (i,j) closes a canonical structure */
1601
1602 /* pop one element from stack */
1603 i = bt_stack[s].i;
1604 j = bt_stack[s].j;
1605 ml = bt_stack[s--].ml;
1606
1607 switch (ml) {
1608 /* backtrack in f5 */
1609 case 0:
1610 {
1611 int p, q;
1612 if (vrna_BT_ext_loop_f5(fc, &j, &p, &q, bp_stack, &b)) {
1613 if (j > 0) {
1614 bt_stack[++s].i = 1;
1615 bt_stack[s].j = j;
1616 bt_stack[s].ml = 0;
1617 }
1618
1619 if (p > 0) {
1620 i = p;
1621 j = q;
1622 goto repeat1;
1623 }
1624
1625 continue;
1626 } else {
1627 vrna_message_warning("backtracking failed in f5, segment [%d,%d]\n", i, j);
1628 ret = 0;
1629 goto backtrack_exit;
1630 }
1631 }
1632 break;
1633
1634 /* trace back in fML array */
1635 case 1:
1636 {
1637 int p, q, comp1, comp2;
1638 if (vrna_BT_mb_loop_split(fc, &i, &j, &p, &q, &comp1, &comp2, bp_stack, &b)) {
1639 if (i > 0) {
1640 bt_stack[++s].i = i;
1641 bt_stack[s].j = j;
1642 bt_stack[s].ml = comp1;
1643 }
1644
1645 if (p > 0) {
1646 bt_stack[++s].i = p;
1647 bt_stack[s].j = q;
1648 bt_stack[s].ml = comp2;
1649 }
1650
1651 continue;
1652 } else {
1653 ret = 0;
1654 goto backtrack_exit;
1655 }
1656 }
1657 break;
1658
1659 /* backtrack in c */
1660 case 2:
1661 bp_stack[++b].i = i;
1662 bp_stack[b].j = j;
1663 goto repeat1;
1664
1665 break;
1666
1667 default:
1668 ret = 0;
1669 goto backtrack_exit;
1670 }
1671
1672 repeat1:
1673
1674 /*----- begin of "repeat:" -----*/
1675 ij = indx[j] + i;
1676
1677 if (canonical)
1678 cij = my_c[ij];
1679
1680 if (noLP) {
1681 if (vrna_BT_stack(fc, &i, &j, &cij, bp_stack, &b)) {
1682 canonical = 0;
1683 goto repeat1;
1684 }
1685 }
1686
1687 canonical = 1;
1688
1689 if (fc->type == VRNA_FC_TYPE_COMPARATIVE)
1690 cij += pscore[indx[j] + i];
1691
1692 if (vrna_BT_hp_loop(fc, i, j, cij, bp_stack, &b))
1693 continue;
1694
1695 if (vrna_BT_int_loop(fc, &i, &j, cij, bp_stack, &b)) {
1696 if (i < 0)
1697 continue;
1698 else
1699 goto repeat1;
1700 }
1701
1702 /* (i.j) must close a multi-loop */
1703 int comp1, comp2;
1704
1705 if (vrna_BT_mb_loop(fc, &i, &j, &k, cij, &comp1, &comp2)) {
1706 bt_stack[++s].i = i;
1707 bt_stack[s].j = k;
1708 bt_stack[s].ml = comp1;
1709 bt_stack[++s].i = k + 1;
1710 bt_stack[s].j = j;
1711 bt_stack[s].ml = comp2;
1712 } else {
1713 vrna_message_warning("backtracking failed in repeat, segment [%d,%d]\n", i, j);
1714 ret = 0;
1715 goto backtrack_exit;
1716 }
1717
1718 /* end of repeat: --------------------------------------------------*/
1719 } /* end of infinite while loop */
1720
1721 backtrack_exit:
1722
1723 bp_stack[0].i = b; /* save the total number of base pairs */
1724
1725 return ret;
1726 }
1727
1728
1729 PRIVATE INLINE int
decompose_pair(vrna_fold_compound_t * fc,int i,int j,struct aux_arrays * aux)1730 decompose_pair(vrna_fold_compound_t *fc,
1731 int i,
1732 int j,
1733 struct aux_arrays *aux)
1734 {
1735 unsigned char hc_decompose;
1736 unsigned int n;
1737 int e, new_c, energy, stackEnergy, ij, dangle_model, noLP,
1738 *DMLi1, *DMLi2, *cc, *cc1;
1739
1740 n = fc->length;
1741 ij = fc->jindx[j] + i;
1742 dangle_model = fc->params->model_details.dangles;
1743 noLP = fc->params->model_details.noLP;
1744 hc_decompose = fc->hc->mx[n * i + j];
1745 DMLi1 = aux->DMLi1;
1746 DMLi2 = aux->DMLi2;
1747 cc = aux->cc;
1748 cc1 = aux->cc1;
1749 e = INF;
1750
1751 /* do we evaluate this pair? */
1752 if (hc_decompose) {
1753 new_c = INF;
1754
1755 /* check for hairpin loop */
1756 energy = vrna_E_hp_loop(fc, i, j);
1757 new_c = MIN2(new_c, energy);
1758
1759 /* check for multibranch loops */
1760 energy = vrna_E_mb_loop_fast(fc, i, j, DMLi1, DMLi2);
1761 new_c = MIN2(new_c, energy);
1762
1763 if (dangle_model == 3) {
1764 /* coaxial stacking */
1765 energy = vrna_E_mb_loop_stack(fc, i, j);
1766 new_c = MIN2(new_c, energy);
1767 }
1768
1769 /* check for interior loops */
1770 energy = vrna_E_int_loop(fc, i, j);
1771 new_c = MIN2(new_c, energy);
1772
1773 /* remember stack energy for --noLP option */
1774 if (noLP) {
1775 stackEnergy = vrna_E_stack(fc, i, j);
1776 new_c = MIN2(new_c, cc1[j - 1] + stackEnergy);
1777 cc[j] = new_c;
1778 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) &&
1779 (cc[j] != INF))
1780 cc[j] -= fc->pscore[ij];
1781
1782 e = cc1[j - 1] + stackEnergy;
1783 } else {
1784 e = new_c;
1785 }
1786
1787 /* finally, check for auxiliary grammar rule(s) */
1788 if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_c)) {
1789 energy = fc->aux_grammar->cb_aux_c(fc, i, j, fc->aux_grammar->data);
1790 new_c = MIN2(new_c, energy);
1791 }
1792
1793 if ((fc->type == VRNA_FC_TYPE_COMPARATIVE) &&
1794 (e != INF))
1795 e -= fc->pscore[ij];
1796 } /* end >> if (pair) << */
1797
1798 return e;
1799 }
1800
1801
1802 PRIVATE INLINE struct aux_arrays *
get_aux_arrays(unsigned int length)1803 get_aux_arrays(unsigned int length)
1804 {
1805 unsigned int j;
1806 struct aux_arrays *aux = (struct aux_arrays *)vrna_alloc(sizeof(struct aux_arrays));
1807
1808 aux->cc = (int *)vrna_alloc(sizeof(int) * (length + 2)); /* auxilary arrays for canonical structures */
1809 aux->cc1 = (int *)vrna_alloc(sizeof(int) * (length + 2)); /* auxilary arrays for canonical structures */
1810 aux->Fmi = (int *)vrna_alloc(sizeof(int) * (length + 1)); /* holds row i of fML (avoids jumps in memory) */
1811 aux->DMLi = (int *)vrna_alloc(sizeof(int) * (length + 1)); /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j]) */
1812 aux->DMLi1 = (int *)vrna_alloc(sizeof(int) * (length + 1)); /* MIN(fML[i+1,k]+fML[k+1,j]) */
1813 aux->DMLi2 = (int *)vrna_alloc(sizeof(int) * (length + 1)); /* MIN(fML[i+2,k]+fML[k+1,j]) */
1814
1815 /* prefill helper arrays */
1816 for (j = 0; j <= length; j++)
1817 aux->Fmi[j] = aux->DMLi[j] = aux->DMLi1[j] = aux->DMLi2[j] = INF;
1818
1819 return aux;
1820 }
1821
1822
1823 PRIVATE INLINE void
rotate_aux_arrays(struct aux_arrays * aux,unsigned int length)1824 rotate_aux_arrays(struct aux_arrays *aux,
1825 unsigned int length)
1826 {
1827 unsigned int j;
1828 int *FF;
1829
1830 FF = aux->DMLi2;
1831 aux->DMLi2 = aux->DMLi1;
1832 aux->DMLi1 = aux->DMLi;
1833 aux->DMLi = FF;
1834 FF = aux->cc1;
1835 aux->cc1 = aux->cc;
1836 aux->cc = FF;
1837 for (j = 1; j <= length; j++)
1838 aux->cc[j] = aux->Fmi[j] = aux->DMLi[j] = INF;
1839 }
1840
1841
1842 PRIVATE INLINE void
free_aux_arrays(struct aux_arrays * aux)1843 free_aux_arrays(struct aux_arrays *aux)
1844 {
1845 free(aux->cc);
1846 free(aux->cc1);
1847 free(aux->Fmi);
1848 free(aux->DMLi);
1849 free(aux->DMLi1);
1850 free(aux->DMLi2);
1851 free(aux);
1852 }
1853