ViennaRNA-bindings-0.233.1.1: C/ViennaRNA/part_func_co.c
/*
partiton function for RNA secondary structures
Ivo L Hofacker
Stephan Bernhart
Ronny Lorenz
Vienna RNA package
*/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h> /* #defines FLT_MAX ... */
#include <limits.h>
#include "ViennaRNA/utils.h"
#include "ViennaRNA/structure_utils.h"
#include "ViennaRNA/energy_par.h"
#include "ViennaRNA/fold_vars.h"
#include "ViennaRNA/PS_dot.h"
#include "ViennaRNA/params.h"
#include "ViennaRNA/loop_energies.h"
#include "ViennaRNA/part_func.h"
#include "ViennaRNA/part_func_co.h"
#ifdef _OPENMP
#include <omp.h>
#endif
/*
#################################
# GLOBAL VARIABLES #
#################################
*/
int mirnatog = 0;
double F_monomer[2] = {0,0}; /* free energies of the two monomers */
/*
#################################
# PRIVATE VARIABLES #
#################################
*/
/* some backward compatibility stuff */
PRIVATE vrna_fold_compound_t *backward_compat_compound = NULL;
PRIVATE int backward_compat = 0;
#ifdef _OPENMP
#pragma omp threadprivate(backward_compat_compound, backward_compat)
#endif
/*
#################################
# PRIVATE FUNCTION DECLARATIONS #
#################################
*/
PRIVATE void pf_co(vrna_fold_compound_t *vc);
PRIVATE void pf_co_bppm(vrna_fold_compound_t *vc, char *structure);
PRIVATE double *Newton_Conc(double ZAB, double ZAA, double ZBB, double concA, double concB,double* ConcVec);
PRIVATE vrna_dimer_pf_t wrap_co_pf_fold(char *sequence,
char *structure,
vrna_exp_param_t *parameters,
int calculate_bppm,
int is_constrained);
/*
#################################
# BEGIN OF FUNCTION DEFINITIONS #
#################################
*/
/*
*****************************************
* BEGIN backward compatibility wrappers *
*****************************************
*/
PRIVATE vrna_dimer_pf_t
wrap_co_pf_fold(char *sequence,
char *structure,
vrna_exp_param_t *parameters,
int calculate_bppm,
int is_constrained){
int length;
char *seq;
vrna_fold_compound_t *vc;
vrna_md_t md;
vc = NULL;
length = strlen(sequence);
/* we need vrna_exp_param_t datastructure to correctly init default hard constraints */
if(parameters)
md = parameters->model_details;
else{
set_model_details(&md); /* get global default parameters */
}
md.compute_bpp = calculate_bppm;
md.min_loop_size = 0;
seq = (char *)vrna_alloc(sizeof(char) * (length + 2));
if(cut_point > -1){
int i;
for(i = 0; i < cut_point-1; i++)
seq[i] = sequence[i];
seq[i] = '&';
for(;i<(int)length;i++)
seq[i+1] = sequence[i];
} else { /* this ensures the allocation of all cofold matrices via vrna_fold_compound_t */
free(seq);
seq = strdup(sequence);
}
vc = vrna_fold_compound(seq, &md, VRNA_OPTION_PF | VRNA_OPTION_HYBRID);
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_constraints_add(vc, (const char *)structure, constraint_options);
}
if(backward_compat_compound)
vrna_fold_compound_free(backward_compat_compound);
backward_compat_compound = vc;
backward_compat = 1;
iindx = backward_compat_compound->iindx;
free(seq);
return vrna_pf_dimer(vc, structure);
}
/*
*****************************************
* END backward compatibility wrappers *
*****************************************
*/
PUBLIC vrna_dimer_pf_t
vrna_pf_dimer(vrna_fold_compound_t *vc,
char *structure){
int n;
FLT_OR_DBL Q;
vrna_dimer_pf_t X;
double free_energy;
char *sequence;
vrna_md_t *md;
vrna_exp_param_t *params;
vrna_mx_pf_t *matrices;
if(!vrna_fold_compound_prepare(vc, VRNA_OPTION_PF | VRNA_OPTION_HYBRID)){
vrna_message_warning("vrna_pf_dimer@part_func_co.c: Failed to prepare vrna_fold_compound");
X.FA = X.FB = X.FAB = X.F0AB = X.FcAB = 0;
return X;
}
params = vc->exp_params;
n = vc->length;
md = &(params->model_details);
matrices = vc->exp_matrices;
sequence = vc->sequence;
#ifdef _OPENMP
/* Explicitly turn off dynamic threads */
omp_set_dynamic(0);
#endif
#ifdef SUN4
nonstandard_arithmetic();
#else
#ifdef HP9
fpsetfastmode(1);
#endif
#endif
/* call user-defined recursion status callback function */
if(vc->stat_cb)
vc->stat_cb(VRNA_STATUS_PF_PRE, vc->auxdata);
pf_co(vc);
/* call user-defined recursion status callback function */
if(vc->stat_cb)
vc->stat_cb(VRNA_STATUS_PF_POST, vc->auxdata);
if (md->backtrack_type=='C')
Q = matrices->qb[vc->iindx[1]-n];
else if (md->backtrack_type=='M')
Q = matrices->qm[vc->iindx[1]-n];
else Q = matrices->q[vc->iindx[1]-n];
/* ensemble free energy in Kcal/mol */
if (Q<=FLT_MIN)
vrna_message_warning("pf_scale too large");
free_energy = (-log(Q)-n*log(params->pf_scale))*params->kT/1000.0;
/* in case we abort because of floating point errors */
if(n>1600)
vrna_message_info(stderr, "free energy = %8.2f", free_energy);
/*probability of molecules being bound together*/
/*Computation of "real" Partition function*/
/*Need that for concentrations*/
if (vc->cutpoint > 0){
double kT, QAB, QToT, Qzero;
kT = params->kT/1000.0;
Qzero = matrices->q[vc->iindx[1] - n];
QAB = (matrices->q[vc->iindx[1] - n]- matrices->q[vc->iindx[1] - (vc->cutpoint - 1)] * matrices->q[vc->iindx[vc->cutpoint] - n]) * params->expDuplexInit;
/*correction for symmetry*/
if((n - (vc->cutpoint - 1) * 2) == 0){
if((strncmp(sequence, sequence + vc->cutpoint - 1, vc->cutpoint - 1)) == 0){
QAB/=2;
}
}
QToT = matrices->q[vc->iindx[1] - (vc->cutpoint - 1)] * matrices->q[vc->iindx[vc->cutpoint] - n] + QAB;
X.FAB = -kT * (log(QToT) + n * log(params->pf_scale));
X.F0AB = -kT * (log(Qzero)+ n * log(params->pf_scale));
X.FcAB = (QAB>1e-17) ? -kT * (log(QAB) + n * log(params->pf_scale)) : 999;
X.FA = -kT * (log(matrices->q[vc->iindx[1] - (vc->cutpoint - 1)]) + (vc->cutpoint - 1) * log(params->pf_scale));
X.FB = -kT * (log(matrices->q[vc->iindx[vc->cutpoint] - n]) + (n - vc->cutpoint + 1) * log(params->pf_scale));
/* printf("QAB=%.9f\tQtot=%.9f\n",QAB/scale[n],QToT/scale[n]);*/
}
else {
X.FA = X.FB = X.FAB = X.F0AB = free_energy;
X.FcAB = 0;
}
/* backtracking to construct binding probabilities of pairs*/
if(md->compute_bpp){
pf_co_bppm(vc, structure);
/*
* Backward compatibility:
* This block may be removed if deprecated functions
* relying on the global variable "pr" vanish from within the package!
*/
pr = vc->exp_matrices->probs;
/*
{
if(pr) free(pr);
pr = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2));
memcpy(pr, probs, sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2));
}
*/
}
#ifdef SUN4
standard_arithmetic();
#else
#ifdef HP9
fpsetfastmode(0);
#endif
#endif
return X;
}
/* forward recursion of pf cofolding */
PRIVATE void
pf_co(vrna_fold_compound_t *vc){
unsigned int *sn;
int n, i,j,k,l, ij, kl, u,u1,u2,ii, type, type_2, tt, cp, turn, maxk, minl;
FLT_OR_DBL *qqm = NULL, *qqm1 = NULL, *qq = NULL, *qq1 = NULL;
FLT_OR_DBL temp, q_temp, Qmax=0;
FLT_OR_DBL qbt1, *tmp;
FLT_OR_DBL *q, *qb, *qm, *qm1;
FLT_OR_DBL *scale;
FLT_OR_DBL *expMLbase;
short *S1;
short s5, s3;
int *my_iindx, *jindx;
char *ptype, *sequence;
vrna_md_t *md;
vrna_hc_t *hc;
vrna_sc_t *sc;
FLT_OR_DBL expMLclosing;
int noGUclosure;
double max_real;
int *rtype;
vrna_exp_param_t *pf_params;
vrna_mx_pf_t *matrices;
int hc_decompose;
char *hard_constraints;
int *hc_up_ext;
int *hc_up_hp;
int *hc_up_int;
int *hc_up_ml;
sequence = vc->sequence;
S1 = vc->sequence_encoding;
n = vc->length;
cp = vc->cutpoint;
my_iindx = vc->iindx;
jindx = vc->jindx;
ptype = vc->ptype;
pf_params = vc->exp_params;
md = &(pf_params->model_details);
sn = vc->strand_number;
rtype = &(md->rtype[0]);
hc = vc->hc;
sc = vc->sc;
expMLclosing = pf_params->expMLclosing;
noGUclosure = md->noGUclosure;
matrices = vc->exp_matrices;
turn = md->min_loop_size;
q = matrices->q;
qb = matrices->qb;
qm = matrices->qm;
qm1 = matrices->qm1;
scale = matrices->scale;
expMLbase = matrices->expMLbase;
hard_constraints = hc->matrix;
hc_up_ext = hc->up_ext;
hc_up_hp = hc->up_hp;
hc_up_int = hc->up_int;
hc_up_ml = hc->up_ml;
max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
/* allocate memory for helper arrays */
qq = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
qq1 = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
qqm = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
qqm1 = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
/* hard code min_loop_size to 0, since we can not be sure yet that this is already the case */
turn = 0;
/*array initialization ; qb,qm,q
qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */
/* for (d=0; d<=TURN; d++) */
for (i=1; i<=n/*-d*/; i++) {
ij = my_iindx[i]-i;
if(hc_up_ext[i]){
q[ij] = scale[1];
if(sc){
if(sc->exp_energy_up)
q[ij] *= sc->exp_energy_up[i][1];
if(sc->exp_f)
q[ij] *= sc->exp_f(i, i, i, i, VRNA_DECOMP_EXT_UP, sc->data);
}
} else {
q[ij] = 0.;
}
qb[ij] = qm[ij] = 0.0;
}
for (i=0; i<=n; i++)
qq[i] = qq1[i] = qqm[i] = qqm1[i] = 0;
for (j = turn + 2; j <= n; j++) {
for (i = j - turn - 1; i >= 1; i--) {
/* construction of partition function of segment i,j */
/* firstly that given i binds j : qb(i,j) */
u = j - i - 1;
ij = my_iindx[i] - j;
type = (unsigned char)ptype[jindx[j] + i];
hc_decompose = hard_constraints[jindx[j] + i];
qbt1 = 0;
q_temp = 0.;
if(hc_decompose){
/* process hairpin loop(s) */
qbt1 += vrna_exp_E_hp_loop(vc, i, j);
qbt1 += vrna_exp_E_int_loop(vc, i, j);
qbt1 += vrna_exp_E_mb_loop_fast(vc, i, j, qqm1);
qb[ij] = qbt1;
} else /* end if allowed to be paired */
qb[ij] = 0.0;
/* construction of qqm matrix containing final stem
contributions to multiple loop partition function
from segment i,j */
qqm[i] = 0.;
if(hc_up_ml[j]){
if (sn[j] == sn[j - 1]) {
q_temp = qqm1[i] * expMLbase[1];
if(sc){
if(sc->exp_energy_up)
q_temp *= sc->exp_energy_up[j][1];
if(sc->exp_f)
q_temp *= sc->exp_f(i, j, i, j-1, VRNA_DECOMP_ML_ML, sc->data);
}
qqm[i] = q_temp;
}
}
if(hc_decompose & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC){
if ((sn[i] == sn[i - 1]) && (sn[j + 1] == sn[j])) {
tt = type;
if(tt == 0)
tt = 7;
qbt1 = qb[ij];
qbt1 *= exp_E_MLstem(tt, (i>1) ? S1[i-1] : -1, (j<n) ? S1[j+1] : -1, pf_params);
if(sc){
if(sc->exp_f)
q_temp *= sc->exp_f(i, j, i, j, VRNA_DECOMP_ML_STEM, sc->data);
}
qqm[i] += qbt1;
}
}
if (qm1) qm1[jindx[j]+i] = qqm[i]; /* for stochastic backtracking */
/*construction of qm matrix containing multiple loop
partition function contributions from segment i,j */
temp = 0.0;
kl = my_iindx[i] - j + 1; /* ii-k=[i,k-1] */
if (sc && sc->exp_f) {
if (j >= cp) {
for (k = j; k > MAX2(i, cp); k--, kl++) {
q_temp = qm[kl] * qqm[k];
q_temp *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data);
temp += q_temp;
}
for (; k > i; k--, kl++) {
q_temp = qm[kl] * qqm[k];
q_temp *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data);
temp += q_temp;
}
} else {
for (k = j; k > i; k--, kl++) {
q_temp = qm[kl] * qqm[k];
q_temp *= sc->exp_f(i, j, k - 1, k, VRNA_DECOMP_ML_ML_ML, sc->data);
temp += q_temp;
}
}
} else { /* without soft-constraints */
if (j >= cp) {
for (k = j; k > MAX2(i, cp); k--, kl++)
temp += qm[kl] * qqm[k];
for (; k > i; k--, kl++)
temp += qm[kl] * qqm[k];
} else {
for (k = j; k > i; k--, kl++)
temp += qm[kl] * qqm[k];
}
}
maxk = MIN2(i + hc_up_ml[i], j);
if (i < cp) { /* we must not have the strand border within the unpaired segment */
maxk = MIN2(maxk, cp - 1);
}
u = 1; /* length of unpaired stretch */
if (sc) {
for (k = i + 1; k <= maxk; k++, u++){
q_temp = expMLbase[u] * qqm[k];
if(sc->exp_energy_up)
q_temp *= sc->exp_energy_up[i][u];
if(sc->exp_f)
q_temp *= sc->exp_f(i, j, k, j, VRNA_DECOMP_ML_ML, sc->data);
temp += q_temp;
}
} else { /* without soft-constraints */
for (k = i + 1; k <= maxk; k++, u++)
temp += expMLbase[u] * qqm[k];
}
qm[ij] = (temp + qqm[i]);
/*auxiliary matrix qq for cubic order q calculation below */
qbt1 = 0.;
if(hc_decompose & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP){
tt = type;
if(tt == 0)
tt = 7;
s5 = ((i > 1) && (sn[i] == sn[i - 1])) ? S1[i - 1] : -1;
s3 = ((j < n) && (sn[j + 1] == sn[j])) ? S1[j + 1] : -1;
qbt1 = qb[ij] * exp_E_ExtLoop(tt, s5, s3, pf_params);
if(sc){
if(sc->exp_f)
qbt1 *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_STEM, sc->data);
}
}
if(hc_up_ext[j]){
q_temp = qq1[i] * scale[1];
if(sc){
if(sc->exp_energy_up)
q_temp *= sc->exp_energy_up[j][1];
if(sc->exp_f)
q_temp *= sc->exp_f(i, j, i, j-1, VRNA_DECOMP_EXT_EXT, sc->data);
}
qbt1 += q_temp;
}
qq[i] = qbt1;
/*construction of partition function for segment i,j */
temp = qq[i];
/* the whole stretch [i,j] is unpaired */
if(hc_up_ext[i] >= (j-i+1)){
q_temp = 1.0 * scale[j-i+1];
if(sc){
if(sc->exp_energy_up)
q_temp *= sc->exp_energy_up[i][j-i+1];
if(sc->exp_f)
q_temp *= sc->exp_f(i, j, i, j, VRNA_DECOMP_EXT_UP, sc->data);
}
temp += q_temp;
}
kl = my_iindx[i] - i;
if (sc && sc->exp_f) {
for (k=i; k<j; k++, kl--){
q_temp = q[kl] * qq[k+1];
q_temp *= sc->exp_f(i, j, k, k+1, VRNA_DECOMP_EXT_EXT_EXT, sc->data);
temp += q_temp;
}
} else {
for (k=i; k<j; k++, kl--)
temp += q[kl] * qq[k+1];
}
q[ij] = temp;
if (temp>Qmax) {
Qmax = temp;
if (Qmax>max_real/10.)
vrna_message_warning("Q close to overflow: %d %d %g", i,j,temp);
}
if (temp>=max_real) {
vrna_message_error("overflow in co_pf_fold while calculating q[%d,%d]\n"
"use larger pf_scale", i,j);
}
}
tmp = qq1; qq1 =qq; qq =tmp;
tmp = qqm1; qqm1=qqm; qqm=tmp;
}
/* clean up */
free(qq);
free(qq1);
free(qqm);
free(qqm1);
}
/* backward recursion of pf cofolding */
PRIVATE void
pf_co_bppm(vrna_fold_compound_t *vc, char *structure){
unsigned int *sn;
int n, i,j,k,l, ij, kl, ii, ll, lj, u1, u2, type, type_2, tt, turn, ov=0, *my_iindx, *jindx, cp;
FLT_OR_DBL temp, Qmax=0, prm_MLb, tmp2, ppp;
FLT_OR_DBL prmt,prmt1, *expMLbase;
FLT_OR_DBL *tmp;
FLT_OR_DBL expMLclosing, *probs, *q1k, *qln, *q, *qb, *qm, *scale;
double max_real;
vrna_exp_param_t *pf_params;
vrna_md_t *md;
short *S, *S1, s5, s3;
char *ptype;
vrna_hc_t *hc;
vrna_sc_t *sc;
vrna_mx_pf_t *matrices;
char *sequence;
char *hard_constraints;
int *hc_up_ext;
int *hc_up_hp;
int *hc_up_int;
int *hc_up_ml;
int *rtype;
sequence = vc->sequence;
n = vc->length;
cp = vc->cutpoint;
pf_params = vc->exp_params;
md = &(pf_params->model_details);
expMLclosing = pf_params->expMLclosing;
S = vc->sequence_encoding2;
S1 = vc->sequence_encoding;
sn = vc->strand_number;
jindx = vc->jindx;
my_iindx = vc->iindx;
ptype = vc->ptype;
rtype = &(md->rtype[0]);
turn = md->min_loop_size;
matrices = vc->exp_matrices;
probs = matrices->probs;
scale = matrices->scale;
q1k = matrices->q1k;
qln = matrices->qln;
q = matrices->q;
qb = matrices->qb;
qm = matrices->qm;
expMLbase = matrices->expMLbase;
hc = vc->hc;
sc = vc->sc;
hard_constraints = hc->matrix;
hc_up_ext = hc->up_ext;
hc_up_hp = hc->up_hp;
hc_up_int = hc->up_int;
hc_up_ml = hc->up_ml;
/* hard code min_loop_size to 0, since we can not be sure yet that this is already the case */
turn = 0;
max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
/* backtracking to construct binding probabilities of pairs*/
if ((S != NULL) && (S1 != NULL)) {
FLT_OR_DBL *Qlout, *Qrout;
FLT_OR_DBL *prm_l = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
FLT_OR_DBL *prm_l1 = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
FLT_OR_DBL *prml = (FLT_OR_DBL *) vrna_alloc(sizeof(FLT_OR_DBL)*(n+2));
Qmax = 0;
Qrout = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (n+2));
Qlout = (FLT_OR_DBL *)vrna_alloc(sizeof(FLT_OR_DBL) * (cp+2));
for (k=1; k<=n; k++) {
q1k[k] = q[my_iindx[1] - k];
qln[k] = q[my_iindx[k] - n];
}
q1k[0] = 1.0;
qln[n+1] = 1.0;
/* 1. exterior pair i,j and initialization of pr array */
for (i=1; i<=n; i++) {
for (j=i; j<=MIN2(i + turn, n); j++)
probs[my_iindx[i]-j] = 0;
for (j = i + turn + 1; j <= n; j++){
ij = my_iindx[i]-j;
if((hard_constraints[jindx[j] + i] & VRNA_CONSTRAINT_CONTEXT_EXT_LOOP) && (qb[ij] > 0.)){
type = ptype[jindx[j] + i];
if(type == 0)
type = 7;
s5 = ((i > 1) && (sn[i] == sn[i - 1])) ? S1[i - 1] : -1;
s3 = ((j < n) && (sn[j + 1] == sn[j])) ? S1[j + 1] : -1;
probs[ij] = q1k[i - 1] * qln[j + 1] / q1k[n];
probs[ij] *= exp_E_ExtLoop(type, s5, s3, pf_params);
if(sc){
if(sc->exp_f){
probs[ij] *= sc->exp_f(1, n, i, j, VRNA_DECOMP_EXT_STEM_OUTSIDE, sc->data);
}
}
} else
probs[ij] = 0;
}
}
for(l = n; l > turn + 1; l--){
/* 2. bonding k,l as substem of 2:loop enclosed by i,j */
for(k = 1; k < l - turn; k++){
kl = my_iindx[k]-l;
type_2 = (unsigned char)ptype[jindx[l] + k];
type_2 = rtype[type_2];
if(qb[kl]==0.) continue;
if(hard_constraints[jindx[l] + k] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP_ENC){
if(type_2 == 0)
type_2 = 7;
for(i = MAX2(1, k - MAXLOOP - 1); i <= k - 1; i++){
u1 = k - i - 1;
if(hc_up_int[i+1] < u1) continue;
for(j = l + 1; j <= MIN2(l + MAXLOOP - k + i + 2, n); j++){
u2 = j-l-1;
if(hc_up_int[l+1] < u2) break;
ij = my_iindx[i] - j;
if(hard_constraints[jindx[j] + i] & VRNA_CONSTRAINT_CONTEXT_INT_LOOP){
if ((sn[k] == sn[i]) && (sn[j] == sn[l])) {
type = (unsigned char)ptype[jindx[j] + i];
if(probs[ij] > 0){
if(type == 0)
type = 7;
tmp2 = probs[ij]
* scale[u1 + u2 + 2]
* exp_E_IntLoop(u1, u2, type, type_2, S1[i+1], S1[j-1], S1[k-1], S1[l+1], pf_params);
if(sc){
if(sc->exp_energy_up)
tmp2 *= sc->exp_energy_up[i+1][u1]
* sc->exp_energy_up[l+1][u2];
if(sc->exp_energy_bp)
tmp2 *= sc->exp_energy_bp[ij];
if(sc->exp_energy_stack){
if((i+1 == k) && (j-1 == l)){
tmp2 *= sc->exp_energy_stack[i]
* sc->exp_energy_stack[k]
* sc->exp_energy_stack[l]
* sc->exp_energy_stack[j];
}
}
if(sc->exp_f)
tmp2 *= sc->exp_f(i, j, k, l, VRNA_DECOMP_PAIR_IL, sc->data);
}
probs[kl] += tmp2;
}
}
}
}
}
}
}
/* 3. bonding k,l as substem of multi-loop enclosed by i,j */
prm_MLb = 0.;
if ((l < n) && (sn[l + 1] == sn[l]))
for (k = 2; k < l - turn; k++) {
kl = my_iindx[k] - l;
i = k - 1;
prmt = prmt1 = 0.0;
ii = my_iindx[i]; /* ii-j=[i,j] */
ll = my_iindx[l+1]; /* ll-j=[l+1,j] */
tt = (unsigned char)ptype[jindx[l+1] + i];
tt = rtype[tt];
if (sn[k] == sn[i]) {
if(hard_constraints[jindx[l+1] + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
if(tt == 0)
tt = 7;
prmt1 = probs[ii-(l+1)]
* expMLclosing
* exp_E_MLstem(tt, S1[l], S1[i+1], pf_params);
if(sc){
/* which decompositions are covered here? => (i, l+1) -> enclosing pair, (k,l) -> enclosed pair, */
if(sc->exp_energy_bp)
prmt1 *= sc->exp_energy_bp[ii - (l+1)];
/*
if(sc->exp_f)
prmt1 *= sc->exp_f(i, l+1, k, l, , sc->data);
*/
}
}
ij = my_iindx[i] - (l+2);
lj = my_iindx[l+1]-(l+1);
for(j = l + 2; j <= n; j++, ij--, lj--){
if(hard_constraints[jindx[j] + i] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
if (sn[j] == sn[j - 1]) { /*??*/
tt = (unsigned char)ptype[jindx[j] + i];
tt = rtype[tt];
if(tt == 0)
tt = 7;
/* which decomposition is covered here? =>
i + 1 = k < l < j:
(i,j) -> enclosing pair
(k, l) -> enclosed pair
(l+1, j-1) -> multiloop part with at least one stem
*/
ppp = probs[ii-j]
* exp_E_MLstem(tt, S1[j-1], S1[i+1], pf_params)
* qm[ll-(j-1)];
if(sc){
if(sc->exp_energy_bp)
ppp *= sc->exp_energy_bp[ij];
/*
if(sc->exp_f)
ppp *= sc->exp_f(i, j, l+1, j-1, , sc->data);
*/
}
prmt += ppp;
}
}
}
}
prmt *= expMLclosing;
tt = ptype[jindx[l] + k];
prml[ i] = prmt;
/* l+1 is unpaired */
if(hc->up_ml[l+1]){
ppp = prm_l1[i] * expMLbase[1];
if(sc){
if(sc->exp_energy_up)
ppp *= sc->exp_energy_up[l+1][1];
/*
if(sc_exp_f)
ppp *= sc->exp_f(, sc->data);
*/
}
prm_l[i] = ppp + prmt1;
} else {
prm_l[i] = prmt1;
}
/* i is unpaired */
if(hc->up_ml[i]){
ppp = prm_MLb*expMLbase[1];
if(sc){
if(sc->exp_energy_up)
ppp *= sc->exp_energy_up[i][1];
/*
if(sc->exp_f)
ppp *= sc->exp_f(, sc->data);
*/
}
prm_MLb = ppp + prml[i];
/* same as: prm_MLb = 0;
for (i=1; i<=k-1; i++) prm_MLb += prml[i]*expMLbase[k-i-1]; */
} else {
prm_MLb = prml[i];
}
prml[i] = prml[ i] + prm_l[i];
if (qb[kl] == 0.) continue;
if(hard_constraints[jindx[l] + k] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP_ENC){
temp = prm_MLb;
for (i=1;i<=k-2; i++) {
if ((sn[i + 1] == sn[i]) && (sn[k] == sn[k - 1])) {
temp += prml[i]*qm[my_iindx[i+1] - (k-1)];
}
}
if(tt == 0)
tt = 7;
s5 = ((k > 1) && (sn[k] == sn[k - 1])) ? S1[k - 1] : -1;
s3 = ((l < n) && (sn[l + 1] == sn[l])) ? S1[l + 1] : -1;
temp *= exp_E_MLstem(tt, s5, s3, pf_params) * scale[2];
probs[kl] += temp;
if (probs[kl]>Qmax) {
Qmax = probs[kl];
if (Qmax>max_real/10.)
vrna_message_warning("P close to overflow: %d %d %g %g",
i, j, probs[kl], qb[kl]);
}
if (probs[kl]>=max_real) {
ov++;
probs[kl]=FLT_MAX;
}
}
} /* end for (k=..) multloop*/
else /* set prm_l to 0 to get prm_l1 to be 0 */
for (i=0; i<=n; i++) prm_l[i]=0;
tmp = prm_l1; prm_l1=prm_l; prm_l=tmp;
/*computation of .(..(...)..&..). type features?*/
if (cp<=0) continue; /* no .(..(...)..&..). type features*/
if ((l==n)||(l<=2)) continue; /* no .(..(...)..&..). type features*/
/*new version with O(n^3)??*/
if (l>cp) {
if (l<n) {
int t,kt;
for (t=n; t>l; t--) {
for (k=1; k<cp; k++) {
int samestrand;
kt = my_iindx[k]-t;
samestrand = (sn[k + 1] == sn[k]) ? 1 : 0;
type = rtype[(unsigned char)ptype[jindx[t] + k]];
if(type == 0)
type = 7;
temp = probs[kt]
* exp_E_ExtLoop(type, S1[t-1], samestrand ? S1[k+1] : -1, pf_params)
* scale[2];
if (l+1<t)
temp *= q[my_iindx[l+1]-(t-1)];
if (samestrand)
temp *= q[my_iindx[k+1]-(cp-1)];
Qrout[l] += temp;
}
}
}
for (k=l-1; k>=cp; k--) {
if (qb[my_iindx[k]-l]) {
kl = my_iindx[k]-l;
type = ptype[jindx[l] + k];
temp = Qrout[l];
if(type == 0)
type = 7;
temp *= exp_E_ExtLoop(type, (k>cp) ? S1[k-1] : -1, (l < n) ? S1[l+1] : -1, pf_params);
if (k>cp)
temp *= q[my_iindx[cp]-(k-1)];
probs[kl] += temp;
}
}
}
else if (l==cp ) {
int t, sk,s;
for (t=2; t<cp;t++) {
for (s=1; s<t; s++) {
for (k=cp; k<=n; k++) {
sk=my_iindx[s]-k;
if (qb[sk]) {
int samestrand;
samestrand = (sn[k] == sn[k - 1]) ? 1 : 0;
type = rtype[(unsigned char)ptype[jindx[k] + s]];
if(type == 0)
type = 7;
temp = probs[sk]
* exp_E_ExtLoop(type, samestrand ? S1[k - 1] : -1, S1[s + 1], pf_params)
* scale[2];
if (s+1<t)
temp *= q[my_iindx[s+1]-(t-1)];
if (samestrand)
temp *= q[my_iindx[cp]-(k-1)];
Qlout[t] += temp;
}
}
}
}
}
else if (l<cp) {
for (k=1; k<l; k++) {
if (qb[my_iindx[k]-l]) {
type = ptype[jindx[l] + k];
temp = Qlout[k];
if(type == 0)
type = 7;
temp *= exp_E_ExtLoop(type, (k>1) ? S1[k-1] : -1, (l<(cp-1)) ? S1[l+1] : -1, pf_params);
if (l+1<cp)
temp *= q[my_iindx[l+1]-(cp-1)];
probs[my_iindx[k]-l] += temp;
}
}
}
} /* end for (l=..) */
free(Qlout);
free(Qrout);
for (i=1; i<=n; i++)
for (j=i+turn+1; j<=n; j++) {
ij = my_iindx[i]-j;
probs[ij] *= qb[ij];
}
if (structure!=NULL){
char *s = vrna_db_from_probs(probs, (unsigned int)n);
memcpy(structure, s, n);
structure[n] = '\0';
free(s);
}
/* clean up */
free(prm_l);
free(prm_l1);
free(prml);
} /* end if (do_backtrack)*/
if(ov > 0)
vrna_message_warning("%d overflows occurred while backtracking;\n"
"you might try a smaller pf_scale than %g\n",
ov, pf_params->pf_scale);
}
PUBLIC void
vrna_pf_dimer_probs(double FAB,
double FA,
double FB,
vrna_plist_t *prAB,
const vrna_plist_t *prA,
const vrna_plist_t *prB,
int Alength,
const vrna_exp_param_t *exp_params) {
/*computes binding probabilities and dimer free energies*/
int i, j;
double pAB;
double mykT;
const vrna_plist_t *lp2;
vrna_plist_t *lp1;
int offset;
mykT = exp_params->kT/1000.;
/* pair probabilities in pr are relative to the null model (without DuplexInit) */
/*Compute probabilities pAB, pAA, pBB*/
pAB = 1. - exp((1/mykT)*(FAB-FA-FB));
/* compute pair probabilities given that it is a dimer */
/* AB dimer */
offset = 0;
lp2 = prA;
if (pAB>0)
for (lp1=prAB; lp1->j>0; lp1++) {
float pp=0;
i = lp1->i;
j = lp1->j;
while (offset+lp2->i < i && lp2->i>0) lp2++;
if (offset+lp2->i == i)
while ((offset+lp2->j) < j && (lp2->j>0)) lp2++;
if (lp2->j == 0) {lp2=prB; offset=Alength;}/* jump to next list */
if ((offset+lp2->i==i) && (offset+lp2->j ==j)) {
pp = lp2->p;
lp2++;
}
lp1->p=(lp1->p-(1-pAB)*pp)/pAB;
if(lp1->p < 0.){
vrna_message_warning("vrna_co_pf_probs: numeric instability detected, probability below zero!");
lp1->p = 0.;
}
}
return;
}
PRIVATE double *
Newton_Conc(double KAB,
double KAA,
double KBB,
double concA,
double concB,
double* ConcVec){
double TOL, EPS, xn, yn, det, cA, cB;
int i;
i = 0;
/*Newton iteration for computing concentrations*/
cA = concA;
cB = concB;
TOL = 1e-6; /*Tolerance for convergence*/
ConcVec = (double*)vrna_alloc(5*sizeof(double)); /* holds concentrations */
do {
/* det = (4.0 * KAA * cA + KAB *cB + 1.0) * (4.0 * KBB * cB + KAB *cA + 1.0) - (KAB *cB) * (KAB *cA); */
det = 1 + 16. *KAA*KBB*cA*cB + KAB*(cA+cB) + 4.*KAA*cA + 4.*KBB*cB + 4.*KAB*(KBB*cB*cB + KAA*cA*cA);
/* xn = ( (2.0 * KBB * cB*cB + KAB *cA *cB + cB - concB) * (KAB *cA) -
(2.0 * KAA * cA*cA + KAB *cA *cB + cA - concA) * (4.0 * KBB * cB + KAB *cA + 1.0) ) /det; */
xn = ( (2.0 * KBB * cB*cB + cB - concB) * (KAB *cA) - KAB*cA*cB*(4. * KBB*cB + 1.) -
(2.0 * KAA * cA*cA + cA - concA) * (4.0 * KBB * cB + KAB *cA + 1.0) ) /det;
/* yn = ( (2.0 * KAA * cA*cA + KAB *cA *cB + cA - concA) * (KAB *cB) -
(2.0 * KBB * cB*cB + KAB *cA *cB + cB - concB) * (4.0 * KAA * cA + KAB *cB + 1.0) ) /det; */
yn = ( (2.0 * KAA * cA*cA + cA - concA) * (KAB *cB) - KAB*cA*cB*(4. * KAA*cA + 1.) -
(2.0 * KBB * cB*cB + cB - concB) * (4.0 * KAA * cA + KAB *cB + 1.0) ) /det;
EPS = fabs(xn/cA) + fabs(yn/cB);
cA += xn;
cB += yn;
i++;
if (i>10000) {
vrna_message_warning("Newton did not converge after %d steps!!",i);
break;
}
} while(EPS>TOL);
ConcVec[0] = cA*cB*KAB ;/*AB concentration*/
ConcVec[1] = cA*cA*KAA ;/*AA concentration*/
ConcVec[2] = cB*cB*KBB ;/*BB concentration*/
ConcVec[3] = cA; /* A concentration*/
ConcVec[4] = cB; /* B concentration*/
return ConcVec;
}
PUBLIC vrna_dimer_conc_t *
vrna_pf_dimer_concentrations(double FcAB,
double FcAA,
double FcBB,
double FEA,
double FEB,
const double *startconc,
const vrna_exp_param_t *exp_params){
/*takes an array of start concentrations, computes equilibrium concentrations of dimers, monomers, returns array of concentrations in strucutre vrna_dimer_conc_t*/
double *ConcVec;
int i;
vrna_dimer_conc_t *Concentration;
double KAA, KAB, KBB, kT;
kT = exp_params->kT/1000.;
Concentration = (vrna_dimer_conc_t *)vrna_alloc(20*sizeof(vrna_dimer_conc_t));
/* Compute equilibrium constants */
/* again note the input free energies are not from the null model (without DuplexInit) */
KAA = exp(( 2.0 * FEA - FcAA)/kT);
KBB = exp(( 2.0 * FEB - FcBB)/kT);
KAB = exp(( FEA + FEB - FcAB)/kT);
/* printf("Kaa..%g %g %g\n", KAA, KBB, KAB); */
for (i=0; ((startconc[i]!=0)||(startconc[i+1]!=0));i+=2) {
ConcVec = Newton_Conc(KAB, KAA, KBB, startconc[i], startconc[i+1], ConcVec);
Concentration[i/2].A0 = startconc[i];
Concentration[i/2].B0 = startconc[i+1];
Concentration[i/2].ABc = ConcVec[0];
Concentration[i/2].AAc = ConcVec[1];
Concentration[i/2].BBc = ConcVec[2];
Concentration[i/2].Ac = ConcVec[3];
Concentration[i/2].Bc = ConcVec[4];
if (!(((i+2)/2)%20)) {
Concentration = (vrna_dimer_conc_t *)vrna_realloc(Concentration,((i+2)/2+20)*sizeof(vrna_dimer_conc_t));
}
free(ConcVec);
}
return Concentration;
}
#if 0
/*
stochastic backtracking in pf_fold arrays
returns random structure S with Boltzman probabilty
p(S) = exp(-E(S)/kT)/Z
*/
PRIVATE void
backtrack_qm1(vrna_fold_compound_t *vc,
int i,
int j,
char *pstruc){
/* i is paired to l, i<l<j; backtrack in qm1 to find l */
int ii, l, type, *jindx, *my_iindx, *rtype, turn;
double qt, r;
FLT_OR_DBL *qm, *qm1, *qb, *expMLbase;
short *S1;
char *ptype;
vrna_md_t *md;
vrna_exp_param_t *pf_params;
vrna_mx_pf_t *matrices;
pf_params = vc->exp_params;
md = &(pf_params->model_details);
S1 = vc->sequence_encoding;
ptype = vc->ptype;
rtype = &(md->rtype[0]);
turn = md->min_loop_size;
matrices = vc->exp_matrices;
qb = matrices->qb;
qm = matrices->qm;
qm1 = matrices->qm1;
expMLbase = matrices->expMLbase;
jindx = vc->jindx;
my_iindx = vc->iindx;
r = vrna_urn() * qm1[jindx[j]+i];
ii = my_iindx[i];
for (qt=0., l=i+turn+1; l<=j; l++) {
type = ptype[jindx[l] + i];
if (type)
qt += qb[ii-l]*exp_E_MLstem(type, S1[i-1], S1[l+1], pf_params) * expMLbase[j-l];
if (qt>=r) break;
}
if (l>j) vrna_message_error("backtrack failed in qm1");
backtrack(vc, i,l, pstruc);
}
PRIVATE void
backtrack(vrna_fold_compound_t *vc,
int i,
int j,
char *pstruc){
int *jindx, *my_iindx, *rtype, turn;
FLT_OR_DBL *qm, *qm1, *qb, *expMLbase, *scale;
vrna_exp_param_t *pf_params;
vrna_mx_pf_t *matrices;
short *S1;
char *ptype, *sequence;
int noGUclosure;
vrna_md_t *md;
sequence = vc->sequence;
pf_params = vc->exp_params;
md = &(pf_params->model_details);
S1 = vc->sequence_encoding;
ptype = vc->ptype;
rtype = &(md->rtype[0]);
turn = md->min_loop_size;
matrices = vc->exp_matrices;
qb = matrices->qb;
qm = matrices->qm;
qm1 = matrices->qm1;
expMLbase = matrices->expMLbase;
scale = matrices->scale;
jindx = vc->jindx;
my_iindx = vc->iindx;
noGUclosure = pf_params->model_details.noGUclosure;
do {
double r, qbt1;
int k, l, type, u, u1;
pstruc[i-1] = '('; pstruc[j-1] = ')';
r = vrna_urn() * qb[my_iindx[i]-j];
type = ptype[jindx[j] + i];
u = j - i - 1;
/*hairpin contribution*/
if (((type==3)||(type==4))&&noGUclosure) qbt1 = 0;
else
qbt1 = exp_E_Hairpin(u, type, S1[i+1], S1[j-1], sequence+i-1, pf_params)*scale[u+2];
if (qbt1>r) return; /* found the hairpin we're done */
for (k=i+1; k<=MIN2(i+MAXLOOP+1,j-turn-2); k++) {
u1 = k-i-1;
for (l=MAX2(k+turn+1,j-1-MAXLOOP+u1); l<j; l++) {
int type_2;
type_2 = ptype[jindx[l] + k];
if (type_2) {
type_2 = rtype[type_2];
qbt1 += qb[my_iindx[k]-l] *
exp_E_IntLoop(u1, j-l-1, type, type_2,
S1[i+1], S1[j-1], S1[k-1], S1[l+1], pf_params)*scale[u1+j-l+1];
}
if (qbt1 > r) break;
}
if (qbt1 > r) break;
}
if (l<j) {
i=k; j=l;
}
else break;
} while (1);
/* backtrack in multi-loop */
{
double r, qt;
int k, ii, jj;
i++; j--;
/* find the first split index */
ii = my_iindx[i]; /* ii-j=[i,j] */
jj = jindx[j]; /* jj+i=[j,i] */
for (qt=0., k=i+1; k<j; k++) qt += qm[ii-(k-1)]*qm1[jj+k];
r = vrna_urn() * qt;
for (qt=0., k=i+1; k<j; k++) {
qt += qm[ii-(k-1)]*qm1[jj+k];
if (qt>=r) break;
}
if (k>=j) vrna_message_error("backtrack failed, can't find split index ");
backtrack_qm1(vc, k, j, pstruc);
j = k-1;
while (j>i) {
/* now backtrack [i ... j] in qm[] */
jj = jindx[j];
ii = my_iindx[i];
r = vrna_urn() * qm[ii - j];
qt = qm1[jj+i]; k=i;
if (qt<r)
for (k=i+1; k<=j; k++) {
qt += (qm[ii-(k-1)]+expMLbase[k-i])*qm1[jj+k];
if (qt >= r) break;
}
if (k>j) vrna_message_error("backtrack failed in qm");
backtrack_qm1(vc, k,j, pstruc);
if (k<i+turn) break; /* no more pairs */
r = vrna_urn() * (qm[ii-(k-1)] + expMLbase[k-i]);
if (expMLbase[k-i] >= r) break; /* no more pairs */
j = k-1;
}
}
}
#endif
/*###########################################*/
/*# deprecated functions below #*/
/*###########################################*/
PUBLIC vrna_dimer_pf_t
co_pf_fold(char *sequence, char *structure){
return wrap_co_pf_fold(sequence, structure, NULL, do_backtrack, fold_constrained);
}
PUBLIC vrna_dimer_pf_t
co_pf_fold_par( char *sequence,
char *structure,
vrna_exp_param_t *parameters,
int calculate_bppm,
int is_constrained){
return wrap_co_pf_fold(sequence, structure, parameters, calculate_bppm, is_constrained);
}
PUBLIC vrna_plist_t *
get_plist(vrna_plist_t *pl,
int length,
double cut_off){
int i, j,n, count, *my_iindx;
my_iindx = backward_compat_compound->iindx;
/*get pair probibilities out of pr array*/
count=0;
n=2;
for (i=1; i<length; i++) {
for (j=i+1; j<=length; j++) {
if (pr[my_iindx[i]-j]<cut_off) continue;
if (count==n*length-1) {
n*=2;
pl=(vrna_plist_t *)vrna_realloc(pl,n*length*sizeof(vrna_plist_t));
}
pl[count].i=i;
pl[count].j=j;
pl[count++].p=pr[my_iindx[i]-j];
/* printf("gpl: %2d %2d %.9f\n",i,j,pr[my_iindx[i]-j]);*/
}
}
pl[count].i=0;
pl[count].j=0; /*->??*/
pl[count++].p=0.;
pl=(vrna_plist_t *)vrna_realloc(pl,(count)*sizeof(vrna_plist_t));
return pl;
}
PUBLIC void
compute_probabilities(double FAB,
double FA,
double FB,
vrna_plist_t *prAB,
vrna_plist_t *prA,
vrna_plist_t *prB,
int Alength) {
if(backward_compat_compound && backward_compat){
vrna_pf_dimer_probs(FAB, FA, FB, prAB, (const vrna_plist_t *)prA, (const vrna_plist_t *)prB, Alength, (const vrna_exp_param_t *)backward_compat_compound->exp_params);
}
}
PUBLIC vrna_dimer_conc_t *
get_concentrations( double FcAB,
double FcAA,
double FcBB,
double FEA,
double FEB,
double *startconc){
return vrna_pf_dimer_concentrations(FcAB, FcAA, FcBB, FEA, FEB, (const double *)startconc, (const vrna_exp_param_t *)backward_compat_compound->exp_params);
}
PUBLIC void
init_co_pf_fold(int length){
/* DO NOTHING */
}
PUBLIC void
free_co_pf_arrays(void){
if(backward_compat_compound && backward_compat){
vrna_fold_compound_free(backward_compat_compound);
backward_compat_compound = NULL;
backward_compat = 0;
}
}
PUBLIC FLT_OR_DBL *
export_co_bppm(void){
if(backward_compat_compound)
return backward_compat_compound->exp_matrices->probs;
else
return NULL;
}
/*----------------------------------------------------------------------*/
PUBLIC void
update_co_pf_params(int length){
if(backward_compat_compound && backward_compat){
vrna_md_t md;
set_model_details(&md);
vrna_exp_params_reset(backward_compat_compound, &md);
/* compatibility with RNAup, may be removed sometime */
pf_scale = backward_compat_compound->exp_params->pf_scale;
}
}
PUBLIC void
update_co_pf_params_par(int length,
vrna_exp_param_t *parameters){
if(backward_compat_compound && backward_compat){
vrna_md_t md;
if(parameters){
vrna_exp_params_subst(backward_compat_compound, parameters);
} else {
set_model_details(&md);
vrna_exp_params_reset(backward_compat_compound, &md);
}
/* compatibility with RNAup, may be removed sometime */
pf_scale = backward_compat_compound->exp_params->pf_scale;
}
}