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 <math.h>
16 #include <float.h> /* #defines FLT_MAX ... */
17 #include <limits.h>
18
19 #include "ViennaRNA/utils/basic.h"
20 #include "ViennaRNA/params/default.h"
21 #include "ViennaRNA/fold_vars.h"
22 #include "ViennaRNA/loops/all.h"
23 #include "ViennaRNA/gquad.h"
24 #include "ViennaRNA/constraints/hard.h"
25 #include "ViennaRNA/constraints/soft.h"
26 #include "ViennaRNA/mfe.h"
27 #include "ViennaRNA/part_func.h"
28
29 #ifdef _OPENMP
30 #include <omp.h>
31 #endif
32
33 /*
34 #################################
35 # GLOBAL VARIABLES #
36 #################################
37 */
38
39 /*
40 #################################
41 # PRIVATE VARIABLES #
42 #################################
43 */
44
45 /*
46 #################################
47 # PRIVATE FUNCTION DECLARATIONS #
48 #################################
49 */
50 PRIVATE int
51 fill_arrays(vrna_fold_compound_t *fc);
52
53
54 PRIVATE void
55 postprocess_circular(vrna_fold_compound_t *fc);
56
57
58 PRIVATE FLT_OR_DBL
59 decompose_pair(vrna_fold_compound_t *fc,
60 int i,
61 int j,
62 vrna_mx_pf_aux_ml_t aux_mx_ml);
63
64
65 /*
66 #################################
67 # BEGIN OF FUNCTION DEFINITIONS #
68 #################################
69 */
70 PUBLIC float
vrna_pf(vrna_fold_compound_t * fc,char * structure)71 vrna_pf(vrna_fold_compound_t *fc,
72 char *structure)
73 {
74 int n;
75 FLT_OR_DBL Q;
76 double free_energy;
77 vrna_md_t *md;
78 vrna_exp_param_t *params;
79 vrna_mx_pf_t *matrices;
80
81 free_energy = (float)(INF / 100.);
82
83 if (fc) {
84 /* make sure, everything is set up properly to start partition function computations */
85 if (!vrna_fold_compound_prepare(fc, VRNA_OPTION_PF)) {
86 vrna_message_warning("vrna_pf@part_func.c: Failed to prepare vrna_fold_compound");
87 return free_energy;
88 }
89
90 n = fc->length;
91 params = fc->exp_params;
92 matrices = fc->exp_matrices;
93 md = &(params->model_details);
94
95 #ifdef _OPENMP
96 /* Explicitly turn off dynamic threads */
97 omp_set_dynamic(0);
98 #endif
99
100 #ifdef SUN4
101 nonstandard_arithmetic();
102 #elif defined(HP9)
103 fpsetfastmode(1);
104 #endif
105
106 /* call user-defined recursion status callback function */
107 if (fc->stat_cb)
108 fc->stat_cb(VRNA_STATUS_PF_PRE, fc->auxdata);
109
110 /* call user-defined grammar pre-condition callback function */
111 if ((fc->aux_grammar) && (fc->aux_grammar->cb_proc))
112 fc->aux_grammar->cb_proc(fc, VRNA_STATUS_PF_PRE, fc->aux_grammar->data);
113
114 if (!fill_arrays(fc)) {
115 #ifdef SUN4
116 standard_arithmetic();
117 #elif defined(HP9)
118 fpsetfastmode(0);
119 #endif
120 return (float)(INF / 100.);
121 }
122
123 if (md->circ)
124 /* do post processing step for circular RNAs */
125 postprocess_circular(fc);
126
127 /* calculate base pairing probability matrix (bppm) */
128 if (md->compute_bpp) {
129 vrna_pairing_probs(fc, structure);
130
131 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
132
133 /*
134 * Backward compatibility:
135 * This block may be removed if deprecated functions
136 * relying on the global variable "pr" vanish from within the package!
137 */
138 pr = matrices->probs;
139
140 #endif
141 }
142
143 /* call user-defined recursion status callback function */
144 if (fc->stat_cb)
145 fc->stat_cb(VRNA_STATUS_PF_POST, fc->auxdata);
146
147 /* call user-defined grammar post-condition callback function */
148 if ((fc->aux_grammar) && (fc->aux_grammar->cb_proc))
149 fc->aux_grammar->cb_proc(fc, VRNA_STATUS_PF_POST, fc->aux_grammar->data);
150
151 switch (md->backtrack_type) {
152 case 'C':
153 Q = matrices->qb[fc->iindx[1] - n];
154 break;
155
156 case 'M':
157 Q = matrices->qm[fc->iindx[1] - n];
158 break;
159
160 default:
161 Q = (md->circ) ? matrices->qo : matrices->q[fc->iindx[1] - n];
162 break;
163 }
164
165 /* ensemble free energy in Kcal/mol */
166 if (Q <= FLT_MIN)
167 vrna_message_warning("pf_scale too large");
168
169 free_energy = (-log(Q) - n * log(params->pf_scale)) *
170 params->kT /
171 1000.0;
172
173 if (fc->type == VRNA_FC_TYPE_COMPARATIVE)
174 free_energy /= fc->n_seq;
175
176 #ifdef SUN4
177 standard_arithmetic();
178 #elif defined(HP9)
179 fpsetfastmode(0);
180 #endif
181 }
182
183 return free_energy;
184 }
185
186
187 PUBLIC vrna_dimer_pf_t
vrna_pf_dimer(vrna_fold_compound_t * fc,char * structure)188 vrna_pf_dimer(vrna_fold_compound_t *fc,
189 char *structure)
190 {
191 unsigned int *so, *se, *ss;
192 int n;
193 FLT_OR_DBL Q;
194 vrna_dimer_pf_t X;
195 double free_energy;
196 char *sequence;
197 vrna_md_t *md;
198 vrna_exp_param_t *params;
199 vrna_mx_pf_t *matrices;
200
201 if (!vrna_fold_compound_prepare(fc, VRNA_OPTION_PF | VRNA_OPTION_HYBRID)) {
202 vrna_message_warning("vrna_pf_dimer@part_func_co.c: Failed to prepare vrna_fold_compound");
203 X.FA = X.FB = X.FAB = X.F0AB = X.FcAB = 0;
204 return X;
205 }
206
207 params = fc->exp_params;
208 n = fc->length;
209 so = fc->strand_order;
210 se = fc->strand_end;
211 ss = fc->strand_start;
212 md = &(params->model_details);
213 matrices = fc->exp_matrices;
214 sequence = fc->sequence;
215
216 #ifdef _OPENMP
217 /* Explicitly turn off dynamic threads */
218 omp_set_dynamic(0);
219 #endif
220
221 #ifdef SUN4
222 nonstandard_arithmetic();
223 #elif defined(HP9)
224 fpsetfastmode(1);
225 #endif
226
227 /* hard code min_loop_size to 0, since we can not be sure yet that this is already the case */
228 md->min_loop_size = 0;
229
230 /* call user-defined recursion status callback function */
231 if (fc->stat_cb)
232 fc->stat_cb(VRNA_STATUS_PF_PRE, fc->auxdata);
233
234 if (!fill_arrays(fc)) {
235 X.FA = X.FB = X.FAB = X.F0AB = (float)(INF / 100.);
236 X.FcAB = 0;
237
238 #ifdef SUN4
239 standard_arithmetic();
240 #elif defined(HP9)
241 fpsetfastmode(0);
242 #endif
243
244 return X;
245 }
246
247 /* call user-defined recursion status callback function */
248 if (fc->stat_cb)
249 fc->stat_cb(VRNA_STATUS_PF_POST, fc->auxdata);
250
251 if (md->backtrack_type == 'C')
252 Q = matrices->qb[fc->iindx[1] - n];
253 else if (md->backtrack_type == 'M')
254 Q = matrices->qm[fc->iindx[1] - n];
255 else
256 Q = matrices->q[fc->iindx[1] - n];
257
258 /* ensemble free energy in Kcal/mol */
259 if (Q <= FLT_MIN)
260 vrna_message_warning("pf_scale too large");
261
262 free_energy = (-log(Q) - n * log(params->pf_scale)) * params->kT / 1000.0;
263 /* in case we abort because of floating point errors */
264 if (n > 1600)
265 vrna_message_info(stderr, "free energy = %8.2f", free_energy);
266
267 /* probability of molecules being bound together */
268
269 /*
270 * Computation of "real" Partition function
271 * Need that for concentrations
272 */
273 if (fc->strands > 1) {
274 double kT, QAB, QToT, Qzero;
275 kT = params->kT / 1000.0;
276 Qzero = matrices->q[fc->iindx[1] - n];
277 QAB =
278 (matrices->q[fc->iindx[1] - n] - matrices->q[fc->iindx[1] - se[so[0]]] *
279 matrices->q[fc->iindx[ss[so[1]]] - n]) * params->expDuplexInit;
280 /*correction for symmetry*/
281 if ((n - 2 * se[so[0]]) == 0)
282 if ((strncmp(sequence, sequence + se[so[0]], se[so[0]])) == 0)
283 QAB /= 2;
284
285 QToT = matrices->q[fc->iindx[1] - se[so[0]]] *
286 matrices->q[fc->iindx[ss[so[1]]] - n] + QAB;
287 X.FAB = -kT * (log(QToT) + n * log(params->pf_scale));
288 X.F0AB = -kT * (log(Qzero) + n * log(params->pf_scale));
289 X.FcAB = (QAB > 1e-17) ? -kT * (log(QAB) + n * log(params->pf_scale)) : 999;
290 X.FA = -kT *
291 (log(matrices->q[fc->iindx[1] - se[so[0]]]) + (se[so[0]]) *
292 log(params->pf_scale));
293 X.FB = -kT *
294 (log(matrices->q[fc->iindx[ss[so[1]]] - n]) + (n - ss[so[1]] + 1) *
295 log(params->pf_scale));
296
297 /* printf("QAB=%.9f\tQtot=%.9f\n",QAB/scale[n],QToT/scale[n]); */
298 } else {
299 X.FA = X.FB = X.FAB = X.F0AB = free_energy;
300 X.FcAB = 0;
301 }
302
303 /* backtracking to construct binding probabilities of pairs */
304 if (md->compute_bpp) {
305 vrna_pairing_probs(fc, structure);
306
307 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
308
309 /*
310 * Backward compatibility:
311 * This block may be removed if deprecated functions
312 * relying on the global variable "pr" vanish from within the package!
313 */
314 pr = fc->exp_matrices->probs;
315
316 #endif
317 }
318
319 #ifdef SUN4
320 standard_arithmetic();
321 #elif defined(HP9)
322 fpsetfastmode(0);
323 #endif
324
325 return X;
326 }
327
328
329 PUBLIC int
vrna_pf_float_precision(void)330 vrna_pf_float_precision(void)
331 {
332 return sizeof(FLT_OR_DBL) == sizeof(float);
333 }
334
335
336 /*
337 #################################
338 # STATIC helper functions below #
339 #################################
340 */
341 PRIVATE int
fill_arrays(vrna_fold_compound_t * fc)342 fill_arrays(vrna_fold_compound_t *fc)
343 {
344 int n, i, j, k, ij, d, *my_iindx, *jindx, with_gquad, turn,
345 with_ud;
346 FLT_OR_DBL temp, Qmax, *q, *qb, *qm, *qm1, *q1k, *qln;
347 double max_real;
348 vrna_ud_t *domains_up;
349 vrna_md_t *md;
350 vrna_mx_pf_t *matrices;
351 vrna_mx_pf_aux_el_t aux_mx_el;
352 vrna_mx_pf_aux_ml_t aux_mx_ml;
353 vrna_exp_param_t *pf_params;
354
355 n = fc->length;
356 my_iindx = fc->iindx;
357 jindx = fc->jindx;
358 matrices = fc->exp_matrices;
359 pf_params = fc->exp_params;
360 domains_up = fc->domains_up;
361 q = matrices->q;
362 qb = matrices->qb;
363 qm = matrices->qm;
364 qm1 = matrices->qm1;
365 q1k = matrices->q1k;
366 qln = matrices->qln;
367 md = &(pf_params->model_details);
368 with_gquad = md->gquad;
369 turn = md->min_loop_size;
370
371 with_ud = (domains_up && domains_up->exp_energy_cb && (!(fc->type == VRNA_FC_TYPE_COMPARATIVE)));
372 Qmax = 0;
373
374 max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
375
376 if (with_ud && domains_up->exp_prod_cb)
377 domains_up->exp_prod_cb(fc, domains_up->data);
378
379 /* no G-Quadruplexes for comparative partition function (yet) */
380 if (with_gquad) {
381 free(fc->exp_matrices->G);
382 fc->exp_matrices->G = NULL;
383
384 switch (fc->type) {
385 case VRNA_FC_TYPE_SINGLE:
386 fc->exp_matrices->G = get_gquad_pf_matrix(fc->sequence_encoding2,
387 fc->exp_matrices->scale,
388 fc->exp_params);
389 break;
390
391 case VRNA_FC_TYPE_COMPARATIVE:
392 fc->exp_matrices->G = get_gquad_pf_matrix_comparative(fc->length,
393 fc->S_cons,
394 fc->S,
395 fc->a2s,
396 fc->exp_matrices->scale,
397 fc->n_seq,
398 fc->exp_params);
399 break;
400 }
401 }
402
403 /* init auxiliary arrays for fast exterior/multibranch loops */
404 aux_mx_el = vrna_exp_E_ext_fast_init(fc);
405 aux_mx_ml = vrna_exp_E_ml_fast_init(fc);
406
407 /*array initialization ; qb,qm,q
408 * qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */
409 for (d = 0; d <= turn; d++)
410 for (i = 1; i <= n - d; i++) {
411 j = i + d;
412 ij = my_iindx[i] - j;
413 qb[ij] = 0.0;
414 }
415
416 for (j = turn + 2; j <= n; j++) {
417 for (i = j - turn - 1; i >= 1; i--) {
418 ij = my_iindx[i] - j;
419
420 qb[ij] = decompose_pair(fc, i, j, aux_mx_ml);
421
422 /* Multibranch loop */
423 qm[ij] = vrna_exp_E_ml_fast(fc, i, j, aux_mx_ml);
424
425 if (qm1) {
426 temp = vrna_exp_E_ml_fast_qqm(aux_mx_ml)[i]; /* for stochastic backtracking and circfold */
427
428 /* apply auxiliary grammar rule for multibranch loop (M1) case */
429 if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_m1))
430 temp += fc->aux_grammar->cb_aux_exp_m1(fc, i, j, fc->aux_grammar->data);
431
432 qm1[jindx[j] + i] = temp;
433 }
434
435 /* Exterior loop */
436 q[ij] = vrna_exp_E_ext_fast(fc, i, j, aux_mx_el);
437
438 /* apply auxiliary grammar rule (storage takes place in user-defined data structure */
439 if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp))
440 fc->aux_grammar->cb_aux_exp(fc, i, j, fc->aux_grammar->data);
441
442 if (q[ij] > Qmax) {
443 Qmax = q[ij];
444 if (Qmax > max_real / 10.)
445 vrna_message_warning("Q close to overflow: %d %d %g", i, j, q[ij]);
446 }
447
448 if (q[ij] >= max_real) {
449 vrna_message_warning("overflow while computing partition function for segment q[%d,%d]\n"
450 "use larger pf_scale", i, j);
451
452 vrna_exp_E_ml_fast_free(aux_mx_ml);
453 vrna_exp_E_ext_fast_free(aux_mx_el);
454
455 return 0; /* failure */
456 }
457 }
458
459 /* rotate auxiliary arrays */
460 vrna_exp_E_ext_fast_rotate(aux_mx_el);
461 vrna_exp_E_ml_fast_rotate(aux_mx_ml);
462 }
463
464 /* prefill linear qln, q1k arrays */
465 if (q1k && qln) {
466 for (k = 1; k <= n; k++) {
467 q1k[k] = q[my_iindx[1] - k];
468 qln[k] = q[my_iindx[k] - n];
469 }
470 q1k[0] = 1.0;
471 qln[n + 1] = 1.0;
472 }
473
474 /* free memory occupied by auxiliary arrays for fast exterior/multibranch loops */
475 vrna_exp_E_ml_fast_free(aux_mx_ml);
476 vrna_exp_E_ext_fast_free(aux_mx_el);
477
478 return 1;
479 }
480
481
482 PRIVATE FLT_OR_DBL
decompose_pair(vrna_fold_compound_t * fc,int i,int j,vrna_mx_pf_aux_ml_t aux_mx_ml)483 decompose_pair(vrna_fold_compound_t *fc,
484 int i,
485 int j,
486 vrna_mx_pf_aux_ml_t aux_mx_ml)
487 {
488 unsigned int n;
489 int *jindx, *pscore;
490 FLT_OR_DBL contribution;
491 double kTn;
492 vrna_hc_t *hc;
493
494 contribution = 0.;
495 n = fc->length;
496 hc = fc->hc;
497
498 if (hc->mx[j * n + i]) {
499 /* process hairpin loop(s) */
500 contribution += vrna_exp_E_hp_loop(fc, i, j);
501 /* process interior loop(s) */
502 contribution += vrna_exp_E_int_loop(fc, i, j);
503 /* process multibranch loop(s) */
504 contribution += vrna_exp_E_mb_loop_fast(fc, i, j, aux_mx_ml);
505
506 if ((fc->aux_grammar) && (fc->aux_grammar->cb_aux_exp_c))
507 contribution += fc->aux_grammar->cb_aux_exp_c(fc, i, j, fc->aux_grammar->data);
508
509 if (fc->type == VRNA_FC_TYPE_COMPARATIVE) {
510 jindx = fc->jindx;
511 pscore = fc->pscore;
512 kTn = fc->exp_params->kT / 10.; /* kT in cal/mol */
513 contribution *= exp(pscore[jindx[j] + i] / kTn);
514 }
515 }
516
517 return contribution;
518 }
519
520
521 /*
522 * calculate partition function for circular case
523 * NOTE: this is the postprocessing step ONLY
524 * You have to call fill_arrays first to calculate
525 * complete circular case!!!
526 */
527 PRIVATE void
postprocess_circular(vrna_fold_compound_t * fc)528 postprocess_circular(vrna_fold_compound_t *fc)
529 {
530 unsigned int **a2s;
531 int u, p, q, k, turn, n, *my_iindx, *jindx, s;
532 FLT_OR_DBL *scale, *qb, *qm, *qm1, *qm2, qo, qho, qio, qmo,
533 qbt1, qot, expMLclosing, n_seq;
534 unsigned char eval;
535 vrna_exp_param_t *pf_params;
536 vrna_mx_pf_t *matrices;
537 vrna_hc_t *hc;
538 vrna_sc_t *sc, **scs;
539
540 n = fc->length;
541 n_seq = (fc->type == VRNA_FC_TYPE_SINGLE) ? 1 : fc->n_seq;
542 matrices = fc->exp_matrices;
543 my_iindx = fc->iindx;
544 jindx = fc->jindx;
545 pf_params = fc->exp_params;
546 hc = fc->hc;
547 qb = matrices->qb;
548 qm = matrices->qm;
549 qm1 = matrices->qm1;
550 qm2 = matrices->qm2;
551 scale = matrices->scale;
552 expMLclosing = pf_params->expMLclosing;
553 turn = pf_params->model_details.min_loop_size;
554 hc = fc->hc;
555 sc = (fc->type == VRNA_FC_TYPE_SINGLE) ? fc->sc : NULL;
556 scs = (fc->type == VRNA_FC_TYPE_COMPARATIVE) ? fc->scs : NULL;
557 a2s = (fc->type == VRNA_FC_TYPE_COMPARATIVE) ? fc->a2s : NULL;
558 qo = qho = qio = qmo = 0.;
559
560 for (p = 1; p < n; p++) {
561 for (q = p + turn + 1; q <= n; q++) {
562 u = n - q + p - 1;
563 if (u < turn)
564 continue;
565
566 /* 1. get exterior hairpin contribution */
567 qho += qb[my_iindx[p] - q] *
568 vrna_exp_E_hp_loop(fc, q, p);
569
570 /* 2. get exterior interior loop contribution */
571 qio += qb[my_iindx[p] - q] *
572 vrna_exp_E_int_loop(fc, q, p);
573 }
574 } /* end of pq double loop */
575
576 /* 3. Multiloops */
577
578 /* construct qm2 matrix for exterior multibranch loop computation */
579 if (hc->f) {
580 switch (fc->type) {
581 case VRNA_FC_TYPE_SINGLE:
582 if ((sc) && (sc->exp_f)) {
583 for (k = 1; k < n - turn - 1; k++) {
584 qot = 0.;
585
586 for (u = k + turn + 1; u < n - turn - 1; u++)
587 if (hc->f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data)) {
588 qot += qm1[jindx[u] + k] *
589 qm1[jindx[n] + (u + 1)] *
590 sc->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
591 }
592
593 qm2[k] = qot;
594 }
595 } else {
596 for (k = 1; k < n - turn - 1; k++) {
597 qot = 0.;
598
599 for (u = k + turn + 1; u < n - turn - 1; u++)
600 if (hc->f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
601 qot += qm1[jindx[u] + k] *
602 qm1[jindx[n] + (u + 1)];
603
604 qm2[k] = qot;
605 }
606 }
607
608 break;
609
610 case VRNA_FC_TYPE_COMPARATIVE:
611 if (scs) {
612 for (k = 1; k < n - turn - 1; k++) {
613 qbt1 = 0.;
614
615 for (u = k + turn + 1; u < n - turn - 1; u++) {
616 if (hc->f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data)) {
617 qot = qm1[jindx[u] + k] *
618 qm1[jindx[n] + (u + 1)];
619
620 for (s = 0; s < n_seq; s++)
621 if ((scs[s]) && scs[s]->exp_f)
622 qot *= scs[s]->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
623
624 qbt1 += qot;
625 }
626 }
627
628 qm2[k] = qbt1;
629 }
630 } else {
631 for (k = 1; k < n - turn - 1; k++) {
632 qot = 0.;
633
634 for (u = k + turn + 1; u < n - turn - 1; u++)
635 if (hc->f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
636 qot += qm1[jindx[u] + k] *
637 qm1[jindx[n] + (u + 1)];
638
639 qm2[k] = qot;
640 }
641 }
642
643 break;
644 }
645 } else {
646 switch (fc->type) {
647 case VRNA_FC_TYPE_SINGLE:
648 if ((sc) && (sc->exp_f)) {
649 for (k = 1; k < n - turn - 1; k++) {
650 qot = 0.;
651
652 for (u = k + turn + 1; u < n - turn - 1; u++)
653 qot += qm1[jindx[u] + k] *
654 qm1[jindx[n] + (u + 1)] *
655 sc->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
656
657 qm2[k] = qot;
658 }
659 } else {
660 for (k = 1; k < n - turn - 1; k++) {
661 qot = 0.;
662
663 for (u = k + turn + 1; u < n - turn - 1; u++)
664 qot += qm1[jindx[u] + k] *
665 qm1[jindx[n] + (u + 1)];
666
667 qm2[k] = qot;
668 }
669 }
670
671 break;
672
673 case VRNA_FC_TYPE_COMPARATIVE:
674 if (scs) {
675 for (k = 1; k < n - turn - 1; k++) {
676 qbt1 = 0.;
677
678 for (u = k + turn + 1; u < n - turn - 1; u++) {
679 qot = qm1[jindx[u] + k] *
680 qm1[jindx[n] + (u + 1)];
681
682 for (s = 0; s < n_seq; s++)
683 if ((scs[s]) && (scs[s]->exp_f))
684 qot *= scs[s]->exp_f(k, n, u, u + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
685
686 qbt1 += qot;
687 }
688
689 qm2[k] = qbt1;
690 }
691 } else {
692 for (k = 1; k < n - turn - 1; k++) {
693 qot = 0.;
694
695 for (u = k + turn + 1; u < n - turn - 1; u++)
696 qot += qm1[jindx[u] + k] *
697 qm1[jindx[n] + (u + 1)];
698
699 qm2[k] = qot;
700 }
701 }
702
703 break;
704 }
705 }
706
707 qbt1 = 0.;
708 /* go through exterior multibranch loop configurations */
709 if (hc->f) {
710 switch (fc->type) {
711 case VRNA_FC_TYPE_SINGLE:
712 if ((sc) && (sc->exp_f)) {
713 for (k = turn + 2; k < n - 2 * turn - 3; k++)
714 if (hc->f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
715 qbt1 += qm[my_iindx[1] - k] *
716 qm2[k + 1] *
717 sc->exp_f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
718 } else {
719 for (k = turn + 2; k < n - 2 * turn - 3; k++)
720 if (hc->f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
721 qbt1 += qm[my_iindx[1] - k] *
722 qm2[k + 1];
723 }
724
725 qbt1 *= expMLclosing;
726 break;
727
728 case VRNA_FC_TYPE_COMPARATIVE:
729 if (scs) {
730 for (k = turn + 2; k < n - 2 * turn - 3; k++)
731 if (hc->f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, hc->data)) {
732 qot = qm[my_iindx[1] - k] *
733 qm2[k + 1];
734
735 for (s = 0; s < n_seq; s++)
736 if ((scs[s]) && (scs[s]->exp_f))
737 qot *= scs[s]->exp_f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
738
739 qbt1 += qot;
740 }
741 } else {
742 for (k = turn + 2; k < n - 2 * turn - 3; k++)
743 if (hc->f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, hc->data))
744 qbt1 += qm[my_iindx[1] - k] *
745 qm2[k + 1];
746 }
747
748 qbt1 *= pow(expMLclosing, fc->n_seq);
749 break;
750 }
751 } else {
752 switch (fc->type) {
753 case VRNA_FC_TYPE_SINGLE:
754 if ((sc) && (sc->exp_f)) {
755 for (k = turn + 2; k < n - 2 * turn - 3; k++)
756 qbt1 += qm[my_iindx[1] - k] *
757 qm2[k + 1] *
758 sc->exp_f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, sc->data);
759 } else {
760 for (k = turn + 2; k < n - 2 * turn - 3; k++)
761 qbt1 += qm[my_iindx[1] - k] *
762 qm2[k + 1];
763 }
764
765 qbt1 *= expMLclosing;
766 break;
767
768 case VRNA_FC_TYPE_COMPARATIVE:
769 if (scs) {
770 for (k = turn + 2; k < n - 2 * turn - 3; k++) {
771 qot = qm[my_iindx[1] - k] *
772 qm2[k + 1];
773
774 for (s = 0; s < n_seq; s++)
775 if ((scs[s]) && (scs[s]->exp_f))
776 qot *= scs[s]->exp_f(1, n, k, k + 1, VRNA_DECOMP_ML_ML_ML, scs[s]->data);
777
778 qbt1 += qot;
779 }
780 } else {
781 for (k = turn + 2; k < n - 2 * turn - 3; k++)
782 qbt1 += qm[my_iindx[1] - k] *
783 qm2[k + 1];
784 }
785
786 qbt1 *= pow(expMLclosing, fc->n_seq);
787 break;
788 }
789 }
790
791 qmo += qbt1;
792
793 /* add an additional pf of 1.0 to take the open chain into account too */
794 eval = (hc->up_ext[1] >= n) ? 1 : 0;
795 if (hc->f)
796 eval = (hc->f(1, n, 1, n, VRNA_DECOMP_EXT_UP, hc->data)) ? eval : 0;
797
798 if (eval) {
799 qbt1 = scale[n];
800
801 switch (fc->type) {
802 case VRNA_FC_TYPE_SINGLE:
803 if (sc) {
804 if (sc->exp_energy_up)
805 qbt1 *= sc->exp_energy_up[1][n];
806
807 if (sc->exp_f)
808 qbt1 *= sc->exp_f(1, n, 1, n, VRNA_DECOMP_EXT_UP, sc->data);
809 }
810
811 break;
812
813 case VRNA_FC_TYPE_COMPARATIVE:
814 if (scs) {
815 for (s = 0; s < fc->n_seq; s++)
816 if (scs[s])
817 if (scs[s]->energy_up)
818 qbt1 *= scs[s]->exp_energy_up[1][a2s[s][n]];
819 }
820
821 break;
822 }
823 qo += qbt1;
824 }
825
826 qo += qho + qio + qmo;
827
828 matrices->qo = qo;
829 matrices->qho = qho;
830 matrices->qio = qio;
831 matrices->qmo = qmo;
832 }
833