/* -*-C-*- */
/* this file contains code for folding circular RNAs */
/* it's #include'd into fold.c */
PRIVATE void
fill_arrays_circ( vrna_fold_compound_t *vc,
sect bt_stack[],
int *bt){
/* variant of fold() for circular RNAs */
/* auxiliarry arrays:
fM2 = multiloop region with exactly two stems, extending to 3' end
for stupid dangles=1 case we also need:
fM_d3 = multiloop region with >= 2 stems, starting at pos 2
(a pair (k,n) will form 3' dangle with pos 1)
fM_d5 = multiloop region with >= 2 stems, extending to pos n-1
(a pair (1,k) will form a 5' dangle with pos n)
*/
int Hi, Hj, Ii, Ij, Ip, Iq, ip, iq, Mi;
int *fM_d3, *fM_d5, Md3i, Md5i, FcMd3, FcMd5;
int FcH, FcI, FcM, Fc, *fM2;
int i,j, ij, u, length, new_c, fm, type;
int *my_c, *my_fML, *indx, dangle_model, turn;
vrna_param_t *P;
char *ptype, *hard_constraints;
short *S1;
vrna_hc_t *hc;
vrna_sc_t *sc;
P = vc->params;
ptype = vc->ptype;
indx = vc->jindx;
S1 = vc->sequence_encoding;
dangle_model = P->model_details.dangles;
turn = P->model_details.min_loop_size;
hc = vc->hc;
sc = vc->sc;
hard_constraints = hc->matrix;
my_c = vc->matrices->c;
my_fML = vc->matrices->fML;
fM2 = vc->matrices->fM2;
length = vc->length;
FcH = FcI = FcM = FcMd3 = FcMd5 = INF;
if(hc->up_ext[1] >= length){
Fc = 0;
if(sc){
if(sc->energy_up)
Fc += sc->energy_up[1][length];
}
} else {
Fc = INF;
}
for(i = 1; i < length; i++)
for(j = i + turn + 1; j <= length; j++){
u = length-j + i-1;
if (u<TURN) continue;
ij = indx[j] + i;
if (!hard_constraints[ij])
continue;
/* exterior hairpin case */
new_c = vrna_E_hp_loop(vc, j, i)
+ my_c[ij];
if (new_c<FcH) {
FcH = new_c;
Hi = i;
Hj = j;
}
/* exterior interior loop case */
ip = iq = 0;
new_c = vrna_E_ext_int_loop(vc, i, j, &ip, &iq)
+ my_c[ij];
if(ip != 0){
if(new_c < FcI){
FcI = new_c;
Ii = i;
Ij = j;
Ip = ip;
Iq = iq;
}
}
} /* end of i,j loop */
Fc = MIN2(Fc, FcH);
Fc = MIN2(Fc, FcI);
/* compute the fM2 array (multi loops with exactly 2 helices) */
/* to get a unique ML decomposition, just use fM1 instead of fML
below. However, that will not work with dangle_model==1 */
for (i=1; i<length-TURN; i++) {
fM2[i] = INF;
for (u=i+TURN; u<length-TURN; u++)
fM2[i] = MIN2(fM2[i], my_fML[indx[u]+i] + my_fML[indx[length]+u+1]);
}
for (i=TURN+1; i<length-2*TURN; i++) {
fm = my_fML[indx[i]+1]+fM2[i+1]+P->MLclosing;
if (fm<FcM) {
FcM=fm; Mi=i;
}
}
Fc = MIN2(Fc, FcM);
if ((dangle_model == 1) || (dangle_model == 3)){
fM_d3 = (int *) vrna_alloc(sizeof(int)*(length+2));
fM_d5 = (int *) vrna_alloc(sizeof(int)*(length+2));
for (i=TURN+1; i<length-TURN; i++) {
fM_d3[i] = INF;
for (u=2+TURN; u<i-TURN; u++)
fM_d3[i] = MIN2(fM_d3[i], my_fML[indx[u]+2] + my_fML[indx[i]+u+1]);
}
for (i=2*TURN+1; i<length-TURN; i++) {
type = ptype[indx[length]+i+1];
if(hard_constraints[indx[length] + i + 1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
if(hc->up_ml[1]){
fm = fM_d3[i]+my_c[indx[length]+i+1]+E_MLstem(type, -1, S1[1], P) + P->MLclosing;
if (fm<FcMd3) {
FcMd3=fm; Md3i=i;
}
if(hc->up_ml[i]){
fm = fM_d3[i-1]+my_c[indx[length]+i+1]+E_MLstem(type, S1[i], S1[1], P) + P->MLclosing;
if (fm<FcMd3) {
FcMd3=fm; Md3i=-i;
}
}
}
}
}
for (i=TURN+1; i<length-TURN; i++) {
fM_d5[i] = INF;
for (u=i+TURN; u<length-TURN; u++)
fM_d5[i] = MIN2(fM_d5[i], my_fML[indx[u]+i] + my_fML[indx[length-1]+u+1]);
}
for (i=TURN+1; i<length-2*TURN; i++) {
if(hard_constraints[indx[i]+1] & VRNA_CONSTRAINT_CONTEXT_MB_LOOP){
type = ptype[indx[i]+1];
if(hc->up_ml[length]){
fm = E_MLstem(type, S1[length], -1, P) + my_c[indx[i]+1] + fM_d5[i+1] + P->MLclosing;
if (fm<FcMd5) {
FcMd5=fm; Md5i=i;
}
if(hc->up_ml[i+1]){
fm = E_MLstem(type, S1[length], S1[i+1], P) + my_c[indx[i]+1] + fM_d5[i+2] + P->MLclosing;
if (fm<FcMd5) {
FcMd5=fm; Md5i=-i;
}
}
}
}
}
if (FcMd5<MIN2(Fc,FcMd3)) {
/* looks like we have to do this ... */
bt_stack[++(*bt)].i = 1;
bt_stack[(*bt)].j = (Md5i>0)?Md5i:-Md5i;
bt_stack[(*bt)].ml = 2;
i = (Md5i>0)?Md5i+1 : -Md5i+2; /* let's backtrack fm_d5[Md5i+1] */
for (u=i+TURN; u<length-TURN; u++)
if (fM_d5[i] == my_fML[indx[u]+i] + my_fML[indx[length-1]+u+1]) {
bt_stack[++(*bt)].i = i;
bt_stack[(*bt)].j = u;
bt_stack[(*bt)].ml = 1;
bt_stack[++(*bt)].i =u+1;
bt_stack[(*bt)].j = length-1;
bt_stack[(*bt)].ml = 1;
break;
}
Fc = FcMd5;
} else if (FcMd3<Fc) {
/* here we go again... */
bt_stack[++(*bt)].i = (Md3i>0)?Md3i+1:-Md3i+1;
bt_stack[(*bt)].j = length;
bt_stack[(*bt)].ml = 2;
i = (Md3i>0)? Md3i : -Md3i-1; /* let's backtrack fm_d3[Md3i] */
for (u=2+TURN; u<i-TURN; u++)
if (fM_d3[i] == my_fML[indx[u]+2] + my_fML[indx[i]+u+1]) {
bt_stack[++(*bt)].i = 2;
bt_stack[(*bt)].j = u;
bt_stack[(*bt)].ml = 1;
bt_stack[++(*bt)].i =u+1;
bt_stack[(*bt)].j = i;
bt_stack[(*bt)].ml = 1;
break;
}
Fc = FcMd3;
}
free(fM_d3);
free(fM_d5);
}
else if(Fc < 0){
if (FcH==Fc) {
bt_stack[++(*bt)].i = Hi;
bt_stack[(*bt)].j = Hj;
bt_stack[(*bt)].ml = 2;
}
else if (FcI==Fc) {
bt_stack[++(*bt)].i = Ii;
bt_stack[(*bt)].j = Ij;
bt_stack[(*bt)].ml = 2;
bt_stack[++(*bt)].i = Ip;
bt_stack[(*bt)].j = Iq;
bt_stack[(*bt)].ml = 2;
}
else if (FcM==Fc) { /* grumpf we found a Multiloop */
/* backtrack in fM2 */
fm = fM2[Mi+1];
for (u=Mi+TURN+1; u<length-TURN; u++)
if (fm == my_fML[indx[u]+Mi+1] + my_fML[indx[length]+u+1]) {
bt_stack[++(*bt)].i=Mi+1;
bt_stack[(*bt)].j=u;
bt_stack[(*bt)].ml = 1;
bt_stack[++(*bt)].i=u+1;
bt_stack[(*bt)].j=length;
bt_stack[(*bt)].ml = 1;
break;
}
bt_stack[++(*bt)].i = 1;
bt_stack[(*bt)].j = Mi;
bt_stack[(*bt)].ml = 1;
}
} else { /* unstructured */
bt_stack[++(*bt)].i = 1;
bt_stack[(*bt)].j = 1;
bt_stack[(*bt)].ml = 0;
}
vc->matrices->FcH = FcH;
vc->matrices->FcI = FcI;
vc->matrices->FcM = FcM;
vc->matrices->Fc = Fc;
}