ViennaRNA-bindings-0.233.1.1: C/ViennaRNA/constraints_hard.c
/* constraints handling */
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#include <limits.h>
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/energy_const.h" /* defines MINPSCORE */
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/utils.h"
#include "ViennaRNA/aln_util.h"
#include "ViennaRNA/file_formats.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/constraints.h"
#include "ViennaRNA/constraints_hard.h"
#ifdef __GNUC__
# define INLINE inline
#else
# define INLINE
#endif
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE void
hc_add_up(vrna_fold_compound_t *vc,
int i,
char option);
PRIVATE INLINE void
hc_cant_pair( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index);
PRIVATE INLINE void
hc_must_pair( unsigned int i,
char c_option,
char *hc,
int *index);
PRIVATE INLINE void
hc_pairs_upstream(unsigned int i,
char c_option,
char *hc,
unsigned int length,
int *index);
PRIVATE INLINE void
hc_pairs_downstream(unsigned int i,
char c_option,
char *hc,
unsigned int length,
int *index);
PRIVATE INLINE void
hc_allow_pair(unsigned int i,
unsigned int j,
char c_option,
char *hc,
int *index);
PRIVATE INLINE void
hc_weak_enforce_pair( unsigned int i,
unsigned int j,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index);
PRIVATE INLINE void
hc_enforce_pair(unsigned int i,
unsigned int j,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index);
PRIVATE INLINE void
hc_intramolecular_only( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int cut,
int *index);
PRIVATE INLINE void
hc_intermolecular_only( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int cut,
int *index);
PRIVATE void
apply_DB_constraint(const char *constraint,
char *ptype,
unsigned int length,
unsigned int min_loop_size,
int cut,
unsigned int options);
PRIVATE void
hc_reset_to_default(vrna_fold_compound_t *vc);
PRIVATE void
hc_update_up(vrna_fold_compound_t *vc);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
PUBLIC void
vrna_message_constraint_options_all(void){
vrna_message_constraint_options( VRNA_CONSTRAINT_DB_PIPE
| VRNA_CONSTRAINT_DB_DOT
| VRNA_CONSTRAINT_DB_X
| VRNA_CONSTRAINT_DB_ANG_BRACK
| VRNA_CONSTRAINT_DB_RND_BRACK);
}
PUBLIC void
vrna_message_constraint_options(unsigned int option){
printf("Input structure constraints using the following notation:\n");
if(option & VRNA_CONSTRAINT_DB_PIPE) printf("| : paired with another base\n");
if(option & VRNA_CONSTRAINT_DB_DOT) printf(". : no constraint at all\n");
if(option & VRNA_CONSTRAINT_DB_X) printf("x : base must not pair\n");
if(option & VRNA_CONSTRAINT_DB_ANG_BRACK) printf("< : base i is paired with a base j<i\n> : base i is paired with a base j>i\n");
if(option & VRNA_CONSTRAINT_DB_RND_BRACK) printf("matching brackets ( ): base i pairs base j\n");
}
PUBLIC void
vrna_hc_init(vrna_fold_compound_t *vc){
unsigned int n;
vrna_hc_t *hc;
n = vc->length;
/* free previous hard constraints */
vrna_hc_free(vc->hc);
/* allocate memory new hard constraints data structure */
hc = (vrna_hc_t *)vrna_alloc(sizeof(vrna_hc_t));
hc->matrix = (char *)vrna_alloc(sizeof(char)*((n*(n+1))/2+2));
hc->up_ext = (int *)vrna_alloc(sizeof(int)*(n+2));
hc->up_hp = (int *)vrna_alloc(sizeof(int)*(n+2));
hc->up_int = (int *)vrna_alloc(sizeof(int)*(n+2));
hc->up_ml = (int *)vrna_alloc(sizeof(int)*(n+2));
/* set new hard constraints */
vc->hc = hc;
/* prefill default values */
hc_reset_to_default(vc);
/* add null pointers for the generalized hard constraint feature */
hc->f = NULL;
hc->data = NULL;
hc->free_data = NULL;
/* update */
hc_update_up(vc);
}
PUBLIC void
vrna_hc_add_up( vrna_fold_compound_t *vc,
int i,
char option){
int j;
if(vc)
if(vc->hc){
if((i <= 0) || (i > vc->length)){
vrna_message_warning("vrna_hc_add_up: position out of range, not doing anything");
return;
}
hc_add_up(vc, i, option);
hc_update_up(vc);
}
}
PUBLIC int
vrna_hc_add_up_batch( vrna_fold_compound_t *vc,
vrna_hc_up_t *constraints){
int i, ret;
ret = 0; /* failure */
if(vc)
if(vc->hc && constraints){
for(i = 0; constraints[i].position != 0; i++){
int pos = constraints[i].position;
char options = constraints[i].options;
if((pos <= 0) || (pos > vc->length)){
vrna_message_warning("vrna_hc_add_up_batch: position out of range, application of hard constraints stops here!");
return ret;
}
hc_add_up(vc, pos, options);
}
hc_update_up(vc);
ret = 1; /* success */
}
return ret;
}
PRIVATE void
hc_add_up(vrna_fold_compound_t *vc,
int i,
char option){
int j;
char type = (char)0;
if(option & VRNA_CONSTRAINT_CONTEXT_ENFORCE){ /* force nucleotide to appear unpaired within a certain type of loop */
/* do not allow i to be paired with any other nucleotide */
if(!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)){
for(j = 1; j < i; j++)
vc->hc->matrix[vc->jindx[i] + j] = (char)0;
for(j = i+1; j <= vc->length; j++)
vc->hc->matrix[vc->jindx[j] + i] = (char)0;
}
type = option & (char)( VRNA_CONSTRAINT_CONTEXT_EXT_LOOP
| VRNA_CONSTRAINT_CONTEXT_HP_LOOP
| VRNA_CONSTRAINT_CONTEXT_INT_LOOP
| VRNA_CONSTRAINT_CONTEXT_MB_LOOP);
vc->hc->matrix[vc->jindx[i] + i] = type;
} else {
type = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
/* do not allow i to be paired with any other nucleotide (in context type) */
if(!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)){
for(j = 1; j < i; j++)
vc->hc->matrix[vc->jindx[i] + j] &= ~type;
for(j = i+1; j <= vc->length; j++)
vc->hc->matrix[vc->jindx[j] + i] &= ~type;
}
vc->hc->matrix[vc->jindx[i] + i] = (char)( VRNA_CONSTRAINT_CONTEXT_EXT_LOOP
| VRNA_CONSTRAINT_CONTEXT_HP_LOOP
| VRNA_CONSTRAINT_CONTEXT_INT_LOOP
| VRNA_CONSTRAINT_CONTEXT_MB_LOOP);
}
}
PUBLIC void
vrna_hc_add_bp_nonspecific( vrna_fold_compound_t *vc,
int i,
int d,
char option){
int p;
char type, t1, t2;
if(vc)
if(vc->hc){
if((i <= 0) || (i > vc->length)){
vrna_message_warning("vrna_hc_add_bp_nonspecific: position out of range, not doing anything");
return;
}
/* position i may pair in provided contexts */
type = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
/* acknowledge pairing direction */
t1 = (d <= 0) ? type : (char)0;
t2 = (d >= 0) ? type : (char)0;
if(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE){
/* only allow for possibly non-canonical pairs, do not enforce them */
for(p = 1; p < i; p++)
vc->hc->matrix[vc->jindx[i] + p] |= t1;
for(p = i+1; p <= vc->length; p++)
vc->hc->matrix[vc->jindx[p] + i] |= t2;
} else {
/* force pairing direction */
for(p = 1; p < i; p++)
vc->hc->matrix[vc->jindx[i] + p] &= t1;
for(p = i+1; p <= vc->length; p++)
vc->hc->matrix[vc->jindx[p] + i] &= t2;
/* nucleotide mustn't be unpaired */
vc->hc->matrix[vc->jindx[i] + i] = (char)0;
}
hc_update_up(vc);
}
}
PUBLIC void
vrna_hc_add_bp( vrna_fold_compound_t *vc,
int i,
int j,
char option){
int k, l;
char type;
if(vc)
if(vc->hc){
if((i <= 0) || (j <= i) || (j > vc->length)){
vrna_message_warning("vrna_hc_add_bp: position out of range, not doing anything");
return;
}
/* reset ptype in case (i,j) is a non-canonical pair */
if(option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS){
if(vc->hc->matrix[vc->jindx[j] + i])
if(vc->ptype[vc->jindx[j] + i] == 0)
vc->ptype[vc->jindx[j] + i] = 7;
}
vc->hc->matrix[vc->jindx[j] + i] = option & VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
if(!(option & VRNA_CONSTRAINT_CONTEXT_NO_REMOVE)){
/*
remove all conflicting base pairs, i.e. do not allow i,j to pair
with any other nucleotide k
*/
for(k = 1; k < i; k++){
vc->hc->matrix[vc->jindx[i] + k] = (char)0;
vc->hc->matrix[vc->jindx[j] + k] = (char)0;
for(l = i+1; l < j; l++)
vc->hc->matrix[vc->jindx[l] + k] = (char)0;
}
for(k = i+1; k < j; k++){
vc->hc->matrix[vc->jindx[k] + i] = (char)0;
vc->hc->matrix[vc->jindx[j] + k] = (char)0;
for(l = j + 1; l <= vc->length; l++)
vc->hc->matrix[vc->jindx[l] + k] = (char)0;
}
for(k = j+1; k <= vc->length; k++){
vc->hc->matrix[vc->jindx[k] + i] = (char)0;
vc->hc->matrix[vc->jindx[k] + j] = (char)0;
}
}
if(option & VRNA_CONSTRAINT_CONTEXT_ENFORCE){
/* do not allow i,j to be unpaired */
vc->hc->matrix[vc->jindx[i] + i] = (char)0;
vc->hc->matrix[vc->jindx[j] + j] = (char)0;
hc_update_up(vc);
}
}
}
PUBLIC void
vrna_hc_free(vrna_hc_t *hc){
if(hc){
free(hc->matrix);
free(hc->up_ext);
free(hc->up_hp);
free(hc->up_int);
free(hc->up_ml);
if(hc->free_data)
hc->free_data(hc->data);
free(hc);
}
}
PUBLIC void
vrna_hc_add_f(vrna_fold_compound_t *vc,
vrna_callback_hc_evaluate *f)
{
if (vc && f) {
if (vc->type == VRNA_FC_TYPE_SINGLE) {
if (!vc->hc)
vrna_hc_init(vc);
vc->hc->f = f;
}
}
}
PUBLIC void
vrna_hc_add_data( vrna_fold_compound_t *vc,
void *data,
vrna_callback_free_auxdata *f)
{
if (vc && data) {
if (vc->type == VRNA_FC_TYPE_SINGLE) {
if (!vc->hc)
vrna_hc_init(vc);
vc->hc->data = data;
vc->hc->free_data = f;
}
}
}
PUBLIC int
vrna_hc_add_from_db(vrna_fold_compound_t *vc,
const char *constraint,
unsigned int options){
int i, d, ret;
vrna_md_t *md;
ret = 0; /* Failure */
if(vc){
if(vc->params)
md = &(vc->params->model_details);
else if(vc->exp_params)
md = &(vc->exp_params->model_details);
else
return ret;
if(!vc->hc)
vrna_hc_init(vc);
/* apply hard constraints from dot-bracket notation */
apply_DB_constraint(constraint,
vc->hc->matrix,
vc->length,
md->min_loop_size,
-1,
options);
hc_update_up(vc);
ret = 1; /* Success */
}
return ret;
}
PRIVATE void
apply_DB_constraint(const char *constraint,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int cut,
unsigned int options){
int n,i,j;
int hx, *stack;
int *index;
char c_option;
if(constraint == NULL) return;
n = (int)strlen(constraint);
stack = (int *) vrna_alloc(sizeof(int)*(n+1));
index = vrna_idx_col_wise(length);
c_option = VRNA_CONSTRAINT_CONTEXT_EXT_LOOP
| VRNA_CONSTRAINT_CONTEXT_HP_LOOP
| VRNA_CONSTRAINT_CONTEXT_INT_LOOP
| VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC
| VRNA_CONSTRAINT_CONTEXT_MB_LOOP
| VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC;
for(hx=0, j=1; j<=n; j++) {
switch (constraint[j-1]) {
/* can't pair */
case 'x': if(options & VRNA_CONSTRAINT_DB_X){
hc_cant_pair(j, c_option, hc, length, min_loop_size, index);
}
break;
/* must pair, i.e. may not be unpaired */
case '|': if(options & VRNA_CONSTRAINT_DB_PIPE){
if(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)
hc_must_pair(j, c_option, hc, index);
}
break;
/* weak enforced pair 'open' */
case '(': if(options & VRNA_CONSTRAINT_DB_RND_BRACK){
stack[hx++]=j;
}
break;
/* weak enforced pair 'close' */
case ')': if(options & VRNA_CONSTRAINT_DB_RND_BRACK){
if (hx<=0) {
vrna_message_error("%s\nunbalanced brackets in constraints", constraint);
}
i = stack[--hx];
if(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)
hc_enforce_pair(i, j, c_option, hc, length, min_loop_size, index);
else
hc_weak_enforce_pair(i, j, c_option, hc, length, min_loop_size, index);
}
break;
/* pairs upstream */
case '<': if(options & VRNA_CONSTRAINT_DB_ANG_BRACK){
hc_pairs_downstream(j, c_option, hc, length, index);
if(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)
hc_must_pair(j, c_option, hc, index);
}
break;
/* pairs downstream */
case '>': if(options & VRNA_CONSTRAINT_DB_ANG_BRACK){
hc_pairs_upstream(j, c_option, hc, length, index);
if(options & VRNA_CONSTRAINT_DB_ENFORCE_BP)
hc_must_pair(j, c_option, hc, index);
}
break;
/* only intramolecular basepairing */
case 'l': if(options & VRNA_CONSTRAINT_DB_INTRAMOL){
hc_intramolecular_only(j, c_option, hc, length, min_loop_size, cut, index);
}
break;
/* only intermolecular bp */
case 'e': if(options & VRNA_CONSTRAINT_DB_INTERMOL){
hc_intermolecular_only(j, c_option, hc, length, min_loop_size, cut, index);
}
break;
case '.': break;
default: vrna_message_warning("Unrecognized character '%c' in pseudo dot-bracket notation constraint string",
constraint[j-1]);
break;
}
}
if (hx!=0) {
vrna_message_error("%s\nunbalanced brackets in constraint string", constraint);
}
/* clean up */
free(index);
free(stack);
}
PRIVATE INLINE void
hc_intramolecular_only( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int cut,
int *index){
unsigned int l;
if(cut > 1){
if(i < cut)
for(l = MAX2(i+min_loop_size, cut); l <= length; l++)
hc[index[l] + i] &= ~c_option;
else
for(l = 1; l < MIN2(cut, i-min_loop_size); l++)
hc[index[i] + l] &= ~c_option;
}
}
PRIVATE INLINE void
hc_intermolecular_only( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int cut,
int *index){
unsigned int l;
if(cut > 1){
if(i < cut){
for(l = 1; l < i; l++)
hc[index[i] + l] &= ~c_option;
for(l = i + 1; l < cut; l++)
hc[index[l] + i] &= ~c_option;
} else {
for(l = cut; l < i; l++)
hc[index[i] + l] &= ~c_option;
for(l = i + 1; l <= length; l++)
hc[index[l] + i] &= ~c_option;
}
}
}
PRIVATE INLINE void
hc_cant_pair( unsigned int i,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index){
hc_pairs_upstream(i, c_option, hc, length, index);
hc_pairs_downstream(i, c_option, hc, length, index);
}
PRIVATE INLINE void
hc_must_pair( unsigned int i,
char c_option,
char *hc,
int *index){
hc[index[i]+i] &= ~c_option;
}
PRIVATE INLINE void
hc_pairs_upstream(unsigned int i,
char c_option,
char *hc,
unsigned int length,
int *index){
unsigned int l;
/* prohibit downstream pairs */
for(l = length; l > i; l--)
hc[index[l] + i] = (char)0;
/* allow upstream pairs of given type */
for(l = i - 1; l >= 1; l--)
hc[index[i] + l] &= c_option;
}
PRIVATE INLINE void
hc_pairs_downstream(unsigned int i,
char c_option,
char *hc,
unsigned int length,
int *index){
unsigned int l;
/* allow downstream pairs of given type */
for(l = length; l > i; l--)
hc[index[l] + i] &= c_option;
/* forbid upstream pairs */
for(l = i - 1; l >= 1; l--)
hc[index[i] + l] = (char)0;
}
PRIVATE INLINE void
hc_allow_pair(unsigned int i,
unsigned int j,
char c_option,
char *hc,
int *index){
hc[index[j] + i] |= c_option;
}
PRIVATE INLINE void
hc_weak_enforce_pair( unsigned int i,
unsigned int j,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index){
unsigned int k, l;
/* don't allow pairs (k,i) 1 <= k < i */
/* don't allow pairs (i,k) i < k <= n */
hc_pairs_upstream(i, (char)0, hc, length, index);
/* don't allow pairs (k,j) 1 <= k < j */
/* don't allow pairs (j,k) j < k <= n */
hc_pairs_upstream(j, (char)0, hc, length, index);
/* don't allow pairs i < k < j < l */
for(k = i+1; k < j; k++)
for(l = j+1; l <= length; l++){
hc[index[l] + k] = 0;
}
/* don't allow pairs k<i<l<j */
for(k = 1; k < i; k++)
for(l = i+1; l < j; l++){
hc[index[l] + k] = 0;
}
/* allow base pair (i,j) */
hc[index[j] + i] |= c_option;
}
PRIVATE INLINE void
hc_enforce_pair(unsigned int i,
unsigned int j,
char c_option,
char *hc,
unsigned int length,
unsigned int min_loop_size,
int *index){
hc_weak_enforce_pair( i,
j,
c_option,
hc,
length,
min_loop_size,
index);
/* forbid i and j to be unpaired */
hc[index[i] + i] = 0;
hc[index[j] + j] = 0;
}
PRIVATE void
hc_reset_to_default(vrna_fold_compound_t *vc){
unsigned int i, j, ij, min_loop_size, n;
int max_span, *idx;
vrna_md_t *md;
vrna_hc_t *hc;
short *S;
md = NULL;
n = vc->length;
hc = vc->hc;
idx = vc->jindx;
S = vc->sequence_encoding;
if(vc->params)
md = &(vc->params->model_details);
else if(vc->exp_params)
md = &(vc->exp_params->model_details);
else
vrna_message_error("missing model_details in fold_compound");
min_loop_size = md->min_loop_size;
max_span = md->max_bp_span;
if((max_span < 5) || (max_span > n))
max_span = n;
/* ######################### */
/* fill with default values */
/* ######################### */
/* 1. unpaired nucleotides are allowed in all contexts */
for(i = 1; i <= n; i++)
hc->matrix[idx[i] + i] = VRNA_CONSTRAINT_CONTEXT_EXT_LOOP
| VRNA_CONSTRAINT_CONTEXT_HP_LOOP
| VRNA_CONSTRAINT_CONTEXT_INT_LOOP
| VRNA_CONSTRAINT_CONTEXT_MB_LOOP;
/* 2. all base pairs with pscore above threshold are allowed in all contexts */
switch(vc->type){
case VRNA_FC_TYPE_COMPARATIVE: for(j = n; j > min_loop_size + 1; j--){
ij = idx[j]+1;
for(i=1; i < j - min_loop_size; i++, ij++){
char opt = (char)0;
if((j-i+1) <= max_span){
if(vc->pscore[idx[j]+i] >= md->cv_fact*MINPSCORE)
opt = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
}
hc->matrix[ij] = opt;
}
}
break;
case VRNA_FC_TYPE_SINGLE: for(j = n; j > min_loop_size + 1; j--){
ij = idx[j]+1;
for(i=1; i < j - min_loop_size; i++, ij++){
char opt = (char)0;
if((j-i+1) <= max_span){
int t = md->pair[S[i]][S[j]];
switch(t){
case 0: break;
case 3: /* fallthrough */
case 4: if(md->noGU){
break;
} else if(md->noGUclosure){
opt = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
opt &= ~(VRNA_CONSTRAINT_CONTEXT_HP_LOOP | VRNA_CONSTRAINT_CONTEXT_MB_LOOP);
break;
} /* else fallthrough */
default: opt = VRNA_CONSTRAINT_CONTEXT_ALL_LOOPS;
break;
}
}
hc->matrix[ij] = opt;
}
}
/* correct for no lonely pairs (assuming that ptypes already incorporate noLP status) */
/* this should be fixed such that ij loses its hard constraint type if it does not
allow for enclosing an interior loop, etc.
*/
/* ???????
Is this necessary? We could leave the noLP option somewhere else, i.e. do not enforce it
on the level of ptype/constraints, but an the level of recursions...
???????
*/
if(md->noLP){
if(!vc->ptype)
vc->ptype = vrna_ptypes(vc->sequence_encoding2, md);
for(i = 1; i < n; i++)
for(j = i + min_loop_size + 1; j <= n; j++){
if(hc->matrix[idx[j] +i]){
if(!vc->ptype[idx[j] + i]){
hc->matrix[idx[j] + i] = (char)0;
}
}
}
}
break;
default: break;
}
/* should we reset the generalized hard constraint feature here? */
if(hc->f || hc->data){
if(hc->free_data)
hc->free_data(hc->data);
hc->f = NULL;
hc->data = NULL;
hc->free_data = NULL;
}
}
PRIVATE void
hc_update_up(vrna_fold_compound_t *vc){
unsigned int i, n;
int *idx;
vrna_hc_t *hc;
n = vc->length;
idx = vc->jindx;
hc = vc->hc;
for(hc->up_ext[n+1] = 0, i = n; i > 0; i--) /* unpaired stretch in exterior loop */
hc->up_ext[i] = (hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) ? 1 + hc->up_ext[i+1] : 0;
for(hc->up_hp[n+1] = 0, i = n; i > 0; i--) /* unpaired stretch in hairpin loop */
hc->up_hp[i] = (hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP) ? 1 + hc->up_hp[i+1] : 0;
for(hc->up_int[n+1] = 0, i = n; i > 0; i--) /* unpaired stretch in interior loop */
hc->up_int[i] = (hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP) ? 1 + hc->up_int[i+1] : 0;
for(hc->up_ml[n+1] = 0, i = n; i > 0; i--) /* unpaired stretch in multibranch loop */
hc->up_ml[i] = (hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) ? 1 + hc->up_ml[i+1] : 0;
/*
* loop arround once more until we find a nucleotide that mustn't
* be unpaired (needed for circular folding)
*/
if(hc->matrix[idx[1]+1] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
hc->up_ext[n+1] = hc->up_ext[1];
for(i = n; i > 0; i--){
if(hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
hc->up_ext[i] = MIN2(n, 1 + hc->up_ext[i+1]);
} else
break;
}
}
if(hc->matrix[idx[1]+1] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP){
hc->up_hp[n+1] = hc->up_hp[1];
for(i = n; i > 0; i--){
if(hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP){
hc->up_hp[i] = MIN2(n, 1 + hc->up_hp[i+1]);
} else
break;
}
}
if(hc->matrix[idx[1]+1] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
hc->up_int[n+1] = hc->up_int[1];
for(i = n; i > 0; i--){
if(hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
hc->up_int[i] = MIN2(n, 1 + hc->up_int[i+1]);
} else
break;
}
}
if(hc->matrix[idx[1]+1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
hc->up_ml[n+1] = hc->up_ml[1];
for(i = n; i > 0; i--){
if(hc->matrix[idx[i]+i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
hc->up_ml[i] = MIN2(n, 1 + hc->up_ml[i+1]);
} else
break;
}
}
}
#ifdef VRNA_BACKWARD_COMPAT
/*###########################################*/
/*# deprecated functions below #*/
/*###########################################*/
PUBLIC void
print_tty_constraint_full(void){
vrna_message_constraint_options_all();
}
PUBLIC void
print_tty_constraint(unsigned int option){
vrna_message_constraint_options(option);
}
PUBLIC void
constrain_ptypes( const char *constraint,
unsigned int length,
char *ptype,
int *BP,
int min_loop_size,
unsigned int idx_type){
int n,i,j,k,l;
int hx, *stack;
char type;
int *index;
if(constraint == NULL) return;
n = (int)strlen(constraint);
stack = vrna_alloc(sizeof(int)*(n+1));
if(!idx_type){ /* index allows access in energy matrices at pos (i,j) via index[j]+i */
index = vrna_idx_col_wise(length);
for(hx=0, j=1; j<=n; j++){
switch(constraint[j-1]){
case '|': if(BP) BP[j] = -1;
break;
case 'x': /* can't pair */
for (l=1; l<j-min_loop_size; l++)
ptype[index[j]+l] = 0;
for (l=j+min_loop_size+1; l<=(int)length; l++)
ptype[index[l]+j] = 0;
break;
case '(': stack[hx++]=j;
/* fallthrough */
case '<': /* pairs upstream */
for (l=1; l<j-min_loop_size; l++)
ptype[index[j]+l] = 0;
break;
case ')': if (hx<=0) {
vrna_message_error("%s\nunbalanced brackets in constraint", constraint);
}
i = stack[--hx];
type = ptype[index[j]+i];
for (k=i+1; k<=(int)length; k++)
ptype[index[k]+i] = 0;
/* don't allow pairs i<k<j<l */
for (l=j; l<=(int)length; l++)
for (k=i+1; k<=j; k++)
ptype[index[l]+k] = 0;
/* don't allow pairs k<i<l<j */
for (l=i; l<=j; l++)
for (k=1; k<=i; k++)
ptype[index[l]+k] = 0;
for (k=1; k<j; k++)
ptype[index[j]+k] = 0;
ptype[index[j]+i] = (type==0) ? 7 : type;
/* fallthrough */
case '>': /* pairs downstream */
for (l=j+min_loop_size+1; l<=(int)length; l++)
ptype[index[l]+j] = 0;
break;
}
}
}
else{ /* index allows access in energy matrices at pos (i,j) via index[i]-j */
index = vrna_idx_row_wise(length);
for(hx=0, j=1; j<=n; j++) {
switch (constraint[j-1]) {
case 'x': /* can't pair */
for (l=1; l<j-min_loop_size; l++)
ptype[index[l]-j] = 0;
for (l=j+min_loop_size+1; l<=(int)length; l++)
ptype[index[j]-l] = 0;
break;
case '(': stack[hx++]=j;
/* fallthrough */
case '<': /* pairs upstream */
for (l=1; l<j-min_loop_size; l++)
ptype[index[l]-j] = 0;
break;
case ')': if (hx<=0) {
vrna_message_error("%s\nunbalanced brackets in constraints", constraint);
}
i = stack[--hx];
type = ptype[index[i]-j];
/* don't allow pairs i<k<j<l */
for (k=i; k<=j; k++)
for (l=j; l<=(int)length; l++)
ptype[index[k]-l] = 0;
/* don't allow pairs k<i<l<j */
for (k=1; k<=i; k++)
for (l=i; l<=j; l++)
ptype[index[k]-l] = 0;
ptype[index[i]-j] = (type==0) ? 7 : type;
/* fallthrough */
case '>': /* pairs downstream */
for (l=j+min_loop_size+1; l<=(int)length; l++)
ptype[index[j]-l] = 0;
break;
}
}
}
if (hx!=0) {
vrna_message_error("%s\nunbalanced brackets in constraint string", constraint);
}
free(index);
free(stack);
}
#endif