1 /*
2 * local pair probabilities for RNA secondary structures
3 *
4 * Stephan Bernhart, Ivo L Hofacker
5 * Vienna RNA package
6 */
7 /*
8 * todo: compute energy z-score for each window
9 *
10 */
11
12 #ifdef HAVE_CONFIG_H
13 #include "config.h"
14 #endif
15
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <math.h>
20 #include <float.h> /* #defines FLT_MAX ... */
21 #include "ViennaRNA/datastructures/basic.h"
22 #include "ViennaRNA/utils/basic.h"
23 #include "ViennaRNA/params/default.h"
24 #include "ViennaRNA/fold_vars.h"
25 #include "ViennaRNA/plotting/probabilities.h"
26 #include "ViennaRNA/part_func.h"
27 #include "ViennaRNA/params/basic.h"
28 #include "ViennaRNA/loops/all.h"
29 #include "ViennaRNA/LPfold.h"
30 #include "ViennaRNA/Lfold.h"
31 #include "ViennaRNA/alphabet.h"
32 #include "ViennaRNA/part_func_window.h"
33
34 /*
35 #################################
36 # GLOBAL VARIABLES #
37 #################################
38 */
39
40 typedef struct {
41 int bpp_print; /* 1 if pairing probabilities should be written to file-handle, 0 if they are returned as vrna_ep_t */
42 int up_print; /* 1 if unpaired probabilities should be written to file-handle, 0 if they are returned as array */
43
44 FILE *fp_pU;
45 double **pU;
46 FLT_OR_DBL bpp_cutoff;
47 FILE *fp_bpp;
48 vrna_ep_t *bpp;
49 unsigned int bpp_max_size;
50 unsigned int bpp_size;
51 vrna_ep_t *stack_prob;
52 unsigned int stack_prob_size;
53 unsigned int stack_prob_max_size;
54 } default_cb_data;
55
56 typedef struct {
57 FLT_OR_DBL *prml;
58 FLT_OR_DBL *prm_l;
59 FLT_OR_DBL *prm_l1;
60 double **pU;
61 double **pUO;
62 double **pUI;
63 double **pUM;
64 double **pUH;
65 } helper_arrays;
66
67 /* soft constraint contributions function (interior-loops) */
68 typedef FLT_OR_DBL (sc_int)(vrna_fold_compound_t *,
69 int,
70 int,
71 int,
72 int);
73
74 /* QI5 contribution function for unpaired probability computations */
75 typedef void (add_QI5)(FLT_OR_DBL **,
76 int,
77 int,
78 FLT_OR_DBL,
79 FLT_OR_DBL);
80
81 /*
82 #################################
83 # PRIVATE VARIABLES #
84 #################################
85 */
86
87 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
88
89 #ifdef _OPENMP
90 #include <omp.h>
91 #endif
92
93 /* some backward compatibility stuff */
94 PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
95 PRIVATE int backward_compat = 0;
96
97 #ifdef _OPENMP
98
99 #pragma omp threadprivate(backward_compat_compound, backward_compat)
100
101 #endif
102
103 #endif
104 /*
105 #################################
106 # PRIVATE FUNCTION DECLARATIONS #
107 #################################
108 */
109
110 PRIVATE void
111 alloc_helper_arrays(vrna_fold_compound_t *vc,
112 int ulength,
113 helper_arrays *aux_arrays,
114 unsigned int options);
115
116
117 PRIVATE void
118 free_helper_arrays(vrna_fold_compound_t *vc,
119 int ulength,
120 helper_arrays *aux_arrays,
121 unsigned int options);
122
123
124 PRIVATE void
125 compute_probs(vrna_fold_compound_t *vc,
126 int j,
127 helper_arrays *aux_arrays,
128 int ulength,
129 vrna_probs_window_callback *cb,
130 void *data,
131 unsigned int options,
132 int *ov);
133
134
135 PRIVATE void
136 make_ptypes(vrna_fold_compound_t *vc,
137 int j);
138
139
140 PRIVATE void
141 probability_correction(vrna_fold_compound_t *vc,
142 int i);
143
144
145 #if 0
146 PRIVATE vrna_ep_t *
147 get_deppp(vrna_fold_compound_t *vc,
148 vrna_ep_t *pl,
149 int start);
150
151
152 #endif
153
154 PRIVATE void
155 compute_pU(vrna_fold_compound_t *vc,
156 int k,
157 int ulength,
158 helper_arrays *aux_arrays,
159 vrna_probs_window_callback *cb,
160 void *data,
161 unsigned int options);
162
163
164 PRIVATE FLT_OR_DBL *
165 compute_stack_probabilities(vrna_fold_compound_t *vc,
166 int start);
167
168
169 PRIVATE void
170 return_pU(int size,
171 int i,
172 int max_size,
173 helper_arrays *aux_arrays,
174 vrna_probs_window_callback *cb,
175 void *data,
176 unsigned int options);
177
178
179 PRIVATE void
180 print_bpp_callback(FLT_OR_DBL *pr,
181 int size,
182 int k,
183 void *data);
184
185
186 PRIVATE void
187 store_bpp_callback(FLT_OR_DBL *pr,
188 int size,
189 int k,
190 void *data);
191
192
193 #if 0
194 PRIVATE void
195 store_stack_prob_callback(FLT_OR_DBL *pr,
196 int size,
197 int k,
198 void *data);
199
200
201 #endif
202
203 PRIVATE void
204 print_pU_callback(double *pU,
205 int size,
206 int k,
207 int ulength,
208 unsigned int type,
209 void *data);
210
211
212 PRIVATE void
213 store_pU_callback(double *pU,
214 int size,
215 int k,
216 int ulength,
217 unsigned int type,
218 void *data);
219
220
221 PRIVATE void
222 backward_compat_callback(FLT_OR_DBL *pr,
223 int pr_size,
224 int i,
225 int max,
226 unsigned int type,
227 void *data);
228
229
230 PRIVATE FLT_OR_DBL
231 sc_contribution(vrna_fold_compound_t *vc,
232 int i,
233 int j,
234 int k,
235 int l);
236
237
238 PRIVATE FLT_OR_DBL
239 sc_dummy(vrna_fold_compound_t *vc,
240 int i,
241 int j,
242 int k,
243 int l);
244
245
246 /*
247 #################################
248 # BEGIN OF FUNCTION DEFINITIONS #
249 #################################
250 */
251 PUBLIC vrna_ep_t *
vrna_pfl_fold(const char * sequence,int window_size,int max_bp_span,float cutoff)252 vrna_pfl_fold(const char *sequence,
253 int window_size,
254 int max_bp_span,
255 float cutoff)
256 {
257 default_cb_data data;
258
259 data.fp_pU = NULL;
260 data.pU = NULL;
261 data.bpp_cutoff = (FLT_OR_DBL)cutoff;
262 data.fp_bpp = NULL;
263 data.bpp = NULL;
264 data.bpp_max_size = 0;
265 data.bpp_size = 0;
266 data.stack_prob = NULL;
267 data.stack_prob_max_size = 0;
268 data.stack_prob_size = 0;
269 data.bpp_print = 0;
270 data.up_print = 0;
271
272 vrna_pfl_fold_cb(sequence, window_size, max_bp_span, &backward_compat_callback, (void *)&data);
273
274 /* resize pair probability list to actual size */
275 data.bpp =
276 (vrna_ep_t *)vrna_realloc(data.bpp, sizeof(vrna_ep_t) * (data.bpp_size + 1));
277 data.bpp[data.bpp_size].i = 0;
278 data.bpp[data.bpp_size].j = 0;
279 data.bpp[data.bpp_size].type = VRNA_PLIST_TYPE_BASEPAIR;
280 data.bpp[data.bpp_size].p = 0;
281
282 return data.bpp;
283 }
284
285
286 PUBLIC double **
vrna_pfl_fold_up(const char * sequence,int ulength,int window_size,int max_bp_span)287 vrna_pfl_fold_up(const char *sequence,
288 int ulength,
289 int window_size,
290 int max_bp_span)
291 {
292 unsigned int i;
293 double **pU;
294 default_cb_data data;
295
296 pU = NULL;
297
298 if (sequence) {
299 i = strlen(sequence);
300 pU = (double **)vrna_alloc(sizeof(double *) * (i + 2));
301
302 data.fp_pU = NULL;
303 data.pU = pU;
304 data.bpp_cutoff = 0.;
305 data.fp_bpp = NULL;
306 data.bpp = NULL;
307 data.bpp_max_size = 0;
308 data.bpp_size = 0;
309 data.stack_prob = NULL;
310 data.stack_prob_max_size = 0;
311 data.stack_prob_size = 0;
312 data.bpp_print = 0;
313 data.up_print = 0;
314
315 vrna_pfl_fold_up_cb(sequence,
316 ulength,
317 window_size,
318 max_bp_span,
319 &backward_compat_callback,
320 (void *)&data);
321 }
322
323 return pU;
324 }
325
326
327 PRIVATE void
alloc_helper_arrays(vrna_fold_compound_t * vc,int ulength,helper_arrays * aux_arrays,unsigned int options)328 alloc_helper_arrays(vrna_fold_compound_t *vc,
329 int ulength,
330 helper_arrays *aux_arrays,
331 unsigned int options)
332 {
333 int i, n;
334
335 n = vc->length;
336
337 aux_arrays->pU = NULL;
338 aux_arrays->pUO = NULL;
339 aux_arrays->pUH = NULL;
340 aux_arrays->pUI = NULL;
341 aux_arrays->pUM = NULL;
342
343 aux_arrays->prm_l = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
344 aux_arrays->prm_l1 = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
345 aux_arrays->prml = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n + 2));
346
347 if ((options & VRNA_PROBS_WINDOW_UP) && (ulength > 0)) {
348 aux_arrays->pU = (double **)vrna_alloc((n + 1) * sizeof(double *));
349 for (i = 1; i <= n; i++)
350 aux_arrays->pU[i] = (double *)vrna_alloc((MAX2(MAXLOOP, ulength) + 2) * sizeof(double));
351
352 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
353 aux_arrays->pUO = (double **)vrna_alloc((n + 1) * sizeof(double *));
354 aux_arrays->pUI = (double **)vrna_alloc((n + 1) * sizeof(double *));
355 aux_arrays->pUM = (double **)vrna_alloc((n + 1) * sizeof(double *));
356 aux_arrays->pUH = (double **)vrna_alloc((n + 1) * sizeof(double *));
357 for (i = 1; i <= n; i++) {
358 aux_arrays->pUH[i] = (double *)vrna_alloc((MAX2(MAXLOOP, ulength) + 2) * sizeof(double));
359 aux_arrays->pUI[i] = (double *)vrna_alloc((MAX2(MAXLOOP, ulength) + 2) * sizeof(double));
360 aux_arrays->pUO[i] = (double *)vrna_alloc((MAX2(MAXLOOP, ulength) + 2) * sizeof(double));
361 aux_arrays->pUM[i] = (double *)vrna_alloc((MAX2(MAXLOOP, ulength) + 2) * sizeof(double));
362 }
363 }
364 }
365 }
366
367
368 PRIVATE void
free_helper_arrays(vrna_fold_compound_t * vc,int ulength,helper_arrays * aux_arrays,unsigned int options)369 free_helper_arrays(vrna_fold_compound_t *vc,
370 int ulength,
371 helper_arrays *aux_arrays,
372 unsigned int options)
373 {
374 int i, n;
375
376 n = vc->length;
377
378 free(aux_arrays->prm_l);
379 free(aux_arrays->prm_l1);
380 free(aux_arrays->prml);
381
382 if ((options & VRNA_PROBS_WINDOW_UP) && (ulength > 0)) {
383 for (i = 1; i <= n; i++)
384 free(aux_arrays->pU[i]);
385 free(aux_arrays->pU);
386
387 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
388 for (i = 1; i <= n; i++) {
389 free(aux_arrays->pUH[i]);
390 free(aux_arrays->pUI[i]);
391 free(aux_arrays->pUO[i]);
392 free(aux_arrays->pUM[i]);
393 }
394 free(aux_arrays->pUH);
395 free(aux_arrays->pUI);
396 free(aux_arrays->pUO);
397 free(aux_arrays->pUM);
398 }
399 }
400 }
401
402
403 PRIVATE void
return_pU(int size,int i,int max_size,helper_arrays * aux_arrays,vrna_probs_window_callback * cb,void * data,unsigned int options)404 return_pU(int size,
405 int i,
406 int max_size,
407 helper_arrays *aux_arrays,
408 vrna_probs_window_callback *cb,
409 void *data,
410 unsigned int options)
411 {
412 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
413 cb(aux_arrays->pUO[i], size, i, max_size, VRNA_PROBS_WINDOW_UP | VRNA_EXT_LOOP, data);
414 cb(aux_arrays->pUH[i], size, i, max_size, VRNA_PROBS_WINDOW_UP | VRNA_HP_LOOP, data);
415 cb(aux_arrays->pUI[i], size, i, max_size, VRNA_PROBS_WINDOW_UP | VRNA_INT_LOOP, data);
416 cb(aux_arrays->pUM[i], size, i, max_size, VRNA_PROBS_WINDOW_UP | VRNA_MB_LOOP, data);
417 } else {
418 cb(aux_arrays->pU[i], size, i, max_size, VRNA_PROBS_WINDOW_UP | VRNA_ANY_LOOP, data);
419 }
420 }
421
422
423 PRIVATE INLINE void
allocate_dp_matrices(vrna_fold_compound_t * vc,int i,unsigned int options)424 allocate_dp_matrices(vrna_fold_compound_t *vc,
425 int i,
426 unsigned int options)
427 {
428 char **ptype;
429 int winSize;
430 FLT_OR_DBL **pR, **q, **qb, **qm, **qm2, **QI5, **qmb, **q2l;
431 vrna_mx_pf_t *mx;
432 vrna_hc_t *hc;
433 vrna_sc_t *sc;
434
435 mx = vc->exp_matrices;
436 pR = mx->pR;
437 q = mx->q_local;
438 qb = mx->qb_local;
439 qm = mx->qm_local;
440 qm2 = mx->qm2_local;
441 QI5 = mx->QI5;
442 qmb = mx->qmb;
443 q2l = mx->q2l;
444 ptype = vc->ptype_local;
445 winSize = vc->window_size;
446 hc = vc->hc;
447
448 /* allocate new part of arrays */
449 pR[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
450 pR[i] -= i;
451 q[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
452 q[i] -= i;
453 qb[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
454 qb[i] -= i;
455 qm[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
456 qm[i] -= i;
457 if (options & VRNA_PROBS_WINDOW_UP) {
458 qm2[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
459 qm2[i] -= i;
460 QI5[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
461 qmb[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
462 q2l[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
463 }
464
465 hc->matrix_local[i] = (unsigned char *)vrna_alloc((winSize + 1) * sizeof(unsigned char));
466 ptype[i] = (char *)vrna_alloc((winSize + 1) * sizeof(char));
467 ptype[i] -= i;
468
469 switch (vc->type) {
470 case VRNA_FC_TYPE_SINGLE:
471 sc = vc->sc;
472 if (sc) {
473 if (sc->exp_energy_bp_local)
474 sc->exp_energy_bp_local[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
475
476 if (sc->exp_energy_up)
477 sc->exp_energy_up[i] = (FLT_OR_DBL *)vrna_alloc((winSize + 1) * sizeof(FLT_OR_DBL));
478
479 vrna_sc_update(vc, i, VRNA_OPTION_PF | VRNA_OPTION_WINDOW);
480 }
481
482 break;
483
484 case VRNA_FC_TYPE_COMPARATIVE:
485 break;
486 }
487 }
488
489
490 PRIVATE INLINE void
free_dp_matrices(vrna_fold_compound_t * vc,unsigned int options)491 free_dp_matrices(vrna_fold_compound_t *vc,
492 unsigned int options)
493 {
494 char **ptype;
495 int i, n, winSize;
496 FLT_OR_DBL **pR, **q, **qb, **qm, **qm2, **QI5, **qmb, **q2l;
497 vrna_mx_pf_t *mx;
498 vrna_hc_t *hc;
499 vrna_sc_t *sc;
500
501 n = (int)vc->length;
502 winSize = vc->window_size;
503 mx = vc->exp_matrices;
504 pR = mx->pR;
505 q = mx->q_local;
506 qb = mx->qb_local;
507 qm = mx->qm_local;
508 ptype = vc->ptype_local;
509 hc = vc->hc;
510 sc = vc->sc;
511
512 for (i = MAX2(1, n - (winSize + MAXLOOP)); i <= n; i++) {
513 free(pR[i] + i);
514 free(q[i] + i);
515 free(qb[i] + i);
516 free(qm[i] + i);
517 pR[i] = NULL;
518 q[i] = NULL;
519 qb[i] = NULL;
520 qm[i] = NULL;
521 if (options & VRNA_PROBS_WINDOW_UP) {
522 qm2 = mx->qm2_local;
523 QI5 = mx->QI5;
524 qmb = mx->qmb;
525 q2l = mx->q2l;
526 free(qm2[i] + i);
527 free(QI5[i]);
528 free(qmb[i]);
529 free(q2l[i]);
530 qm2[i] = NULL;
531 QI5[i] = NULL;
532 qmb[i] = NULL;
533 q2l[i] = NULL;
534 }
535
536 free(hc->matrix_local[i]);
537 hc->matrix_local[i] = NULL;
538 free(ptype[i] + i);
539 ptype[i] = NULL;
540
541 if (sc) {
542 if (sc->exp_energy_up)
543 free(sc->exp_energy_up[i]);
544
545 if (sc->exp_energy_bp_local)
546 free(sc->exp_energy_bp_local[i]);
547 }
548 }
549 }
550
551
552 PRIVATE INLINE void
init_dp_matrices(vrna_fold_compound_t * vc,unsigned int options)553 init_dp_matrices(vrna_fold_compound_t *vc,
554 unsigned int options)
555 {
556 int j, max_j, winSize;
557
558 winSize = vc->window_size;
559 max_j = MIN2((int)vc->length, 2 * winSize + MAXLOOP + 2);
560
561 for (j = 1; j <= max_j; j++)
562 allocate_dp_matrices(vc, j, options);
563 }
564
565
566 PRIVATE INLINE void
rotate_dp_matrices(vrna_fold_compound_t * vc,int j,unsigned int options)567 rotate_dp_matrices(vrna_fold_compound_t *vc,
568 int j,
569 unsigned int options)
570 {
571 char **ptype;
572 int i, winSize, length;
573 FLT_OR_DBL **pR, **q, **qb, **qm, **qm2, **QI5, **qmb, **q2l;
574 vrna_mx_pf_t *mx;
575 vrna_hc_t *hc;
576 vrna_sc_t *sc;
577
578 length = vc->length;
579 winSize = vc->window_size;
580 mx = vc->exp_matrices;
581 pR = mx->pR;
582 q = mx->q_local;
583 qb = mx->qb_local;
584 qm = mx->qm_local;
585 ptype = vc->ptype_local;
586 hc = vc->hc;
587 sc = vc->sc;
588
589 if (j > 2 * winSize + MAXLOOP + 1) {
590 i = j - (2 * winSize + MAXLOOP + 1);
591 /* free arrays may be faster than pointer rotation and reset to 0 values */
592 free(pR[i] + i);
593 free(q[i] + i);
594 free(qb[i] + i);
595 free(qm[i] + i);
596 pR[i] = NULL;
597 q[i] = NULL;
598 qb[i] = NULL;
599 qm[i] = NULL;
600 if (options & VRNA_PROBS_WINDOW_UP) {
601 qm2 = mx->qm2_local;
602 QI5 = mx->QI5;
603 qmb = mx->qmb;
604 q2l = mx->q2l;
605 free(qm2[i] + i);
606 free(QI5[i]);
607 free(qmb[i]);
608 free(q2l[i]);
609 qm2[i] = NULL;
610 QI5[i] = NULL;
611 qmb[i] = NULL;
612 q2l[i] = NULL;
613 }
614
615 free(hc->matrix_local[i]);
616 hc->matrix_local[i] = NULL;
617 free(ptype[i] + i);
618 ptype[i] = NULL;
619
620 if (sc) {
621 if (sc->exp_energy_up) {
622 free(sc->exp_energy_up[i]);
623 sc->exp_energy_up[i] = NULL;
624 }
625
626 if (sc->exp_energy_bp_local) {
627 free(sc->exp_energy_bp_local[i]);
628 sc->exp_energy_bp_local[i] = NULL;
629 }
630 }
631
632 if (j + 1 <= length)
633 /* get arrays for next round */
634 allocate_dp_matrices(vc, j + 1, options);
635 }
636 }
637
638
639 PRIVATE INLINE void
init_constraints(vrna_fold_compound_t * fc,unsigned int options)640 init_constraints(vrna_fold_compound_t *fc,
641 unsigned int options)
642 {
643 int j, max_j, winSize;
644
645 winSize = fc->window_size;
646 max_j = MIN2((int)fc->length, 2 * winSize + MAXLOOP + 2);
647
648 for (j = 1; j <= max_j; j++) {
649 make_ptypes(fc, j);
650 vrna_hc_update(fc, j, VRNA_CONSTRAINT_WINDOW_UPDATE_5);
651 vrna_sc_update(fc, j, VRNA_OPTION_PF | VRNA_OPTION_WINDOW);
652 }
653 }
654
655
656 PRIVATE INLINE void
rotate_constraints(vrna_fold_compound_t * fc,int j,unsigned int options)657 rotate_constraints(vrna_fold_compound_t *fc,
658 int j,
659 unsigned int options)
660 {
661 if (j + 1 <= fc->length) {
662 make_ptypes(fc, j + 1);
663 vrna_hc_update(fc, j + 1, VRNA_CONSTRAINT_WINDOW_UPDATE_5);
664 vrna_sc_update(fc, j + 1, VRNA_OPTION_PF | VRNA_OPTION_WINDOW);
665 }
666 }
667
668
669 PUBLIC int
vrna_probs_window(vrna_fold_compound_t * vc,int ulength,unsigned int options,vrna_probs_window_callback * cb,void * data)670 vrna_probs_window(vrna_fold_compound_t *vc,
671 int ulength,
672 unsigned int options,
673 vrna_probs_window_callback *cb,
674 void *data)
675 {
676 unsigned char hc_decompose;
677 int n, i, j, k, maxl, ov, winSize, pairSize, turn;
678 FLT_OR_DBL temp, Qmax, qbt1, **q, **qb, **qm, **qm2, **pR;
679 double max_real, *Fwindow;
680 vrna_exp_param_t *pf_params;
681 vrna_md_t *md;
682 vrna_mx_pf_t *matrices;
683 vrna_hc_t *hc;
684 helper_arrays aux_arrays;
685 vrna_mx_pf_aux_el_t aux_mx_el;
686 vrna_mx_pf_aux_ml_t aux_mx_ml;
687
688 ov = 0;
689 Qmax = 0;
690
691 if ((!vc) || (!cb))
692 return 0; /* failure */
693
694 if (!vrna_fold_compound_prepare(vc, VRNA_OPTION_PF | VRNA_OPTION_WINDOW)) {
695 vrna_message_warning("vrna_probs_window: "
696 "Failed to prepare vrna_fold_compound");
697 return 0; /* failure */
698 }
699
700 /* here space for initializing everything */
701
702 n = vc->length;
703 pf_params = vc->exp_params;
704 md = &(pf_params->model_details);
705 matrices = vc->exp_matrices;
706 winSize = vc->window_size;
707 pairSize = md->max_bp_span;
708 turn = md->min_loop_size;
709
710 q = matrices->q_local;
711 qb = matrices->qb_local;
712 qm = matrices->qm_local;
713 qm2 = matrices->qm2_local;
714 pR = matrices->pR;
715
716 hc = vc->hc;
717
718 alloc_helper_arrays(vc, ulength, &aux_arrays, options);
719
720 Fwindow =
721 (options & VRNA_PROBS_WINDOW_PF) ? (double *)vrna_alloc(sizeof(double) * (winSize + 1)) : NULL;
722
723 /* very short molecule ? */
724 if (n < turn + 2) {
725 if ((options & VRNA_PROBS_WINDOW_UP) && (ulength > 0)) {
726 for (i = 1; i <= n; i++) {
727 maxl = MIN2(MAX2(MAXLOOP, ulength), n);
728 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
729 for (j = 0; j <= maxl; j++) {
730 aux_arrays.pUO[i][j] = 1.;
731 aux_arrays.pUH[i][j] = 0.;
732 aux_arrays.pUI[i][j] = 0.;
733 aux_arrays.pUM[i][j] = 0.;
734 }
735 } else {
736 for (j = 0; j <= maxl; j++)
737 aux_arrays.pU[i][j] = 1.;
738 }
739
740 return_pU(maxl, i, ulength, &aux_arrays, cb, data, options);
741 }
742 }
743
744 free_helper_arrays(vc, ulength, &aux_arrays, options);
745
746 return 1; /* success */
747 }
748
749 init_dp_matrices(vc, options);
750 init_constraints(vc, options);
751
752 /* init auxiliary arrays for fast exterior/multibranch loops */
753 aux_mx_el = vrna_exp_E_ext_fast_init(vc);
754 aux_mx_ml = vrna_exp_E_ml_fast_init(vc);
755
756 max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
757
758 /* start recursions */
759 for (j = turn + 2; j <= n + winSize; j++) {
760 if (j <= n) {
761 vrna_exp_E_ext_fast_update(vc, j, aux_mx_el);
762 for (i = j - turn - 1; i >= MAX2(1, (j - winSize + 1)); i--) {
763 hc_decompose = hc->matrix_local[i][j - i];
764 qbt1 = 0.;
765
766 /*
767 * construction of partition function of segment i,j
768 * firstly that given i bound to j : qb(i,j)
769 */
770 if (hc_decompose) {
771 /* process hairpin loop(s) */
772 qbt1 += vrna_exp_E_hp_loop(vc, i, j);
773 /* process interior loop(s) */
774 qbt1 += vrna_exp_E_int_loop(vc, i, j);
775 /* process multibranch loop(s) */
776 qbt1 += vrna_exp_E_mb_loop_fast(vc, i, j, aux_mx_ml);
777 }
778
779 qb[i][j] = qbt1;
780
781 /* Multibranch loop */
782 qm[i][j] = vrna_exp_E_ml_fast(vc, i, j, aux_mx_ml);
783 if ((options & VRNA_PROBS_WINDOW_UP) && (ulength > 0)) {
784 /* new qm2 computation done here */
785 const FLT_OR_DBL *qqm = vrna_exp_E_ml_fast_qqm(aux_mx_ml);
786 temp = 0.0;
787 for (k = i + 1; k <= j; k++)
788 temp += qm[i][k - 1] *
789 qqm[k];
790 qm2[i][j] = temp;
791 }
792
793 /* Exterior loop */
794 q[i][j] = temp = vrna_exp_E_ext_fast(vc, i, j, aux_mx_el);
795
796 if (temp > Qmax) {
797 Qmax = temp;
798 if (Qmax > max_real / 10.) {
799 vrna_message_warning("vrna_probs_window: "
800 "Q close to overflow: %d %d %g\n",
801 i,
802 j,
803 temp);
804 }
805 }
806
807 if (temp >= max_real) {
808 vrna_message_warning("vrna_probs_window: "
809 "overflow while computing partition function for segment q[%d,%d]\n"
810 "use larger pf_scale",
811 i,
812 j);
813
814 vrna_exp_E_ml_fast_free(aux_mx_ml);
815 vrna_exp_E_ext_fast_free(aux_mx_el);
816 free_helper_arrays(vc, ulength, &aux_arrays, options);
817
818 return 0; /* failure */
819 }
820 } /* end for i */
821
822 /*
823 * here we return the partition function for subsegments [i...j] in terms
824 * of ensemble free energies G_ij = -RT * ln(Q_ij) in kcal/mol
825 */
826 if (options & VRNA_PROBS_WINDOW_PF) {
827 int start = MAX2(1, j - winSize + 1);
828 Fwindow -= start;
829 for (i = start; i <= j; i++)
830 Fwindow[i] = (double)(-log(q[i][j]) - (j - i + 1) * log(pf_params->pf_scale)) *
831 pf_params->kT / 1000.0;
832
833 cb(Fwindow, j, start, winSize, VRNA_PROBS_WINDOW_PF, data);
834 Fwindow += start;
835 }
836
837 /*
838 * just as a general service, I save here the free energy of the windows
839 * no output is generated, however,...
840 */
841 if ((j >= winSize) && (options & VRNA_PROBS_WINDOW_UP)) {
842 FLT_OR_DBL eee = 0.;
843 eee = (FLT_OR_DBL)(-log(q[j - winSize + 1][j]) - winSize * log(pf_params->pf_scale)) *
844 pf_params->kT / 1000.0;
845
846 /* we could return this to the user via callback cb() if we were nice */
847
848 aux_arrays.pU[j][0] = eee;
849 }
850
851 /* rotate auxiliary arrays */
852 vrna_exp_E_ext_fast_rotate(aux_mx_el);
853 vrna_exp_E_ml_fast_rotate(aux_mx_ml);
854 }
855
856 if (j > winSize) {
857 compute_probs(vc, j, &aux_arrays, ulength, cb, data, options, &ov);
858
859 if ((options & VRNA_PROBS_WINDOW_UP) && (j > winSize + MAXLOOP + 1))
860 compute_pU(vc, j - winSize - MAXLOOP - 1, ulength, &aux_arrays, cb, data, options);
861
862 if (j > 2 * winSize + MAXLOOP + 1) {
863 int start = j - (2 * winSize + MAXLOOP + 1);
864 probability_correction(vc, start);
865 if (options & VRNA_PROBS_WINDOW_BPP) {
866 cb(pR[start],
867 MIN2(start + winSize, n),
868 start,
869 winSize,
870 VRNA_PROBS_WINDOW_BPP,
871 data);
872 }
873
874 if (options & VRNA_PROBS_WINDOW_STACKP) {
875 int start = j - (2 * winSize - MAXLOOP);
876 if (start > 1) {
877 FLT_OR_DBL *stack_probs = compute_stack_probabilities(vc, start);
878 stack_probs -= start + 1;
879 cb(stack_probs,
880 MIN2(n - start + turn, pairSize),
881 start,
882 winSize,
883 VRNA_PROBS_WINDOW_STACKP,
884 data);
885 stack_probs += start + 1;
886 free(stack_probs);
887 }
888 }
889
890 rotate_dp_matrices(vc, j, options);
891 rotate_constraints(vc, j, options);
892 }
893 } /* end if (do_backtrack) */
894 } /* end for j */
895
896 /* finish output */
897 if (options & VRNA_PROBS_WINDOW_UP)
898 for (j = MAX2(1, n - MAXLOOP); j <= n; j++)
899 compute_pU(vc, j, ulength, &aux_arrays, cb, data, options);
900
901 for (j = MAX2(n - winSize - MAXLOOP, 1); j <= n; j++) {
902 probability_correction(vc, j);
903 if (options & VRNA_PROBS_WINDOW_BPP) {
904 cb(pR[j],
905 MIN2(j + winSize, n),
906 j,
907 winSize,
908 VRNA_PROBS_WINDOW_BPP,
909 data);
910 }
911
912 if ((options & VRNA_PROBS_WINDOW_STACKP) && j < n) {
913 int start = j;
914 if (start > 1) {
915 FLT_OR_DBL *stack_probs = compute_stack_probabilities(vc, start);
916 stack_probs -= start + 1;
917 cb(stack_probs,
918 MIN2(n - start + turn, pairSize),
919 start,
920 winSize,
921 VRNA_PROBS_WINDOW_STACKP,
922 data);
923 stack_probs += start + 1;
924 free(stack_probs);
925 }
926 }
927 }
928
929 if (ov > 0) {
930 vrna_message_warning("vrna_probs_window: "
931 "%d overflows occurred while backtracking;\n"
932 "you might try a smaller pf_scale than %g\n",
933 ov,
934 pf_params->pf_scale);
935 }
936
937 free_dp_matrices(vc, options);
938 free_helper_arrays(vc, ulength, &aux_arrays, options);
939
940 /* free memory occupied by auxiliary arrays for fast exterior/multibranch loops */
941 vrna_exp_E_ml_fast_free(aux_mx_ml);
942 vrna_exp_E_ext_fast_free(aux_mx_el);
943
944 free(Fwindow);
945
946 return 1; /* success */
947 }
948
949
950 PRIVATE FLT_OR_DBL
sc_contribution(vrna_fold_compound_t * vc,int i,int j,int k,int l)951 sc_contribution(vrna_fold_compound_t *vc,
952 int i,
953 int j,
954 int k,
955 int l)
956 {
957 FLT_OR_DBL q;
958 vrna_sc_t *sc;
959
960 q = 1.;
961 sc = vc->sc;
962
963 if (sc->exp_energy_up)
964 q *= sc->exp_energy_up[i + 1][k - i - 1] *
965 sc->exp_energy_up[l + 1][j - l - 1];
966
967 if (sc->exp_energy_bp_local)
968 q *= sc->exp_energy_bp_local[i][j - i];
969
970 if ((sc->exp_energy_stack) && (i + 1 == k) && (l + 1 == j)) {
971 q *= sc->exp_energy_stack[i] *
972 sc->exp_energy_stack[k] *
973 sc->exp_energy_stack[l] *
974 sc->exp_energy_stack[j];
975 }
976
977 if (sc->f)
978 q *= sc->f(i, j, k, l, VRNA_DECOMP_PAIR_IL, sc->data);
979
980 return q;
981 }
982
983
984 PRIVATE FLT_OR_DBL
sc_dummy(vrna_fold_compound_t * vc,int i,int j,int k,int l)985 sc_dummy(vrna_fold_compound_t *vc,
986 int i,
987 int j,
988 int k,
989 int l)
990 {
991 return 1.;
992 }
993
994
995 PRIVATE void
add_QI5_contribution(FLT_OR_DBL ** QI5,int i,int j,FLT_OR_DBL q,FLT_OR_DBL qkl)996 add_QI5_contribution(FLT_OR_DBL **QI5,
997 int i,
998 int j,
999 FLT_OR_DBL q,
1000 FLT_OR_DBL qkl)
1001 {
1002 QI5[i][j] += q * qkl;
1003 }
1004
1005
1006 PRIVATE void
add_QI5_dummy(FLT_OR_DBL ** QI5,int i,int j,FLT_OR_DBL q,FLT_OR_DBL qkl)1007 add_QI5_dummy(FLT_OR_DBL **QI5,
1008 int i,
1009 int j,
1010 FLT_OR_DBL q,
1011 FLT_OR_DBL qkl)
1012 {
1013 return;
1014 }
1015
1016
1017 PRIVATE void
compute_probs(vrna_fold_compound_t * vc,int j,helper_arrays * aux_arrays,int ulength,vrna_probs_window_callback * cb,void * data,unsigned int options,int * ov)1018 compute_probs(vrna_fold_compound_t *vc,
1019 int j,
1020 helper_arrays *aux_arrays,
1021 int ulength,
1022 vrna_probs_window_callback *cb,
1023 void *data,
1024 unsigned int options,
1025 int *ov)
1026 {
1027 char **ptype;
1028 short *S1;
1029 int start_i, i, k, l, n, m, winSize, turn, type, type_2, tt, *rtype;
1030 FLT_OR_DBL *prml, *prm_l, *prm_l1, **pR, **QI5, **qmb, **q2l, **qb, **q, **qm,
1031 *scale, *expMLbase, expMLclosing, temp, prm_MLb, prmt1, prmt, *tmp,
1032 Qmax;
1033 double max_real;
1034 vrna_exp_param_t *pf_params;
1035 vrna_md_t *md;
1036 vrna_hc_t *hc;
1037 vrna_sc_t *sc;
1038 sc_int *sc_int_f;
1039 add_QI5 *add_QI5_f;
1040
1041 max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
1042
1043 prml = aux_arrays->prml;
1044 prm_l = aux_arrays->prm_l;
1045 prm_l1 = aux_arrays->prm_l1;
1046
1047 n = vc->length;
1048 winSize = vc->window_size;
1049 S1 = vc->sequence_encoding;
1050 ptype = vc->ptype_local;
1051 pf_params = vc->exp_params;
1052 md = &(pf_params->model_details);
1053 turn = md->min_loop_size;
1054 rtype = &(md->rtype[0]);
1055 expMLclosing = pf_params->expMLclosing;
1056 scale = vc->exp_matrices->scale;
1057 expMLbase = vc->exp_matrices->expMLbase;
1058 hc = vc->hc;
1059 sc = vc->sc;
1060
1061 pR = vc->exp_matrices->pR;
1062 QI5 = vc->exp_matrices->QI5;
1063 qmb = vc->exp_matrices->qmb;
1064 q2l = vc->exp_matrices->q2l;
1065 q = vc->exp_matrices->q_local;
1066 qb = vc->exp_matrices->qb_local;
1067 qm = vc->exp_matrices->qm_local;
1068
1069 Qmax = 0;
1070
1071 /* assign helper functions */
1072 if (sc)
1073 sc_int_f = &sc_contribution;
1074 else
1075 sc_int_f = &sc_dummy;
1076
1077 if (options & VRNA_PROBS_WINDOW_UP)
1078 add_QI5_f = &add_QI5_contribution;
1079 else
1080 add_QI5_f = &add_QI5_dummy;
1081
1082 /* start recursion */
1083
1084 /*
1085 * i=j-winSize;
1086 * initialize multiloopfs
1087 */
1088 for (k = j - winSize; k <= MIN2(n, j); k++) {
1089 prml[k] = 0;
1090 prm_l[k] = 0;
1091 /* prm_l1[k]=0; others stay */
1092 }
1093 k = j - winSize;
1094 prm_l1[k] = 0;
1095
1096 for (l = k + turn + 1; l <= MIN2(n, k + winSize - 1); l++) {
1097 int a;
1098 pR[k][l] = 0; /* set zero at start */
1099 type = vrna_get_ptype_window(k, l + k, ptype);
1100
1101 if (qb[k][l] == 0)
1102 continue;
1103
1104 /* Exterior loop cases */
1105 if (hc->matrix_local[k][l - k] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
1106 for (a = MAX2(1, l - winSize + 2); a < MIN2(k, n - winSize + 2); a++)
1107 pR[k][l] += q[a][k - 1] *
1108 q[l + 1][a + winSize - 1] /
1109 q[a][a + winSize - 1];
1110
1111 if (l - k + 1 == winSize) {
1112 pR[k][l] += 1. / q[k][l];
1113 } else {
1114 if (k + winSize - 1 <= n) /* k outermost */
1115 pR[k][l] += q[l + 1][k + winSize - 1] /
1116 q[k][k + winSize - 1];
1117
1118 if (l - winSize + 1 >= 1) /* l outermost */
1119 pR[k][l] += q[l - winSize + 1][k - 1] /
1120 q[l - winSize + 1][l];
1121 }
1122
1123 pR[k][l] *= vrna_exp_E_ext_stem(type,
1124 (k > 1) ? S1[k - 1] : -1,
1125 (l < n) ? S1[l + 1] : -1,
1126 pf_params);
1127 }
1128
1129 if (hc->matrix_local[k][l - k] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC) {
1130 FLT_OR_DBL ppp;
1131
1132 type_2 = rtype[vrna_get_ptype_window(k, l + k, ptype)];
1133 ppp = 0.;
1134 start_i = k - MAXLOOP - 1;
1135
1136 if (start_i < l - winSize + 1)
1137 start_i = l - winSize + 1;
1138
1139 if (start_i < 1)
1140 start_i = 1;
1141
1142 int u1 = 0;
1143 short sk1, sl1, si1;
1144
1145 sk1 = S1[k - 1];
1146 sl1 = S1[l + 1];
1147 for (i = k - 1; i >= start_i; i--, u1++) {
1148 int max_m = i + winSize - 1;
1149
1150 if (hc->up_int[i + 1] < u1)
1151 break;
1152
1153 si1 = S1[i + 1];
1154
1155 if (max_m > l + MAXLOOP - u1 + 1)
1156 max_m = l + MAXLOOP - u1 + 1;
1157
1158 if (max_m > n)
1159 max_m = n;
1160
1161 for (m = l + 1; m <= max_m; m++) {
1162 int u2 = m - l - 1;
1163
1164 if (hc->up_int[l + 1] < u2)
1165 break;
1166
1167 if (hc->matrix_local[i][m - i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
1168 type = vrna_get_ptype_window(i, m + i, ptype);
1169 if (pR[i][m] > 0) {
1170 temp = pR[i][m] *
1171 exp_E_IntLoop(u1,
1172 u2,
1173 type,
1174 type_2,
1175 si1,
1176 S1[m - 1],
1177 sk1,
1178 sl1,
1179 pf_params) *
1180 sc_int_f(vc, i, m, k, l) *
1181 scale[u1 + u2 + 2];
1182
1183 add_QI5_f(QI5, i, k - i - 1, temp, qb[k][l]);
1184 add_QI5_f(QI5, l, m - l - 1, temp, qb[k][l]);
1185
1186 ppp += temp;
1187 }
1188 }
1189 }
1190 }
1191
1192 pR[k][l] += ppp;
1193 }
1194 }
1195
1196 /* 3. bonding k,l as substem of multi-loop enclosed by i,m */
1197 prm_MLb = 0.;
1198 if (k > 1) {
1199 /* sonst nix! */
1200 for (l = MIN2(n - 1, k + winSize - 2); l >= k + turn + 1; l--) {
1201 FLT_OR_DBL ppp;
1202
1203 /* opposite direction */
1204 m = l + 1;
1205 prmt = prmt1 = 0.0;
1206
1207 for (i = MAX2(1, l - winSize + 2); i < k - 1 /* turn */; i++) {
1208 if (hc->matrix_local[i][m - i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1209 tt = rtype[vrna_get_ptype_window(i, m + i, ptype)];
1210 ppp = pR[i][m] *
1211 exp_E_MLstem(tt, S1[m - 1], S1[i + 1], pf_params) *
1212 qm[i + 1][k - 1];
1213
1214 if (sc)
1215 if (sc->exp_energy_bp_local)
1216 ppp *= sc->exp_energy_bp_local[i][m - i];
1217
1218 prmt += ppp;
1219 }
1220 }
1221 prmt *= expMLclosing;
1222
1223 prml[m] = prmt;
1224
1225 if (hc->matrix_local[k - 1][m - k + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1226 tt = rtype[vrna_get_ptype_window(k - 1, m + k - 1, ptype)];
1227 prmt1 = pR[k - 1][m] *
1228 expMLclosing *
1229 exp_E_MLstem(tt,
1230 S1[l],
1231 S1[k],
1232 pf_params);
1233
1234
1235 if (sc)
1236 if (sc->exp_energy_bp_local)
1237 prmt1 *= sc->exp_energy_bp_local[k - 1][m - k + 1];
1238 }
1239
1240 /* k-1 is unpaired */
1241 if (hc->up_ml[k - 1]) {
1242 ppp = prm_l1[m] * expMLbase[1];
1243
1244 if (sc)
1245 if (sc->exp_energy_up)
1246 ppp *= sc->exp_energy_up[k - 1][1];
1247
1248 prm_l[m] = ppp + prmt1;
1249 } else {
1250 /* skip configuration where k-1 is unpaired */
1251 prm_l[m] = prmt1;
1252 }
1253
1254 /* m is unpaired */
1255 if (hc->up_ml[m]) {
1256 ppp = prm_MLb * expMLbase[1];
1257
1258 if (sc)
1259 if (sc->exp_energy_up)
1260 ppp *= sc->exp_energy_up[m][1];
1261
1262 prm_MLb = ppp + prml[m];
1263 } else {
1264 prm_MLb = prml[m];
1265 }
1266
1267 /*
1268 * same as: prm_MLb = 0;
1269 * for (i=n; i>k; i--) prm_MLb += prml[i]*expMLbase[k-i-1];
1270 */
1271 prml[m] = prml[m] + prm_l[m];
1272
1273 if (qb[k][l] == 0.)
1274 continue;
1275
1276 if (hc->matrix_local[k][l - k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) {
1277 tt = vrna_get_ptype_window(k, l + k, ptype);
1278
1279 if (options & VRNA_PROBS_WINDOW_UP) {
1280 double dang;
1281 /* coefficient for computations of unpairedarrays */
1282 dang = qb[k][l] *
1283 exp_E_MLstem(tt,
1284 (k > 1) ? S1[k - 1] : -1,
1285 (l < n) ? S1[l + 1] : -1,
1286 pf_params) *
1287 scale[2];
1288
1289 for (m = MIN2(k + winSize - 2, n); m >= l + 2; m--) {
1290 qmb[l][m - l - 1] += prml[m] * dang;
1291 q2l[l][m - l - 1] += (prml[m] - prm_l[m]) * dang;
1292 }
1293 }
1294
1295 temp = prm_MLb;
1296
1297 for (m = MIN2(k + winSize - 2, n); m >= l + 2; m--)
1298 temp += prml[m] * qm[l + 1][m - 1];
1299
1300
1301 temp *= exp_E_MLstem(tt,
1302 (k > 1) ? S1[k - 1] : -1,
1303 (l < n) ? S1[l + 1] : -1,
1304 pf_params) * scale[2];
1305 pR[k][l] += temp;
1306 }
1307
1308 if (pR[k][l] > Qmax) {
1309 Qmax = pR[k][l];
1310 if (Qmax > max_real / 10.)
1311 vrna_message_warning("P close to overflow: %d %d %g %g\n",
1312 i, m, pR[k][l], qb[k][l]);
1313 }
1314
1315 if (pR[k][l] >= max_real) {
1316 (*ov)++;
1317 pR[k][l] = FLT_MAX;
1318 }
1319 } /* end for (l=..) */
1320 }
1321
1322 tmp = prm_l1;
1323 aux_arrays->prm_l1 = prm_l;
1324 aux_arrays->prm_l = tmp;
1325 }
1326
1327
1328 PRIVATE void
probability_correction(vrna_fold_compound_t * vc,int i)1329 probability_correction(vrna_fold_compound_t *vc,
1330 int i)
1331 {
1332 int j, howoften, pairdist, turn, n, winSize;
1333 FLT_OR_DBL **qb, **pR;
1334
1335 n = vc->length;
1336 winSize = vc->window_size;
1337 turn = vc->exp_params->model_details.min_loop_size;
1338 howoften = 0; /* how many samples do we have for this pair */
1339
1340 qb = vc->exp_matrices->qb_local;
1341 pR = vc->exp_matrices->pR;
1342
1343 for (j = i + turn; j < MIN2(i + winSize, n + 1); j++) {
1344 pairdist = (j - i + 1);
1345 /* 4cases */
1346 howoften = MIN2(winSize - pairdist + 1, i); /* pairdist,start */
1347 howoften = MIN2(howoften, n - j + 1); /* end */
1348 howoften = MIN2(howoften, n - winSize + 1); /* windowsize */
1349 pR[i][j] *= qb[i][j] / howoften;
1350 }
1351 return;
1352 }
1353
1354
1355 PRIVATE void
make_ptypes(vrna_fold_compound_t * vc,int i)1356 make_ptypes(vrna_fold_compound_t *vc,
1357 int i)
1358 {
1359 /* make new entries in ptype array */
1360 char **ptype;
1361 const short *S;
1362 int j, type, pairSize, n;
1363 vrna_md_t *md;
1364
1365 ptype = vc->ptype_local;
1366 md = &(vc->exp_params->model_details);
1367 pairSize = md->max_bp_span;
1368 S = vc->sequence_encoding2;
1369 n = vc->length;
1370
1371 for (j = i; j <= MIN2(i + pairSize, n); j++) {
1372 type = md->pair[S[i]][S[j]];
1373 ptype[i][j] = (char)type;
1374 }
1375 return;
1376 }
1377
1378
1379 #if 0
1380 PRIVATE vrna_ep_t *
1381 get_deppp(vrna_fold_compound_t *vc,
1382 vrna_ep_t *pl,
1383 int start)
1384 {
1385 /* compute dependent pair probabilities */
1386 int i, j, count = 0;
1387 double tmp;
1388 vrna_ep_t *temp;
1389 char **ptype;
1390 short *S1;
1391 FLT_OR_DBL **qb, *scale;
1392 int *rtype, turn, pairsize, length;
1393
1394 vrna_exp_param_t *pf_params;
1395
1396 S1 = vc->sequence_encoding;
1397 pf_params = vc->exp_params;
1398 ptype = vc->ptype_local;
1399 qb = vc->exp_matrices->qb_local;
1400 scale = vc->exp_matrices->scale;
1401 rtype = &(pf_params->model_details.rtype[0]);
1402 turn = pf_params->model_details.min_loop_size;
1403 pairsize = pf_params->model_details.max_bp_span;
1404 length = vc->length;
1405
1406 temp = (vrna_ep_t *)vrna_alloc(pairsize * sizeof(vrna_ep_t)); /* holds temporary deppp */
1407 for (j = start + turn; j < MIN2(start + pairsize, length); j++) {
1408 if ((qb[start][j] * qb[start - 1][(j + 1)]) > 10e-200) {
1409 int type = ptype[start - 1][j + 1];
1410 int type_2 = rtype[(unsigned char)ptype[start][j]];
1411 tmp = qb[start][j] / qb[start - 1][(j + 1)] * exp_E_IntLoop(0,
1412 0,
1413 type,
1414 type_2,
1415 S1[start],
1416 S1[j],
1417 S1[start - 1],
1418 S1[j + 1],
1419 pf_params) * scale[2];
1420 temp[count].i = start;
1421 temp[count].j = j;
1422 temp[count++].p = tmp;
1423 }
1424 }
1425 /* write it to list of deppps */
1426 for (i = 0; pl[i].i != 0; i++);
1427 pl = (vrna_ep_t *)vrna_realloc(pl, (i + count + 1) * sizeof(vrna_ep_t));
1428 for (j = 0; j < count; j++) {
1429 pl[i + j].i = temp[j].i;
1430 pl[i + j].j = temp[j].j;
1431 pl[i + j].p = temp[j].p;
1432 }
1433 pl[i + count].i = 0;
1434 pl[i + count].j = 0;
1435 pl[i + count].p = 0;
1436 free(temp);
1437 return pl;
1438 }
1439
1440
1441 #endif
1442
1443 PRIVATE FLT_OR_DBL *
compute_stack_probabilities(vrna_fold_compound_t * vc,int start)1444 compute_stack_probabilities(vrna_fold_compound_t *vc,
1445 int start)
1446 {
1447 /* compute dependent pair probabilities */
1448 char **ptype;
1449 short *S1;
1450 int j, max_j, *rtype, turn, pairsize, length, type, type_2;
1451 FLT_OR_DBL **qb, *scale, *probs;
1452 double tmp;
1453 vrna_exp_param_t *pf_params;
1454
1455 length = vc->length;
1456 S1 = vc->sequence_encoding;
1457 pf_params = vc->exp_params;
1458 ptype = vc->ptype_local;
1459 qb = vc->exp_matrices->qb_local;
1460 scale = vc->exp_matrices->scale;
1461 rtype = &(pf_params->model_details.rtype[0]);
1462 turn = pf_params->model_details.min_loop_size;
1463 pairsize = pf_params->model_details.max_bp_span;
1464
1465 max_j = MIN2(start + pairsize, length) - 1;
1466
1467 probs = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (max_j - start + 1));
1468
1469 for (j = start + turn + 1; j <= max_j; j++) {
1470 if ((qb[start][j] * qb[start - 1][(j + 1)]) > 10e-200) {
1471 type = vrna_get_ptype_window(start - 1, j + 1 + start - 1, ptype);
1472 type_2 = rtype[vrna_get_ptype_window(start, j + start, ptype)];
1473 tmp = qb[start][j] /
1474 qb[start - 1][(j + 1)] *
1475 exp_E_IntLoop(0,
1476 0,
1477 type,
1478 type_2,
1479 S1[start],
1480 S1[j],
1481 S1[start - 1],
1482 S1[j + 1],
1483 pf_params) *
1484 scale[2];
1485 probs[j - start - 1] = tmp;
1486 }
1487 }
1488 return probs;
1489 }
1490
1491
1492 /*
1493 * Here: Space for questions...
1494 */
1495 PRIVATE void
compute_pU(vrna_fold_compound_t * vc,int k,int ulength,helper_arrays * aux_arrays,vrna_probs_window_callback * cb,void * data,unsigned int options)1496 compute_pU(vrna_fold_compound_t *vc,
1497 int k,
1498 int ulength,
1499 helper_arrays *aux_arrays,
1500 vrna_probs_window_callback *cb,
1501 void *data,
1502 unsigned int options)
1503 {
1504 /*
1505 * here, we try to add a function computing all unpaired probabilities starting at some i,
1506 * going down to $unpaired, to be unpaired, i.e. a list with entries from 1 to unpaired for
1507 * every i, with the probability of a stretch of length x, starting at i-x+1, to be unpaired
1508 */
1509 char **ptype;
1510 short *S1;
1511 int startu, i5, j3, len, obp, *rtype, turn, winSize, n, leftmost,
1512 rightmost, tt;
1513 FLT_OR_DBL expMLclosing, *expMLbase, **q, **qm, **qm2, *scale, **pR, **QI5,
1514 **q2l, **qmb;
1515 double qqq, temp, *QBE, *QBI, *QBM, *QBH, **pU, **pUO, **pUH, **pUI, **pUM;
1516 vrna_exp_param_t *pf_params;
1517 vrna_hc_t *hc;
1518 vrna_sc_t *sc;
1519
1520 n = vc->length;
1521 winSize = vc->window_size;
1522 S1 = vc->sequence_encoding;
1523 pf_params = vc->exp_params;
1524 ptype = vc->ptype_local;
1525 rtype = &(pf_params->model_details.rtype[0]);
1526 scale = vc->exp_matrices->scale;
1527 q = vc->exp_matrices->q_local;
1528 qm = vc->exp_matrices->qm_local;
1529 qm2 = vc->exp_matrices->qm2_local;
1530 expMLbase = vc->exp_matrices->expMLbase;
1531 expMLclosing = pf_params->expMLclosing;
1532 pR = vc->exp_matrices->pR;
1533 QI5 = vc->exp_matrices->QI5;
1534 q2l = vc->exp_matrices->q2l;
1535 qmb = vc->exp_matrices->qmb;
1536 turn = pf_params->model_details.min_loop_size;
1537 hc = vc->hc;
1538 sc = vc->sc;
1539
1540 pU = aux_arrays->pU;
1541 pUO = aux_arrays->pUO;
1542 pUH = aux_arrays->pUH;
1543 pUI = aux_arrays->pUI;
1544 pUM = aux_arrays->pUM;
1545
1546 QBE = (double *)vrna_alloc((MAX2(ulength, MAXLOOP) + 2) * sizeof(double));
1547 QBM = (double *)vrna_alloc((MAX2(ulength, MAXLOOP) + 2) * sizeof(double));
1548 QBI = (double *)vrna_alloc((MAX2(ulength, MAXLOOP) + 2) * sizeof(double));
1549 QBH = (double *)vrna_alloc((MAX2(ulength, MAXLOOP) + 2) * sizeof(double));
1550
1551 /*
1552 * first, we will
1553 * for k<=ulength, pU[k][k]=0, because no bp can enclose it
1554 */
1555
1556 /* compute pu[k+ulength][ulength] */
1557 for (i5 = MAX2(k + ulength - winSize + 1, 1); i5 <= k; i5++) {
1558 for (j3 = k + ulength + 1; j3 <= MIN2(n, i5 + winSize - 1); j3++) {
1559 /* Multiloops */
1560 if (hc->matrix_local[i5][j3 - i5] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1561 tt = rtype[vrna_get_ptype_window(i5, j3 + i5, ptype)];
1562 temp = 0.;
1563 /*
1564 * (.. >-----|..........)
1565 * i5 j j+ulength j3
1566 */
1567 /* (..{}{}-----|......) */
1568 if ((hc->up_ml[k + 1] >= j3 - k - 1) && (i5 < k)) {
1569 qqq = qm2[i5 + 1][k] * expMLbase[j3 - k - 1];
1570
1571 if (sc) {
1572 if (sc->exp_energy_up)
1573 qqq *= sc->exp_energy_up[k + 1][j3 - k - 1];
1574
1575 if (sc->f)
1576 qqq *= sc->f(i5, j3, i5 + 1, k, VRNA_DECOMP_PAIR_ML, sc->data);
1577 }
1578
1579 temp += qqq;
1580 }
1581
1582 /* (..|-----|{}{}) */
1583 if ((hc->up_ml[i5 + 1] >= k + ulength - i5) && (j3 - 1 > k + ulength)) {
1584 qqq = qm2[k + ulength + 1][j3 - 1] *
1585 expMLbase[k + ulength - i5];
1586
1587 if (sc) {
1588 if (sc->exp_energy_up)
1589 qqq *= sc->exp_energy_up[i5 + 1][k + ulength - i5];
1590
1591 if (sc->f)
1592 qqq *= sc->f(i5, j3, k + ulength + 1, j3, VRNA_DECOMP_PAIR_ML, sc->data);
1593 }
1594
1595 temp += qqq;
1596 }
1597
1598 /* ({}|-----|{}) */
1599 if ((hc->up_ml[k + 1] >= ulength) && (i5 < k) && (j3 - 1 > k + ulength)) {
1600 qqq = qm[i5 + 1][k] *
1601 qm[k + ulength + 1][j3 - 1] *
1602 expMLbase[ulength];
1603
1604 if (sc) {
1605 if (sc->exp_energy_up)
1606 qqq *= sc->exp_energy_up[k + 1][ulength];
1607
1608 if (sc->f)
1609 qqq *= sc->f(i5, j3, k, k + ulength + 1, VRNA_DECOMP_PAIR_ML_OUTSIDE, sc->data);
1610 }
1611
1612 temp += qqq;
1613 }
1614
1615 /* add dangles, multloopclosing etc. */
1616 qqq = exp_E_MLstem(tt,
1617 S1[j3 - 1],
1618 S1[i5 + 1],
1619 pf_params) *
1620 scale[2] *
1621 expMLclosing;
1622
1623 if (sc)
1624 if (sc->exp_energy_bp_local)
1625 qqq *= sc->exp_energy_bp_local[i5][j3 - i5];
1626
1627 temp *= qqq;
1628
1629 pU[k + ulength][ulength] += temp * pR[i5][j3];
1630
1631 if (options & VRNA_PROBS_WINDOW_UP_SPLIT)
1632 pUM[k + ulength][ulength] += temp * pR[i5][j3];
1633 }
1634
1635 /* add hairpins */
1636 if (hc->matrix_local[i5][j3 - i5] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) {
1637 temp = vrna_exp_E_hp_loop(vc, i5, j3);
1638 pU[k + ulength][ulength] += temp * pR[i5][j3];
1639
1640 if (options & VRNA_PROBS_WINDOW_UP_SPLIT)
1641 pUH[k + ulength][ulength] += temp * pR[i5][j3];
1642 }
1643 }
1644 }
1645
1646 /* Add Interior loop contribution to QBE (and QBI) */
1647 temp = 0.;
1648 for (len = winSize; len > MAX2(ulength, MAXLOOP); len--)
1649 temp += QI5[k][len];
1650
1651 for (; len > 0; len--) {
1652 temp += QI5[k][len];
1653 QBI[len] += temp;
1654 QBE[len] += temp;
1655 }
1656
1657 /* Add Hairpin loop contribution to QBE (and QBH) */
1658 temp = 0.;
1659 for (obp = MIN2(n, k + winSize - 1); obp > k + ulength; obp--)
1660 temp += pR[k][obp] *
1661 vrna_exp_E_hp_loop(vc, k, obp);
1662
1663 for (obp = MIN2(n, MIN2(k + winSize - 1, k + ulength)); obp > k + 1; obp--) {
1664 temp += pR[k][obp] *
1665 vrna_exp_E_hp_loop(vc, k, obp);
1666 QBH[obp - k - 1] += temp;
1667 QBE[obp - k - 1] += temp;
1668 }
1669
1670 /*
1671 * Add up Multiloopterms qmb[l][m]+=prml[m]*dang;
1672 * q2l[l][m]+=(prml[m]-prm_l[m])*dang;
1673 */
1674
1675 temp = 0.;
1676
1677 /* add (()()____) type cont. to I3 */
1678 if (sc && sc->exp_energy_up) {
1679 for (len = winSize; len >= ulength; len--)
1680 if (hc->up_ml[k + 1] >= len) {
1681 temp += q2l[k][len] *
1682 expMLbase[len] *
1683 sc->exp_energy_up[k + 1][len];
1684 }
1685
1686 for (; len > 0; len--) {
1687 if (hc->up_ml[k + 1] >= len) {
1688 temp += q2l[k][len] *
1689 expMLbase[len] *
1690 sc->exp_energy_up[k + 1][len];
1691 }
1692
1693 QBM[len] += temp;
1694 QBE[len] += temp;
1695 }
1696 } else {
1697 for (len = winSize; len >= ulength; len--)
1698 if (hc->up_ml[k + 1] >= len)
1699 temp += q2l[k][len] *
1700 expMLbase[len];
1701
1702 for (; len > 0; len--) {
1703 if (hc->up_ml[k + 1] >= len)
1704 temp += q2l[k][len] *
1705 expMLbase[len];
1706
1707 QBM[len] += temp;
1708 QBE[len] += temp;
1709 }
1710 }
1711
1712 /* add (()___()) */
1713 for (len = 1; len < ulength; len++) {
1714 if (hc->up_ml[k + 1] >= len) {
1715 for (obp = k + len + turn; obp <= MIN2(n, k + winSize - 1); obp++) {
1716 temp = qmb[k][obp - k - 1] *
1717 qm[k + len + 1 /*2*/][obp - 1] *
1718 expMLbase[len];
1719
1720 if (sc)
1721 if (sc->exp_energy_up)
1722 temp *= sc->exp_energy_up[k + 1][len];
1723
1724 QBM[len] += temp;
1725 QBE[len] += temp;
1726 }
1727 }
1728 }
1729
1730 /* add (___()()) */
1731 for (len = 1; len < ulength; len++) {
1732 if (hc->up_ml[k + 1] >= len) {
1733 for (obp = k + len + turn + turn; obp <= MIN2(n, k + winSize - 1); obp++) {
1734 if (hc->matrix_local[k][obp - k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1735 tt = rtype[vrna_get_ptype_window(k, obp + k, ptype)];
1736 temp = exp_E_MLstem(tt, S1[obp - 1], S1[k + 1], pf_params) *
1737 scale[2] *
1738 expMLbase[len] *
1739 expMLclosing *
1740 pR[k][obp] *
1741 qm2[k + len + 1][obp - 1]; /* k:obp */
1742
1743 if (sc) {
1744 if (sc->exp_energy_up)
1745 temp *= sc->exp_energy_up[k + 1][len];
1746
1747 if (sc->exp_energy_bp)
1748 temp *= sc->exp_energy_bp_local[k][obp - k];
1749 }
1750
1751 QBM[len] += temp;
1752 QBE[len] += temp;
1753 }
1754 }
1755 }
1756 }
1757
1758 /*
1759 * After computing all these contributions in QBE[len], that k is paired
1760 * and the unpaired stretch is AT LEAST len long, we start to add that to
1761 * the old unpaired thingies;
1762 */
1763 for (len = 1; len <= MIN2(MAX2(ulength, MAXLOOP), n - k); len++)
1764 pU[k + len][len] += pU[k + len][len + 1] + QBE[len];
1765
1766 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
1767 for (len = 1; len <= MIN2(MAX2(ulength, MAXLOOP), n - k); len++) {
1768 pUH[k + len][len] += pUH[k + len][len + 1] + QBH[len];
1769 pUM[k + len][len] += pUM[k + len][len + 1] + QBM[len];
1770 pUI[k + len][len] += pUI[k + len][len + 1] + QBI[len];
1771 }
1772 /* open chain */
1773 if ((ulength >= winSize) && (k >= ulength) && (hc->up_ext[k - winSize + 1] >= winSize))
1774 pUO[k][winSize] = scale[winSize] / q[k - winSize + 1][k];
1775 }
1776
1777 /* open chain */
1778 if ((ulength >= winSize) && (k >= ulength) && (hc->up_ext[k - winSize + 1] >= winSize)) {
1779 if (sc && sc->exp_energy_up) {
1780 pU[k][winSize] = scale[winSize] *
1781 sc->exp_energy_up[k][winSize] /
1782 q[k - winSize + 1][k];
1783 } else {
1784 pU[k][winSize] = scale[winSize] / q[k - winSize + 1][k];
1785 }
1786 }
1787
1788 /*
1789 * now the not enclosed by any base pair terms for whatever it is we do not need anymore...
1790 * ... which should be e.g; k, again
1791 */
1792 for (startu = MIN2(ulength, k); startu > 0; startu--) {
1793 temp = 0.;
1794 /* check whether soft constraint unpaired contributions available */
1795 if (sc && sc->exp_energy_up) {
1796 if (hc->up_ext[k - startu + 1] >= startu) {
1797 for (i5 = MAX2(1, k - winSize + 2); i5 <= MIN2(k - startu, n - winSize + 1); i5++)
1798 temp += q[i5][k - startu] *
1799 q[k + 1][i5 + winSize - 1] *
1800 scale[startu] *
1801 sc->exp_energy_up[k - startu + 1][startu] /
1802 q[i5][i5 + winSize - 1];
1803
1804 /* the 2 Cases where the borders are on the edge of the interval */
1805 if ((k >= winSize) && (startu + 1 <= winSize)) {
1806 temp += q[k - winSize + 1][k - startu] *
1807 scale[startu] *
1808 sc->exp_energy_up[k - startu + 1][startu] /
1809 q[k - winSize + 1][k];
1810 }
1811
1812 if ((k <= n - winSize + startu) && (k - startu >= 0) && (k < n) &&
1813 (startu + 1 <= winSize)) {
1814 temp += q[k + 1][k - startu + winSize] *
1815 scale[startu] *
1816 sc->exp_energy_up[k - startu + 1][startu] /
1817 q[k - startu + 1][k - startu + winSize];
1818 }
1819 }
1820 } else {
1821 if (hc->up_ext[k - startu + 1] >= startu) {
1822 for (i5 = MAX2(1, k - winSize + 2); i5 <= MIN2(k - startu, n - winSize + 1); i5++)
1823 temp += q[i5][k - startu] *
1824 q[k + 1][i5 + winSize - 1] *
1825 scale[startu] /
1826 q[i5][i5 + winSize - 1];
1827
1828 /* the 2 Cases where the borders are on the edge of the interval */
1829 if ((k >= winSize) && (startu + 1 <= winSize))
1830 temp += q[k - winSize + 1][k - startu] *
1831 scale[startu] /
1832 q[k - winSize + 1][k];
1833
1834 if ((k <= n - winSize + startu) && (k - startu >= 0) && (k < n) && (startu + 1 <= winSize))
1835 temp += q[k + 1][k - startu + winSize] *
1836 scale[startu] /
1837 q[k - startu + 1][k - startu + winSize];
1838 }
1839 }
1840
1841 /* Divide by number of possible windows */
1842 leftmost = MAX2(1, k - winSize + 1);
1843 rightmost = MIN2(n - winSize + 1, k - startu + 1);
1844
1845 pU[k][startu] += temp;
1846 pU[k][startu] /= (rightmost - leftmost + 1);
1847
1848 if (options & VRNA_PROBS_WINDOW_UP_SPLIT) {
1849 pUO[k][startu] += temp;
1850
1851 /* Do we want to make a distinction between those? */
1852 pUO[k][startu] /= (rightmost - leftmost + 1);
1853 pUH[k][startu] /= (rightmost - leftmost + 1);
1854 pUI[k][startu] /= (rightmost - leftmost + 1);
1855 pUM[k][startu] /= (rightmost - leftmost + 1);
1856 }
1857 }
1858 free(QBE);
1859 free(QBI);
1860 free(QBH);
1861 free(QBM);
1862
1863 /* call return callback */
1864 return_pU(MIN2(ulength, k), k, ulength, aux_arrays, cb, data, options);
1865
1866 return;
1867 }
1868
1869
1870 PRIVATE void
print_bpp_callback(FLT_OR_DBL * pr,int size,int k,void * data)1871 print_bpp_callback(FLT_OR_DBL *pr,
1872 int size,
1873 int k,
1874 void *data)
1875 {
1876 int j;
1877 FILE *fp = ((default_cb_data *)data)->fp_bpp;
1878 FLT_OR_DBL cutoff = ((default_cb_data *)data)->bpp_cutoff;
1879
1880 for (j = k + 1; j <= size; j++) {
1881 if (pr[j] < cutoff)
1882 continue;
1883
1884 fprintf(fp, "%d %d %g\n", k, j, pr[j]);
1885 }
1886 }
1887
1888
1889 PRIVATE void
store_bpp_callback(FLT_OR_DBL * pr,int size,int k,void * data)1890 store_bpp_callback(FLT_OR_DBL *pr,
1891 int size,
1892 int k,
1893 void *data)
1894 {
1895 int j;
1896 vrna_ep_t *pl = ((default_cb_data *)data)->bpp;
1897 unsigned int pl_size = ((default_cb_data *)data)->bpp_size;
1898 unsigned int pl_max_size = ((default_cb_data *)data)->bpp_max_size;
1899 FLT_OR_DBL cutoff = ((default_cb_data *)data)->bpp_cutoff;
1900
1901 if (pl_max_size == 0) {
1902 /* init if necessary */
1903 pl_max_size = 100;
1904 pl = (vrna_ep_t *)vrna_realloc(pl, sizeof(vrna_ep_t) * pl_max_size);
1905 }
1906
1907 for (j = k + 1; j <= size; j++) {
1908 if (pr[j] < cutoff)
1909 continue;
1910
1911 /* resize vrna_ep_t memory if necessary */
1912 if (pl_size >= pl_max_size - 1) {
1913 pl_max_size *= 1.5;
1914 pl = (vrna_ep_t *)vrna_realloc(pl, sizeof(vrna_ep_t) * pl_max_size);
1915 }
1916
1917 pl[pl_size].i = k;
1918 pl[pl_size].j = j;
1919 pl[pl_size].type = VRNA_PLIST_TYPE_BASEPAIR;
1920 pl[pl_size++].p = pr[j];
1921 }
1922
1923 /* mark end of vrna_ep_t */
1924 pl[pl_size].i = 0;
1925 pl[pl_size].j = 0;
1926 pl[pl_size].type = VRNA_PLIST_TYPE_BASEPAIR;
1927 pl[pl_size].p = 0.;
1928
1929 /* update data */
1930 ((default_cb_data *)data)->bpp = pl;
1931 ((default_cb_data *)data)->bpp_size = pl_size;
1932 ((default_cb_data *)data)->bpp_max_size = pl_max_size;
1933 }
1934
1935
1936 #if 0
1937 PRIVATE void
1938 store_stack_prob_callback(FLT_OR_DBL *pr,
1939 int size,
1940 int k,
1941 void *data)
1942 {
1943 int j;
1944 vrna_ep_t *pl = ((default_cb_data *)data)->stack_prob;
1945 unsigned int pl_size = ((default_cb_data *)data)->stack_prob_size;
1946 unsigned int pl_max_size = ((default_cb_data *)data)->stack_prob_max_size;
1947 FLT_OR_DBL cutoff = ((default_cb_data *)data)->bpp_cutoff;
1948
1949 if (pl_max_size == 0) {
1950 /* init if necessary */
1951 pl_max_size = 100;
1952 pl = (vrna_ep_t *)vrna_realloc(pl, sizeof(vrna_ep_t) * pl_max_size);
1953 }
1954
1955 for (j = k + 1; j <= size; j++) {
1956 if (pr[j] < cutoff)
1957 continue;
1958
1959 /* resize vrna_ep_t memory if necessary */
1960 if (pl_size >= pl_max_size - 1) {
1961 pl_max_size *= 1.5;
1962 pl = (vrna_ep_t *)vrna_realloc(pl, sizeof(vrna_ep_t) * pl_max_size);
1963 }
1964
1965 pl[pl_size].i = k;
1966 pl[pl_size].j = j;
1967 pl[pl_size].type = VRNA_PLIST_TYPE_BASEPAIR;
1968 pl[pl_size++].p = pr[j];
1969 }
1970
1971 /* mark end of vrna_ep_t */
1972 pl[pl_size].i = 0;
1973 pl[pl_size].j = 0;
1974 pl[pl_size].type = VRNA_PLIST_TYPE_BASEPAIR;
1975 pl[pl_size].p = 0.;
1976
1977 /* update data */
1978 ((default_cb_data *)data)->stack_prob = pl;
1979 ((default_cb_data *)data)->stack_prob_size = pl_size;
1980 ((default_cb_data *)data)->stack_prob_max_size = pl_max_size;
1981 }
1982
1983
1984 #endif
1985
1986
1987 PRIVATE void
print_pU_callback(double * pU,int size,int k,int ulength,unsigned int type,void * data)1988 print_pU_callback(double *pU,
1989 int size,
1990 int k,
1991 int ulength,
1992 unsigned int type,
1993 void *data)
1994 {
1995 if (type & VRNA_PROBS_WINDOW_UP) {
1996 int i;
1997 FILE *fp = ((default_cb_data *)data)->fp_pU;
1998
1999 fprintf(fp, "%d\t", k);
2000
2001 for (i = 1; i < size; i++)
2002 fprintf(fp, "%.7g\t", pU[i]);
2003 fprintf(fp, "%.7g", pU[size]);
2004
2005 if ((type & VRNA_ANY_LOOP) == VRNA_ANY_LOOP)
2006 fprintf(fp, "\n");
2007 else if (type & VRNA_EXT_LOOP)
2008 fprintf(fp, "\tE\n");
2009 else if (type & VRNA_HP_LOOP)
2010 fprintf(fp, "\tH\n");
2011 else if (type & VRNA_INT_LOOP)
2012 fprintf(fp, "\tI\n");
2013 else if (type & VRNA_MB_LOOP)
2014 fprintf(fp, "\tM\n");
2015 else
2016 vrna_message_warning("unknown loop type");
2017 }
2018 }
2019
2020
2021 PRIVATE void
store_pU_callback(double * pU,int size,int k,int ulength,unsigned int type,void * data)2022 store_pU_callback(double *pU,
2023 int size,
2024 int k,
2025 int ulength,
2026 unsigned int type,
2027 void *data)
2028 {
2029 int i;
2030 double **pU_storage = ((default_cb_data *)data)->pU;
2031
2032 if ((type & VRNA_PROBS_WINDOW_UP) && ((type & VRNA_ANY_LOOP) == VRNA_ANY_LOOP)) {
2033 pU_storage[k] = (double *)vrna_alloc(sizeof(double) * (ulength + 1));
2034 for (i = 1; i <= size; i++)
2035 pU_storage[k][i] = pU[i];
2036 }
2037 }
2038
2039
2040 PRIVATE void
backward_compat_callback(FLT_OR_DBL * pr,int pr_size,int i,int max,unsigned int type,void * data)2041 backward_compat_callback(FLT_OR_DBL *pr,
2042 int pr_size,
2043 int i,
2044 int max,
2045 unsigned int type,
2046 void *data)
2047 {
2048 default_cb_data *d = (default_cb_data *)data;
2049
2050 if (type & VRNA_PROBS_WINDOW_BPP) {
2051 if (d->bpp_print)
2052 print_bpp_callback(pr, pr_size, i, data);
2053 else
2054 store_bpp_callback(pr, pr_size, i, data);
2055 } else if (type & VRNA_PROBS_WINDOW_UP) {
2056 if (d->up_print)
2057 print_pU_callback(pr, pr_size, i, max, type, data);
2058 else
2059 store_pU_callback(pr, pr_size, i, max, type, data);
2060 }
2061 }
2062
2063
2064 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
2065
2066 /*
2067 *###########################################
2068 *# deprecated functions below #
2069 *###########################################
2070 */
2071
2072 PRIVATE void
2073 putoutpU_prob_old(double **pU,
2074 int length,
2075 int ulength,
2076 FILE *fp,
2077 int energies,
2078 vrna_exp_param_t *parameters);
2079
2080
2081 PRIVATE void
2082 putoutpU_prob_bin_old(double **pU,
2083 int length,
2084 int ulength,
2085 FILE *fp,
2086 int energies,
2087 vrna_exp_param_t *parameters);
2088
2089
2090 PRIVATE vrna_ep_t *
wrap_pf_foldLP(char * sequence,int winSize,int pairSize,float cutoffb,double ** pU,vrna_ep_t ** dpp2,FILE * pUfp,FILE * spup,vrna_exp_param_t * parameters)2091 wrap_pf_foldLP(char *sequence,
2092 int winSize,
2093 int pairSize,
2094 float cutoffb,
2095 double **pU,
2096 vrna_ep_t **dpp2,
2097 FILE *pUfp,
2098 FILE *spup,
2099 vrna_exp_param_t *parameters)
2100 {
2101 int ulength, r;
2102 vrna_fold_compound_t *vc;
2103 vrna_md_t md;
2104 default_cb_data data;
2105
2106 vc = NULL;
2107 ulength = 0;
2108
2109 /*
2110 * if present, extract model details from provided parameters variable,
2111 * to properly initialize the fold compound. Otherwise use default
2112 * settings taken from deprecated global variables
2113 */
2114 if (parameters)
2115 vrna_md_copy(&md, &(parameters->model_details));
2116 else
2117 set_model_details(&md);
2118
2119 md.compute_bpp = 1; /* turn on base pair probability computations */
2120 md.window_size = winSize; /* set size of sliding window */
2121 md.max_bp_span = pairSize; /* set maximum base pair span */
2122
2123 vc = vrna_fold_compound(sequence, &md, VRNA_OPTION_DEFAULT | VRNA_OPTION_WINDOW);
2124
2125 /*
2126 * if present, attach a copy of the parameters structure instead of the
2127 * default parameters but take care of re-setting it to (initialized)
2128 * model details
2129 */
2130 free(vc->exp_params);
2131 if (parameters) {
2132 vrna_md_copy(&(parameters->model_details), &(vc->params->model_details));
2133 vc->exp_params = vrna_exp_params_copy(parameters);
2134 } else {
2135 vc->exp_params = vrna_exp_params(&(vc->params->model_details));
2136 }
2137
2138 /* propagate global pf_scale into vc->exp_params */
2139 vc->exp_params->pf_scale = pf_scale;
2140
2141 if (backward_compat_compound && backward_compat)
2142 vrna_fold_compound_free(backward_compat_compound);
2143
2144 backward_compat_compound = vc;
2145 backward_compat = 1;
2146 iindx = backward_compat_compound->iindx; /* for backward compatibility and Perl wrapper */
2147
2148 if (pU)
2149 ulength = (int)pU[0][0] + 0.49;
2150
2151 data.fp_pU = pUfp;
2152 data.pU = pU;
2153 data.bpp_cutoff = (FLT_OR_DBL)cutoffb;
2154 data.fp_bpp = spup;
2155 data.bpp = NULL;
2156 data.bpp_max_size = 0;
2157 data.bpp_size = 0;
2158 data.stack_prob = NULL;
2159 data.stack_prob_max_size = 0;
2160 data.stack_prob_size = 0;
2161
2162 data.bpp_print = (spup) ? 1 : 0;
2163 data.up_print = (pUfp) ? 1 : 0;
2164
2165 unsigned int options = VRNA_PROBS_WINDOW_BPP; /* always compute base pair probabilities */
2166
2167 if (dpp2 && (*dpp2))
2168 options |= VRNA_PROBS_WINDOW_STACKP;
2169
2170 if (ulength > 0)
2171 options |= VRNA_PROBS_WINDOW_UP;
2172
2173 r = vrna_probs_window(vc, ulength, options, &backward_compat_callback, (void *)&data);
2174
2175 if (!r)
2176 return NULL;
2177
2178 if (dpp2 && (*dpp2)) {
2179 data.stack_prob = (vrna_ep_t *)vrna_realloc(data.stack_prob,
2180 sizeof(vrna_ep_t) *
2181 (data.stack_prob_size + 1));
2182 data.stack_prob[data.stack_prob_size].i = 0;
2183 data.stack_prob[data.stack_prob_size].j = 0;
2184 data.stack_prob[data.stack_prob_size].type = VRNA_PLIST_TYPE_BASEPAIR;
2185 data.stack_prob[data.stack_prob_size].p = 0;
2186 free(*dpp2); /* free already occupied memory */
2187 *dpp2 = data.stack_prob;
2188 }
2189
2190 if (!spup) {
2191 data.bpp =
2192 (vrna_ep_t *)vrna_realloc(data.bpp, sizeof(vrna_ep_t) * (data.bpp_size + 1));
2193 data.bpp[data.bpp_size].i = 0;
2194 data.bpp[data.bpp_size].j = 0;
2195 data.bpp[data.bpp_size].type = VRNA_PLIST_TYPE_BASEPAIR;
2196 data.bpp[data.bpp_size].p = 0;
2197 return data.bpp;
2198 } else {
2199 return NULL;
2200 }
2201 }
2202
2203
2204 PUBLIC void
init_pf_foldLP(int length)2205 init_pf_foldLP(int length)
2206 {
2207 /* DO NOTHING */
2208 }
2209
2210
2211 PUBLIC void
update_pf_paramsLP(int length)2212 update_pf_paramsLP(int length)
2213 {
2214 if (backward_compat_compound && backward_compat) {
2215 vrna_md_t md;
2216 set_model_details(&md);
2217 vrna_exp_params_reset(backward_compat_compound, &md);
2218
2219 /* compatibility with RNAup, may be removed sometime */
2220 pf_scale = backward_compat_compound->exp_params->pf_scale;
2221 }
2222 }
2223
2224
2225 PUBLIC void
update_pf_paramsLP_par(int length,vrna_exp_param_t * parameters)2226 update_pf_paramsLP_par(int length,
2227 vrna_exp_param_t *parameters)
2228 {
2229 if (backward_compat_compound && backward_compat) {
2230 vrna_md_t md;
2231 if (parameters) {
2232 vrna_exp_params_subst(backward_compat_compound, parameters);
2233 } else {
2234 set_model_details(&md);
2235 vrna_exp_params_reset(backward_compat_compound, &md);
2236 }
2237
2238 /* compatibility with RNAup, may be removed sometime */
2239 pf_scale = backward_compat_compound->exp_params->pf_scale;
2240 }
2241 }
2242
2243
2244 PUBLIC vrna_ep_t *
pfl_fold(char * sequence,int winSize,int pairSize,float cutoffb,double ** pU,vrna_ep_t ** dpp2,FILE * pUfp,FILE * spup)2245 pfl_fold(char *sequence,
2246 int winSize,
2247 int pairSize,
2248 float cutoffb,
2249 double **pU,
2250 vrna_ep_t **dpp2,
2251 FILE *pUfp,
2252 FILE *spup)
2253 {
2254 return wrap_pf_foldLP(sequence, winSize, pairSize, cutoffb, pU, dpp2, pUfp, spup, NULL);
2255 }
2256
2257
2258 PUBLIC vrna_ep_t *
pfl_fold_par(char * sequence,int winSize,int pairSize,float cutoffb,double ** pU,vrna_ep_t ** dpp2,FILE * pUfp,FILE * spup,vrna_exp_param_t * parameters)2259 pfl_fold_par(char *sequence,
2260 int winSize,
2261 int pairSize,
2262 float cutoffb,
2263 double **pU,
2264 vrna_ep_t **dpp2,
2265 FILE *pUfp,
2266 FILE *spup,
2267 vrna_exp_param_t *parameters)
2268 {
2269 return wrap_pf_foldLP(sequence, winSize, pairSize, cutoffb, pU, dpp2, pUfp, spup, parameters);
2270 }
2271
2272
2273 PUBLIC void
putoutpU_prob(double ** pU,int length,int ulength,FILE * fp,int energies)2274 putoutpU_prob(double **pU,
2275 int length,
2276 int ulength,
2277 FILE *fp,
2278 int energies)
2279 {
2280 if (backward_compat_compound && backward_compat)
2281 putoutpU_prob_old(pU, length, ulength, fp, energies, backward_compat_compound->exp_params);
2282 else
2283 vrna_message_warning("putoutpU_prob: Not doing anything! First, run pfl_fold()!");
2284 }
2285
2286
2287 PUBLIC void
putoutpU_prob_par(double ** pU,int length,int ulength,FILE * fp,int energies,vrna_exp_param_t * parameters)2288 putoutpU_prob_par(double **pU,
2289 int length,
2290 int ulength,
2291 FILE *fp,
2292 int energies,
2293 vrna_exp_param_t *parameters)
2294 {
2295 if ((pU) && (fp) && (parameters))
2296 putoutpU_prob_old(pU, length, ulength, fp, energies, parameters);
2297 }
2298
2299
2300 PRIVATE void
putoutpU_prob_old(double ** pU,int length,int ulength,FILE * fp,int energies,vrna_exp_param_t * parameters)2301 putoutpU_prob_old(double **pU,
2302 int length,
2303 int ulength,
2304 FILE *fp,
2305 int energies,
2306 vrna_exp_param_t *parameters)
2307 {
2308 /* put out unpaireds */
2309 int i, k;
2310 double temp, kT = parameters->kT / 1000.0;
2311
2312 if (energies)
2313 fprintf(fp, "#opening energies\n #i$\tl=");
2314 else
2315 fprintf(fp, "#unpaired probabilities\n #i$\tl=");
2316
2317 for (i = 1; i <= ulength; i++)
2318 fprintf(fp, "%d\t", i);
2319 fprintf(fp, "\n");
2320
2321 for (k = 1; k <= length; k++) {
2322 fprintf(fp, "%d\t", k);
2323 for (i = 1; i <= ulength; i++) {
2324 if (i > k) {
2325 fprintf(fp, "NA\t");
2326 continue;
2327 }
2328
2329 if (energies)
2330 temp = -log(pU[k][i]) * kT;
2331 else
2332 temp = pU[k][i];
2333
2334 fprintf(fp, "%.7g\t", temp);
2335 }
2336 fprintf(fp, "\n");
2337 free(pU[k]);
2338 }
2339 fflush(fp);
2340 }
2341
2342
2343 PUBLIC void
putoutpU_prob_bin(double ** pU,int length,int ulength,FILE * fp,int energies)2344 putoutpU_prob_bin(double **pU,
2345 int length,
2346 int ulength,
2347 FILE *fp,
2348 int energies)
2349 {
2350 if (backward_compat_compound && backward_compat)
2351 putoutpU_prob_bin_old(pU, length, ulength, fp, energies, backward_compat_compound->exp_params);
2352 else
2353 vrna_message_warning("putoutpU_prob_bin: Not doing anything! First, run pfl_fold()!");
2354 }
2355
2356
2357 PUBLIC void
putoutpU_prob_bin_par(double ** pU,int length,int ulength,FILE * fp,int energies,vrna_exp_param_t * parameters)2358 putoutpU_prob_bin_par(double **pU,
2359 int length,
2360 int ulength,
2361 FILE *fp,
2362 int energies,
2363 vrna_exp_param_t *parameters)
2364 {
2365 if ((pU) && (fp) && (parameters))
2366 putoutpU_prob_bin_old(pU, length, ulength, fp, energies, parameters);
2367 }
2368
2369
2370 PRIVATE void
putoutpU_prob_bin_old(double ** pU,int length,int ulength,FILE * fp,int energies,vrna_exp_param_t * parameters)2371 putoutpU_prob_bin_old(double **pU,
2372 int length,
2373 int ulength,
2374 FILE *fp,
2375 int energies,
2376 vrna_exp_param_t *parameters)
2377 {
2378 /* put out unpaireds */
2379 int i, k, *p;
2380 double kT = parameters->kT / 1000.0;
2381
2382 p = (int *)vrna_alloc(sizeof(int) * 1);
2383 /* write first line */
2384 p[0] = ulength; /* u length */
2385 fwrite(p, sizeof(int), 1, fp);
2386 p[0] = length; /* seq length */
2387 fwrite(p, sizeof(int), 1, fp);
2388 for (k = 3; k <= (length + 20); k++) {
2389 /* all the other lines are set to 1000000 because we are at ulength=0 */
2390 p[0] = 1000000;
2391 fwrite(p, sizeof(int), 1, fp);
2392 }
2393 /* data */
2394 for (i = 1; i <= ulength; i++) {
2395 for (k = 1; k <= 11; k++) {
2396 /* write first ten entries to 1000000 */
2397 p[0] = 1000000;
2398 fwrite(p, sizeof(int), 1, fp);
2399 }
2400 for (k = 1; k <= length; k++) {
2401 /* write data now */
2402 if (i > k) {
2403 p[0] = 1000000; /* check if u > pos */
2404 fwrite(p, sizeof(int), 1, fp);
2405 continue;
2406 } else {
2407 p[0] = (int)rint(100 * (-log(pU[k][i]) * kT));
2408 fwrite(p, sizeof(int), 1, fp);
2409 }
2410 }
2411 for (k = 1; k <= 9; k++) {
2412 /* finish by writing the last 10 entries */
2413 p[0] = 1000000;
2414 fwrite(p, sizeof(int), 1, fp);
2415 }
2416 }
2417 /* free pU array; */
2418 for (k = 1; k <= length; k++)
2419 free(pU[k]);
2420 free(p);
2421 fflush(fp);
2422 }
2423
2424
2425 #endif
2426