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