/*
suboptimal folding - Stefan Wuchty, Walter Fontana & Ivo Hofacker
Vienna RNA package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include "ViennaRNA/fold.h"
#include "ViennaRNA/constraints.h"
#include "ViennaRNA/utils.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/fold_vars.h"
#include "list.h"
#include "ViennaRNA/eval.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/loop_energies.h"
#include "ViennaRNA/cofold.h"
#include "ViennaRNA/gquad.h"
#include "ViennaRNA/subopt.h"
/* hack */
#include "ViennaRNA/color_output.inc"
#ifdef _OPENMP
#include <omp.h>
#endif
#define true 1
#define false 0
#define ON_SAME_STRAND(I,J,C) (((I)>=(C))||((J)<(C)))
/**
* @brief Sequence interval stack element used in subopt.c
*/
typedef struct INTERVAL {
int i;
int j;
int array_flag;
} INTERVAL;
typedef struct {
char *structure;
LIST *Intervals;
int partial_energy;
int is_duplex;
/* int best_energy; */ /* best attainable energy */
} STATE;
typedef struct {
LIST *Intervals;
LIST *Stack;
int nopush;
} subopt_env;
struct old_subopt_dat {
unsigned long max_sol;
unsigned long n_sol;
SOLUTION *SolutionList;
FILE *fp;
};
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
PUBLIC int subopt_sorted=0; /* output sorted by energy */
PUBLIC int density_of_states[MAXDOS+1];
PUBLIC double print_energy = 9999; /* printing threshold for use with logML */
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
/* some backward compatibility stuff */
PRIVATE int backward_compat = 0;
PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
#ifdef _OPENMP
#pragma omp threadprivate(backward_compat_compound, backward_compat)
#endif
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
#ifdef VRNA_BACKWARD_COMPAT
PRIVATE SOLUTION *
wrap_subopt(char *seq,
char *structure,
vrna_param_t *parameters,
int delta,
int is_constrained,
int is_circular,
FILE *fp);
#endif
PRIVATE void make_pair(int i, int j, STATE *state);
/* mark a gquadruplex in the resulting dot-bracket structure */
PRIVATE void make_gquad(int i, int L, int l[3], STATE *state);
PRIVATE INTERVAL *make_interval (int i, int j, int ml);
/*@out@*/ PRIVATE STATE *make_state(/*@only@*/LIST *Intervals,
/*@only@*/ /*@null@*/ char *structure,
int partial_energy, int is_duplex, int length);
PRIVATE STATE *copy_state(STATE * state);
PRIVATE void print_state(STATE * state);
//PRIVATE void UNUSED print_stack(LIST * list);
/*@only@*/ PRIVATE LIST *make_list(void);
PRIVATE void push(LIST * list, /*@only@*/ void *data);
PRIVATE void *pop(LIST * list);
PRIVATE int best_attainable_energy(vrna_fold_compound_t *vc, STATE * state);
PRIVATE void scan_interval(vrna_fold_compound_t *vc, int i, int j, int array_flag, int threshold, STATE * state, subopt_env *env);
PRIVATE void free_interval_node(/*@only@*/ INTERVAL * node);
PRIVATE void free_state_node(/*@only@*/ STATE * node);
PRIVATE void push_back(LIST *Stack, STATE * state);
PRIVATE char* get_structure(STATE * state);
PRIVATE int compare(const void *solution1, const void *solution2);
PRIVATE void make_output(SOLUTION *SL, int cp, FILE *fp);
PRIVATE void repeat(vrna_fold_compound_t *vc, int i, int j, STATE * state, int part_energy, int temp_energy, int best_energy, int threshold, subopt_env *env);
PRIVATE void repeat_gquad(vrna_fold_compound_t *vc, int i, int j, STATE *state, int part_energy, int temp_energy, int best_energy, int threshold, subopt_env *env);
PRIVATE void old_subopt_print( const char *structure, float energy, void *data);
PRIVATE void old_subopt_store( const char *structure, float energy, void *data);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
/*---------------------------------------------------------------------------*/
/*List routines--------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
PRIVATE void
make_pair(int i, int j, STATE *state)
{
state->structure[i-1] = '(';
state->structure[j-1] = ')';
}
PRIVATE void
make_gquad(int i, int L, int l[3], STATE *state)
{
int x;
for(x = 0; x < L; x++){
state->structure[i - 1 + x] = '+';
state->structure[i - 1 + x + L + l[0]] = '+';
state->structure[i - 1 + x + 2*L + l[0] + l[1]] = '+';
state->structure[i - 1 + x + 3*L + l[0] + l[1] + l[2]] = '+';
}
}
/*---------------------------------------------------------------------------*/
PRIVATE INTERVAL *
make_interval(int i, int j, int array_flag)
{
INTERVAL *interval;
interval = lst_newnode(sizeof(INTERVAL));
interval->i = i;
interval->j = j;
interval->array_flag = array_flag;
return interval;
}
/*---------------------------------------------------------------------------*/
PRIVATE void
free_interval_node(INTERVAL * node)
{
lst_freenode(node);
}
/*---------------------------------------------------------------------------*/
PRIVATE void
free_state_node(STATE * node)
{
free(node->structure);
if (node->Intervals)
lst_kill(node->Intervals, lst_freenode);
lst_freenode(node);
}
/*---------------------------------------------------------------------------*/
PRIVATE STATE *
make_state(LIST * Intervals,
char *structure,
int partial_energy,
int is_duplex,
int length)
{
STATE *state;
state = lst_newnode(sizeof(STATE));
if (Intervals)
state->Intervals = Intervals;
else
state->Intervals = lst_init();
if (structure)
state->structure = structure;
else {
int i;
state->structure = (char *) vrna_alloc(length+1);
for (i=0; i<length; i++)
state->structure[i] = '.';
}
state->partial_energy = partial_energy;
return state;
}
/*---------------------------------------------------------------------------*/
PRIVATE STATE *
copy_state(STATE * state)
{
STATE *new_state;
void *after;
INTERVAL *new_interval, *next;
new_state = lst_newnode(sizeof(STATE));
new_state->Intervals = lst_init();
new_state->partial_energy = state->partial_energy;
/* new_state->best_energy = state->best_energy; */
if (state->Intervals->count) {
after = LST_HEAD(new_state->Intervals);
for ( next = lst_first(state->Intervals); next; next = lst_next(next))
{
new_interval = lst_newnode(sizeof(INTERVAL));
*new_interval = *next;
lst_insertafter(new_state->Intervals, new_interval, after);
after = new_interval;
}
}
new_state->structure = strdup(state->structure);
if (!new_state->structure) vrna_message_error("out of memory");
return new_state;
}
/*---------------------------------------------------------------------------*/
/*@unused @*/ PRIVATE void
print_state(STATE * state)
{
INTERVAL *next;
if (state->Intervals->count)
{
printf("%d intervals:\n", state->Intervals->count);
for (next = lst_first(state->Intervals); next; next = lst_next(next))
{
printf("[%d,%d],%d ", next->i, next->j, next->array_flag);
}
printf("\n");
}
printf("partial structure: %s\n", state->structure);
printf("\n");
printf(" partial_energy: %d\n", state->partial_energy);
/* printf(" best_energy: %d\n", state->best_energy); */
(void) fflush(stdout);
}
/*---------------------------------------------------------------------------*/
/*@unused @*/ PRIVATE void
print_stack(LIST * list)
{
void *rec;
printf("================\n");
printf("%d states\n", list->count);
for (rec = lst_first(list); rec; rec = lst_next(rec))
{
printf("state-----------\n");
print_state(rec);
}
printf("================\n");
}
/*---------------------------------------------------------------------------*/
PRIVATE LIST *
make_list(void)
{
return lst_init();
}
/*---------------------------------------------------------------------------*/
PRIVATE void
push(LIST * list, void *data)
{
lst_insertafter(list, data, LST_HEAD(list));
}
/* PRIVATE void */
/* push_stack(STATE *state) { */ /* keep the stack sorted by energy */
/* STATE *after, *next; */
/* nopush = false; */
/* next = after = LST_HEAD(Stack); */
/* while ( next = lst_next(next)) { */
/* if ( next->best_energy >= state->best_energy ) break; */
/* after = next; */
/* } */
/* lst_insertafter(Stack, state, after); */
/* } */
/*---------------------------------------------------------------------------*/
PRIVATE void *
pop(LIST * list)
{
void *data;
data = lst_deletenext(list, LST_HEAD(list));
return data;
}
/*---------------------------------------------------------------------------*/
/*auxiliary routines---------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
PRIVATE int
best_attainable_energy( vrna_fold_compound_t *vc,
STATE *state){
/* evaluation of best possible energy attainable within remaining intervals */
register int sum;
INTERVAL *next;
vrna_md_t *md;
vrna_mx_mfe_t *matrices;
int *indx;
md = &(vc->params->model_details);
matrices = vc->matrices;
indx = vc->jindx;
sum = state->partial_energy; /* energy of already found elements */
for (next = lst_first(state->Intervals); next; next = lst_next(next))
{
if (next->array_flag == 0)
sum += (md->circ) ? matrices->Fc : matrices->f5[next->j];
else if (next->array_flag == 1)
sum += matrices->fML[indx[next->j] + next->i];
else if (next->array_flag == 2)
sum += matrices->c[indx[next->j] + next->i];
else if (next->array_flag == 3)
sum += matrices->fM1[indx[next->j] + next->i];
else if (next->array_flag == 4)
sum += matrices->fc[next->i];
else if (next->array_flag == 5)
sum += matrices->fc[next->j];
else if (next->array_flag == 6)
sum += matrices->ggg[indx[next->j] + next->i];
}
return sum;
}
/*---------------------------------------------------------------------------*/
PRIVATE void
push_back(LIST *Stack, STATE * state)
{
push(Stack, copy_state(state));
return;
}
/*---------------------------------------------------------------------------*/
PRIVATE char*
get_structure(STATE * state)
{
char* structure;
structure = strdup(state->structure);
return structure;
}
/*---------------------------------------------------------------------------*/
PRIVATE int
compare(const void *solution1, const void *solution2)
{
if (((SOLUTION *) solution1)->energy > ((SOLUTION *) solution2)->energy)
return 1;
if (((SOLUTION *) solution1)->energy < ((SOLUTION *) solution2)->energy)
return -1;
return strcmp(((SOLUTION *) solution1)->structure,
((SOLUTION *) solution2)->structure);
}
/*---------------------------------------------------------------------------*/
PRIVATE void make_output(SOLUTION *SL, int cp, FILE *fp) /* prints stuff */
{
SOLUTION *sol;
for (sol = SL; sol->structure!=NULL; sol++){
char *e_string = vrna_strdup_printf(" %6.2f", sol->energy);
print_structure(fp, sol->structure, e_string);
free(e_string);
}
}
PRIVATE STATE *
derive_new_state( int i,
int j,
STATE *s,
int e,
int flag){
STATE *s_new = copy_state(s);
INTERVAL *ival = make_interval(i, j, flag);
push(s_new->Intervals, ival);
s_new->partial_energy += e;
return s_new;
}
PRIVATE void
fork_state( int i,
int j,
STATE *s,
int e,
int flag,
subopt_env *env){
STATE *s_new = derive_new_state(i, j, s, e, flag);
push(env->Stack, s_new);
env->nopush = false;
}
PRIVATE void
fork_int_state( int i, int j,
int p, int q,
STATE *s,
int e,
subopt_env *env){
STATE *s_new = derive_new_state(p, q, s, e, 2);
make_pair(i, j, s_new);
make_pair(p, q, s_new);
push(env->Stack, s_new);
env->nopush = false;
}
PRIVATE void
fork_state_pair(int i,
int j,
STATE *s,
int e,
subopt_env *env){
STATE *new_state;
new_state = copy_state(s);
make_pair(i, j, new_state);
new_state->partial_energy += e;
push(env->Stack, new_state);
env->nopush = false;
}
PRIVATE void
fork_two_states_pair( int i,
int j,
int k,
STATE *s,
int e,
int flag1,
int flag2,
subopt_env *env){
INTERVAL *interval1, *interval2;
STATE *new_state;
new_state = copy_state(s);
interval1 = make_interval(i+1, k-1, flag1);
interval2 = make_interval(k, j-1, flag2);
if (k-i < j-k) { /* push larger interval first */
push(new_state->Intervals, interval1);
push(new_state->Intervals, interval2);
} else {
push(new_state->Intervals, interval2);
push(new_state->Intervals, interval1);
}
make_pair(i, j, new_state);
new_state->partial_energy += e;
push(env->Stack, new_state);
env->nopush = false;
}
PRIVATE void
fork_two_states(int i,
int j,
int p,
int q,
STATE *s,
int e,
int flag1,
int flag2,
subopt_env *env){
INTERVAL *interval1, *interval2;
STATE *new_state;
new_state = copy_state(s);
interval1 = make_interval(i, j, flag1);
interval2 = make_interval(p, q, flag2);
if((j - i) < (q - p)){
push(new_state->Intervals, interval1);
push(new_state->Intervals, interval2);
} else {
push(new_state->Intervals, interval2);
push(new_state->Intervals, interval1);
}
new_state->partial_energy += e;
push(env->Stack, new_state);
env->nopush = false;
}
/*---------------------------------------------------------------------------*/
/* start of subopt backtracking ---------------------------------------------*/
/*---------------------------------------------------------------------------*/
PUBLIC SOLUTION *
vrna_subopt(vrna_fold_compound_t *vc,
int delta,
int sorted,
FILE *fp){
SOLUTION *SolutionList;
unsigned long max_sol, n_sol;
struct old_subopt_dat data;
data.SolutionList = NULL;
data.max_sol = 128;
data.n_sol = 0;
data.fp = fp;
if(vc){
/* SolutionList stores the suboptimal structures found */
data.SolutionList = (SOLUTION *) vrna_alloc(data.max_sol*sizeof(SOLUTION));
/* end initialize ------------------------------------------------------- */
if (fp) {
float min_en;
char *SeQ, *energies = NULL;
if(vc->cutpoint > 0)
min_en = vrna_mfe_dimer(vc, NULL);
else
min_en = vrna_mfe(vc, NULL);
SeQ = vrna_cut_point_insert(vc->sequence, vc->cutpoint);
energies = vrna_strdup_printf(" %6.2f %6.2f", min_en, (float)delta/100.);
print_structure(fp, SeQ, energies);
free(SeQ);
free(energies);
vrna_mx_mfe_free(vc);
}
/* call subopt() */
vrna_subopt_cb(vc, delta, (!sorted && fp) ? &old_subopt_print : &old_subopt_store, (void *)&data);
if(sorted){
/* sort structures by energy */
if(data.n_sol > 0)
qsort(data.SolutionList, data.n_sol - 1, sizeof(SOLUTION), compare);
if(fp)
make_output(data.SolutionList, vc->cutpoint, fp);
}
if(fp){ /* we've printed everything -- free solutions */
SOLUTION *sol;
for(sol = data.SolutionList; sol->structure != NULL; sol++)
free(sol->structure);
free(data.SolutionList);
data.SolutionList = NULL;
}
}
return data.SolutionList;
}
PUBLIC void
vrna_subopt_cb( vrna_fold_compound_t *vc,
int delta,
vrna_subopt_callback *cb,
void *data){
subopt_env *env;
STATE *state;
INTERVAL *interval;
int maxlevel, count, partial_energy, old_dangles, logML, dangle_model, length, circular, threshold, cp;
double structure_energy, min_en, eprint;
char *struc, *structure, *sequence;
float correction;
vrna_param_t *P;
vrna_md_t *md;
int minimal_energy;
int Fc;
int *f5;
vrna_fold_compound_prepare(vc, VRNA_OPTION_MFE | VRNA_OPTION_HYBRID);
sequence = vc->sequence;
length = vc->length;
cp = vc->cutpoint;
P = vc->params;
md = &(P->model_details);
/* do mfe folding to get fill arrays and get ground state energy */
/* in case dangles is neither 0 or 2, set dangles=2 while folding */
circular = md->circ;
logML = md->logML;
old_dangles = dangle_model = md->dangles;
if(md->uniq_ML != 1) /* failsafe mechanism to enforce valid fM1 array */
md->uniq_ML = 1;
/* temporarily set dangles to 2 if necessary */
if((md->dangles != 0) && (md->dangles != 2))
md->dangles = 2;
struc = (char *)vrna_alloc(sizeof(char) * (length + 1));
if(circular){
min_en = vrna_mfe(vc, struc);
Fc = vc->matrices->Fc;
f5 = vc->matrices->f5;
/* restore dangle model */
md->dangles = old_dangles;
/* re-evaluate in case we're using logML etc */
min_en = vrna_eval_structure(vc, struc);
} else {
min_en = vrna_mfe_dimer(vc, struc);
f5 = vc->matrices->f5;
/* restore dangle model */
md->dangles = old_dangles;
/* re-evaluate in case we're using logML etc */
min_en = vrna_eval_structure(vc, struc);
}
free(struc);
eprint = print_energy + min_en;
correction = (min_en < 0) ? -0.1 : 0.1;
/* Initialize ------------------------------------------------------------ */
maxlevel = 0;
count = 0;
partial_energy = 0;
/* Initialize the stack ------------------------------------------------- */
minimal_energy = (circular) ? Fc : f5[length];
threshold = minimal_energy + delta;
if(threshold > INF){
vrna_message_warning("Energy range too high, limiting to reasonable value");
threshold = INF-EMAX;
}
/* init env data structure */
env = (subopt_env *)vrna_alloc(sizeof(subopt_env));
env->Stack = NULL;
env->nopush = true;
env->Stack = make_list(); /* anchor */
env->Intervals = make_list(); /* initial state: */
interval = make_interval(1, length, 0); /* interval [1,length,0] */
push(env->Intervals, interval);
env->nopush = false;
state = make_state(env->Intervals, NULL, partial_energy,0, length);
/* state->best_energy = minimal_energy; */
push(env->Stack, state);
env->nopush = false;
/* end initialize ------------------------------------------------------- */
while (1) { /* forever, til nothing remains on stack */
maxlevel = (env->Stack->count > maxlevel ? env->Stack->count : maxlevel);
if (LST_EMPTY (env->Stack)) /* we are done! clean up and quit */
{
/* fprintf(stderr, "maxlevel: %d\n", maxlevel); */
lst_kill(env->Stack, free_state_node);
cb(NULL, 0, data); /* NULL (last time to call callback function */
break;
}
/* pop the last element ---------------------------------------------- */
state = pop(env->Stack); /* current state to work with */
if (LST_EMPTY(state->Intervals))
{
int e;
/* state has no intervals left: we got a solution */
count++;
structure = get_structure(state);
structure_energy = state->partial_energy / 100.;
#ifdef CHECK_ENERGY
structure_energy = vrna_eval_structure(vc, structure);
if (!logML)
if ((double) (state->partial_energy / 100.) != structure_energy) {
vrna_message_error("%s %6.2f %6.2f",
structure,
state->partial_energy / 100.,
structure_energy );
exit(1);
}
#endif
if (logML || (dangle_model==1) || (dangle_model==3)) { /* recalc energy */
structure_energy = vrna_eval_structure(vc, structure);
}
e = (int) ((structure_energy-min_en)*10. - correction); /* avoid rounding errors */
if (e>MAXDOS) e=MAXDOS;
density_of_states[e]++;
if(structure_energy <= eprint){
char *outstruct = vrna_cut_point_insert(structure, cp);
cb((const char *)outstruct, structure_energy, data);
free(outstruct);
}
free(structure);
}
else {
/* get (and remove) next interval of state to analyze */
interval = pop(state->Intervals);
scan_interval(vc, interval->i, interval->j, interval->array_flag, threshold, state, env);
free_interval_node(interval); /* free the current interval */
}
free_state_node(state); /* free the current state */
} /* end of while (1) */
/* cleanup memory */
free(env);
}
PRIVATE void
scan_interval(vrna_fold_compound_t *vc,
int i,
int j,
int array_flag,
int threshold,
STATE * state,
subopt_env *env){
/* real backtrack routine */
/* array_flag = 0: trace back in f5-array */
/* array_flag = 1: trace back in fML-array */
/* array_flag = 2: trace back in repeat() */
/* array_flag = 3: trace back in fM1-array */
STATE *new_state, *temp_state;
INTERVAL *new_interval;
vrna_param_t *P;
vrna_md_t *md;
register int k, fi, cij, ij;
register int type;
register int dangle_model;
register int noLP;
int element_energy, best_energy;
int *fc, *f5, *c, *fML, *fM1, *ggg;
int FcH, FcI, FcM, *fM2;
int length, *indx, *rtype, circular, with_gquad, turn, cp;
char *ptype;
short *S1;
char *hard_constraints, hc_decompose;
vrna_hc_t *hc;
vrna_sc_t *sc;
length = vc->length;
cp = vc->cutpoint;
indx = vc->jindx;
ptype = vc->ptype;
S1 = vc->sequence_encoding;
P = vc->params;
md = &(P->model_details);
rtype = &(md->rtype[0]);
dangle_model = md->dangles;
noLP = md->noLP;
circular = md->circ;
with_gquad = md->gquad;
turn = md->min_loop_size;
fc = vc->matrices->fc;
f5 = vc->matrices->f5;
c = vc->matrices->c;
fML = vc->matrices->fML;
fM1 = vc->matrices->fM1;
ggg = vc->matrices->ggg;
FcH = vc->matrices->FcH;
FcI = vc->matrices->FcI;
FcM = vc->matrices->FcM;
fM2 = vc->matrices->fM2;
hc = vc->hc;
hard_constraints = hc->matrix;
sc = vc->sc;
best_energy = best_attainable_energy(vc, state); /* .. on remaining intervals */
env->nopush = true;
if ((i > 1) && (!array_flag))
vrna_message_error ("Error while backtracking!");
if (j < i + turn + 1 && ON_SAME_STRAND(i,j,cp)) { /* minimal structure element */
if(array_flag == 0){
/* do not forget to add f5[j], since it may contain pseudo energies from soft constraining */
state->partial_energy += f5[j];
}
if (env->nopush){
push_back(env->Stack, state);
env->nopush = false;
}
return;
}
ij = indx[j] + i;
/* 13131313131313131313131313131313131313131313131313131313131313131313131 */
if (array_flag == 3 || array_flag == 1) {
/* array_flag = 3: interval i,j was generated during */
/* a multiloop decomposition using array fM1 in repeat() */
/* or in this block */
/* array_flag = 1: interval i,j was generated from a */
/* stack, bulge, or internal loop in repeat() */
/* or in this block */
if(hc->up_ml[j]){
if (array_flag == 3)
fi = fM1[indx[j-1] + i] + P->MLbase;
else
fi = fML[indx[j-1] + i] + P->MLbase;
if(sc){
if(sc->energy_up)
fi += sc->energy_up[j][1];
if(sc->f)
fi += sc->f(i, j, i, j - 1, VRNA_DECOMP_ML_ML, sc->data);
}
if ((fi + best_energy <= threshold)&&(ON_SAME_STRAND(j-1,j, cp))) {
/* no basepair, nibbling of 3'-end */
fork_state(i, j-1, state, P->MLbase, array_flag, env);
}
}
hc_decompose = hard_constraints[ij];
if (hc_decompose & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC) { /* i,j may pair */
cij = c[ij];
type = ptype[ij];
if(type == 0)
type = 7;
switch(dangle_model){
case 0: element_energy = E_MLstem(type, -1, -1, P);
break;
default: element_energy = E_MLstem(type,
(((i > 1)&&(ON_SAME_STRAND(i-1,i,cp))) || circular) ? S1[i-1] : -1,
(((j < length)&&(ON_SAME_STRAND(j,j+1,cp))) || circular) ? S1[j+1] : -1,
P);
break;
}
if(sc){
/* should be unnecessary
if(sc->energy_bp)
element_energy += sc->energy_bp[ij];
*/
if(sc->f)
element_energy += sc->f(i, j, i, j, VRNA_DECOMP_ML_STEM, sc->data);
}
cij += element_energy;
if (cij + best_energy <= threshold)
repeat(vc, i, j, state, element_energy, 0, best_energy, threshold, env);
} else if (with_gquad){
element_energy = E_MLstem(0, -1, -1, P);
cij = ggg[ij] + element_energy;
if(cij + best_energy <= threshold)
repeat_gquad(vc, i, j, state, element_energy, 0, best_energy, threshold, env);
}
} /* array_flag == 3 || array_flag == 1 */
/* 11111111111111111111111111111111111111111111111111111111111111111111111 */
if (array_flag == 1) {
/* array_flag = 1: interval i,j was generated from a */
/* stack, bulge, or internal loop in repeat() */
/* or in this block */
int stopp, k1j;
if ((ON_SAME_STRAND(i-1,i,cp))&&(ON_SAME_STRAND(j,j+1,cp))) { /*backtrack in FML only if multiloop is possible*/
for ( k = i+turn+1 ; k <= j-1-turn ; k++) {
/* Multiloop decomposition if i,j contains more than 1 stack */
if(with_gquad){
if(ON_SAME_STRAND(k, k+1, cp)){
element_energy = E_MLstem(0, -1, -1, P);
if(fML[indx[k]+i] + ggg[indx[j] + k + 1] + element_energy + best_energy <= threshold){
temp_state = derive_new_state(i, k, state, 0, array_flag);
env->nopush = false;
repeat_gquad(vc, k+1, j, temp_state, element_energy, fML[indx[k]+i], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
k1j = indx[j] + k + 1;
if(hard_constraints[k1j] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC){
short s5, s3;
type = ptype[k1j];
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s5 = s3 = -1;
break;
default: s5 = (ON_SAME_STRAND(i-1,i,cp)) ? S1[k] : -1;
s3 = (ON_SAME_STRAND(j,j+1,cp)) ? S1[j+1] : -1;
break;
}
element_energy = E_MLstem(type, s5, s3, P);
if(sc){
/* should be unnecessary
if(sc->energy_bp)
element_energy += sc->energy_bp[k1j];
*/
if(sc->f)
element_energy += sc->f(i, j, k, k + 1, VRNA_DECOMP_ML_ML_STEM, sc->data);
}
if(ON_SAME_STRAND(k, k+1, cp)){
if(fML[indx[k]+i] + c[k1j] + element_energy + best_energy <= threshold){
temp_state = derive_new_state(i, k, state, 0, array_flag);
env->nopush = false;
repeat(vc, k+1, j, temp_state, element_energy, fML[indx[k]+i], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
}
}
stopp=(cp>0)? (cp-2):(length); /*if cp -1: k on cut, => no ml*/
stopp=MIN2(stopp, j-1-turn);
if (i>cp) stopp=j-1-turn;
else if (i==cp) stopp=0; /*not a multi loop*/
int up = 1;
for(k = i; k <= stopp; k++, up++){
if(hc->up_ml[i] >= up){
k1j = indx[j] + k + 1;
/* Multiloop decomposition if i,j contains only 1 stack */
if(with_gquad){
element_energy = E_MLstem(0, -1, -1, P) + P->MLbase * up;
if(sc){
if(sc->energy_up)
element_energy += sc->energy_up[i][up];
}
if(ggg[k1j] + element_energy + best_energy <= threshold)
repeat_gquad(vc, k+1, j, state, element_energy, 0, best_energy, threshold, env);
}
if(hard_constraints[k1j] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC){
int s5, s3;
type = ptype[k1j];
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s5 = s3 = -1;
break;
default: s5 = (ON_SAME_STRAND(k-1,k,cp)) ? S1[k] : -1;
s3 = (ON_SAME_STRAND(j,j+1,cp)) ? S1[j+1] : -1;
break;
}
element_energy = E_MLstem(type, s5, s3, P);
element_energy += P->MLbase * up;
if(sc){
if(sc->energy_up)
element_energy += sc->energy_up[i][up];
/* should be unnecessary
if(sc->energy_bp)
element_energy += sc->energy_bp[k1j];
*/
}
if (c[k1j] + element_energy + best_energy <= threshold)
repeat(vc, k+1, j, state, element_energy, 0, best_energy, threshold, env);
}
}
}
} /* array_flag == 1 */
/* 22222222222222222222222222222222222222222222222222 */
/* */
/* array_flag = 2: interval i,j was generated from a */
/* stack, bulge, or internal loop in repeat() */
/* */
/* 22222222222222222222222222222222222222222222222222 */
if(array_flag == 2){
repeat(vc, i, j, state, 0, 0, best_energy, threshold, env);
if (env->nopush){
if (!noLP){
vrna_message_warning("%d,%d\nOops, no solution in repeat!", i, j);
}
}
return;
}
/* 00000000000000000000000000000000000000000000000000 */
/* */
/* array_flag = 0: interval i,j was found while */
/* tracing back through f5-array and c-array */
/* or within this block */
/* */
/* 00000000000000000000000000000000000000000000000000 */
if((array_flag == 0) && !circular){
int s5, s3, kj, tmp_en;
if(hc->up_ext[j]){
tmp_en = 0;
if(sc){
if(sc->energy_up)
tmp_en += sc->energy_up[j][1];
}
if (f5[j-1] + tmp_en + best_energy <= threshold) {
/* no basepair, nibbling of 3'-end */
fork_state(i, j-1, state, tmp_en, 0, env);
}
}
for (k = j-turn-1; k > 1; k--) {
kj = indx[j] + k;
if(with_gquad){
if(ON_SAME_STRAND(k,j,cp)){
element_energy = 0;
if(f5[k-1] + ggg[kj] + element_energy + best_energy <= threshold){
temp_state = derive_new_state(1, k-1, state, 0, 0);
env->nopush = false;
/* backtrace the quadruplex */
repeat_gquad(vc, k, j, temp_state, element_energy, f5[k-1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
if(hard_constraints[kj] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[kj];
if(type == 0)
type = 7;
/* k and j pair */
switch(dangle_model){
case 0: s5 = s3 = -1;
break;
default: s5 = (ON_SAME_STRAND(k-1,k,cp)) ? S1[k-1] : -1;
s3 = ((j < length)&&(ON_SAME_STRAND(j,j+1,cp))) ? S1[j+1] : -1;
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
if (!(ON_SAME_STRAND(k,j,cp)))/*&&(state->is_duplex==0))*/ {
element_energy+=P->DuplexInit;
/*state->is_duplex=1;*/
}
if (f5[k-1] + c[kj] + element_energy + best_energy <= threshold){
temp_state = derive_new_state(1, k-1, state, 0, 0);
env->nopush = false;
repeat(vc, k, j, temp_state, element_energy, f5[k-1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
kj = indx[j] + 1;
if(with_gquad){
if(ON_SAME_STRAND(k,j,cp)){
element_energy = 0;
if(ggg[kj] + element_energy + best_energy <= threshold){
/* backtrace the quadruplex */
repeat_gquad(vc, 1, j, state, element_energy, 0, best_energy, threshold, env);
}
}
}
if(hard_constraints[kj] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[kj];
s5 = -1;
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s3 = -1;
break;
default: s3 = (j < length) && (ON_SAME_STRAND(j,j+1,cp)) ? S1[j+1] : -1;
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
if(!(ON_SAME_STRAND(1,j,cp)))
element_energy += P->DuplexInit;
if (c[kj] + element_energy + best_energy <= threshold)
repeat(vc, 1, j, state, element_energy, 0, best_energy, threshold, env);
}
} /* end array_flag == 0 && !circular*/
/* or do we subopt circular? */
else if(array_flag == 0){
int k, l, p, q, tmp_en;
/* if we've done everything right, we will never reach this case more than once */
/* right after the initilization of the stack with ([1,n], empty, 0) */
/* lets check, if we can have an open chain without breaking the threshold */
/* this is an ugly work-arround cause in case of an open chain we do not have to */
/* backtrack anything further... */
if(hc->up_ext[1] >= length){
tmp_en = 0;
if(sc){
if(sc->energy_up)
tmp_en += sc->energy_up[1][length];
}
if(tmp_en <= threshold){
new_state = derive_new_state(1,2,state,0,0);
new_state->partial_energy = 0;
push(env->Stack, new_state);
env->nopush = false;
}
}
/* ok, lets check if we can do an exterior hairpin without breaking the threshold */
/* best energy should be 0 if we are here */
if(FcH + best_energy <= threshold){
/* lets search for all exterior hairpin cases, that fit into our threshold barrier */
/* we use index k,l to avoid confusion with i,j index of our state... */
/* if we reach here, i should be 1 and j should be n respectively */
for(k=i; k<j; k++){
if(hc->up_hp[1] < k)
break;
for (l=j; l >= k + turn + 1; l--){
int kl, tmpE;
kl = indx[l] + k;
tmpE = vrna_E_hp_loop(vc, l, k);
if(c[kl] + tmpE + best_energy <= threshold){
/* what we really have to do is something like this, isn't it? */
/* we have to create a new state, with interval [k,l], then we */
/* add our loop energy as initial energy of this state and put */
/* the state onto the stack R... for further refinement... */
/* we also denote this new interval to be scanned in C */
fork_state(k, l, state, tmpE, 2, env);
}
}
}
}
/* now lets see, if we can do an exterior interior loop without breaking the threshold */
if(FcI + best_energy <= threshold){
/* now we search for our exterior interior loop possibilities */
for(k=i; k<j; k++){
for (l=j; l >= k + turn + 1; l--){
int kl, type, tmpE;
kl = indx[l]+k; /* just confusing these indices ;-) */
if(hard_constraints[kl] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
type = ptype[kl];
type = rtype[type];
if(type == 0)
type = 7;
for (p = l+1; p < j ; p++){
int u1, qmin;
u1 = p-l-1;
if (u1+k-1>MAXLOOP) break;
if (hc->up_int[l+1] < u1) break;
qmin = u1+k-1+j-MAXLOOP;
if(qmin<p+turn+1) qmin = p+turn+1;
for(q = j; q >= qmin; q--){
int u2, type_2;
if(hc->up_int[q+1] < (j - q + k - 1))
break;
if(hard_constraints[indx[q]+p] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
type_2 = rtype[ptype[indx[q]+p]];
if(type_2 == 0)
type_2 = 7;
u2 = k-1 + j-q;
if(u1+u2>MAXLOOP) continue;
tmpE = E_IntLoop(u1, u2, type, type_2, S1[l+1], S1[k-1], S1[p-1], S1[q+1], P);
if(sc){
if(sc->energy_up)
tmpE += sc->energy_up[l+1][p-l-1]
+ sc->energy_up[q+1][j-q]
+ sc->energy_up[1][k-1];
if(sc->energy_stack)
if(u1 + u2 == 0)
tmpE += sc->energy_stack[k]
+ sc->energy_stack[l]
+ sc->energy_stack[p]
+ sc->energy_stack[q];
}
if(c[kl] + c[indx[q]+p] + tmpE + best_energy <= threshold){
/* ok, similar to the hairpin stuff, we add new states onto the stack R */
/* but in contrast to the hairpin decomposition, we have to add two new */
/* intervals, enclosed by k,l and p,q respectively and we also have to */
/* add the partial energy, that comes from the exterior interior loop */
fork_two_states(k, l, p, q, state, tmpE, 2, 2, env);
}
}
}
}
}
}
}
}
/* and last but not least, we have a look, if we can do an exterior multiloop within the energy threshold */
if(FcM <= threshold){
/* this decomposition will be somehow more complicated...so lets see what we do here... */
/* first we want to find out which split inidices we can use without exceeding the threshold */
int tmpE2;
for (k=turn+1; k<j-2*turn; k++){
tmpE2 = fML[indx[k]+1]+fM2[k+1]+P->MLclosing;
if(tmpE2 + best_energy <= threshold){
/* grmpfh, we have found a possible split index k so we have to split fM2 and fML now */
/* lets do it first in fM2 anyway */
for(l=k+turn+2; l<j-turn-1; l++){
tmpE2 = fM1[indx[l]+k+1] + fM1[indx[j]+l+1];
if(tmpE2 + fML[indx[k]+1] + P->MLclosing <= threshold){
/* we've (hopefully) found a valid decomposition of fM2 and therefor we have all */
/* three intervals for our new state to be pushed on stack R */
new_state = copy_state(state);
/* first interval leads for search in fML array */
new_interval = make_interval(1, k, 1);
push(new_state->Intervals, new_interval);
env->nopush = false;
/* next, we have the first interval that has to be traced in fM1 */
new_interval = make_interval(k+1, l, 3);
push(new_state->Intervals, new_interval);
env->nopush = false;
/* and the last of our three intervals is also one to be traced within fM1 array... */
new_interval = make_interval(l+1, j, 3);
push(new_state->Intervals, new_interval);
env->nopush = false;
/* mmh, we add the energy for closing the multiloop now... */
new_state->partial_energy += P->MLclosing;
/* next we push our state onto the R stack */
push(env->Stack, new_state);
env->nopush = false;
}
/* else we search further... */
}
/* ok, we have to decompose fML now... */
}
}
}
} /* thats all folks for the circular case... */
/* 44444444444444444444444444444444444444444444444444 */
/* */
/* array_flag = 4: interval i,j was found while */
/* tracing back through fc-array smaller than than cp */
/* or within this block */
/* */
/* 44444444444444444444444444444444444444444444444444 */
if (array_flag == 4) {
int ik, s5, s3, tmp_en;
if(hc->up_ext[i]){
tmp_en = 0;
if(sc){
if(sc->energy_up)
tmp_en += sc->energy_up[i][1];
}
if (fc[i+1] + tmp_en + best_energy <= threshold) {
/* no basepair, nibbling of 5'-end */
fork_state(i+1, j, state, tmp_en, 4, env);
}
}
for (k = i+TURN+1; k < j; k++) {
ik = indx[k] + i;
if(with_gquad){
if(fc[k+1] + ggg[ik] + best_energy <= threshold){
temp_state = derive_new_state(k+1, j, state, 0, 4);
env->nopush = false;
repeat_gquad(vc, i, k, temp_state, 0, fc[k+1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
if(hard_constraints[ik] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[ik];
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s5 = s3 = -1;
break;
default: s5 = (i > 1) ? S1[i-1]: -1;
s3 = S1[k+1];
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
/* should be unnecessary
if(sc){
if(sc->energy_bp)
element_energy += sc->energy_bp[ik];
}
*/
if (fc[k+1] + c[ik] + element_energy + best_energy <= threshold){
temp_state = derive_new_state(k+1, j, state, 0, 4);
env->nopush = false;
repeat(vc, i, k, temp_state, element_energy, fc[k+1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
ik = indx[cp -1] + i; /* indx[j] + i; */
if(with_gquad){
if(ggg[ik] + best_energy <= threshold)
repeat_gquad(vc, i, cp - 1, state, 0, 0, best_energy, threshold, env);
}
if(hard_constraints[ik] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[ik];
s3 = -1;
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s5 = -1;
break;
default: s5 = (i>1) ? S1[i-1] : -1;
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
/* should be unnecessary
if(sc){
if(sc->energy_bp)
element_energy += sc->energy_bp[ik];
}
*/
if(c[ik] + element_energy + best_energy <= threshold)
repeat(vc, i, cp-1, state, element_energy, 0, best_energy, threshold, env);
}
} /* array_flag == 4 */
/* 55555555555555555555555555555555555555555555555555 */
/* */
/* array_flag = 5: interval cp=i,j was found while */
/* tracing back through fc-array greater than cp */
/* or within this block */
/* */
/* 55555555555555555555555555555555555555555555555555 */
if (array_flag == 5) {
int kj, s5, s3, tmp_en;
if(hc->up_ext[j]){
tmp_en = 0;
if(sc){
if(sc->energy_up)
tmp_en += sc->energy_up[j][1];
}
if (fc[j-1] + tmp_en + best_energy <= threshold) {
/* no basepair, nibbling of 3'-end */
fork_state(i, j-1, state, tmp_en, 5, env);
}
}
for (k = j-TURN-1; k > i; k--) {
kj = indx[j] + k;
if(with_gquad){
if(fc[k-1] + ggg[kj] + best_energy <= threshold){
temp_state = derive_new_state(i, k-1, state, 0, 5);
env->nopush = false;
repeat_gquad(vc, k, j, temp_state, 0, fc[k-1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
if(hard_constraints[kj] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[kj];
element_energy = 0;
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s3 = s5 = -1;
break;
default: s5 = S1[k-1];
s3 = (j < length) ? S1[j+1] : -1;
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
/* should be unnecessary
if(sc){
if(sc->energy_bp)
element_energy += sc->energy_bp[kj];
}
*/
if (fc[k-1] + c[kj] + element_energy + best_energy <= threshold) {
temp_state = derive_new_state(i, k-1, state, 0, 5);
env->nopush = false;
repeat(vc, k, j, temp_state, element_energy, fc[k-1], best_energy, threshold, env);
free_state_node(temp_state);
}
}
}
kj = indx[j] + cp; /* indx[j] + i; */
if(with_gquad){
if(ggg[kj] + best_energy <= threshold)
repeat_gquad(vc, cp, j, state, 0, 0, best_energy, threshold, env);
}
if(hard_constraints[kj] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
type = ptype[kj];
s5 = -1;
if(type == 0)
type = 7;
switch(dangle_model){
case 0: s3 = -1;
break;
default: s3 = (j<length) ? S1[j+1] : -1;
break;
}
element_energy = E_ExtLoop(type, s5, s3, P);
if (c[kj] + element_energy + best_energy <= threshold)
repeat(vc, cp, j, state, element_energy, 0, best_energy, threshold, env);
}
} /* array_flag == 5 */
if (array_flag == 6) { /* we have a gquad */
repeat_gquad(vc, i, j, state, 0, 0, best_energy, threshold, env);
if (env->nopush){
vrna_message_warning("%d,%d\nOops, no solution in gquad-repeat!", i, j);
}
return;
}
if (env->nopush){
push_back(env->Stack, state);
env->nopush = false;
}
return;
}
/*---------------------------------------------------------------------------*/
PRIVATE void
repeat_gquad( vrna_fold_compound_t *vc,
int i,
int j,
STATE *state,
int part_energy,
int temp_energy,
int best_energy,
int threshold,
subopt_env *env){
int *ggg, *indx, element_energy, cp;
short *S1;
vrna_param_t *P;
indx = vc->jindx;
cp = vc->cutpoint;
ggg = vc->matrices->ggg;
S1 = vc->sequence_encoding;
P = vc->params;
/* find all gquads that fit into the energy range and the interval [i,j] */
STATE *new_state;
best_energy += part_energy; /* energy of current structural element */
best_energy += temp_energy; /* energy from unpushed interval */
if(ON_SAME_STRAND(i,j,cp)){
element_energy = ggg[indx[j] + i];
if(element_energy + best_energy <= threshold){
int cnt;
int *L;
int *l;
/* find out how many gquads we might expect in the interval [i,j] */
int num_gquads = get_gquad_count(S1, i, j);
num_gquads++;
L = (int *)vrna_alloc(sizeof(int) * num_gquads);
l = (int *)vrna_alloc(sizeof(int) * num_gquads * 3);
L[0] = -1;
get_gquad_pattern_exhaustive(S1, i, j, P, L, l, threshold - best_energy);
for(cnt = 0; L[cnt] != -1; cnt++){
new_state = copy_state(state);
make_gquad(i, L[cnt], &(l[3*cnt]), new_state);
new_state->partial_energy += part_energy;
new_state->partial_energy += element_energy;
/* new_state->best_energy =
hairpin[unpaired] + element_energy + best_energy; */
push(env->Stack, new_state);
env->nopush = false;
}
free(L);
free(l);
}
}
best_energy -= part_energy;
best_energy -= temp_energy;
return;
}
PRIVATE void
repeat( vrna_fold_compound_t *vc,
int i,
int j,
STATE * state,
int part_energy,
int temp_energy,
int best_energy,
int threshold,
subopt_env *env){
/* routine to find stacks, bulges, internal loops and multiloops */
/* within interval closed by basepair i,j */
STATE *new_state;
vrna_param_t *P;
vrna_md_t *md;
register int ij, k, p, q, energy, new;
register int mm;
register int no_close, type, type_2;
char *ptype;
int element_energy;
int *fc, *c, *fML, *fM1, *ggg;
int rt, *indx, *rtype, noGUclosure, noLP, with_gquad, dangle_model, turn, cp;
short *S1;
vrna_hc_t *hc;
vrna_sc_t *sc;
S1 = vc->sequence_encoding;
ptype = vc->ptype;
indx = vc->jindx;
cp = vc->cutpoint;
P = vc->params;
md = &(P->model_details);
rtype = &(md->rtype[0]);
noGUclosure = md->noGUclosure;
noLP = md->noLP;
with_gquad = md->gquad;
dangle_model = md->dangles;
turn = md->min_loop_size;
fc = vc->matrices->fc;
c = vc->matrices->c;
fML = vc->matrices->fML;
fM1 = vc->matrices->fM1;
ggg = vc->matrices->ggg;
hc = vc->hc;
sc = vc->sc;
ij = indx[j]+i;
type = ptype[ij];
/*
if (type==0) fprintf(stderr, "repeat: Warning: %d %d can't pair\n", i,j);
*/
if(type == 0)
type = 7;
no_close = (((type == 3) || (type == 4)) && noGUclosure);
if(hc->matrix[ij] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
if (noLP) /* always consider the structure with additional stack */
if(i + turn + 2 < j){
if(hc->matrix[indx[j-1]+i+1] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC){
type_2 = rtype[ptype[indx[j-1]+i+1]];
if(type_2 == 0)
type_2 = 7;
energy = 0;
if(ON_SAME_STRAND(i,i+1,cp) && ON_SAME_STRAND(j-1,j, cp)){
energy = E_IntLoop(0, 0, type, type_2,S1[i+1],S1[j-1],S1[i+1],S1[j-1], P);
if(sc){
if(sc->energy_bp)
energy += sc->energy_bp[ij];
if(sc->energy_stack)
energy += sc->energy_stack[i]
+ sc->energy_stack[i+1]
+ sc->energy_stack[j-1]
+ sc->energy_stack[j];
if(sc->f)
energy += sc->f(i, j, i+1, j-1, VRNA_DECOMP_PAIR_IL, sc->data);
}
new_state = derive_new_state(i+1, j-1, state, part_energy + energy, 2);
make_pair(i, j, new_state);
make_pair(i+1, j-1, new_state);
/* new_state->best_energy = new + best_energy; */
push(env->Stack, new_state);
env->nopush = false;
if (i==1 || state->structure[i-2]!='(' || state->structure[j]!=')')
/* adding a stack is the only possible structure */
return;
}
}
}
}
best_energy += part_energy; /* energy of current structural element */
best_energy += temp_energy; /* energy from unpushed interval */
if(hc->matrix[ij] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
for (p = i + 1; p <= MIN2 (j-2-turn, i+MAXLOOP+1); p++) {
int minq = j-i+p-MAXLOOP-2;
if (minq<p+1+turn) minq = p+1+turn;
if(hc->up_int[i+1] < (p - i - 1))
break;
for (q = j - 1; q >= minq; q--) {
if(hc->up_int[q+1] < (j - q - 1))
break;
/* skip stack if noLP, since we've already processed it above */
if((noLP) && (p==i+1) && (q==j-1))
continue;
if(!(hc->matrix[indx[q]+p] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC))
continue;
type_2 = ptype[indx[q]+p];
if(type_2 == 0)
type_2 = 7;
if (noGUclosure)
if (no_close||(type_2==3)||(type_2==4))
if ((p>i+1)||(q<j-1)) continue; /* continue unless stack */
if (ON_SAME_STRAND(i,p,cp) && ON_SAME_STRAND(q,j,cp)) {
energy = E_IntLoop(p-i-1, j-q-1, type, rtype[type_2],
S1[i+1],S1[j-1],S1[p-1],S1[q+1], P);
new = energy + c[indx[q]+p];
if(sc){
if(sc->energy_up)
energy += sc->energy_up[i+1][p-i-1]
+ sc->energy_up[q+1][j-q-1];
if(sc->energy_bp)
energy += sc->energy_bp[ij];
if(sc->energy_stack)
if((p == i+1) && (q == j-1))
energy += sc->energy_stack[i]
+ sc->energy_stack[p]
+ sc->energy_stack[q]
+ sc->energy_stack[j];
if(sc->f)
energy += sc->f(i, j, p, q, VRNA_DECOMP_PAIR_IL, sc->data);
}
new = energy + c[indx[q]+p];
if (new + best_energy <= threshold) {
/* stack, bulge, or interior loop */
fork_int_state(i, j, p, q, state, part_energy + energy, env);
}
}/*end of if block */
} /* end of q-loop */
} /* end of p-loop */
}
if (!ON_SAME_STRAND(i,j,cp)) { /*look in fc*/
if(hc->matrix[ij] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
rt = rtype[type];
if(rt == 0)
rt = 7;
element_energy=0;
switch(dangle_model){
case 0: element_energy = E_ExtLoop(rt, -1, -1, P);
break;
default: element_energy = E_ExtLoop(rt, (ON_SAME_STRAND(j-1,j,cp)) ? S1[j-1] : -1, (ON_SAME_STRAND(i,i+1,cp)) ? S1[i+1] : -1, P);
break;
}
if (fc[i+1] + fc[j-1] +element_energy + best_energy <= threshold)
{
fork_two_states_pair(i, j, cp, state, part_energy + element_energy, 4, 5, env);
}
}
}
mm = P->MLclosing;
rt = rtype[type];
if((hc->matrix[ij] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP) && (i != cp-1) && (j != cp)){
if(rt == 0)
rt = 7;
element_energy = mm;
switch(dangle_model){
case 0: element_energy = E_MLstem(rt, -1, -1, P) + mm;
break;
default: element_energy = E_MLstem(rt, S1[j-1], S1[i+1], P) + mm;
break;
}
if(sc){
if(sc->energy_bp)
element_energy += sc->energy_bp[ij];
}
for (k = i + turn + 2; k <= j - turn - 2; k++) {
/* multiloop decomposition */
if ((fML[indx[k-1] + i+1] + fM1[indx[j-1] + k] +
element_energy + best_energy) <= threshold)
{
fork_two_states_pair(i, j, k, state, part_energy + element_energy, 1, 3, env);
}
}
}
if (ON_SAME_STRAND(i,j,cp)) {
if(hc->matrix[ij] & VRNA_CONSTRAINT_CONTEXT_HP_LOOP){
if(no_close)
element_energy = FORBIDDEN;
else
element_energy = vrna_E_hp_loop(vc, i, j);
if (element_energy + best_energy <= threshold) {
/* hairpin structure */
fork_state_pair(i, j, state, part_energy + element_energy, env);
}
}
if(with_gquad){
/* now we have to find all loops where (i,j) encloses a gquad in an interior loops style */
int cnt, *p, *q, *en, tmp_en;
p = q = en = NULL;
en = E_GQuad_IntLoop_exhaustive(i, j, &p, &q, type, S1, ggg, threshold - best_energy, indx, P);
for(cnt = 0; p[cnt] != -1; cnt++){
if((hc->up_int[i+1] >= p[cnt] - i - 1) && (hc->up_int[q[cnt]+1] >= j - q[cnt] - 1)){
tmp_en = en[cnt];
if(sc){
if(sc->energy_bp)
tmp_en += sc->energy_bp[ij];
if(sc->energy_up)
tmp_en += sc->energy_up[i+1][p[cnt] - i - 1]
+ sc->energy_up[q[cnt]+1][j - q[cnt] - 1];
}
new_state = derive_new_state(p[cnt], q[cnt], state, tmp_en + part_energy, 6);
make_pair(i, j, new_state);
/* new_state->best_energy = new + best_energy; */
push(env->Stack, new_state);
env->nopush = false;
}
}
free(en);
free(p);
free(q);
}
}
best_energy -= part_energy;
best_energy -= temp_energy;
return;
}
PRIVATE void
old_subopt_print( const char *structure,
float energy,
void *data){
struct old_subopt_dat *d = (struct old_subopt_dat *)data;
if(structure && d->fp){
char *e_string = vrna_strdup_printf(" %6.2f", energy);
print_structure(d->fp, structure, e_string);
free(e_string);
}
}
PRIVATE void
old_subopt_store( const char *structure,
float energy,
void *data){
struct old_subopt_dat *d = (struct old_subopt_dat *)data;
/* store solution */
if(d->n_sol + 1 == d->max_sol){
d->max_sol *= 2;
d->SolutionList = (SOLUTION *)vrna_realloc(d->SolutionList, d->max_sol*sizeof(SOLUTION));
}
if(structure){
d->SolutionList[d->n_sol].energy = energy;
d->SolutionList[d->n_sol++].structure = strdup(structure);
} else {
d->SolutionList[d->n_sol].energy = 0;
d->SolutionList[d->n_sol++].structure = NULL;
}
}
/*###########################################*/
/*# deprecated functions below #*/
/*###########################################*/
#ifdef VRNA_BACKWARD_COMPAT
PUBLIC SOLUTION *
subopt( char *seq,
char *structure,
int delta,
FILE *fp){
return wrap_subopt(seq, structure, NULL, delta, fold_constrained, 0, fp);
}
PUBLIC SOLUTION *
subopt_circ(char *seq,
char *structure,
int delta,
FILE *fp){
return wrap_subopt(seq, structure, NULL, delta, fold_constrained, 1, fp);
}
PUBLIC SOLUTION *subopt_par(char *seq,
char *structure,
vrna_param_t *parameters,
int delta,
int is_constrained,
int is_circular,
FILE *fp){
return wrap_subopt(seq, structure, parameters, delta, is_constrained, is_circular, fp);
}
PRIVATE SOLUTION *
wrap_subopt(char *string,
char *structure,
vrna_param_t *parameters,
int delta,
int is_constrained,
int is_circular,
FILE *fp){
vrna_fold_compound_t *vc;
vrna_param_t *P;
char *seq;
#ifdef _OPENMP
/* Explicitly turn off dynamic threads */
omp_set_dynamic(0);
#endif
/* we need the parameter structure for hard constraints */
if(parameters){
P = vrna_params_copy(parameters);
} else {
vrna_md_t md;
set_model_details(&md);
md.temperature = temperature;
P = vrna_params(&md);
}
P->model_details.circ = is_circular;
P->model_details.uniq_ML = uniq_ML = 1;
/* what about cofold sequences here? Is it safe to call the below cut_point_insert() ? */
/* dirty hack to reinsert the '&' according to the global variable 'cut_point' */
seq = vrna_cut_point_insert(string, cut_point);
vc = vrna_fold_compound(seq, &(P->model_details), ((is_circular == 0) ? VRNA_OPTION_HYBRID : VRNA_OPTION_DEFAULT));
if(parameters){ /* replace params if necessary */
free(vc->params);
vc->params = P;
} else {
free(P);
}
/* handle hard constraints in pseudo dot-bracket format if passed via simple interface */
if(is_constrained && structure){
unsigned int constraint_options = 0;
constraint_options |= VRNA_CONSTRAINT_DB
| VRNA_CONSTRAINT_DB_PIPE
| VRNA_CONSTRAINT_DB_DOT
| VRNA_CONSTRAINT_DB_X
| VRNA_CONSTRAINT_DB_ANG_BRACK
| VRNA_CONSTRAINT_DB_RND_BRACK
| VRNA_CONSTRAINT_DB_INTRAMOL
| VRNA_CONSTRAINT_DB_INTERMOL;
vrna_constraints_add(vc, (const char *)structure, constraint_options);
}
if(backward_compat_compound && backward_compat)
vrna_fold_compound_free(backward_compat_compound);
backward_compat_compound = vc;
backward_compat = 1;
/* cleanup */
free(seq);
return vrna_subopt(vc, delta, subopt_sorted, fp);
}
#endif
/*---------------------------------------------------------------------------*/
/* Well, that is the end!----------------------------------------------------*/
/*---------------------------------------------------------------------------*/