1 /* constraints handling */
2
3 #ifdef HAVE_CONFIG_H
4 #include "config.h"
5 #endif
6
7 #include <assert.h>
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11 #include <ctype.h>
12 #include <string.h>
13 #include <limits.h>
14
15 #include "ViennaRNA/params/default.h"
16 #include "ViennaRNA/params/constants.h" /* defines MINPSCORE */
17 #include "ViennaRNA/fold_vars.h"
18 #include "ViennaRNA/utils/basic.h"
19 #include "ViennaRNA/utils/alignments.h"
20 #include "ViennaRNA/io/file_formats.h"
21 #include "ViennaRNA/params/basic.h"
22 #include "ViennaRNA/constraints/basic.h"
23 #include "ViennaRNA/constraints/hard.h"
24
25
26 #ifdef __GNUC__
27 # define INLINE inline
28 #else
29 # define INLINE
30 #endif
31
32 #define STATE_CLEAN (unsigned char)0
33 #define STATE_DIRTY_UP (unsigned char)1
34 #define STATE_DIRTY_BP (unsigned char)2
35 #define STATE_UNINITIALIZED (unsigned char)4
36
37 #include "hc_depot.inc"
38
39 /*
40 #################################
41 # GLOBAL VARIABLES #
42 #################################
43 */
44
45 /*
46 #################################
47 # PRIVATE VARIABLES #
48 #################################
49 */
50
51 /*
52 #################################
53 # PRIVATE FUNCTION DECLARATIONS #
54 #################################
55 */
56 PRIVATE unsigned char
57 default_pair_constraint(vrna_fold_compound_t *fc,
58 int i,
59 int j);
60
61
62 PRIVATE INLINE void
63 populate_hc_up(vrna_fold_compound_t *fc,
64 unsigned int i,
65 unsigned int options);
66
67
68 PRIVATE INLINE void
69 populate_hc_bp(vrna_fold_compound_t *fc,
70 unsigned int i,
71 unsigned int options);
72
73
74 PRIVATE void
75 apply_DB_constraint(vrna_fold_compound_t *vc,
76 const char *constraint,
77 unsigned int options);
78
79
80 PRIVATE void
81 hc_reset_to_default(vrna_fold_compound_t *vc);
82
83
84 PRIVATE void
85 hc_update_up(vrna_fold_compound_t *vc);
86
87
88 PRIVATE void
89 hc_update_up_window(vrna_fold_compound_t *vc,
90 int i,
91 unsigned int options);
92
93
94 PRIVATE void
95 default_hc_up(vrna_fold_compound_t *fc,
96 unsigned int options);
97
98 PRIVATE void
99 default_hc_bp(vrna_fold_compound_t *fc,
100 unsigned int options);
101
102 PRIVATE void
103 prepare_hc_up(vrna_fold_compound_t *fc,
104 unsigned int options);
105
106
107 PRIVATE void
108 prepare_hc_bp(vrna_fold_compound_t *fc,
109 unsigned int options);
110
111 /*
112 #################################
113 # BEGIN OF FUNCTION DEFINITIONS #
114 #################################
115 */
116 PUBLIC void
vrna_message_constraint_options_all(void)117 vrna_message_constraint_options_all(void)
118 {
119 vrna_message_constraint_options(VRNA_CONSTRAINT_DB_PIPE
120 | VRNA_CONSTRAINT_DB_DOT
121 | VRNA_CONSTRAINT_DB_X
122 | VRNA_CONSTRAINT_DB_ANG_BRACK
123 | VRNA_CONSTRAINT_DB_RND_BRACK);
124 }
125
126
127 PUBLIC void
vrna_message_constraint_options(unsigned int option)128 vrna_message_constraint_options(unsigned int option)
129 {
130 printf("Input structure constraints using the following notation:\n");
131 if (option & VRNA_CONSTRAINT_DB_PIPE)
132 printf("| : paired with another base\n");
133
134 if (option & VRNA_CONSTRAINT_DB_DOT)
135 printf(". : no constraint at all\n");
136
137 if (option & VRNA_CONSTRAINT_DB_X)
138 printf("x : base must not pair\n");
139
140 if (option & VRNA_CONSTRAINT_DB_ANG_BRACK)
141 printf("< : base i is paired downstream with a base i < j\n"
142 "> : base i is paired upstream with a base j < i\n");
143
144 if (option & VRNA_CONSTRAINT_DB_RND_BRACK)
145 printf("matching brackets ( ): base i pairs base j\n");
146 }
147
148
149 PUBLIC void
vrna_hc_init(vrna_fold_compound_t * vc)150 vrna_hc_init(vrna_fold_compound_t *vc)
151 {
152 unsigned int n;
153 vrna_hc_t *hc;
154
155 n = vc->length;
156
157 /* free previous hard constraints */
158 vrna_hc_free(vc->hc);
159
160 /* allocate memory new hard constraints data structure */
161 hc = (vrna_hc_t *)vrna_alloc(sizeof(vrna_hc_t));
162 hc->type = VRNA_HC_DEFAULT;
163 hc->n = n;
164 hc->mx = (unsigned char *)vrna_alloc(sizeof(unsigned char) * ((n + 1) * (n + 1)));
165 hc->up_ext = (int *)vrna_alloc(sizeof(int) * (n + 2));
166 hc->up_hp = (int *)vrna_alloc(sizeof(int) * (n + 2));
167 hc->up_int = (int *)vrna_alloc(sizeof(int) * (n + 2));
168 hc->up_ml = (int *)vrna_alloc(sizeof(int) * (n + 2));
169 hc->depot = NULL;
170 hc->state = STATE_UNINITIALIZED;
171
172 /* set new hard constraints */
173 vc->hc = hc;
174
175 /* prefill default values */
176 hc_reset_to_default(vc);
177
178 /* add null pointers for the generalized hard constraint feature */
179 hc->f = NULL;
180 hc->data = NULL;
181 hc->free_data = NULL;
182
183 /* update */
184 hc_update_up(vc);
185 }
186
187
188 PUBLIC void
vrna_hc_init_window(vrna_fold_compound_t * vc)189 vrna_hc_init_window(vrna_fold_compound_t *vc)
190 {
191 unsigned int n;
192 vrna_hc_t *hc;
193
194 n = vc->length;
195
196 /* free previous hard constraints */
197 vrna_hc_free(vc->hc);
198
199 /* allocate memory new hard constraints data structure */
200 hc = (vrna_hc_t *)vrna_alloc(sizeof(vrna_hc_t));
201 hc->type = VRNA_HC_WINDOW;
202 hc->n = n;
203 hc->matrix_local = (unsigned char **)vrna_alloc(sizeof(unsigned char *) * (n + 2));
204 hc->up_ext = NULL;
205 hc->up_hp = NULL;
206 hc->up_int = NULL;
207 hc->up_ml = NULL;
208 hc->depot = NULL;
209 hc->state = STATE_UNINITIALIZED;
210
211 /* set new hard constraints */
212 vc->hc = hc;
213
214 /* add null pointers for the generalized hard constraint feature */
215 hc->f = NULL;
216 hc->data = NULL;
217 hc->free_data = NULL;
218 }
219
220
221 PUBLIC void
vrna_hc_update(vrna_fold_compound_t * fc,unsigned int i,unsigned int options)222 vrna_hc_update(vrna_fold_compound_t *fc,
223 unsigned int i,
224 unsigned int options)
225 {
226 unsigned int n, maxdist;
227 vrna_hc_t *hc;
228
229 if (fc) {
230 n = fc->length;
231 maxdist = fc->window_size;
232 hc = fc->hc;
233
234 if (i > n) {
235 vrna_message_warning("vrna_hc_update(): Position %u out of range!",
236 " (Sequence length: %u)",
237 i, n);
238 } else {
239 /* init up_xx arrays if necessary */
240 if (!hc->up_ext) {
241 hc->up_ext = (int *)vrna_alloc(sizeof(int) * (n + 2));
242 hc->up_hp = (int *)vrna_alloc(sizeof(int) * (n + 2));
243 hc->up_int = (int *)vrna_alloc(sizeof(int) * (n + 2));
244 hc->up_ml = (int *)vrna_alloc(sizeof(int) * (n + 2));
245
246 hc_update_up(fc);
247 }
248
249 populate_hc_up(fc, i, options);
250 populate_hc_bp(fc, i, options);
251 }
252 }
253 }
254
255
256 PUBLIC int
vrna_hc_prepare(vrna_fold_compound_t * fc,unsigned int options)257 vrna_hc_prepare(vrna_fold_compound_t *fc,
258 unsigned int options)
259 {
260 int ret = 0;
261
262
263 if (fc) {
264 if (options & VRNA_OPTION_WINDOW) {
265 /* check for minimal hard constraints structure */
266 if ((!fc->hc) || (fc->hc->type != VRNA_HC_WINDOW) || (!fc->hc->matrix_local))
267 vrna_hc_init_window(fc);
268 } else {
269 if (fc->hc->state & STATE_UNINITIALIZED) {
270 default_hc_up(fc, options);
271 default_hc_bp(fc, options);
272 }
273
274 if (fc->hc->state & STATE_DIRTY_UP)
275 prepare_hc_up(fc, options);
276
277 if (fc->hc->state & STATE_DIRTY_BP)
278 prepare_hc_bp(fc, options);
279
280 if (fc->hc->state & ~STATE_CLEAN)
281 hc_update_up(fc);
282 }
283
284 fc->hc->state = STATE_CLEAN;
285 ret = 1;
286 }
287
288 return ret;
289 }
290
291
292 PUBLIC void
vrna_hc_add_up(vrna_fold_compound_t * fc,int i,unsigned char option)293 vrna_hc_add_up(vrna_fold_compound_t *fc,
294 int i,
295 unsigned char option)
296 {
297 unsigned int actual_i, strand_i;
298
299 if (fc) {
300 if (fc->hc) {
301 if ((i <= 0) || (i > fc->length)) {
302 vrna_message_warning("vrna_hc_add_up: position out of range, not doing anything");
303 return;
304 }
305
306 strand_i = fc->strand_number[i];
307 actual_i = (unsigned int)i - fc->strand_start[strand_i] + 1;
308
309 hc_depot_store_up(fc,
310 actual_i,
311 strand_i,
312 option);
313
314 fc->hc->state |= STATE_DIRTY_UP;
315 }
316 }
317 }
318
319
320 PUBLIC int
vrna_hc_add_up_strand(vrna_fold_compound_t * fc,unsigned int i,unsigned int strand,unsigned char option)321 vrna_hc_add_up_strand(vrna_fold_compound_t *fc,
322 unsigned int i,
323 unsigned int strand,
324 unsigned char option)
325 {
326 int ret = 0;
327
328 if ((fc) &&
329 (fc->hc) &&
330 (strand < fc->strands) &&
331 (i > 0)) {
332 unsigned int n_i = (fc->type == VRNA_FC_TYPE_SINGLE) ?
333 fc->nucleotides[strand].length :
334 fc->alignment[strand].sequences[0].length;
335
336 if (i > n_i)
337 return ret;
338
339 hc_depot_store_up(fc,
340 i,
341 strand,
342 option);
343
344 fc->hc->state |= STATE_DIRTY_UP;
345
346 ret = 1;
347 }
348
349 return ret;
350 }
351
352
353 PUBLIC int
vrna_hc_add_up_batch(vrna_fold_compound_t * fc,vrna_hc_up_t * constraints)354 vrna_hc_add_up_batch(vrna_fold_compound_t *fc,
355 vrna_hc_up_t *constraints)
356 {
357 unsigned char options;
358 unsigned int strand_i, actual_i, *ss, *sn;
359 int i, ret, pos;
360
361 ret = 0; /* failure */
362
363 if ((fc) &&
364 (constraints)) {
365 if (fc->hc) {
366 sn = fc->strand_number;
367 ss = fc->strand_start;
368
369 for (i = 0; constraints[i].position != 0; i++) {
370 pos = constraints[i].position;
371 options = constraints[i].options;
372
373 if ((pos <= 0) || (pos > fc->length))
374 break;
375
376 strand_i = sn[pos];
377 actual_i = pos - ss[strand_i] + 1;
378
379 hc_depot_store_up(fc,
380 actual_i,
381 strand_i,
382 options);
383
384 ret++;
385 }
386 }
387 }
388
389 if (ret)
390 fc->hc->state |= STATE_DIRTY_UP;
391
392 return ret;
393 }
394
395
396 PUBLIC int
vrna_hc_add_up_strand_batch(vrna_fold_compound_t * fc,vrna_hc_up_t * constraints)397 vrna_hc_add_up_strand_batch(vrna_fold_compound_t *fc,
398 vrna_hc_up_t *constraints)
399 {
400 unsigned char options;
401 unsigned int i, strand, pos, n_pos, *ss, *sn;
402 int ret;
403
404 ret = 0; /* failure */
405
406 if ((fc) &&
407 (constraints)) {
408 if (fc->hc) {
409 sn = fc->strand_number;
410 ss = fc->strand_start;
411
412 for (i = 0; constraints[i].position != 0; i++) {
413 pos = (unsigned int)constraints[i].position;
414 strand = (unsigned int)constraints[i].strand;
415 options = constraints[i].options;
416
417 if (strand < fc->strands) {
418 n_pos = (fc->type == VRNA_FC_TYPE_SINGLE) ?
419 fc->nucleotides[strand].length :
420 fc->alignment[strand].sequences[0].length;
421
422 if (pos > n_pos)
423 break;
424
425 hc_depot_store_up(fc,
426 pos,
427 strand,
428 options);
429
430 ret++;
431 } else {
432 break;
433 }
434 }
435 }
436 }
437
438 if (ret)
439 fc->hc->state |= STATE_DIRTY_UP;
440
441 return ret;
442 }
443
444
445 PUBLIC void
vrna_hc_add_bp_nonspecific(vrna_fold_compound_t * vc,int i,int d,unsigned char option)446 vrna_hc_add_bp_nonspecific(vrna_fold_compound_t *vc,
447 int i,
448 int d,
449 unsigned char option)
450 {
451 unsigned int n, strand, actual_i, *sn, *ss, *se;
452 int p;
453 vrna_hc_t *hc;
454
455 if (vc) {
456
457 if (vc->hc) {
458 if ((i <= 0) || (i > vc->length)) {
459 vrna_message_warning("vrna_hc_add_bp_nonspecific: position out of range, not doing anything");
460 return;
461 }
462
463 hc = vc->hc;
464 n = hc->n;
465 sn = vc->strand_number;
466 ss = vc->strand_start;
467 se = vc->strand_end;
468 strand = sn[i];
469 actual_i = i - ss[strand] + 1;
470
471 hc_depot_store_nonspec(vc,
472 actual_i,
473 strand,
474 d,
475 option);
476
477 hc->state |= STATE_DIRTY_UP;
478 }
479 }
480 }
481
482
483 PUBLIC int
vrna_hc_add_bp_strand(vrna_fold_compound_t * fc,unsigned int i,unsigned int strand_i,unsigned int j,unsigned int strand_j,unsigned char option)484 vrna_hc_add_bp_strand(vrna_fold_compound_t *fc,
485 unsigned int i,
486 unsigned int strand_i,
487 unsigned int j,
488 unsigned int strand_j,
489 unsigned char option)
490 {
491 unsigned int n, *sn, *se, *ss, n_i, n_j, k, l;
492 int ret, turn;
493 vrna_hc_t *hc;
494
495 ret = 0;
496
497 if ((fc) &&
498 (fc->hc) &&
499 (strand_i < fc->strands) &&
500 (strand_j < fc->strands) &&
501 (i > 0) &&
502 (j > 0)) {
503 sn = fc->strand_number;
504 ss = fc->strand_start;
505 se = fc->strand_end;
506 hc = fc->hc;
507 n_i = (fc->type == VRNA_FC_TYPE_SINGLE) ?
508 fc->nucleotides[strand_i].length :
509 fc->alignment[strand_i].sequences[0].length;
510 n_j = (fc->type == VRNA_FC_TYPE_SINGLE) ?
511 fc->nucleotides[strand_j].length :
512 fc->alignment[strand_j].sequences[0].length;
513 turn = fc->params->model_details.min_loop_size;
514
515 /* check for strand length ranges */
516 if ((i > n_i) || (j > n_j))
517 return ret;
518
519 /* check for minimum hairpin loop length for intra-molecular pairs */
520 if ((strand_i == strand_j) &&
521 ((j - i - 1) < turn))
522 return ret;
523
524 /* store the constraint */
525 hc_depot_store_bp(fc,
526 i,
527 strand_i,
528 j,
529 strand_j,
530 option);
531
532 hc->state |= STATE_DIRTY_BP;
533
534 ret = 1;
535 }
536
537 return ret;
538 }
539
540
541 PUBLIC int
vrna_hc_add_bp(vrna_fold_compound_t * vc,int i,int j,unsigned char option)542 vrna_hc_add_bp(vrna_fold_compound_t *vc,
543 int i,
544 int j,
545 unsigned char option)
546 {
547 unsigned int n, *sn, *se, *ss, strand_i, strand_j, actual_i, actual_j;
548 int k, l, ret;
549 vrna_hc_t *hc;
550
551 ret = 0;
552
553 if (vc) {
554 sn = vc->strand_number;
555 ss = vc->strand_start;
556 se = vc->strand_end;
557
558 if (vc->hc) {
559 if ((i <= 0) || (j <= i) || (j > vc->length)) {
560 vrna_message_warning("vrna_hc_add_bp: position out of range, omitting constraint");
561 } else if ((sn[i] == sn[j]) && ((j - i - 1) < vc->params->model_details.min_loop_size)) {
562 vrna_message_warning(
563 "vrna_hc_add_bp: Pairing partners (%d, %d) violate minimum loop size settings of %dnt, omitting constraint",
564 i,
565 j,
566 vc->params->model_details.min_loop_size);
567 } else {
568
569 hc = vc->hc;
570 n = hc->n;
571
572 /*
573 determine the corresponding strand numbers and the actual
574 position (relative to the strand)
575 */
576 strand_i = sn[i];
577 strand_j = sn[j];
578 actual_i = i - ss[strand_i] + 1;
579 actual_j = j - ss[strand_j] + 1;
580
581 ret = vrna_hc_add_bp_strand(vc,
582 actual_i,
583 strand_i,
584 actual_j,
585 strand_j,
586 option);
587 }
588 }
589 }
590
591 return ret;
592 }
593
594
595 PUBLIC void
vrna_hc_free(vrna_hc_t * hc)596 vrna_hc_free(vrna_hc_t *hc)
597 {
598 if (hc) {
599 if (hc->type == VRNA_HC_DEFAULT) {
600 free(hc->mx);
601 } else if (hc->type == VRNA_HC_WINDOW) {
602 unsigned int i;
603 free(hc->matrix_local);
604 }
605
606 hc_depot_free(hc);
607
608 free(hc->up_ext);
609 free(hc->up_hp);
610 free(hc->up_int);
611 free(hc->up_ml);
612
613 if (hc->free_data)
614 hc->free_data(hc->data);
615
616 free(hc);
617 }
618 }
619
620
621 PUBLIC void
vrna_hc_add_f(vrna_fold_compound_t * vc,vrna_callback_hc_evaluate * f)622 vrna_hc_add_f(vrna_fold_compound_t *vc,
623 vrna_callback_hc_evaluate *f)
624 {
625 if (vc && f) {
626 if (vc->type == VRNA_FC_TYPE_SINGLE) {
627 if (!vc->hc)
628 vrna_hc_init(vc);
629
630 vc->hc->f = f;
631 }
632 }
633 }
634
635
636 PUBLIC void
vrna_hc_add_data(vrna_fold_compound_t * vc,void * data,vrna_callback_free_auxdata * f)637 vrna_hc_add_data(vrna_fold_compound_t *vc,
638 void *data,
639 vrna_callback_free_auxdata *f)
640 {
641 if (vc && data) {
642 if (vc->type == VRNA_FC_TYPE_SINGLE) {
643 if (!vc->hc)
644 vrna_hc_init(vc);
645
646 vc->hc->data = data;
647 vc->hc->free_data = f;
648 }
649 }
650 }
651
652
653 PUBLIC int
vrna_hc_add_from_db(vrna_fold_compound_t * vc,const char * constraint,unsigned int options)654 vrna_hc_add_from_db(vrna_fold_compound_t *vc,
655 const char *constraint,
656 unsigned int options)
657 {
658 const char *structure_constraint;
659 char *tmp;
660 int ret;
661
662 ret = 0; /* Failure */
663
664 if (vc) {
665 tmp = NULL;
666 if ((!vc->params) && (!vc->exp_params))
667 return ret;
668
669 if (!vc->hc)
670 vrna_hc_init(vc);
671
672 if (options & VRNA_CONSTRAINT_DB_WUSS) {
673 tmp = vrna_db_from_WUSS(constraint);
674 structure_constraint = (const char *)tmp;
675 } else {
676 structure_constraint = constraint;
677 }
678
679 /* apply hard constraints from dot-bracket notation */
680 apply_DB_constraint(vc, structure_constraint, options);
681
682 ret = 1; /* Success */
683
684 free(tmp);
685 }
686
687 return ret;
688 }
689
690
691 /*
692 #####################################
693 # BEGIN OF STATIC HELPER FUNCTIONS #
694 #####################################
695 */
696 PRIVATE unsigned char
default_pair_constraint(vrna_fold_compound_t * fc,int i,int j)697 default_pair_constraint(vrna_fold_compound_t *fc,
698 int i,
699 int j)
700 {
701 unsigned char constraint, can_stack;
702 short *S;
703 unsigned int *sn;
704 int type;
705 vrna_md_t *md;
706
707 sn = fc->strand_number;
708 md = &(fc->params->model_details);
709 constraint = VRNA_CONSTRAINT_CONTEXT_NONE;
710
711 switch (fc->type) {
712 case VRNA_FC_TYPE_SINGLE:
713 S = fc->sequence_encoding2;
714 if (((j - i) < md->max_bp_span) &&
715 ((sn[i] != sn[j]) || ((j - i) > md->min_loop_size))) {
716 type = md->pair[S[i]][S[j]];
717 switch (type) {
718 case 0:
719 break;
720 case 3:
721 /* fallthrough */
722 case 4:
723 if (md->noGU) {
724 break;
725 } else if (md->noGUclosure) {
726 constraint = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
727 constraint &= ~(VRNA_CONSTRAINT_CONTEXT_HP_LOOP | VRNA_CONSTRAINT_CONTEXT_MB_LOOP);
728 break;
729 }
730
731 /* else fallthrough */
732 default:
733 constraint = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
734 break;
735 }
736
737 if (md->noLP) {
738 /* check, whether this nucleotide can actually stack with anything or only forms isolated pairs */
739 can_stack = VRNA_CONSTRAINT_CONTEXT_NONE;
740
741 /* can it be enclosed by another base pair? */
742 if ((i > 1) &&
743 (j < fc->length) &&
744 (((j - i + 2) < md->max_bp_span) || (sn[i - 1] != sn[j + 1])) &&
745 (md->pair[S[i - 1]][S[j + 1]]))
746 can_stack = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
747
748 /* can it enclose another base pair? */
749 if ((i + 2 < j) &&
750 (((j - i - 2) > md->min_loop_size) || (sn[i + 1] != sn[j - 1])) &&
751 (md->pair[S[i + 1]][S[j - 1]]))
752 can_stack = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
753
754 constraint &= can_stack;
755 }
756 }
757
758 break;
759
760 case VRNA_FC_TYPE_COMPARATIVE:
761 if ((sn[i] != sn[j]) ||
762 (((j - i + 1) <= md->max_bp_span) && ((j - i - 1) >= md->min_loop_size))) {
763 int min_score = md->cv_fact * MINPSCORE;
764 int act_score = (fc->hc->type == VRNA_HC_WINDOW) ?
765 fc->pscore_local[i][j - i] :
766 fc->pscore[fc->jindx[j] + i];
767 if (act_score >= min_score)
768 constraint = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
769
770 if (md->noLP) {
771 /* check, whether this nucleotide can actually stack with anything or only forms isolated pairs */
772 can_stack = VRNA_CONSTRAINT_CONTEXT_NONE;
773
774 /* can it be enclosed by another base pair? */
775 if ((i > 1) &&
776 (j < fc->length) &&
777 (((j - i + 2) < md->max_bp_span) || (sn[i - 1] != sn[j + 1]))) {
778 int outer_pscore = (fc->hc->type == VRNA_HC_WINDOW) ?
779 fc->pscore_local[i - 1][j - i + 2] :
780 fc->pscore[fc->jindx[j + 1] + i - 1];
781 if (outer_pscore >= min_score)
782 can_stack = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
783 }
784
785 /* can it enclose another base pair? */
786 if ((i + 2 < j) &&
787 (((j - i - 2) > md->min_loop_size) || (sn[i + 1] != sn[j - 1]))) {
788 int inner_pscore = (fc->hc->type == VRNA_HC_WINDOW) ?
789 fc->pscore_local[i + 1][j - i - 2] :
790 fc->pscore[fc->jindx[j - 1] + i + 1];
791 if (inner_pscore >= min_score)
792 can_stack = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
793 }
794
795 constraint &= can_stack;
796 }
797 }
798
799 break;
800 }
801
802 return constraint;
803 }
804
805
806 PRIVATE INLINE void
populate_hc_up(vrna_fold_compound_t * fc,unsigned int i,unsigned int options)807 populate_hc_up(vrna_fold_compound_t *fc,
808 unsigned int i,
809 unsigned int options)
810 {
811 unsigned char context;
812 unsigned int actual_i, strand;
813
814 vrna_hc_t *hc = fc->hc;
815
816 if (hc->type == VRNA_HC_WINDOW) {
817 strand = fc->strand_number[i];
818 actual_i = fc->strand_start[strand] + i - 1;
819
820 if ((hc->depot) &&
821 (hc->depot->up) &&
822 (hc->depot->up_size[strand] >= i)) {
823 /* We use user-defined constraints for unpaired nucleotides */
824 context = hc->depot->up[strand][i].context;
825
826 if (hc->depot->up[strand][i].nonspec) {
827 /* this nucleotide must not stay unpaired */
828 hc->matrix_local[actual_i][0] = VRNA_CONSTRAINT_CONTEXT_NONE;
829
830 /* pairing direction will be enforced by populate_hc_bp() */
831 } else if (context & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
832 hc->matrix_local[i][0] = context & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
833 /* restriction of base pairing contexts and removal of possible pairing will be handled by populate_hc_bp() */
834 } else {
835 hc->matrix_local[i][0] = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
836 /* restriction of base pairing contexts and removal of possible pairing will be handled by populate_hc_bp() */
837 }
838 } else {
839 /* if no constraints are available, we simply allow it to be unpaired in all contexts */
840 hc->matrix_local[i][0] = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
841 }
842 } else {
843 /* do something reasonable here... */
844 }
845 }
846
847
848 PRIVATE void
default_hc_up(vrna_fold_compound_t * fc,unsigned int options)849 default_hc_up(vrna_fold_compound_t *fc,
850 unsigned int options)
851 {
852 unsigned int i, n;
853 vrna_hc_t *hc;
854 vrna_hc_depot_t *depot;
855
856 hc = fc->hc;
857
858 if (options & VRNA_OPTION_WINDOW) {
859
860 } else {
861 n = fc->length;
862
863 for (i = 1; i <= n; i++) {
864 hc->mx[n * i + i] = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
865 }
866 }
867 }
868
869 PRIVATE void
prepare_hc_up(vrna_fold_compound_t * fc,unsigned int options)870 prepare_hc_up(vrna_fold_compound_t *fc,
871 unsigned int options)
872 {
873 unsigned char option, type, t1, t2;
874 unsigned int i, j, k, n, s, *ss;
875 vrna_hc_t *hc;
876 vrna_hc_depot_t *depot;
877
878 hc = fc->hc;
879
880 if (options & VRNA_OPTION_WINDOW) {
881
882 } else {
883 n = fc->length;
884 ss = fc->strand_start;
885 depot = hc->depot;
886
887 /* 2. apply constraints as stored in depot */
888 if ((depot) && (depot->up)) {
889 for (s = 0; s < depot->strands; s++) {
890 for (k = 1; k <= depot->up_size[s]; k++) {
891 /* process nucleotide-specific constraint */
892 option = depot->up[s][k].context;
893 i = ss[s] + k - 1; /* constraint position in current strand order */
894
895 if (depot->up[s][k].nonspec) {
896 /* this is actually a must-pair constraint */
897 type = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
898 /* acknowledge pairing direction */
899 t1 = (depot->up[s][k].direction <= 0) ? type : VRNA_CONSTRAINT_CONTEXT_NONE;
900 t2 = (depot->up[s][k].direction >= 0) ? type : VRNA_CONSTRAINT_CONTEXT_NONE;
901
902 if (option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE) {
903 /* only allow for possibly non-canonical pairs, do not enforce them */
904 for (j = 1; j < i; j++) {
905 hc->mx[n * i + j] |= t1;
906 hc->mx[n * j + i] |= t1;
907 }
908 for (j = i + 1; j <= n; j++) {
909 hc->mx[n * i + j] |= t2;
910 hc->mx[n * j + i] |= t2;
911 }
912 } else {
913 /* force pairing direction */
914 for (j = 1; j < i; j++) {
915 hc->mx[n * i + j] &= t1;
916 hc->mx[n * j + i] &= t1;
917 }
918 for (j = i + 1; j <= n; j++) {
919 hc->mx[n * i + j] &= t2;
920 hc->mx[n * j + i] &= t2;
921 }
922 /* nucleotide mustn't be unpaired */
923 hc->mx[n * i + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
924 }
925 } else {
926 /* 'regular' nucleotide-specific constraint */
927 if (option & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
928 /* force nucleotide to appear unpaired within a certain type of loop */
929 /* do not allow i to be paired with any other nucleotide */
930 if (!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
931 for (j = 1; j < i; j++) {
932 hc->mx[n * i + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
933 hc->mx[n * j + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
934 }
935 for (j = i + 1; j <= n; j++) {
936 hc->mx[n * i + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
937 hc->mx[n * j + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
938 }
939 }
940
941 type = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
942
943 hc->mx[n * i + i] = type;
944 } else {
945 type = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
946
947 /* do not allow i to be paired with any other nucleotide (in context type) */
948 if (!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
949 for (j = 1; j < i; j++) {
950 hc->mx[n * i + j] &= ~type;
951 hc->mx[n * j + i] &= ~type;
952 }
953 for (j = i + 1; j <= n; j++) {
954 hc->mx[n * i + j] &= ~type;
955 hc->mx[n * j + i] &= ~type;
956 }
957 }
958
959 hc->mx[n * i + i] = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
960 }
961 }
962 }
963 }
964 }
965 }
966 }
967
968
969 PRIVATE void
default_hc_bp(vrna_fold_compound_t * fc,unsigned int options)970 default_hc_bp(vrna_fold_compound_t *fc,
971 unsigned int options)
972 {
973 unsigned int i, j, n;
974 vrna_hc_t *hc;
975
976 hc = fc->hc;
977
978 if (options & VRNA_OPTION_WINDOW) {
979
980 } else {
981 n = fc->length;
982
983 for (j = n; j > 1; j--) {
984 for (i = 1; i < j; i++) {
985 hc->mx[n * i + j] = default_pair_constraint(fc, i, j);
986 hc->mx[n * j + i] = hc->mx[n * i + j];
987 }
988 }
989 }
990 }
991
992
993 PRIVATE void
prepare_hc_bp(vrna_fold_compound_t * fc,unsigned int options)994 prepare_hc_bp(vrna_fold_compound_t *fc,
995 unsigned int options)
996 {
997 unsigned char option, type;
998 unsigned int i, j, k, p, q, n, actual_i, actual_j, strand_j, start, end, t_start, t_end, s, *sn, *ss, *se;
999 int *idx;
1000 vrna_hc_t *hc;
1001 vrna_hc_depot_t *depot;
1002
1003 hc = fc->hc;
1004 depot = hc->depot;
1005 sn = fc->strand_number;
1006 ss = fc->strand_start;
1007 se = fc->strand_end;
1008
1009 if ((!depot) || (!depot->bp))
1010 return;
1011
1012 if (options & VRNA_OPTION_WINDOW) {
1013
1014 } else {
1015 n = fc->length;
1016 idx = fc->jindx;
1017
1018 /* 2. apply constraints stored in depot */
1019 for (s = 0; s < depot->strands; s++) {
1020 for (actual_i = 1; actual_i <= depot->bp_size[s]; actual_i++) {
1021 for (k = 0; k < depot->bp[s][actual_i].list_size; k++) {
1022 option = depot->bp[s][actual_i].context[k];
1023 strand_j = depot->bp[s][actual_i].strand_j[k];
1024 actual_j = depot->bp[s][actual_i].j[k];
1025 i = ss[s] + actual_i - 1; /* constraint position in current strand order */
1026 j = ss[strand_j] + actual_j - 1; /* constraint position in current strand order */
1027
1028 if (i < j) {
1029 /* apply the constraint */
1030 hc->mx[n * i + j] = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1031 hc->mx[n * j + i] = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1032
1033 /* is the ptype reset actually required??? */
1034 if ((fc->type == VRNA_FC_TYPE_SINGLE) &&
1035 (option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS)) {
1036 /* reset ptype in case (i,j) is a non-canonical pair */
1037 if (fc->ptype[idx[j] + i] == 0)
1038 fc->ptype[idx[j] + i] = 7;
1039 }
1040
1041 if (!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1042 /*
1043 * remove all conflicting base pairs, i.e. do not allow i,j to pair
1044 * with any other nucleotide k
1045 */
1046 for (p = 1; p < i; p++) {
1047 hc->mx[n * i + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1048 hc->mx[n * p + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1049 hc->mx[n * j + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1050 hc->mx[n * p + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
1051
1052 for (q = i + 1; q < j; q++) {
1053 hc->mx[n * p + q] = VRNA_CONSTRAINT_CONTEXT_NONE;
1054 hc->mx[n * q + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1055 }
1056 }
1057 for (p = i + 1; p < j; p++) {
1058 hc->mx[n * i + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1059 hc->mx[n * p + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1060 hc->mx[n * j + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1061 hc->mx[n * p + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
1062
1063 for (q = j + 1; q <= n; q++) {
1064 hc->mx[n * p + q] = VRNA_CONSTRAINT_CONTEXT_NONE;
1065 hc->mx[n * q + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1066 }
1067 }
1068 for (p = j + 1; p <= n; p++) {
1069 hc->mx[n * i + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1070 hc->mx[n * p + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1071 hc->mx[n * j + p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1072 hc->mx[n * p + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
1073 }
1074 }
1075
1076 if (option & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
1077 /* do not allow i,j to be unpaired */
1078 hc->mx[n * i + i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1079 hc->mx[n * j + j] = VRNA_CONSTRAINT_CONTEXT_NONE;
1080 }
1081 }
1082 }
1083 }
1084 }
1085 }
1086 }
1087
1088
1089 PRIVATE INLINE void
populate_hc_bp(vrna_fold_compound_t * fc,unsigned int i,unsigned int options)1090 populate_hc_bp(vrna_fold_compound_t *fc,
1091 unsigned int i,
1092 unsigned int options)
1093 {
1094 unsigned char constraint, type, t1, t2;
1095 unsigned int max_span, maxdist, j, k, p, n, l, turn, strand, sj, sl, actual_i, actual_j, actual_l, *sn, *ss;
1096 vrna_hc_t *hc;
1097
1098 n = fc->length;
1099 maxdist = fc->window_size;
1100 max_span = fc->params->model_details.max_bp_span;
1101 sn = fc->strand_number;
1102 ss = fc->strand_start;
1103 strand = sn[i];
1104 actual_i = i - ss[strand] + 1;
1105 turn = fc->params->model_details.min_loop_size;
1106 hc = fc->hc;
1107
1108 if (options & VRNA_CONSTRAINT_WINDOW_UPDATE_3) {
1109 /* the sliding window moves from 3' to 5' side */
1110
1111 /* apply default constraints first */
1112 for (k = turn + 1; k < max_span; k++) {
1113 j = i + k;
1114 if (j > n)
1115 break;
1116
1117 constraint = default_pair_constraint(fc, i, j);
1118
1119 hc->matrix_local[i][j - i] = constraint;
1120 }
1121
1122 if (hc->depot) {
1123 /* 1. apply remainder of (partly) applied nucleotide-specific constraints */
1124 if (hc->depot->up) {
1125
1126 /* 1.a apply nucleotide specific constraints for i */
1127 if ((hc->depot->up[strand]) &&
1128 (hc->depot->up_size[strand] >= actual_i)) {
1129
1130 constraint = hc->depot->up[strand][actual_i].context;
1131
1132 if (hc->depot->up[strand][actual_i].nonspec) {
1133 /* handle unspecific pairing contraint */
1134
1135 /* remove downstream pairs if necessary */
1136 if (hc->depot->up[strand][actual_i].direction < 0) {
1137 for (k = i + turn + 1; k < MIN2(i + max_span, n + 1); k++)
1138 hc->matrix_local[i][k - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1139 } else {
1140 /* enforce the base pair type */
1141 for (k = i + turn + 1; k < MIN2(i + max_span, n + 1); k++)
1142 hc->matrix_local[i][k - i] &= constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1143 }
1144 } else {
1145 /* handle 'regular' unpairedness constraint */
1146
1147 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
1148 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1149 /* do not allow i to be paired with any other nucleotide */
1150 for (k = i + turn + 1; k < MIN2(i + max_span, n + 1); k++)
1151 hc->matrix_local[i][k - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1152 }
1153 } else {
1154 type = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1155 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1156 /* only allow i to be paired for particular types */
1157 for (k = i + turn + 1; k < MIN2(i + max_span, n + 1); k++)
1158 hc->matrix_local[i][k - i] &= ~type;
1159 }
1160 }
1161 }
1162 }
1163
1164
1165 /*
1166 1.b apply remainder of (partly) applied nucleotide-specific
1167 constraints for j with i < j < i + maxdist
1168 */
1169 for (k = 1; k < max_span; k++) {
1170 j = i + k;
1171 sj = sn[j];
1172 actual_j = j - ss[sj] + 1;
1173
1174 if (j > n)
1175 break;
1176
1177 if ((hc->depot->up[sj]) &&
1178 (hc->depot->up_size[sj] >= actual_j)) {
1179 constraint = hc->depot->up[sj][actual_j].context;
1180
1181 if (hc->depot->up[sj][actual_j].nonspec) {
1182 /* j must be paired, check preferred direction */
1183 if (hc->depot->up[sj][actual_j].direction > 0) {
1184 /* j must pair downstream, so we remove base pair (i, j) */
1185 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1186 } else {
1187 /* enforce base pair type */
1188 hc->matrix_local[i][j - i] &= constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1189 }
1190 } else if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1191 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
1192 /* remove (i, j) base pair */
1193 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1194 } else {
1195 constraint = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1196 hc->matrix_local[i][j - i] &= ~constraint;
1197 }
1198 }
1199 }
1200 }
1201 }
1202
1203 /* 2. now for the actual base pair constraints */
1204 if (hc->depot->bp) {
1205
1206 /* 2.a apply base pair specific constraint for nucleotide i */
1207 if ((hc->depot->bp[strand]) &&
1208 (hc->depot->bp_size[strand] >= actual_i) &&
1209 (hc->depot->bp[strand][actual_i].list_size > 0)) {
1210 /* go through list of all constraints for this nucleotide */
1211 for (size_t cnt = 0; cnt < hc->depot->bp[strand][actual_i].list_size; cnt++) {
1212 constraint = hc->depot->bp[strand][actual_i].context[cnt];
1213 actual_j = hc->depot->bp[strand][actual_i].j[cnt];
1214 sj = hc->depot->bp[strand][actual_i].strand_j[cnt];
1215 j = ss[sj] + actual_j - 1;
1216
1217 /* apply the constraint */
1218
1219 /* do not allow i to stay unpaired if necessary */
1220 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE)
1221 hc->matrix_local[i][0] = VRNA_CONSTRAINT_CONTEXT_NONE;
1222
1223 if (j > i) {
1224 /* j is 3' base of base pair (i, j) */
1225 if (i + max_span > j)
1226 hc->matrix_local[i][j - i] = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1227
1228 /* remove other base pairs violating the constraint */
1229 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1230 /* remove base pairs (i, k) with k != j */
1231 for (k = i + turn + 1; k < j; k++)
1232 hc->matrix_local[i][k - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1233
1234 for (k = j + 1; k < MIN2(i + max_span, n + 1); k++)
1235 hc->matrix_local[i][k - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1236 }
1237 } else {
1238 /* j is 5' base of base pair (j, i) */
1239
1240 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1241 /* this base pairs upstream, so remove all downstream pairing partners */
1242 for (k = turn + 1; k < max_span; k++) {
1243 j = i + k;
1244 if (j > n)
1245 break;
1246
1247 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1248 }
1249 }
1250 }
1251 }
1252 }
1253
1254 /* 2.b acknowledge base pair constraints that involve j with i < j < i + max_span */
1255 for (k = 1; k < max_span; k++) {
1256 j = i + k;
1257 sj = sn[j];
1258 actual_j = j - ss[sj] + 1;
1259
1260 if (j > n)
1261 break;
1262
1263 if ((hc->depot->bp[sj]) &&
1264 (hc->depot->bp_size[sj] >= actual_j) &&
1265 (hc->depot->bp[sj][actual_j].list_size > 0)) {
1266 /* go through list of all constraints for this nucleotide */
1267 for (size_t cnt = 0; cnt < hc->depot->bp[sj][actual_j].list_size; cnt++) {
1268 constraint = hc->depot->bp[sj][actual_j].context[cnt];
1269 actual_l = hc->depot->bp[sj][actual_j].j[cnt];
1270 sl = hc->depot->bp[sj][actual_j].strand_j[cnt];
1271 l = ss[sl] + actual_l - 1;
1272
1273 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1274 if (l > j) {
1275 /* remove base pairs (i, p) with i < j <= p <= l */
1276 for (p = j; p < MIN2(i + max_span, l + 1); p++)
1277 hc->matrix_local[i][p - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1278 } else if (i < l) {
1279 /* remove base pairs (i, p) with i < l <= p < j */
1280 for (p = l; p <= j; p++)
1281 hc->matrix_local[i][p - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1282 } else if (l < i) {
1283 /* remove base pairs (i, p) with i < j <= p */
1284 for (p = j; p < MIN2(i + max_span, n + 1); p++)
1285 hc->matrix_local[i][p - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1286 }
1287 }
1288 }
1289 }
1290 }
1291 }
1292 }
1293
1294 hc_update_up_window(fc, i, options);
1295 } else if (options & VRNA_CONSTRAINT_WINDOW_UPDATE_5) {
1296 /* the sliding window moves from 5' to 3' side (i is 3' nucleotide) */
1297
1298 j = i;
1299 sj = strand;
1300 actual_j = actual_i;
1301
1302 if ((j <= n) && (j > turn)) {
1303 unsigned int start_i, min_i;
1304
1305 start_i = (j > max_span) ? j - max_span + 1 : 1;
1306 min_i = (j > maxdist) ? j - maxdist + 1 : 1;
1307
1308 /* apply default constraints first */
1309 for (i = start_i; i < j - turn; i++)
1310 hc->matrix_local[i][j - i] = default_pair_constraint(fc, i, j);
1311
1312 /* next are user-defined constraints */
1313 if (hc->depot) {
1314 /* 1. apply remainder of (partly) applied nucleotide-specific constraints */
1315 if (hc->depot->up) {
1316
1317 /* 1.a apply nucleotide-specific constraints for j */
1318 if ((hc->depot->up[sj]) &&
1319 (hc->depot->up_size[sj] >= actual_j)) {
1320
1321 constraint = hc->depot->up[sj][actual_j].context;
1322
1323 if (hc->depot->up[sj][actual_j].nonspec) {
1324 /* handle unspecific pairing constraint */
1325
1326 /* remove upstream pairs if necessary */
1327 if (hc->depot->up[sj][actual_j].direction > 0) {
1328 /* j is only allowed to pair downstream */
1329 /* remove all base pairs (i, j) with j - maxdist < i < j */
1330 for (i = start_i; i < j - turn; i++)
1331 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1332 } else {
1333 /* enforce the base pair type */
1334 for (i = start_i; i < j - turn; i++)
1335 hc->matrix_local[i][j - i] &= constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1336 }
1337 } else {
1338 /* handle 'regular' unpairedness constraint */
1339
1340 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
1341 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1342 /* do not allow j to be paired with any other nucleotide */
1343 for (i = start_i; i < j - turn; i++)
1344 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1345 }
1346 } else {
1347 type = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1348 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1349 /* only allow j to be paired for particular types */
1350 for (i = start_i; i < j - turn; i++)
1351 hc->matrix_local[i][j - i] &= ~type;
1352 }
1353 }
1354 }
1355 }
1356
1357 /* 1.b apply remainder of (partly) applied nucleotide-specific
1358 constraints for i with j - maxdist < i < j
1359 */
1360
1361 for (i = start_i; i < j; i++) {
1362 strand = sn[i];
1363 actual_i = i - ss[strand] + 1;
1364
1365 if ((hc->depot->up[strand]) &&
1366 (hc->depot->up_size[strand] >= actual_i)) {
1367 constraint = hc->depot->up[strand][actual_i].context;
1368
1369 if (hc->depot->up[strand][actual_i].nonspec) {
1370 /* i must be paired, check preferred direction */
1371 if (hc->depot->up[strand][actual_i].direction < 0) {
1372 /* remove (i, j) base pair */
1373 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1374 } else {
1375 /* enforce base pair type */
1376 hc->matrix_local[i][j - i] &= constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1377 }
1378 } else if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1379 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE) {
1380 /* remove (i, j) base pair */
1381 hc->matrix_local[i][j - i] = VRNA_CONSTRAINT_CONTEXT_NONE;
1382 } else {
1383 constraint = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1384 hc->matrix_local[i][j - i] &= ~constraint;
1385 }
1386 }
1387 }
1388 }
1389 }
1390
1391 /* 2. now for the actual base pair constraints */
1392 if (hc->depot->bp) {
1393 /* 2.a apply base pair specific constraint for nucleotide j */
1394 if ((hc->depot->bp[sj]) &&
1395 (hc->depot->bp_size[sj] >= actual_j) &&
1396 (hc->depot->bp[sj][actual_j].list_size > 0)) {
1397 for (size_t cnt = 0; cnt < hc->depot->bp[sj][actual_j].list_size; cnt++) {
1398 constraint = hc->depot->bp[sj][actual_j].context[cnt];
1399 actual_i = hc->depot->bp[sj][actual_j].j[cnt];
1400 strand = hc->depot->bp[sj][actual_j].strand_j[cnt];
1401 i = ss[strand] + actual_i - 1;
1402
1403 /* apply the constraint */
1404 /* do not allow j to stay unpaired if necessary */
1405 if (constraint & VRNA_CONSTRAINT_CONTEXT_ENFORCE)
1406 hc->matrix_local[j][0] = VRNA_CONSTRAINT_CONTEXT_NONE;
1407
1408 if (i < j) {
1409 /* j is 3' base of base pair (i, j) */
1410 if (i + max_span > j)
1411 hc->matrix_local[i][j - i] = constraint & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1412
1413 /* remove other base pairs violating the constraint */
1414 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1415 /* remove base pairs (k, j) with k != i */
1416 for (k = start_i; k < i; k++)
1417 hc->matrix_local[k][j - k] = VRNA_CONSTRAINT_CONTEXT_NONE;
1418
1419 for (k = i + 1; k < j - turn; k++)
1420 hc->matrix_local[k][j - k] = VRNA_CONSTRAINT_CONTEXT_NONE;
1421 }
1422 } else {
1423 /* j is 5' base of base pair (j, i) */
1424
1425 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1426 /* this base pair pairs downstream, so remove all upstream pairing partners */
1427 for (k = start_i; k < j - turn; k++)
1428 hc->matrix_local[k][j - k] = VRNA_CONSTRAINT_CONTEXT_NONE;
1429 }
1430 }
1431 }
1432 }
1433
1434 /* 2.b acknowledge base pair constraints that involve i with j - maxdist < i < j */
1435 for (i = min_i; i < j; i++) {
1436 strand = sn[i];
1437 actual_i = i - ss[strand] + 1;
1438
1439 if ((hc->depot->bp[strand]) &&
1440 (hc->depot->bp_size[strand] >= actual_i) &&
1441 (hc->depot->bp[strand][actual_i].list_size > 0)) {
1442 for (size_t cnt = 0; cnt < hc->depot->bp[strand][actual_i].list_size; cnt++) {
1443 constraint = hc->depot->bp[strand][actual_i].context[cnt];
1444 actual_l = hc->depot->bp[strand][actual_i].j[cnt];
1445 sl = hc->depot->bp[strand][actual_i].strand_j[cnt];
1446 l = ss[sl] + actual_l - 1;
1447
1448 if (!(constraint & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)) {
1449 if (l < i) {
1450 /* remove base pairs (p, j) with l <= p <= i < j */
1451 for (p = MAX2(l, start_i); p <= i; p++)
1452 hc->matrix_local[p][j - p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1453 } else if (l < j) {
1454 /* remove base pairs (p, j) with i <= p <= l < j */
1455 for (p = MAX2(i, start_i); p <= l; p++)
1456 hc->matrix_local[p][j - p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1457 } else if (l > j) {
1458 /* remove base pairs (p, j) with p <= i */
1459 for (p = start_i; p <= i; p++)
1460 hc->matrix_local[p][j - p] = VRNA_CONSTRAINT_CONTEXT_NONE;
1461 }
1462 }
1463 }
1464 }
1465 }
1466 }
1467 }
1468 }
1469
1470 hc_update_up_window(fc, j, options);
1471 }
1472 }
1473
1474
1475 struct hc_bp {
1476 int i;
1477 int j;
1478 unsigned char options;
1479 };
1480
1481
1482 PRIVATE void
apply_DB_constraint(vrna_fold_compound_t * vc,const char * constraint,unsigned int options)1483 apply_DB_constraint(vrna_fold_compound_t *vc,
1484 const char *constraint,
1485 unsigned int options)
1486 {
1487 char *sequence;
1488 short *S;
1489 unsigned int length, min_loop_size, num_up, num_bp, num_bp_unspecific,
1490 size_up, size_bp, size_bp_unspecific;
1491 int n, i, j, hx, *stack, cut;
1492 vrna_md_t *md;
1493 vrna_hc_up_t *up;
1494 struct hc_bp *bp, *bp_unspecific;
1495
1496 if (constraint == NULL)
1497 return;
1498
1499 sequence = vc->sequence;
1500 length = (int)vc->length;
1501 S = vc->sequence_encoding2;
1502 md = &(vc->params->model_details);
1503 min_loop_size = md->min_loop_size;
1504 cut = vc->cutpoint;
1505 n = (int)strlen(constraint);
1506 stack = (int *)vrna_alloc(sizeof(int) * (n + 1));
1507 size_up = size_bp = size_bp_unspecific = 10;
1508 num_up = num_bp = num_bp_unspecific = 0;
1509 up = (vrna_hc_up_t *)vrna_alloc(sizeof(vrna_hc_up_t) * size_up);
1510 bp = (struct hc_bp *)vrna_alloc(sizeof(struct hc_bp) * size_bp);
1511 bp_unspecific = (struct hc_bp *)vrna_alloc(sizeof(struct hc_bp) * size_bp_unspecific);
1512
1513
1514 for (hx = 0, j = 1; j <= n; j++) {
1515 switch (constraint[j - 1]) {
1516 /* can't pair */
1517 case 'x':
1518 if (options & VRNA_CONSTRAINT_DB_X) {
1519 up[num_up].position = j;
1520 up[num_up].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1521
1522 num_up++;
1523
1524 if (num_up == size_up) {
1525 size_up *= 1.4;
1526 up = (vrna_hc_up_t *)vrna_realloc(up, sizeof(vrna_hc_up_t) * size_up);
1527 }
1528 }
1529
1530 break;
1531
1532 /* must pair, i.e. may not be unpaired */
1533 case '|':
1534 if (options & VRNA_CONSTRAINT_DB_PIPE) {
1535 bp_unspecific[num_bp_unspecific].i = j; /* position */
1536 bp_unspecific[num_bp_unspecific].j = 0; /* direction */
1537 bp_unspecific[num_bp_unspecific].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1538
1539 num_bp_unspecific++;
1540
1541 if (num_bp_unspecific == size_bp_unspecific) {
1542 size_bp_unspecific *= 1.4;
1543 bp_unspecific = (struct hc_bp *)vrna_realloc(bp_unspecific,
1544 sizeof(struct hc_bp) *
1545 size_bp_unspecific);
1546 }
1547 }
1548
1549 break;
1550
1551 /* weak enforced pair 'open' */
1552 case '(':
1553 if (options & VRNA_CONSTRAINT_DB_RND_BRACK)
1554 stack[hx++] = j;
1555
1556 break;
1557
1558 /* weak enforced pair 'close' */
1559 case ')':
1560 if (options & VRNA_CONSTRAINT_DB_RND_BRACK) {
1561 if (hx <= 0) {
1562 vrna_message_warning("vrna_hc_add_from_db: "
1563 "Unbalanced brackets in constraint string\n%s\n"
1564 "No constraints will be applied!",
1565 constraint);
1566 goto db_constraints_exit;
1567 }
1568
1569 i = stack[--hx];
1570
1571 if (options & VRNA_CONSTRAINT_DB_CANONICAL_BP) {
1572 /* check whether this pair forms a non-canoncial base pair */
1573 if (md->pair[S[i]][S[j]] == 0) {
1574 vrna_message_warning("Removing non-canonical base pair %c%c (%d,%d) from constraint",
1575 sequence[i - 1], sequence[j - 1],
1576 i, j);
1577 break;
1578 }
1579 }
1580
1581 if ((j - i - 1) < md->min_loop_size) {
1582 vrna_message_warning("vrna_hc_add_from_db: "
1583 "Pairing partners (%d, %d) violate minimum loop size settings of %dnt, omitting constraint",
1584 i,
1585 j,
1586 md->min_loop_size);
1587 break;
1588 }
1589
1590 bp[num_bp].i = i;
1591 bp[num_bp].j = j;
1592 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1593
1594 if (options & VRNA_CONSTRAINT_DB_ENFORCE_BP)
1595 bp[num_bp].options |= VRNA_CONSTRAINT_CONTEXT_ENFORCE;
1596
1597 num_bp++;
1598
1599 if (num_bp == size_bp) {
1600 size_bp *= 1.4;
1601 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1602 }
1603 }
1604
1605 break;
1606
1607 /* pairs downstream */
1608 case '<':
1609 if (options & VRNA_CONSTRAINT_DB_ANG_BRACK) {
1610 bp_unspecific[num_bp_unspecific].i = j;
1611 bp_unspecific[num_bp_unspecific].j = 1;
1612 bp_unspecific[num_bp_unspecific].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1613
1614 num_bp_unspecific++;
1615
1616 if (num_bp_unspecific == size_bp_unspecific) {
1617 size_bp_unspecific *= 1.4;
1618 bp_unspecific = (struct hc_bp *)vrna_realloc(bp_unspecific,
1619 sizeof(struct hc_bp) *
1620 size_bp_unspecific);
1621 }
1622
1623 if (!(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)) {
1624 /* (re-)allow this nucleotide to stay unpaired for nostalgic reasons */
1625 up[num_up].position = j;
1626 up[num_up].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS |
1627 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1628
1629 num_up++;
1630
1631 if (num_up == size_up) {
1632 size_up *= 1.4;
1633 up = (vrna_hc_up_t *)vrna_realloc(up, sizeof(vrna_hc_up_t) * size_up);
1634 }
1635 }
1636 }
1637
1638 break;
1639
1640 /* pairs upstream */
1641 case '>':
1642 if (options & VRNA_CONSTRAINT_DB_ANG_BRACK) {
1643 bp_unspecific[num_bp_unspecific].i = j;
1644 bp_unspecific[num_bp_unspecific].j = -1;
1645 bp_unspecific[num_bp_unspecific].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
1646
1647 num_bp_unspecific++;
1648
1649 if (num_bp_unspecific == size_bp_unspecific) {
1650 size_bp_unspecific *= 1.4;
1651 bp_unspecific = (struct hc_bp *)vrna_realloc(bp_unspecific,
1652 sizeof(struct hc_bp) *
1653 size_bp_unspecific);
1654 }
1655
1656 if (!(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)) {
1657 /* (re-)allow this nucleotide to stay unpaired for nostalgic reasons */
1658 up[num_up].position = j;
1659 up[num_up].options = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS |
1660 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1661
1662 num_up++;
1663
1664 if (num_up == size_up) {
1665 size_up *= 1.4;
1666 up = (vrna_hc_up_t *)vrna_realloc(up, sizeof(vrna_hc_up_t) * size_up);
1667 }
1668 }
1669 }
1670
1671 break;
1672
1673 /* only intramolecular basepairing */
1674 case 'l':
1675 if (options & VRNA_CONSTRAINT_DB_INTRAMOL) {
1676 unsigned int l;
1677 if (cut > 1) {
1678 if (j < cut) {
1679 for (l = MAX2(j + min_loop_size, cut); l <= length; l++) {
1680 bp[num_bp].i = j;
1681 bp[num_bp].j = l;
1682 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1683 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1684
1685 num_bp++;
1686
1687 if (num_bp == size_bp) {
1688 size_bp *= 1.4;
1689 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1690 }
1691 }
1692 } else {
1693 for (l = 1; l < MIN2(cut, j - min_loop_size); l++) {
1694 bp[num_bp].i = l;
1695 bp[num_bp].j = j;
1696 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1697 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1698
1699 num_bp++;
1700
1701 if (num_bp == size_bp) {
1702 size_bp *= 1.4;
1703 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1704 }
1705 }
1706 }
1707 }
1708 }
1709
1710 break;
1711
1712 /* only intermolecular bp */
1713 case 'e':
1714 if (options & VRNA_CONSTRAINT_DB_INTERMOL) {
1715 unsigned int l;
1716 if (cut > 1) {
1717 if (j < cut) {
1718 for (l = 1; l < j; l++) {
1719 bp[num_bp].i = l;
1720 bp[num_bp].j = j;
1721 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1722 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1723
1724 num_bp++;
1725
1726 if (num_bp == size_bp) {
1727 size_bp *= 1.4;
1728 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1729 }
1730 }
1731
1732 for (l = j + 1; l < cut; l++) {
1733 bp[num_bp].i = j;
1734 bp[num_bp].j = l;
1735 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1736 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1737
1738 num_bp++;
1739
1740 if (num_bp == size_bp) {
1741 size_bp *= 1.4;
1742 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1743 }
1744 }
1745 } else {
1746 for (l = cut; l < j; l++) {
1747 bp[num_bp].i = l;
1748 bp[num_bp].j = j;
1749 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1750 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1751
1752 num_bp++;
1753
1754 if (num_bp == size_bp) {
1755 size_bp *= 1.4;
1756 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1757 }
1758 }
1759
1760 for (l = j + 1; l <= length; l++) {
1761 bp[num_bp].i = j;
1762 bp[num_bp].j = l;
1763 bp[num_bp].options = VRNA_CONSTRAINT_CONTEXT_NONE |
1764 VRNA_CONSTRAINT_CONTEXT_NO_REMOVE;
1765
1766 num_bp++;
1767
1768 if (num_bp == size_bp) {
1769 size_bp *= 1.4;
1770 bp = (struct hc_bp *)vrna_realloc(bp, sizeof(struct hc_bp) * size_bp);
1771 }
1772 }
1773 }
1774 }
1775 }
1776
1777 break;
1778
1779 case '.':
1780 break;
1781
1782 default:
1783 vrna_message_warning(
1784 "vrna_hc_add_from_db: "
1785 "Unrecognized character '%c' in constraint string",
1786 constraint[j - 1]);
1787 break;
1788 }
1789 }
1790
1791 if (hx != 0) {
1792 vrna_message_warning("vrna_hc_add_from_db: "
1793 "Unbalanced brackets in constraint string\n%s\n"
1794 "No constraints will be applied!",
1795 constraint);
1796 goto db_constraints_exit;
1797 }
1798
1799 /* finally, apply constraints */
1800
1801 /* 1st, unspecific pairing states */
1802 for (i = 0; i < num_bp_unspecific; i++)
1803 vrna_hc_add_bp_nonspecific(vc,
1804 bp_unspecific[i].i, /* nucleotide position */
1805 bp_unspecific[i].j, /* pairing direction */
1806 bp_unspecific[i].options);
1807
1808 /* 2nd, specific base pairs */
1809 for (i = 0; i < num_bp; i++)
1810 vrna_hc_add_bp(vc,
1811 bp[i].i,
1812 bp[i].j,
1813 bp[i].options);
1814
1815 /* 3rd, unpaired constraints */
1816 if (num_up > 0) {
1817 up[num_up].position = 0; /* end of list marker */
1818 vrna_hc_add_up_batch(vc, up);
1819 }
1820
1821 db_constraints_exit:
1822
1823 /* clean up */
1824 free(up);
1825 free(bp);
1826 free(bp_unspecific);
1827 free(stack);
1828 }
1829
1830
1831 PRIVATE void
hc_reset_to_default(vrna_fold_compound_t * vc)1832 hc_reset_to_default(vrna_fold_compound_t *vc)
1833 {
1834 vrna_hc_t *hc = vc->hc;
1835
1836 /* ######################### */
1837 /* fill with default values */
1838 /* ######################### */
1839
1840 default_hc_up(vc, 0);
1841
1842 default_hc_bp(vc, 0);
1843
1844 /* should we reset the generalized hard constraint feature here? */
1845 if (hc->f || hc->data) {
1846 if (hc->free_data)
1847 hc->free_data(hc->data);
1848
1849 hc->f = NULL;
1850 hc->data = NULL;
1851 hc->free_data = NULL;
1852 }
1853 }
1854
1855
1856 PRIVATE void
hc_update_up(vrna_fold_compound_t * vc)1857 hc_update_up(vrna_fold_compound_t *vc)
1858 {
1859 unsigned int i, n, u;
1860 vrna_hc_t *hc;
1861
1862 n = vc->length;
1863 hc = vc->hc;
1864
1865 if (hc->type == VRNA_HC_WINDOW) {
1866 /* do nothing for now! */
1867 } else {
1868 for (hc->up_ext[n + 1] = 0, i = n; i > 0; i--) /* unpaired stretch in exterior loop */
1869 hc->up_ext[i] = (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) ? 1 +
1870 hc->up_ext[i + 1] : 0;
1871
1872 for (hc->up_hp[n + 1] = 0, i = n; i > 0; i--) /* unpaired stretch in hairpin loop */
1873 hc->up_hp[i] = (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) ? 1 +
1874 hc->up_hp[i + 1] : 0;
1875
1876 for (hc->up_int[n + 1] = 0, i = n; i > 0; i--) /* unpaired stretch in interior loop */
1877 hc->up_int[i] = (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) ? 1 +
1878 hc->up_int[i + 1] : 0;
1879
1880 for (hc->up_ml[n + 1] = 0, i = n; i > 0; i--) /* unpaired stretch in multibranch loop */
1881 hc->up_ml[i] = (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) ? 1 +
1882 hc->up_ml[i + 1] : 0;
1883
1884 /*
1885 * loop arround once more until we find a nucleotide that mustn't
1886 * be unpaired (needed for circular folding)
1887 * Note, circular fold is only possible for single strand predictions
1888 */
1889 if (vc->strands < 2) {
1890 if (hc->mx[n + 1] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) {
1891 hc->up_ext[n + 1] = hc->up_ext[1];
1892 for (i = n; i > 0; i--) {
1893 if (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP)
1894 hc->up_ext[i] = MIN2(n, 1 + hc->up_ext[i + 1]);
1895 else
1896 break;
1897 }
1898 }
1899
1900 if (hc->mx[n + 1] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) {
1901 hc->up_hp[n + 1] = hc->up_hp[1];
1902 for (i = n; i > 0; i--) {
1903 if (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP)
1904 hc->up_hp[i] = MIN2(n, 1 + hc->up_hp[i + 1]);
1905 else
1906 break;
1907 }
1908 }
1909
1910 if (hc->mx[n + 1] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) {
1911 hc->up_int[n + 1] = hc->up_int[1];
1912 for (i = n; i > 0; i--) {
1913 if (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP)
1914 hc->up_int[i] = MIN2(n, 1 + hc->up_int[i + 1]);
1915 else
1916 break;
1917 }
1918 }
1919
1920 if (hc->mx[n + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) {
1921 hc->up_ml[n + 1] = hc->up_ml[1];
1922 for (i = n; i > 0; i--) {
1923 if (hc->mx[n * i + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP)
1924 hc->up_ml[i] = MIN2(n, 1 + hc->up_ml[i + 1]);
1925 else
1926 break;
1927 }
1928 }
1929 }
1930 }
1931 }
1932
1933
1934 PRIVATE void
hc_update_up_window(vrna_fold_compound_t * vc,int i,unsigned int options)1935 hc_update_up_window(vrna_fold_compound_t *vc,
1936 int i,
1937 unsigned int options)
1938 {
1939 vrna_hc_t *hc;
1940 int k, winsize, up_ext, up_hp, up_int, up_ml;
1941
1942 hc = vc->hc;
1943 winsize = vc->window_size;
1944
1945 if (options & VRNA_CONSTRAINT_WINDOW_UPDATE_5) {
1946 up_ext = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) ?
1947 1 : 0;
1948 up_hp = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) ?
1949 1 : 0;
1950 up_int = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) ?
1951 1 : 0;
1952 up_ml = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) ?
1953 1 : 0;
1954 } else {
1955 up_ext = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) ?
1956 1 + hc->up_ext[i + 1] :
1957 0;
1958 up_hp = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) ?
1959 1 + hc->up_hp[i + 1] :
1960 0;
1961 up_int = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) ?
1962 1 + hc->up_int[i + 1] :
1963 0;
1964 up_ml = (hc->matrix_local[i][0] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) ?
1965 1 + hc->up_ml[i + 1] :
1966 0;
1967 }
1968
1969 hc->up_ext[i] = up_ext;
1970 hc->up_hp[i] = up_hp;
1971 hc->up_int[i] = up_int;
1972 hc->up_ml[i] = up_ml;
1973
1974 if (options & VRNA_CONSTRAINT_WINDOW_UPDATE_5) {
1975 /* the sliding window proceeds from 5' to 3' so we update constraints 3' to 5' */
1976 if (up_ext > 0)
1977 for (k = i - 1; k >= MAX2(1, i - winsize); k--) {
1978 if (hc->up_ext[k] < 1)
1979 break;
1980 hc->up_ext[k] += up_ext;
1981 }
1982
1983 if (up_hp > 0)
1984 for (k = i - 1; k >= MAX2(1, i - winsize); k--) {
1985 if (hc->up_hp[k] < 1)
1986 break;
1987 hc->up_hp[k] += up_hp;
1988 }
1989
1990 if (up_int > 0)
1991 for (k = i - 1; k >= MAX2(1, i - winsize); k--) {
1992 if (hc->up_int[k] < 1)
1993 break;
1994 hc->up_int[k] += up_int;
1995 }
1996
1997 if (up_ml > 0)
1998 for (k = i - 1; k >= MAX2(1, i - winsize); k--) {
1999 if (hc->up_ml[k] < 1)
2000 break;
2001 hc->up_ml[k] += up_ml;
2002 }
2003 }
2004 }
2005
2006
2007 #ifndef VRNA_DISABLE_BACKWARD_COMPATIBILITY
2008
2009 /*###########################################*/
2010 /*# deprecated functions below #*/
2011 /*###########################################*/
2012
2013 PUBLIC void
print_tty_constraint_full(void)2014 print_tty_constraint_full(void)
2015 {
2016 vrna_message_constraint_options_all();
2017 }
2018
2019
2020 PUBLIC void
print_tty_constraint(unsigned int option)2021 print_tty_constraint(unsigned int option)
2022 {
2023 vrna_message_constraint_options(option);
2024 }
2025
2026
2027 PUBLIC void
constrain_ptypes(const char * constraint,unsigned int length,char * ptype,int * BP,int min_loop_size,unsigned int idx_type)2028 constrain_ptypes(const char *constraint,
2029 unsigned int length,
2030 char *ptype,
2031 int *BP,
2032 int min_loop_size,
2033 unsigned int idx_type)
2034 {
2035 int n, i, j, k, l;
2036 int hx, *stack;
2037 char type;
2038 int *index;
2039
2040 if (constraint == NULL)
2041 return;
2042
2043 n = (int)strlen(constraint);
2044
2045 stack = vrna_alloc(sizeof(int) * (n + 1));
2046
2047 if (!idx_type) {
2048 /* index allows access in energy matrices at pos (i,j) via index[j]+i */
2049 index = vrna_idx_col_wise(length);
2050
2051 for (hx = 0, j = 1; j <= n; j++) {
2052 switch (constraint[j - 1]) {
2053 case '|':
2054 if (BP)
2055 BP[j] = -1;
2056
2057 break;
2058 case 'x': /* can't pair */
2059 for (l = 1; l < j - min_loop_size; l++)
2060 ptype[index[j] + l] = 0;
2061 for (l = j + min_loop_size + 1; l <= (int)length; l++)
2062 ptype[index[l] + j] = 0;
2063 break;
2064 case '(':
2065 stack[hx++] = j;
2066 /* fallthrough */
2067 case '<': /* pairs upstream */
2068 for (l = 1; l < j - min_loop_size; l++)
2069 ptype[index[j] + l] = 0;
2070 break;
2071 case ')':
2072 if (hx <= 0)
2073 vrna_message_error("%s\nunbalanced brackets in constraint", constraint);
2074
2075 i = stack[--hx];
2076 type = ptype[index[j] + i];
2077 for (k = i + 1; k <= (int)length; k++)
2078 ptype[index[k] + i] = 0;
2079 /* don't allow pairs i<k<j<l */
2080 for (l = j; l <= (int)length; l++)
2081 for (k = i + 1; k <= j; k++)
2082 ptype[index[l] + k] = 0;
2083 /* don't allow pairs k<i<l<j */
2084 for (l = i; l <= j; l++)
2085 for (k = 1; k <= i; k++)
2086 ptype[index[l] + k] = 0;
2087 for (k = 1; k < j; k++)
2088 ptype[index[j] + k] = 0;
2089 ptype[index[j] + i] = (type == 0) ? 7 : type;
2090 /* fallthrough */
2091 case '>': /* pairs downstream */
2092 for (l = j + min_loop_size + 1; l <= (int)length; l++)
2093 ptype[index[l] + j] = 0;
2094 break;
2095 }
2096 }
2097 } else {
2098 /* index allows access in energy matrices at pos (i,j) via index[i]-j */
2099 index = vrna_idx_row_wise(length);
2100
2101 for (hx = 0, j = 1; j <= n; j++) {
2102 switch (constraint[j - 1]) {
2103 case 'x': /* can't pair */
2104 for (l = 1; l < j - min_loop_size; l++)
2105 ptype[index[l] - j] = 0;
2106 for (l = j + min_loop_size + 1; l <= (int)length; l++)
2107 ptype[index[j] - l] = 0;
2108 break;
2109 case '(':
2110 stack[hx++] = j;
2111 /* fallthrough */
2112 case '<': /* pairs upstream */
2113 for (l = 1; l < j - min_loop_size; l++)
2114 ptype[index[l] - j] = 0;
2115 break;
2116 case ')':
2117 if (hx <= 0)
2118 vrna_message_error("%s\nunbalanced brackets in constraints", constraint);
2119
2120 i = stack[--hx];
2121 type = ptype[index[i] - j];
2122 /* don't allow pairs i<k<j<l */
2123 for (k = i; k <= j; k++)
2124 for (l = j; l <= (int)length; l++)
2125 ptype[index[k] - l] = 0;
2126 /* don't allow pairs k<i<l<j */
2127 for (k = 1; k <= i; k++)
2128 for (l = i; l <= j; l++)
2129 ptype[index[k] - l] = 0;
2130 ptype[index[i] - j] = (type == 0) ? 7 : type;
2131 /* fallthrough */
2132 case '>': /* pairs downstream */
2133 for (l = j + min_loop_size + 1; l <= (int)length; l++)
2134 ptype[index[j] - l] = 0;
2135 break;
2136 }
2137 }
2138 }
2139
2140 if (hx != 0)
2141 vrna_message_error("%s\nunbalanced brackets in constraint string", constraint);
2142
2143 free(index);
2144 free(stack);
2145 }
2146
2147
2148 #endif
2149