binaries/src/ViennaRNA/lib/2Dfold.c

   1 /*
   2       minimum free energy
   3       RNA secondary structure with
   4       basepair distance d_1 to reference structure 1 and distance d_2 to reference structure 2
   5
   6 */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <math.h>
  11 #include <ctype.h>
  12 #include <string.h>
  13 #include "utils.h"
  14 #include "energy_par.h"
  15 #include "fold_vars.h"
  16 #include "fold.h"
  17 #include "pair_mat.h"
  18 #include "loop_energies.h"
  19 #include "mm.h"
  20 #include "params.h"
  21 #ifdef _OPENMP
  22 #include <omp.h>
  23 #endif
  24 #include "2Dfold.h"
  25
  26 /*
  27 #################################
  28 # GLOBAL VARIABLES              #
  29 #################################
  30 */
  31 int compute_2Dfold_F3 = 0;
  32
  33 /*
  34 #################################
  35 # PRIVATE VARIABLES             #
  36 #################################
  37 */
  38
  39 /*
  40 #################################
  41 # PRIVATE FUNCTION DECLARATIONS #
  42 #################################
  43 */
  44 PRIVATE void  mfe_linear(TwoDfold_vars *vars);
  45 PRIVATE void  mfe_circ(TwoDfold_vars *vars);
  46
  47 PRIVATE void  initialize_TwoDfold_vars(TwoDfold_vars *vars);
  48 PUBLIC  void  update_TwoDfold_params(TwoDfold_vars *vars);
  49 PRIVATE void  make_ptypes(TwoDfold_vars *vars);
  50
  51 PRIVATE void  backtrack_f5(unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars);
  52 PRIVATE void  backtrack_c(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars);
  53 PRIVATE void  backtrack_m(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars);
  54 PRIVATE void  backtrack_m1(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars);
  55 PRIVATE void  backtrack_fc(int k, int l, char *structure, TwoDfold_vars *vars);
  56 PRIVATE void  backtrack_m2(unsigned int i, int k, int l, char *structure, TwoDfold_vars *vars);
  57
  58 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);
  59 INLINE  PRIVATE void  preparePosteriorBoundaries(int size, int shift, int *min_k, int *max_k, int **min_l, int **max_l);
  60 INLINE  PRIVATE void  updatePosteriorBoundaries(int d1, int d2, int *min_k, int *max_k, int **min_l, int **max_l);
  61 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);
  62 INLINE  PRIVATE void  prepareArray(int ***array, int min_k, int max_k, int *min_l, int *max_l);
  63 INLINE  PRIVATE void  prepareArray2(unsigned long ***array, int min_k, int max_k, int *min_l, int *max_l);
  64
  65
  66
  67 /*
  68 #################################
  69 # BEGIN OF FUNCTION DEFINITIONS #
  70 #################################
  71 */
  72
  73 PUBLIC TwoDfold_vars *get_TwoDfold_variables(const char *seq, const char *structure1, const char *structure2, int circ){
  74   unsigned int size, length, i;
  75   int *index;
  76   TwoDfold_vars *vars;
  77   length = strlen(seq);
  78   vars = (TwoDfold_vars *)malloc(sizeof(TwoDfold_vars));
  79   vars->sequence     = (char *)space(length + 1);
  80   strcpy(vars->sequence, seq);
  81   vars->seq_length   = length;
  82   if(vars->seq_length < 1) nrerror("get_TwoDfold_variables: sequence must be longer than 0");
  83   size                        = ((length + 1) * (length + 2)/2);
  84
  85   vars->reference_pt1   = make_pair_table(structure1);
  86   vars->reference_pt2   = make_pair_table(structure2);
  87   vars->referenceBPs1   = make_referenceBP_array(vars->reference_pt1, TURN);
  88   vars->referenceBPs2   = make_referenceBP_array(vars->reference_pt2, TURN);
  89   vars->bpdist          = compute_BPdifferences(vars->reference_pt1, vars->reference_pt2, TURN);
  90   vars->do_backtrack    = 1;
  91   vars->dangles         = dangles;
  92   vars->circ            = circ;
  93   vars->temperature     = temperature;
  94   vars->ptype           = space(sizeof(char) * size);
  95   vars->P               = NULL;
  96   vars->S               = NULL;
  97   vars->S1              = NULL;
  98   vars->my_iindx        = get_iindx(length);
  99   index                 = vars->my_iindx;
 100   /* compute maximum matching with reference structure 1 disallowed */
 101   vars->mm1          = maximumMatchingConstraint(vars->sequence, vars->reference_pt1);
 102   /* compute maximum matching with reference structure 2 disallowed */
 103   vars->mm2          = maximumMatchingConstraint(vars->sequence, vars->reference_pt2);
 104
 105   vars->maxD1        = vars->mm1[index[1]-length] + vars->referenceBPs1[index[1]-length];
 106   vars->maxD2        = vars->mm2[index[1]-length] + vars->referenceBPs2[index[1]-length];
 107
 108   /* allocate memory for the energy matrices and min-/max-index helper arrays */
 109   vars->E_C              = (int ***) space(sizeof(int **)  * size);
 110   vars->l_min_values     = (int **)  space(sizeof(int *)   * size);
 111   vars->l_max_values     = (int **)  space(sizeof(int *)   * size);
 112   vars->k_min_values     = (int *)   space(sizeof(int)     * size);
 113   vars->k_max_values     = (int *)   space(sizeof(int)     * size);
 114
 115   vars->E_F5             = (int ***) space(sizeof(int **)  * (length + 1));
 116   vars->l_min_values_f   = (int **)  space(sizeof(int *)   * (length + 1));
 117   vars->l_max_values_f   = (int **)  space(sizeof(int *)   * (length + 1));
 118   vars->k_min_values_f   = (int *)   space(sizeof(int)     * (length + 1));
 119   vars->k_max_values_f   = (int *)   space(sizeof(int)     * (length + 1));
 120
 121   if(compute_2Dfold_F3){
 122     vars->E_F3             = (int ***) space(sizeof(int **)  * (length + 1));
 123     vars->l_min_values_f3  = (int **)  space(sizeof(int *)   * (length + 1));
 124     vars->l_max_values_f3  = (int **)  space(sizeof(int *)   * (length + 1));
 125     vars->k_min_values_f3  = (int *)   space(sizeof(int)     * (length + 1));
 126     vars->k_max_values_f3  = (int *)   space(sizeof(int)     * (length + 1));
 127   }
 128   else vars->E_F3 = NULL;
 129
 130   vars->E_M              = (int ***) space(sizeof(int **)  * size);
 131   vars->l_min_values_m   = (int **)  space(sizeof(int *)   * size);
 132   vars->l_max_values_m   = (int **)  space(sizeof(int *)   * size);
 133   vars->k_min_values_m   = (int *)   space(sizeof(int)     * size);
 134   vars->k_max_values_m   = (int *)   space(sizeof(int)     * size);
 135
 136   vars->E_M1             = (int ***) space(sizeof(int **)  * size);
 137   vars->l_min_values_m1  = (int **)  space(sizeof(int *)   * size);
 138   vars->l_max_values_m1  = (int **)  space(sizeof(int *)   * size);
 139   vars->k_min_values_m1  = (int *)   space(sizeof(int)     * size);
 140   vars->k_max_values_m1  = (int *)   space(sizeof(int)     * size);
 141
 142 #ifdef COUNT_STATES
 143   vars->N_C              = (unsigned long ***) space(sizeof(unsigned long **)  * size);
 144   vars->N_F5             = (unsigned long ***) space(sizeof(unsigned long **)  * (length + 1));
 145   vars->N_M              = (unsigned long ***) space(sizeof(unsigned long **)  * size);
 146   vars->N_M1             = (unsigned long ***) space(sizeof(unsigned long **)  * size);
 147 #endif
 148
 149
 150   if(circ){
 151     vars->E_M2_rem       = (int *)   space(sizeof(int)     * (length + 1));
 152     vars->E_M2            = (int ***) space(sizeof(int **)  * (length + 1));
 153     vars->l_min_values_m2 = (int **)  space(sizeof(int *)   * (length + 1));
 154     vars->l_max_values_m2 = (int **)  space(sizeof(int *)   * (length + 1));
 155     vars->k_min_values_m2 = (int *)   space(sizeof(int)     * (length + 1));
 156     vars->k_max_values_m2 = (int *)   space(sizeof(int)     * (length + 1));
 157   }
 158   else{
 159     vars->E_M2_rem       = NULL;
 160     vars->E_M2            = NULL;
 161     vars->l_min_values_m2 = NULL;
 162     vars->l_max_values_m2 = NULL;
 163     vars->k_min_values_m2 = NULL;
 164     vars->k_max_values_m2 = NULL;
 165   }
 166
 167   vars->E_Fc              = NULL;
 168   vars->E_FcH             = NULL;
 169   vars->E_FcI             = NULL;
 170   vars->E_FcM             = NULL;
 171
 172   vars->E_Fc_rem         = INF;
 173   vars->E_FcH_rem        = INF;
 174   vars->E_FcI_rem        = INF;
 175   vars->E_FcM_rem        = INF;
 176
 177   vars->E_C_rem          = (int *) space(sizeof(int) * size);
 178   vars->E_M_rem          = (int *) space(sizeof(int) * size);
 179   vars->E_M1_rem         = (int *) space(sizeof(int) * size);
 180   vars->E_F5_rem         = (int *) space(sizeof(int) * (length+1));
 181   /* init rest arrays */
 182   for(i=0;i<size;i++){
 183     vars->E_C_rem[i] = vars->E_M_rem[i] = vars->E_M1_rem[i] = INF;
 184   }
 185   for(i=0;i<=length;i++)
 186     vars->E_F5_rem[i] = INF;
 187   if(vars->E_M2_rem)
 188     for(i=0;i<=length;i++)
 189       vars->E_M2_rem[i] = INF;
 190
 191   return vars;
 192 }
 193
 194 PUBLIC void destroy_TwoDfold_variables(TwoDfold_vars *vars){
 195   unsigned int i, j, ij;
 196   int cnt1;
 197   if(vars == NULL) return;
 198
 199   free(vars->E_C_rem);
 200   free(vars->E_M_rem);
 201   free(vars->E_M1_rem);
 202   free(vars->E_F5_rem);
 203   if(vars->E_M2_rem) free(vars->E_M2_rem);
 204
 205 #ifdef _OPENMP
 206   #pragma omp sections private(i,j,ij,cnt1)
 207   {
 208
 209   #pragma omp section
 210   {
 211 #endif
 212
 213 #ifdef COUNT_STATES
 214   if(vars->N_C != NULL){
 215     for(i = 1; i < vars->seq_length; i++){
 216       for(j = i; j <= vars->seq_length; j++){
 217         ij = vars->my_iindx[i] - j;
 218         if(!vars->N_C[ij]) continue;
 219         for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++)
 220           if(vars->l_min_values[ij][cnt1] < INF){
 221             vars->N_C[ij][cnt1] += vars->l_min_values[ij][cnt1]/2;
 222             free(vars->N_C[ij][cnt1]);
 223           }
 224         if(vars->k_min_values[ij] < INF){
 225           vars->N_C[ij] += vars->k_min_values[ij];
 226           free(vars->N_C[ij]);
 227         }
 228       }
 229     }
 230     free(vars->N_C);
 231   }
 232 #endif
 233
 234   if(vars->E_C != NULL){
 235     for(i = 1; i < vars->seq_length; i++){
 236       for(j = i; j <= vars->seq_length; j++){
 237         ij = vars->my_iindx[i] - j;
 238         if(!vars->E_C[ij]) continue;
 239         for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++)
 240           if(vars->l_min_values[ij][cnt1] < INF){
 241             vars->E_C[ij][cnt1] += vars->l_min_values[ij][cnt1]/2;
 242             free(vars->E_C[ij][cnt1]);
 243           }
 244         if(vars->k_min_values[ij] < INF){
 245           vars->E_C[ij] += vars->k_min_values[ij];
 246           free(vars->E_C[ij]);
 247           vars->l_min_values[ij] += vars->k_min_values[ij];
 248           vars->l_max_values[ij] += vars->k_min_values[ij];
 249           free(vars->l_min_values[ij]);
 250           free(vars->l_max_values[ij]);
 251         }
 252       }
 253     }
 254     free(vars->E_C);
 255     free(vars->l_min_values);
 256     free(vars->l_max_values);
 257     free(vars->k_min_values);
 258     free(vars->k_max_values);
 259   }
 260 #ifdef _OPENMP
 261   }
 262   #pragma omp section
 263   {
 264 #endif
 265
 266 #ifdef COUNT_STATES
 267   if(vars->N_M != NULL){
 268     for(i = 1; i < vars->seq_length; i++){
 269       for(j = i; j <= vars->seq_length; j++){
 270         ij = vars->my_iindx[i] - j;
 271         if(!vars->N_M[ij]) continue;
 272         for(cnt1 = vars->k_min_values_m[ij]; cnt1 <= vars->k_max_values_m[ij]; cnt1++)
 273           if(vars->l_min_values_m[ij][cnt1] < INF){
 274             vars->N_M[ij][cnt1] += vars->l_min_values_m[ij][cnt1]/2;
 275             free(vars->N_M[ij][cnt1]);
 276           }
 277         if(vars->k_min_values_m[ij] < INF){
 278           vars->N_M[ij] += vars->k_min_values_m[ij];
 279           free(vars->N_M[ij]);
 280         }
 281       }
 282     }
 283     free(vars->N_M);
 284   }
 285 #endif
 286
 287   if(vars->E_M != NULL){
 288     for(i = 1; i < vars->seq_length; i++){
 289       for(j = i; j <= vars->seq_length; j++){
 290         ij = vars->my_iindx[i] - j;
 291         if(!vars->E_M[ij]) continue;
 292         for(cnt1 = vars->k_min_values_m[ij]; cnt1 <= vars->k_max_values_m[ij]; cnt1++)
 293           if(vars->l_min_values_m[ij][cnt1] < INF){
 294             vars->E_M[ij][cnt1] += vars->l_min_values_m[ij][cnt1]/2;
 295             free(vars->E_M[ij][cnt1]);
 296           }
 297         if(vars->k_min_values_m[ij] < INF){
 298           vars->E_M[ij] += vars->k_min_values_m[ij];
 299           free(vars->E_M[ij]);
 300           vars->l_min_values_m[ij] += vars->k_min_values_m[ij];
 301           vars->l_max_values_m[ij] += vars->k_min_values_m[ij];
 302           free(vars->l_min_values_m[ij]);
 303           free(vars->l_max_values_m[ij]);
 304         }
 305       }
 306     }
 307     free(vars->E_M);
 308     free(vars->l_min_values_m);
 309     free(vars->l_max_values_m);
 310     free(vars->k_min_values_m);
 311     free(vars->k_max_values_m);
 312   }
 313 #ifdef _OPENMP
 314   }
 315   #pragma omp section
 316   {
 317 #endif
 318
 319 #ifdef COUNT_STATES
 320   if(vars->N_M1 != NULL){
 321     for(i = 1; i < vars->seq_length; i++){
 322       for(j = i; j <= vars->seq_length; j++){
 323         ij = vars->my_iindx[i] - j;
 324         if(!vars->N_M1[ij]) continue;
 325         for(cnt1 = vars->k_min_values_m1[ij]; cnt1 <= vars->k_max_values_m1[ij]; cnt1++)
 326           if(vars->l_min_values_m1[ij][cnt1] < INF){
 327             vars->N_M1[ij][cnt1] += vars->l_min_values_m1[ij][cnt1]/2;
 328             free(vars->N_M1[ij][cnt1]);
 329           }
 330         if(vars->k_min_values_m1[ij] < INF){
 331           vars->N_M1[ij] += vars->k_min_values_m1[ij];
 332           free(vars->N_M1[ij]);
 333         }
 334       }
 335     }
 336     free(vars->N_M1);
 337   }
 338 #endif
 339
 340   if(vars->E_M1 != NULL){
 341     for(i = 1; i < vars->seq_length; i++){
 342       for(j = i; j <= vars->seq_length; j++){
 343         ij = vars->my_iindx[i] - j;
 344         if(!vars->E_M1[ij]) continue;
 345         for(cnt1 = vars->k_min_values_m1[ij]; cnt1 <= vars->k_max_values_m1[ij]; cnt1++)
 346           if(vars->l_min_values_m1[ij][cnt1] < INF){
 347             vars->E_M1[ij][cnt1] += vars->l_min_values_m1[ij][cnt1]/2;
 348             free(vars->E_M1[ij][cnt1]);
 349           }
 350         if(vars->k_min_values_m1[ij] < INF){
 351           vars->E_M1[ij] += vars->k_min_values_m1[ij];
 352           free(vars->E_M1[ij]);
 353           vars->l_min_values_m1[ij] += vars->k_min_values_m1[ij];
 354           vars->l_max_values_m1[ij] += vars->k_min_values_m1[ij];
 355           free(vars->l_min_values_m1[ij]);
 356           free(vars->l_max_values_m1[ij]);
 357         }
 358       }
 359     }
 360     free(vars->E_M1);
 361     free(vars->l_min_values_m1);
 362     free(vars->l_max_values_m1);
 363     free(vars->k_min_values_m1);
 364     free(vars->k_max_values_m1);
 365   }
 366 #ifdef _OPENMP
 367   }
 368   #pragma omp section
 369   {
 370 #endif
 371   if(vars->E_M2 != NULL){
 372     for(i = 1; i < vars->seq_length-TURN-1; i++){
 373       if(!vars->E_M2[i]) continue;
 374       for(cnt1 = vars->k_min_values_m2[i]; cnt1 <= vars->k_max_values_m2[i]; cnt1++)
 375         if(vars->l_min_values_m2[i][cnt1] < INF){
 376           vars->E_M2[i][cnt1] += vars->l_min_values_m2[i][cnt1]/2;
 377           free(vars->E_M2[i][cnt1]);
 378         }
 379       if(vars->k_min_values_m2[i] < INF){
 380         vars->E_M2[i] += vars->k_min_values_m2[i];
 381         free(vars->E_M2[i]);
 382         vars->l_min_values_m2[i] += vars->k_min_values_m2[i];
 383         vars->l_max_values_m2[i] += vars->k_min_values_m2[i];
 384         free(vars->l_min_values_m2[i]);
 385         free(vars->l_max_values_m2[i]);
 386       }
 387     }
 388     free(vars->E_M2);
 389     free(vars->l_min_values_m2);
 390     free(vars->l_max_values_m2);
 391     free(vars->k_min_values_m2);
 392     free(vars->k_max_values_m2);
 393   }
 394 #ifdef _OPENMP
 395   }
 396   #pragma omp section
 397   {
 398 #endif
 399
 400 #ifdef COUNT_STATES
 401   if(vars->N_F5 != NULL){
 402     for(i = 1; i <= vars->seq_length; i++){
 403       if(!vars->N_F5[i]) continue;
 404       for(cnt1 = vars->k_min_values_f[i]; cnt1 <= vars->k_max_values_f[i]; cnt1++)
 405         if(vars->l_min_values_f[i][cnt1] < INF){
 406           vars->N_F5[i][cnt1] += vars->l_min_values_f[i][cnt1]/2;
 407           free(vars->N_F5[i][cnt1]);
 408         }
 409       if(vars->k_min_values_f[i] < INF){
 410         vars->N_F5[i] += vars->k_min_values_f[i];
 411         free(vars->N_F5[i]);
 412       }
 413     }
 414     free(vars->N_F5);
 415   }
 416 #endif
 417
 418   if(vars->E_F5 != NULL){
 419     for(i = 1; i <= vars->seq_length; i++){
 420       if(!vars->E_F5[i]) continue;
 421       for(cnt1 = vars->k_min_values_f[i]; cnt1 <= vars->k_max_values_f[i]; cnt1++)
 422         if(vars->l_min_values_f[i][cnt1] < INF){
 423           vars->E_F5[i][cnt1] += vars->l_min_values_f[i][cnt1]/2;
 424           free(vars->E_F5[i][cnt1]);
 425         }
 426       if(vars->k_min_values_f[i] < INF){
 427         vars->E_F5[i] += vars->k_min_values_f[i];
 428         free(vars->E_F5[i]);
 429         vars->l_min_values_f[i] += vars->k_min_values_f[i];
 430         vars->l_max_values_f[i] += vars->k_min_values_f[i];
 431         free(vars->l_min_values_f[i]);
 432         free(vars->l_max_values_f[i]);
 433       }
 434     }
 435     free(vars->E_F5);
 436     free(vars->l_min_values_f);
 437     free(vars->l_max_values_f);
 438     free(vars->k_min_values_f);
 439     free(vars->k_max_values_f);
 440   }
 441
 442   if(vars->E_F3 != NULL){
 443     for(i = 1; i <= vars->seq_length; i++){
 444       if(!vars->E_F3[i]) continue;
 445       for(cnt1 = vars->k_min_values_f3[i]; cnt1 <= vars->k_max_values_f3[i]; cnt1++)
 446         if(vars->l_min_values_f3[i][cnt1] < INF){
 447           vars->E_F3[i][cnt1] += vars->l_min_values_f3[i][cnt1]/2;
 448           free(vars->E_F3[i][cnt1]);
 449         }
 450       if(vars->k_min_values_f3[i] < INF){
 451         vars->E_F3[i] += vars->k_min_values_f3[i];
 452         free(vars->E_F3[i]);
 453         vars->l_min_values_f3[i] += vars->k_min_values_f3[i];
 454         vars->l_max_values_f3[i] += vars->k_min_values_f3[i];
 455         free(vars->l_min_values_f3[i]);
 456         free(vars->l_max_values_f3[i]);
 457       }
 458     }
 459     free(vars->E_F3);
 460     free(vars->l_min_values_f3);
 461     free(vars->l_max_values_f3);
 462     free(vars->k_min_values_f3);
 463     free(vars->k_max_values_f3);
 464   }
 465
 466 #ifdef _OPENMP
 467   }
 468   #pragma omp section
 469   {
 470 #endif
 471   if(vars->E_Fc != NULL){
 472     for(cnt1 = vars->k_min_values_fc; cnt1 <= vars->k_max_values_fc; cnt1++)
 473       if(vars->l_min_values_fc[cnt1] < INF){
 474         vars->E_Fc[cnt1] += vars->l_min_values_fc[cnt1]/2;
 475         free(vars->E_Fc[cnt1]);
 476       }
 477     if(vars->k_min_values_fc < INF){
 478       vars->E_Fc += vars->k_min_values_fc;
 479       free(vars->E_Fc);
 480       vars->l_min_values_fc += vars->k_min_values_fc;
 481       vars->l_max_values_fc += vars->k_min_values_fc;
 482       free(vars->l_min_values_fc);
 483       free(vars->l_max_values_fc);
 484     }
 485   }
 486 #ifdef _OPENMP
 487   }
 488   #pragma omp section
 489   {
 490 #endif
 491   if(vars->E_FcI != NULL){
 492     for(cnt1 = vars->k_min_values_fcI; cnt1 <= vars->k_max_values_fcI; cnt1++)
 493       if(vars->l_min_values_fcI[cnt1] < INF){
 494         vars->E_FcI[cnt1] += vars->l_min_values_fcI[cnt1]/2;
 495         free(vars->E_FcI[cnt1]);
 496       }
 497     if(vars->k_min_values_fcI < INF){
 498       vars->E_FcI += vars->k_min_values_fcI;
 499       free(vars->E_FcI);
 500       vars->l_min_values_fcI += vars->k_min_values_fcI;
 501       vars->l_max_values_fcI += vars->k_min_values_fcI;
 502       free(vars->l_min_values_fcI);
 503       free(vars->l_max_values_fcI);
 504     }
 505   }
 506 #ifdef _OPENMP
 507   }
 508   #pragma omp section
 509   {
 510 #endif
 511   if(vars->E_FcH != NULL){
 512     for(cnt1 = vars->k_min_values_fcH; cnt1 <= vars->k_max_values_fcH; cnt1++)
 513       if(vars->l_min_values_fcH[cnt1] < INF){
 514         vars->E_FcH[cnt1] += vars->l_min_values_fcH[cnt1]/2;
 515         free(vars->E_FcH[cnt1]);
 516       }
 517     if(vars->k_min_values_fcH < INF){
 518       vars->E_FcH += vars->k_min_values_fcH;
 519       free(vars->E_FcH);
 520       vars->l_min_values_fcH += vars->k_min_values_fcH;
 521       vars->l_max_values_fcH += vars->k_min_values_fcH;
 522       free(vars->l_min_values_fcH);
 523       free(vars->l_max_values_fcH);
 524     }
 525   }
 526 #ifdef _OPENMP
 527   }
 528   #pragma omp section
 529   {
 530 #endif
 531   if(vars->E_FcM != NULL){
 532     for(cnt1 = vars->k_min_values_fcM; cnt1 <= vars->k_max_values_fcM; cnt1++)
 533       if(vars->l_min_values_fcM[cnt1] < INF){
 534         vars->E_FcM[cnt1] += vars->l_min_values_fcM[cnt1]/2;
 535         free(vars->E_FcM[cnt1]);
 536       }
 537     if(vars->k_min_values_fcM < INF){
 538       vars->E_FcM += vars->k_min_values_fcM;
 539       free(vars->E_FcM);
 540       vars->l_min_values_fcM += vars->k_min_values_fcM;
 541       vars->l_max_values_fcM += vars->k_min_values_fcM;
 542       free(vars->l_min_values_fcM);
 543       free(vars->l_max_values_fcM);
 544     }
 545   }
 546
 547 #ifdef _OPENMP
 548   }
 549   #pragma omp section
 550   {
 551 #endif
 552
 553   if(vars->P != NULL)             free(vars->P);
 554   if(vars->sequence != NULL)      free(vars->sequence);
 555   if(vars->reference_pt1 != NULL) free(vars->reference_pt1);
 556   if(vars->reference_pt2 != NULL) free(vars->reference_pt2);
 557   if(vars->referenceBPs1 != NULL) free(vars->referenceBPs1);
 558   if(vars->referenceBPs2 != NULL) free(vars->referenceBPs2);
 559   if(vars->ptype != NULL)         free(vars->ptype);
 560   if(vars->S != NULL)             free(vars->S);
 561   if(vars->S1 != NULL)            free(vars->S1);
 562
 563   if(vars->mm1 != NULL)           free(vars->mm1);
 564   if(vars->mm2 != NULL)           free(vars->mm2);
 565   if(vars->bpdist != NULL)        free(vars->bpdist);
 566 #ifdef _OPENMP
 567   }
 568   }
 569 #endif
 570
 571   if(vars->my_iindx != NULL)      free(vars->my_iindx);
 572
 573   free(vars);
 574 }
 575
 576 PRIVATE void initialize_TwoDfold_vars(TwoDfold_vars *vars){
 577   update_TwoDfold_params(vars);
 578   /* this call updates the params in the ViennaRNA fold.o which is a global, so be careful
 579   *  whith calling it parallel... need a workarround or fix of ViennaRNA fold stuff
 580   */
 581   update_fold_params();
 582 }
 583
 584
 585 PUBLIC TwoDfold_solution **TwoDfold(TwoDfold_vars *vars, int distance1, int distance2){
 586   unsigned int i, d1, d2;
 587   unsigned int maxD1;
 588   unsigned int maxD2;
 589   unsigned int length;
 590   TwoDfold_solution **output;
 591
 592   initialize_TwoDfold_vars(vars);
 593   if(fabs(vars->P->temperature - temperature)>1e-6) update_TwoDfold_params(vars);
 594   vars->S   = encode_sequence(vars->sequence, 0);
 595   vars->S1  = encode_sequence(vars->sequence, 1);
 596
 597   make_ptypes(vars);
 598
 599   maxD1 = vars->maxD1;
 600   maxD2 = vars->maxD2;
 601
 602   if(distance1 >= 0){
 603     if((unsigned int)distance1 > maxD1)
 604       fprintf(stderr,
 605               "limiting maximum basepair distance 1 to %u\n",
 606               maxD1);
 607     else
 608       maxD1 = (unsigned int)distance1;
 609   }
 610
 611   if(distance2 >= 0){
 612     if((unsigned int)distance2 > maxD2)
 613       fprintf(stderr,
 614               "limiting maximum basepair distance 2 to %u\n",
 615               maxD2);
 616     else
 617       maxD2 = (unsigned int)distance2;
 618   }
 619
 620   vars->maxD1 = maxD1;
 621   vars->maxD2 = maxD2;
 622   output = (TwoDfold_solution **)space((vars->maxD1+1) * sizeof(TwoDfold_solution *));
 623
 624   mfe_linear(vars);
 625   if(vars->circ) mfe_circ(vars);
 626
 627   length = vars->seq_length;
 628
 629   for(d1=0; d1<=maxD1;d1++){
 630     output[d1] = (TwoDfold_solution *)space((vars->maxD2+1)*sizeof(TwoDfold_solution));
 631 #ifdef _OPENMP
 632   #pragma omp parallel for private(d2)
 633 #endif
 634     for(d2=0; d2<=maxD2;d2++){
 635       output[d1][d2].en = (float)INF/(float)100.;
 636       output[d1][d2].s  = NULL;
 637     }
 638     if(     (d1 >= ((vars->circ) ? vars->k_min_values_fc : vars->k_min_values_f[length]))
 639         &&  (d1 <= ((vars->circ) ? vars->k_max_values_fc : vars->k_max_values_f[length]))){
 640 #ifdef _OPENMP
 641   #pragma omp parallel for private(d2, i)
 642 #endif
 643       for(  d2  = ((vars->circ) ? vars->l_min_values_fc[d1] : vars->l_min_values_f[length][d1]);
 644             d2 <= ((vars->circ) ? vars->l_max_values_fc[d1] : vars->l_max_values_f[length][d1]);
 645             d2 += 2){
 646         output[d1][d2].en = (float)((vars->circ) ? vars->E_Fc[d1][d2/2] : vars->E_F5[length][d1][d2/2])/(float)100.;
 647         if(vars->do_backtrack && (output[d1][d2].en != (float)INF/(float)100.)){
 648           char *mfe_structure = (char *)space(length+1);
 649           for(i=0;i<length;i++) mfe_structure[i] = '.';
 650           mfe_structure[i] = '\0';
 651           (vars->circ) ? backtrack_fc(d1, d2, mfe_structure, vars) : backtrack_f5(length, d1, d2, mfe_structure, vars);
 652           output[d1][d2].s = mfe_structure;
 653         }
 654       }
 655     }
 656
 657   }
 658   return output;
 659 }
 660
 661 PUBLIC TwoDfold_solution *TwoDfoldList(TwoDfold_vars *vars, int distance1, int distance2){
 662   unsigned int  i, d1, d2;
 663   unsigned int  maxD1;
 664   unsigned int  maxD2;
 665   unsigned int  length;
 666   unsigned int  counter = 0;
 667   int           en = 0;
 668   TwoDfold_solution *output;
 669
 670   initialize_TwoDfold_vars(vars);
 671   if(fabs(vars->P->temperature - temperature)>1e-6) update_TwoDfold_params(vars);
 672   vars->S   = encode_sequence(vars->sequence, 0);
 673   vars->S1  = encode_sequence(vars->sequence, 1);
 674
 675   make_ptypes(vars);
 676
 677   maxD1 = vars->maxD1;
 678   maxD2 = vars->maxD2;
 679
 680   if(distance1 >= 0){
 681     if((unsigned int)distance1 > maxD1)
 682       fprintf(stderr,
 683               "TwoDfoldList@2Dfold.c: limiting maximum basepair distance 1 to %u\n",
 684               maxD1);
 685     else
 686       maxD1 = (unsigned int)distance1;
 687   }
 688
 689   if(distance2 >= 0){
 690     if((unsigned int)distance2 > maxD2)
 691       fprintf(stderr,
 692               "TwoDfoldList@2Dfold.c: limiting maximum basepair distance 2 to %u\n",
 693               maxD2);
 694     else
 695       maxD2 = (unsigned int)distance2;
 696   }
 697
 698   vars->maxD1 = maxD1;
 699   vars->maxD2 = maxD2;
 700   output = (TwoDfold_solution *)space((((vars->maxD1+1)*(vars->maxD2+2))/2 + 2) * sizeof(TwoDfold_solution));
 701
 702   mfe_linear(vars);
 703   if(vars->circ) mfe_circ(vars);
 704
 705   length = vars->seq_length;
 706
 707   for(d1=0; d1<=maxD1;d1++){
 708     if((d1 >= ((vars->circ) ? vars->k_min_values_fc : vars->k_min_values_f[length]))
 709         &&  (d1 <= ((vars->circ) ? vars->k_max_values_fc : vars->k_max_values_f[length]))){
 710       for(d2  = ((vars->circ) ? vars->l_min_values_fc[d1] : vars->l_min_values_f[length][d1]);
 711           d2 <= ((vars->circ) ? vars->l_max_values_fc[d1] : vars->l_max_values_f[length][d1]);
 712           d2 += 2){
 713         en = ((vars->circ) ? vars->E_Fc[d1][d2/2] : vars->E_F5[length][d1][d2/2]);
 714         if(en == INF) continue;
 715         output[counter].k   = d1;
 716         output[counter].l   = d2;
 717         output[counter].en  = (float)en/(float)100.;
 718         if(vars->do_backtrack){
 719           char *mfe_structure = (char *)space(length+1);
 720           for(i=0;i<length;i++) mfe_structure[i] = '.';
 721           mfe_structure[i] = '\0';
 722           (vars->circ) ? backtrack_fc((int)d1, (int)d2, mfe_structure, vars) : backtrack_f5(length, (int)d1, (int)d2, mfe_structure, vars);
 723           output[counter].s = mfe_structure;
 724         }
 725         else output[counter].s = NULL;
 726         counter++;
 727       }
 728     }
 729   }
 730
 731   /* store entry for remaining partition if it exists */
 732   en = ((vars->circ) ? vars->E_Fc_rem : vars->E_F5_rem[length]);
 733   if(en != INF){
 734     output[counter].k   = -1;
 735     output[counter].l   = -1;
 736     output[counter].en  =  (float)en/(float)100.;
 737     if(vars->do_backtrack){
 738       char *mfe_structure = (char *)space(length+1);
 739       for(i=0;i<length;i++) mfe_structure[i] = '.';
 740       mfe_structure[i] = '\0';
 741       (vars->circ) ? backtrack_fc(-1, -1, mfe_structure, vars) : backtrack_f5(length, -1, -1, mfe_structure, vars);
 742       output[counter].s = mfe_structure;
 743     }
 744     else output[counter].s = NULL;
 745     counter++;
 746   }
 747
 748   /* insert end-marker entry */
 749   output[counter].k = output[counter].l = INF;
 750   counter++;
 751
 752   /* resize to actual dataset amount */
 753   output = (TwoDfold_solution*)xrealloc(output, sizeof(TwoDfold_solution) * counter);
 754   return output;
 755 }
 756
 757
 758 PUBLIC char *TwoDfold_backtrack_f5(unsigned int j, int k, int l, TwoDfold_vars *vars){
 759   unsigned int i;
 760   char *mfe_structure = (char *)space(j+1);
 761   if(j < TURN + 2) return NULL;
 762
 763   for(i=0; i < j; i++) mfe_structure[i] = '.';
 764   mfe_structure[i] = '\0';
 765
 766   backtrack_f5(j, k, l, mfe_structure, vars);
 767   return mfe_structure;
 768 }
 769
 770 PRIVATE void mfe_linear(TwoDfold_vars *vars){
 771
 772   unsigned int  d, i, j, ij, maxD1, maxD2, seq_length, dia, dib, dja, djb, *referenceBPs1, *referenceBPs2, *mm1, *mm2, *bpdist;
 773   int           cnt1, cnt2, cnt3, cnt4, d1, d2, energy, dangles, temp2, type, additional_en, *my_iindx, circ;
 774   short         *S1, *reference_pt1, *reference_pt2;
 775   char          *sequence, *ptype;
 776   paramT        *P;
 777
 778   /* dereferenciate things we often need */
 779   P               = vars->P;
 780   sequence        = vars->sequence;
 781   seq_length      = vars->seq_length;
 782   maxD1           = vars->maxD1;
 783   maxD2           = vars->maxD2;
 784   S1              = vars->S1;
 785   ptype           = vars->ptype;
 786   reference_pt1   = vars->reference_pt1;
 787   reference_pt2   = vars->reference_pt2;
 788   my_iindx        = vars->my_iindx;
 789   referenceBPs1   = vars->referenceBPs1;
 790   referenceBPs2   = vars->referenceBPs2;
 791   dangles         = vars->dangles;
 792   mm1             = vars->mm1;
 793   mm2             = vars->mm2;
 794   bpdist          = vars->bpdist;
 795   circ            = vars->circ;
 796
 797   for (d = TURN+2; d <= seq_length; d++) { /* i,j in [1..length] */
 798 #ifdef _OPENMP
 799   #pragma omp parallel for private(additional_en, j, energy, temp2, i, ij, dia,dib,dja,djb,cnt1,cnt2,cnt3,cnt4, d1, d2)
 800 #endif
 801     for (j = d; j <= seq_length; j++) {
 802       unsigned int p, q, pq, u, maxp, dij;
 803       int type_2, type, tt, no_close, base_d1, base_d2;
 804
 805       i = j-d+1;
 806       dij = j - i - 1;
 807       ij = my_iindx[i]-j;
 808       type = ptype[ij];
 809
 810       no_close = (((type==3)||(type==4))&&no_closingGU);
 811
 812       if (type) {   /* we have a pair */
 813         /* increase or decrease distance-to-reference value depending whether (i,j) is included in
 814         *  reference or has to be introduced
 815         */
 816         base_d1 = ((unsigned int)reference_pt1[i] != j) ? 1 : -1;
 817         base_d2 = ((unsigned int)reference_pt2[i] != j) ? 1 : -1;
 818
 819         /* HAIRPIN STRUCTURES */
 820
 821         /* get distance to reference if closing the hairpin
 822         *  d = dbp(T_{i,j}, {i,j})
 823         */
 824         d1 = base_d1 + referenceBPs1[ij];
 825         d2 = base_d2 + referenceBPs2[ij];
 826
 827         int min_k, max_k, min_l, max_l;
 828         int real_min_k, real_max_k, *min_l_real, *max_l_real;
 829
 830         min_l = min_k = 0;
 831         max_k = mm1[ij] + referenceBPs1[ij];
 832         max_l = mm2[ij] + referenceBPs2[ij];
 833
 834         prepareBoundaries(min_k,
 835                           max_k,
 836                           min_l,
 837                           max_l,
 838                           bpdist[ij],
 839                           &vars->k_min_values[ij],
 840                           &vars->k_max_values[ij],
 841                           &vars->l_min_values[ij],
 842                           &vars->l_max_values[ij]
 843                           );
 844
 845         preparePosteriorBoundaries( vars->k_max_values[ij] - vars->k_min_values[ij] + 1,
 846                                     vars->k_min_values[ij],
 847                                     &real_min_k,
 848                                     &real_max_k,
 849                                     &min_l_real,
 850                                     &max_l_real
 851                                   );
 852
 853         prepareArray( &vars->E_C[ij],
 854                       vars->k_min_values[ij],
 855                       vars->k_max_values[ij],
 856                       vars->l_min_values[ij],
 857                       vars->l_max_values[ij]
 858                     );
 859
 860 #ifdef COUNT_STATES
 861         prepareArray2( &vars->N_C[ij],
 862                       vars->k_min_values[ij],
 863                       vars->k_max_values[ij],
 864                       vars->l_min_values[ij],
 865                       vars->l_max_values[ij]
 866                     );
 867 #endif
 868
 869         /* d1 and d2 are the distancies to both references introduced by closing a hairpin structure at (i,j) */
 870         if((d1 >= 0) && (d2 >= 0)){
 871           if(((unsigned int)d1<=maxD1) && ((unsigned int)d2 <= maxD2)){
 872             vars->E_C[ij][d1][d2/2] = (no_close) ? FORBIDDEN : E_Hairpin(dij, type, S1[i+1], S1[j-1], sequence+i-1, P);
 873             updatePosteriorBoundaries(d1,
 874                                       d2,
 875                                       &real_min_k,
 876                                       &real_max_k,
 877                                       &min_l_real,
 878                                       &max_l_real
 879                                       );
 880 #ifdef COUNT_STATES
 881             vars->N_C[ij][d1][d2/2] = 1;
 882 #endif
 883           }
 884           else{
 885             vars->E_C_rem[ij] = (no_close) ? FORBIDDEN : E_Hairpin(dij, type, S1[i+1], S1[j-1], sequence+i-1, P);
 886           }
 887         }
 888         /* INTERIOR LOOP STRUCTURES */
 889         maxp = MIN2(j-2-TURN,i+MAXLOOP+1);
 890         for(p = i+1; p <= maxp; p++){
 891           unsigned int minq = p + TURN + 1;
 892           unsigned int ln_pre = dij + p;
 893           if(ln_pre > minq + MAXLOOP) minq = ln_pre - MAXLOOP - 1;
 894           for(q = minq; q < j; q++){
 895             pq = my_iindx[p]-q;
 896             /* set distance to reference structure... */
 897             type_2 = ptype[pq];
 898
 899             if (type_2==0) continue;
 900             type_2 = rtype[type_2];
 901
 902             /* get distance to reference if closing the interior loop
 903             *  d2 = dbp(S_{i,j}, S_{p.q} + {i,j})
 904             */
 905             d1 = base_d1 + referenceBPs1[ij] - referenceBPs1[pq];
 906             d2 = base_d2 + referenceBPs2[ij] - referenceBPs2[pq];
 907
 908             if(no_closingGU)
 909               if(no_close||(type_2==3)||(type_2==4))
 910                 if((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
 911
 912             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);
 913
 914             if(vars->E_C[pq] != NULL){
 915               for(cnt1 = vars->k_min_values[pq]; cnt1 <= vars->k_max_values[pq]; cnt1++){
 916                 for(cnt2 = vars->l_min_values[pq][cnt1]; cnt2 <= vars->l_max_values[pq][cnt1]; cnt2+=2){
 917                   if(vars->E_C[pq][cnt1][cnt2/2] != INF){
 918                     if(((cnt1 + d1) <= maxD1) && ((cnt2+d2) <= maxD2)){
 919                         vars->E_C[ij][cnt1 + d1][(cnt2 + d2)/2] = MIN2( vars->E_C[ij][cnt1 + d1][(cnt2 + d2)/2],
 920                                                                         vars->E_C[pq][cnt1][cnt2/2] + energy
 921                                                                       );
 922                         updatePosteriorBoundaries(cnt1 + d1,
 923                                                   cnt2 + d2,
 924                                                   &real_min_k,
 925                                                   &real_max_k,
 926                                                   &min_l_real,
 927                                                   &max_l_real
 928                                                   );
 929 #ifdef COUNT_STATES
 930                        vars->N_C[ij][cnt1 + d1][(cnt2 + d2)/2] += vars->N_C[pq][cnt1][cnt2/2];
 931 #endif
 932                     }
 933                     /* collect all cases where d1+cnt1 or d2+cnt2 exceeds maxD1, maxD2, respectively */
 934                     else{
 935                       vars->E_C_rem[ij] = MIN2(vars->E_C_rem[ij], vars->E_C[pq][cnt1][cnt2/2] + energy);
 936                     }
 937                   }
 938                 }
 939               }
 940             }
 941             /* collect all contributions where C[pq] already lies outside k_max, l_max boundary */
 942             if(vars->E_C_rem[pq] != INF){
 943               vars->E_C_rem[ij] = MIN2(vars->E_C_rem[ij], vars->E_C_rem[pq] + energy);
 944             }
 945           } /* end q-loop */
 946         } /* end p-loop */
 947
 948
 949
 950         /* MULTI LOOP STRUCTURES */
 951         if(!no_close){
 952
 953           /* dangle energies for multiloop closing stem */
 954           tt = rtype[type];
 955           temp2 = P->MLclosing;
 956           if(dangles == 2)
 957             temp2 += E_MLstem(tt, S1[j-1], S1[i+1], P);
 958           else
 959             temp2 += E_MLstem(tt, -1, -1, P);
 960
 961           for(u=i+TURN+2; u<j-TURN-2;u++){
 962             int i1u   = my_iindx[i+1]-u;
 963             int u1j1  = my_iindx[u+1]-j+1;
 964             /* check all cases where either M or M1 are already out of scope of maxD1 and/or maxD2 */
 965             if(vars->E_M_rem[i1u] != INF){
 966               for(cnt3 = vars->k_min_values_m1[u1j1];
 967                   cnt3 <= vars->k_max_values_m1[u1j1];
 968                   cnt3++)
 969                 for(cnt4 = vars->l_min_values_m1[u1j1][cnt3];
 970                     cnt4 <= vars->l_max_values_m1[u1j1][cnt3];
 971                     cnt4+=2){
 972                   if(vars->E_M1[u1j1][cnt3][cnt4/2]!= INF){
 973                     vars->E_C_rem[ij] = MIN2(vars->E_C_rem[ij],
 974                                               vars->E_M_rem[i1u]
 975                                             + vars->E_M1[u1j1][cnt3][cnt4/2]
 976                                             + temp2
 977                                               );
 978                   }
 979                 }
 980               if(vars->E_M1_rem[u1j1] != INF){
 981                 vars->E_C_rem[ij] = MIN2(vars->E_C_rem[ij],
 982                                           vars->E_M_rem[i1u]
 983                                         + vars->E_M1_rem[u1j1]
 984                                         + temp2
 985                                           );
 986               }
 987             }
 988             if(vars->E_M1_rem[u1j1] != INF){
 989               for(cnt1 = vars->k_min_values_m[i1u];
 990                   cnt1 <= vars->k_max_values_m[i1u];
 991                   cnt1++)
 992                 for(cnt2 = vars->l_min_values_m[i1u][cnt1];
 993                     cnt2 <= vars->l_max_values_m[i1u][cnt1];
 994                     cnt2+=2)
 995                   if(vars->E_M[i1u][cnt1][cnt2/2] != INF){
 996                     vars->E_C_rem[ij] = MIN2(vars->E_C_rem[ij],
 997                                               vars->E_M[i1u][cnt1][cnt2/2]
 998                                             + vars->E_M1_rem[u1j1]
 999                                             + temp2
1000                                               );
1001                   }
1002             }
1003             /* get distance to reference if closing the multiloop
1004             *  d = dbp(S_{i,j}, {i,j} + S_{i+1,u} + S_{u+1,j-1})
1005             */
1006             if(!vars->E_M[i1u]) continue;
1007             if(!vars->E_M1[u1j1]) continue;
1008
1009             d1 = base_d1 + referenceBPs1[ij] - referenceBPs1[i1u] - referenceBPs1[u1j1];
1010             d2 = base_d2 + referenceBPs2[ij] - referenceBPs2[i1u] - referenceBPs2[u1j1];
1011
1012             for(cnt1 = vars->k_min_values_m[i1u];
1013                 cnt1 <= vars->k_max_values_m[i1u];
1014                 cnt1++)
1015               for(cnt2 = vars->l_min_values_m[i1u][cnt1];
1016                   cnt2 <= vars->l_max_values_m[i1u][cnt1];
1017                   cnt2+=2)
1018                 for(cnt3 = vars->k_min_values_m1[u1j1];
1019                     cnt3 <= vars->k_max_values_m1[u1j1];
1020                     cnt3++)
1021                   for(cnt4 = vars->l_min_values_m1[u1j1][cnt3];
1022                       cnt4 <= vars->l_max_values_m1[u1j1][cnt3];
1023                       cnt4+=2){
1024                     if((vars->E_M[i1u][cnt1][cnt2/2] != INF) && (vars->E_M1[u1j1][cnt3][cnt4/2]!= INF)){
1025                       if(((cnt1+cnt3+d1) <= maxD1) && ((cnt2+cnt4+d2) <= maxD2)){
1026                         vars->E_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2] = MIN2( vars->E_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2],
1027                                                                               vars->E_M[i1u][cnt1][cnt2/2]
1028                                                                             + vars->E_M1[u1j1][cnt3][cnt4/2]
1029                                                                             + temp2
1030                                                                             );
1031                         updatePosteriorBoundaries(cnt1 + cnt3 + d1,
1032                                                   cnt2 + cnt4 + d2,
1033                                                   &real_min_k,
1034                                                   &real_max_k,
1035                                                   &min_l_real,
1036                                                   &max_l_real
1037                                                 );
1038 #ifdef COUNT_STATES
1039                         vars->N_C[ij][cnt1+cnt3+d1][(cnt2+cnt4+d2)/2] += vars->N_M[i1u][cnt1][cnt2/2] * vars->N_M1[u1j1][cnt3][cnt4/2];
1040 #endif
1041                       }
1042                       /* collect all cases where d1+cnt1+cnt3 or d2+cnt2+cnt4 exceeds maxD1, maxD2, respectively */
1043                       else{
1044                         vars->E_C_rem[ij] = MIN2(  vars->E_C_rem[ij],
1045                                                     vars->E_M[i1u][cnt1][cnt2/2]
1046                                                   + vars->E_M1[u1j1][cnt3][cnt4/2]
1047                                                   + temp2
1048                                                   );
1049                       }
1050                     }
1051                   }
1052           }
1053         }
1054
1055         /* resize and move memory portions of energy matrix E_C */
1056         adjustArrayBoundaries(&vars->E_C[ij],
1057                               &vars->k_min_values[ij],
1058                               &vars->k_max_values[ij],
1059                               &vars->l_min_values[ij],
1060                               &vars->l_max_values[ij],
1061                               real_min_k,
1062                               real_max_k,
1063                               min_l_real,
1064                               max_l_real
1065                               );
1066 #ifdef COUNT_STATES
1067         /* actually we should adjust the array boundaries here but we might never use the count states option more than once so what....*/
1068 #endif
1069       } /* end >> if (pair) << */
1070
1071
1072
1073       /* done with c[i,j], now compute fML[i,j] */
1074       /* free ends ? -----------------------------------------*/
1075
1076
1077       dia = referenceBPs1[ij] - referenceBPs1[my_iindx[i+1]-j];
1078       dib = referenceBPs2[ij] - referenceBPs2[my_iindx[i+1]-j];
1079       dja = referenceBPs1[ij] - referenceBPs1[ij+1];
1080       djb = referenceBPs2[ij] - referenceBPs2[ij+1];
1081
1082       if(dangles==2)
1083         temp2 = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P);
1084       else
1085         temp2 = E_MLstem(type, -1, -1, P);
1086
1087       int min_k_guess, max_k_guess, min_l_guess, max_l_guess;
1088       int min_k_real_m, max_k_real_m, *min_l_real_m, *max_l_real_m;
1089       int min_k_real_m1, max_k_real_m1, *min_l_real_m1, *max_l_real_m1;
1090
1091       min_k_guess = min_l_guess = 0;
1092       max_k_guess = mm1[ij] + referenceBPs1[ij];
1093       max_l_guess = mm2[ij] + referenceBPs2[ij];
1094
1095       prepareBoundaries(min_k_guess,
1096                         max_k_guess,
1097                         min_l_guess,
1098                         max_l_guess,
1099                         bpdist[ij],
1100                         &vars->k_min_values_m[ij],
1101                         &vars->k_max_values_m[ij],
1102                         &vars->l_min_values_m[ij],
1103                         &vars->l_max_values_m[ij]
1104                         );
1105
1106       prepareBoundaries(min_k_guess,
1107                         max_k_guess,
1108                         min_l_guess,
1109                         max_l_guess,
1110                         bpdist[ij],
1111                         &vars->k_min_values_m1[ij],
1112                         &vars->k_max_values_m1[ij],
1113                         &vars->l_min_values_m1[ij],
1114                         &vars->l_max_values_m1[ij]
1115                         );
1116
1117       preparePosteriorBoundaries( vars->k_max_values_m[ij] - vars->k_min_values_m[ij] + 1,
1118                                   vars->k_min_values_m[ij],
1119                                   &min_k_real_m,
1120                                   &max_k_real_m,
1121                                   &min_l_real_m,
1122                                   &max_l_real_m
1123                                 );
1124       preparePosteriorBoundaries( vars->k_max_values_m1[ij] - vars->k_min_values_m1[ij] + 1,
1125                                   vars->k_min_values_m1[ij],
1126                                   &min_k_real_m1,
1127                                   &max_k_real_m1,
1128                                   &min_l_real_m1,
1129                                   &max_l_real_m1
1130                                 );
1131
1132       prepareArray( &vars->E_M[ij],
1133                     vars->k_min_values_m[ij],
1134                     vars->k_max_values_m[ij],
1135                     vars->l_min_values_m[ij],
1136                     vars->l_max_values_m[ij]
1137                   );
1138
1139       prepareArray( &vars->E_M1[ij],
1140                     vars->k_min_values_m1[ij],
1141                     vars->k_max_values_m1[ij],
1142                     vars->l_min_values_m1[ij],
1143                     vars->l_max_values_m1[ij]
1144                   );
1145 #ifdef COUNT_STATES
1146       prepareArray2( &vars->N_M[ij],
1147                     vars->k_min_values_m[ij],
1148                     vars->k_max_values_m[ij],
1149                     vars->l_min_values_m[ij],
1150                     vars->l_max_values_m[ij]
1151                   );
1152       prepareArray2( &vars->N_M1[ij],
1153                     vars->k_min_values_m1[ij],
1154                     vars->k_max_values_m1[ij],
1155                     vars->l_min_values_m1[ij],
1156                     vars->l_max_values_m1[ij]
1157                   );
1158 #endif
1159
1160       /* now to the actual computations... */
1161       /* 1st E_M[ij] = E_M1[ij] = E_C[ij] + b */
1162       if(vars->E_C_rem[ij] != INF){
1163         vars->E_M_rem[ij] = vars->E_M1_rem[ij] = temp2 + vars->E_C_rem[ij];
1164       }
1165       if(vars->E_C[ij])
1166         for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++){
1167           for(cnt2 = vars->l_min_values[ij][cnt1]; cnt2 <= vars->l_max_values[ij][cnt1]; cnt2+=2){
1168             if(vars->E_C[ij][cnt1][cnt2/2] != INF){
1169               vars->E_M[ij][cnt1][cnt2/2] = vars->E_M1[ij][cnt1][cnt2/2] = temp2 + vars->E_C[ij][cnt1][cnt2/2];
1170               updatePosteriorBoundaries(cnt1,
1171                                         cnt2,
1172                                         &min_k_real_m,
1173                                         &max_k_real_m,
1174                                         &min_l_real_m,
1175                                         &max_l_real_m
1176                                         );
1177               updatePosteriorBoundaries(cnt1,
1178                                         cnt2,
1179                                         &min_k_real_m1,
1180                                         &max_k_real_m1,
1181                                         &min_l_real_m1,
1182                                         &max_l_real_m1
1183                                         );
1184 #ifdef COUNT_STATES
1185              vars->N_M[ij][cnt1][cnt2/2] = vars->N_M1[ij][cnt1][cnt2/2] = vars->N_C[ij][cnt1][cnt2/2];
1186 #endif
1187             }
1188           }
1189         }
1190
1191       /* 2nd E_M[ij] = MIN(E_M[ij], E_M[i+1,j] + c) */
1192       if(vars->E_M_rem[my_iindx[i+1]-j] != INF){
1193         vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1194                                   vars->E_M_rem[my_iindx[i+1]-j] + P->MLbase
1195                                   );
1196       }
1197       if(vars->E_M[my_iindx[i+1]-j])
1198         for(cnt1 = vars->k_min_values_m[my_iindx[i+1]-j];
1199             cnt1 <= vars->k_max_values_m[my_iindx[i+1]-j];
1200             cnt1++){
1201           for(cnt2 = vars->l_min_values_m[my_iindx[i+1]-j][cnt1];
1202               cnt2 <= vars->l_max_values_m[my_iindx[i+1]-j][cnt1];
1203               cnt2+=2){
1204             if(vars->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] != INF){
1205               if(((cnt1 + dia) <= maxD1) && ((cnt2 + dib) <= maxD2)){
1206                 vars->E_M[ij][cnt1+dia][(cnt2+dib)/2] = MIN2( vars->E_M[ij][cnt1+dia][(cnt2+dib)/2],
1207                                                               vars->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase
1208                                                             );
1209                 updatePosteriorBoundaries(cnt1 + dia,
1210                                           cnt2 + dib,
1211                                           &min_k_real_m,
1212                                           &max_k_real_m,
1213                                           &min_l_real_m,
1214                                           &max_l_real_m
1215                                           );
1216 #ifdef COUNT_STATES
1217                 vars->N_M[ij][cnt1+dia][(cnt2+dib)/2] += vars->N_M[my_iindx[i+1]-j][cnt1][cnt2/2];
1218 #endif
1219               }
1220               /* collect all cases where dia+cnt1 or dib+cnt2 exceeds maxD1, maxD2, respectively */
1221               else{
1222                 vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1223                                           vars->E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase
1224                                           );
1225               }
1226             }
1227           }
1228         }
1229
1230       /* 3rd E_M[ij] = MIN(E_M[ij], E_M[i,j-1] + c) */
1231       if(vars->E_M_rem[ij+1] != INF){
1232         vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1233                                   vars->E_M_rem[ij+1] + P->MLbase
1234                                   );
1235       }
1236       if(vars->E_M[ij+1])
1237         for(cnt1 = vars->k_min_values_m[ij+1];
1238             cnt1 <= vars->k_max_values_m[ij+1];
1239             cnt1++){
1240           for(cnt2 = vars->l_min_values_m[ij+1][cnt1];
1241               cnt2 <= vars->l_max_values_m[ij+1][cnt1];
1242               cnt2+=2){
1243             if(vars->E_M[ij+1][cnt1][cnt2/2] != INF){
1244               if(((cnt1 + dja) <= maxD1) && ((cnt2 + djb) <= maxD2)){
1245                 vars->E_M[ij][cnt1+dja][(cnt2+djb)/2] = MIN2( vars->E_M[ij][cnt1+dja][(cnt2+djb)/2],
1246                                                               vars->E_M[ij+1][cnt1][cnt2/2] + P->MLbase
1247                                                             );
1248                 updatePosteriorBoundaries(cnt1 + dja,
1249                                           cnt2 + djb,
1250                                           &min_k_real_m,
1251                                           &max_k_real_m,
1252                                           &min_l_real_m,
1253                                           &max_l_real_m
1254                                           );
1255 #ifdef COUNT_STATES
1256                 vars->N_M[ij][cnt1+dja][(cnt2+djb)/2] += vars->N_M[ij+1][cnt1][cnt2/2];
1257 #endif
1258               }
1259               /* collect all cases where dja+cnt1 or djb+cnt2 exceeds maxD1, maxD2, respectively */
1260               else{
1261                 vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1262                                           vars->E_M[ij+1][cnt1][cnt2/2] + P->MLbase
1263                                           );
1264               }
1265             }
1266           }
1267         }
1268
1269       /* 4th E_M1[ij] = MIN(E_M1[ij], E_M1[i,j-1] + c) */
1270       if(vars->E_M1_rem[ij+1] != INF){
1271         vars->E_M1_rem[ij] = MIN2( vars->E_M1_rem[ij],
1272                                     vars->E_M1_rem[ij+1] + P->MLbase
1273                                   );
1274       }
1275       if(vars->E_M1[ij+1])
1276         for(cnt1 = vars->k_min_values_m1[ij+1];
1277             cnt1 <= vars->k_max_values_m1[ij+1];
1278             cnt1++){
1279           for(cnt2 = vars->l_min_values_m1[ij+1][cnt1];
1280               cnt2 <= vars->l_max_values_m1[ij+1][cnt1];
1281               cnt2+=2){
1282             if(vars->E_M1[ij+1][cnt1][cnt2/2] != INF){
1283               if(((cnt1 + dja) <= maxD1) && ((cnt2 + djb) <= maxD2)){
1284                 vars->E_M1[ij][cnt1+dja][(cnt2+djb)/2]  = MIN2( vars->E_M1[ij][cnt1+dja][(cnt2+djb)/2],
1285                                                                 vars->E_M1[ij+1][cnt1][cnt2/2] + P->MLbase
1286                                                               );
1287                 updatePosteriorBoundaries(cnt1 + dja,
1288                                           cnt2 + djb,
1289                                           &min_k_real_m1,
1290                                           &max_k_real_m1,
1291                                           &min_l_real_m1,
1292                                           &max_l_real_m1
1293                                           );
1294 #ifdef COUNT_STATES
1295                 vars->N_M1[ij][cnt1+dja][(cnt2+djb)/2]  += vars->N_M1[ij+1][cnt1][cnt2/2];
1296 #endif
1297               }
1298               /* collect all cases where dja+cnt1 or djb+cnt2 exceeds maxD1, maxD2, respectively */
1299               else{
1300                 vars->E_M1_rem[ij] = MIN2( vars->E_M1_rem[ij],
1301                                             vars->E_M1[ij+1][cnt1][cnt2/2] + P->MLbase
1302                                           );
1303               }
1304             }
1305           }
1306         }
1307
1308
1309       /* 5th E_M[ij] = MIN(E_M[ij], min(E_M[i,k] + E_M[k+1,j])) */
1310       if(j > TURN + 2)
1311       for (u = i+1+TURN; u <= j-2-TURN; u++){
1312         /* check all cases where M(i,u) and/or M(u+1,j) are already out of scope of maxD1 and/or maxD2 */
1313         if(vars->E_M_rem[my_iindx[i]-u] != INF){
1314           for(cnt3 = vars->k_min_values_m[my_iindx[u+1]-j];
1315               cnt3 <= vars->k_max_values_m[my_iindx[u+1]-j];
1316               cnt3++){
1317             for(cnt4 = vars->l_min_values_m[my_iindx[u+1]-j][cnt3];
1318                 cnt4 <= vars->l_max_values_m[my_iindx[u+1]-j][cnt3];
1319                 cnt4+=2){
1320               if(vars->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] != INF){
1321                   vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1322                                             vars->E_M_rem[my_iindx[i]-u] + vars->E_M[my_iindx[u+1]-j][cnt3][cnt4/2]
1323                                             );
1324               }
1325             }
1326           }
1327           if(vars->E_M_rem[my_iindx[u+1]-j] != INF){
1328             vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1329                                       vars->E_M_rem[my_iindx[i]-u] + vars->E_M_rem[my_iindx[u+1]-j]
1330                                       );
1331           }
1332         }
1333         if(vars->E_M_rem[my_iindx[u+1]-j] != INF){
1334           for(cnt1 = vars->k_min_values_m[my_iindx[i]-u];
1335               cnt1 <= vars->k_max_values_m[my_iindx[i]-u];
1336               cnt1++){
1337             for(cnt2 = vars->l_min_values_m[my_iindx[i]-u][cnt1];
1338                 cnt2 <= vars->l_max_values_m[my_iindx[i]-u][cnt1];
1339                 cnt2+=2){
1340               if(vars->E_M[my_iindx[i]-u][cnt1][cnt2/2] != INF){
1341                 vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1342                                           vars->E_M[my_iindx[i]-u][cnt1][cnt2/2] + vars->E_M_rem[my_iindx[u+1]-j]
1343                                           );
1344               }
1345             }
1346           }
1347         }
1348         if(!vars->E_M[my_iindx[i]-u]) continue;
1349         if(!vars->E_M[my_iindx[u+1]-j]) continue;
1350
1351         dia = referenceBPs1[ij] - referenceBPs1[my_iindx[i]-u] - referenceBPs1[my_iindx[u+1]-j];
1352         dib = referenceBPs2[ij] - referenceBPs2[my_iindx[i]-u] - referenceBPs2[my_iindx[u+1]-j];
1353
1354         for(cnt1 = vars->k_min_values_m[my_iindx[i]-u];
1355             cnt1 <= vars->k_max_values_m[my_iindx[i]-u];
1356             cnt1++){
1357           for(cnt2 = vars->l_min_values_m[my_iindx[i]-u][cnt1];
1358               cnt2 <= vars->l_max_values_m[my_iindx[i]-u][cnt1];
1359               cnt2+=2){
1360             for(cnt3 = vars->k_min_values_m[my_iindx[u+1]-j];
1361                 cnt3 <= vars->k_max_values_m[my_iindx[u+1]-j];
1362                 cnt3++){
1363               for(cnt4 = vars->l_min_values_m[my_iindx[u+1]-j][cnt3];
1364                   cnt4 <= vars->l_max_values_m[my_iindx[u+1]-j][cnt3];
1365                   cnt4+=2){
1366                 if((vars->E_M[my_iindx[i]-u][cnt1][cnt2/2] != INF) && (vars->E_M[my_iindx[u+1]-j][cnt3][cnt4/2] != INF)){
1367                   if(((cnt1 + cnt3 + dia) <= maxD1) && ((cnt2 + cnt4 + dib) <= maxD2)){
1368                     vars->E_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2] = MIN2( vars->E_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2],
1369                                                                             vars->E_M[my_iindx[i]-u][cnt1][cnt2/2]
1370                                                                           + vars->E_M[my_iindx[u+1]-j][cnt3][cnt4/2]
1371                                                                           );
1372                     updatePosteriorBoundaries(cnt1 + cnt3 + dia,
1373                                               cnt2 + cnt4 + dib,
1374                                               &min_k_real_m,
1375                                               &max_k_real_m,
1376                                               &min_l_real_m,
1377                                               &max_l_real_m
1378                                               );
1379 #ifdef COUNT_STATES
1380                     vars->N_M[ij][cnt1+cnt3+dia][(cnt2+cnt4+dib)/2] += vars->N_M[my_iindx[i]-u][cnt1][cnt2/2] * vars->N_M1[my_iindx[u+1]-j][cnt3][cnt4/2];
1381 #endif
1382                   }
1383                   /* collect all cases where dia+cnt1+cnt3 or dib+cnt2+cnt4 exceeds maxD1, maxD2, respectively */
1384                   else{
1385                     vars->E_M_rem[ij] = MIN2(vars->E_M_rem[ij],
1386                                               vars->E_M[my_iindx[i]-u][cnt1][cnt2/2] + vars->E_M[my_iindx[u+1]-j][cnt3][cnt4/2]
1387                                               );
1388                   }
1389                 }
1390               }
1391             }
1392           }
1393         }
1394       }
1395
1396       /* thats all folks for the multiloop decomposition... */
1397
1398       adjustArrayBoundaries(&vars->E_M[ij],
1399                             &vars->k_min_values_m[ij],
1400                             &vars->k_max_values_m[ij],
1401                             &vars->l_min_values_m[ij],
1402                             &vars->l_max_values_m[ij],
1403                             min_k_real_m,
1404                             max_k_real_m,
1405                             min_l_real_m,
1406                             max_l_real_m
1407                             );
1408
1409       adjustArrayBoundaries(&vars->E_M1[ij],
1410                             &vars->k_min_values_m1[ij],
1411                             &vars->k_max_values_m1[ij],
1412                             &vars->l_min_values_m1[ij],
1413                             &vars->l_max_values_m1[ij],
1414                             min_k_real_m1,
1415                             max_k_real_m1,
1416                             min_l_real_m1,
1417                             max_l_real_m1
1418                             );
1419
1420 #ifdef COUNT_STATES
1421         /* actually we should adjust the array boundaries here but we might never use the count states option more than once so what....*/
1422 #endif
1423     } /* end of j-loop */
1424   }
1425
1426   /* calculate energies of 5' and 3' fragments */
1427
1428   /* prepare first entries in E_F5 */
1429   for(cnt1 = 1; cnt1 <= TURN+1; cnt1++){
1430     vars->E_F5[cnt1] = (int **)space(sizeof(int *));
1431     vars->E_F5[cnt1][0] = (int *)space(sizeof(int));
1432     vars->E_F5[cnt1][0][0] = 0;
1433     vars->E_F5_rem[cnt1] = INF;
1434     vars->k_min_values_f[cnt1] = vars->k_max_values_f[cnt1] = 0;
1435     vars->l_min_values_f[cnt1] = (int *)space(sizeof(int));
1436     vars->l_max_values_f[cnt1] = (int *)space(sizeof(int));
1437     vars->l_min_values_f[cnt1][0] = vars->l_max_values_f[cnt1][0] = 0;
1438 #ifdef COUNT_STATES
1439     vars->N_F5[cnt1] = (unsigned long **)space(sizeof(unsigned long *));
1440     vars->N_F5[cnt1][0] = (unsigned long *)space(sizeof(unsigned long));
1441     vars->N_F5[cnt1][0][0] = 1;
1442 #endif
1443
1444   }
1445
1446
1447
1448   for (j=TURN+2; j <= seq_length; j++) {
1449
1450     unsigned int da = referenceBPs1[my_iindx[1]-j] - referenceBPs1[my_iindx[1]-j+1];
1451     unsigned int db = referenceBPs2[my_iindx[1]-j] - referenceBPs2[my_iindx[1]-j+1];
1452
1453     type=ptype[my_iindx[1]-j];
1454     additional_en = 0;
1455     if(type){
1456       if(dangles == 2)
1457         additional_en += E_ExtLoop(type, -1, j < seq_length ? S1[j+1] : -1, P);
1458       else
1459         additional_en += E_ExtLoop(type, -1, -1, P);
1460     }
1461
1462     /* make min and max k guess for memory allocation */
1463     int min_k_guess, max_k_guess, min_l_guess, max_l_guess;
1464     int *min_l_real, *max_l_real, min_k_real, max_k_real;
1465
1466     min_k_guess = min_l_guess = 0;
1467     max_k_guess = referenceBPs1[my_iindx[1]-j] + mm1[my_iindx[1]-j];
1468     max_l_guess = referenceBPs2[my_iindx[1]-j] + mm2[my_iindx[1]-j];
1469
1470     prepareBoundaries(min_k_guess,
1471                       max_k_guess,
1472                       min_l_guess,
1473                       max_l_guess,
1474                       bpdist[my_iindx[1]-j],
1475                       &vars->k_min_values_f[j],
1476                       &vars->k_max_values_f[j],
1477                       &vars->l_min_values_f[j],
1478                       &vars->l_max_values_f[j]
1479                       );
1480
1481     preparePosteriorBoundaries( vars->k_max_values_f[j] - vars->k_min_values_f[j] + 1,
1482                                 vars->k_min_values_f[j],
1483                                 &min_k_real,
1484                                 &max_k_real,
1485                                 &min_l_real,
1486                                 &max_l_real
1487                               );
1488
1489     prepareArray( &vars->E_F5[j],
1490                   vars->k_min_values_f[j],
1491                   vars->k_max_values_f[j],
1492                   vars->l_min_values_f[j],
1493                   vars->l_max_values_f[j]
1494                 );
1495 #ifdef COUNT_STATES
1496     prepareArray2( &vars->N_F5[j],
1497                   vars->k_min_values_f[j],
1498                   vars->k_max_values_f[j],
1499                   vars->l_min_values_f[j],
1500                   vars->l_max_values_f[j]
1501                 );
1502 #endif
1503
1504     /* begin the actual computation of 5' end energies */
1505
1506     /* j-1 is unpaired ... */
1507     vars->E_F5_rem[j] = vars->E_F5_rem[j-1];
1508     for(cnt1 = vars->k_min_values_f[j-1]; cnt1 <= vars->k_max_values_f[j-1]; cnt1++){
1509       for(cnt2 = vars->l_min_values_f[j-1][cnt1]; cnt2 <= vars->l_max_values_f[j-1][cnt1]; cnt2+=2){
1510         if(((cnt1 + da) <= maxD1) && ((cnt2 + db) <= maxD2)){
1511           vars->E_F5[j][cnt1+da][(cnt2+db)/2] = MIN2( vars->E_F5[j][cnt1+da][(cnt2+db)/2],
1512                                                       vars->E_F5[j-1][cnt1][cnt2/2]
1513                                                     );
1514           updatePosteriorBoundaries(cnt1 + da,
1515                                     cnt2 + db,
1516                                     &min_k_real,
1517                                     &max_k_real,
1518                                     &min_l_real,
1519                                     &max_l_real
1520                                     );
1521 #ifdef COUNT_STATES
1522           vars->N_F5[j][cnt1+da][(cnt2+db)/2] += vars->N_F5[j-1][cnt1][cnt2/2];
1523 #endif
1524         }
1525         /* collect all cases where da+cnt1 or db+cnt2 exceeds maxD1, maxD2, respectively */
1526         else{
1527           vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j], vars->E_F5[j-1][cnt1][cnt2/2]);
1528         }
1529       }
1530     }
1531     /* j pairs with 1 */
1532     if(vars->E_C_rem[my_iindx[1]-j] != INF){
1533       vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j], vars->E_C_rem[my_iindx[1]-j] + additional_en);
1534     }
1535     if(vars->E_C[my_iindx[1]-j])
1536       for(cnt1 = vars->k_min_values[my_iindx[1]-j]; cnt1 <= vars->k_max_values[my_iindx[1]-j]; cnt1++)
1537         for(cnt2 = vars->l_min_values[my_iindx[1]-j][cnt1]; cnt2 <= vars->l_max_values[my_iindx[1]-j][cnt1]; cnt2+=2){
1538           if(vars->E_C[my_iindx[1]-j][cnt1][cnt2/2] != INF){
1539             vars->E_F5[j][cnt1][cnt2/2] = MIN2( vars->E_F5[j][cnt1][cnt2/2],
1540                                                 vars->E_C[my_iindx[1]-j][cnt1][cnt2/2]+ additional_en
1541                                               );
1542             updatePosteriorBoundaries(cnt1,
1543                                       cnt2,
1544                                       &min_k_real,
1545                                       &max_k_real,
1546                                       &min_l_real,
1547                                       &max_l_real
1548                                       );
1549 #ifdef COUNT_STATES
1550             vars->N_F5[j][cnt1][cnt2/2] += vars->N_C[my_iindx[1]-j][cnt1][cnt2/2];
1551 #endif
1552           }
1553         }
1554
1555     /* j pairs with some other nucleotide -> see below */
1556     for (i=j-TURN-1; i>1; i--) {
1557       ij = my_iindx[i]-j;
1558       type = ptype[ij];
1559       if (type) {
1560         if(dangles == 2)
1561           additional_en = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P);
1562         else
1563           additional_en = E_ExtLoop(type, -1, -1, P);
1564
1565         if(vars->E_C_rem[ij] != INF){
1566           for(cnt3 = vars->k_min_values_f[i-1]; cnt3 <= vars->k_max_values_f[i-1]; cnt3++)
1567             for(cnt4 = vars->l_min_values_f[i-1][cnt3]; cnt4 <= vars->l_max_values_f[i-1][cnt3]; cnt4+=2){
1568               if(vars->E_F5[i-1][cnt3][cnt4/2] != INF){
1569                 vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j],
1570                                           vars->E_F5[i-1][cnt3][cnt4/2] + vars->E_C_rem[ij] + additional_en
1571                                           );
1572               }
1573             }
1574           if(vars->E_F5_rem[i-1] != INF){
1575             vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j],
1576                                       vars->E_F5_rem[i-1] + vars->E_C_rem[ij] + additional_en
1577                                       );
1578           }
1579         }
1580         if((vars->E_F5_rem[i-1] != INF) && (vars->E_C[ij])){
1581           for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++)
1582             for(cnt2 = vars->l_min_values[ij][cnt1]; cnt2 <= vars->l_max_values[ij][cnt1]; cnt2+=2)
1583               if(vars->E_C[ij][cnt1][cnt2/2]!= INF){
1584                 vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j],
1585                                           vars->E_F5_rem[i-1] + vars->E_C[ij][cnt1][cnt2/2] + additional_en
1586                                           );
1587               }
1588         }
1589         if(!vars->E_C[ij]) continue;
1590
1591         unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1];
1592         unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1];
1593
1594         for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++)
1595           for(cnt2 = vars->l_min_values[ij][cnt1]; cnt2 <= vars->l_max_values[ij][cnt1]; cnt2+=2)
1596             for(cnt3 = vars->k_min_values_f[i-1]; cnt3 <= vars->k_max_values_f[i-1]; cnt3++)
1597               for(cnt4 = vars->l_min_values_f[i-1][cnt3]; cnt4 <= vars->l_max_values_f[i-1][cnt3]; cnt4+=2){
1598                 if(vars->E_F5[i-1][cnt3][cnt4/2] != INF && vars->E_C[ij][cnt1][cnt2/2]!= INF){
1599                   if(((cnt1 + cnt3 + d1a) <= maxD1) && ((cnt2 + cnt4 + d1b) <= maxD2)){
1600                     vars->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] = MIN2( vars->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2],
1601                                                                             vars->E_F5[i-1][cnt3][cnt4/2] + vars->E_C[ij][cnt1][cnt2/2] + additional_en
1602                                                                           );
1603                     updatePosteriorBoundaries(cnt1 + cnt3 + d1a,
1604                                               cnt2 + cnt4 + d1b,
1605                                               &min_k_real,
1606                                               &max_k_real,
1607                                               &min_l_real,
1608                                               &max_l_real
1609                                               );
1610 #ifdef COUNT_STATES
1611                     vars->N_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] += vars->N_F5[i-1][cnt3][cnt4/2] * vars->N_C[ij][cnt1][cnt2/2];
1612 #endif
1613                   }
1614                   /* collect all cases where d1a+cnt1+cnt3 or d1b+cnt2+cnt4 exceeds maxD1, maxD2, respectively */
1615                   else{
1616                     vars->E_F5_rem[j] = MIN2(vars->E_F5_rem[j],
1617                                               vars->E_F5[i-1][cnt3][cnt4/2] + vars->E_C[ij][cnt1][cnt2/2] + additional_en
1618                                               );
1619                   }
1620                 }
1621               }
1622       }
1623     }
1624
1625     /* resize and move memory portions of energy matrix E_F5 */
1626     adjustArrayBoundaries(&vars->E_F5[j],
1627                           &vars->k_min_values_f[j],
1628                           &vars->k_max_values_f[j],
1629                           &vars->l_min_values_f[j],
1630                           &vars->l_max_values_f[j],
1631                           min_k_real,
1632                           max_k_real,
1633                           min_l_real,
1634                           max_l_real
1635                           );
1636
1637   } /* end of j-loop */
1638
1639
1640   if(compute_2Dfold_F3){
1641
1642     /* prepare first entries in E_F3 */
1643     for(cnt1 = seq_length; cnt1 >= seq_length-TURN-1; cnt1--){
1644       vars->E_F3[cnt1]        = (int **)space(sizeof(int *));
1645       vars->E_F3[cnt1][0]     = (int *) space(sizeof(int));
1646       vars->E_F3[cnt1][0][0]  = 0;
1647       vars->k_min_values_f3[cnt1]     = vars->k_max_values_f3[cnt1] = 0;
1648       vars->l_min_values_f3[cnt1]     = (int *)space(sizeof(int));
1649       vars->l_max_values_f3[cnt1]     = (int *)space(sizeof(int));
1650       vars->l_min_values_f3[cnt1][0]  = vars->l_max_values_f3[cnt1][0] = 0;
1651     }
1652     /* begin calculations */
1653     for (j=seq_length-TURN-2; j >= 1; j--){
1654
1655       unsigned int da = referenceBPs1[my_iindx[j]-seq_length] - referenceBPs1[my_iindx[j+1]-seq_length];
1656       unsigned int db = referenceBPs2[my_iindx[j]-seq_length] - referenceBPs2[my_iindx[j+1]-seq_length];
1657
1658       type=ptype[my_iindx[j]-seq_length];
1659       additional_en = 0;
1660       if(type){
1661         if(dangles == 2)
1662           additional_en += E_ExtLoop(type, j > 1 ? S1[j-1] : -1, -1, P);
1663         else
1664           additional_en += E_ExtLoop(type, -1, -1, P);
1665       }
1666
1667       /* make min and max k guess for memory allocation */
1668       int min_k_guess, max_k_guess, min_l_guess, max_l_guess;
1669       int *min_l_real, *max_l_real, min_k_real, max_k_real;
1670
1671       min_k_guess = min_l_guess = 0;
1672       max_k_guess = referenceBPs1[my_iindx[j]-seq_length] + mm1[my_iindx[j]-seq_length];
1673       max_l_guess = referenceBPs2[my_iindx[j]-seq_length] + mm2[my_iindx[j]-seq_length];
1674
1675       prepareBoundaries(min_k_guess,
1676                         max_k_guess,
1677                         min_l_guess,
1678                         max_l_guess,
1679                         bpdist[my_iindx[j]-seq_length],
1680                         &vars->k_min_values_f3[j],
1681                         &vars->k_max_values_f3[j],
1682                         &vars->l_min_values_f3[j],
1683                         &vars->l_max_values_f3[j]
1684                         );
1685
1686       preparePosteriorBoundaries( vars->k_max_values_f3[j] - vars->k_min_values_f3[j] + 1,
1687                                   vars->k_min_values_f3[j],
1688                                   &min_k_real,
1689                                   &max_k_real,
1690                                   &min_l_real,
1691                                   &max_l_real
1692                                 );
1693
1694       prepareArray( &vars->E_F3[j],
1695                     vars->k_min_values_f3[j],
1696                     vars->k_max_values_f3[j],
1697                     vars->l_min_values_f3[j],
1698                     vars->l_max_values_f3[j]
1699                   );
1700       /* begin the actual computation of 5' end energies */
1701
1702       /* j is unpaired ... */
1703       for(cnt1 = vars->k_min_values_f3[j+1]; cnt1 <= vars->k_max_values_f3[j+1]; cnt1++){
1704         for(cnt2 = vars->l_min_values_f3[j+1][cnt1]; cnt2 <= vars->l_max_values_f3[j+1][cnt1]; cnt2+=2){
1705           vars->E_F3[j][cnt1+da][(cnt2+db)/2] = MIN2( vars->E_F3[j][cnt1+da][(cnt2+db)/2],
1706                                                       vars->E_F3[j+1][cnt1][cnt2/2]
1707                                                     );
1708           updatePosteriorBoundaries(cnt1 + da,
1709                                     cnt2 + db,
1710                                     &min_k_real,
1711                                     &max_k_real,
1712                                     &min_l_real,
1713                                     &max_l_real
1714                                     );
1715         }
1716       }
1717       /* j pairs with n */
1718       if(vars->E_C[my_iindx[j]-seq_length])
1719         for(cnt1 = vars->k_min_values[my_iindx[j]-seq_length]; cnt1 <= vars->k_max_values[my_iindx[j]-seq_length]; cnt1++)
1720           for(cnt2 = vars->l_min_values[my_iindx[j]-seq_length][cnt1]; cnt2 <= vars->l_max_values[my_iindx[j]-seq_length][cnt1]; cnt2+=2){
1721             if(vars->E_C[my_iindx[j]-seq_length][cnt1][cnt2/2] != INF){
1722               vars->E_F3[j][cnt1][cnt2/2] = MIN2( vars->E_F3[j][cnt1][cnt2/2],
1723                                                   vars->E_C[my_iindx[j]-seq_length][cnt1][cnt2/2]+ additional_en
1724                                                 );
1725               updatePosteriorBoundaries(cnt1,
1726                                         cnt2,
1727                                         &min_k_real,
1728                                         &max_k_real,
1729                                         &min_l_real,
1730                                         &max_l_real
1731                                         );
1732             }
1733           }
1734
1735       /* j pairs with some other nucleotide -> see below */
1736       for (i=j-TURN-1; i>1; i--) {
1737         ij = my_iindx[i]-j;
1738         if(!vars->E_C[ij]) continue;
1739         type = ptype[ij];
1740         if (type) {
1741           unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1];
1742           unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1];
1743
1744           if(dangles == 2)
1745             additional_en = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P);
1746           else
1747             additional_en = E_ExtLoop(type, -1, -1, P);
1748
1749           for(cnt1 = vars->k_min_values[ij]; cnt1 <= vars->k_max_values[ij]; cnt1++)
1750             for(cnt2 = vars->l_min_values[ij][cnt1]; cnt2 <= vars->l_max_values[ij][cnt1]; cnt2+=2)
1751               for(cnt3 = vars->k_min_values_f[i-1]; cnt3 <= vars->k_max_values_f[i-1]; cnt3++)
1752                 for(cnt4 = vars->l_min_values_f[i-1][cnt3]; cnt4 <= vars->l_max_values_f[i-1][cnt3]; cnt4+=2){
1753                   if(vars->E_F5[i-1][cnt3][cnt4/2] != INF && vars->E_C[ij][cnt1][cnt2/2]!= INF){
1754                     vars->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] = MIN2( vars->E_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2],
1755                                                                             vars->E_F5[i-1][cnt3][cnt4/2] + vars->E_C[ij][cnt1][cnt2/2] + additional_en
1756                                                                           );
1757                     updatePosteriorBoundaries(cnt1 + cnt3 + d1a,
1758                                               cnt2 + cnt4 + d1b,
1759                                               &min_k_real,
1760                                               &max_k_real,
1761                                               &min_l_real,
1762                                               &max_l_real
1763                                               );
1764 #ifdef COUNT_STATES
1765                     vars->N_F5[j][cnt1+cnt3+d1a][(cnt2+cnt4+d1b)/2] += vars->N_F5[i-1][cnt3][cnt4/2] * vars->N_C[ij][cnt1][cnt2/2];
1766 #endif
1767                   }
1768                 }
1769         }
1770       }
1771
1772       /* resize and move memory portions of energy matrix E_F5 */
1773       adjustArrayBoundaries(&vars->E_F5[j],
1774                             &vars->k_min_values_f[j],
1775                             &vars->k_max_values_f[j],
1776                             &vars->l_min_values_f[j],
1777                             &vars->l_max_values_f[j],
1778                             min_k_real,
1779                             max_k_real,
1780                             min_l_real,
1781                             max_l_real
1782                             );
1783
1784     } /* end of j-loop */
1785
1786
1787
1788
1789   }
1790 }
1791
1792
1793 /*---------------------------------------------------------------------------*/
1794
1795 PUBLIC void update_TwoDfold_params(TwoDfold_vars *vars){
1796   if(vars->P) free(vars->P);
1797   vars->P = scale_parameters();
1798   make_pair_matrix();
1799 }
1800
1801 /*---------------------------------------------------------------------------*/
1802 PRIVATE void make_ptypes(TwoDfold_vars *vars) {
1803   int n,i,j,k,l;
1804
1805   n=vars->S[0];
1806   for (k=1; k<n-TURN; k++)
1807     for (l=1; l<=2; l++) {
1808       int type,ntype=0,otype=0;
1809       i=k; j = i+TURN+l; if (j>n) continue;
1810       type = pair[vars->S[i]][vars->S[j]];
1811       while ((i>=1)&&(j<=n)) {
1812         if ((i>1)&&(j<n)) ntype = pair[vars->S[i-1]][vars->S[j+1]];
1813         if (noLonelyPairs && (!otype) && (!ntype))
1814           type = 0; /* i.j can only form isolated pairs */
1815         vars->ptype[vars->my_iindx[i]-j] = (char) type;
1816         otype =  type;
1817         type  = ntype;
1818         i--; j++;
1819       }
1820     }
1821 }
1822
1823 PRIVATE void backtrack_f5(unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars){
1824   int           *my_iindx, energy, type, dangles, cnt1, cnt2, cnt3, cnt4;
1825   int           **l_min_values, **l_max_values,**l_min_values_f, **l_max_values_f;
1826   int           *k_min_values, *k_max_values,*k_min_values_f, *k_max_values_f;
1827   int           ***E_C, ***E_F5;
1828   int           *E_C_rem, *E_F5_rem;
1829   unsigned int   i, ij, seq_length, maxD1, maxD2;
1830   short *S1;
1831   unsigned int   *referenceBPs1, *referenceBPs2;
1832   char  *ptype;
1833   paramT   *P;
1834   unsigned int da, db;
1835
1836   P               = vars->P;
1837   seq_length      = vars->seq_length;
1838   S1              = vars->S1;
1839   ptype           = vars->ptype;
1840   my_iindx        = vars->my_iindx;
1841   referenceBPs1   = vars->referenceBPs1;
1842   referenceBPs2   = vars->referenceBPs2;
1843   dangles         = vars->dangles;
1844   E_F5            = vars->E_F5;
1845   l_min_values_f  = vars->l_min_values_f;
1846   l_max_values_f  = vars->l_max_values_f;
1847   k_min_values_f  = vars->k_min_values_f;
1848   k_max_values_f  = vars->k_max_values_f;
1849
1850   E_C             = vars->E_C;
1851   l_min_values    = vars->l_min_values;
1852   l_max_values    = vars->l_max_values;
1853   k_min_values    = vars->k_min_values;
1854   k_max_values    = vars->k_max_values;
1855
1856   E_F5_rem       = vars->E_F5_rem;
1857   E_C_rem        = vars->E_C_rem;
1858   maxD1           = vars->maxD1;
1859   maxD2           = vars->maxD2;
1860
1861   da = referenceBPs1[my_iindx[1]-j] - referenceBPs1[my_iindx[1]-j+1];
1862   db = referenceBPs2[my_iindx[1]-j] - referenceBPs2[my_iindx[1]-j+1];
1863
1864   if(j<TURN+2) return;
1865
1866   /* F5[j] == F5[j-1] ? */
1867   if(k == -1){
1868     if(E_F5_rem[j]==INF)
1869       return;
1870     else if(E_F5_rem[j] == E_F5_rem[j-1]){
1871       backtrack_f5(j-1,k,l,structure, vars);
1872       return;
1873     }
1874     else if(E_F5[j-1]){
1875       for(cnt1 =  k_min_values_f[j-1];
1876         cnt1 <= k_max_values_f[j-1];
1877         cnt1++){
1878         for(cnt2 = l_min_values_f[j-1][cnt1];
1879             cnt2 <= l_max_values_f[j-1][cnt1];
1880             cnt2+=2){
1881           if(((cnt1 + da) > maxD1) || ((cnt2 + db) > maxD2)){
1882             if(E_F5_rem[j] == E_F5[j-1][cnt1][cnt2/2]){
1883               backtrack_f5(j-1, cnt1, cnt2, structure, vars);
1884               return;
1885             }
1886           }
1887         }
1888       }
1889     }
1890   }
1891   else if((k >= da) && (l >= db)){
1892     if(E_F5[j-1]){
1893       if((k - da >= k_min_values_f[j-1]) && (k - da <= k_max_values_f[j-1])){
1894         if((l - db >= l_min_values_f[j-1][k-da]) && (l - db <= l_max_values_f[j-1][k-da]))
1895           if(E_F5[j-1][k-da][(l-db)/2] == E_F5[j][k][l/2]){
1896             backtrack_f5(j-1, k-da, l-db, structure, vars);
1897             return;
1898           }
1899       }
1900     }
1901   }
1902
1903   type = ptype[my_iindx[1]-j];
1904   if(type){
1905     if(dangles == 2)
1906       energy = E_ExtLoop(type, -1, j < seq_length ? S1[j+1] : -1, P);
1907     else
1908       energy = E_ExtLoop(type, -1, -1, P);
1909
1910     if(k == -1){
1911       if(E_C_rem[my_iindx[1]-j] + energy == E_F5_rem[j]){
1912           backtrack_c(1, j, -1, -1, structure, vars);
1913           return;
1914       }
1915     }
1916     else if(k >= k_min_values[my_iindx[1]-j] && (k <= k_max_values[my_iindx[1]-j])){
1917
1918       if((l >= l_min_values[my_iindx[1]-j][k]) && (l <= l_max_values[my_iindx[1]-j][k]))
1919         if(E_C[my_iindx[1]-j][k][l/2] + energy == E_F5[j][k][l/2]){
1920           backtrack_c(1, j, k, l, structure, vars);
1921           return;
1922         }
1923     }
1924   }
1925
1926   for (i=j-TURN-1; i>1; i--) {
1927     ij = my_iindx[i]-j;
1928     type = ptype[ij];
1929     if (type) {
1930       unsigned int d1a = referenceBPs1[my_iindx[1]-j] - referenceBPs1[ij] - referenceBPs1[my_iindx[1]-i+1];
1931       unsigned int d1b = referenceBPs2[my_iindx[1]-j] - referenceBPs2[ij] - referenceBPs2[my_iindx[1]-i+1];
1932
1933       if(dangles == 2)
1934         energy = E_ExtLoop(type, S1[i-1], j < seq_length ? S1[j+1] : -1, P);
1935       else
1936         energy = E_ExtLoop(type, -1, -1, P);
1937
1938       if(k == -1){
1939         if(E_C_rem[ij] != INF){
1940           for(cnt1 = k_min_values_f[i-1];
1941               cnt1 <= k_max_values_f[i-1];
1942               cnt1++){
1943             for(cnt2 = l_min_values_f[i-1][cnt1];
1944                 cnt2 <= l_max_values_f[i-1][cnt1];
1945                 cnt2+=2){
1946               if(E_F5_rem[j] == (E_F5[i-1][cnt1][cnt2/2] + E_C_rem[ij] + energy)){
1947                 backtrack_f5(i-1, cnt1, cnt2, structure, vars);
1948                 backtrack_c(i,j,-1,-1,structure, vars);
1949                 return;
1950               }
1951             }
1952           }
1953           if(E_F5_rem[j] == (E_F5_rem[i-1] + E_C_rem[ij] + energy)){
1954             backtrack_f5(i-1, -1, -1, structure, vars);
1955             backtrack_c(i,j,-1,-1,structure,vars);
1956             return;
1957           }
1958         }
1959         if(E_F5_rem[i-1] != INF){
1960           for(cnt1 = k_min_values[ij];
1961               cnt1 <= k_max_values[ij];
1962               cnt1++){
1963             for(cnt2 = l_min_values[ij][cnt1];
1964                 cnt2 <= l_max_values[ij][cnt1];
1965                 cnt2 += 2){
1966                if(E_F5_rem[j] == (E_F5_rem[i-1] + E_C[ij][cnt1][cnt2/2] + energy)){
1967                 backtrack_f5(i-1,-1,-1,structure,vars);
1968                 backtrack_c(i,j,cnt1,cnt2,structure,vars);
1969                 return;
1970               }
1971             }
1972           }
1973         }
1974         for(cnt1 = k_min_values_f[i-1];
1975             cnt1 <= k_max_values_f[i-1];
1976             cnt1++)
1977           for(cnt2 = l_min_values_f[i-1][cnt1];
1978               cnt2 <= l_max_values_f[i-1][cnt1];
1979               cnt2 += 2)
1980             for(cnt3 = k_min_values[ij];
1981                 cnt3 <= k_max_values[ij];
1982                 cnt3++)
1983               for(cnt4 = l_min_values[ij][cnt3];
1984                   cnt4 <= l_max_values[ij][cnt3];
1985                   cnt4 += 2){
1986                 if(((cnt1 + cnt3 + d1a)>maxD1) || ((cnt2+cnt4+d1b)>maxD2)){
1987                   if(E_F5_rem[j] == (E_F5[i-1][cnt1][cnt2/2] + E_C[ij][cnt3][cnt4/2] + energy)){
1988                     backtrack_f5(i-1,cnt1,cnt2,structure,vars);
1989                     backtrack_c(i,j,cnt3,cnt4,structure,vars);
1990                     return;
1991                   }
1992                 }
1993               }
1994       }
1995       else if((k >= d1a) && (l >= d1b)){
1996         int k_f_max = MIN2(k-d1a, k_max_values_f[i-1]);
1997
1998         for(cnt1 = k_min_values_f[i-1]; cnt1 <= k_f_max; cnt1++){
1999           int l_f_max = MIN2(l - d1b, l_max_values_f[i-1][cnt1]);
2000           for(cnt2 = l_min_values_f[i-1][cnt1]; cnt2 <= l_f_max; cnt2+=2){
2001             int k_c = k - d1a - cnt1;
2002             if((k_c >= k_min_values[ij]) && (k_c <= k_max_values[ij])){
2003               int l_c = l - d1b - cnt2;
2004               if((l_c >= l_min_values[ij][k_c]) && (l_c <= l_max_values[ij][k_c])){
2005                 if(E_F5[j][k][l/2] == (E_F5[i-1][cnt1][cnt2/2] + E_C[ij][k_c][l_c/2] + energy)){
2006                   backtrack_f5(i-1, cnt1, cnt2, structure, vars);
2007                   backtrack_c(i, j, k_c, l_c, structure, vars);
2008                   return;
2009                 }
2010               }
2011             }
2012           }
2013
2014         }
2015       }
2016     }
2017   }
2018   nrerror("backtracking failed in f5");
2019 }
2020
2021 PRIVATE void backtrack_c(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars){
2022   unsigned int p, q, pq, ij, maxp, maxD1, maxD2;
2023   int *my_iindx, type, type_2, energy, no_close, dangles, base_d1, base_d2, d1, d2, cnt1, cnt2, cnt3, cnt4;
2024   int           **l_min_values, **l_max_values,**l_min_values_m, **l_max_values_m,**l_min_values_m1, **l_max_values_m1;
2025   int           *k_min_values, *k_max_values,*k_min_values_m, *k_max_values_m,*k_min_values_m1, *k_max_values_m1;
2026   int           ***E_C, ***E_M, ***E_M1, *E_C_rem, *E_M_rem, *E_M1_rem;
2027   short *S1;
2028   unsigned int   *referenceBPs1, *referenceBPs2;
2029   char  *ptype, *sequence;
2030   paramT   *P;
2031
2032   P               = vars->P;
2033   sequence        = vars->sequence;
2034   S1              = vars->S1;
2035   ptype           = vars->ptype;
2036   my_iindx        = vars->my_iindx;
2037   referenceBPs1   = vars->referenceBPs1;
2038   referenceBPs2   = vars->referenceBPs2;
2039   dangles         = vars->dangles;
2040
2041   E_C             = vars->E_C;
2042   l_min_values    = vars->l_min_values;
2043   l_max_values    = vars->l_max_values;
2044   k_min_values    = vars->k_min_values;
2045   k_max_values    = vars->k_max_values;
2046
2047   E_M             = vars->E_M;
2048   l_min_values_m  = vars->l_min_values_m;
2049   l_max_values_m  = vars->l_max_values_m;
2050   k_min_values_m  = vars->k_min_values_m;
2051   k_max_values_m  = vars->k_max_values_m;
2052
2053   E_M1            = vars->E_M1;
2054   l_min_values_m1 = vars->l_min_values_m1;
2055   l_max_values_m1 = vars->l_max_values_m1;
2056   k_min_values_m1 = vars->k_min_values_m1;
2057   k_max_values_m1 = vars->k_max_values_m1;
2058
2059   E_C_rem        = vars->E_C_rem;
2060   E_M_rem        = vars->E_M_rem;
2061   E_M1_rem       = vars->E_M1_rem;
2062   maxD1           = vars->maxD1;
2063   maxD2           = vars->maxD2;
2064
2065
2066   ij = my_iindx[i]-j;
2067
2068   int e = (k==-1) ? E_C_rem[ij] : E_C[ij][k][l/2];
2069
2070   type = ptype[ij];
2071
2072   no_close = (((type==3)||(type==4))&&no_closingGU);
2073   structure[i-1] = '(';
2074   structure[j-1] = ')';
2075
2076   base_d1 = ((unsigned int)vars->reference_pt1[i] != j) ? 1 : -1;
2077   base_d2 = ((unsigned int)vars->reference_pt2[i] != j) ? 1 : -1;
2078
2079   base_d1 += referenceBPs1[ij];
2080   base_d2 += referenceBPs2[ij];
2081
2082   if(k == -1){
2083     if(((unsigned int)base_d1 > maxD1) || ((unsigned int)base_d2 > maxD2)){
2084       if(e == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], sequence+i-1, P)) return;
2085     }
2086   }
2087   else{
2088     if((unsigned int)base_d1 == k)
2089       if((unsigned int)base_d2 == l)
2090         if(E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], sequence+i-1, P) == e) return;
2091   }
2092   maxp = MIN2(j-2-TURN,i+MAXLOOP+1);
2093   for(p = i+1; p <= maxp; p++){
2094     unsigned int minq, ln_pre;
2095     minq = p + TURN + 1;
2096     ln_pre = j - i - 1;
2097     if(ln_pre > minq + MAXLOOP) minq = ln_pre - MAXLOOP - 1;
2098     for (q = minq; q < j; q++) {
2099       pq = my_iindx[p]-q;
2100       type_2 = ptype[pq];
2101       if (type_2==0) continue;
2102       type_2 = rtype[type_2];
2103
2104       /* d2 = dbp(S_{i,j}, S_{p.q} + {i,j}) */
2105       d1 = base_d1 - referenceBPs1[pq];
2106       d2 = base_d2 - referenceBPs2[pq];
2107
2108       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);
2109
2110
2111       if(k == -1){
2112         if(E_C_rem[pq] != INF)
2113           if(e == (E_C_rem[pq] + energy)){
2114             backtrack_c(p,q,-1,-1,structure,vars);
2115             return;
2116           }
2117         if(E_C[pq])
2118         for(cnt1 = k_min_values[pq];
2119             cnt1 <= k_max_values[pq];
2120             cnt1++)
2121           for(cnt2 = l_min_values[pq][cnt1];
2122               cnt2 <= l_max_values[pq][cnt1];
2123               cnt2 += 2){
2124             if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){
2125               if(e == (E_C[pq][cnt1][cnt2/2] + energy)){
2126                 backtrack_c(p,q,cnt1,cnt2,structure,vars);
2127                 return;
2128               }
2129             }
2130           }
2131       }
2132       else{
2133         if(!E_C[pq]) continue;
2134         if(d1 <= k && d2 <= l){
2135           if((k-d1 >= k_min_values[pq]) && (k-d1) <= k_max_values[pq])
2136             if((l - d2 >= l_min_values[pq][k-d1]) && (l-d2 <= l_max_values[pq][k-d1]))
2137               if(E_C[pq][k-d1][(l-d2)/2] + energy == e){
2138                 backtrack_c(p, q, k-d1, l-d2, structure, vars);
2139                 return;
2140               }
2141         }
2142       }
2143     } /* end q-loop */
2144   } /* end p-loop */
2145
2146   /* multi-loop decomposition ------------------------*/
2147   if(!no_close){
2148     unsigned int u;
2149     int tt;
2150     if(k==-1){
2151       for(u=i+TURN+2; u<j-TURN-2;u++){
2152         int i1u, u1j1;
2153         i1u   = my_iindx[i+1]-u;
2154         u1j1  = my_iindx[u+1]-j+1;
2155         tt = rtype[type];
2156         energy = P->MLclosing;
2157         if(dangles == 2)
2158           energy += E_MLstem(tt, S1[j-1], S1[i+1], P);
2159         else
2160           energy += E_MLstem(tt, -1, -1, P);
2161
2162
2163         if(E_M_rem[i1u] != INF){
2164           if(E_M1[u1j1])
2165           for(cnt1 = k_min_values_m1[u1j1];
2166               cnt1 <= k_max_values_m1[u1j1];
2167               cnt1++)
2168             for(cnt2 = l_min_values_m1[u1j1][cnt1];
2169                 cnt2 <= l_max_values_m1[u1j1][cnt1];
2170                 cnt2 += 2){
2171               if(e == (E_M_rem[i1u] + E_M1[u1j1][cnt1][cnt2/2] + energy)){
2172                 backtrack_m(i+1,u,-1,-1,structure,vars);
2173                 backtrack_m1(u+1,j-1,cnt1,cnt2,structure,vars);
2174                 return;
2175               }
2176             }
2177           if(E_M1_rem[u1j1] != INF){
2178             if(e == (E_M_rem[i1u] + E_M1_rem[u1j1] + energy)){
2179               backtrack_m(i+1, u, -1, -1, structure, vars);
2180               backtrack_m1(u+1, j-1, -1, -1, structure, vars);
2181               return;
2182             }
2183           }
2184         }
2185         if(E_M1_rem[u1j1] != INF){
2186           if(E_M[i1u])
2187           for(cnt1 = k_min_values_m[i1u];
2188               cnt1 <= k_max_values_m[i1u];
2189               cnt1++)
2190             for(cnt2 = l_min_values_m[i1u][cnt1];
2191                 cnt2 <= l_max_values_m[i1u][cnt1];
2192                 cnt2 += 2)
2193               if(e == (E_M[i1u][cnt1][cnt2/2] + E_M1_rem[u1j1] + energy)){
2194                 backtrack_m(i+1,u,cnt1,cnt2,structure,vars);
2195                 backtrack_m1(u+1,j-1,-1,-1,structure,vars);
2196                 return;
2197               }
2198         }
2199
2200         /* now all cases where we exceed the maxD1/D2 scope by combination of E_M and E_M1 */
2201         if(!E_M[i1u]) continue;
2202         if(!E_M1[u1j1]) continue;
2203         /* get distance to reference if closing this multiloop
2204         *  dist3 = dbp(S_{i,j}, {i,j} + S_{i+1.u} + S_{u+1,j-1})
2205         */
2206         d1 = base_d1 - referenceBPs1[i1u] - referenceBPs1[u1j1];
2207         d2 = base_d2 - referenceBPs2[i1u] - referenceBPs2[u1j1];
2208
2209         for(cnt1 = vars->k_min_values_m[i1u];
2210             cnt1 <= vars->k_max_values_m[i1u];
2211             cnt1++)
2212           for(cnt2 = vars->l_min_values_m[i1u][cnt1];
2213               cnt2 <= vars->l_max_values_m[i1u][cnt1];
2214               cnt2+=2)
2215             for(cnt3 = vars->k_min_values_m1[u1j1];
2216                 cnt3 <= vars->k_max_values_m1[u1j1];
2217                 cnt3++)
2218               for(cnt4 = vars->l_min_values_m1[u1j1][cnt3];
2219                   cnt4 <= vars->l_max_values_m1[u1j1][cnt3];
2220                   cnt4+=2){
2221                 if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){
2222                   if(e == (E_M[i1u][cnt1][cnt2/2] + E_M1[u1j1][cnt3][cnt4/2] + energy)){
2223                     backtrack_m(i+1,u,cnt1,cnt2,structure,vars);
2224                     backtrack_m1(u+1,j-1,cnt3,cnt4,structure,vars);
2225                     return;
2226                   }
2227                 }
2228               }
2229       }
2230     }
2231     else{
2232       for(u=i+TURN+2; u<j-TURN-2;u++){
2233         int i1u, u1j1;
2234         i1u   = my_iindx[i+1]-u;
2235         u1j1  = my_iindx[u+1]-j+1;
2236         if(!E_M[i1u]) continue;
2237         if(!E_M1[u1j1]) continue;
2238
2239         /* get distance to reference if closing this multiloop
2240         *  dist3 = dbp(S_{i,j}, {i,j} + S_{i+1.u} + S_{u+1,j-1})
2241         */
2242         d1 = base_d1 - referenceBPs1[i1u] - referenceBPs1[u1j1];
2243         d2 = base_d2 - referenceBPs2[i1u] - referenceBPs2[u1j1];
2244
2245         tt = rtype[type];
2246         energy = P->MLclosing;
2247         if(dangles == 2)
2248           energy += E_MLstem(tt, S1[j-1], S1[i+1], P);
2249         else
2250           energy += E_MLstem(tt, -1, -1, P);
2251
2252         if((d1 <= k) && (d2 <= l))
2253           for(cnt1 = k_min_values_m[i1u];
2254               cnt1 <= MIN2(k-d1, k_max_values_m[i1u]);
2255               cnt1++)
2256             for(cnt2 = l_min_values_m[i1u][cnt1];
2257                 cnt2 <= MIN2(l-d2, l_max_values_m[i1u][cnt1]);
2258                 cnt2+=2)
2259               if(     ((k-d1-cnt1) >= k_min_values_m1[u1j1])
2260                   &&  ((k-d1-cnt1) <= k_max_values_m1[u1j1]))
2261                 if(     ((l-d2-cnt2) >= l_min_values_m1[u1j1][k-d1-cnt1])
2262                     &&  ((l-d2-cnt2) <= l_max_values_m1[u1j1][k-d1-cnt1]))
2263                   if(e == (energy + E_M[i1u][cnt1][cnt2/2] + E_M1[u1j1][k-d1-cnt1][(l-d2-cnt2)/2])){
2264                     backtrack_m(i+1, u, cnt1, cnt2, structure, vars);
2265                     backtrack_m1(u+1, j-1, k-d1-cnt1, l-d2-cnt2, structure, vars);
2266                     return;
2267                   }
2268       }
2269     }
2270   }
2271   nrerror("backtracking failed in c");
2272 }
2273
2274 PRIVATE void backtrack_m(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars){
2275   unsigned int u, ij, seq_length, base_d1, base_d2, d1, d2, maxD1, maxD2;
2276   int *my_iindx, type, energy, dangles,circ, cnt1, cnt2, cnt3, cnt4;
2277   int           **l_min_values, **l_max_values,**l_min_values_m, **l_max_values_m;
2278   int           *k_min_values, *k_max_values,*k_min_values_m, *k_max_values_m;
2279   int           ***E_C, ***E_M, *E_C_rem, *E_M_rem;
2280   short *S1;
2281   unsigned int   *referenceBPs1, *referenceBPs2;
2282   char  *ptype;
2283   paramT   *P;
2284
2285   P           = vars->P;
2286   seq_length  = vars->seq_length;
2287   S1          = vars->S1;
2288   circ        = vars->circ;
2289   ptype       = vars->ptype;
2290   my_iindx    = vars->my_iindx;
2291   referenceBPs1  = vars->referenceBPs1;
2292   referenceBPs2  = vars->referenceBPs2;
2293   dangles     = vars->dangles;
2294
2295   E_C             = vars->E_C;
2296   l_min_values    = vars->l_min_values;
2297   l_max_values    = vars->l_max_values;
2298   k_min_values    = vars->k_min_values;
2299   k_max_values    = vars->k_max_values;
2300
2301   E_M             = vars->E_M;
2302   l_min_values_m  = vars->l_min_values_m;
2303   l_max_values_m  = vars->l_max_values_m;
2304   k_min_values_m  = vars->k_min_values_m;
2305   k_max_values_m  = vars->k_max_values_m;
2306
2307   E_C_rem        = vars->E_C_rem;
2308   E_M_rem        = vars->E_M_rem;
2309   maxD1           = vars->maxD1;
2310   maxD2           = vars->maxD2;
2311
2312   ij = my_iindx[i]-j;
2313   int e = (k == -1) ? E_M_rem[ij] : E_M[ij][k][l/2];
2314
2315   base_d1 = referenceBPs1[ij];
2316   base_d2 = referenceBPs2[ij];
2317
2318   if(k == -1){
2319     /* new_fML = ML(i+1,j)+c */
2320     d1 = base_d1 - referenceBPs1[my_iindx[i+1]-j];
2321     d2 = base_d2 - referenceBPs2[my_iindx[i+1]-j];
2322     if(E_M_rem[my_iindx[i+1]-j] != INF){
2323       if(e == (E_M_rem[my_iindx[i+1]-j] + P->MLbase)){
2324         backtrack_m(i+1,j,-1,-1,structure,vars);
2325         return;
2326       }
2327     }
2328     if(E_M[my_iindx[i+1]-j])
2329     for(cnt1 = k_min_values_m[my_iindx[i+1]-j];
2330         cnt1 <= k_max_values_m[my_iindx[i+1]-j];
2331         cnt1++)
2332       for(cnt2 = l_min_values_m[my_iindx[i+1]-j][cnt1];
2333           cnt2 <= l_max_values_m[my_iindx[i+1]-j][cnt1];
2334           cnt2 += 2)
2335         if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){
2336           if(e == (E_M[my_iindx[i+1]-j][cnt1][cnt2/2] + P->MLbase)){
2337             backtrack_m(i+1,j,cnt1,cnt2,structure,vars);
2338             return;
2339           }
2340         }
2341
2342     /* new_fML = min(ML(i,j-1) + c, new_fML) */
2343     d1 = base_d1 - referenceBPs1[ij+1];
2344     d2 = base_d2 - referenceBPs2[ij+1];
2345     if(E_M_rem[ij+1] != INF){
2346       if(e == (E_M_rem[ij+1] + P->MLbase)){
2347         backtrack_m(i,j-1,-1,-1,structure,vars);
2348         return;
2349       }
2350     }
2351     if(E_M[ij+1])
2352     for(cnt1 = k_min_values_m[ij+1];
2353         cnt1 <= k_max_values_m[ij+1];
2354         cnt1++)
2355       for(cnt2 = l_min_values_m[ij+1][cnt1];
2356           cnt2 <= l_max_values_m[ij+1][cnt1];
2357           cnt2 += 2)
2358         if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){
2359           if(e == (E_M[ij+1][cnt1][cnt2/2] + P->MLbase)){
2360             backtrack_m(i,j-1,cnt1,cnt2,structure,vars);
2361             return;
2362           }
2363         }
2364
2365     /* new_fML = min(new_fML, C(i,j)+b) */
2366     if(E_C_rem[ij] != INF){
2367       type = ptype[ij];
2368       if(dangles == 2)
2369         energy = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P);
2370       else
2371         energy = E_MLstem(type, -1, -1, P);
2372       if(e == (E_C_rem[ij] + energy)){
2373         backtrack_c(i,j,-1,-1,structure,vars);
2374         return;
2375       }
2376     }
2377
2378     /* modular decomposition -------------------------------*/
2379     for(u = i+1+TURN; u <= j-2-TURN; u++){
2380       int iu, uj;
2381       iu = my_iindx[i]-u;
2382       uj = my_iindx[u+1]-j;
2383       type = ptype[uj];
2384
2385       d1 = base_d1 - referenceBPs1[iu] - referenceBPs1[uj];
2386       d2 = base_d2 - referenceBPs2[iu] - referenceBPs2[uj];
2387
2388       if(dangles == 2)
2389         energy = E_MLstem(type, S1[u], (j < seq_length) || circ ? S1[j+1] : -1, P);
2390       else
2391         energy = E_MLstem(type, -1, -1, P);
2392
2393       if(E_M_rem[iu] != INF){
2394         if(E_C[uj])
2395         for(cnt1 = k_min_values[uj];
2396             cnt1 <= k_max_values[uj];
2397             cnt1++)
2398           for(cnt2 = l_min_values[uj][cnt1];
2399               cnt2 <= l_max_values[uj][cnt1];
2400               cnt2 += 2)
2401             if(e == (E_M_rem[iu] + E_C[uj][cnt1][cnt2/2] + energy)){
2402               backtrack_m(i,u,-1,-1,structure,vars);
2403               backtrack_c(u+1,j,cnt1,cnt2,structure, vars);
2404               return;
2405             }
2406         if(E_C_rem[uj] != INF){
2407           if(e == (E_M_rem[iu] + E_C_rem[uj] + energy)){
2408             backtrack_m(i,u,-1,-1,structure,vars);
2409             backtrack_c(u+1,j,-1,-1,structure,vars);
2410             return;
2411           }
2412         }
2413       }
2414       if(E_C_rem[uj] != INF){
2415         if(E_M[iu])
2416         for(cnt1 = k_min_values_m[iu];
2417             cnt1 <= k_max_values_m[iu];
2418             cnt1++)
2419           for(cnt2 = l_min_values_m[iu][cnt1];
2420               cnt2 <= l_max_values_m[iu][cnt1];
2421               cnt2 += 2)
2422             if(e == (E_M[iu][cnt1][cnt2/2] + E_C_rem[uj] + energy)){
2423               backtrack_m(i,u,cnt1,cnt2,structure,vars);
2424               backtrack_c(u+1,j,-1,-1,structure,vars);
2425               return;
2426             }
2427       }
2428
2429       if(!E_M[iu]) continue;
2430       if(!E_C[uj]) continue;
2431
2432       for(cnt1 = k_min_values_m[iu];
2433           cnt1 <= k_max_values_m[iu];
2434           cnt1++)
2435         for(cnt2 = l_min_values_m[iu][cnt1];
2436             cnt2 <= l_max_values_m[iu][cnt1];
2437             cnt2 += 2)
2438           for(cnt3 = k_min_values[uj];
2439               cnt3 <= k_max_values[uj];
2440               cnt3++){
2441             for(cnt4 = l_min_values[uj][cnt3];
2442                 cnt4 <= l_max_values[uj][cnt3];
2443                 cnt4 += 2)
2444               if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2))
2445                 if(e == (E_M[iu][cnt1][cnt2/2] + E_C[uj][cnt3][cnt4/2] + energy)){
2446                   backtrack_m(i, u, cnt1, cnt2, structure, vars);
2447                   backtrack_c(u+1, j, cnt3, cnt4, structure, vars);
2448                   return;
2449                 }
2450           }
2451     }
2452
2453   } /* end if (k == -1) */
2454   else{
2455     d1 = base_d1 - referenceBPs1[my_iindx[i+1]-j];
2456     d2 = base_d2 - referenceBPs2[my_iindx[i+1]-j];
2457     /* new_fML = ML(i+1,j)+c */
2458     if(d1 <= k && d2 <= l)
2459       if((k-d1 >= k_min_values_m[my_iindx[i+1]-j]) && (k-d1 <= k_max_values_m[my_iindx[i+1]-j]))
2460         if((l-d2 >= l_min_values_m[my_iindx[i+1]-j][k-d1]) && (l-d2 <= l_max_values_m[my_iindx[i+1]-j][k-d1])){
2461           if(E_M[my_iindx[i+1]-j][k-d1][(l-d2)/2] + P->MLbase == e){
2462             backtrack_m(i+1, j, k-d1, l-d2, structure, vars);
2463             return;
2464           }
2465         }
2466
2467     d1 = base_d1 - referenceBPs1[ij+1];
2468     d2 = base_d2 - referenceBPs2[ij+1];
2469
2470     /* new_fML = min(ML(i,j-1) + c, new_fML) */
2471     if(E_M[ij+1])
2472       if(d1 <= k && d2 <= l)
2473         if((k-d1 >= k_min_values_m[ij+1]) && (k-d1 <= k_max_values_m[ij+1]))
2474           if((l-d2 >= l_min_values_m[ij+1][k-d1]) && (l-d2 <= l_max_values_m[ij+1][k-d1]))
2475             if(E_M[ij+1][k-d1][(l-d2)/2] + P->MLbase == e){
2476               backtrack_m(i, j-1, k-d1, l-d2, structure, vars);
2477               return;
2478             }
2479
2480     /* new_fML = min(new_fML, C(i,j)+b) */
2481     if(E_C[ij]){
2482       type = ptype[ij];
2483
2484       if(dangles == 2)
2485         energy = E_MLstem(type, ((i > 1) || circ) ? S1[i-1] : -1, ((j < seq_length) || circ) ? S1[j+1] : -1, P);
2486       else
2487         energy = E_MLstem(type, -1, -1, P);
2488
2489       if((k >= k_min_values[ij]) && (k <= k_max_values[ij]))
2490         if((l >= l_min_values[ij][k]) && (l <= l_max_values[ij][k])){
2491           if(E_C[ij][k][l/2] + energy == e){
2492             backtrack_c(i, j, k, l, structure, vars);
2493             return;
2494           }
2495         }
2496     }
2497
2498       /* modular decomposition -------------------------------*/
2499
2500     for(u = i+1+TURN; u <= j-2-TURN; u++){
2501       if(!E_M[my_iindx[i]-u]) continue;
2502       if(!E_C[my_iindx[u+1]-j]) continue;
2503       type = ptype[my_iindx[u+1]-j];
2504
2505       d1 = base_d1 - referenceBPs1[my_iindx[i]-u] - referenceBPs1[my_iindx[u+1]-j];
2506       d2 = base_d2 - referenceBPs2[my_iindx[i]-u] - referenceBPs2[my_iindx[u+1]-j];
2507
2508       if(dangles == 2)
2509         energy = E_MLstem(type, S1[u], ((j < seq_length) || circ) ? S1[j+1] : -1, P);
2510       else
2511         energy = E_MLstem(type, -1, -1, P);
2512
2513       if(d1 <= k && d2 <= l)
2514         for(cnt1 = k_min_values_m[my_iindx[i]-u]; cnt1 <= MIN2(k-d1, k_max_values_m[my_iindx[i]-u]); cnt1++)
2515           for(cnt2 = l_min_values_m[my_iindx[i]-u][cnt1]; cnt2 <= MIN2(l-d2, l_max_values_m[my_iindx[i]-u][cnt1]); cnt2+=2)
2516             if((k-d1-cnt1 >= k_min_values[my_iindx[u+1]-j]) && (k-d1-cnt1 <= k_max_values[my_iindx[u+1]-j]))
2517               if((l-d2-cnt2 >= l_min_values[my_iindx[u+1]-j][k-d1-cnt1]) && (l-d2-cnt2 <= l_max_values[my_iindx[u+1]-j][k-d1-cnt1]))
2518                 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){
2519                   backtrack_m(i, u, cnt1, cnt2, structure, vars);
2520                   backtrack_c(u+1, j, k-d1-cnt1, l-d2-cnt2, structure, vars);
2521                   return;
2522                 }
2523     }
2524   }
2525   nrerror("backtracking failed in fML\n");
2526 }
2527
2528 PRIVATE void backtrack_m1(unsigned int i, unsigned int j, int k, int l, char *structure, TwoDfold_vars *vars){
2529   unsigned int  ij, seq_length, d1, d2, *referenceBPs1, *referenceBPs2, maxD1, maxD2;
2530   int           *my_iindx, **l_min_values, **l_max_values,**l_min_values_m1, **l_max_values_m1;
2531   int           *k_min_values, *k_max_values,*k_min_values_m1, *k_max_values_m1, cnt1, cnt2;
2532   int           ***E_C, ***E_M1, *E_C_rem, *E_M1_rem, type, dangles, circ, energy, e_m1;
2533
2534   short         *S1;
2535   char          *ptype;
2536   paramT        *P;
2537
2538   P               = vars->P;
2539   seq_length      = vars->seq_length;
2540   S1              = vars->S1;
2541   ptype           = vars->ptype;
2542   circ            = vars->circ;
2543   my_iindx        = vars->my_iindx;
2544   referenceBPs1   = vars->referenceBPs1;
2545   referenceBPs2   = vars->referenceBPs2;
2546   dangles         = vars->dangles;
2547
2548   E_C             = vars->E_C;
2549   l_min_values    = vars->l_min_values;
2550   l_max_values    = vars->l_max_values;
2551   k_min_values    = vars->k_min_values;
2552   k_max_values    = vars->k_max_values;
2553
2554   E_M1            = vars->E_M1;
2555   l_min_values_m1 = vars->l_min_values_m1;
2556   l_max_values_m1 = vars->l_max_values_m1;
2557   k_min_values_m1 = vars->k_min_values_m1;
2558   k_max_values_m1 = vars->k_max_values_m1;
2559
2560   E_C_rem        = vars->E_C_rem;
2561   E_M1_rem       = vars->E_M1_rem;
2562   maxD1           = vars->maxD1;
2563   maxD2           = vars->maxD2;
2564
2565   ij    = my_iindx[i]-j;
2566   e_m1  = (k == -1) ? E_M1_rem[ij] : E_M1[ij][k][l/2];
2567
2568   type = ptype[ij];
2569   d1 = referenceBPs1[ij] - referenceBPs1[ij+1];
2570   d2 = referenceBPs2[ij] - referenceBPs2[ij+1];
2571
2572   if(dangles == 2)
2573     energy = E_MLstem(type, (i > 1) || circ ? S1[i-1] : -1, (j < seq_length) || circ ? S1[j+1] : -1, P);
2574   else
2575     energy = E_MLstem(type, -1, -1, P);
2576
2577   if(k == -1){
2578     if(E_C_rem[ij] != INF){
2579       if(e_m1 == (E_C_rem[ij] + energy)){
2580         backtrack_c(i,j,-1,-1,structure,vars);
2581         return;
2582       }
2583     }
2584     if(E_M1_rem[ij+1] != INF){
2585       if(e_m1 == (E_M1_rem[ij+1] + P->MLbase)){
2586         backtrack_m1(i,j-1,-1,-1,structure,vars);
2587         return;
2588       }
2589     }
2590     for(cnt1 = k_min_values_m1[ij+1];
2591         cnt1 <= k_max_values_m1[ij+1];
2592         cnt1++)
2593       for(cnt2 = l_min_values_m1[ij+1][cnt1];
2594           cnt2 <= l_max_values_m1[ij+1][cnt1];
2595           cnt2 += 2)
2596         if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2)){
2597           if(e_m1  == (E_M1[ij+1][cnt1][cnt2/2] + P->MLbase)){
2598             backtrack_m1(i,j-1,cnt1,cnt2,structure,vars);
2599             return;
2600           }
2601         }
2602   }
2603   else{
2604     if(E_C[ij])
2605       if((k >= k_min_values[ij]) && (k <= k_max_values[ij]))
2606         if((l >= l_min_values[ij][k]) && (l <= l_max_values[ij][k]))
2607           if(E_C[ij][k][l/2] + energy == e_m1){
2608             backtrack_c(i, j, k, l, structure, vars);
2609             return;
2610           }
2611
2612     if(d1 <= k && d2 <= l)
2613       if((k-d1 >= k_min_values_m1[ij+1]) && (k-d1 <= k_max_values_m1[ij+1]))
2614         if((l-d2 >= l_min_values_m1[ij+1][k-d1]) && (l-d2 <= l_max_values_m1[ij+1][k-d1]))
2615           if(E_M1[ij+1][k-d1][(l-d2)/2] + P->MLbase == e_m1){
2616             backtrack_m1(i, j-1, k-d1, l-d2, structure, vars);
2617             return;
2618           }
2619   }
2620   nrerror("backtack failed in m1\n");
2621 }
2622
2623 PRIVATE void backtrack_fc(int k, int l, char *structure, TwoDfold_vars *vars){
2624   unsigned int   d, i, j, seq_length, base_d1, base_d2, d1, d2, maxD1, maxD2;
2625   int   *my_iindx, energy, cnt1, cnt2, cnt3, cnt4;
2626   short *S1;
2627   unsigned int   *referenceBPs1, *referenceBPs2;
2628   char  *sequence, *ptype;
2629   int **E_Fc, **E_FcH, **E_FcI, **E_FcM, ***E_C, ***E_M, ***E_M2;
2630   int *E_C_rem, *E_M_rem, *E_M2_rem, E_Fc_rem, E_FcH_rem, E_FcI_rem, E_FcM_rem;
2631   int **l_min_values, **l_max_values, *k_min_values, *k_max_values;
2632   int **l_min_values_m, **l_max_values_m, *k_min_values_m, *k_max_values_m;
2633   int **l_min_values_m2, **l_max_values_m2, *k_min_values_m2, *k_max_values_m2;
2634   int *l_min_values_fcH, *l_max_values_fcH, k_min_values_fcH, k_max_values_fcH;
2635   int *l_min_values_fcI, *l_max_values_fcI, k_min_values_fcI, k_max_values_fcI;
2636   int *l_min_values_fcM, *l_max_values_fcM, k_min_values_fcM, k_max_values_fcM;
2637   paramT   *P;
2638   P                 = vars->P;
2639   sequence          = vars->sequence;
2640   seq_length        = vars->seq_length;
2641   S1                = vars->S1;
2642   ptype             = vars->ptype;
2643   my_iindx          = vars->my_iindx;
2644   referenceBPs1     = vars->referenceBPs1;
2645   referenceBPs2     = vars->referenceBPs2;
2646
2647   base_d1           = referenceBPs1[my_iindx[1]-seq_length];
2648   base_d2           = referenceBPs2[my_iindx[1]-seq_length];
2649
2650   E_C               = vars->E_C;
2651   l_min_values      = vars->l_min_values;
2652   l_max_values      = vars->l_max_values;
2653   k_min_values      = vars->k_min_values;
2654   k_max_values      = vars->k_max_values;
2655
2656   E_M               = vars->E_M;
2657   l_min_values_m    = vars->l_min_values_m;
2658   l_max_values_m    = vars->l_max_values_m;
2659   k_min_values_m    = vars->k_min_values_m;
2660   k_max_values_m    = vars->k_max_values_m;
2661
2662   E_M2              = vars->E_M2;
2663   l_min_values_m2   = vars->l_min_values_m2;
2664   l_max_values_m2   = vars->l_max_values_m2;
2665   k_min_values_m2   = vars->k_min_values_m2;
2666   k_max_values_m2   = vars->k_max_values_m2;
2667
2668   E_Fc              = vars->E_Fc;
2669
2670   E_FcI             = vars->E_FcI;
2671   l_min_values_fcI  = vars->l_min_values_fcI;
2672   l_max_values_fcI  = vars->l_max_values_fcI;
2673   k_min_values_fcI  = vars->k_min_values_fcI;
2674   k_max_values_fcI  = vars->k_max_values_fcI;
2675
2676   E_FcH             = vars->E_FcH;
2677   l_min_values_fcH  = vars->l_min_values_fcH;
2678   l_max_values_fcH  = vars->l_max_values_fcH;
2679   k_min_values_fcH  = vars->k_min_values_fcH;
2680   k_max_values_fcH  = vars->k_max_values_fcH;
2681
2682   E_FcM             = vars->E_FcM;
2683   l_min_values_fcM  = vars->l_min_values_fcM;
2684   l_max_values_fcM  = vars->l_max_values_fcM;
2685   k_min_values_fcM  = vars->k_min_values_fcM;
2686   k_max_values_fcM  = vars->k_max_values_fcM;
2687
2688   E_C_rem          = vars->E_C_rem;
2689   E_M_rem          = vars->E_M_rem;
2690   E_M2_rem         = vars->E_M2_rem;
2691   E_Fc_rem         = vars->E_Fc_rem;
2692   E_FcH_rem        = vars->E_FcH_rem;
2693   E_FcI_rem        = vars->E_FcI_rem;
2694   E_FcM_rem        = vars->E_FcM_rem;
2695   maxD1             = vars->maxD1;
2696   maxD2             = vars->maxD2;
2697
2698
2699   if(k==-1){
2700     /* check if mfe might be open chain */
2701     if(E_Fc_rem == 0)
2702       if((referenceBPs1[my_iindx[1]-seq_length] > maxD1) || (referenceBPs2[my_iindx[1]-seq_length] > maxD2))
2703         return;
2704
2705     /* check for hairpin configurations */
2706     if(E_Fc_rem == E_FcH_rem){
2707       for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
2708         for (j = d; j <= seq_length; j++) {
2709           unsigned int u, ij;
2710           int type, no_close;
2711           char loopseq[10];
2712           i = j-d+1;
2713           ij = my_iindx[i]-j;
2714           u = seq_length-j + i-1;
2715           if (u<TURN) continue;
2716           type = ptype[ij];
2717           no_close = (((type==3)||(type==4))&&no_closingGU);
2718           type=rtype[type];
2719           if (!type) continue;
2720           if(no_close) continue;
2721
2722           d1 = base_d1 - referenceBPs1[ij];
2723           d2 = base_d2 - referenceBPs2[ij];
2724           if (u<7) {
2725             strcpy(loopseq , sequence+j-1);
2726             strncat(loopseq, sequence, i);
2727           }
2728           energy = E_Hairpin(u, type, S1[j+1], S1[i-1],  loopseq, P);
2729
2730           if(E_C_rem[ij] != INF){
2731             if(E_Fc_rem == (E_C_rem[ij] + energy)){
2732               backtrack_c(i,j,-1,-1,structure,vars);
2733               return;
2734             }
2735           }
2736           if(E_C[ij])
2737             for(cnt1 = k_min_values[ij];
2738                 cnt1 <= k_max_values[ij];
2739                 cnt1++)
2740               for(cnt2 = l_min_values[ij][cnt1];
2741                   cnt2 <= l_max_values[ij][cnt1];
2742                   cnt2 += 2)
2743                 if(((cnt1 + d1) > maxD1) || ((cnt2 + d2) > maxD2))
2744                   if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + energy)){
2745                     backtrack_c(i,j,cnt1,cnt2,structure,vars);
2746                     return;
2747                   }
2748         }
2749     }
2750
2751     /* check for interior loop configurations */
2752     if(E_Fc_rem == E_FcI_rem){
2753       for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
2754         for (j = d; j <= seq_length; j++) {
2755           unsigned int u, ij, p, q, pq;
2756           int type, type_2;
2757           i = j-d+1;
2758           ij = my_iindx[i]-j;
2759           u = seq_length-j + i-1;
2760           if (u<TURN) continue;
2761           type = rtype[(unsigned int)ptype[ij]];
2762           if (!type) continue;
2763
2764           for(p = j+1; p < seq_length ; p++){
2765             unsigned int u1, qmin, ln_pre;
2766             u1 = p-j-1;
2767             if (u1+i-1>MAXLOOP) break;
2768             qmin = p + TURN + 1;
2769             ln_pre = u1 + i + seq_length;
2770             if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1;
2771             for(q = qmin; q <= seq_length; q++){
2772               unsigned int u2;
2773               pq = my_iindx[p]-q;
2774               type_2 = rtype[(unsigned int)ptype[pq]];
2775               if (type_2==0) continue;
2776               u2 = i-1 + seq_length-q;
2777               if (u1+u2>MAXLOOP) continue;
2778               energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P);
2779               if(E_C_rem[ij] != INF){
2780                 if(E_C[pq])
2781                   for(cnt1 = k_min_values[pq];
2782                       cnt1 <= k_max_values[pq];
2783                       cnt1++)
2784                     for(cnt2 = l_min_values[pq][cnt1];
2785                         cnt2 <= l_max_values[pq][cnt1];
2786                         cnt2 += 2)
2787                       if(E_Fc_rem == (E_C_rem[ij] + E_C[pq][cnt1][cnt2/2] + energy)){
2788                         backtrack_c(i,j,-1,-1,structure,vars);
2789                         backtrack_c(p,q,cnt1,cnt2,structure,vars);
2790                         return;
2791                       }
2792                 if(E_C_rem[pq] != INF){
2793                   if(E_Fc_rem == (E_C_rem[ij] + E_C_rem[pq] + energy)){
2794                     backtrack_c(i,j,-1,-1,structure,vars);
2795                     backtrack_c(p,q,-1,-1,structure,vars);
2796                     return;
2797                   }
2798                 }
2799               }
2800               if(E_C_rem[pq] != INF){
2801                 if(E_C[ij])
2802                   for(cnt1 = k_min_values[ij];
2803                       cnt1 <= k_max_values[ij];
2804                       cnt1++)
2805                     for(cnt2 = l_min_values[ij][cnt1];
2806                         cnt2 <= l_max_values[ij][cnt1];
2807                         cnt2 += 2)
2808                       if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + E_C_rem[pq] + energy)){
2809                         backtrack_c(i,j,cnt1,cnt2,structure,vars);
2810                         backtrack_c(p,q,-1,-1,structure,vars);
2811                         return;
2812                       }
2813               }
2814
2815               if(!(E_C[ij])) continue;
2816               if(!(E_C[pq])) continue;
2817
2818               /* get distance to reference if closing the interior loop
2819               *  d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j})
2820               *  d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j})
2821               */
2822               d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq];
2823               d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq];
2824               for(cnt1 = k_min_values[ij];
2825                   cnt1 <= k_max_values[ij];
2826                   cnt1++)
2827                 for(cnt2 = l_min_values[ij][cnt1];
2828                     cnt2 <= l_max_values[ij][cnt1];
2829                     cnt2 += 2)
2830                   for(cnt3 = k_min_values[pq];
2831                       cnt3 <= k_max_values[pq];
2832                       cnt3++)
2833                     for(cnt4 = l_min_values[pq][cnt3];
2834                         cnt4 <= l_max_values[pq][cnt3];
2835                         cnt4 += 2)
2836                       if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2))
2837                         if(E_Fc_rem == (E_C[ij][cnt1][cnt2/2] + E_C[pq][cnt3][cnt4/2] + energy)){
2838                           backtrack_c(i, j, cnt1, cnt2, structure, vars);
2839                           backtrack_c(p, q, cnt3, cnt4, structure, vars);
2840                           return;
2841                         }
2842             } /* end for p */
2843           } /* end for q */
2844         }
2845
2846     }
2847
2848     /* check for multi loop configurations */
2849     if(E_Fc_rem == E_FcM_rem){
2850       if(seq_length > 2*TURN)
2851         for (i=TURN+1; i<seq_length-2*TURN; i++) {
2852           /* get distancies to references
2853           * d3a = dbp(T1_[1,n}, T1_{1,k} + T1_{k+1, n})
2854           * d3b = dbp(T2_[1,n}, T2_{1,k} + T2_{k+1, n})
2855           */
2856           if(E_M_rem[my_iindx[1]-i] != INF){
2857             if(E_M2[i+1])
2858               for(cnt1 = k_min_values_m2[i+1];
2859                   cnt1 <= k_max_values_m2[i+1];
2860                   cnt1++)
2861                 for(cnt2 = l_min_values_m2[i+1][cnt1];
2862                     cnt2 <= l_max_values_m2[i+1][cnt1];
2863                     cnt2 += 2)
2864                   if(E_Fc_rem == (E_M_rem[my_iindx[1]-i] + E_M2[i+1][cnt1][cnt2/2] + P->MLclosing)){
2865                     backtrack_m(1,i,-1,-1,structure,vars);
2866                     backtrack_m2(i+1,cnt1,cnt2,structure,vars);
2867                     return;
2868                   }
2869             if(E_M2_rem[i+1] != INF){
2870               if(E_Fc_rem == (E_M_rem[my_iindx[1]-i] + E_M2_rem[i+1] + P->MLclosing)){
2871                 backtrack_m(1,i,-1,-1,structure,vars);
2872                 backtrack_m2(i+1,-1,-1,structure,vars);
2873                 return;
2874               }
2875             }
2876           }
2877           if(E_M2_rem[i+1] != INF){
2878             if(E_M[my_iindx[1]-i])
2879               for(cnt1 = k_min_values_m[my_iindx[1]-i];
2880                   cnt1 <= k_max_values_m[my_iindx[1]-i];
2881                   cnt1++)
2882                 for(cnt2 = l_min_values_m[my_iindx[1]-i][cnt1];
2883                     cnt2 <= l_max_values_m[my_iindx[1]-i][cnt1];
2884                     cnt2 += 2)
2885                   if(E_Fc_rem == (E_M[my_iindx[1]-i][cnt1][cnt2/2] + E_M2_rem[i+1] + P->MLclosing)){
2886                     backtrack_m(1,i,cnt1,cnt2,structure,vars);
2887                     backtrack_m2(i+1,-1,-1,structure,vars);
2888                     return;
2889                   }
2890           }
2891
2892           if(!(E_M[my_iindx[1]-i])) continue;
2893           if(!(E_M2[i+1])) continue;
2894
2895           d1 = base_d1 - referenceBPs1[my_iindx[1]-i] - referenceBPs1[my_iindx[i+1]-seq_length];
2896           d2 = base_d2 - referenceBPs2[my_iindx[1]-i] - referenceBPs2[my_iindx[i+1]-seq_length];
2897           for(cnt1 = k_min_values_m[my_iindx[1]-i];
2898               cnt1 <= k_max_values_m[my_iindx[1]-i];
2899               cnt1++)
2900             for(cnt2 = l_min_values_m[my_iindx[1]-i][cnt1];
2901                 cnt2 <= l_max_values_m[my_iindx[1]-i][cnt1];
2902                 cnt2 += 2)
2903               for(cnt3 = k_min_values_m2[i+1];
2904                   cnt3 <= k_max_values_m2[i+1];
2905                   cnt3++)
2906                 for(cnt4 = l_min_values_m2[i+1][cnt3];
2907                     cnt4 <= l_max_values_m2[i+1][cnt3];
2908                     cnt4 += 2)
2909                   if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){
2910                     if(E_Fc_rem == (E_M[my_iindx[1]-i][cnt1][cnt2/2] + E_M2[i+1][cnt3][cnt4/2] + P->MLclosing)){
2911                       backtrack_m(1, i, cnt1, cnt2, structure, vars);
2912                       backtrack_m2(i+1, cnt3, cnt4, structure, vars);
2913                       return;
2914                     }
2915                   }
2916         }
2917     }
2918   }
2919   else{
2920     /* open chain ? */
2921     if(E_Fc[k][l/2] == 0)
2922       if((k == referenceBPs1[my_iindx[1]-seq_length]) && (l == referenceBPs2[my_iindx[1]-seq_length])){
2923         return;
2924       }
2925     if((k >= k_min_values_fcH) && (k <= k_max_values_fcH)){
2926       if((l >= l_min_values_fcH[k]) && (l <= l_max_values_fcH[k]))
2927         if(E_Fc[k][l/2] == E_FcH[k][l/2]){
2928           for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
2929             for (j = d; j <= seq_length; j++) {
2930               unsigned int u, ij;
2931               int type, no_close;
2932               char loopseq[10];
2933               i = j-d+1;
2934               ij = my_iindx[i]-j;
2935               if (!E_C[ij]) continue;
2936               u = seq_length-j + i-1;
2937               if (u<TURN) continue;
2938
2939               type = ptype[ij];
2940
2941               no_close = (((type==3)||(type==4))&&no_closingGU);
2942
2943               type=rtype[type];
2944
2945               if (!type) continue;
2946               if(no_close) continue;
2947
2948               d1 = base_d1 - referenceBPs1[ij];
2949               d2 = base_d2 - referenceBPs2[ij];
2950               if (u<7) {
2951                 strcpy(loopseq , sequence+j-1);
2952                 strncat(loopseq, sequence, i);
2953               }
2954               energy = E_Hairpin(u, type, S1[j+1], S1[i-1],  loopseq, P);
2955               if((k >= d1) && (l >= d2))
2956                 if((k-d1 >= k_min_values[ij]) && (k-d1 <= k_max_values[ij]))
2957                   if((l-d2 >= l_min_values[ij][k-d1]) && (l-d2 <= l_max_values[ij][k-d1])){
2958                     if(E_Fc[k][l/2] == E_C[ij][k-d1][(l-d2)/2] + energy){
2959                       backtrack_c(i, j, k-d1, l-d2, structure, vars);
2960                       return;
2961                     }
2962                   }
2963             }
2964         }
2965     }
2966
2967     if((k >= k_min_values_fcI) && (k <= k_max_values_fcI)){
2968       if((l >= l_min_values_fcI[k]) && (l <= l_max_values_fcI[k]))
2969         if(E_Fc[k][l/2] == E_FcI[k][l/2]){
2970           for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
2971             for (j = d; j <= seq_length; j++) {
2972               unsigned int u, ij, p, q, pq;
2973               int type, type_2;
2974               i = j-d+1;
2975               ij = my_iindx[i]-j;
2976               if(!E_C[ij]) continue;
2977               u = seq_length-j + i-1;
2978               if (u<TURN) continue;
2979
2980               type = ptype[ij];
2981
2982               type=rtype[type];
2983
2984               if (!type) continue;
2985
2986               for(p = j+1; p < seq_length ; p++){
2987                 unsigned int u1, qmin, ln_pre;
2988                 u1 = p-j-1;
2989                 if (u1+i-1>MAXLOOP) break;
2990                 qmin = p + TURN + 1;
2991                 ln_pre = u1 + i + seq_length;
2992                 if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1;
2993                 for(q = qmin; q <= seq_length; q++){
2994                   unsigned int u2;
2995                   pq = my_iindx[p]-q;
2996                   if(!E_C[pq]) continue;
2997                   type_2 = rtype[(unsigned int)ptype[pq]];
2998                   if (type_2==0) continue;
2999                   u2 = i-1 + seq_length-q;
3000                   if (u1+u2>MAXLOOP) continue;
3001                   /* get distance to reference if closing the interior loop
3002                   *  d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j})
3003                   *  d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j})
3004                   */
3005                   d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq];
3006                   d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq];
3007                   energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P);
3008                   if((k >= d1) && (l >= d2))
3009                     for(cnt1 = k_min_values[ij]; cnt1 <= MIN2(k_max_values[ij], k - d1); cnt1++)
3010                       for(cnt2 = l_min_values[ij][cnt1]; cnt2 <= MIN2(l_max_values[ij][cnt1], l - d2); cnt2+=2)
3011                         if((k - d1 - cnt1 >= k_min_values[pq]) && (k - d1 - cnt1 <= k_max_values[pq]))
3012                           if((l - d2 - cnt2 >= l_min_values[pq][k-d1-cnt1]) && (l - d2 - cnt2 <= l_max_values[pq][k-d1-cnt1])){
3013                             if((E_C[ij][cnt1][cnt2/2] + E_C[pq][k-d1-cnt1][(l-d2-cnt2)/2] + energy) == E_Fc[k][l/2]){
3014                               backtrack_c(i, j, cnt1, cnt2, structure, vars);
3015                               backtrack_c(p, q, k - d1 - cnt1, l - d2 - cnt2, structure, vars);
3016                               return;
3017                             }
3018                           }
3019                 }
3020               }
3021             }
3022         }
3023     }
3024
3025     if((k >= k_min_values_fcM) && (k <= k_max_values_fcM)){
3026       if((l >= l_min_values_fcM[k]) && (l <= l_max_values_fcM[k]))
3027         if(E_Fc[k][l/2] == E_FcM[k][l/2]){
3028           if(seq_length > 2*TURN)
3029             for (i=TURN+1; i<seq_length-2*TURN; i++) {
3030               /* get distancies to references
3031               * d3a = dbp(T1_[1,n}, T1_{1,k} + T1_{k+1, n})
3032               * d3b = dbp(T2_[1,n}, T2_{1,k} + T2_{k+1, n})
3033               */
3034               if(!E_M[my_iindx[1]-i]) continue;
3035               if(!E_M2[i+1]) continue;
3036               d1 = base_d1 - referenceBPs1[my_iindx[1]-i] - referenceBPs1[my_iindx[i+1]-seq_length];
3037               d2 = base_d2 - referenceBPs2[my_iindx[1]-i] - referenceBPs2[my_iindx[i+1]-seq_length];
3038               if((k >= d1) && (l >= d2))
3039                 for(cnt1 = k_min_values_m[my_iindx[1]-i]; cnt1 <= MIN2(k_max_values_m[my_iindx[1]-i], k-d1); cnt1++)
3040                   for(cnt2 = l_min_values_m[my_iindx[1]-i][cnt1]; cnt2 <= MIN2(l_max_values_m[my_iindx[1]-i][cnt1], l-d2); cnt2+=2)
3041                     if((k - d1 - cnt1 >= k_min_values_m2[i+1]) && (k - d1 - cnt1 <= k_max_values_m2[i+1]))
3042                       if((l - d2 - cnt2 >= l_min_values_m2[i+1][k-d1-cnt1]) && (l - d2 - cnt2 <= l_max_values_m2[i+1][k-d1-cnt1]))
3043                         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]){
3044                           backtrack_m(1, i, cnt1, cnt2, structure, vars);
3045                           backtrack_m2(i+1, k - d1 - cnt1, l - d2 - cnt2, structure, vars);
3046                           return;
3047                         }
3048             }
3049         }
3050     }
3051   }
3052   nrerror("backtack failed in fc\n");
3053 }
3054
3055
3056 PRIVATE void backtrack_m2(unsigned int i, int k, int l, char *structure, TwoDfold_vars *vars){
3057   unsigned int   j, ij, j3, n;
3058   unsigned int   *referenceBPs1, *referenceBPs2;
3059   unsigned int d1, d2, base_d1, base_d2, maxD1, maxD2;
3060   int *my_iindx, cnt1, cnt2, cnt3, cnt4;
3061   int ***E_M1, ***E_M2, *E_M2_rem, *E_M1_rem, e;
3062   int **l_min_values_m1, **l_max_values_m1, *k_min_values_m1, *k_max_values_m1;
3063
3064   n               = vars->seq_length;
3065   my_iindx        = vars->my_iindx;
3066   referenceBPs1   = vars->referenceBPs1;
3067   referenceBPs2   = vars->referenceBPs2;
3068
3069   E_M1            = vars->E_M1;
3070   l_min_values_m1 = vars->l_min_values_m1;
3071   l_max_values_m1 = vars->l_max_values_m1;
3072   k_min_values_m1 = vars->k_min_values_m1;
3073   k_max_values_m1 = vars->k_max_values_m1;
3074
3075   E_M1_rem        = vars->E_M1_rem;
3076
3077   E_M2            = vars->E_M2;
3078
3079   E_M2_rem        = vars->E_M2_rem;
3080
3081   maxD1           = vars->maxD1;
3082   maxD2           = vars->maxD2;
3083
3084   base_d1 = referenceBPs1[my_iindx[i]-n];
3085   base_d2 = referenceBPs2[my_iindx[i]-n];
3086
3087   if(k == -1){
3088     e = E_M2_rem[i];
3089     for (j=i+TURN+1; j<n-TURN-1; j++){
3090       if(E_M1_rem[my_iindx[i]-j] != INF){
3091         if(E_M1[my_iindx[j+1]-n])
3092           for(cnt1 = k_min_values_m1[my_iindx[j+1]-n];
3093               cnt1 <= k_max_values_m1[my_iindx[j+1]-n];
3094               cnt1++)
3095             for(cnt2 = l_min_values_m1[my_iindx[j+1]-n][cnt1];
3096                 cnt2 <= l_max_values_m1[my_iindx[j+1]-n][cnt1];
3097                 cnt2++)
3098               if(e == E_M1_rem[my_iindx[i]-j] + E_M1[my_iindx[j+1]-n][cnt1][cnt2/2]){
3099                 backtrack_m1(i, j, k, l, structure, vars);
3100                 backtrack_m1(j+1, n, cnt1, cnt2, structure, vars);
3101                 return;
3102               }
3103         if(E_M1_rem[my_iindx[j+1]-n] != INF){
3104           if(e == E_M1_rem[my_iindx[i]-j] + E_M1_rem[my_iindx[j+1]-n]){
3105             backtrack_m1(i, j, k, l, structure, vars);
3106             backtrack_m1(j+1, n, k, l, structure, vars);
3107             return;
3108           }
3109         }
3110       }
3111       if(E_M1_rem[my_iindx[j+1]-n] != INF){
3112         if(E_M1[my_iindx[i]-j])
3113           for(cnt1 = k_min_values_m1[my_iindx[i]-j];
3114               cnt1 <= k_max_values_m1[my_iindx[i]-j];
3115               cnt1++)
3116             for(cnt2 = l_min_values_m1[my_iindx[i]-j][cnt1];
3117                 cnt2 <= l_max_values_m1[my_iindx[i]-j][cnt1];
3118                 cnt2 += 2)
3119               if(e == E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1_rem[my_iindx[j+1]-n]){
3120                 backtrack_m1(i, j, cnt1, cnt2, structure, vars);
3121                 backtrack_m1(j+1, n, k, l, structure, vars);
3122                 return;
3123               }
3124       }
3125
3126
3127       if(!E_M1[my_iindx[i]-j]) continue;
3128       if(!E_M1[my_iindx[j+1]-n]) continue;
3129
3130       d1 = referenceBPs1[my_iindx[i]-n] - referenceBPs1[my_iindx[i]-j] - referenceBPs1[my_iindx[j+1]-n];
3131       d2 = referenceBPs2[my_iindx[i]-n] - referenceBPs2[my_iindx[i]-j] - referenceBPs2[my_iindx[j+1]-n];
3132
3133       for(cnt1 = k_min_values_m1[my_iindx[i]-j]; cnt1 <= k_max_values_m1[my_iindx[i]-j]; cnt1++)
3134         for(cnt2 = l_min_values_m1[my_iindx[i]-j][cnt1]; cnt2 <= l_max_values_m1[my_iindx[i]-j][cnt1]; cnt2+=2){
3135           for(cnt3 = k_min_values_m1[my_iindx[j+1]-n]; cnt3 <= k_max_values_m1[my_iindx[j+1]-n]; cnt3++)
3136             for(cnt4 = l_min_values_m1[my_iindx[j+1]-n][cnt3]; cnt4 <= l_max_values_m1[my_iindx[j+1]-n][cnt3]; cnt4+=2){
3137               if(((cnt1 + cnt3 + d1) > maxD1) || ((cnt2 + cnt4 + d2) > maxD2)){
3138                 if(e == E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-n][cnt3][cnt4/2]){
3139                   backtrack_m1(i, j, cnt1, cnt2, structure, vars);
3140                   backtrack_m1(j+1, n, cnt3, cnt4, structure, vars);
3141                   return;
3142                 }
3143               }
3144             }
3145         }
3146     }
3147   }
3148   else{
3149     for(j=i+TURN+1; j<n-TURN-1; j++){
3150       if(!E_M1[my_iindx[i]-j]) continue;
3151       if(!E_M1[my_iindx[j+1]-n]) continue;
3152
3153       ij = my_iindx[i]-j;
3154       j3 = my_iindx[j+1]-n;
3155       d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[j3];
3156       d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[j3];
3157
3158       for(cnt1 = k_min_values_m1[ij]; cnt1 <= MIN2(k_max_values_m1[ij], k - d1); cnt1++)
3159         for(cnt2 = l_min_values_m1[ij][cnt1]; cnt2 <= MIN2(l_max_values_m1[ij][cnt1], l-d2); cnt2+=2)
3160           if((k - d1 - cnt1 >= k_min_values_m1[j3]) && (k - d1 - cnt1 <= k_max_values_m1[j3]))
3161             if((l - d2 - cnt2 >= l_min_values_m1[j3][k - d1 - cnt1]) && (l - d2 - cnt2 <= l_max_values_m1[j3][k-d1-cnt1]))
3162               if(E_M1[ij][cnt1][cnt2/2] + E_M1[j3][k-d1-cnt1][(l-d2-cnt2)/2] == E_M2[i][k][l/2]){
3163                 backtrack_m1(i, j, cnt1, cnt2, structure, vars);
3164                 backtrack_m1(j+1, n, k-d1-cnt1, l-d2-cnt2, structure, vars);
3165                 return;
3166               }
3167     }
3168   }
3169   nrerror("backtack failed in m2\n");
3170 }
3171
3172 PRIVATE void mfe_circ(TwoDfold_vars *vars){
3173   unsigned int  d, i, j, maxD1, maxD2, seq_length, *referenceBPs1, *referenceBPs2, d1, d2, base_d1, base_d2, *mm1, *mm2, *bpdist;
3174   int           *my_iindx, energy, cnt1, cnt2, cnt3, cnt4;
3175   short         *S1;
3176   char          *sequence, *ptype;
3177   int           ***E_C, ***E_M, ***E_M1;
3178   int           *E_C_rem, *E_M_rem, *E_M1_rem;
3179   int           **l_min_values, **l_max_values, **l_min_values_m, **l_max_values_m, **l_min_values_m1, **l_max_values_m1;
3180   int           *k_min_values, *k_max_values,*k_min_values_m, *k_max_values_m,*k_min_values_m1, *k_max_values_m1;
3181
3182   paramT        *P;
3183
3184   P               = vars->P;
3185   sequence        = vars->sequence;
3186   seq_length      = vars->seq_length;
3187   maxD1           = vars->maxD1;
3188   maxD2           = vars->maxD2;
3189   S1              = vars->S1;
3190   ptype           = vars->ptype;
3191   my_iindx        = vars->my_iindx;
3192   referenceBPs1   = vars->referenceBPs1;
3193   referenceBPs2   = vars->referenceBPs2;
3194   mm1             = vars->mm1;
3195   mm2             = vars->mm2;
3196   bpdist          = vars->bpdist;
3197
3198   E_C             = vars->E_C;
3199   l_min_values    = vars->l_min_values;
3200   l_max_values    = vars->l_max_values;
3201   k_min_values    = vars->k_min_values;
3202   k_max_values    = vars->k_max_values;
3203
3204   E_M             = vars->E_M;
3205   l_min_values_m  = vars->l_min_values_m;
3206   l_max_values_m  = vars->l_max_values_m;
3207   k_min_values_m  = vars->k_min_values_m;
3208   k_max_values_m  = vars->k_max_values_m;
3209
3210   E_M1            = vars->E_M1;
3211   l_min_values_m1 = vars->l_min_values_m1;
3212   l_max_values_m1 = vars->l_max_values_m1;
3213   k_min_values_m1 = vars->k_min_values_m1;
3214   k_max_values_m1 = vars->k_max_values_m1;
3215
3216   E_C_rem        = vars->E_C_rem;
3217   E_M_rem        = vars->E_M_rem;
3218   E_M1_rem       = vars->E_M1_rem;
3219
3220 #ifdef _OPENMP
3221   #pragma omp parallel for private(d1,d2,cnt1,cnt2,cnt3,cnt4,j, i)
3222 #endif
3223   for(i=1; i<seq_length-TURN-1; i++){
3224     /* guess memory requirements for M2 */
3225
3226     int min_k, max_k, max_l, min_l;
3227     int min_k_real, max_k_real, *min_l_real, *max_l_real;
3228
3229     min_k = min_l = 0;
3230     max_k = mm1[my_iindx[i]-seq_length] + referenceBPs1[my_iindx[i] - seq_length];
3231     max_l = mm2[my_iindx[i]-seq_length] + referenceBPs2[my_iindx[i] - seq_length];
3232
3233     prepareBoundaries(min_k,
3234                       max_k,
3235                       min_l,
3236                       max_l,
3237                       bpdist[my_iindx[i] - seq_length],
3238                       &vars->k_min_values_m2[i],
3239                       &vars->k_max_values_m2[i],
3240                       &vars->l_min_values_m2[i],
3241                       &vars->l_max_values_m2[i]
3242                       );
3243
3244     prepareArray( &vars->E_M2[i],
3245                   vars->k_min_values_m2[i],
3246                   vars->k_max_values_m2[i],
3247                   vars->l_min_values_m2[i],
3248                   vars->l_max_values_m2[i]
3249                  );
3250
3251     preparePosteriorBoundaries( vars->k_max_values_m2[i] - vars->k_min_values_m2[i] + 1,
3252                                 vars->k_min_values_m2[i],
3253                                 &min_k_real,
3254                                 &max_k_real,
3255                                 &min_l_real,
3256                                 &max_l_real
3257                               );
3258
3259     /* begin filling of M2 array */
3260     for (j=i+TURN+1; j<seq_length-TURN-1; j++){
3261       if(E_M1_rem[my_iindx[i]-j] != INF){
3262         if(E_M1[my_iindx[j+1]-seq_length])
3263           for(cnt1 = k_min_values_m1[my_iindx[j+1]-seq_length];
3264               cnt1 <= k_max_values_m1[my_iindx[j+1]-seq_length];
3265               cnt1++)
3266             for(cnt2 = l_min_values_m1[my_iindx[j+1]-seq_length][cnt1];
3267                 cnt2 <= l_max_values_m1[my_iindx[j+1]-seq_length][cnt1];
3268                 cnt2++)
3269               vars->E_M2_rem[i] = MIN2(vars->E_M2_rem[i],
3270                                   E_M1_rem[my_iindx[i]-j] + E_M1[my_iindx[j+1]-seq_length][cnt1][cnt2/2]
3271                                   );
3272         if(E_M1_rem[my_iindx[j+1]-seq_length] != INF)
3273           vars->E_M2_rem[i] = MIN2(vars->E_M2_rem[i], E_M1_rem[my_iindx[i]-j] + E_M1_rem[my_iindx[j+1]-seq_length]);
3274       }
3275       if(E_M1_rem[my_iindx[j+1]-seq_length] != INF){
3276         if(E_M1[my_iindx[i]-j])
3277           for(cnt1 = k_min_values_m1[my_iindx[i]-j];
3278               cnt1 <= k_max_values_m1[my_iindx[i]-j];
3279               cnt1++)
3280             for(cnt2 = l_min_values_m1[my_iindx[i]-j][cnt1];
3281                 cnt2 <= l_max_values_m1[my_iindx[i]-j][cnt1];
3282                 cnt2 += 2)
3283               vars->E_M2_rem[i] = MIN2(vars->E_M2_rem[i],
3284                                   E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1_rem[my_iindx[j+1]-seq_length]
3285                                   );
3286       }
3287
3288
3289       if(!E_M1[my_iindx[i]-j]) continue;
3290       if(!E_M1[my_iindx[j+1]-seq_length]) continue;
3291
3292       d1 = referenceBPs1[my_iindx[i]-seq_length] - referenceBPs1[my_iindx[i]-j] - referenceBPs1[my_iindx[j+1]-seq_length];
3293       d2 = referenceBPs2[my_iindx[i]-seq_length] - referenceBPs2[my_iindx[i]-j] - referenceBPs2[my_iindx[j+1]-seq_length];
3294
3295       for(cnt1 = k_min_values_m1[my_iindx[i]-j]; cnt1 <= k_max_values_m1[my_iindx[i]-j]; cnt1++)
3296         for(cnt2 = l_min_values_m1[my_iindx[i]-j][cnt1]; cnt2 <= l_max_values_m1[my_iindx[i]-j][cnt1]; cnt2+=2){
3297           for(cnt3 = k_min_values_m1[my_iindx[j+1]-seq_length]; cnt3 <= k_max_values_m1[my_iindx[j+1]-seq_length]; cnt3++)
3298             for(cnt4 = l_min_values_m1[my_iindx[j+1]-seq_length][cnt3]; cnt4 <= l_max_values_m1[my_iindx[j+1]-seq_length][cnt3]; cnt4+=2){
3299               if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){
3300                 vars->E_M2[i][cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2( vars->E_M2[i][cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2],
3301                                                                         E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-seq_length][cnt3][cnt4/2]
3302                                                                       );
3303                 updatePosteriorBoundaries(cnt1+cnt3+d1,
3304                                           cnt2+cnt4+d2,
3305                                           &min_k_real,
3306                                           &max_k_real,
3307                                           &min_l_real,
3308                                           &max_l_real
3309                                           );
3310               }
3311               else{
3312                 vars->E_M2_rem[i] = MIN2(vars->E_M2_rem[i],
3313                                     E_M1[my_iindx[i]-j][cnt1][cnt2/2] + E_M1[my_iindx[j+1]-seq_length][cnt3][cnt4/2]
3314                                     );
3315               }
3316             }
3317         }
3318     }
3319
3320     /* resize and move memory portions of energy matrix E_M2 */
3321     adjustArrayBoundaries(&vars->E_M2[i],
3322                           &vars->k_min_values_m2[i],
3323                           &vars->k_max_values_m2[i],
3324                           &vars->l_min_values_m2[i],
3325                           &vars->l_max_values_m2[i],
3326                           min_k_real,
3327                           max_k_real,
3328                           min_l_real,
3329                           max_l_real
3330                           );
3331   } /* end for i */
3332
3333   base_d1 = referenceBPs1[my_iindx[1]-seq_length];
3334   base_d2 = referenceBPs2[my_iindx[1]-seq_length];
3335
3336   /* guess memory requirements for E_FcH, E_FcI and E_FcM */
3337
3338   int min_k, max_k, max_l, min_l;
3339   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;
3340   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;
3341
3342   min_k = min_l = 0;
3343
3344   max_k = mm1[my_iindx[1] - seq_length] + referenceBPs1[my_iindx[1] - seq_length];
3345   max_l = mm2[my_iindx[1] - seq_length] + referenceBPs2[my_iindx[1] - seq_length];
3346
3347 #ifdef _OPENMP
3348   #pragma omp sections
3349   {
3350
3351   #pragma omp section
3352   {
3353 #endif
3354   prepareBoundaries(min_k,
3355                     max_k,
3356                     min_l,
3357                     max_l,
3358                     bpdist[my_iindx[1] - seq_length],
3359                     &vars->k_min_values_fc,
3360                     &vars->k_max_values_fc,
3361                     &vars->l_min_values_fc,
3362                     &vars->l_max_values_fc
3363                     );
3364   prepareArray( &vars->E_Fc,
3365                 vars->k_min_values_fc,
3366                 vars->k_max_values_fc,
3367                 vars->l_min_values_fc,
3368                 vars->l_max_values_fc
3369               );
3370 #ifdef _OPENMP
3371   }
3372   #pragma omp section
3373   {
3374 #endif
3375   prepareBoundaries(min_k,
3376                     max_k,
3377                     min_l,
3378                     max_l,
3379                     bpdist[my_iindx[1] - seq_length],
3380                     &vars->k_min_values_fcH,
3381                     &vars->k_max_values_fcH,
3382                     &vars->l_min_values_fcH,
3383                     &vars->l_max_values_fcH
3384                     );
3385   prepareArray( &vars->E_FcH,
3386                 vars->k_min_values_fcH,
3387                 vars->k_max_values_fcH,
3388                 vars->l_min_values_fcH,
3389                 vars->l_max_values_fcH
3390               );
3391 #ifdef _OPENMP
3392   }
3393   #pragma omp section
3394   {
3395 #endif
3396   prepareBoundaries(min_k,
3397                     max_k,
3398                     min_l,
3399                     max_l,
3400                     bpdist[my_iindx[1] - seq_length],
3401                     &vars->k_min_values_fcI,
3402                     &vars->k_max_values_fcI,
3403                     &vars->l_min_values_fcI,
3404                     &vars->l_max_values_fcI
3405                     );
3406   prepareArray( &vars->E_FcI,
3407                 vars->k_min_values_fcI,
3408                 vars->k_max_values_fcI,
3409                 vars->l_min_values_fcI,
3410                 vars->l_max_values_fcI
3411               );
3412 #ifdef _OPENMP
3413   }
3414   #pragma omp section
3415   {
3416 #endif
3417   prepareBoundaries(min_k,
3418                     max_k,
3419                     min_l,
3420                     max_l,
3421                     bpdist[my_iindx[1] - seq_length],
3422                     &vars->k_min_values_fcM,
3423                     &vars->k_max_values_fcM,
3424                     &vars->l_min_values_fcM,
3425                     &vars->l_max_values_fcM
3426                     );
3427   prepareArray( &vars->E_FcM,
3428                 vars->k_min_values_fcM,
3429                 vars->k_max_values_fcM,
3430                 vars->l_min_values_fcM,
3431                 vars->l_max_values_fcM
3432               );
3433 #ifdef _OPENMP
3434   }
3435   #pragma omp section
3436   {
3437 #endif
3438   preparePosteriorBoundaries( max_k - min_k + 1,
3439                               min_k,
3440                               &min_k_real,
3441                               &max_k_real,
3442                               &min_l_real,
3443                               &max_l_real
3444                             );
3445 #ifdef _OPENMP
3446   }
3447   #pragma omp section
3448   {
3449 #endif
3450   preparePosteriorBoundaries( max_k - min_k + 1,
3451                               min_k,
3452                               &min_k_real_fcH,
3453                               &max_k_real_fcH,
3454                               &min_l_real_fcH,
3455                               &max_l_real_fcH
3456                             );
3457 #ifdef _OPENMP
3458   }
3459   #pragma omp section
3460   {
3461 #endif
3462   preparePosteriorBoundaries( max_k - min_k + 1,
3463                               min_k,
3464                               &min_k_real_fcI,
3465                               &max_k_real_fcI,
3466                               &min_l_real_fcI,
3467                               &max_l_real_fcI
3468                             );
3469
3470 #ifdef _OPENMP
3471   }
3472   #pragma omp section
3473   {
3474 #endif
3475   preparePosteriorBoundaries( max_k - min_k + 1,
3476                               min_k,
3477                               &min_k_real_fcM,
3478                               &max_k_real_fcM,
3479                               &min_l_real_fcM,
3480                               &max_l_real_fcM
3481                             );
3482 #ifdef _OPENMP
3483   }
3484   }
3485 #endif
3486
3487   /* begin actual energy calculations */
3488 #ifdef _OPENMP
3489   #pragma omp sections private(d, d1,d2,cnt1,cnt2,cnt3,cnt4,j, i, energy)
3490   {
3491
3492   #pragma omp section
3493   {
3494 #endif
3495   for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
3496     for (j = d; j <= seq_length; j++) {
3497       unsigned int  u, ij;
3498       int           type, no_close;
3499       char          loopseq[10];
3500       i = j-d+1;
3501       ij = my_iindx[i]-j;
3502       u = seq_length-j + i-1;
3503       if (u<TURN) continue;
3504
3505       type = ptype[ij];
3506
3507       no_close = (((type==3)||(type==4))&&no_closingGU);
3508
3509       type=rtype[type];
3510
3511       if (!type) continue;
3512       if(no_close) continue;
3513
3514       d1 = base_d1 - referenceBPs1[ij];
3515       d2 = base_d2 - referenceBPs2[ij];
3516       if (u<7) {
3517         strcpy(loopseq , sequence+j-1);
3518         strncat(loopseq, sequence, i);
3519       }
3520       energy = E_Hairpin(u, type, S1[j+1], S1[i-1],  loopseq, P);
3521
3522       if(E_C_rem[ij] != INF)
3523         vars->E_FcH_rem = MIN2(vars->E_FcH_rem, E_C_rem[ij] + energy);
3524
3525       if (!E_C[ij]) continue;
3526       for(cnt1 = k_min_values[ij]; cnt1 <= k_max_values[ij]; cnt1++)
3527         for(cnt2 = l_min_values[ij][cnt1]; cnt2 <= l_max_values[ij][cnt1]; cnt2 += 2){
3528           if(((cnt1 + d1) <= maxD1) && ((cnt2 + d2) <= maxD2)){
3529             vars->E_FcH[cnt1 + d1][(cnt2+d2)/2] = MIN2( vars->E_FcH[cnt1 + d1][(cnt2+d2)/2],
3530                                                   energy + E_C[ij][cnt1][cnt2/2]
3531                                                 );
3532             updatePosteriorBoundaries(cnt1 + d1,
3533                                       cnt2 + d2,
3534                                       &min_k_real_fcH,
3535                                       &max_k_real_fcH,
3536                                       &min_l_real_fcH,
3537                                       &max_l_real_fcH
3538                                     );
3539           }
3540           else
3541             vars->E_FcH_rem = MIN2(vars->E_FcH_rem, energy + E_C[ij][cnt1][cnt2/2]);
3542         }
3543     }
3544   /* end of i-j loop */
3545
3546   /* resize and move memory portions of energy matrix E_FcH */
3547   adjustArrayBoundaries(&vars->E_FcH,
3548                         &vars->k_min_values_fcH,
3549                         &vars->k_max_values_fcH,
3550                         &vars->l_min_values_fcH,
3551                         &vars->l_max_values_fcH,
3552                         min_k_real_fcH,
3553                         max_k_real_fcH,
3554                         min_l_real_fcH,
3555                         max_l_real_fcH
3556                       );
3557 #ifdef _OPENMP
3558   }
3559   #pragma omp section
3560   {
3561 #endif
3562   for (d = TURN+2; d <= seq_length; d++) /* i,j in [1..length] */
3563     for (j = d; j <= seq_length; j++) {
3564       unsigned int  u, ij, p, q, pq;
3565       int           type, type_2, no_close;
3566       i = j-d+1;
3567       ij = my_iindx[i]-j;
3568       u = seq_length-j + i-1;
3569       if (u<TURN) continue;
3570
3571       type = ptype[ij];
3572
3573       no_close = (((type==3)||(type==4))&&no_closingGU);
3574
3575       type=rtype[type];
3576
3577       if (!type) continue;
3578       if(no_close) continue;
3579
3580       if(E_C_rem[ij] != INF){
3581         for(p = j+1; p < seq_length ; p++){
3582           unsigned int u1, qmin, ln_pre;
3583           u1 = p-j-1;
3584           if (u1+i-1>MAXLOOP) break;
3585           qmin = p + TURN + 1;
3586           ln_pre = u1 + i + seq_length;
3587           if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1;
3588           for(q = qmin; q <= seq_length; q++){
3589             unsigned int u2;
3590             pq = my_iindx[p]-q;
3591             type_2 = rtype[(unsigned int)ptype[pq]];
3592             if (type_2==0) continue;
3593             u2 = i-1 + seq_length-q;
3594             if (u1+u2>MAXLOOP) continue;
3595             /* get distance to reference if closing the interior loop
3596             *  d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j})
3597             *  d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j})
3598             */
3599             d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq];
3600             d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq];
3601             energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P);
3602
3603             if(E_C_rem[pq] != INF)
3604               vars->E_FcI_rem = MIN2(vars->E_FcI_rem, E_C_rem[ij] + E_C_rem[pq] + energy);
3605
3606             if(E_C[pq])
3607               for(cnt1 = k_min_values[pq];
3608                   cnt1 <= k_max_values[pq];
3609                   cnt1++)
3610                 for(cnt2 = l_min_values[pq][cnt1];
3611                     cnt2 <= l_max_values[pq][cnt1];
3612                     cnt2 += 2)
3613                   vars->E_FcI_rem = MIN2(vars->E_FcI_rem, E_C_rem[ij] + E_C[pq][cnt1][cnt2/2] + energy);
3614           }
3615         }
3616       }
3617
3618       if(E_C[ij]){
3619         for(p = j+1; p < seq_length ; p++){
3620           unsigned int u1, qmin, ln_pre;
3621           u1 = p-j-1;
3622           if (u1+i-1>MAXLOOP) break;
3623           qmin = p + TURN + 1;
3624           ln_pre = u1 + i + seq_length;
3625           if(ln_pre > qmin + MAXLOOP) qmin = ln_pre - MAXLOOP - 1;
3626           for(q = qmin; q <= seq_length; q++){
3627             unsigned int u2;
3628             pq = my_iindx[p]-q;
3629             type_2 = rtype[(unsigned int)ptype[pq]];
3630             if (type_2==0) continue;
3631             u2 = i-1 + seq_length-q;
3632             if (u1+u2>MAXLOOP) continue;
3633             /* get distance to reference if closing the interior loop
3634             *  d2a = dbp(T1_[1,n}, T1_{p,q} + T1_{i,j})
3635             *  d2b = dbp(T2_[1,n}, T2_{p,q} + T2_{i,j})
3636             */
3637             d1 = base_d1 - referenceBPs1[ij] - referenceBPs1[pq];
3638             d2 = base_d2 - referenceBPs2[ij] - referenceBPs2[pq];
3639             energy = E_IntLoop(u1, u2, type, type_2, S1[j+1], S1[i-1], S1[p-1], S1[q+1], P);
3640             if(E_C_rem[pq] != INF){
3641               for(cnt1 = k_min_values[ij];
3642                   cnt1 <= k_max_values[ij];
3643                   cnt1++)
3644                 for(cnt2 = l_min_values[ij][cnt1];
3645                     cnt2 <= l_max_values[ij][cnt1];
3646                     cnt2 += 2)
3647                   vars->E_FcI_rem = MIN2(vars->E_FcI_rem, E_C[ij][cnt1][cnt2/2] + E_C_rem[pq] + energy);
3648             }
3649
3650             if(E_C[pq])
3651               for(cnt1 = k_min_values[ij];
3652                   cnt1 <= k_max_values[ij];
3653                   cnt1++)
3654                 for(cnt2 = l_min_values[ij][cnt1];
3655                     cnt2 <= l_max_values[ij][cnt1];
3656                     cnt2 += 2)
3657                   for(cnt3 = k_min_values[pq];
3658                       cnt3 <= k_max_values[pq];
3659                       cnt3++)
3660                     for(cnt4 = l_min_values[pq][cnt3];
3661                         cnt4 <= l_max_values[pq][cnt3];
3662                         cnt4 += 2){
3663                       if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){
3664                         vars->E_FcI[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2(
3665                                                                           vars->E_FcI[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2],
3666                                                                           E_C[ij][cnt1][cnt2/2]
3667                                                                         + E_C[pq][cnt3][cnt4/2]
3668                                                                         + energy
3669                                                                             );
3670                         updatePosteriorBoundaries(cnt1 + cnt3 + d1,
3671                                                   cnt2 + cnt4 + d2,
3672                                                   &min_k_real_fcI,
3673                                                   &max_k_real_fcI,
3674                                                   &min_l_real_fcI,
3675                                                   &max_l_real_fcI
3676                                                 );
3677                       }
3678                       else{
3679                         vars->E_FcI_rem = MIN2(
3680                                       vars->E_FcI_rem,
3681                                       E_C[ij][cnt1][cnt2/2]
3682                                     + E_C[pq][cnt3][cnt4/2]
3683                                     + energy
3684                                           );
3685                       }
3686                     }
3687           }
3688         }
3689       }
3690     }
3691   /* end of i-j loop */
3692
3693   /* resize and move memory portions of energy matrix E_FcI */
3694   adjustArrayBoundaries(&vars->E_FcI,
3695                         &vars->k_min_values_fcI,
3696                         &vars->k_max_values_fcI,
3697                         &vars->l_min_values_fcI,
3698                         &vars->l_max_values_fcI,
3699                         min_k_real_fcI,
3700                         max_k_real_fcI,
3701                         min_l_real_fcI,
3702                         max_l_real_fcI
3703                       );
3704 #ifdef _OPENMP
3705   }
3706   #pragma omp section
3707   {
3708 #endif
3709   if(seq_length > 2*TURN){
3710     for (i=TURN+1; i<seq_length-2*TURN; i++) {
3711       /* get distancies to references
3712       * d3a = dbp(T1_[1,n}, T1_{1,k} + T1_{k+1, n})
3713       * d3b = dbp(T2_[1,n}, T2_{1,k} + T2_{k+1, n})
3714       */
3715       d1 = base_d1 - referenceBPs1[my_iindx[1]-i] - referenceBPs1[my_iindx[i+1]-seq_length];
3716       d2 = base_d2 - referenceBPs2[my_iindx[1]-i] - referenceBPs2[my_iindx[i+1]-seq_length];
3717
3718       if(E_M_rem[my_iindx[1]-i] != INF){
3719         if(vars->E_M2[i+1])
3720           for(cnt1 = vars->k_min_values_m2[i+1];
3721               cnt1 <= vars->k_max_values_m2[i+1];
3722               cnt1++)
3723             for(cnt2 = vars->l_min_values_m2[i+1][cnt1];
3724                 cnt2 <= vars->l_max_values_m2[i+1][cnt1];
3725                 cnt2 += 2)
3726               vars->E_FcM_rem = MIN2(vars->E_FcM_rem, E_M_rem[my_iindx[1]-i] + vars->E_M2[i+1][cnt1][cnt2/2] + P->MLclosing);
3727         if(vars->E_M2_rem[i+1] != INF)
3728           vars->E_FcM_rem = MIN2(vars->E_FcM_rem, E_M_rem[my_iindx[1]-i] + vars->E_M2_rem[i+1] + P->MLclosing);
3729       }
3730       if(vars->E_M2_rem[i+1] != INF){
3731         if(E_M[my_iindx[1]-i])
3732           for(cnt1 = k_min_values_m[my_iindx[1]-i];
3733               cnt1 <= k_max_values_m[my_iindx[1]-i];
3734               cnt1++)
3735             for(cnt2 = l_min_values_m[my_iindx[1]-i][cnt1];
3736                 cnt2 <= l_max_values_m[my_iindx[1]-i][cnt1];
3737                 cnt2 += 2)
3738               vars->E_FcM_rem = MIN2(vars->E_FcM_rem, E_M[my_iindx[1]-i][cnt1][cnt2/2] + vars->E_M2_rem[i+1] + P->MLclosing);
3739       }
3740
3741       if(!E_M[my_iindx[1]-i]) continue;
3742       if(!vars->E_M2[i+1]) continue;
3743       for(cnt1 = k_min_values_m[my_iindx[1]-i]; cnt1 <= k_max_values_m[my_iindx[1]-i]; cnt1++)
3744         for(cnt2 = l_min_values_m[my_iindx[1]-i][cnt1]; cnt2 <= l_max_values_m[my_iindx[1]-i][cnt1]; cnt2 += 2)
3745           for(cnt3 = vars->k_min_values_m2[i+1]; cnt3 <= vars->k_max_values_m2[i+1]; cnt3++)
3746             for(cnt4 = vars->l_min_values_m2[i+1][cnt3]; cnt4 <= vars->l_max_values_m2[i+1][cnt3]; cnt4 += 2){
3747               if(((cnt1 + cnt3 + d1) <= maxD1) && ((cnt2 + cnt4 + d2) <= maxD2)){
3748                 vars->E_FcM[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2] = MIN2(
3749                                                                   vars->E_FcM[cnt1 + cnt3 + d1][(cnt2 + cnt4 + d2)/2],
3750                                                                   E_M[my_iindx[1]-i][cnt1][cnt2/2]
3751                                                                 + vars->E_M2[i+1][cnt3][cnt4/2]
3752                                                                 + P->MLclosing
3753                                                                     );
3754                 updatePosteriorBoundaries(cnt1 + cnt3 + d1,
3755                                           cnt2 + cnt4 + d2,
3756                                           &min_k_real_fcM,
3757                                           &max_k_real_fcM,
3758                                           &min_l_real_fcM,
3759                                           &max_l_real_fcM
3760                                         );
3761               }
3762               else{
3763                 vars->E_FcM_rem = MIN2(
3764                               vars->E_FcM_rem,
3765                               E_M[my_iindx[1]-i][cnt1][cnt2/2]
3766                             + vars->E_M2[i+1][cnt3][cnt4/2]
3767                             + P->MLclosing
3768                                   );
3769               }
3770             }
3771     }
3772   }
3773   /* resize and move memory portions of energy matrix E_FcM */
3774   adjustArrayBoundaries(&vars->E_FcM,
3775                         &vars->k_min_values_fcM,
3776                         &vars->k_max_values_fcM,
3777                         &vars->l_min_values_fcM,
3778                         &vars->l_max_values_fcM,
3779                         min_k_real_fcM,
3780                         max_k_real_fcM,
3781                         min_l_real_fcM,
3782                         max_l_real_fcM
3783                       );
3784 #ifdef _OPENMP
3785   }
3786   }
3787 #endif
3788
3789
3790
3791   /* compute E_Fc_rem */
3792   vars->E_Fc_rem = MIN2(vars->E_FcH_rem, vars->E_FcI_rem);
3793   vars->E_Fc_rem = MIN2(vars->E_Fc_rem, vars->E_FcM_rem);
3794   /* add the case were structure is unfolded chain */
3795   if((referenceBPs1[my_iindx[1]-seq_length] > maxD1) || (referenceBPs2[my_iindx[1]-seq_length] > maxD2))
3796     vars->E_Fc_rem = MIN2(vars->E_Fc_rem, 0);
3797
3798
3799   /* compute all E_Fc */
3800   for(cnt1 = vars->k_min_values_fcH; cnt1 <= vars->k_max_values_fcH; cnt1++)
3801     for(cnt2 = vars->l_min_values_fcH[cnt1]; cnt2 <= vars->l_max_values_fcH[cnt1]; cnt2 += 2){
3802       vars->E_Fc[cnt1][cnt2/2] = MIN2(vars->E_Fc[cnt1][cnt2/2],
3803                                       vars->E_FcH[cnt1][cnt2/2]
3804                                       );
3805       updatePosteriorBoundaries(cnt1,
3806                                 cnt2,
3807                                 &min_k_real,
3808                                 &max_k_real,
3809                                 &min_l_real,
3810                                 &max_l_real
3811                                 );
3812     }
3813   for(cnt1 = vars->k_min_values_fcI; cnt1 <= vars->k_max_values_fcI; cnt1++)
3814     for(cnt2 = vars->l_min_values_fcI[cnt1]; cnt2 <= vars->l_max_values_fcI[cnt1]; cnt2 += 2){
3815       vars->E_Fc[cnt1][cnt2/2] = MIN2(vars->E_Fc[cnt1][cnt2/2],
3816                                       vars->E_FcI[cnt1][cnt2/2]
3817                                       );
3818       updatePosteriorBoundaries(cnt1,
3819                                 cnt2,
3820                                 &min_k_real,
3821                                 &max_k_real,
3822                                 &min_l_real,
3823                                 &max_l_real
3824                                 );
3825     }
3826   for(cnt1 = vars->k_min_values_fcM; cnt1 <= vars->k_max_values_fcM; cnt1++)
3827     for(cnt2 = vars->l_min_values_fcM[cnt1]; cnt2 <= vars->l_max_values_fcM[cnt1]; cnt2 += 2){
3828       vars->E_Fc[cnt1][cnt2/2] = MIN2(vars->E_Fc[cnt1][cnt2/2],
3829                                       vars->E_FcM[cnt1][cnt2/2]
3830                                       );
3831       updatePosteriorBoundaries(cnt1,
3832                                 cnt2,
3833                                 &min_k_real,
3834                                 &max_k_real,
3835                                 &min_l_real,
3836                                 &max_l_real
3837                                 );
3838     }
3839   /* add the case were structure is unfolded chain */
3840   vars->E_Fc[referenceBPs1[my_iindx[1]-seq_length]][referenceBPs2[my_iindx[1]-seq_length]/2] = MIN2(vars->E_Fc[referenceBPs1[my_iindx[1]-seq_length]][referenceBPs2[my_iindx[1]-seq_length]/2],
3841                                                                                                     0);
3842   updatePosteriorBoundaries(referenceBPs1[my_iindx[1]-seq_length],
3843                             referenceBPs2[my_iindx[1]-seq_length],
3844                             &min_k_real,
3845                             &max_k_real,
3846                             &min_l_real,
3847                             &max_l_real
3848                             );
3849
3850
3851   adjustArrayBoundaries(&vars->E_Fc,
3852                         &vars->k_min_values_fc,
3853                         &vars->k_max_values_fc,
3854                         &vars->l_min_values_fc,
3855                         &vars->l_max_values_fc,
3856                         min_k_real,
3857                         max_k_real,
3858                         min_l_real,
3859                         max_l_real
3860                       );
3861
3862 }
3863
3864
3865
3866
3867 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){
3868   int cnt1;
3869   int k_diff_pre  = k_min_post - *k_min;
3870   int mem_size    = k_max_post - k_min_post + 1;
3871
3872   if(k_min_post < INF){
3873     /* free all the unused memory behind actual data */
3874     for(cnt1 = k_max_post + 1; cnt1 <= *k_max; cnt1++){
3875       (*array)[cnt1] += (*l_min)[cnt1]/2;
3876       free((*array)[cnt1]);
3877     }
3878
3879     /* free unused memory before actual data */
3880     for(cnt1 = *k_min; cnt1 < k_min_post; cnt1++){
3881       (*array)[cnt1] += (*l_min)[cnt1]/2;
3882       free((*array)[cnt1]);
3883     }
3884     /* move data to front and thereby eliminating unused memory in front of actual data */
3885     if(k_diff_pre > 0){
3886       memmove((int **)(*array),((int **)(*array)) + k_diff_pre, sizeof(int *) * mem_size);
3887       memmove((int *) (*l_min),((int *) (*l_min)) + k_diff_pre, sizeof(int)   * mem_size);
3888       memmove((int *) (*l_max),((int *) (*l_max)) + k_diff_pre, sizeof(int)   * mem_size);
3889     }
3890
3891     /* reallocating memory to actual size used */
3892     *array  += *k_min;
3893     *array = (int **)realloc(*array, sizeof(int *) * mem_size);
3894     *array -= k_min_post;
3895
3896     *l_min  += *k_min;
3897     *l_min = (int *)realloc(*l_min, sizeof(int) * mem_size);
3898     *l_min -= k_min_post;
3899
3900     *l_max  += *k_min;
3901     *l_max = (int *)realloc(*l_max, sizeof(int) * mem_size);
3902     *l_max -= k_min_post;
3903
3904     /* adjust l dimension of array */
3905     for(cnt1 = k_min_post; cnt1 <= k_max_post; cnt1++){
3906       if(l_min_post[cnt1] < INF){
3907         /* new memsize */
3908         mem_size        = (l_max_post[cnt1] - l_min_post[cnt1] + 1)/2 + 1;
3909         /* reshift the pointer */
3910         (*array)[cnt1]  += (*l_min)[cnt1]/2;
3911
3912         int shift       = (l_min_post[cnt1]%2 == (*l_min)[cnt1]%2) ? 0 : 1;
3913         /* eliminate unused memory in front of actual data */
3914         unsigned int    start = (l_min_post[cnt1] - (*l_min)[cnt1])/2 + shift;
3915         if(start > 0)
3916           memmove((int *)((*array)[cnt1]), (int *)((*array)[cnt1])+start, sizeof(int) * mem_size);
3917         (*array)[cnt1]  = (int *) realloc((*array)[cnt1], sizeof(int) * mem_size);
3918
3919         (*array)[cnt1]  -= l_min_post[cnt1]/2;
3920       }
3921       else{
3922         /* free according memory */
3923         (*array)[cnt1] += (*l_min)[cnt1]/2;
3924         free((*array)[cnt1]);
3925       }
3926
3927       (*l_min)[cnt1] = l_min_post[cnt1];
3928       (*l_max)[cnt1] = l_max_post[cnt1];
3929     }
3930   }
3931   else{
3932     /* we have to free all unused memory */
3933     for(cnt1 = *k_min; cnt1 <= *k_max; cnt1++){
3934       (*array)[cnt1] += (*l_min)[cnt1]/2;
3935       free((*array)[cnt1]);
3936     }
3937     (*l_min) += *k_min;
3938     (*l_max) += *k_min;
3939     free(*l_min);
3940     free(*l_max);
3941     (*array) += *k_min;
3942     free(*array);
3943     *array = NULL;
3944   }
3945
3946   l_min_post  += *k_min;
3947   l_max_post  += *k_min;
3948   free(l_min_post);
3949   free(l_max_post);
3950   *k_min      = k_min_post;
3951   *k_max      = k_max_post;
3952 }
3953
3954 INLINE PRIVATE void preparePosteriorBoundaries(int size, int shift, int *min_k, int *max_k, int **min_l, int **max_l){
3955   int i;
3956   *min_k  = INF;
3957   *max_k  = 0;
3958
3959   *min_l  = (int *)space(sizeof(int) * size);
3960   *max_l  = (int *)space(sizeof(int) * size);
3961
3962   for(i = 0; i < size; i++){
3963     (*min_l)[i] = INF;
3964     (*max_l)[i] = 0;
3965   }
3966
3967   *min_l  -= shift;
3968   *max_l  -= shift;
3969 }
3970
3971 INLINE PRIVATE void updatePosteriorBoundaries(int d1, int d2, int *min_k, int *max_k, int **min_l, int **max_l){
3972   (*min_l)[d1]  = MIN2((*min_l)[d1], d2);
3973   (*max_l)[d1]  = MAX2((*max_l)[d1], d2);
3974   *min_k        = MIN2(*min_k, d1);
3975   *max_k        = MAX2(*max_k, d1);
3976 }
3977
3978 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){
3979   int cnt;
3980   int mem = max_k_pre - min_k_pre + 1;
3981
3982   *min_k  = min_k_pre;
3983   *max_k  = max_k_pre;
3984   *min_l  = (int *) space(sizeof(int) * mem);
3985   *max_l  = (int *) space(sizeof(int) * mem);
3986
3987   *min_l  -= min_k_pre;
3988   *max_l  -= min_k_pre;
3989
3990   /* for each k guess the according minimum l*/
3991   for(cnt = min_k_pre; cnt <= max_k_pre; cnt++){
3992     (*min_l)[cnt] = min_l_pre;
3993     (*max_l)[cnt] = max_l_pre;
3994     while((*min_l)[cnt] + cnt < bpdist) (*min_l)[cnt]++;
3995     if((bpdist % 2) != (((*min_l)[cnt] + cnt) % 2)) (*min_l)[cnt]++;
3996   }
3997 }
3998
3999 INLINE  PRIVATE void  prepareArray(int ***array, int min_k, int max_k, int *min_l, int *max_l){
4000   int i, j, mem;
4001   *array  = (int **)space(sizeof(int *) * (max_k - min_k + 1));
4002   *array  -= min_k;
4003
4004   for(i = min_k; i <= max_k; i++){
4005     mem = (max_l[i] - min_l[i] + 1)/2 + 1;
4006     (*array)[i] = (int *)space(sizeof(int) * mem);
4007     for(j = 0; j < mem; j++)
4008       (*array)[i][j] = INF;
4009     (*array)[i]  -= min_l[i]/2;
4010   }
4011 }
4012
4013 INLINE  PRIVATE void  prepareArray2(unsigned long ***array, int min_k, int max_k, int *min_l, int *max_l){
4014   int i, mem;
4015   *array  = (unsigned long **)space(sizeof(unsigned long *) * (max_k - min_k + 1));
4016   *array  -= min_k;
4017
4018   for(i = min_k; i <= max_k; i++){
4019     mem = (max_l[i] - min_l[i] + 1)/2 + 1;
4020     (*array)[i] = (unsigned long *)space(sizeof(unsigned long) * mem);
4021     (*array)[i] -= min_l[i]/2;
4022   }
4023 }
4024