/* minimum free energy RNA secondary structure with basepair distance d_1 to reference structure 1 and distance d_2 to reference structure 2 */ #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include #include #include #include "ViennaRNA/utils.h" #include "ViennaRNA/energy_par.h" #include "ViennaRNA/fold_vars.h" #include "ViennaRNA/fold.h" #include "ViennaRNA/loop_energies.h" #include "ViennaRNA/params.h" #ifdef _OPENMP #include #endif #include "ViennaRNA/2Dfold.h" /* ################################# # GLOBAL VARIABLES # ################################# */ int compute_2Dfold_F3 = 0; /* ################################# # PRIVATE VARIABLES # ################################# */ /* ################################# # PRIVATE FUNCTION DECLARATIONS # ################################# */ PRIVATE void mfe_linear(vrna_fold_compound_t *vc); PRIVATE void mfe_circ(vrna_fold_compound_t *vc); PUBLIC void update_TwoDfold_params(TwoDfold_vars *vars); PRIVATE void backtrack_f5(unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void backtrack_c(unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void backtrack_m(unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void backtrack_m1(unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void backtrack_fc(int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void backtrack_m2(unsigned int i, int k, int l, char *structure, vrna_fold_compound_t *vc); PRIVATE void adjustArrayBoundaries(int ***array, int *k_min, int *k_max, int **l_min, int **l_max, int k_min_real, int k_max_real, int *l_min_real, int *l_max_real); INLINE PRIVATE void preparePosteriorBoundaries(int size, int shift, int *min_k, int *max_k, int **min_l, int **max_l); INLINE PRIVATE void updatePosteriorBoundaries(int d1, int d2, int *min_k, int *max_k, int **min_l, int **max_l); INLINE PRIVATE void prepareBoundaries(int min_k_pre, int max_k_pre, int min_l_pre, int max_l_pre, int bpdist, int *min_k, int *max_k, int **min_l, int **max_l); INLINE PRIVATE void prepareArray(int ***array, int min_k, int max_k, int *min_l, int *max_l); INLINE PRIVATE void prepareArray2(unsigned long ***array, int min_k, int max_k, int *min_l, int *max_l); /* ################################# # BEGIN OF FUNCTION DEFINITIONS # ################################# */ #if 0 PRIVATE void initialize_TwoDfold_vars(TwoDfold_vars *vars){ update_TwoDfold_params(vars); /* this call updates the params in the ViennaRNA fold.o which is a global, so be careful * whith calling it parallel... need a workarround or fix of ViennaRNA fold stuff */ update_fold_params(); } PUBLIC TwoDfold_solution **TwoDfold(TwoDfold_vars *vars, int distance1, int distance2){ unsigned int i, d1, d2; unsigned int maxD1; unsigned int maxD2; unsigned int length; TwoDfold_solution **output; initialize_TwoDfold_vars(vars); if(fabs(vars->P->temperature - temperature)>1e-6) update_TwoDfold_params(vars); vars->S = encode_sequence(vars->sequence, 0); vars->S1 = encode_sequence(vars->sequence, 1); make_ptypes(vars); maxD1 = vars->maxD1; maxD2 = vars->maxD2; if(distance1 >= 0){ if((unsigned int)distance1 > maxD1) fprintf(stderr, "limiting maximum basepair distance 1 to %u\n", maxD1); else maxD1 = (unsigned int)distance1; } if(distance2 >= 0){ if((unsigned int)distance2 > maxD2) fprintf(stderr, "limiting maximum basepair distance 2 to %u\n", maxD2); else maxD2 = (unsigned int)distance2; } vars->maxD1 = maxD1; vars->maxD2 = maxD2; output = (TwoDfold_solution **)vrna_alloc((vars->maxD1+1) * sizeof(TwoDfold_solution *)); mfe_linear(vars); if(vars->circ) mfe_circ(vars); length = vars->seq_length; for(d1=0; d1<=maxD1;d1++){ output[d1] = (TwoDfold_solution *)vrna_alloc((vars->maxD2+1)*sizeof(TwoDfold_solution)); #ifdef _OPENMP #pragma omp parallel for private(d2) #endif for(d2=0; d2<=maxD2;d2++){ output[d1][d2].en = (float)INF/(float)100.; output[d1][d2].s = NULL; } if( (d1 >= ((vars->circ) ? vars->k_min_values_fc : vars->k_min_values_f[length])) && (d1 <= ((vars->circ) ? vars->k_max_values_fc : vars->k_max_values_f[length]))){ #ifdef _OPENMP #pragma omp parallel for private(d2, i) #endif for( d2 = ((vars->circ) ? vars->l_min_values_fc[d1] : vars->l_min_values_f[length][d1]); d2 <= ((vars->circ) ? vars->l_max_values_fc[d1] : vars->l_max_values_f[length][d1]); d2 += 2){ output[d1][d2].en = (float)((vars->circ) ? vars->E_Fc[d1][d2/2] : vars->E_F5[length][d1][d2/2])/(float)100.; if(vars->do_backtrack && (output[d1][d2].en != (float)INF/(float)100.)){ char *mfe_structure = (char *)vrna_alloc(length+1); for(i=0;icirc) ? backtrack_fc(d1, d2, mfe_structure, vars) : backtrack_f5(length, d1, d2, mfe_structure, vars); output[d1][d2].s = mfe_structure; } } } } return output; } #endif PUBLIC vrna_sol_TwoD_t * vrna_mfe_TwoD(vrna_fold_compound_t *vars, int distance1, int distance2){ unsigned int i, d1, d2; unsigned int maxD1; unsigned int maxD2; unsigned int length; unsigned int counter = 0; int en = 0; vrna_sol_TwoD_t *output; vrna_md_t *md; vrna_mx_mfe_t *matrices; maxD1 = vars->maxD1; maxD2 = vars->maxD2; matrices = vars->matrices; md = &(vars->params->model_details); if(distance1 >= 0){ if((unsigned int)distance1 > maxD1) vrna_message_warning("vrna_mfe_TwoD@2Dfold.c: limiting maximum basepair distance 1 to %u\n", maxD1); else maxD1 = (unsigned int)distance1; } if(distance2 >= 0){ if((unsigned int)distance2 > maxD2) vrna_message_warning("vrna_mfe_TwoD@2Dfold.c: limiting maximum basepair distance 2 to %u\n", maxD2); else maxD2 = (unsigned int)distance2; } vars->maxD1 = maxD1; vars->maxD2 = maxD2; output = (vrna_sol_TwoD_t *)vrna_alloc((((vars->maxD1+1)*(vars->maxD2+2))/2 + 2) * sizeof(vrna_sol_TwoD_t)); mfe_linear(vars); if(md->circ) mfe_circ(vars); length = vars->length; for(d1=0; d1<=maxD1;d1++){ if((d1 >= ((md->circ) ? matrices->k_min_Fc : matrices->k_min_F5[length])) && (d1 <= ((md->circ) ? matrices->k_max_Fc : matrices->k_max_F5[length]))){ for(d2 = ((md->circ) ? matrices->l_min_Fc[d1] : matrices->l_min_F5[length][d1]); d2 <= ((md->circ) ? matrices->l_max_Fc[d1] : matrices->l_max_F5[length][d1]); d2 += 2){ en = ((md->circ) ? matrices->E_Fc[d1][d2/2] : matrices->E_F5[length][d1][d2/2]); if(en == INF) continue; output[counter].k = d1; output[counter].l = d2; output[counter].en = (float)en/(float)100.; if(md->backtrack){ char *mfe_structure = (char *)vrna_alloc(length+1); for(i=0;icirc) ? backtrack_fc((int)d1, (int)d2, mfe_structure, vars) : backtrack_f5(length, (int)d1, (int)d2, mfe_structure, vars); output[counter].s = mfe_structure; } else output[counter].s = NULL; counter++; } } } /* store entry for remaining partition if it exists */ en = ((md->circ) ? matrices->E_Fc_rem : matrices->E_F5_rem[length]); if(en != INF){ output[counter].k = -1; output[counter].l = -1; output[counter].en = (float)en/(float)100.; if(md->backtrack){ char *mfe_structure = (char *)vrna_alloc(length+1); for(i=0;icirc) ? backtrack_fc(-1, -1, mfe_structure, vars) : backtrack_f5(length, -1, -1, mfe_structure, vars); output[counter].s = mfe_structure; } else output[counter].s = NULL; counter++; } /* insert end-marker entry */ output[counter].k = output[counter].l = INF; counter++; /* resize to actual dataset amount */ output = (vrna_sol_TwoD_t *)vrna_realloc(output, sizeof(vrna_sol_TwoD_t) * counter); return output; } PUBLIC char * vrna_backtrack5_TwoD( vrna_fold_compound_t *vc, int k, int l, unsigned int j){ unsigned int i; char *mfe_structure = (char *)vrna_alloc(j+1); if(j < TURN + 2) return NULL; for(i=0; i < j; i++) mfe_structure[i] = '.'; mfe_structure[i] = '\0'; backtrack_f5(j, k, l, mfe_structure, vc); return mfe_structure; } PRIVATE void mfe_linear(vrna_fold_compound_t *vc){ unsigned int d, i, j, ij, maxD1, maxD2, seq_length, dia, dib, dja, djb, *referenceBPs1, *referenceBPs2, *mm1, *mm2, *bpdist; int cnt1, cnt2, cnt3, cnt4, d1, d2, energy, dangles, temp2, type, additional_en, *my_iindx, *jindx, circ, *rtype; short *S1, *reference_pt1, *reference_pt2; char *sequence, *ptype; vrna_param_t *P; vrna_mx_mfe_t *matrices; vrna_md_t *md; /* dereferenciate things we often need */ P = vc->params; md = &(P->model_details); matrices = vc->matrices; sequence = vc->sequence; seq_length = vc->length; maxD1 = vc->maxD1; maxD2 = vc->maxD2; S1 = vc->sequence_encoding; ptype = vc->ptype; rtype = &(md->rtype[0]); reference_pt1 = vc->reference_pt1; reference_pt2 = vc->reference_pt2; my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; mm1 = vc->mm1; mm2 = vc->mm2; bpdist = vc->bpdist; dangles = md->dangles; circ = md->circ; for (d = TURN+2; d <= seq_length; d++) { /* i,j in [1..length] */ #ifdef _OPENMP #pragma omp parallel for private(additional_en, j, energy, temp2, i, ij, dia,dib,dja,djb,cnt1,cnt2,cnt3,cnt4, d1, d2) #endif for (j = d; j <= seq_length; j++) { unsigned int p, q, pq, u, maxp, dij; int type_2, type, tt, no_close, base_d1, base_d2; i = j-d+1; dij = j - i - 1; ij = my_iindx[i]-j; type = ptype[jindx[j] + i]; no_close = (((type==3)||(type==4))&&no_closingGU); if (type) { /* we have a pair */ /* increase or decrease distance-to-reference value depending whether (i,j) is included in * reference or has to be introduced */ base_d1 = ((unsigned int)reference_pt1[i] != j) ? 1 : -1; base_d2 = ((unsigned int)reference_pt2[i] != j) ? 1 : -1; /* HAIRPIN STRUCTURES */ /* get distance to reference if closing the hairpin * d = dbp(T_{i,j}, {i,j}) */ d1 = base_d1 + referenceBPs1[ij]; d2 = base_d2 + referenceBPs2[ij]; int min_k, max_k, min_l, max_l; int real_min_k, real_max_k, *min_l_real, *max_l_real; min_l = min_k = 0; max_k = mm1[ij] + referenceBPs1[ij]; max_l = mm2[ij] + referenceBPs2[ij]; prepareBoundaries(min_k, max_k, min_l, max_l, bpdist[ij], &matrices->k_min_C[ij], &matrices->k_max_C[ij], &matrices->l_min_C[ij], &matrices->l_max_C[ij] ); preparePosteriorBoundaries( matrices->k_max_C[ij] - matrices->k_min_C[ij] + 1, matrices->k_min_C[ij], &real_min_k, &real_max_k, &min_l_real, &max_l_real ); prepareArray( &matrices->E_C[ij], matrices->k_min_C[ij], matrices->k_max_C[ij], matrices->l_min_C[ij], matrices->l_max_C[ij] ); #ifdef COUNT_STATES prepareArray2( &matrices->N_C[ij], matrices->k_min_C[ij], matrices->k_max_C[ij], matrices->l_min_C[ij], matrices->l_max_C[ij] ); #endif /* d1 and d2 are the distancies to both references introduced by closing a hairpin structure at (i,j) */ if((d1 >= 0) && (d2 >= 0)){ if(((unsigned int)d1<=maxD1) && ((unsigned int)d2 <= maxD2)){ matrices->E_C[ij][d1][d2/2] = (no_close) ? FORBIDDEN : E_Hairpin(dij, type, S1[i+1], S1[j-1], sequence+i-1, P); updatePosteriorBoundaries(d1, d2, &real_min_k, &real_max_k, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_C[ij][d1][d2/2] = 1; #endif } else{ matrices->E_C_rem[ij] = (no_close) ? FORBIDDEN : E_Hairpin(dij, type, S1[i+1], S1[j-1], sequence+i-1, P); } } /* INTERIOR LOOP STRUCTURES */ maxp = MIN2(j-2-TURN,i+MAXLOOP+1); for(p = i+1; p <= maxp; p++){ unsigned int minq = p + TURN + 1; unsigned int ln_pre = dij + p; if(ln_pre > minq + MAXLOOP) minq = ln_pre - MAXLOOP - 1; for(q = minq; q < j; q++){ pq = my_iindx[p]-q; /* set distance to reference structure... */ type_2 = ptype[jindx[q] + p]; if (type_2==0) continue; type_2 = rtype[type_2]; /* get distance to reference if closing the interior loop * d2 = dbp(S_{i,j}, S_{p.q} + {i,j}) */ d1 = base_d1 + referenceBPs1[ij] - referenceBPs1[pq]; d2 = base_d2 + referenceBPs2[ij] - referenceBPs2[pq]; if(no_closingGU) if(no_close||(type_2==3)||(type_2==4)) if((p>i+1)||(qE_C[pq] != NULL){ for(cnt1 = matrices->k_min_C[pq]; cnt1 <= matrices->k_max_C[pq]; cnt1++){ for(cnt2 = matrices->l_min_C[pq][cnt1]; cnt2 <= matrices->l_max_C[pq][cnt1]; cnt2+=2){ if(matrices->E_C[pq][cnt1][cnt2/2] != INF){ if(((cnt1 + d1) <= maxD1) && ((cnt2+d2) <= maxD2)){ matrices->E_C[ij][cnt1 + d1][(cnt2 + d2)/2] = MIN2( matrices->E_C[ij][cnt1 + d1][(cnt2 + d2)/2], matrices->E_C[pq][cnt1][cnt2/2] + energy ); updatePosteriorBoundaries(cnt1 + d1, cnt2 + d2, &real_min_k, &real_max_k, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_C[ij][cnt1 + d1][(cnt2 + d2)/2] += matrices->N_C[pq][cnt1][cnt2/2]; #endif } /* collect all cases where d1+cnt1 or d2+cnt2 exceeds maxD1, maxD2, respectively */ else{ matrices->E_C_rem[ij] = MIN2(matrices->E_C_rem[ij], matrices->E_C[pq][cnt1][cnt2/2] + energy); } } } } } /* collect all contributions where C[pq] already lies outside k_max, l_max boundary */ if(matrices->E_C_rem[pq] != INF){ matrices->E_C_rem[ij] = MIN2(matrices->E_C_rem[ij], matrices->E_C_rem[pq] + energy); } } /* end q-loop */ } /* end p-loop */ /* MULTI LOOP STRUCTURES */ if(!no_close){ /* dangle energies for multiloop closing stem */ tt = rtype[type]; temp2 = P->MLclosing; if(dangles == 2) temp2 += E_MLstem(tt, S1[j-1], S1[i+1], P); else temp2 += E_MLstem(tt, -1, -1, P); for(u=i+TURN+2; uE_M_rem[i1u] != INF){ for(cnt3 = matrices->k_min_M1[u1j1]; cnt3 <= matrices->k_max_M1[u1j1]; cnt3++) for(cnt4 = matrices->l_min_M1[u1j1][cnt3]; cnt4 <= matrices->l_max_M1[u1j1][cnt3]; cnt4+=2){ if(matrices->E_M1[u1j1][cnt3][cnt4/2]!= INF){ matrices->E_C_rem[ij] = MIN2(matrices->E_C_rem[ij], matrices->E_M_rem[i1u] + matrices->E_M1[u1j1][cnt3][cnt4/2] + temp2 ); } } if(matrices->E_M1_rem[u1j1] != INF){ matrices->E_C_rem[ij] = MIN2(matrices->E_C_rem[ij], matrices->E_M_rem[i1u] + matrices->E_M1_rem[u1j1] + temp2 ); } } if(matrices->E_M1_rem[u1j1] != INF){ for(cnt1 = matrices->k_min_M[i1u]; cnt1 <= matrices->k_max_M[i1u]; cnt1++) for(cnt2 = matrices->l_min_M[i1u][cnt1]; cnt2 <= matrices->l_max_M[i1u][cnt1]; cnt2+=2) if(matrices->E_M[i1u][cnt1][cnt2/2] != INF){ matrices->E_C_rem[ij] = MIN2(matrices->E_C_rem[ij], matrices->E_M[i1u][cnt1][cnt2/2] + matrices->E_M1_rem[u1j1] + temp2 ); } } /* get distance to reference if closing the multiloop * d = dbp(S_{i,j}, {i,j} + S_{i+1,u} + S_{u+1,j-1}) */ if(!matrices->E_M[i1u]) continue; if(!matrices->E_M1[u1j1]) continue; d1 = base_d1 + referenceBPs1[ij] - referenceBPs1[i1u] - referenceBPs1[u1j1]; d2 = base_d2 + referenceBPs2[ij] - referenceBPs2[i1u] - referenceBPs2[u1j1]; for(cnt1 = matrices->k_min_M[i1u]; cnt1 <= matrices->k_max_M[i1u]; cnt1++) for(cnt2 = matrices->l_min_M[i1u][cnt1]; cnt2 <= matrices->l_max_M[i1u][cnt1]; cnt2+=2) for(cnt3 = matrices->k_min_M1[u1j1]; cnt3 <= matrices->k_max_M1[u1j1]; cnt3++) for(cnt4 = matrices->l_min_M1[u1j1][cnt3]; cnt4 <= matrices->l_max_M1[u1j1][cnt3]; cnt4+=2){ if((matrices->E_M[i1u][cnt1][cnt2/2] != INF) && (matrices->E_M1[u1j1][cnt3][cnt4/2]!= INF)){ if(((cnt1+cnt3+d1) <= maxD1) && ((cnt2+cnt4+d2) <= maxD2)){ matrices->E_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2] = MIN2( matrices->E_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2], matrices->E_M[i1u][cnt1][cnt2/2] + matrices->E_M1[u1j1][cnt3][cnt4/2] + temp2 ); updatePosteriorBoundaries(cnt1 + cnt3 + d1, cnt2 + cnt4 + d2, &real_min_k, &real_max_k, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2] += matrices->N_M[i1u][cnt1][cnt2/2] * matrices->N_M1[u1j1][cnt3][cnt4/2]; #endif } /* collect all cases where d1+cnt1+cnt3 or d2+cnt2+cnt4 exceeds maxD1, maxD2, respectively */ else{ matrices->E_C_rem[ij] = MIN2( matrices->E_C_rem[ij], matrices->E_M[i1u][cnt1][cnt2/2] + matrices->E_M1[u1j1][cnt3][cnt4/2] + temp2 ); } } } } } /* resize and move memory portions of energy matrix E_C */ adjustArrayBoundaries(&matrices->E_C[ij], &matrices->k_min_C[ij], &matrices->k_max_C[ij], &matrices->l_min_C[ij], &matrices->l_max_C[ij], real_min_k, real_max_k, min_l_real, max_l_real ); #ifdef COUNT_STATES /* actually we should adjust the array boundaries here but we might never use the count states option more than once so what....*/ #endif } /* end >> if (pair) << */ /* done with c[i,j], now compute fML[i,j] */ /* free ends ? -----------------------------------------*/ dia = referenceBPs1[ij] - referenceBPs1[my_iindx[i+1]-j]; dib = referenceBPs2[ij] - referenceBPs2[my_iindx[i+1]-j]; dja = referenceBPs1[ij] - referenceBPs1[ij+1]; djb = referenceBPs2[ij] - referenceBPs2[ij+1]; if(dangles==2) temp2 = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P); else temp2 = E_MLstem(type, -1, -1, P); int min_k_guess, max_k_guess, min_l_guess, max_l_guess; int min_k_real_m, max_k_real_m, *min_l_real_m, *max_l_real_m; int min_k_real_m1, max_k_real_m1, *min_l_real_m1, *max_l_real_m1; min_k_guess = min_l_guess = 0; max_k_guess = mm1[ij] + referenceBPs1[ij]; max_l_guess = mm2[ij] + referenceBPs2[ij]; prepareBoundaries(min_k_guess, max_k_guess, min_l_guess, max_l_guess, bpdist[ij], &matrices->k_min_M[ij], &matrices->k_max_M[ij], &matrices->l_min_M[ij], &matrices->l_max_M[ij] ); prepareBoundaries(min_k_guess, max_k_guess, min_l_guess, max_l_guess, bpdist[ij], &matrices->k_min_M1[ij], &matrices->k_max_M1[ij], &matrices->l_min_M1[ij], &matrices->l_max_M1[ij] ); preparePosteriorBoundaries( matrices->k_max_M[ij] - matrices->k_min_M[ij] + 1, matrices->k_min_M[ij], &min_k_real_m, &max_k_real_m, &min_l_real_m, &max_l_real_m ); preparePosteriorBoundaries( matrices->k_max_M1[ij] - matrices->k_min_M1[ij] + 1, matrices->k_min_M1[ij], &min_k_real_m1, &max_k_real_m1, &min_l_real_m1, &max_l_real_m1 ); prepareArray( &matrices->E_M[ij], matrices->k_min_M[ij], matrices->k_max_M[ij], matrices->l_min_M[ij], matrices->l_max_M[ij] ); prepareArray( &matrices->E_M1[ij], matrices->k_min_M1[ij], matrices->k_max_M1[ij], matrices->l_min_M1[ij], matrices->l_max_M1[ij] ); #ifdef COUNT_STATES prepareArray2( &matrices->N_M[ij], matrices->k_min_M[ij], matrices->k_max_M[ij], matrices->l_min_M[ij], matrices->l_max_M[ij] ); prepareArray2( &matrices->N_M1[ij], matrices->k_min_M1[ij], matrices->k_max_M1[ij], matrices->l_min_M1[ij], matrices->l_max_M1[ij] ); #endif /* now to the actual computations... */ /* 1st E_M[ij] = E_M1[ij] = E_C[ij] + b */ if(matrices->E_C_rem[ij] != INF){ matrices->E_M_rem[ij] = matrices->E_M1_rem[ij] = temp2 + matrices->E_C_rem[ij]; } if(matrices->E_C[ij]) for(cnt1 = matrices->k_min_C[ij]; cnt1 <= matrices->k_max_C[ij]; cnt1++){ for(cnt2 = matrices->l_min_C[ij][cnt1]; cnt2 <= matrices->l_max_C[ij][cnt1]; cnt2+=2){ if(matrices->E_C[ij][cnt1][cnt2/2] != INF){ matrices->E_M[ij][cnt1][cnt2/2] = matrices->E_M1[ij][cnt1][cnt2/2] = temp2 + matrices->E_C[ij][cnt1][cnt2/2]; updatePosteriorBoundaries(cnt1, cnt2, &min_k_real_m, &max_k_real_m, &min_l_real_m, &max_l_real_m ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real_m1, &max_k_real_m1, &min_l_real_m1, &max_l_real_m1 ); #ifdef COUNT_STATES matrices->N_M[ij][cnt1][cnt2/2] = matrices->N_M1[ij][cnt1][cnt2/2] = matrices->N_C[ij][cnt1][cnt2/2]; #endif } } } /* 2nd E_M[ij] = MIN(E_M[ij], E_M[i+1,j] + c) */ if(matrices->E_M_rem[my_iindx[i+1]-j] != INF){ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M_rem[my_iindx[i+1]-j] + P->MLbase ); } if(matrices->E_M[my_iindx[i+1]-j]) for(cnt1 = matrices->k_min_M[my_iindx[i+1]-j]; cnt1 <= matrices->k_max_M[my_iindx[i+1]-j]; cnt1++){ for(cnt2 = matrices->l_min_M[my_iindx[i+1]-j][cnt1]; cnt2 <= matrices->l_max_M[my_iindx[i+1]-j][cnt1]; cnt2+=2){ if(matrices->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] != INF){ if(((cnt1 + dia) <= maxD1) && ((cnt2 + dib) <= maxD2)){ matrices->E_M[ij][cnt1+dia][(cnt2+dib)/2] = MIN2( matrices->E_M[ij][cnt1+dia][(cnt2+dib)/2], matrices->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase ); updatePosteriorBoundaries(cnt1 + dia, cnt2 + dib, &min_k_real_m, &max_k_real_m, &min_l_real_m, &max_l_real_m ); #ifdef COUNT_STATES matrices->N_M[ij][cnt1+dia][(cnt2+dib)/2] += matrices->N_M[my_iindx[i+1]-j][cnt1][cnt2/2]; #endif } /* collect all cases where dia+cnt1 or dib+cnt2 exceeds maxD1, maxD2, respectively */ else{ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase ); } } } } /* 3rd E_M[ij] = MIN(E_M[ij], E_M[i,j-1] + c) */ if(matrices->E_M_rem[ij+1] != INF){ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M_rem[ij+1] + P->MLbase ); } if(matrices->E_M[ij+1]) for(cnt1 = matrices->k_min_M[ij+1]; cnt1 <= matrices->k_max_M[ij+1]; cnt1++){ for(cnt2 = matrices->l_min_M[ij+1][cnt1]; cnt2 <= matrices->l_max_M[ij+1][cnt1]; cnt2+=2){ if(matrices->E_M[ij+1][cnt1][cnt2/2] != INF){ if(((cnt1 + dja) <= maxD1) && ((cnt2 + djb) <= maxD2)){ matrices->E_M[ij][cnt1+dja][(cnt2+djb)/2] = MIN2( matrices->E_M[ij][cnt1+dja][(cnt2+djb)/2], matrices->E_M[ij+1][cnt1][cnt2/2] + P->MLbase ); updatePosteriorBoundaries(cnt1 + dja, cnt2 + djb, &min_k_real_m, &max_k_real_m, &min_l_real_m, &max_l_real_m ); #ifdef COUNT_STATES matrices->N_M[ij][cnt1+dja][(cnt2+djb)/2] += matrices->N_M[ij+1][cnt1][cnt2/2]; #endif } /* collect all cases where dja+cnt1 or djb+cnt2 exceeds maxD1, maxD2, respectively */ else{ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M[ij+1][cnt1][cnt2/2] + P->MLbase ); } } } } /* 4th E_M1[ij] = MIN(E_M1[ij], E_M1[i,j-1] + c) */ if(matrices->E_M1_rem[ij+1] != INF){ matrices->E_M1_rem[ij] = MIN2( matrices->E_M1_rem[ij], matrices->E_M1_rem[ij+1] + P->MLbase ); } if(matrices->E_M1[ij+1]) for(cnt1 = matrices->k_min_M1[ij+1]; cnt1 <= matrices->k_max_M1[ij+1]; cnt1++){ for(cnt2 = matrices->l_min_M1[ij+1][cnt1]; cnt2 <= matrices->l_max_M1[ij+1][cnt1]; cnt2+=2){ if(matrices->E_M1[ij+1][cnt1][cnt2/2] != INF){ if(((cnt1 + dja) <= maxD1) && ((cnt2 + djb) <= maxD2)){ matrices->E_M1[ij][cnt1+dja][(cnt2+djb)/2] = MIN2( matrices->E_M1[ij][cnt1+dja][(cnt2+djb)/2], matrices->E_M1[ij+1][cnt1][cnt2/2] + P->MLbase ); updatePosteriorBoundaries(cnt1 + dja, cnt2 + djb, &min_k_real_m1, &max_k_real_m1, &min_l_real_m1, &max_l_real_m1 ); #ifdef COUNT_STATES matrices->N_M1[ij][cnt1+dja][(cnt2+djb)/2] += matrices->N_M1[ij+1][cnt1][cnt2/2]; #endif } /* collect all cases where dja+cnt1 or djb+cnt2 exceeds maxD1, maxD2, respectively */ else{ matrices->E_M1_rem[ij] = MIN2( matrices->E_M1_rem[ij], matrices->E_M1[ij+1][cnt1][cnt2/2] + P->MLbase ); } } } } /* 5th E_M[ij] = MIN(E_M[ij], min(E_M[i,k] + E_M[k+1,j])) */ if(j > TURN + 2) for (u = i+1+TURN; u <= j-2-TURN; u++){ /* check all cases where M(i,u) and/or M(u+1,j) are already out of scope of maxD1 and/or maxD2 */ if(matrices->E_M_rem[my_iindx[i]-u] != INF){ for(cnt3 = matrices->k_min_M[my_iindx[u+1]-j]; cnt3 <= matrices->k_max_M[my_iindx[u+1]-j]; cnt3++){ for(cnt4 = matrices->l_min_M[my_iindx[u+1]-j][cnt3]; cnt4 <= matrices->l_max_M[my_iindx[u+1]-j][cnt3]; cnt4+=2){ if(matrices->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] != INF){ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M_rem[my_iindx[i]-u] + matrices->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] ); } } } if(matrices->E_M_rem[my_iindx[u+1]-j] != INF){ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M_rem[my_iindx[i]-u] + matrices->E_M_rem[my_iindx[u+1]-j] ); } } if(matrices->E_M_rem[my_iindx[u+1]-j] != INF){ for(cnt1 = matrices->k_min_M[my_iindx[i]-u]; cnt1 <= matrices->k_max_M[my_iindx[i]-u]; cnt1++){ for(cnt2 = matrices->l_min_M[my_iindx[i]-u][cnt1]; cnt2 <= matrices->l_max_M[my_iindx[i]-u][cnt1]; cnt2+=2){ if(matrices->E_M[my_iindx[i]-u][cnt1][cnt2/2] != INF){ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M[my_iindx[i]-u][cnt1][cnt2/2] + matrices->E_M_rem[my_iindx[u+1]-j] ); } } } } if(!matrices->E_M[my_iindx[i]-u]) continue; if(!matrices->E_M[my_iindx[u+1]-j]) continue; dia = referenceBPs1[ij] - referenceBPs1[my_iindx[i]-u] - referenceBPs1[my_iindx[u+1]-j]; dib = referenceBPs2[ij] - referenceBPs2[my_iindx[i]-u] - referenceBPs2[my_iindx[u+1]-j]; for(cnt1 = matrices->k_min_M[my_iindx[i]-u]; cnt1 <= matrices->k_max_M[my_iindx[i]-u]; cnt1++){ for(cnt2 = matrices->l_min_M[my_iindx[i]-u][cnt1]; cnt2 <= matrices->l_max_M[my_iindx[i]-u][cnt1]; cnt2+=2){ for(cnt3 = matrices->k_min_M[my_iindx[u+1]-j]; cnt3 <= matrices->k_max_M[my_iindx[u+1]-j]; cnt3++){ for(cnt4 = matrices->l_min_M[my_iindx[u+1]-j][cnt3]; cnt4 <= matrices->l_max_M[my_iindx[u+1]-j][cnt3]; cnt4+=2){ if((matrices->E_M[my_iindx[i]-u][cnt1][cnt2/2] != INF) && (matrices->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] != INF)){ if(((cnt1 + cnt3 + dia) <= maxD1) && ((cnt2 + cnt4 + dib) <= maxD2)){ matrices->E_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2] = MIN2( matrices->E_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2], matrices->E_M[my_iindx[i]-u][cnt1][cnt2/2] + matrices->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] ); updatePosteriorBoundaries(cnt1 + cnt3 + dia, cnt2 + cnt4 + dib, &min_k_real_m, &max_k_real_m, &min_l_real_m, &max_l_real_m ); #ifdef COUNT_STATES matrices->N_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2] += matrices->N_M[my_iindx[i]-u][cnt1][cnt2/2] * matrices->N_M1[my_iindx[u+1]-j][cnt3][cnt4/2]; #endif } /* collect all cases where dia+cnt1+cnt3 or dib+cnt2+cnt4 exceeds maxD1, maxD2, respectively */ else{ matrices->E_M_rem[ij] = MIN2(matrices->E_M_rem[ij], matrices->E_M[my_iindx[i]-u][cnt1][cnt2/2] + matrices->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] ); } } } } } } } /* thats all folks for the multiloop decomposition... */ adjustArrayBoundaries(&matrices->E_M[ij], &matrices->k_min_M[ij], &matrices->k_max_M[ij], &matrices->l_min_M[ij], &matrices->l_max_M[ij], min_k_real_m, max_k_real_m, min_l_real_m, max_l_real_m ); adjustArrayBoundaries(&matrices->E_M1[ij], &matrices->k_min_M1[ij], &matrices->k_max_M1[ij], &matrices->l_min_M1[ij], &matrices->l_max_M1[ij], min_k_real_m1, max_k_real_m1, min_l_real_m1, max_l_real_m1 ); #ifdef COUNT_STATES /* actually we should adjust the array boundaries here but we might never use the count states option more than once so what....*/ #endif } /* end of j-loop */ } /* calculate energies of 5' and 3' fragments */ /* prepare first entries in E_F5 */ for(cnt1 = 1; cnt1 <= TURN+1; cnt1++){ matrices->E_F5[cnt1] = (int **)vrna_alloc(sizeof(int *)); matrices->E_F5[cnt1][0] = (int *)vrna_alloc(sizeof(int)); matrices->E_F5[cnt1][0][0] = 0; matrices->E_F5_rem[cnt1] = INF; matrices->k_min_F5[cnt1] = matrices->k_max_F5[cnt1] = 0; matrices->l_min_F5[cnt1] = (int *)vrna_alloc(sizeof(int)); matrices->l_max_F5[cnt1] = (int *)vrna_alloc(sizeof(int)); matrices->l_min_F5[cnt1][0] = matrices->l_max_F5[cnt1][0] = 0; #ifdef COUNT_STATES matrices->N_F5[cnt1] = (unsigned long **)vrna_alloc(sizeof(unsigned long *)); matrices->N_F5[cnt1][0] = (unsigned long *)vrna_alloc(sizeof(unsigned long)); matrices->N_F5[cnt1][0][0] = 1; #endif } for (j=TURN+2; j <= seq_length; j++) { unsigned int da = referenceBPs1[my_iindx[1]-j] - referenceBPs1[my_iindx[1]-j+1]; unsigned int db = referenceBPs2[my_iindx[1]-j] - referenceBPs2[my_iindx[1]-j+1]; type=ptype[jindx[j] + 1]; additional_en = 0; if(type){ if(dangles == 2) additional_en += E_ExtLoop(type, -1, j < seq_length ? S1[j+1] : -1, P); else additional_en += E_ExtLoop(type, -1, -1, P); } /* make min and max k guess for memory allocation */ int min_k_guess, max_k_guess, min_l_guess, max_l_guess; int *min_l_real, *max_l_real, min_k_real, max_k_real; min_k_guess = min_l_guess = 0; max_k_guess = referenceBPs1[my_iindx[1]-j] + mm1[my_iindx[1]-j]; max_l_guess = referenceBPs2[my_iindx[1]-j] + mm2[my_iindx[1]-j]; prepareBoundaries(min_k_guess, max_k_guess, min_l_guess, max_l_guess, bpdist[my_iindx[1]-j], &matrices->k_min_F5[j], &matrices->k_max_F5[j], &matrices->l_min_F5[j], &matrices->l_max_F5[j] ); preparePosteriorBoundaries( matrices->k_max_F5[j] - matrices->k_min_F5[j] + 1, matrices->k_min_F5[j], &min_k_real, &max_k_real, &min_l_real, &max_l_real ); prepareArray( &matrices->E_F5[j], matrices->k_min_F5[j], matrices->k_max_F5[j], matrices->l_min_F5[j], matrices->l_max_F5[j] ); #ifdef COUNT_STATES prepareArray2( &matrices->N_F5[j], matrices->k_min_F5[j], matrices->k_max_F5[j], matrices->l_min_F5[j], matrices->l_max_F5[j] ); #endif /* begin the actual computation of 5' end energies */ /* j-1 is unpaired ... */ matrices->E_F5_rem[j] = matrices->E_F5_rem[j-1]; for(cnt1 = matrices->k_min_F5[j-1]; cnt1 <= matrices->k_max_F5[j-1]; cnt1++){ for(cnt2 = matrices->l_min_F5[j-1][cnt1]; cnt2 <= matrices->l_max_F5[j-1][cnt1]; cnt2+=2){ if(((cnt1 + da) <= maxD1) && ((cnt2 + db) <= maxD2)){ matrices->E_F5[j][cnt1+da][(cnt2+db)/2] = MIN2( matrices->E_F5[j][cnt1+da][(cnt2+db)/2], matrices->E_F5[j-1][cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1 + da, cnt2 + db, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_F5[j][cnt1+da][(cnt2+db)/2] += matrices->N_F5[j-1][cnt1][cnt2/2]; #endif } /* collect all cases where da+cnt1 or db+cnt2 exceeds maxD1, maxD2, respectively */ else{ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_F5[j-1][cnt1][cnt2/2]); } } } /* j pairs with 1 */ if(matrices->E_C_rem[my_iindx[1]-j] != INF){ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_C_rem[my_iindx[1]-j] + additional_en); } if(matrices->E_C[my_iindx[1]-j]) for(cnt1 = matrices->k_min_C[my_iindx[1]-j]; cnt1 <= matrices->k_max_C[my_iindx[1]-j]; cnt1++) for(cnt2 = matrices->l_min_C[my_iindx[1]-j][cnt1]; cnt2 <= matrices->l_max_C[my_iindx[1]-j][cnt1]; cnt2+=2){ if(matrices->E_C[my_iindx[1]-j][cnt1][cnt2/2] != INF){ matrices->E_F5[j][cnt1][cnt2/2] = MIN2( matrices->E_F5[j][cnt1][cnt2/2], matrices->E_C[my_iindx[1]-j][cnt1][cnt2/2]+ additional_en ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_F5[j][cnt1][cnt2/2] += matrices->N_C[my_iindx[1]-j][cnt1][cnt2/2]; #endif } } /* j pairs with some other nucleotide -> see below */ for (i=j-TURN-1; i>1; i--) { ij = my_iindx[i]-j; type = ptype[jindx[j] + i]; if (type) { if(dangles == 2) additional_en = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P); else additional_en = E_ExtLoop(type, -1, -1, P); if(matrices->E_C_rem[ij] != INF){ for(cnt3 = matrices->k_min_F5[i-1]; cnt3 <= matrices->k_max_F5[i-1]; cnt3++) for(cnt4 = matrices->l_min_F5[i-1][cnt3]; cnt4 <= matrices->l_max_F5[i-1][cnt3]; cnt4+=2){ if(matrices->E_F5[i-1][cnt3][cnt4/2] != INF){ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_F5[i-1][cnt3][cnt4/2] + matrices->E_C_rem[ij] + additional_en ); } } if(matrices->E_F5_rem[i-1] != INF){ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_F5_rem[i-1] + matrices->E_C_rem[ij] + additional_en ); } } if((matrices->E_F5_rem[i-1] != INF) && (matrices->E_C[ij])){ for(cnt1 = matrices->k_min_C[ij]; cnt1 <= matrices->k_max_C[ij]; cnt1++) for(cnt2 = matrices->l_min_C[ij][cnt1]; cnt2 <= matrices->l_max_C[ij][cnt1]; cnt2+=2) if(matrices->E_C[ij][cnt1][cnt2/2]!= INF){ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_F5_rem[i-1] + matrices->E_C[ij][cnt1][cnt2/2] + additional_en ); } } if(!matrices->E_C[ij]) continue; unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1]; unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1]; for(cnt1 = matrices->k_min_C[ij]; cnt1 <= matrices->k_max_C[ij]; cnt1++) for(cnt2 = matrices->l_min_C[ij][cnt1]; cnt2 <= matrices->l_max_C[ij][cnt1]; cnt2+=2) for(cnt3 = matrices->k_min_F5[i-1]; cnt3 <= matrices->k_max_F5[i-1]; cnt3++) for(cnt4 = matrices->l_min_F5[i-1][cnt3]; cnt4 <= matrices->l_max_F5[i-1][cnt3]; cnt4+=2){ if(matrices->E_F5[i-1][cnt3][cnt4/2] != INF && matrices->E_C[ij][cnt1][cnt2/2]!= INF){ if(((cnt1 + cnt3 + d1a) <= maxD1) && ((cnt2 + cnt4 + d1b) <= maxD2)){ matrices->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] = MIN2( matrices->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2], matrices->E_F5[i-1][cnt3][cnt4/2] + matrices->E_C[ij][cnt1][cnt2/2] + additional_en ); updatePosteriorBoundaries(cnt1 + cnt3 + d1a, cnt2 + cnt4 + d1b, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] += matrices->N_F5[i-1][cnt3][cnt4/2] * matrices->N_C[ij][cnt1][cnt2/2]; #endif } /* collect all cases where d1a+cnt1+cnt3 or d1b+cnt2+cnt4 exceeds maxD1, maxD2, respectively */ else{ matrices->E_F5_rem[j] = MIN2(matrices->E_F5_rem[j], matrices->E_F5[i-1][cnt3][cnt4/2] + matrices->E_C[ij][cnt1][cnt2/2] + additional_en ); } } } } } /* resize and move memory portions of energy matrix E_F5 */ adjustArrayBoundaries(&matrices->E_F5[j], &matrices->k_min_F5[j], &matrices->k_max_F5[j], &matrices->l_min_F5[j], &matrices->l_max_F5[j], min_k_real, max_k_real, min_l_real, max_l_real ); } /* end of j-loop */ if(compute_2Dfold_F3){ /* prepare first entries in E_F3 */ for(cnt1 = seq_length; cnt1 >= seq_length-TURN-1; cnt1--){ matrices->E_F3[cnt1] = (int **)vrna_alloc(sizeof(int *)); matrices->E_F3[cnt1][0] = (int *) vrna_alloc(sizeof(int)); matrices->E_F3[cnt1][0][0] = 0; matrices->k_min_F3[cnt1] = matrices->k_max_F3[cnt1] = 0; matrices->l_min_F3[cnt1] = (int *)vrna_alloc(sizeof(int)); matrices->l_max_F3[cnt1] = (int *)vrna_alloc(sizeof(int)); matrices->l_min_F3[cnt1][0] = matrices->l_max_F3[cnt1][0] = 0; } /* begin calculations */ for (j=seq_length-TURN-2; j >= 1; j--){ unsigned int da = referenceBPs1[my_iindx[j]-seq_length] - referenceBPs1[my_iindx[j+1]-seq_length]; unsigned int db = referenceBPs2[my_iindx[j]-seq_length] - referenceBPs2[my_iindx[j+1]-seq_length]; type=ptype[jindx[seq_length] + j]; additional_en = 0; if(type){ if(dangles == 2) additional_en += E_ExtLoop(type, j > 1 ? S1[j-1] : -1, -1, P); else additional_en += E_ExtLoop(type, -1, -1, P); } /* make min and max k guess for memory allocation */ int min_k_guess, max_k_guess, min_l_guess, max_l_guess; int *min_l_real, *max_l_real, min_k_real, max_k_real; min_k_guess = min_l_guess = 0; max_k_guess = referenceBPs1[my_iindx[j]-seq_length] + mm1[my_iindx[j]-seq_length]; max_l_guess = referenceBPs2[my_iindx[j]-seq_length] + mm2[my_iindx[j]-seq_length]; prepareBoundaries(min_k_guess, max_k_guess, min_l_guess, max_l_guess, bpdist[my_iindx[j]-seq_length], &matrices->k_min_F3[j], &matrices->k_max_F3[j], &matrices->l_min_F3[j], &matrices->l_max_F3[j] ); preparePosteriorBoundaries( matrices->k_max_F3[j] - matrices->k_min_F3[j] + 1, matrices->k_min_F3[j], &min_k_real, &max_k_real, &min_l_real, &max_l_real ); prepareArray( &matrices->E_F3[j], matrices->k_min_F3[j], matrices->k_max_F3[j], matrices->l_min_F3[j], matrices->l_max_F3[j] ); /* begin the actual computation of 5' end energies */ /* j is unpaired ... */ for(cnt1 = matrices->k_min_F3[j+1]; cnt1 <= matrices->k_max_F3[j+1]; cnt1++){ for(cnt2 = matrices->l_min_F3[j+1][cnt1]; cnt2 <= matrices->l_max_F3[j+1][cnt1]; cnt2+=2){ matrices->E_F3[j][cnt1+da][(cnt2+db)/2] = MIN2( matrices->E_F3[j][cnt1+da][(cnt2+db)/2], matrices->E_F3[j+1][cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1 + da, cnt2 + db, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } } /* j pairs with n */ if(matrices->E_C[my_iindx[j]-seq_length]) for(cnt1 = matrices->k_min_C[my_iindx[j]-seq_length]; cnt1 <= matrices->k_max_C[my_iindx[j]-seq_length]; cnt1++) for(cnt2 = matrices->l_min_C[my_iindx[j]-seq_length][cnt1]; cnt2 <= matrices->l_max_C[my_iindx[j]-seq_length][cnt1]; cnt2+=2){ if(matrices->E_C[my_iindx[j]-seq_length][cnt1][cnt2/2] != INF){ matrices->E_F3[j][cnt1][cnt2/2] = MIN2( matrices->E_F3[j][cnt1][cnt2/2], matrices->E_C[my_iindx[j]-seq_length][cnt1][cnt2/2]+ additional_en ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } } /* j pairs with some other nucleotide -> see below */ for (i=j-TURN-1; i>1; i--) { ij = my_iindx[i]-j; if(!matrices->E_C[ij]) continue; type = ptype[jindx[j] + i]; if (type) { unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1]; unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1]; if(dangles == 2) additional_en = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P); else additional_en = E_ExtLoop(type, -1, -1, P); for(cnt1 = matrices->k_min_C[ij]; cnt1 <= matrices->k_max_C[ij]; cnt1++) for(cnt2 = matrices->l_min_C[ij][cnt1]; cnt2 <= matrices->l_max_C[ij][cnt1]; cnt2+=2) for(cnt3 = matrices->k_min_F5[i-1]; cnt3 <= matrices->k_max_F5[i-1]; cnt3++) for(cnt4 = matrices->l_min_F5[i-1][cnt3]; cnt4 <= matrices->l_max_F5[i-1][cnt3]; cnt4+=2){ if(matrices->E_F5[i-1][cnt3][cnt4/2] != INF && matrices->E_C[ij][cnt1][cnt2/2]!= INF){ matrices->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] = MIN2( matrices->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2], matrices->E_F5[i-1][cnt3][cnt4/2] + matrices->E_C[ij][cnt1][cnt2/2] + additional_en ); updatePosteriorBoundaries(cnt1 + cnt3 + d1a, cnt2 + cnt4 + d1b, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); #ifdef COUNT_STATES matrices->N_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] += matrices->N_F5[i-1][cnt3][cnt4/2] * matrices->N_C[ij][cnt1][cnt2/2]; #endif } } } } /* resize and move memory portions of energy matrix E_F5 */ adjustArrayBoundaries(&matrices->E_F5[j], &matrices->k_min_F5[j], &matrices->k_max_F5[j], &matrices->l_min_F5[j], &matrices->l_max_F5[j], min_k_real, max_k_real, min_l_real, max_l_real ); } /* end of j-loop */ } } /*---------------------------------------------------------------------------*/ /*---------------------------------------------------------------------------*/ PRIVATE void backtrack_f5( unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc){ int *my_iindx, *jindx, energy, type, dangles, cnt1, cnt2, cnt3, cnt4; int **l_min_C, **l_max_C,**l_min_F5, **l_max_F5; int *k_min_C, *k_max_C,*k_min_F5, *k_max_F5; int ***E_C, ***E_F5; int *E_C_rem, *E_F5_rem; unsigned int i, ij, seq_length, maxD1, maxD2; short *S1; unsigned int *referenceBPs1, *referenceBPs2; char *ptype; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; unsigned int da, db; P = vc->params; md = &(P->model_details); matrices = vc->matrices; seq_length = vc->length; S1 = vc->sequence_encoding; ptype = vc->ptype; my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; dangles = md->dangles; E_F5 = matrices->E_F5; l_min_F5 = matrices->l_min_F5; l_max_F5 = matrices->l_max_F5; k_min_F5 = matrices->k_min_F5; k_max_F5 = matrices->k_max_F5; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_F5_rem = matrices->E_F5_rem; E_C_rem = matrices->E_C_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; da = referenceBPs1[my_iindx[1]-j] - referenceBPs1[my_iindx[1]-j+1]; db = referenceBPs2[my_iindx[1]-j] - referenceBPs2[my_iindx[1]-j+1]; if(j maxD1) || ((cnt2 + db) > maxD2)){ if(E_F5_rem[j] == E_F5[j-1][cnt1][cnt2/2]){ backtrack_f5(j-1, cnt1, cnt2, structure, vc); return; } } } } } } else if((k >= da) && (l >= db)){ if(E_F5[j-1]){ if((k - da >= k_min_F5[j-1]) && (k - da <= k_max_F5[j-1])){ if((l - db >= l_min_F5[j-1][k-da]) && (l - db <= l_max_F5[j-1][k-da])) if(E_F5[j-1][k-da][(l-db)/2] == E_F5[j][k][l/2]){ backtrack_f5(j-1, k-da, l-db, structure, vc); return; } } } } type = ptype[jindx[j] + 1]; if(type){ if(dangles == 2) energy = E_ExtLoop(type, -1, j < seq_length ? S1[j+1] : -1, P); else energy = E_ExtLoop(type, -1, -1, P); if(k == -1){ if(E_C_rem[my_iindx[1]-j] + energy == E_F5_rem[j]){ backtrack_c(1, j, -1, -1, structure, vc); return; } } else if(k >= k_min_C[my_iindx[1]-j] && (k <= k_max_C[my_iindx[1]-j])){ if((l >= l_min_C[my_iindx[1]-j][k]) && (l <= l_max_C[my_iindx[1]-j][k])) if(E_C[my_iindx[1]-j][k][l/2] + energy == E_F5[j][k][l/2]){ backtrack_c(1, j, k, l, structure, vc); return; } } } for (i=j-TURN-1; i>1; i--) { ij = my_iindx[i]-j; type = ptype[jindx[j] + i]; if (type) { unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1]; unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1]; if(dangles == 2) energy = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P); else energy = E_ExtLoop(type, -1, -1, P); if(k == -1){ if(E_C_rem[ij] != INF){ for(cnt1 = k_min_F5[i-1]; cnt1 <= k_max_F5[i-1]; cnt1++){ for(cnt2 = l_min_F5[i-1][cnt1]; cnt2 <= l_max_F5[i-1][cnt1]; cnt2+=2){ if(E_F5_rem[j] == (E_F5[i-1][cnt1][cnt2/2] + E_C_rem[ij] + energy)){ backtrack_f5(i-1, cnt1, cnt2, structure, vc); backtrack_c(i,j,-1,-1,structure, vc); return; } } } if(E_F5_rem[j] == (E_F5_rem[i-1] + E_C_rem[ij] + energy)){ backtrack_f5(i-1, -1, -1, structure, vc); backtrack_c(i,j,-1,-1,structure,vc); return; } } if(E_F5_rem[i-1] != INF){ for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++){ for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2){ if(E_F5_rem[j] == (E_F5_rem[i-1] + E_C[ij][cnt1][cnt2/2] + energy)){ backtrack_f5(i-1,-1,-1,structure,vc); backtrack_c(i,j,cnt1,cnt2,structure,vc); return; } } } } for(cnt1 = k_min_F5[i-1]; cnt1 <= k_max_F5[i-1]; cnt1++) for(cnt2 = l_min_F5[i-1][cnt1]; cnt2 <= l_max_F5[i-1][cnt1]; cnt2 += 2) for(cnt3 = k_min_C[ij]; cnt3 <= k_max_C[ij]; cnt3++) for(cnt4 = l_min_C[ij][cnt3]; cnt4 <= l_max_C[ij][cnt3]; cnt4 += 2){ if(((cnt1 + cnt3 + d1a)>maxD1) || ((cnt2+cnt4+d1b)>maxD2)){ if(E_F5_rem[j] == (E_F5[i-1][cnt1][cnt2/2] + E_C[ij][cnt3][cnt4/2] + energy)){ backtrack_f5(i-1,cnt1,cnt2,structure,vc); backtrack_c(i,j,cnt3,cnt4,structure,vc); return; } } } } else if((k >= d1a) && (l >= d1b)){ int k_f_max = MIN2(k-d1a, k_max_F5[i-1]); for(cnt1 = k_min_F5[i-1]; cnt1 <= k_f_max; cnt1++){ int l_f_max = MIN2(l - d1b, l_max_F5[i-1][cnt1]); for(cnt2 = l_min_F5[i-1][cnt1]; cnt2 <= l_f_max; cnt2+=2){ int k_c = k - d1a - cnt1; if((k_c >= k_min_C[ij]) && (k_c <= k_max_C[ij])){ int l_c = l - d1b - cnt2; if((l_c >= l_min_C[ij][k_c]) && (l_c <= l_max_C[ij][k_c])){ if(E_F5[j][k][l/2] == (E_F5[i-1][cnt1][cnt2/2] + E_C[ij][k_c][l_c/2] + energy)){ backtrack_f5(i-1, cnt1, cnt2, structure, vc); backtrack_c(i, j, k_c, l_c, structure, vc); return; } } } } } } } } vrna_message_error("backtracking failed in f5"); } PRIVATE void backtrack_c(unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc){ unsigned int p, q, pq, ij, maxp, maxD1, maxD2; int *my_iindx, *jindx, type, type_2, energy, no_close, dangles, base_d1, base_d2, d1, d2, cnt1, cnt2, cnt3, cnt4, *rtype; int **l_min_C, **l_max_C,**l_min_M, **l_max_M,**l_min_M1, **l_max_M1; int *k_min_C, *k_max_C,*k_min_M, *k_max_M,*k_min_M1, *k_max_M1; int ***E_C, ***E_M, ***E_M1, *E_C_rem, *E_M_rem, *E_M1_rem; short *S1; unsigned int *referenceBPs1, *referenceBPs2; char *ptype, *sequence; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; P = vc->params; md = &(P->model_details); matrices = vc->matrices; sequence = vc->sequence; S1 = vc->sequence_encoding; ptype = vc->ptype; rtype = &(md->rtype[0]); my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; dangles = md->dangles; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_M = matrices->E_M; l_min_M = matrices->l_min_M; l_max_M = matrices->l_max_M; k_min_M = matrices->k_min_M; k_max_M = matrices->k_max_M; E_M1 = matrices->E_M1; l_min_M1 = matrices->l_min_M1; l_max_M1 = matrices->l_max_M1; k_min_M1 = matrices->k_min_M1; k_max_M1 = matrices->k_max_M1; E_C_rem = matrices->E_C_rem; E_M_rem = matrices->E_M_rem; E_M1_rem = matrices->E_M1_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; ij = my_iindx[i]-j; int e = (k==-1) ? E_C_rem[ij] : E_C[ij][k][l/2]; type = ptype[jindx[j] + i]; no_close = (((type==3)||(type==4))&&no_closingGU); structure[i-1] = '('; structure[j-1] = ')'; base_d1 = ((unsigned int)vc->reference_pt1[i] != j) ? 1 : -1; base_d2 = ((unsigned int)vc->reference_pt2[i] != j) ? 1 : -1; base_d1 += referenceBPs1[ij]; base_d2 += referenceBPs2[ij]; if(k == -1){ if(((unsigned int)base_d1 > maxD1) || ((unsigned int)base_d2 > maxD2)){ if(e == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], sequence+i-1, P)) return; } } else{ if((unsigned int)base_d1 == k) if((unsigned int)base_d2 == l) if(E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], sequence+i-1, P) == e) return; } maxp = MIN2(j-2-TURN,i+MAXLOOP+1); for(p = i+1; p <= maxp; p++){ unsigned int minq, ln_pre; minq = p + TURN + 1; ln_pre = j - i - 1; if(ln_pre > minq + MAXLOOP) minq = ln_pre - MAXLOOP - 1; for (q = minq; q < j; q++) { pq = my_iindx[p]-q; type_2 = ptype[jindx[q] + p]; if (type_2==0) continue; type_2 = rtype[type_2]; /* d2 = dbp(S_{i,j}, S_{p.q} + {i,j}) */ d1 = base_d1 - referenceBPs1[pq]; d2 = base_d2 - referenceBPs2[pq]; energy = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1], P); if(k == -1){ if(E_C_rem[pq] != INF) if(e == (E_C_rem[pq] + energy)){ backtrack_c(p,q,-1,-1,structure,vc); return; } if(E_C[pq]) for(cnt1 = k_min_C[pq]; cnt1 <= k_max_C[pq]; cnt1++) for(cnt2 = l_min_C[pq][cnt1]; cnt2 <= l_max_C[pq][cnt1]; cnt2 += 2){ if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){ if(e == (E_C[pq][cnt1][cnt2/2] + energy)){ backtrack_c(p,q,cnt1,cnt2,structure,vc); return; } } } } else{ if(!E_C[pq]) continue; if(d1 <= k && d2 <= l){ if((k-d1 >= k_min_C[pq]) && (k-d1) <= k_max_C[pq]) if((l - d2 >= l_min_C[pq][k-d1]) && (l-d2 <= l_max_C[pq][k-d1])) if(E_C[pq][k-d1][(l-d2)/2] + energy == e){ backtrack_c(p, q, k-d1, l-d2, structure, vc); return; } } } } /* end q-loop */ } /* end p-loop */ /* multi-loop decomposition ------------------------*/ if(!no_close){ unsigned int u; int tt; if(k==-1){ for(u=i+TURN+2; uMLclosing; if(dangles == 2) energy += E_MLstem(tt, S1[j-1], S1[i+1], P); else energy += E_MLstem(tt, -1, -1, P); if(E_M_rem[i1u] != INF){ if(E_M1[u1j1]) for(cnt1 = k_min_M1[u1j1]; cnt1 <= k_max_M1[u1j1]; cnt1++) for(cnt2 = l_min_M1[u1j1][cnt1]; cnt2 <= l_max_M1[u1j1][cnt1]; cnt2 += 2){ if(e == (E_M_rem[i1u] + E_M1[u1j1][cnt1][cnt2/2] + energy)){ backtrack_m(i+1,u,-1,-1,structure,vc); backtrack_m1(u+1,j-1,cnt1,cnt2,structure,vc); return; } } if(E_M1_rem[u1j1] != INF){ if(e == (E_M_rem[i1u] + E_M1_rem[u1j1] + energy)){ backtrack_m(i+1, u, -1, -1, structure, vc); backtrack_m1(u+1, j-1, -1, -1, structure, vc); return; } } } if(E_M1_rem[u1j1] != INF){ if(E_M[i1u]) for(cnt1 = k_min_M[i1u]; cnt1 <= k_max_M[i1u]; cnt1++) for(cnt2 = l_min_M[i1u][cnt1]; cnt2 <= l_max_M[i1u][cnt1]; cnt2 += 2) if(e == (E_M[i1u][cnt1][cnt2/2] + E_M1_rem[u1j1] + energy)){ backtrack_m(i+1,u,cnt1,cnt2,structure,vc); backtrack_m1(u+1,j-1,-1,-1,structure,vc); return; } } /* now all cases where we exceed the maxD1/D2 scope by combination of E_M and E_M1 */ if(!E_M[i1u]) continue; if(!E_M1[u1j1]) continue; /* get distance to reference if closing this multiloop * dist3 = dbp(S_{i,j}, {i,j} + S_{i+1.u} + S_{u+1,j-1}) */ d1 = base_d1 - referenceBPs1[i1u] - referenceBPs1[u1j1]; d2 = base_d2 - referenceBPs2[i1u] - referenceBPs2[u1j1]; for(cnt1 = matrices->k_min_M[i1u]; cnt1 <= matrices->k_max_M[i1u]; cnt1++) for(cnt2 = matrices->l_min_M[i1u][cnt1]; cnt2 <= matrices->l_max_M[i1u][cnt1]; cnt2+=2) for(cnt3 = matrices->k_min_M1[u1j1]; cnt3 <= matrices->k_max_M1[u1j1]; cnt3++) for(cnt4 = matrices->l_min_M1[u1j1][cnt3]; cnt4 <= matrices->l_max_M1[u1j1][cnt3]; cnt4+=2){ if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){ if(e == (E_M[i1u][cnt1][cnt2/2] + E_M1[u1j1][cnt3][cnt4/2] + energy)){ backtrack_m(i+1,u,cnt1,cnt2,structure,vc); backtrack_m1(u+1,j-1,cnt3,cnt4,structure,vc); return; } } } } } else{ for(u=i+TURN+2; uMLclosing; if(dangles == 2) energy += E_MLstem(tt, S1[j-1], S1[i+1], P); else energy += E_MLstem(tt, -1, -1, P); if((d1 <= k) && (d2 <= l)) for(cnt1 = k_min_M[i1u]; cnt1 <= MIN2(k-d1, k_max_M[i1u]); cnt1++) for(cnt2 = l_min_M[i1u][cnt1]; cnt2 <= MIN2(l-d2, l_max_M[i1u][cnt1]); cnt2+=2) if( ((k-d1-cnt1) >= k_min_M1[u1j1]) && ((k-d1-cnt1) <= k_max_M1[u1j1])) if( ((l-d2-cnt2) >= l_min_M1[u1j1][k-d1-cnt1]) && ((l-d2-cnt2) <= l_max_M1[u1j1][k-d1-cnt1])) if(e == (energy + E_M[i1u][cnt1][cnt2/2] + E_M1[u1j1][k-d1-cnt1][(l-d2-cnt2)/2])){ backtrack_m(i+1, u, cnt1, cnt2, structure, vc); backtrack_m1(u+1, j-1, k-d1-cnt1, l-d2-cnt2, structure, vc); return; } } } } vrna_message_error("backtracking failed in c"); } PRIVATE void backtrack_m(unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc){ unsigned int u, ij, seq_length, base_d1, base_d2, d1, d2, maxD1, maxD2; int *my_iindx, *jindx, type, energy, dangles,circ, cnt1, cnt2, cnt3, cnt4; int **l_min_C, **l_max_C,**l_min_M, **l_max_M; int *k_min_C, *k_max_C,*k_min_M, *k_max_M; int ***E_C, ***E_M, *E_C_rem, *E_M_rem; short *S1; unsigned int *referenceBPs1, *referenceBPs2; char *ptype; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; P = vc->params; md = &(P->model_details); matrices = vc->matrices; seq_length = vc->length; S1 = vc->sequence_encoding; circ = md->circ; ptype = vc->ptype; my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; dangles = md->dangles; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_M = matrices->E_M; l_min_M = matrices->l_min_M; l_max_M = matrices->l_max_M; k_min_M = matrices->k_min_M; k_max_M = matrices->k_max_M; E_C_rem = matrices->E_C_rem; E_M_rem = matrices->E_M_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; ij = my_iindx[i]-j; int e = (k == -1) ? E_M_rem[ij] : E_M[ij][k][l/2]; base_d1 = referenceBPs1[ij]; base_d2 = referenceBPs2[ij]; if(k == -1){ /* new_fML = ML(i+1,j)+c */ d1 = base_d1 - referenceBPs1[my_iindx[i+1]-j]; d2 = base_d2 - referenceBPs2[my_iindx[i+1]-j]; if(E_M_rem[my_iindx[i+1]-j] != INF){ if(e == (E_M_rem[my_iindx[i+1]-j] + P->MLbase)){ backtrack_m(i+1,j,-1,-1,structure,vc); return; } } if(E_M[my_iindx[i+1]-j]) for(cnt1 = k_min_M[my_iindx[i+1]-j]; cnt1 <= k_max_M[my_iindx[i+1]-j]; cnt1++) for(cnt2 = l_min_M[my_iindx[i+1]-j][cnt1]; cnt2 <= l_max_M[my_iindx[i+1]-j][cnt1]; cnt2 += 2) if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){ if(e == (E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase)){ backtrack_m(i+1,j,cnt1,cnt2,structure,vc); return; } } /* new_fML = min(ML(i,j-1) + c, new_fML) */ d1 = base_d1 - referenceBPs1[ij+1]; d2 = base_d2 - referenceBPs2[ij+1]; if(E_M_rem[ij+1] != INF){ if(e == (E_M_rem[ij+1] + P->MLbase)){ backtrack_m(i,j-1,-1,-1,structure,vc); return; } } if(E_M[ij+1]) for(cnt1 = k_min_M[ij+1]; cnt1 <= k_max_M[ij+1]; cnt1++) for(cnt2 = l_min_M[ij+1][cnt1]; cnt2 <= l_max_M[ij+1][cnt1]; cnt2 += 2) if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){ if(e == (E_M[ij+1][cnt1][cnt2/2] + P->MLbase)){ backtrack_m(i,j-1,cnt1,cnt2,structure,vc); return; } } /* new_fML = min(new_fML, C(i,j)+b) */ if(E_C_rem[ij] != INF){ type = ptype[jindx[j] + i]; if(dangles == 2) energy = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P); else energy = E_MLstem(type, -1, -1, P); if(e == (E_C_rem[ij] + energy)){ backtrack_c(i,j,-1,-1,structure,vc); return; } } /* modular decomposition -------------------------------*/ for(u = i+1+TURN; u <= j-2-TURN; u++){ int iu, uj; iu = my_iindx[i]-u; uj = my_iindx[u+1]-j; type = ptype[jindx[j] + u + 1]; d1 = base_d1 - referenceBPs1[iu] - referenceBPs1[uj]; d2 = base_d2 - referenceBPs2[iu] - referenceBPs2[uj]; if(dangles == 2) energy = E_MLstem(type, S1[u], (j < seq_length) || circ ? S1[j+1] : -1, P); else energy = E_MLstem(type, -1, -1, P); if(E_M_rem[iu] != INF){ if(E_C[uj]) for(cnt1 = k_min_C[uj]; cnt1 <= k_max_C[uj]; cnt1++) for(cnt2 = l_min_C[uj][cnt1]; cnt2 <= l_max_C[uj][cnt1]; cnt2 += 2) if(e == (E_M_rem[iu] + E_C[uj][cnt1][cnt2/2] + energy)){ backtrack_m(i,u,-1,-1,structure,vc); backtrack_c(u+1,j,cnt1,cnt2,structure, vc); return; } if(E_C_rem[uj] != INF){ if(e == (E_M_rem[iu] + E_C_rem[uj] + energy)){ backtrack_m(i,u,-1,-1,structure,vc); backtrack_c(u+1,j,-1,-1,structure,vc); return; } } } if(E_C_rem[uj] != INF){ if(E_M[iu]) for(cnt1 = k_min_M[iu]; cnt1 <= k_max_M[iu]; cnt1++) for(cnt2 = l_min_M[iu][cnt1]; cnt2 <= l_max_M[iu][cnt1]; cnt2 += 2) if(e == (E_M[iu][cnt1][cnt2/2] + E_C_rem[uj] + energy)){ backtrack_m(i,u,cnt1,cnt2,structure,vc); backtrack_c(u+1,j,-1,-1,structure,vc); return; } } if(!E_M[iu]) continue; if(!E_C[uj]) continue; for(cnt1 = k_min_M[iu]; cnt1 <= k_max_M[iu]; cnt1++) for(cnt2 = l_min_M[iu][cnt1]; cnt2 <= l_max_M[iu][cnt1]; cnt2 += 2) for(cnt3 = k_min_C[uj]; cnt3 <= k_max_C[uj]; cnt3++){ for(cnt4 = l_min_C[uj][cnt3]; cnt4 <= l_max_C[uj][cnt3]; cnt4 += 2) if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)) if(e == (E_M[iu][cnt1][cnt2/2] + E_C[uj][cnt3][cnt4/2] + energy)){ backtrack_m(i, u, cnt1, cnt2, structure, vc); backtrack_c(u+1, j, cnt3, cnt4, structure, vc); return; } } } } /* end if (k == -1) */ else{ d1 = base_d1 - referenceBPs1[my_iindx[i+1]-j]; d2 = base_d2 - referenceBPs2[my_iindx[i+1]-j]; /* new_fML = ML(i+1,j)+c */ if(d1 <= k && d2 <= l) if((k-d1 >= k_min_M[my_iindx[i+1]-j]) && (k-d1 <= k_max_M[my_iindx[i+1]-j])) if((l-d2 >= l_min_M[my_iindx[i+1]-j][k-d1]) && (l-d2 <= l_max_M[my_iindx[i+1]-j][k-d1])){ if(E_M[my_iindx[i+1]-j][k-d1][(l-d2)/2] + P->MLbase == e){ backtrack_m(i+1, j, k-d1, l-d2, structure, vc); return; } } d1 = base_d1 - referenceBPs1[ij+1]; d2 = base_d2 - referenceBPs2[ij+1]; /* new_fML = min(ML(i,j-1) + c, new_fML) */ if(E_M[ij+1]) if(d1 <= k && d2 <= l) if((k-d1 >= k_min_M[ij+1]) && (k-d1 <= k_max_M[ij+1])) if((l-d2 >= l_min_M[ij+1][k-d1]) && (l-d2 <= l_max_M[ij+1][k-d1])) if(E_M[ij+1][k-d1][(l-d2)/2] + P->MLbase == e){ backtrack_m(i, j-1, k-d1, l-d2, structure, vc); return; } /* new_fML = min(new_fML, C(i,j)+b) */ if(E_C[ij]){ type = ptype[jindx[j] + i]; if(dangles == 2) energy = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P); else energy = E_MLstem(type, -1, -1, P); if((k >= k_min_C[ij]) && (k <= k_max_C[ij])) if((l >= l_min_C[ij][k]) && (l <= l_max_C[ij][k])){ if(E_C[ij][k][l/2] + energy == e){ backtrack_c(i, j, k, l, structure, vc); return; } } } /* modular decomposition -------------------------------*/ for(u = i+1+TURN; u <= j-2-TURN; u++){ if(!E_M[my_iindx[i]-u]) continue; if(!E_C[my_iindx[u+1]-j]) continue; type = ptype[jindx[j] + u + 1]; d1 = base_d1 - referenceBPs1[my_iindx[i]-u] - referenceBPs1[my_iindx[u+1]-j]; d2 = base_d2 - referenceBPs2[my_iindx[i]-u] - referenceBPs2[my_iindx[u+1]-j]; if(dangles == 2) energy = E_MLstem(type, S1[u], ((j < seq_length) || circ) ? S1[j+1] : -1, P); else energy = E_MLstem(type, -1, -1, P); if(d1 <= k && d2 <= l) for(cnt1 = k_min_M[my_iindx[i]-u]; cnt1 <= MIN2(k-d1, k_max_M[my_iindx[i]-u]); cnt1++) for(cnt2 = l_min_M[my_iindx[i]-u][cnt1]; cnt2 <= MIN2(l-d2, l_max_M[my_iindx[i]-u][cnt1]); cnt2+=2) if((k-d1-cnt1 >= k_min_C[my_iindx[u+1]-j]) && (k-d1-cnt1 <= k_max_C[my_iindx[u+1]-j])) if((l-d2-cnt2 >= l_min_C[my_iindx[u+1]-j][k-d1-cnt1]) && (l-d2-cnt2 <= l_max_C[my_iindx[u+1]-j][k-d1-cnt1])) if(E_M[my_iindx[i]-u][cnt1][cnt2/2] + E_C[my_iindx[u+1]-j][k-d1-cnt1][(l-d2-cnt2)/2] + energy == e){ backtrack_m(i, u, cnt1, cnt2, structure, vc); backtrack_c(u+1, j, k-d1-cnt1, l-d2-cnt2, structure, vc); return; } } } vrna_message_error("backtracking failed in fML\n"); } PRIVATE void backtrack_m1( unsigned int i, unsigned int j, int k, int l, char *structure, vrna_fold_compound_t *vc){ unsigned int ij, seq_length, d1, d2, *referenceBPs1, *referenceBPs2, maxD1, maxD2; int *my_iindx, *jindx, **l_min_C, **l_max_C,**l_min_M1, **l_max_M1; int *k_min_C, *k_max_C,*k_min_M1, *k_max_M1, cnt1, cnt2; int ***E_C, ***E_M1, *E_C_rem, *E_M1_rem, type, dangles, circ, energy, e_m1; short *S1; char *ptype; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; P = vc->params; md = &(P->model_details); matrices = vc->matrices; seq_length = vc->length; S1 = vc->sequence_encoding; ptype = vc->ptype; circ = md->circ; my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; dangles = md->dangles; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_M1 = matrices->E_M1; l_min_M1 = matrices->l_min_M1; l_max_M1 = matrices->l_max_M1; k_min_M1 = matrices->k_min_M1; k_max_M1 = matrices->k_max_M1; E_C_rem = matrices->E_C_rem; E_M1_rem = matrices->E_M1_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; ij = my_iindx[i]-j; e_m1 = (k == -1) ? E_M1_rem[ij] : E_M1[ij][k][l/2]; type = ptype[jindx[j] + i]; d1 = referenceBPs1[ij] - referenceBPs1[ij+1]; d2 = referenceBPs2[ij] - referenceBPs2[ij+1]; if(dangles == 2) energy = E_MLstem(type, (i > 1) || circ ? S1[i-1] : -1, (j < seq_length) || circ ? S1[j+1] : -1, P); else energy = E_MLstem(type, -1, -1, P); if(k == -1){ if(E_C_rem[ij] != INF){ if(e_m1 == (E_C_rem[ij] + energy)){ backtrack_c(i,j,-1,-1,structure,vc); return; } } if(E_M1_rem[ij+1] != INF){ if(e_m1 == (E_M1_rem[ij+1] + P->MLbase)){ backtrack_m1(i,j-1,-1,-1,structure,vc); return; } } for(cnt1 = k_min_M1[ij+1]; cnt1 <= k_max_M1[ij+1]; cnt1++) for(cnt2 = l_min_M1[ij+1][cnt1]; cnt2 <= l_max_M1[ij+1][cnt1]; cnt2 += 2) if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){ if(e_m1 == (E_M1[ij+1][cnt1][cnt2/2] + P->MLbase)){ backtrack_m1(i,j-1,cnt1,cnt2,structure,vc); return; } } } else{ if(E_C[ij]) if((k >= k_min_C[ij]) && (k <= k_max_C[ij])) if((l >= l_min_C[ij][k]) && (l <= l_max_C[ij][k])) if(E_C[ij][k][l/2] + energy == e_m1){ backtrack_c(i, j, k, l, structure, vc); return; } if(d1 <= k && d2 <= l) if((k-d1 >= k_min_M1[ij+1]) && (k-d1 <= k_max_M1[ij+1])) if((l-d2 >= l_min_M1[ij+1][k-d1]) && (l-d2 <= l_max_M1[ij+1][k-d1])) if(E_M1[ij+1][k-d1][(l-d2)/2] + P->MLbase == e_m1){ backtrack_m1(i, j-1, k-d1, l-d2, structure, vc); return; } } vrna_message_error("backtack failed in m1\n"); } PRIVATE void backtrack_fc( int k, int l, char *structure, vrna_fold_compound_t *vc){ unsigned int d, i, j, seq_length, base_d1, base_d2, d1, d2, maxD1, maxD2; int *my_iindx, *jindx, energy, cnt1, cnt2, cnt3, cnt4, *rtype; short *S1; unsigned int *referenceBPs1, *referenceBPs2; char *sequence, *ptype; int **E_Fc, **E_FcH, **E_FcI, **E_FcM, ***E_C, ***E_M, ***E_M2; int *E_C_rem, *E_M_rem, *E_M2_rem, E_Fc_rem, E_FcH_rem, E_FcI_rem, E_FcM_rem; int **l_min_C, **l_max_C, *k_min_C, *k_max_C; int **l_min_M, **l_max_M, *k_min_M, *k_max_M; int **l_min_M2, **l_max_M2, *k_min_M2, *k_max_M2; int *l_min_FcH, *l_max_FcH, k_min_FcH, k_max_FcH; int *l_min_FcI, *l_max_FcI, k_min_FcI, k_max_FcI; int *l_min_FcM, *l_max_FcM, k_min_FcM, k_max_FcM; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; P = vc->params; md = &(P->model_details); matrices = vc->matrices; sequence = vc->sequence; seq_length = vc->length; S1 = vc->sequence_encoding; ptype = vc->ptype; rtype = &(md->rtype[0]); my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; base_d1 = referenceBPs1[my_iindx[1]-seq_length]; base_d2 = referenceBPs2[my_iindx[1]-seq_length]; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_M = matrices->E_M; l_min_M = matrices->l_min_M; l_max_M = matrices->l_max_M; k_min_M = matrices->k_min_M; k_max_M = matrices->k_max_M; E_M2 = matrices->E_M2; l_min_M2 = matrices->l_min_M2; l_max_M2 = matrices->l_max_M2; k_min_M2 = matrices->k_min_M2; k_max_M2 = matrices->k_max_M2; E_Fc = matrices->E_Fc; E_FcI = matrices->E_FcI; l_min_FcI = matrices->l_min_FcI; l_max_FcI = matrices->l_max_FcI; k_min_FcI = matrices->k_min_FcI; k_max_FcI = matrices->k_max_FcI; E_FcH = matrices->E_FcH; l_min_FcH = matrices->l_min_FcH; l_max_FcH = matrices->l_max_FcH; k_min_FcH = matrices->k_min_FcH; k_max_FcH = matrices->k_max_FcH; E_FcM = matrices->E_FcM; l_min_FcM = matrices->l_min_FcM; l_max_FcM = matrices->l_max_FcM; k_min_FcM = matrices->k_min_FcM; k_max_FcM = matrices->k_max_FcM; E_C_rem = matrices->E_C_rem; E_M_rem = matrices->E_M_rem; E_M2_rem = matrices->E_M2_rem; E_Fc_rem = matrices->E_Fc_rem; E_FcH_rem = matrices->E_FcH_rem; E_FcI_rem = matrices->E_FcI_rem; E_FcM_rem = matrices->E_FcM_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; if(k==-1){ /* check if mfe might be open chain */ if(E_Fc_rem == 0) if((referenceBPs1[my_iindx[1]-seq_length] > maxD1) || (referenceBPs2[my_iindx[1]-seq_length] > maxD2)) return; /* check for hairpin configurations */ if(E_Fc_rem == E_FcH_rem){ for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij; int type, no_close; char loopseq[10]; i = j-d+1; ij = my_iindx[i]-j; u = seq_length-j + i-1; if (u maxD1) || ((cnt2 + d2) > maxD2)) if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + energy)){ backtrack_c(i,j,cnt1,cnt2,structure,vc); return; } } } /* check for interior loop configurations */ if(E_Fc_rem == E_FcI_rem){ for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij, p, q, pq; int type, type_2; i = j-d+1; ij = my_iindx[i]-j; u = seq_length-j + i-1; if (uMAXLOOP) break; qmin = p + TURN + 1; ln_pre = u1 + i + seq_length; if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1; for(q = qmin; q <= seq_length; q++){ unsigned int u2; pq = my_iindx[p]-q; type_2 = rtype[(unsigned int)ptype[jindx[q] + p]]; if (type_2==0) continue; u2 = i-1 + seq_length-q; if (u1+u2>MAXLOOP) continue; energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P); if(E_C_rem[ij] != INF){ if(E_C[pq]) for(cnt1 = k_min_C[pq]; cnt1 <= k_max_C[pq]; cnt1++) for(cnt2 = l_min_C[pq][cnt1]; cnt2 <= l_max_C[pq][cnt1]; cnt2 += 2) if(E_Fc_rem == (E_C_rem[ij] + E_C[pq][cnt1][cnt2/2] + energy)){ backtrack_c(i,j,-1,-1,structure,vc); backtrack_c(p,q,cnt1,cnt2,structure,vc); return; } if(E_C_rem[pq] != INF){ if(E_Fc_rem == (E_C_rem[ij] + E_C_rem[pq] + energy)){ backtrack_c(i,j,-1,-1,structure,vc); backtrack_c(p,q,-1,-1,structure,vc); return; } } } if(E_C_rem[pq] != INF){ if(E_C[ij]) for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2) if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + E_C_rem[pq] + energy)){ backtrack_c(i,j,cnt1,cnt2,structure,vc); backtrack_c(p,q,-1,-1,structure,vc); return; } } if(!(E_C[ij])) continue; if(!(E_C[pq])) continue; /* get distance to reference if closing the interior loop * d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j}) * d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j}) */ d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq]; d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq]; for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2) for(cnt3 = k_min_C[pq]; cnt3 <= k_max_C[pq]; cnt3++) for(cnt4 = l_min_C[pq][cnt3]; cnt4 <= l_max_C[pq][cnt3]; cnt4 += 2) if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)) if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + E_C[pq][cnt3][cnt4/2] + energy)){ backtrack_c(i, j, cnt1, cnt2, structure, vc); backtrack_c(p, q, cnt3, cnt4, structure, vc); return; } } /* end for p */ } /* end for q */ } } /* check for multi loop configurations */ if(E_Fc_rem == E_FcM_rem){ if(seq_length > 2*TURN) for (i=TURN+1; iMLclosing)){ backtrack_m(1,i,-1,-1,structure,vc); backtrack_m2(i+1,cnt1,cnt2,structure,vc); return; } if(E_M2_rem[i+1] != INF){ if(E_Fc_rem == (E_M_rem[my_iindx[1]-i] + E_M2_rem[i+1] + P->MLclosing)){ backtrack_m(1,i,-1,-1,structure,vc); backtrack_m2(i+1,-1,-1,structure,vc); return; } } } if(E_M2_rem[i+1] != INF){ if(E_M[my_iindx[1]-i]) for(cnt1 = k_min_M[my_iindx[1]-i]; cnt1 <= k_max_M[my_iindx[1]-i]; cnt1++) for(cnt2 = l_min_M[my_iindx[1]-i][cnt1]; cnt2 <= l_max_M[my_iindx[1]-i][cnt1]; cnt2 += 2) if(E_Fc_rem == (E_M[my_iindx[1]-i][cnt1][cnt2/2] + E_M2_rem[i+1] + P->MLclosing)){ backtrack_m(1,i,cnt1,cnt2,structure,vc); backtrack_m2(i+1,-1,-1,structure,vc); return; } } if(!(E_M[my_iindx[1]-i])) continue; if(!(E_M2[i+1])) continue; d1 = base_d1 - referenceBPs1[my_iindx[1]-i] - referenceBPs1[my_iindx[i+1]-seq_length]; d2 = base_d2 - referenceBPs2[my_iindx[1]-i] - referenceBPs2[my_iindx[i+1]-seq_length]; for(cnt1 = k_min_M[my_iindx[1]-i]; cnt1 <= k_max_M[my_iindx[1]-i]; cnt1++) for(cnt2 = l_min_M[my_iindx[1]-i][cnt1]; cnt2 <= l_max_M[my_iindx[1]-i][cnt1]; cnt2 += 2) for(cnt3 = k_min_M2[i+1]; cnt3 <= k_max_M2[i+1]; cnt3++) for(cnt4 = l_min_M2[i+1][cnt3]; cnt4 <= l_max_M2[i+1][cnt3]; cnt4 += 2) if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){ if(E_Fc_rem == (E_M[my_iindx[1]-i][cnt1][cnt2/2] + E_M2[i+1][cnt3][cnt4/2] + P->MLclosing)){ backtrack_m(1, i, cnt1, cnt2, structure, vc); backtrack_m2(i+1, cnt3, cnt4, structure, vc); return; } } } } } else{ /* open chain ? */ if(E_Fc[k][l/2] == 0) if((k == referenceBPs1[my_iindx[1]-seq_length]) && (l == referenceBPs2[my_iindx[1]-seq_length])){ return; } if((k >= k_min_FcH) && (k <= k_max_FcH)){ if((l >= l_min_FcH[k]) && (l <= l_max_FcH[k])) if(E_Fc[k][l/2] == E_FcH[k][l/2]){ for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij; int type, no_close; char loopseq[10]; i = j-d+1; ij = my_iindx[i]-j; if (!E_C[ij]) continue; u = seq_length-j + i-1; if (u= d1) && (l >= d2)) if((k-d1 >= k_min_C[ij]) && (k-d1 <= k_max_C[ij])) if((l-d2 >= l_min_C[ij][k-d1]) && (l-d2 <= l_max_C[ij][k-d1])){ if(E_Fc[k][l/2] == E_C[ij][k-d1][(l-d2)/2] + energy){ backtrack_c(i, j, k-d1, l-d2, structure, vc); return; } } } } } if((k >= k_min_FcI) && (k <= k_max_FcI)){ if((l >= l_min_FcI[k]) && (l <= l_max_FcI[k])) if(E_Fc[k][l/2] == E_FcI[k][l/2]){ for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij, p, q, pq; int type, type_2; i = j-d+1; ij = my_iindx[i]-j; if(!E_C[ij]) continue; u = seq_length-j + i-1; if (uMAXLOOP) break; qmin = p + TURN + 1; ln_pre = u1 + i + seq_length; if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1; for(q = qmin; q <= seq_length; q++){ unsigned int u2; pq = my_iindx[p]-q; if(!E_C[pq]) continue; type_2 = rtype[(unsigned int)ptype[jindx[q] + p]]; if (type_2==0) continue; u2 = i-1 + seq_length-q; if (u1+u2>MAXLOOP) continue; /* get distance to reference if closing the interior loop * d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j}) * d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j}) */ d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq]; d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq]; energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P); if((k >= d1) && (l >= d2)) for(cnt1 = k_min_C[ij]; cnt1 <= MIN2(k_max_C[ij], k - d1); cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= MIN2(l_max_C[ij][cnt1], l - d2); cnt2+=2) if((k - d1 - cnt1 >= k_min_C[pq]) && (k - d1 - cnt1 <= k_max_C[pq])) if((l - d2 - cnt2 >= l_min_C[pq][k-d1-cnt1]) && (l - d2 - cnt2 <= l_max_C[pq][k-d1-cnt1])){ if((E_C[ij][cnt1][cnt2/2] + E_C[pq][k-d1-cnt1][(l-d2-cnt2)/2] + energy) == E_Fc[k][l/2]){ backtrack_c(i, j, cnt1, cnt2, structure, vc); backtrack_c(p, q, k - d1 - cnt1, l - d2 - cnt2, structure, vc); return; } } } } } } } if((k >= k_min_FcM) && (k <= k_max_FcM)){ if((l >= l_min_FcM[k]) && (l <= l_max_FcM[k])) if(E_Fc[k][l/2] == E_FcM[k][l/2]){ if(seq_length > 2*TURN) for (i=TURN+1; i= d1) && (l >= d2)) for(cnt1 = k_min_M[my_iindx[1]-i]; cnt1 <= MIN2(k_max_M[my_iindx[1]-i], k-d1); cnt1++) for(cnt2 = l_min_M[my_iindx[1]-i][cnt1]; cnt2 <= MIN2(l_max_M[my_iindx[1]-i][cnt1], l-d2); cnt2+=2) if((k - d1 - cnt1 >= k_min_M2[i+1]) && (k - d1 - cnt1 <= k_max_M2[i+1])) if((l - d2 - cnt2 >= l_min_M2[i+1][k-d1-cnt1]) && (l - d2 - cnt2 <= l_max_M2[i+1][k-d1-cnt1])) if((E_M[my_iindx[1]-i][cnt1][cnt2/2] + E_M2[i+1][k-d1-cnt1][(l-d2-cnt2)/2] + P->MLclosing) == E_FcM[k][l/2]){ backtrack_m(1, i, cnt1, cnt2, structure, vc); backtrack_m2(i+1, k - d1 - cnt1, l - d2 - cnt2, structure, vc); return; } } } } } vrna_message_error("backtack failed in fc\n"); } PRIVATE void backtrack_m2( unsigned int i, int k, int l, char *structure, vrna_fold_compound_t *vc){ unsigned int j, ij, j3, n; unsigned int *referenceBPs1, *referenceBPs2; unsigned int d1, d2, base_d1, base_d2, maxD1, maxD2; int *my_iindx, cnt1, cnt2, cnt3, cnt4; int ***E_M1, ***E_M2, *E_M2_rem, *E_M1_rem, e; int **l_min_M1, **l_max_M1, *k_min_M1, *k_max_M1; vrna_mx_mfe_t *matrices; matrices = vc->matrices; n = vc->length; my_iindx = vc->iindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; E_M1 = matrices->E_M1; l_min_M1 = matrices->l_min_M1; l_max_M1 = matrices->l_max_M1; k_min_M1 = matrices->k_min_M1; k_max_M1 = matrices->k_max_M1; E_M1_rem = matrices->E_M1_rem; E_M2 = matrices->E_M2; E_M2_rem = matrices->E_M2_rem; maxD1 = vc->maxD1; maxD2 = vc->maxD2; base_d1 = referenceBPs1[my_iindx[i]-n]; base_d2 = referenceBPs2[my_iindx[i]-n]; if(k == -1){ e = E_M2_rem[i]; for (j=i+TURN+1; j maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){ if(e == E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-n][cnt3][cnt4/2]){ backtrack_m1(i, j, cnt1, cnt2, structure, vc); backtrack_m1(j+1, n, cnt3, cnt4, structure, vc); return; } } } } } } else{ for(j=i+TURN+1; j= k_min_M1[j3]) && (k - d1 - cnt1 <= k_max_M1[j3])) if((l - d2 - cnt2 >= l_min_M1[j3][k - d1 - cnt1]) && (l - d2 - cnt2 <= l_max_M1[j3][k-d1-cnt1])) if(E_M1[ij][cnt1][cnt2/2] + E_M1[j3][k-d1-cnt1][(l-d2-cnt2)/2] == E_M2[i][k][l/2]){ backtrack_m1(i, j, cnt1, cnt2, structure, vc); backtrack_m1(j+1, n, k-d1-cnt1, l-d2-cnt2, structure, vc); return; } } } vrna_message_error("backtack failed in m2\n"); } PRIVATE void mfe_circ(vrna_fold_compound_t *vc){ unsigned int d, i, j, maxD1, maxD2, seq_length, *referenceBPs1, *referenceBPs2, d1, d2, base_d1, base_d2, *mm1, *mm2, *bpdist; int *my_iindx, *jindx, energy, cnt1, cnt2, cnt3, cnt4, *rtype; short *S1; char *sequence, *ptype; int ***E_C, ***E_M, ***E_M1; int *E_C_rem, *E_M_rem, *E_M1_rem; int **l_min_C, **l_max_C, **l_min_M, **l_max_M, **l_min_M1, **l_max_M1; int *k_min_C, *k_max_C,*k_min_M, *k_max_M,*k_min_M1, *k_max_M1; vrna_param_t *P; vrna_md_t *md; vrna_mx_mfe_t *matrices; P = vc->params; md = &(P->model_details); matrices = vc->matrices; sequence = vc->sequence; seq_length = vc->length; maxD1 = vc->maxD1; maxD2 = vc->maxD2; S1 = vc->sequence_encoding; ptype = vc->ptype; rtype = &(md->rtype[0]); my_iindx = vc->iindx; jindx = vc->jindx; referenceBPs1 = vc->referenceBPs1; referenceBPs2 = vc->referenceBPs2; mm1 = vc->mm1; mm2 = vc->mm2; bpdist = vc->bpdist; E_C = matrices->E_C; l_min_C = matrices->l_min_C; l_max_C = matrices->l_max_C; k_min_C = matrices->k_min_C; k_max_C = matrices->k_max_C; E_M = matrices->E_M; l_min_M = matrices->l_min_M; l_max_M = matrices->l_max_M; k_min_M = matrices->k_min_M; k_max_M = matrices->k_max_M; E_M1 = matrices->E_M1; l_min_M1 = matrices->l_min_M1; l_max_M1 = matrices->l_max_M1; k_min_M1 = matrices->k_min_M1; k_max_M1 = matrices->k_max_M1; E_C_rem = matrices->E_C_rem; E_M_rem = matrices->E_M_rem; E_M1_rem = matrices->E_M1_rem; #ifdef _OPENMP #pragma omp parallel for private(d1,d2,cnt1,cnt2,cnt3,cnt4,j, i) #endif for(i=1; ik_min_M2[i], &matrices->k_max_M2[i], &matrices->l_min_M2[i], &matrices->l_max_M2[i] ); prepareArray( &matrices->E_M2[i], matrices->k_min_M2[i], matrices->k_max_M2[i], matrices->l_min_M2[i], matrices->l_max_M2[i] ); preparePosteriorBoundaries( matrices->k_max_M2[i] - matrices->k_min_M2[i] + 1, matrices->k_min_M2[i], &min_k_real, &max_k_real, &min_l_real, &max_l_real ); /* begin filling of M2 array */ for (j=i+TURN+1; jE_M2_rem[i] = MIN2(matrices->E_M2_rem[i], E_M1_rem[my_iindx[i]-j] + E_M1[my_iindx[j+1]-seq_length][cnt1][cnt2/2] ); if(E_M1_rem[my_iindx[j+1]-seq_length] != INF) matrices->E_M2_rem[i] = MIN2(matrices->E_M2_rem[i], E_M1_rem[my_iindx[i]-j] + E_M1_rem[my_iindx[j+1]-seq_length]); } if(E_M1_rem[my_iindx[j+1]-seq_length] != INF){ if(E_M1[my_iindx[i]-j]) for(cnt1 = k_min_M1[my_iindx[i]-j]; cnt1 <= k_max_M1[my_iindx[i]-j]; cnt1++) for(cnt2 = l_min_M1[my_iindx[i]-j][cnt1]; cnt2 <= l_max_M1[my_iindx[i]-j][cnt1]; cnt2 += 2) matrices->E_M2_rem[i] = MIN2(matrices->E_M2_rem[i], E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1_rem[my_iindx[j+1]-seq_length] ); } if(!E_M1[my_iindx[i]-j]) continue; if(!E_M1[my_iindx[j+1]-seq_length]) continue; d1 = referenceBPs1[my_iindx[i]-seq_length] - referenceBPs1[my_iindx[i]-j] - referenceBPs1[my_iindx[j+1]-seq_length]; d2 = referenceBPs2[my_iindx[i]-seq_length] - referenceBPs2[my_iindx[i]-j] - referenceBPs2[my_iindx[j+1]-seq_length]; for(cnt1 = k_min_M1[my_iindx[i]-j]; cnt1 <= k_max_M1[my_iindx[i]-j]; cnt1++) for(cnt2 = l_min_M1[my_iindx[i]-j][cnt1]; cnt2 <= l_max_M1[my_iindx[i]-j][cnt1]; cnt2+=2){ for(cnt3 = k_min_M1[my_iindx[j+1]-seq_length]; cnt3 <= k_max_M1[my_iindx[j+1]-seq_length]; cnt3++) for(cnt4 = l_min_M1[my_iindx[j+1]-seq_length][cnt3]; cnt4 <= l_max_M1[my_iindx[j+1]-seq_length][cnt3]; cnt4+=2){ if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){ matrices->E_M2[i][cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2( matrices->E_M2[i][cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2], E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-seq_length][cnt3][cnt4/2] ); updatePosteriorBoundaries(cnt1+cnt3+d1, cnt2+cnt4+d2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } else{ matrices->E_M2_rem[i] = MIN2(matrices->E_M2_rem[i], E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-seq_length][cnt3][cnt4/2] ); } } } } /* resize and move memory portions of energy matrix E_M2 */ adjustArrayBoundaries(&matrices->E_M2[i], &matrices->k_min_M2[i], &matrices->k_max_M2[i], &matrices->l_min_M2[i], &matrices->l_max_M2[i], min_k_real, max_k_real, min_l_real, max_l_real ); } /* end for i */ base_d1 = referenceBPs1[my_iindx[1]-seq_length]; base_d2 = referenceBPs2[my_iindx[1]-seq_length]; /* guess memory requirements for E_FcH, E_FcI and E_FcM */ int min_k, max_k, max_l, min_l; int min_k_real, max_k_real, min_k_real_fcH, max_k_real_fcH, min_k_real_fcI, max_k_real_fcI, min_k_real_fcM, max_k_real_fcM; int *min_l_real, *max_l_real, *min_l_real_fcH, *max_l_real_fcH, *min_l_real_fcI, *max_l_real_fcI,*min_l_real_fcM, *max_l_real_fcM; min_k = min_l = 0; max_k = mm1[my_iindx[1] - seq_length] + referenceBPs1[my_iindx[1] - seq_length]; max_l = mm2[my_iindx[1] - seq_length] + referenceBPs2[my_iindx[1] - seq_length]; #ifdef _OPENMP #pragma omp sections { #pragma omp section { #endif prepareBoundaries(min_k, max_k, min_l, max_l, bpdist[my_iindx[1] - seq_length], &matrices->k_min_Fc, &matrices->k_max_Fc, &matrices->l_min_Fc, &matrices->l_max_Fc ); prepareArray( &matrices->E_Fc, matrices->k_min_Fc, matrices->k_max_Fc, matrices->l_min_Fc, matrices->l_max_Fc ); #ifdef _OPENMP } #pragma omp section { #endif prepareBoundaries(min_k, max_k, min_l, max_l, bpdist[my_iindx[1] - seq_length], &matrices->k_min_FcH, &matrices->k_max_FcH, &matrices->l_min_FcH, &matrices->l_max_FcH ); prepareArray( &matrices->E_FcH, matrices->k_min_FcH, matrices->k_max_FcH, matrices->l_min_FcH, matrices->l_max_FcH ); #ifdef _OPENMP } #pragma omp section { #endif prepareBoundaries(min_k, max_k, min_l, max_l, bpdist[my_iindx[1] - seq_length], &matrices->k_min_FcI, &matrices->k_max_FcI, &matrices->l_min_FcI, &matrices->l_max_FcI ); prepareArray( &matrices->E_FcI, matrices->k_min_FcI, matrices->k_max_FcI, matrices->l_min_FcI, matrices->l_max_FcI ); #ifdef _OPENMP } #pragma omp section { #endif prepareBoundaries(min_k, max_k, min_l, max_l, bpdist[my_iindx[1] - seq_length], &matrices->k_min_FcM, &matrices->k_max_FcM, &matrices->l_min_FcM, &matrices->l_max_FcM ); prepareArray( &matrices->E_FcM, matrices->k_min_FcM, matrices->k_max_FcM, matrices->l_min_FcM, matrices->l_max_FcM ); #ifdef _OPENMP } #pragma omp section { #endif preparePosteriorBoundaries( max_k - min_k + 1, min_k, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); #ifdef _OPENMP } #pragma omp section { #endif preparePosteriorBoundaries( max_k - min_k + 1, min_k, &min_k_real_fcH, &max_k_real_fcH, &min_l_real_fcH, &max_l_real_fcH ); #ifdef _OPENMP } #pragma omp section { #endif preparePosteriorBoundaries( max_k - min_k + 1, min_k, &min_k_real_fcI, &max_k_real_fcI, &min_l_real_fcI, &max_l_real_fcI ); #ifdef _OPENMP } #pragma omp section { #endif preparePosteriorBoundaries( max_k - min_k + 1, min_k, &min_k_real_fcM, &max_k_real_fcM, &min_l_real_fcM, &max_l_real_fcM ); #ifdef _OPENMP } } #endif /* begin actual energy calculations */ #ifdef _OPENMP #pragma omp sections private(d, d1,d2,cnt1,cnt2,cnt3,cnt4,j, i, energy) { #pragma omp section { #endif for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij; int type, no_close; char loopseq[10]; i = j-d+1; ij = my_iindx[i]-j; u = seq_length-j + i-1; if (uE_FcH_rem = MIN2(matrices->E_FcH_rem, E_C_rem[ij] + energy); if (!E_C[ij]) continue; for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2){ if(((cnt1 + d1) <= maxD1) && ((cnt2 + d2) <= maxD2)){ matrices->E_FcH[cnt1 + d1][(cnt2+d2)/2] = MIN2( matrices->E_FcH[cnt1 + d1][(cnt2+d2)/2], energy + E_C[ij][cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1 + d1, cnt2 + d2, &min_k_real_fcH, &max_k_real_fcH, &min_l_real_fcH, &max_l_real_fcH ); } else matrices->E_FcH_rem = MIN2(matrices->E_FcH_rem, energy + E_C[ij][cnt1][cnt2/2]); } } /* end of i-j loop */ /* resize and move memory portions of energy matrix E_FcH */ adjustArrayBoundaries(&matrices->E_FcH, &matrices->k_min_FcH, &matrices->k_max_FcH, &matrices->l_min_FcH, &matrices->l_max_FcH, min_k_real_fcH, max_k_real_fcH, min_l_real_fcH, max_l_real_fcH ); #ifdef _OPENMP } #pragma omp section { #endif for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */ for (j = d; j <= seq_length; j++) { unsigned int u, ij, p, q, pq; int type, type_2, no_close; i = j-d+1; ij = my_iindx[i]-j; u = seq_length-j + i-1; if (uMAXLOOP) break; qmin = p + TURN + 1; ln_pre = u1 + i + seq_length; if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1; for(q = qmin; q <= seq_length; q++){ unsigned int u2; pq = my_iindx[p]-q; type_2 = rtype[(unsigned int)ptype[jindx[q] + p]]; if (type_2==0) continue; u2 = i-1 + seq_length-q; if (u1+u2>MAXLOOP) continue; /* get distance to reference if closing the interior loop * d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j}) * d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j}) */ d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq]; d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq]; energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P); if(E_C_rem[pq] != INF) matrices->E_FcI_rem = MIN2(matrices->E_FcI_rem, E_C_rem[ij] + E_C_rem[pq] + energy); if(E_C[pq]) for(cnt1 = k_min_C[pq]; cnt1 <= k_max_C[pq]; cnt1++) for(cnt2 = l_min_C[pq][cnt1]; cnt2 <= l_max_C[pq][cnt1]; cnt2 += 2) matrices->E_FcI_rem = MIN2(matrices->E_FcI_rem, E_C_rem[ij] + E_C[pq][cnt1][cnt2/2] + energy); } } } if(E_C[ij]){ for(p = j+1; p < seq_length ; p++){ unsigned int u1, qmin, ln_pre; u1 = p-j-1; if (u1+i-1>MAXLOOP) break; qmin = p + TURN + 1; ln_pre = u1 + i + seq_length; if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1; for(q = qmin; q <= seq_length; q++){ unsigned int u2; pq = my_iindx[p]-q; type_2 = rtype[(unsigned int)ptype[jindx[q] + p]]; if (type_2==0) continue; u2 = i-1 + seq_length-q; if (u1+u2>MAXLOOP) continue; /* get distance to reference if closing the interior loop * d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j}) * d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j}) */ d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq]; d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq]; energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P); if(E_C_rem[pq] != INF){ for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2) matrices->E_FcI_rem = MIN2(matrices->E_FcI_rem, E_C[ij][cnt1][cnt2/2] + E_C_rem[pq] + energy); } if(E_C[pq]) for(cnt1 = k_min_C[ij]; cnt1 <= k_max_C[ij]; cnt1++) for(cnt2 = l_min_C[ij][cnt1]; cnt2 <= l_max_C[ij][cnt1]; cnt2 += 2) for(cnt3 = k_min_C[pq]; cnt3 <= k_max_C[pq]; cnt3++) for(cnt4 = l_min_C[pq][cnt3]; cnt4 <= l_max_C[pq][cnt3]; cnt4 += 2){ if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){ matrices->E_FcI[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2( matrices->E_FcI[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2], E_C[ij][cnt1][cnt2/2] + E_C[pq][cnt3][cnt4/2] + energy ); updatePosteriorBoundaries(cnt1 + cnt3 + d1, cnt2 + cnt4 + d2, &min_k_real_fcI, &max_k_real_fcI, &min_l_real_fcI, &max_l_real_fcI ); } else{ matrices->E_FcI_rem = MIN2( matrices->E_FcI_rem, E_C[ij][cnt1][cnt2/2] + E_C[pq][cnt3][cnt4/2] + energy ); } } } } } } /* end of i-j loop */ /* resize and move memory portions of energy matrix E_FcI */ adjustArrayBoundaries(&matrices->E_FcI, &matrices->k_min_FcI, &matrices->k_max_FcI, &matrices->l_min_FcI, &matrices->l_max_FcI, min_k_real_fcI, max_k_real_fcI, min_l_real_fcI, max_l_real_fcI ); #ifdef _OPENMP } #pragma omp section { #endif if(seq_length > 2*TURN){ for (i=TURN+1; iE_M2[i+1]) for(cnt1 = matrices->k_min_M2[i+1]; cnt1 <= matrices->k_max_M2[i+1]; cnt1++) for(cnt2 = matrices->l_min_M2[i+1][cnt1]; cnt2 <= matrices->l_max_M2[i+1][cnt1]; cnt2 += 2) matrices->E_FcM_rem = MIN2(matrices->E_FcM_rem, E_M_rem[my_iindx[1]-i] + matrices->E_M2[i+1][cnt1][cnt2/2] + P->MLclosing); if(matrices->E_M2_rem[i+1] != INF) matrices->E_FcM_rem = MIN2(matrices->E_FcM_rem, E_M_rem[my_iindx[1]-i] + matrices->E_M2_rem[i+1] + P->MLclosing); } if(matrices->E_M2_rem[i+1] != INF){ if(E_M[my_iindx[1]-i]) for(cnt1 = k_min_M[my_iindx[1]-i]; cnt1 <= k_max_M[my_iindx[1]-i]; cnt1++) for(cnt2 = l_min_M[my_iindx[1]-i][cnt1]; cnt2 <= l_max_M[my_iindx[1]-i][cnt1]; cnt2 += 2) matrices->E_FcM_rem = MIN2(matrices->E_FcM_rem, E_M[my_iindx[1]-i][cnt1][cnt2/2] + matrices->E_M2_rem[i+1] + P->MLclosing); } if(!E_M[my_iindx[1]-i]) continue; if(!matrices->E_M2[i+1]) continue; for(cnt1 = k_min_M[my_iindx[1]-i]; cnt1 <= k_max_M[my_iindx[1]-i]; cnt1++) for(cnt2 = l_min_M[my_iindx[1]-i][cnt1]; cnt2 <= l_max_M[my_iindx[1]-i][cnt1]; cnt2 += 2) for(cnt3 = matrices->k_min_M2[i+1]; cnt3 <= matrices->k_max_M2[i+1]; cnt3++) for(cnt4 = matrices->l_min_M2[i+1][cnt3]; cnt4 <= matrices->l_max_M2[i+1][cnt3]; cnt4 += 2){ if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){ matrices->E_FcM[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2( matrices->E_FcM[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2], E_M[my_iindx[1]-i][cnt1][cnt2/2] + matrices->E_M2[i+1][cnt3][cnt4/2] + P->MLclosing ); updatePosteriorBoundaries(cnt1 + cnt3 + d1, cnt2 + cnt4 + d2, &min_k_real_fcM, &max_k_real_fcM, &min_l_real_fcM, &max_l_real_fcM ); } else{ matrices->E_FcM_rem = MIN2( matrices->E_FcM_rem, E_M[my_iindx[1]-i][cnt1][cnt2/2] + matrices->E_M2[i+1][cnt3][cnt4/2] + P->MLclosing ); } } } } /* resize and move memory portions of energy matrix E_FcM */ adjustArrayBoundaries(&matrices->E_FcM, &matrices->k_min_FcM, &matrices->k_max_FcM, &matrices->l_min_FcM, &matrices->l_max_FcM, min_k_real_fcM, max_k_real_fcM, min_l_real_fcM, max_l_real_fcM ); #ifdef _OPENMP } } #endif /* compute E_Fc_rem */ matrices->E_Fc_rem = MIN2(matrices->E_FcH_rem, matrices->E_FcI_rem); matrices->E_Fc_rem = MIN2(matrices->E_Fc_rem, matrices->E_FcM_rem); /* add the case were structure is unfolded chain */ if((referenceBPs1[my_iindx[1]-seq_length] > maxD1) || (referenceBPs2[my_iindx[1]-seq_length] > maxD2)) matrices->E_Fc_rem = MIN2(matrices->E_Fc_rem, 0); /* compute all E_Fc */ for(cnt1 = matrices->k_min_FcH; cnt1 <= matrices->k_max_FcH; cnt1++) for(cnt2 = matrices->l_min_FcH[cnt1]; cnt2 <= matrices->l_max_FcH[cnt1]; cnt2 += 2){ matrices->E_Fc[cnt1][cnt2/2] = MIN2(matrices->E_Fc[cnt1][cnt2/2], matrices->E_FcH[cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } for(cnt1 = matrices->k_min_FcI; cnt1 <= matrices->k_max_FcI; cnt1++) for(cnt2 = matrices->l_min_FcI[cnt1]; cnt2 <= matrices->l_max_FcI[cnt1]; cnt2 += 2){ matrices->E_Fc[cnt1][cnt2/2] = MIN2(matrices->E_Fc[cnt1][cnt2/2], matrices->E_FcI[cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } for(cnt1 = matrices->k_min_FcM; cnt1 <= matrices->k_max_FcM; cnt1++) for(cnt2 = matrices->l_min_FcM[cnt1]; cnt2 <= matrices->l_max_FcM[cnt1]; cnt2 += 2){ matrices->E_Fc[cnt1][cnt2/2] = MIN2(matrices->E_Fc[cnt1][cnt2/2], matrices->E_FcM[cnt1][cnt2/2] ); updatePosteriorBoundaries(cnt1, cnt2, &min_k_real, &max_k_real, &min_l_real, &max_l_real ); } /* add the case were structure is unfolded chain */ matrices->E_Fc[referenceBPs1[my_iindx[1]-seq_length]][referenceBPs2[my_iindx[1]-seq_length]/2] = MIN2(matrices->E_Fc[referenceBPs1[my_iindx[1]-seq_length]][referenceBPs2[my_iindx[1]-seq_length]/2], 0); updatePosteriorBoundaries(referenceBPs1[my_iindx[1]-seq_length], referenceBPs2[my_iindx[1]-seq_length], &min_k_real, &max_k_real, &min_l_real, &max_l_real ); adjustArrayBoundaries(&matrices->E_Fc, &matrices->k_min_Fc, &matrices->k_max_Fc, &matrices->l_min_Fc, &matrices->l_max_Fc, min_k_real, max_k_real, min_l_real, max_l_real ); } PRIVATE void adjustArrayBoundaries(int ***array, int *k_min, int *k_max, int **l_min, int **l_max, int k_min_post, int k_max_post, int *l_min_post, int *l_max_post){ int cnt1; int k_diff_pre = k_min_post - *k_min; int mem_size = k_max_post - k_min_post + 1; if(k_min_post < INF){ /* free all the unused memory behind actual data */ for(cnt1 = k_max_post + 1; cnt1 <= *k_max; cnt1++){ (*array)[cnt1] += (*l_min)[cnt1]/2; free((*array)[cnt1]); } /* free unused memory before actual data */ for(cnt1 = *k_min; cnt1 < k_min_post; cnt1++){ (*array)[cnt1] += (*l_min)[cnt1]/2; free((*array)[cnt1]); } /* move data to front and thereby eliminating unused memory in front of actual data */ if(k_diff_pre > 0){ memmove((int **)(*array),((int **)(*array)) + k_diff_pre, sizeof(int *) * mem_size); memmove((int *) (*l_min),((int *) (*l_min)) + k_diff_pre, sizeof(int) * mem_size); memmove((int *) (*l_max),((int *) (*l_max)) + k_diff_pre, sizeof(int) * mem_size); } /* reallocating memory to actual size used */ *array += *k_min; *array = (int **)realloc(*array, sizeof(int *) * mem_size); *array -= k_min_post; *l_min += *k_min; *l_min = (int *)realloc(*l_min, sizeof(int) * mem_size); *l_min -= k_min_post; *l_max += *k_min; *l_max = (int *)realloc(*l_max, sizeof(int) * mem_size); *l_max -= k_min_post; /* adjust l dimension of array */ for(cnt1 = k_min_post; cnt1 <= k_max_post; cnt1++){ if(l_min_post[cnt1] < INF){ /* new memsize */ mem_size = (l_max_post[cnt1] - l_min_post[cnt1] + 1)/2 + 1; /* reshift the pointer */ (*array)[cnt1] += (*l_min)[cnt1]/2; int shift = (l_min_post[cnt1]%2 == (*l_min)[cnt1]%2) ? 0 : 1; /* eliminate unused memory in front of actual data */ unsigned int start = (l_min_post[cnt1] - (*l_min)[cnt1])/2 + shift; if(start > 0) memmove((int *)((*array)[cnt1]), (int *)((*array)[cnt1])+start, sizeof(int) * mem_size); (*array)[cnt1] = (int *) realloc((*array)[cnt1], sizeof(int) * mem_size); (*array)[cnt1] -= l_min_post[cnt1]/2; } else{ /* free according memory */ (*array)[cnt1] += (*l_min)[cnt1]/2; free((*array)[cnt1]); } (*l_min)[cnt1] = l_min_post[cnt1]; (*l_max)[cnt1] = l_max_post[cnt1]; } } else{ /* we have to free all unused memory */ for(cnt1 = *k_min; cnt1 <= *k_max; cnt1++){ (*array)[cnt1] += (*l_min)[cnt1]/2; free((*array)[cnt1]); } (*l_min) += *k_min; (*l_max) += *k_min; free(*l_min); free(*l_max); (*array) += *k_min; free(*array); *array = NULL; } l_min_post += *k_min; l_max_post += *k_min; free(l_min_post); free(l_max_post); *k_min = k_min_post; *k_max = k_max_post; } PRIVATE INLINE void preparePosteriorBoundaries(int size, int shift, int *min_k, int *max_k, int **min_l, int **max_l){ int i; *min_k = INF; *max_k = 0; *min_l = (int *)vrna_alloc(sizeof(int) * size); *max_l = (int *)vrna_alloc(sizeof(int) * size); for(i = 0; i < size; i++){ (*min_l)[i] = INF; (*max_l)[i] = 0; } *min_l -= shift; *max_l -= shift; } PRIVATE INLINE void updatePosteriorBoundaries(int d1, int d2, int *min_k, int *max_k, int **min_l, int **max_l){ (*min_l)[d1] = MIN2((*min_l)[d1], d2); (*max_l)[d1] = MAX2((*max_l)[d1], d2); *min_k = MIN2(*min_k, d1); *max_k = MAX2(*max_k, d1); } INLINE PRIVATE void prepareBoundaries(int min_k_pre, int max_k_pre, int min_l_pre, int max_l_pre, int bpdist, int *min_k, int *max_k, int **min_l, int **max_l){ int cnt; int mem = max_k_pre - min_k_pre + 1; *min_k = min_k_pre; *max_k = max_k_pre; *min_l = (int *) vrna_alloc(sizeof(int) * mem); *max_l = (int *) vrna_alloc(sizeof(int) * mem); *min_l -= min_k_pre; *max_l -= min_k_pre; /* for each k guess the according minimum l*/ for(cnt = min_k_pre; cnt <= max_k_pre; cnt++){ (*min_l)[cnt] = min_l_pre; (*max_l)[cnt] = max_l_pre; while((*min_l)[cnt] + cnt < bpdist) (*min_l)[cnt]++; if((bpdist % 2) != (((*min_l)[cnt] + cnt) % 2)) (*min_l)[cnt]++; } } INLINE PRIVATE void prepareArray(int ***array, int min_k, int max_k, int *min_l, int *max_l){ int i, j, mem; *array = (int **)vrna_alloc(sizeof(int *) * (max_k - min_k + 1)); *array -= min_k; for(i = min_k; i <= max_k; i++){ mem = (max_l[i] - min_l[i] + 1)/2 + 1; (*array)[i] = (int *)vrna_alloc(sizeof(int) * mem); for(j = 0; j < mem; j++) (*array)[i][j] = INF; (*array)[i] -= min_l[i]/2; } } INLINE PRIVATE void prepareArray2(unsigned long ***array, int min_k, int max_k, int *min_l, int *max_l){ int i, mem; *array = (unsigned long **)vrna_alloc(sizeof(unsigned long *) * (max_k - min_k + 1)); *array -= min_k; for(i = min_k; i <= max_k; i++){ mem = (max_l[i] - min_l[i] + 1)/2 + 1; (*array)[i] = (unsigned long *)vrna_alloc(sizeof(unsigned long) * mem); (*array)[i] -= min_l[i]/2; } } /* ################################# # OLD API support # ################################# */ /* crosslink data from vars->compatibility to TwoDfold_vars structure */ PRIVATE INLINE void crosslink(TwoDfold_vars *vars){ vrna_fold_compound_t *c; vrna_mx_mfe_t *m; c = vars->compatibility; m = c->matrices; vars->sequence = c->sequence; vars->seq_length = c->length; vars->reference_pt1 = c->reference_pt1; vars->reference_pt2 = c->reference_pt2; vars->referenceBPs1 = c->referenceBPs1; vars->referenceBPs2 = c->referenceBPs2; vars->bpdist = c->bpdist; vars->do_backtrack = 1; vars->dangles = c->params->model_details.dangles; vars->circ = c->params->model_details.circ; vars->temperature = c->params->model_details.temperature; vars->ptype = c->ptype_pf_compat; vars->P = c->params; vars->S = c->sequence_encoding2; vars->S1 = c->sequence_encoding; vars->my_iindx = c->iindx; vars->mm1 = c->mm1; vars->mm2 = c->mm2; vars->maxD1 = c->maxD1; vars->maxD2 = c->maxD2; vars->E_C = m->E_C; vars->l_min_values = m->l_min_C; vars->l_max_values = m->l_max_C; vars->k_min_values = m->k_min_C; vars->k_max_values = m->k_max_C; vars->E_F5 = m->E_F5; vars->l_min_values_f = m->l_min_F5; vars->l_max_values_f = m->l_max_F5; vars->k_min_values_f = m->k_min_F5; vars->k_max_values_f = m->k_max_F5; vars->E_F3 = m->E_F3; vars->l_min_values_f3 = m->l_min_F3; vars->l_max_values_f3 = m->l_max_F3; vars->k_min_values_f3 = m->k_min_F3; vars->k_max_values_f3 = m->k_max_F3; vars->E_M = m->E_M; vars->l_min_values_m = m->l_min_M; vars->l_max_values_m = m->l_max_M; vars->k_min_values_m = m->k_min_M; vars->k_max_values_m = m->k_max_M; vars->E_M1 = m->E_M1; vars->l_min_values_m1 = m->l_min_M1; vars->l_max_values_m1 = m->l_max_M1; vars->k_min_values_m1 = m->k_min_M1; vars->k_max_values_m1 = m->k_max_M1; #ifdef COUNT_STATES vars->N_C = m->N_C; vars->N_F5 = m->N_F5; vars->N_M = m->N_M; vars->N_M1 = m->N_M1; #endif vars->E_M2_rem = m->E_M2_rem; vars->E_M2 = m->E_M2; vars->l_min_values_m2 = m->l_min_M2; vars->l_max_values_m2 = m->l_max_M2; vars->k_min_values_m2 = m->k_min_M2; vars->k_max_values_m2 = m->k_max_M2; vars->E_Fc = m->E_Fc; vars->E_FcH = m->E_FcH; vars->E_FcI = m->E_FcI; vars->E_FcM = m->E_FcM; vars->E_Fc_rem = m->E_Fc_rem; vars->E_FcH_rem = m->E_FcH_rem; vars->E_FcI_rem = m->E_FcI_rem; vars->E_FcM_rem = m->E_FcM_rem; vars->E_C_rem = m->E_C_rem; vars->E_M_rem = m->E_M_rem; vars->E_M1_rem = m->E_M1_rem; vars->E_F5_rem = m->E_F5_rem; } PUBLIC TwoDfold_vars * get_TwoDfold_variables( const char *seq, const char *structure1, const char *structure2, int circ){ vrna_md_t md; TwoDfold_vars *vars; vrna_fold_compound_t *c; vrna_mx_mfe_t *m; set_model_details(&md); md.circ = circ; vars = (TwoDfold_vars *)vrna_alloc(sizeof(TwoDfold_vars)); vars->compatibility = vrna_fold_compound_TwoD(seq, structure1, structure2, &md, VRNA_OPTION_MFE); crosslink(vars); return vars; } PUBLIC char * TwoDfold_backtrack_f5(unsigned int j, int k, int l, TwoDfold_vars *vars){ return vrna_backtrack5_TwoD(vars->compatibility, k, l, j); } PUBLIC void destroy_TwoDfold_variables(TwoDfold_vars *vars){ if(vars == NULL) return; vrna_fold_compound_free(vars->compatibility); free(vars); } PUBLIC vrna_sol_TwoD_t * TwoDfoldList( TwoDfold_vars *vars, int distance1, int distance2){ vrna_sol_TwoD_t *sol; sol = vrna_mfe_TwoD(vars->compatibility, distance1, distance2); crosslink(vars); return sol; } PUBLIC void update_TwoDfold_params(TwoDfold_vars *vars){ vrna_md_t md; set_model_details(&md); free(vars->compatibility->params); vars->compatibility->params = vrna_params(&md); crosslink(vars); }