1 /* constraints handling */
2
3 #ifdef HAVE_CONFIG_H
4 #include "config.h"
5 #endif
6
7 #include <assert.h>
8 #include <config.h>
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <math.h>
12 #include <ctype.h>
13 #include <string.h>
14 #include <limits.h>
15
16 #include "ViennaRNA/params/default.h"
17 #include "ViennaRNA/params/constants.h" /* defines MINPSCORE */
18 #include "ViennaRNA/fold_vars.h"
19 #include "ViennaRNA/utils/basic.h"
20 #include "ViennaRNA/utils/alignments.h"
21 #include "ViennaRNA/io/file_formats.h"
22 #include "ViennaRNA/params/basic.h"
23 #include "ViennaRNA/constraints/SHAPE.h"
24 #include "ViennaRNA/constraints/soft.h"
25
26
27 #ifndef INLINE
28 # ifdef __GNUC__
29 # define INLINE inline
30 # else
31 # define INLINE
32 # endif
33 #endif
34
35 #define STATE_CLEAN (unsigned char)0
36 #define STATE_DIRTY_UP_MFE (unsigned char)1
37 #define STATE_DIRTY_UP_PF (unsigned char)2
38 #define STATE_DIRTY_BP_MFE (unsigned char)4
39 #define STATE_DIRTY_BP_PF (unsigned char)8
40
41 /*
42 #################################
43 # GLOBAL VARIABLES #
44 #################################
45 */
46
47 /*
48 #################################
49 # PRIVATE VARIABLES #
50 #################################
51 */
52
53 /*
54 #################################
55 # PRIVATE FUNCTION DECLARATIONS #
56 #################################
57 */
58 PRIVATE void
59 sc_reset_up(vrna_fold_compound_t *fc,
60 const FLT_OR_DBL *constraints,
61 unsigned int options);
62
63
64 PRIVATE void
65 sc_reset_bp(vrna_fold_compound_t *fc,
66 const FLT_OR_DBL **constraints,
67 unsigned int options);
68
69
70 PRIVATE void
71 sc_add_up(vrna_fold_compound_t *fc,
72 unsigned int i,
73 FLT_OR_DBL energy,
74 unsigned int options);
75
76
77 PRIVATE void
78 sc_add_bp(vrna_fold_compound_t *fc,
79 unsigned int i,
80 unsigned int j,
81 FLT_OR_DBL energy,
82 unsigned int options);
83
84
85 PRIVATE void
86 prepare_sc_up_mfe(vrna_fold_compound_t *fc,
87 unsigned int options);
88
89
90 PRIVATE void
91 prepare_sc_bp_mfe(vrna_fold_compound_t *fc,
92 unsigned int options);
93
94
95 PRIVATE void
96 prepare_sc_up_pf(vrna_fold_compound_t *fc,
97 unsigned int options);
98
99
100 PRIVATE void
101 prepare_sc_bp_pf(vrna_fold_compound_t *fc,
102 unsigned int options);
103
104
105 PRIVATE void
106 prepare_sc_stack_pf(vrna_fold_compound_t *fc);
107
108
109 PRIVATE INLINE void
110 sc_init_up_storage(vrna_sc_t *sc);
111
112
113 PRIVATE INLINE void
114 populate_sc_up_mfe(vrna_fold_compound_t *fc,
115 unsigned int i,
116 unsigned int n);
117
118
119 PRIVATE INLINE void
120 populate_sc_up_pf(vrna_fold_compound_t *fc,
121 unsigned int i,
122 unsigned int n);
123
124
125 PRIVATE INLINE void
126 sc_init_bp_storage(vrna_sc_t *sc);
127
128
129 PRIVATE INLINE void
130 sc_store_bp(vrna_sc_bp_storage_t **container,
131 unsigned int i,
132 unsigned int start,
133 unsigned int end,
134 int e);
135
136
137 PRIVATE INLINE int
138 get_stored_bp_contributions(vrna_sc_bp_storage_t *container,
139 unsigned int j);
140
141
142 PRIVATE INLINE void
143 populate_sc_bp_mfe(vrna_fold_compound_t *fc,
144 unsigned int i,
145 unsigned int n);
146
147
148 PRIVATE INLINE void
149 populate_sc_bp_pf(vrna_fold_compound_t *fc,
150 unsigned int i,
151 unsigned int n);
152
153
154 PRIVATE INLINE void
155 free_sc_up(vrna_sc_t *sc);
156
157
158 PRIVATE INLINE void
159 free_sc_bp(vrna_sc_t *sc);
160
161
162 PRIVATE vrna_sc_t *
163 init_sc_default(unsigned int n);
164
165
166 PRIVATE vrna_sc_t *
167 init_sc_window(unsigned int n);
168
169
170 PRIVATE void
171 nullify(vrna_sc_t *sc);
172
173
174 /*
175 #################################
176 # BEGIN OF FUNCTION DEFINITIONS #
177 #################################
178 */
179 PUBLIC void
vrna_sc_init(vrna_fold_compound_t * fc)180 vrna_sc_init(vrna_fold_compound_t *fc)
181 {
182 unsigned int s, n, N;
183
184 if (fc) {
185 vrna_sc_remove(fc);
186
187 n = fc->length;
188
189 switch (fc->type) {
190 case VRNA_FC_TYPE_SINGLE:
191 fc->sc = init_sc_default(n);
192 break;
193
194 case VRNA_FC_TYPE_COMPARATIVE:
195 N = fc->n_seq;
196 fc->scs = (vrna_sc_t **)vrna_alloc(sizeof(vrna_sc_t *) * (N + 1));
197
198 for (s = 0; s < N; s++)
199 fc->scs[s] = init_sc_default(n);
200
201 break;
202 }
203 }
204 }
205
206
207 PUBLIC void
vrna_sc_init_window(vrna_fold_compound_t * fc)208 vrna_sc_init_window(vrna_fold_compound_t *fc)
209 {
210 unsigned int s, n, N;
211
212 if (fc) {
213 vrna_sc_remove(fc);
214
215 n = fc->length;
216
217 switch (fc->type) {
218 case VRNA_FC_TYPE_SINGLE:
219 fc->sc = init_sc_window(n);
220 break;
221
222 case VRNA_FC_TYPE_COMPARATIVE:
223 N = fc->n_seq;
224 fc->scs = (vrna_sc_t **)vrna_alloc(sizeof(vrna_sc_t *) * (N + 1));
225
226 for (s = 0; s < N; s++)
227 fc->scs[s] = init_sc_window(n);
228
229 break;
230 }
231 }
232 }
233
234
235 PUBLIC void
vrna_sc_prepare(vrna_fold_compound_t * fc,unsigned int options)236 vrna_sc_prepare(vrna_fold_compound_t *fc,
237 unsigned int options)
238 {
239 if (fc) {
240 if (options & VRNA_OPTION_MFE) {
241 prepare_sc_up_mfe(fc, options);
242 prepare_sc_bp_mfe(fc, options);
243 }
244
245 if (options & VRNA_OPTION_PF) {
246 prepare_sc_up_pf(fc, options);
247 prepare_sc_bp_pf(fc, options);
248 prepare_sc_stack_pf(fc);
249 }
250 }
251 }
252
253
254 PUBLIC int
vrna_sc_update(vrna_fold_compound_t * fc,unsigned int i,unsigned int options)255 vrna_sc_update(vrna_fold_compound_t *fc,
256 unsigned int i,
257 unsigned int options)
258 {
259 unsigned int n, maxdist;
260 vrna_sc_t *sc;
261
262 if (fc) {
263 n = fc->length;
264 maxdist = (unsigned int)fc->window_size;
265
266 if (i > n) {
267 vrna_message_warning("vrna_sc_update(): Position %u out of range!"
268 " (Sequence length: %u)",
269 i, n);
270 } else if (i > 0) {
271 maxdist = MIN2(maxdist, n - i + 1);
272
273 switch (fc->type) {
274 case VRNA_FC_TYPE_SINGLE:
275 sc = fc->sc;
276
277 if ((sc) &&
278 (options & VRNA_OPTION_WINDOW)) {
279 /* sliding-window mode, i.e. local structure prediction */
280 if (sc->up_storage) {
281 if (options & VRNA_OPTION_MFE)
282 populate_sc_up_mfe(fc, i, maxdist);
283
284 if (options & VRNA_OPTION_PF)
285 populate_sc_up_pf(fc, i, maxdist);
286 }
287
288 if (sc->bp_storage) {
289 if (options & VRNA_OPTION_MFE)
290 populate_sc_bp_mfe(fc, i, maxdist);
291
292 if (options & VRNA_OPTION_PF)
293 populate_sc_bp_pf(fc, i, maxdist);
294 }
295
296 return 1;
297 }
298
299 break;
300
301 case VRNA_FC_TYPE_COMPARATIVE:
302 break;
303 }
304 }
305 }
306
307 return 0;
308 }
309
310
311 PUBLIC void
vrna_sc_remove(vrna_fold_compound_t * fc)312 vrna_sc_remove(vrna_fold_compound_t *fc)
313 {
314 unsigned int s;
315
316 if (fc) {
317 switch (fc->type) {
318 case VRNA_FC_TYPE_SINGLE:
319 vrna_sc_free(fc->sc);
320 fc->sc = NULL;
321 break;
322
323 case VRNA_FC_TYPE_COMPARATIVE:
324 if (fc->scs) {
325 for (s = 0; s < fc->n_seq; s++)
326 vrna_sc_free(fc->scs[s]);
327 free(fc->scs);
328 }
329
330 fc->scs = NULL;
331 break;
332 }
333 }
334 }
335
336
337 PUBLIC void
vrna_sc_free(vrna_sc_t * sc)338 vrna_sc_free(vrna_sc_t *sc)
339 {
340 if (sc) {
341 free_sc_up(sc);
342 free_sc_bp(sc);
343
344 free(sc->energy_stack);
345 free(sc->exp_energy_stack);
346
347 if (sc->free_data)
348 sc->free_data(sc->data);
349
350 free(sc);
351 }
352 }
353
354
355 PUBLIC int
vrna_sc_set_bp(vrna_fold_compound_t * fc,const FLT_OR_DBL ** constraints,unsigned int options)356 vrna_sc_set_bp(vrna_fold_compound_t *fc,
357 const FLT_OR_DBL **constraints,
358 unsigned int options)
359 {
360 if (fc && (fc->type == VRNA_FC_TYPE_SINGLE)) {
361 sc_reset_bp(fc, constraints, options);
362
363 if (options & VRNA_OPTION_MFE)
364 prepare_sc_bp_mfe(fc, options);
365
366 if (options & VRNA_OPTION_PF) /* prepare Boltzmann factors for the BP soft constraints */
367 prepare_sc_bp_pf(fc, options);
368
369 return 1;
370 }
371
372 return 0;
373 }
374
375
376 PUBLIC int
vrna_sc_add_bp(vrna_fold_compound_t * fc,int i,int j,FLT_OR_DBL energy,unsigned int options)377 vrna_sc_add_bp(vrna_fold_compound_t *fc,
378 int i,
379 int j,
380 FLT_OR_DBL energy,
381 unsigned int options)
382 {
383 if (fc && (fc->type == VRNA_FC_TYPE_SINGLE)) {
384 if ((i < 1) ||
385 (i > fc->length) ||
386 (j < i) ||
387 (j > fc->length)) {
388 vrna_message_warning("vrna_sc_add_bp(): Base pair (%d, %d) out of range!"
389 " (Sequence length: %d)",
390 i, j, fc->length);
391 } else {
392 sc_add_bp(fc, (unsigned int)i, (unsigned int)j, energy, options);
393
394 if (options & VRNA_OPTION_MFE)
395 prepare_sc_bp_mfe(fc, options);
396
397 if (options & VRNA_OPTION_PF) /* prepare Boltzmann factors for the BP soft constraints */
398 prepare_sc_bp_pf(fc, options);
399
400 return 1;
401 }
402 }
403
404 return 0;
405 }
406
407
408 PUBLIC int
vrna_sc_set_up(vrna_fold_compound_t * fc,const FLT_OR_DBL * constraints,unsigned int options)409 vrna_sc_set_up(vrna_fold_compound_t *fc,
410 const FLT_OR_DBL *constraints,
411 unsigned int options)
412 {
413 if (fc && (fc->type == VRNA_FC_TYPE_SINGLE)) {
414 sc_reset_up(fc, constraints, options);
415
416 if (options & VRNA_OPTION_MFE)
417 prepare_sc_up_mfe(fc, options);
418
419 if (options & VRNA_OPTION_PF)
420 prepare_sc_up_pf(fc, options);
421
422 return 1;
423 }
424
425 return 0;
426 }
427
428
429 PUBLIC int
vrna_sc_add_up(vrna_fold_compound_t * fc,int i,FLT_OR_DBL energy,unsigned int options)430 vrna_sc_add_up(vrna_fold_compound_t *fc,
431 int i,
432 FLT_OR_DBL energy,
433 unsigned int options)
434 {
435 if (fc && (fc->type == VRNA_FC_TYPE_SINGLE)) {
436 if ((i < 1) || (i > fc->length)) {
437 vrna_message_warning("vrna_sc_add_up(): Nucleotide position %d out of range!"
438 " (Sequence length: %d)",
439 i, fc->length);
440 } else {
441 sc_add_up(fc, (unsigned int)i, energy, options);
442
443 if (options & VRNA_OPTION_MFE)
444 prepare_sc_up_mfe(fc, options);
445
446 if (options & VRNA_OPTION_PF)
447 prepare_sc_up_pf(fc, options);
448
449 return 1;
450 }
451 }
452
453 return 0;
454 }
455
456
457 PUBLIC int
vrna_sc_set_stack(vrna_fold_compound_t * fc,const FLT_OR_DBL * constraints,unsigned int options)458 vrna_sc_set_stack(vrna_fold_compound_t *fc,
459 const FLT_OR_DBL *constraints,
460 unsigned int options)
461 {
462 unsigned int i;
463
464 if ((fc) &&
465 (constraints) &&
466 (fc->type == VRNA_FC_TYPE_SINGLE)) {
467 if (!fc->sc) {
468 if (options & VRNA_OPTION_WINDOW)
469 vrna_sc_init_window(fc);
470 else
471 vrna_sc_init(fc);
472 }
473
474 free(fc->sc->energy_stack);
475 fc->sc->energy_stack = (int *)vrna_alloc(sizeof(int) * (fc->length + 1));
476
477 for (i = 1; i <= fc->length; ++i)
478 fc->sc->energy_stack[i] = (int)roundf(constraints[i] * 100.);
479
480 return 1;
481 }
482
483 return 0;
484 }
485
486
487 PUBLIC int
vrna_sc_set_stack_comparative(vrna_fold_compound_t * fc,const FLT_OR_DBL ** constraints,unsigned int options)488 vrna_sc_set_stack_comparative(vrna_fold_compound_t *fc,
489 const FLT_OR_DBL **constraints,
490 unsigned int options)
491 {
492 unsigned int i, s;
493
494 if ((fc) &&
495 (constraints) &&
496 (fc->type == VRNA_FC_TYPE_COMPARATIVE)) {
497 if (!fc->scs) {
498 if (options & VRNA_OPTION_WINDOW)
499 vrna_sc_init_window(fc);
500 else
501 vrna_sc_init(fc);
502 }
503
504 for (s = 0; s < fc->n_seq; s++) {
505 free(fc->scs[s]->energy_stack);
506 fc->scs[s]->energy_stack = NULL;
507
508 if (constraints[s]) {
509 fc->scs[s]->energy_stack = (int *)vrna_alloc(sizeof(int) * (fc->length + 1));
510
511 for (i = 1; i <= fc->length; ++i)
512 fc->scs[s]->energy_stack[i] = (int)roundf(constraints[s][i] * 100.);
513 }
514 }
515
516 return 1;
517 }
518
519 return 0;
520 }
521
522
523 PUBLIC int
vrna_sc_add_stack(vrna_fold_compound_t * fc,int i,FLT_OR_DBL energy,unsigned int options)524 vrna_sc_add_stack(vrna_fold_compound_t *fc,
525 int i,
526 FLT_OR_DBL energy,
527 unsigned int options)
528 {
529 if ((fc) &&
530 (fc->type == VRNA_FC_TYPE_SINGLE)) {
531 if ((i < 1) || (i > fc->length)) {
532 vrna_message_warning("vrna_sc_add_stack*(): Nucleotide position %d out of range!"
533 " (Sequence length: %d)",
534 i, fc->length);
535 } else {
536 if (!fc->sc) {
537 if (options & VRNA_OPTION_WINDOW)
538 vrna_sc_init_window(fc);
539 else
540 vrna_sc_init(fc);
541 }
542
543 if (!fc->sc->energy_stack)
544 fc->sc->energy_stack = (int *)vrna_alloc(sizeof(int) * (fc->length + 1));
545
546 fc->sc->energy_stack[i] += (int)roundf(energy * 100.);
547
548 return 1;
549 }
550 }
551
552 return 0;
553 }
554
555
556 PUBLIC int
vrna_sc_add_stack_comparative(vrna_fold_compound_t * fc,int i,const FLT_OR_DBL * energies,unsigned int options)557 vrna_sc_add_stack_comparative(vrna_fold_compound_t *fc,
558 int i,
559 const FLT_OR_DBL *energies,
560 unsigned int options)
561 {
562 unsigned int s;
563
564 if ((fc) &&
565 (fc->type == VRNA_FC_TYPE_COMPARATIVE)) {
566 if ((i < 1) || (i > fc->length)) {
567 vrna_message_warning("vrna_sc_add_stack*(): Nucleotide position %d out of range!"
568 " (Alignment length: %d)",
569 i, fc->length);
570 } else {
571 if (!fc->scs) {
572 if (options & VRNA_OPTION_WINDOW)
573 vrna_sc_init_window(fc);
574 else
575 vrna_sc_init(fc);
576 }
577
578 for (s = 0; s < fc->n_seq; s++) {
579 if (!fc->scs[s]->energy_stack)
580 fc->scs[s]->energy_stack = (int *)vrna_alloc(sizeof(int) * (fc->length + 1));
581
582 fc->scs[s]->energy_stack[i] += (int)roundf(energies[s] * 100.);
583 }
584
585 return 1;
586 }
587 }
588
589 return 0;
590 }
591
592
593 PUBLIC int
vrna_sc_add_data(vrna_fold_compound_t * fc,void * data,vrna_callback_free_auxdata * free_data)594 vrna_sc_add_data(vrna_fold_compound_t *fc,
595 void *data,
596 vrna_callback_free_auxdata *free_data)
597 {
598 if ((fc) &&
599 (fc->type == VRNA_FC_TYPE_SINGLE)) {
600 if (!fc->sc)
601 vrna_sc_init(fc);
602
603 fc->sc->data = data;
604 fc->sc->free_data = free_data;
605 return 1;
606 }
607
608 return 0;
609 }
610
611
612 PUBLIC int
vrna_sc_add_data_comparative(vrna_fold_compound_t * fc,void ** data,vrna_callback_free_auxdata ** free_data)613 vrna_sc_add_data_comparative(vrna_fold_compound_t *fc,
614 void **data,
615 vrna_callback_free_auxdata **free_data)
616 {
617 unsigned int s;
618
619 if ((fc) &&
620 (fc->type == VRNA_FC_TYPE_COMPARATIVE)) {
621 if (!fc->scs)
622 vrna_sc_init(fc);
623
624 if (data)
625 for (s = 0; s < fc->n_seq; s++)
626 fc->scs[s]->data = data[s];
627
628 if (free_data)
629 for (s = 0; s < fc->n_seq; s++)
630 fc->scs[s]->free_data = free_data[s];
631
632 return 1;
633 }
634
635 return 0;
636 }
637
638
639 PUBLIC int
vrna_sc_add_f(vrna_fold_compound_t * fc,vrna_callback_sc_energy * f)640 vrna_sc_add_f(vrna_fold_compound_t *fc,
641 vrna_callback_sc_energy *f)
642 {
643 if ((fc) &&
644 (f) &&
645 (fc->type == VRNA_FC_TYPE_SINGLE)) {
646 if (!fc->sc)
647 vrna_sc_init(fc);
648
649 fc->sc->f = f;
650 return 1;
651 }
652
653 return 0;
654 }
655
656
657 PUBLIC int
vrna_sc_add_f_comparative(vrna_fold_compound_t * fc,vrna_callback_sc_energy ** f)658 vrna_sc_add_f_comparative(vrna_fold_compound_t *fc,
659 vrna_callback_sc_energy **f)
660 {
661 unsigned int s;
662
663 if ((fc) &&
664 (f) &&
665 (fc->type == VRNA_FC_TYPE_COMPARATIVE)) {
666 if (!fc->scs)
667 vrna_sc_init(fc);
668
669 for (s = 0; s < fc->n_seq; s++)
670 fc->scs[s]->f = f[s];
671
672 return 1;
673 }
674
675 return 0;
676 }
677
678
679 PUBLIC int
vrna_sc_add_bt(vrna_fold_compound_t * fc,vrna_callback_sc_backtrack * f)680 vrna_sc_add_bt(vrna_fold_compound_t *fc,
681 vrna_callback_sc_backtrack *f)
682 {
683 if ((fc) &&
684 (f) &&
685 (fc->type == VRNA_FC_TYPE_SINGLE)) {
686 if (!fc->sc)
687 vrna_sc_init(fc);
688
689 fc->sc->bt = f;
690 return 1;
691 }
692
693 return 0;
694 }
695
696
697 PUBLIC int
vrna_sc_add_exp_f(vrna_fold_compound_t * fc,vrna_callback_sc_exp_energy * exp_f)698 vrna_sc_add_exp_f(vrna_fold_compound_t *fc,
699 vrna_callback_sc_exp_energy *exp_f)
700 {
701 if ((fc) &&
702 (exp_f) &&
703 (fc->type == VRNA_FC_TYPE_SINGLE)) {
704 if (!fc->sc)
705 vrna_sc_init(fc);
706
707 fc->sc->exp_f = exp_f;
708 return 1;
709 }
710
711 return 0;
712 }
713
714
715 PUBLIC int
vrna_sc_add_exp_f_comparative(vrna_fold_compound_t * fc,vrna_callback_sc_exp_energy ** exp_f)716 vrna_sc_add_exp_f_comparative(vrna_fold_compound_t *fc,
717 vrna_callback_sc_exp_energy **exp_f)
718 {
719 unsigned int s;
720
721 if ((fc) &&
722 (exp_f) &&
723 (fc->type == VRNA_FC_TYPE_COMPARATIVE)) {
724 if (!fc->scs)
725 vrna_sc_init(fc);
726
727 for (s = 0; s < fc->n_seq; s++)
728 fc->scs[s]->exp_f = exp_f[s];
729
730 return 1;
731 }
732
733 return 0;
734 }
735
736
737 /*
738 #####################################
739 # BEGIN OF STATIC HELPER FUNCTIONS #
740 #####################################
741 */
742 PRIVATE INLINE void
sc_init_up_storage(vrna_sc_t * sc)743 sc_init_up_storage(vrna_sc_t *sc)
744 {
745 if (!sc->up_storage)
746 sc->up_storage = (int *)vrna_alloc(sizeof(int) * (sc->n + 2));
747 }
748
749
750 /* pupulate sc->energy_up array at position i from sc->up_storage data */
751 PRIVATE INLINE void
populate_sc_up_mfe(vrna_fold_compound_t * fc,unsigned int i,unsigned int n)752 populate_sc_up_mfe(vrna_fold_compound_t *fc,
753 unsigned int i,
754 unsigned int n)
755 {
756 unsigned int j;
757 vrna_sc_t *sc = fc->sc;
758
759 sc->energy_up[i][0] = 0;
760 for (j = 1; j <= n; j++)
761 sc->energy_up[i][j] = sc->energy_up[i][j - 1]
762 + sc->up_storage[i + j - 1];
763 }
764
765
766 PRIVATE INLINE void
populate_sc_up_pf(vrna_fold_compound_t * fc,unsigned int i,unsigned int n)767 populate_sc_up_pf(vrna_fold_compound_t *fc,
768 unsigned int i,
769 unsigned int n)
770 {
771 unsigned int j;
772 double GT, kT;
773 vrna_sc_t *sc = fc->sc;
774
775 kT = fc->exp_params->kT;
776
777 sc->exp_energy_up[i][0] = 1.;
778
779 for (j = 1; j <= n; j++) {
780 GT = (double)sc->up_storage[i + j - 1] * 10.; /* convert deka-cal/mol to cal/mol */
781 sc->exp_energy_up[i][j] = sc->exp_energy_up[i][j - 1]
782 * (FLT_OR_DBL)exp(-GT / kT);
783 }
784 }
785
786
787 PRIVATE INLINE void
sc_init_bp_storage(vrna_sc_t * sc)788 sc_init_bp_storage(vrna_sc_t *sc)
789 {
790 unsigned int i;
791
792 if (sc->bp_storage == NULL) {
793 sc->bp_storage = (vrna_sc_bp_storage_t **)vrna_alloc(
794 sizeof(vrna_sc_bp_storage_t *) * (sc->n + 2));
795
796 for (i = 1; i <= sc->n; i++)
797 /* by default we do not change energy contributions for any base pairs */
798 sc->bp_storage[i] = NULL;
799 }
800 }
801
802
803 PRIVATE INLINE void
sc_store_bp(vrna_sc_bp_storage_t ** container,unsigned int i,unsigned int start,unsigned int end,int e)804 sc_store_bp(vrna_sc_bp_storage_t **container,
805 unsigned int i,
806 unsigned int start,
807 unsigned int end,
808 int e)
809 {
810 unsigned int size, cnt = 0;
811
812 if (!container[i]) {
813 container[i] = (vrna_sc_bp_storage_t *)vrna_alloc(sizeof(vrna_sc_bp_storage_t) * 2);
814 } else {
815 /* find out total size of container */
816 for (size = 0; container[i][size].interval_start != 0; size++);
817
818 /* find position where we want to insert the new constraint */
819 for (cnt = 0; cnt < size; cnt++) {
820 if (container[i][cnt].interval_start > start)
821 break; /* want to insert before current constraint */
822
823 if (container[i][cnt].interval_end < end)
824 continue; /* want to insert after current constraint */
825 }
826 /* increase memory for bp constraints */
827 container[i] = (vrna_sc_bp_storage_t *)vrna_realloc(container[i],
828 sizeof(vrna_sc_bp_storage_t) * (size + 2));
829 /* shift trailing constraints by 1 entry */
830 memmove(container[i] + cnt + 1, container[i] + cnt,
831 sizeof(vrna_sc_bp_storage_t) * (size - cnt + 1));
832 }
833
834 container[i][cnt].interval_start = start;
835 container[i][cnt].interval_end = end;
836 container[i][cnt].e = e;
837 }
838
839
840 PRIVATE INLINE int
get_stored_bp_contributions(vrna_sc_bp_storage_t * container,unsigned int j)841 get_stored_bp_contributions(vrna_sc_bp_storage_t *container,
842 unsigned int j)
843 {
844 unsigned int cnt;
845 int e;
846
847 e = 0;
848
849 /* go through list of constraints for current position i */
850 for (cnt = 0; container[cnt].interval_start != 0; cnt++) {
851 if (container[cnt].interval_start > j)
852 break; /* only constraints for pairs (i,q) with q > j left */
853
854 if (container[cnt].interval_end < j)
855 continue; /* constraint for pairs (i,q) with q < j */
856
857 /* constraint has interval [p,q] with p <= j <= q */
858 e += container[cnt].e;
859 }
860
861 return e;
862 }
863
864
865 PRIVATE INLINE void
populate_sc_bp_mfe(vrna_fold_compound_t * fc,unsigned int i,unsigned int maxdist)866 populate_sc_bp_mfe(vrna_fold_compound_t *fc,
867 unsigned int i,
868 unsigned int maxdist)
869 {
870 unsigned int j, k, turn, n;
871 int e, *idx;
872 vrna_sc_t *sc;
873
874 n = fc->length;
875 turn = fc->params->model_details.min_loop_size;
876 sc = fc->sc;
877 idx = fc->jindx;
878
879 if (sc->bp_storage[i]) {
880 for (k = turn + 1; k < maxdist; k++) {
881 j = i + k;
882
883 if (j > n)
884 break;
885
886 e = get_stored_bp_contributions(sc->bp_storage[i], j);
887
888 switch (sc->type) {
889 case VRNA_SC_DEFAULT:
890 sc->energy_bp[idx[j] + i] = e;
891 break;
892
893 case VRNA_SC_WINDOW:
894 sc->energy_bp_local[i][j - i] = e;
895 break;
896 }
897 }
898 } else {
899 for (k = turn + 1; k < maxdist; k++) {
900 j = i + k;
901 if (j > n)
902 break;
903
904 switch (sc->type) {
905 case VRNA_SC_DEFAULT:
906 sc->energy_bp[idx[j] + i] = 0;
907 break;
908
909 case VRNA_SC_WINDOW:
910 sc->energy_bp_local[i][j - i] = 0;
911 break;
912 }
913 }
914 }
915 }
916
917
918 PRIVATE INLINE void
populate_sc_bp_pf(vrna_fold_compound_t * fc,unsigned int i,unsigned int maxdist)919 populate_sc_bp_pf(vrna_fold_compound_t *fc,
920 unsigned int i,
921 unsigned int maxdist)
922 {
923 unsigned int j, k, turn, n;
924 int e, *idx;
925 FLT_OR_DBL q;
926 double GT, kT;
927 vrna_exp_param_t *exp_params;
928 vrna_sc_t *sc;
929
930 n = fc->length;
931 exp_params = fc->exp_params;
932 kT = exp_params->kT;
933 turn = exp_params->model_details.min_loop_size;
934 sc = fc->sc;
935 idx = fc->jindx;
936
937 if (sc->bp_storage[i]) {
938 for (k = turn + 1; k < maxdist; k++) {
939 j = i + k;
940
941 if (j > n)
942 break;
943
944 e = get_stored_bp_contributions(sc->bp_storage[i], j);
945
946 GT = e * 10.;
947 q = (FLT_OR_DBL)exp(-GT / kT);
948
949 switch (sc->type) {
950 case VRNA_SC_DEFAULT:
951 sc->exp_energy_bp[idx[j] + i] = q;
952 break;
953
954 case VRNA_SC_WINDOW:
955 sc->exp_energy_bp_local[i][j - i] = q;
956 break;
957 }
958 }
959 } else {
960 for (k = turn + 1; k < maxdist; k++) {
961 j = i + k;
962 if (j > n)
963 break;
964
965 switch (sc->type) {
966 case VRNA_SC_DEFAULT:
967 sc->exp_energy_bp[idx[j] + i] = 1.;
968 break;
969
970 case VRNA_SC_WINDOW:
971 sc->exp_energy_bp_local[i][j - i] = 1.;
972 break;
973 }
974 }
975 }
976 }
977
978
979 PRIVATE void
sc_add_bp(vrna_fold_compound_t * fc,unsigned int i,unsigned int j,FLT_OR_DBL energy,unsigned int options)980 sc_add_bp(vrna_fold_compound_t *fc,
981 unsigned int i,
982 unsigned int j,
983 FLT_OR_DBL energy,
984 unsigned int options)
985 {
986 vrna_sc_t *sc;
987
988 if ((options & VRNA_OPTION_WINDOW) && (!fc->sc))
989 vrna_sc_init_window(fc);
990 else if (!fc->sc)
991 vrna_sc_init(fc);
992
993 sc = fc->sc;
994 sc_init_bp_storage(sc);
995 sc_store_bp(sc->bp_storage, i, j, j, (int)roundf(energy * 100.));
996 sc->state |= STATE_DIRTY_BP_MFE | STATE_DIRTY_BP_PF;
997 }
998
999
1000 PRIVATE INLINE void
free_sc_up(vrna_sc_t * sc)1001 free_sc_up(vrna_sc_t *sc)
1002 {
1003 unsigned int i;
1004
1005 free(sc->up_storage);
1006
1007 sc->up_storage = NULL;
1008
1009 if (sc->type == VRNA_SC_DEFAULT) {
1010 if (sc->energy_up)
1011 for (i = 0; i <= sc->n + 1; i++)
1012 free(sc->energy_up[i]);
1013
1014 if (sc->exp_energy_up)
1015 for (i = 0; i <= sc->n + 1; i++)
1016 free(sc->exp_energy_up[i]);
1017 }
1018
1019 free(sc->energy_up);
1020 sc->energy_up = NULL;
1021
1022 free(sc->exp_energy_up);
1023 sc->exp_energy_up = NULL;
1024
1025 sc->state &= ~(STATE_DIRTY_UP_MFE | STATE_DIRTY_UP_PF);
1026 }
1027
1028
1029 PRIVATE INLINE void
free_sc_bp(vrna_sc_t * sc)1030 free_sc_bp(vrna_sc_t *sc)
1031 {
1032 unsigned int i;
1033
1034 if (sc->bp_storage) {
1035 for (i = 1; i <= sc->n; i++)
1036 free(sc->bp_storage[i]);
1037 free(sc->bp_storage);
1038 sc->bp_storage = NULL;
1039 }
1040
1041 switch (sc->type) {
1042 case VRNA_SC_DEFAULT:
1043 free(sc->energy_bp);
1044 sc->energy_bp = NULL;
1045
1046 free(sc->exp_energy_bp);
1047 sc->energy_bp = NULL;
1048
1049 break;
1050
1051 case VRNA_SC_WINDOW:
1052 free(sc->energy_bp_local);
1053 sc->energy_bp_local = NULL;
1054
1055 free(sc->exp_energy_bp_local);
1056 sc->exp_energy_bp_local = NULL;
1057
1058 break;
1059 }
1060 sc->state &= ~(STATE_DIRTY_BP_MFE | STATE_DIRTY_BP_PF);
1061 }
1062
1063
1064 PRIVATE void
sc_reset_up(vrna_fold_compound_t * fc,const FLT_OR_DBL * constraints,unsigned int options)1065 sc_reset_up(vrna_fold_compound_t *fc,
1066 const FLT_OR_DBL *constraints,
1067 unsigned int options)
1068 {
1069 unsigned int i, n;
1070 vrna_sc_t *sc;
1071
1072 n = fc->length;
1073
1074 if (!fc->sc) {
1075 if (options & VRNA_OPTION_WINDOW)
1076 vrna_sc_init_window(fc);
1077 else
1078 vrna_sc_init(fc);
1079 }
1080
1081 sc = fc->sc;
1082
1083 if (constraints) {
1084 free_sc_up(sc);
1085
1086 /* initialize container for unpaired probabilities */
1087 sc_init_up_storage(sc);
1088
1089 /* add contributions to storage container */
1090 for (i = 1; i <= n; i++)
1091 sc->up_storage[i] = (int)roundf(constraints[i] * 100.); /* convert to 10kal/mol */
1092
1093 sc->state |= STATE_DIRTY_UP_MFE | STATE_DIRTY_UP_PF;
1094 } else {
1095 free_sc_up(sc);
1096 }
1097 }
1098
1099
1100 PRIVATE void
sc_reset_bp(vrna_fold_compound_t * fc,const FLT_OR_DBL ** constraints,unsigned int options)1101 sc_reset_bp(vrna_fold_compound_t *fc,
1102 const FLT_OR_DBL **constraints,
1103 unsigned int options)
1104 {
1105 unsigned int i, j, n;
1106 vrna_sc_t *sc;
1107
1108 n = fc->length;
1109
1110 if (!fc->sc) {
1111 if (options & VRNA_OPTION_WINDOW)
1112 vrna_sc_init_window(fc);
1113 else
1114 vrna_sc_init(fc);
1115 }
1116
1117 sc = fc->sc;
1118
1119 if (constraints) {
1120 free_sc_bp(sc);
1121
1122 /* initialize container for base pair constraints */
1123 sc_init_bp_storage(sc);
1124
1125 /* add contributions to storage container */
1126 for (i = 1; i < n; i++)
1127 for (j = i + 1; j <= n; j++)
1128 sc_store_bp(sc->bp_storage, i, j, j, (int)roundf(constraints[i][j] * 100.));
1129
1130 sc->state |= STATE_DIRTY_BP_MFE | STATE_DIRTY_BP_PF;
1131 } else {
1132 free_sc_bp(sc);
1133 }
1134 }
1135
1136
1137 PRIVATE void
sc_add_up(vrna_fold_compound_t * fc,unsigned int i,FLT_OR_DBL energy,unsigned int options)1138 sc_add_up(vrna_fold_compound_t *fc,
1139 unsigned int i,
1140 FLT_OR_DBL energy,
1141 unsigned int options)
1142 {
1143 vrna_sc_t *sc;
1144
1145 if ((options & VRNA_OPTION_WINDOW) && (!fc->sc))
1146 vrna_sc_init_window(fc);
1147 else if (!fc->sc)
1148 vrna_sc_init(fc);
1149
1150 sc = fc->sc;
1151 sc_init_up_storage(sc);
1152 sc->up_storage[i] += (int)roundf(energy * 100.);
1153 sc->state |= STATE_DIRTY_UP_MFE | STATE_DIRTY_UP_PF;
1154 }
1155
1156
1157 /* populate sc->energy_up arrays for usage in MFE computations */
1158 PRIVATE void
prepare_sc_up_mfe(vrna_fold_compound_t * fc,unsigned int options)1159 prepare_sc_up_mfe(vrna_fold_compound_t *fc,
1160 unsigned int options)
1161 {
1162 unsigned int i, n;
1163 vrna_sc_t *sc;
1164
1165 n = fc->length;
1166 switch (fc->type) {
1167 case VRNA_FC_TYPE_SINGLE:
1168 sc = fc->sc;
1169 if (sc) {
1170 /* prepare sc for unpaired nucleotides only if we actually have some to apply */
1171 if (sc->up_storage) {
1172 if (sc->state & STATE_DIRTY_UP_MFE) {
1173 /* allocate memory such that we can access the soft constraint
1174 * energies of a subsequence of length j starting at position i
1175 * via sc->energy_up[i][j]
1176 */
1177 sc->energy_up = (int **)vrna_realloc(sc->energy_up, sizeof(int *) * (n + 2));
1178
1179 if (options & VRNA_OPTION_WINDOW) {
1180 /*
1181 * simply init with NULL pointers, since the sliding-window implementation must take
1182 * care of allocating the required memory and filling in appropriate energy contributions
1183 */
1184 for (i = 0; i <= n + 1; i++)
1185 sc->energy_up[i] = NULL;
1186 } else {
1187 for (i = 1; i <= n; i++)
1188 sc->energy_up[i] = (int *)vrna_realloc(sc->energy_up[i], sizeof(int) * (n - i + 2));
1189
1190
1191 sc->energy_up[0] = (int *)vrna_realloc(sc->energy_up[0], sizeof(int));
1192 sc->energy_up[n + 1] = (int *)vrna_realloc(sc->energy_up[n + 1], sizeof(int));
1193
1194 /* now add soft constraints as stored in container for unpaired sc */
1195 for (i = 1; i <= n; i++)
1196 populate_sc_up_mfe(fc, i, (n - i + 1));
1197
1198 sc->energy_up[0][0] = 0;
1199 sc->energy_up[n + 1][0] = 0;
1200 }
1201
1202 sc->state &= ~STATE_DIRTY_UP_MFE;
1203 }
1204 } else if (sc->energy_up) {
1205 /* remove any unpaired sc if storage container is empty */
1206 free_sc_up(sc);
1207 }
1208 }
1209
1210 break;
1211
1212 case VRNA_FC_TYPE_COMPARATIVE:
1213 break; /* do nothing for now */
1214 }
1215 }
1216
1217
1218 /* populate sc->exp_energy_up arrays for usage in partition function computations */
1219 PRIVATE void
prepare_sc_up_pf(vrna_fold_compound_t * fc,unsigned int options)1220 prepare_sc_up_pf(vrna_fold_compound_t *fc,
1221 unsigned int options)
1222 {
1223 unsigned int i, n;
1224 vrna_sc_t *sc;
1225
1226 n = fc->length;
1227
1228 switch (fc->type) {
1229 case VRNA_FC_TYPE_SINGLE:
1230 sc = fc->sc;
1231 if (sc) {
1232 /* prepare sc for unpaired nucleotides only if we actually have some to apply */
1233 if (sc->up_storage) {
1234 if (sc->state & STATE_DIRTY_UP_PF) {
1235 /* allocate memory such that we can access the soft constraint
1236 * energies of a subsequence of length j starting at position i
1237 * via sc->exp_energy_up[i][j]
1238 */
1239 sc->exp_energy_up = (FLT_OR_DBL **)vrna_realloc(sc->exp_energy_up,
1240 sizeof(FLT_OR_DBL *) * (n + 2));
1241
1242 if (options & VRNA_OPTION_WINDOW) {
1243 /*
1244 * simply init with NULL pointers, since the sliding-window implementation must take
1245 * care of allocating the required memory and filling in appropriate Boltzmann factors
1246 */
1247 for (i = 0; i <= n + 1; i++)
1248 sc->exp_energy_up[i] = NULL;
1249 } else {
1250 for (i = 1; i <= n; i++)
1251 sc->exp_energy_up[i] =
1252 (FLT_OR_DBL *)vrna_realloc(sc->exp_energy_up[i],
1253 sizeof(FLT_OR_DBL) * (n - i + 2));
1254
1255 sc->exp_energy_up[0] = (FLT_OR_DBL *)vrna_realloc(sc->exp_energy_up[0], sizeof(FLT_OR_DBL));
1256 sc->exp_energy_up[n + 1] = (FLT_OR_DBL *)vrna_realloc(sc->exp_energy_up[n + 1], sizeof(FLT_OR_DBL));
1257
1258 for (i = 1; i <= n; i++)
1259 populate_sc_up_pf(fc, i, (n - i + 1));
1260
1261 sc->exp_energy_up[0][0] = 1.;
1262 sc->exp_energy_up[n + 1][0] = 1.;
1263 }
1264
1265 sc->state &= ~STATE_DIRTY_UP_PF;
1266 }
1267 }
1268 }
1269
1270 break;
1271
1272 case VRNA_FC_TYPE_COMPARATIVE:
1273 break; /* do nothing for now */
1274 }
1275 }
1276
1277
1278 /* populate sc->energy_bp arrays for usage in MFE computations */
1279 PRIVATE void
prepare_sc_bp_mfe(vrna_fold_compound_t * fc,unsigned int options)1280 prepare_sc_bp_mfe(vrna_fold_compound_t *fc,
1281 unsigned int options)
1282 {
1283 unsigned int i, n;
1284 vrna_sc_t *sc;
1285
1286 n = fc->length;
1287
1288 switch (fc->type) {
1289 case VRNA_FC_TYPE_SINGLE:
1290 sc = fc->sc;
1291 if (sc) {
1292 /* prepare sc for base paired positions only if we actually have some to apply */
1293 if (sc->bp_storage) {
1294 if (sc->state & STATE_DIRTY_BP_MFE) {
1295 if (options & VRNA_OPTION_WINDOW) {
1296 sc->energy_bp_local =
1297 (int **)vrna_realloc(sc->energy_bp_local, sizeof(int *) * (n + 2));
1298 } else {
1299 sc->energy_bp =
1300 (int *)vrna_realloc(sc->energy_bp, sizeof(int) * (((n + 1) * (n + 2)) / 2));
1301
1302 for (i = 1; i < n; i++)
1303 populate_sc_bp_mfe(fc, i, n);
1304 }
1305
1306 sc->state &= ~STATE_DIRTY_BP_MFE;
1307 }
1308 } else {
1309 free_sc_bp(sc);
1310 }
1311 }
1312
1313 break;
1314
1315 case VRNA_FC_TYPE_COMPARATIVE:
1316 break; /* do nothing for now */
1317 }
1318 }
1319
1320
1321 /* populate sc->exp_energy_bp arrays for usage in partition function computations */
1322 PRIVATE void
prepare_sc_bp_pf(vrna_fold_compound_t * fc,unsigned int options)1323 prepare_sc_bp_pf(vrna_fold_compound_t *fc,
1324 unsigned int options)
1325 {
1326 unsigned int i, n;
1327 vrna_sc_t *sc;
1328
1329 n = fc->length;
1330
1331 switch (fc->type) {
1332 case VRNA_FC_TYPE_SINGLE:
1333 sc = fc->sc;
1334 if (sc) {
1335 /* prepare sc for base paired positions only if we actually have some to apply */
1336 if (sc->bp_storage) {
1337 if (sc->state & STATE_DIRTY_BP_PF) {
1338 if (options & VRNA_OPTION_WINDOW) {
1339 sc->exp_energy_bp_local =
1340 (FLT_OR_DBL **)vrna_realloc(sc->exp_energy_bp_local,
1341 sizeof(FLT_OR_DBL *) * (n + 2));
1342 } else {
1343 sc->exp_energy_bp =
1344 (FLT_OR_DBL *)vrna_realloc(sc->exp_energy_bp,
1345 sizeof(FLT_OR_DBL) * (((n + 1) * (n + 2)) / 2));
1346
1347 for (i = 1; i < n; i++)
1348 populate_sc_bp_pf(fc, i, n);
1349 }
1350
1351 sc->state &= ~STATE_DIRTY_BP_PF;
1352 }
1353 }
1354 }
1355
1356 break;
1357
1358 case VRNA_FC_TYPE_COMPARATIVE:
1359 break; /* do nothing for now */
1360 }
1361 }
1362
1363
1364 PRIVATE void
prepare_sc_stack_pf(vrna_fold_compound_t * fc)1365 prepare_sc_stack_pf(vrna_fold_compound_t *fc)
1366 {
1367 unsigned int s, n_seq;
1368 int i;
1369 vrna_sc_t *sc, **scs;
1370
1371 switch (fc->type) {
1372 case VRNA_FC_TYPE_SINGLE:
1373 sc = fc->sc;
1374 if ((sc) && (sc->energy_stack)) {
1375 if (!sc->exp_energy_stack) {
1376 sc->exp_energy_stack = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (fc->length + 1));
1377 for (i = 0; i <= fc->length; ++i)
1378 sc->exp_energy_stack[i] = 1.;
1379 }
1380
1381 for (i = 1; i <= fc->length; ++i)
1382 sc->exp_energy_stack[i] = (FLT_OR_DBL)exp(
1383 -(sc->energy_stack[i] * 10.) / fc->exp_params->kT);
1384 }
1385
1386 break;
1387
1388 case VRNA_FC_TYPE_COMPARATIVE:
1389 scs = fc->scs;
1390 n_seq = fc->n_seq;
1391 if (scs) {
1392 for (s = 0; s < n_seq; s++) {
1393 if (scs[s] && scs[s]->energy_stack) {
1394 if (!scs[s]->exp_energy_stack) {
1395 scs[s]->exp_energy_stack =
1396 (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (fc->a2s[s][fc->length] + 1));
1397 for (i = 0; i <= fc->a2s[s][fc->length]; i++)
1398 scs[s]->exp_energy_stack[i] = 1.;
1399 }
1400
1401 for (i = 1; i <= fc->a2s[s][fc->length]; ++i)
1402 scs[s]->exp_energy_stack[i] = (FLT_OR_DBL)exp(
1403 -(scs[s]->energy_stack[i] * 10.) / fc->exp_params->kT);
1404 }
1405 }
1406 }
1407
1408 break;
1409 }
1410 }
1411
1412
1413 PRIVATE vrna_sc_t *
init_sc_default(unsigned int n)1414 init_sc_default(unsigned int n)
1415 {
1416 vrna_sc_t *sc, init = {
1417 .type = VRNA_SC_DEFAULT
1418 };
1419
1420 sc = vrna_alloc(sizeof(vrna_sc_t));
1421
1422 if (sc) {
1423 memcpy(sc, &init, sizeof(vrna_sc_t));
1424 nullify(sc);
1425 sc->n = n;
1426 }
1427
1428 return sc;
1429 }
1430
1431
1432 PRIVATE vrna_sc_t *
init_sc_window(unsigned int n)1433 init_sc_window(unsigned int n)
1434 {
1435 vrna_sc_t *sc, init = {
1436 .type = VRNA_SC_WINDOW
1437 };
1438
1439 sc = vrna_alloc(sizeof(vrna_sc_t));
1440
1441 if (sc) {
1442 memcpy(sc, &init, sizeof(vrna_sc_t));
1443 nullify(sc);
1444 sc->n = n;
1445 }
1446
1447 return sc;
1448 }
1449
1450
1451 PRIVATE void
nullify(vrna_sc_t * sc)1452 nullify(vrna_sc_t *sc)
1453 {
1454 if (sc) {
1455 sc->state = STATE_CLEAN;
1456 sc->up_storage = NULL;
1457 sc->bp_storage = NULL;
1458 sc->energy_up = NULL;
1459 sc->energy_stack = NULL;
1460 sc->exp_energy_stack = NULL;
1461 sc->exp_energy_up = NULL;
1462
1463 sc->f = NULL;
1464 sc->exp_f = NULL;
1465 sc->data = NULL;
1466 sc->free_data = NULL;
1467
1468 switch (sc->type) {
1469 case VRNA_SC_DEFAULT:
1470 sc->energy_bp = NULL;
1471 sc->exp_energy_bp = NULL;
1472 break;
1473
1474 case VRNA_SC_WINDOW:
1475 sc->energy_bp_local = NULL;
1476 sc->exp_energy_bp_local = NULL;
1477 break;
1478 }
1479 }
1480 }
1481