binaries/src/ViennaRNA/lib/cofold.c

   1 /* Last changed Time-stamp: <2008-12-03 17:44:38 ivo> */
   2 /*
   3                   minimum free energy
   4                   RNA secondary structure prediction
   5
   6                   c Ivo Hofacker, Chrisoph Flamm
   7                   original implementation by
   8                   Walter Fontana
   9
  10                   Vienna RNA package
  11 */
  12
  13 #include <config.h>
  14 #include <stdio.h>
  15 #include <stdlib.h>
  16 #include <math.h>
  17 #include <ctype.h>
  18 #include <string.h>
  19 #include <limits.h>
  20
  21 #include "utils.h"
  22 #include "energy_par.h"
  23 #include "fold_vars.h"
  24 #include "pair_mat.h"
  25 #include "params.h"
  26 #include "subopt.h"
  27 #include "fold.h"
  28 #include "loop_energies.h"
  29 #include "gquad.h"
  30 #include "cofold.h"
  31
  32 #ifdef _OPENMP
  33 #include <omp.h>
  34 #endif
  35
  36 #define PAREN
  37
  38 #define STACK_BULGE1      1       /* stacking energies for bulges of size 1 */
  39 #define NEW_NINIO         1       /* new asymetry penalty */
  40 #define MAXSECTORS        500     /* dimension for a backtrack array */
  41 #define LOCALITY          0.      /* locality parameter for base-pairs */
  42 #undef TURN
  43 #define TURN              0       /* reset minimal base pair span for intermolecular pairings */
  44 #define TURN2             3       /* used by zukersubopt */
  45 #define SAME_STRAND(I,J)  (((I)>=cut_point)||((J)<cut_point))
  46
  47 /*
  48 #################################
  49 # GLOBAL VARIABLES              #
  50 #################################
  51 */
  52
  53
  54 /*
  55 #################################
  56 # PRIVATE VARIABLES             #
  57 #################################
  58 */
  59
  60 PRIVATE float   mfe1, mfe2;       /* minimum free energies of the monomers */
  61 PRIVATE int     *indx   = NULL;   /* index for moving in the triangle matrices c[] and fMl[]*/
  62 PRIVATE int     *c      = NULL;   /* energy array, given that i-j pair */
  63 PRIVATE int     *cc     = NULL;   /* linear array for calculating canonical structures */
  64 PRIVATE int     *cc1    = NULL;   /*   "     "        */
  65 PRIVATE int     *f5     = NULL;   /* energy of 5' end */
  66 PRIVATE int     *fc     = NULL;   /* energy from i to cutpoint (and vice versa if i>cut) */
  67 PRIVATE int     *fML    = NULL;   /* multi-loop auxiliary energy array */
  68 PRIVATE int     *fM1    = NULL;   /* second ML array, only for subopt */
  69 PRIVATE int     *Fmi    = NULL;   /* holds row i of fML (avoids jumps in memory) */
  70 PRIVATE int     *DMLi   = NULL;   /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
  71 PRIVATE int     *DMLi1  = NULL;   /*             MIN(fML[i+1,k]+fML[k+1,j])  */
  72 PRIVATE int     *DMLi2  = NULL;   /*             MIN(fML[i+2,k]+fML[k+1,j])  */
  73 PRIVATE char    *ptype  = NULL;   /* precomputed array of pair types */
  74 PRIVATE short   *S = NULL, *S1 = NULL;
  75 PRIVATE paramT  *P          = NULL;
  76 PRIVATE int     init_length = -1;
  77 PRIVATE int     zuker       = 0; /* Do Zuker style suboptimals? */
  78 PRIVATE sect    sector[MAXSECTORS];   /* stack for backtracking */
  79 PRIVATE int     length;
  80 PRIVATE bondT   *base_pair2 = NULL;
  81 PRIVATE int     *BP; /* contains the structure constrainsts: BP[i]
  82                         -1: | = base must be paired
  83                         -2: < = base must be paired with j<i
  84                         -3: > = base must be paired with j>i
  85                         -4: x = base must not pair
  86                         positive int: base is paired with int      */
  87 PRIVATE int     struct_constrained  = 0;
  88 PRIVATE int     with_gquad          = 0;
  89
  90 PRIVATE int     *ggg = NULL;    /* minimum free energies of the gquadruplexes */
  91
  92 #ifdef _OPENMP
  93
  94 #pragma omp threadprivate(mfe1, mfe2, indx, c, cc, cc1, f5, fc, fML, fM1, Fmi, DMLi, DMLi1, DMLi2,\
  95                           ptype, S, S1, P, zuker, sector, length, base_pair2, BP, struct_constrained,\
  96                           ggg, with_gquad)
  97
  98 #endif
  99
 100 /*
 101 #################################
 102 # PRIVATE FUNCTION DECLARATIONS #
 103 #################################
 104 */
 105
 106 PRIVATE void  init_cofold(int length, paramT *parameters);
 107 PRIVATE void  get_arrays(unsigned int size);
 108 /* PRIVATE void  scale_parameters(void); */
 109 PRIVATE void  make_ptypes(const short *S, const char *structure);
 110 PRIVATE void  backtrack(const char *sequence);
 111 PRIVATE int   fill_arrays(const char *sequence);
 112 PRIVATE void  free_end(int *array, int i, int start);
 113
 114 /*
 115 #################################
 116 # BEGIN OF FUNCTION DEFINITIONS #
 117 #################################
 118 */
 119
 120 /*--------------------------------------------------------------------------*/
 121 PRIVATE void init_cofold(int length, paramT *parameters){
 122
 123 #ifdef _OPENMP
 124 /* Explicitly turn off dynamic threads */
 125   omp_set_dynamic(0);
 126 #endif
 127
 128   if (length<1) nrerror("init_cofold: argument must be greater 0");
 129   free_co_arrays();
 130   get_arrays((unsigned) length);
 131   init_length=length;
 132
 133   indx = get_indx((unsigned) length);
 134
 135   update_cofold_params_par(parameters);
 136 }
 137
 138 /*--------------------------------------------------------------------------*/
 139
 140 PRIVATE void get_arrays(unsigned int size){
 141   if(size >= (unsigned int)sqrt((double)INT_MAX))
 142     nrerror("get_arrays@cofold.c: sequence length exceeds addressable range");
 143
 144   c     = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 145   fML   = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 146   if (uniq_ML)
 147     fM1    = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 148
 149   ptype = (char *) space(sizeof(char)*((size*(size+1))/2+2));
 150   f5    = (int *) space(sizeof(int)*(size+2));
 151   fc    = (int *) space(sizeof(int)*(size+2));
 152   cc    = (int *) space(sizeof(int)*(size+2));
 153   cc1   = (int *) space(sizeof(int)*(size+2));
 154   Fmi   = (int *) space(sizeof(int)*(size+1));
 155   DMLi  = (int *) space(sizeof(int)*(size+1));
 156   DMLi1 = (int *) space(sizeof(int)*(size+1));
 157   DMLi2 = (int *) space(sizeof(int)*(size+1));
 158
 159   base_pair2 = (bondT *) space(sizeof(bondT)*(1+size/2+200)); /* add a guess of how many G's may be involved in a G quadruplex */
 160 }
 161
 162 /*--------------------------------------------------------------------------*/
 163
 164 PUBLIC void free_co_arrays(void){
 165   if(indx)        free(indx);
 166   if(c)           free(c);
 167   if(fML)         free(fML);
 168   if(f5)          free(f5);
 169   if(cc)          free(cc);
 170   if(cc1)         free(cc1);
 171   if(fc)          free(fc);
 172   if(ptype)       free(ptype);
 173   if(fM1)         free(fM1);
 174   if(base_pair2)  free(base_pair2);
 175   if(Fmi)         free(Fmi);
 176   if(DMLi)        free(DMLi);
 177   if(DMLi1)       free(DMLi1);
 178   if(DMLi2)       free(DMLi2);
 179   if(P)           free(P);
 180   if(ggg)         free(ggg);
 181
 182   indx = c = fML = f5 = cc = cc1 = fc = fM1 = Fmi = DMLi = DMLi1 = DMLi2 = ggg = NULL;
 183   ptype       = NULL;
 184   base_pair2  = NULL;
 185   P           = NULL;
 186   init_length = 0;
 187 }
 188
 189
 190 /*--------------------------------------------------------------------------*/
 191
 192 PUBLIC void export_cofold_arrays_gq(  int **f5_p,
 193                                       int **c_p,
 194                                       int **fML_p,
 195                                       int **fM1_p,
 196                                       int **fc_p,
 197                                       int **ggg_p,
 198                                       int **indx_p,
 199                                       char **ptype_p){
 200
 201   /* make the DP arrays available to routines such as subopt() */
 202   *f5_p = f5; *c_p = c;
 203   *fML_p = fML; *fM1_p = fM1;
 204   *ggg_p = ggg;
 205   *indx_p = indx; *ptype_p = ptype;
 206   *fc_p =fc;
 207 }
 208
 209 PUBLIC void export_cofold_arrays( int **f5_p,
 210                                   int **c_p,
 211                                   int **fML_p,
 212                                   int **fM1_p,
 213                                   int **fc_p,
 214                                   int **indx_p,
 215                                   char **ptype_p){
 216
 217   /* make the DP arrays available to routines such as subopt() */
 218   *f5_p = f5; *c_p = c;
 219   *fML_p = fML; *fM1_p = fM1;
 220   *indx_p = indx; *ptype_p = ptype;
 221   *fc_p =fc;
 222 }
 223
 224 /*--------------------------------------------------------------------------*/
 225
 226 PUBLIC float cofold(const char *string, char *structure) {
 227   return cofold_par(string, structure, NULL, fold_constrained);
 228 }
 229
 230 PUBLIC float cofold_par(const char *string,
 231                         char *structure,
 232                         paramT *parameters,
 233                         int is_constrained){
 234
 235   int i, length, energy, bonus=0, bonus_cnt=0;
 236
 237   zuker = 0;
 238
 239   struct_constrained  = is_constrained;
 240   length              = (int) strlen(string);
 241
 242 #ifdef _OPENMP
 243   /* always init everything since all global static variables are uninitialized when entering a thread */
 244   init_cofold(length, parameters);
 245 #else
 246   if(parameters) init_cofold(length, parameters);
 247   else if (length>init_length) init_cofold(length, parameters);
 248   else if (fabs(P->temperature - temperature)>1e-6) update_cofold_params_par(parameters);
 249 #endif
 250
 251   with_gquad  = P->model_details.gquad;
 252   S           = encode_sequence(string, 0);
 253   S1          = encode_sequence(string, 1);
 254   S1[0]       = S[0]; /* store length at pos. 0 */
 255
 256   BP = (int *)space(sizeof(int)*(length+2));
 257   make_ptypes(S, structure);
 258
 259   energy = fill_arrays(string);
 260
 261   backtrack(string);
 262
 263 #ifdef PAREN
 264   parenthesis_structure(structure, base_pair2, length);
 265 #else
 266   letter_structure(structure, base_pair2, length);
 267 #endif
 268
 269   /*
 270   *  Backward compatibility:
 271   *  This block may be removed if deprecated functions
 272   *  relying on the global variable "base_pair" vanish from within the package!
 273   */
 274   base_pair = base_pair2;
 275   /*
 276   {
 277     if(base_pair) free(base_pair);
 278     base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
 279     memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
 280   }
 281   */
 282
 283   /* check constraints */
 284   for(i=1;i<=length;i++) {
 285     if((BP[i]<0)&&(BP[i]>-4)) {
 286       bonus_cnt++;
 287       if((BP[i]==-3)&&(structure[i-1]==')')) bonus++;
 288       if((BP[i]==-2)&&(structure[i-1]=='(')) bonus++;
 289       if((BP[i]==-1)&&(structure[i-1]!='.')) bonus++;
 290     }
 291
 292     if(BP[i]>i) {
 293       int l;
 294       bonus_cnt++;
 295       for(l=1; l<=base_pair2[0].i; l++)
 296         if(base_pair2[l].i != base_pair2[l].j)
 297           if((i==base_pair2[l].i)&&(BP[i]==base_pair2[l].j)) bonus++;
 298     }
 299   }
 300
 301   if (bonus_cnt>bonus) fprintf(stderr,"\ncould not enforce all constraints\n");
 302   bonus*=BONUS;
 303
 304   free(S); free(S1); free(BP);
 305
 306   energy += bonus;      /*remove bonus energies from result */
 307
 308   if (backtrack_type=='C')
 309     return (float) c[indx[length]+1]/100.;
 310   else if (backtrack_type=='M')
 311     return (float) fML[indx[length]+1]/100.;
 312   else
 313     return (float) energy/100.;
 314 }
 315
 316 PRIVATE int fill_arrays(const char *string) {
 317   /* fill "c", "fML" and "f5" arrays and return  optimal energy */
 318
 319   int   i, j, k, length, energy;
 320   int   decomp, new_fML, max_separation;
 321   int   no_close, type, type_2, tt, maxj;
 322   int   bonus=0;
 323   int   dangle_model  = P->model_details.dangles;
 324   int   noGUclosure   = P->model_details.noGUclosure;
 325   int   noLP          = P->model_details.noLP;
 326
 327   length = (int) strlen(string);
 328
 329   max_separation = (int) ((1.-LOCALITY)*(double)(length-2)); /* not in use */
 330
 331   if(with_gquad)
 332     ggg = get_gquad_matrix(S, P);
 333
 334   for (j=1; j<=length; j++) {
 335     Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
 336     fc[j]=0;
 337   }
 338
 339   for (j = 1; j<=length; j++)
 340     for (i=1; i<=j; i++) {
 341       c[indx[j]+i] = fML[indx[j]+i] = INF;
 342       if (uniq_ML) fM1[indx[j]+i] = INF;
 343     }
 344
 345   for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
 346
 347     maxj=(zuker)? (MIN2(i+cut_point-1,length)):length;
 348     for (j = i+TURN+1; j <= maxj; j++) {
 349       int p, q, ij;
 350       ij = indx[j]+i;
 351       bonus = 0;
 352       type = ptype[ij];
 353
 354       /* enforcing structure constraints */
 355       if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
 356       if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
 357       if ((BP[i]==-4)||(BP[j]==-4)) type=0;
 358
 359       no_close = (((type==3)||(type==4))&&noGUclosure&&(bonus==0));
 360
 361       if (j-i-1 > max_separation) type = 0;  /* forces locality degree */
 362
 363       if (type) {   /* we have a pair */
 364         int new_c=0, stackEnergy=INF;
 365         short si, sj;
 366         si  = SAME_STRAND(i, i+1) ? S1[i+1] : -1;
 367         sj  = SAME_STRAND(j-1, j) ? S1[j-1] : -1;
 368         /* hairpin ----------------------------------------------*/
 369
 370         if (SAME_STRAND(i,j)) {
 371           if (no_close) new_c = FORBIDDEN;
 372           else
 373             new_c = E_Hairpin(j-i-1, type, si, sj, string+i-1, P);
 374         }
 375         else {
 376           if (dangle_model)
 377             new_c += E_ExtLoop(rtype[type], sj, si, P);
 378           else
 379             new_c += E_ExtLoop(rtype[type], -1, -1, P);
 380         }
 381         /*--------------------------------------------------------
 382           check for elementary structures involving more than one
 383           closing pair.
 384           --------------------------------------------------------*/
 385
 386         for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
 387           int minq = j-i+p-MAXLOOP-2;
 388           if (minq<p+1+TURN) minq = p+1+TURN;
 389           for (q = minq; q < j; q++) {
 390             type_2 = ptype[indx[q]+p];
 391
 392             if (type_2==0) continue;
 393             type_2 = rtype[type_2];
 394
 395             if (noGUclosure)
 396               if (no_close||(type_2==3)||(type_2==4))
 397                 if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
 398
 399             if (SAME_STRAND(i,p) && SAME_STRAND(q,j))
 400               energy = E_IntLoop(p-i-1, j-q-1, type, type_2, si, sj, S1[p-1], S1[q+1], P);
 401             else
 402               energy = E_IntLoop_Co(rtype[type], rtype[type_2],
 403                                     i, j, p, q,
 404                                     cut_point,
 405                                     si, sj,
 406                                     S1[p-1], S1[q+1],
 407                                     dangle_model,
 408                                     P);
 409
 410             new_c = MIN2(energy+c[indx[q]+p], new_c);
 411             if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
 412
 413           } /* end q-loop */
 414         } /* end p-loop */
 415
 416         /* multi-loop decomposition ------------------------*/
 417
 418
 419         if (!no_close) {
 420           int MLenergy;
 421
 422           if((si >= 0) && (sj >= 0)){
 423             decomp = DMLi1[j-1];
 424             tt = rtype[type];
 425             MLenergy = P->MLclosing;
 426             switch(dangle_model){
 427               case 0:   MLenergy += decomp + E_MLstem(tt, -1, -1, P);
 428                         break;
 429               case 2:   MLenergy += decomp + E_MLstem(tt, sj, si, P);
 430                         break;
 431               default:  decomp += E_MLstem(tt, -1, -1, P);
 432                         decomp = MIN2(decomp, DMLi1[j-2] + E_MLstem(tt, sj, -1, P) + P->MLbase);
 433                         decomp = MIN2(decomp, DMLi2[j-1] + E_MLstem(tt, -1, si, P) + P->MLbase);
 434                         decomp = MIN2(decomp, DMLi2[j-2] + E_MLstem(tt, sj, si, P) + 2*P->MLbase);
 435                         MLenergy += decomp;
 436                         break;
 437             }
 438             new_c = MIN2(new_c, MLenergy);
 439           }
 440
 441           if (!SAME_STRAND(i,j)) { /* cut is somewhere in the multiloop*/
 442             decomp = fc[i+1] + fc[j-1];
 443             tt = rtype[type];
 444             switch(dangle_model){
 445               case 0:   decomp += E_ExtLoop(tt, -1, -1, P);
 446                         break;
 447               case 2:   decomp += E_ExtLoop(tt, sj, si, P);
 448                         break;
 449               default:  decomp += E_ExtLoop(tt, -1, -1, P);
 450                         decomp = MIN2(decomp, fc[i+2] + fc[j-2] + E_ExtLoop(tt, sj, si, P));
 451                         decomp = MIN2(decomp, fc[i+2] + fc[j-1] + E_ExtLoop(tt, -1, si, P));
 452                         decomp = MIN2(decomp, fc[i+1] + fc[j-2] + E_ExtLoop(tt, sj, -1, P));
 453                         break;
 454             }
 455             new_c = MIN2(new_c, decomp);
 456           }
 457         } /* end >> if (!no_close) << */
 458
 459         /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
 460
 461         if (dangle_model==3) {
 462           decomp = INF;
 463           for (k = i+2+TURN; k < j-2-TURN; k++) {
 464             type_2 = ptype[indx[k]+i+1]; type_2 = rtype[type_2];
 465             if (type_2)
 466               decomp = MIN2(decomp, c[indx[k]+i+1]+P->stack[type][type_2]+
 467                             fML[indx[j-1]+k+1]);
 468             type_2 = ptype[indx[j-1]+k+1]; type_2 = rtype[type_2];
 469             if (type_2)
 470               decomp = MIN2(decomp, c[indx[j-1]+k+1]+P->stack[type][type_2]+
 471                             fML[indx[k]+i+1]);
 472           }
 473           /* no TermAU penalty if coax stack */
 474           decomp += 2*P->MLintern[1] + P->MLclosing;
 475           new_c = MIN2(new_c, decomp);
 476         }
 477
 478         if(with_gquad){
 479           /* include all cases where a g-quadruplex may be enclosed by base pair (i,j) */
 480           if (!no_close && SAME_STRAND(i,j)) {
 481             tt = rtype[type];
 482             energy = E_GQuad_IntLoop(i, j, type, S1, ggg, indx, P);
 483             new_c = MIN2(new_c, energy);
 484           }
 485         }
 486
 487         new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
 488         cc[j] = new_c + bonus;
 489         if (noLP){
 490           if (SAME_STRAND(i,i+1) && SAME_STRAND(j-1,j))
 491             c[ij] = cc1[j-1]+stackEnergy+bonus;
 492           else /* currently we don't allow stacking over the cut point */
 493             c[ij] = FORBIDDEN;
 494         }
 495         else
 496           c[ij] = cc[j];
 497
 498       } /* end >> if (pair) << */
 499
 500       else c[ij] = INF;
 501
 502
 503       /* done with c[i,j], now compute fML[i,j] */
 504       /* free ends ? -----------------------------------------*/
 505       new_fML=INF;
 506       if (SAME_STRAND(i-1,i)) {
 507         if (SAME_STRAND(i,i+1)) new_fML = fML[ij+1]+P->MLbase;
 508         if (SAME_STRAND(j-1,j)) new_fML = MIN2(fML[indx[j-1]+i]+P->MLbase, new_fML);
 509         if (SAME_STRAND(j,j+1)) {
 510           energy = c[ij];
 511           if(dangle_model == 2) energy += E_MLstem(type,(i>1) ? S1[i-1] : -1, (j<length) ? S1[j+1] : -1, P);
 512           else energy += E_MLstem(type, -1, -1, P);
 513           new_fML = MIN2(new_fML, energy);
 514           if(uniq_ML){
 515             fM1[ij] = energy;
 516             if(SAME_STRAND(j-1,j)) fM1[ij] = MIN2(energy, fM1[indx[j-1]+i] + P->MLbase);
 517           }
 518         }
 519         if (dangle_model%2==1) {  /* normal dangles */
 520           if (SAME_STRAND(i,i+1)) {
 521             tt = ptype[ij+1]; /* i+1,j */
 522             new_fML = MIN2(new_fML, c[ij+1] + P->MLbase + E_MLstem(tt, S1[i], -1, P));
 523           }
 524           if (SAME_STRAND(j-1,j)) {
 525             tt = ptype[indx[j-1]+i]; /* i,j-1 */
 526             new_fML = MIN2(new_fML, c[indx[j-1]+i] + P->MLbase + E_MLstem(tt, -1, S1[j], P));
 527           }
 528           if ((SAME_STRAND(j-1,j))&&(SAME_STRAND(i,i+1))) {
 529             tt = ptype[indx[j-1]+i+1]; /* i+1,j-1 */
 530             new_fML = MIN2(new_fML, c[indx[j-1]+i+1] + 2*P->MLbase + E_MLstem(tt, S1[i], S1[j], P));
 531           }
 532         }
 533       }
 534
 535       if(with_gquad){
 536         if(SAME_STRAND(i, j))
 537           new_fML = MIN2(new_fML, ggg[indx[j] + i] + E_MLstem(0, -1, -1, P));
 538       }
 539
 540       /* modular decomposition -------------------------------*/
 541
 542       {
 543         int stopp;     /*loop 1 up to cut, then loop 2*/
 544         stopp=(cut_point>0)? (cut_point):(j-2-TURN);
 545         for (decomp=INF, k = i+1+TURN; k<stopp; k++)
 546           decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);
 547         k++;
 548         for (;k <= j-2-TURN;k++)
 549           decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);
 550       }
 551       DMLi[j] = decomp;               /* store for use in ML decompositon */
 552       new_fML = MIN2(new_fML,decomp);
 553
 554       /* coaxial stacking */
 555       if (dangle_model==3) {
 556         int stopp;
 557         stopp=(cut_point>0)? (cut_point):(j-2-TURN);
 558         /* additional ML decomposition as two coaxially stacked helices */
 559         for (decomp = INF, k = i+1+TURN; k<stopp; k++) {
 560           type = ptype[indx[k]+i]; type = rtype[type];
 561           type_2 = ptype[indx[j]+k+1]; type_2 = rtype[type_2];
 562           if (type && type_2)
 563             decomp = MIN2(decomp,
 564                           c[indx[k]+i]+c[indx[j]+k+1]+P->stack[type][type_2]);
 565         }
 566         k++;
 567         for (;k <= j-2-TURN; k++) {
 568           type = ptype[indx[k]+i]; type = rtype[type];
 569           type_2 = ptype[indx[j]+k+1]; type_2 = rtype[type_2];
 570           if (type && type_2)
 571             decomp = MIN2(decomp,
 572                           c[indx[k]+i]+c[indx[j]+k+1]+P->stack[type][type_2]);
 573         }
 574
 575         decomp += 2*P->MLintern[1];
 576
 577 #if 0
 578         /* This is needed for Y shaped ML loops with coax stacking of
 579            interior pairs, but backtracking will fail if activated */
 580         DMLi[j] = MIN2(DMLi[j], decomp);
 581         if (SAME_STRAND(j-1,j)) DMLi[j] = MIN2(DMLi[j], DMLi[j-1]+P->MLbase);
 582         if (SAME_STRAND(i,i+1)) DMLi[j] = MIN2(DMLi[j], DMLi1[j]+P->MLbase);
 583         new_fML = MIN2(new_fML, DMLi[j]);
 584 #endif
 585         new_fML = MIN2(new_fML, decomp);
 586       }
 587
 588       fML[ij] = Fmi[j] = new_fML;     /* substring energy */
 589
 590     }
 591
 592     if (i==cut_point)
 593       for (j=i; j<=maxj; j++)
 594         free_end(fc, j, cut_point);
 595     if (i<cut_point)
 596       free_end(fc,i,cut_point-1);
 597
 598
 599     {
 600       int *FF; /* rotate the auxilliary arrays */
 601       FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
 602       FF = cc1; cc1=cc; cc=FF;
 603       for (j=1; j<=maxj; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
 604     }
 605   }
 606
 607   /* calculate energies of 5' and 3' fragments */
 608
 609   for (i=1; i<=length; i++)
 610     free_end(f5, i, 1);
 611
 612   if (cut_point>0) {
 613     mfe1=f5[cut_point-1];
 614     mfe2=fc[length];
 615     /* add DuplexInit, check whether duplex*/
 616     for (i=cut_point; i<=length; i++) {
 617       f5[i]=MIN2(f5[i]+P->DuplexInit, fc[i]+fc[1]);
 618     }
 619   }
 620
 621   energy = f5[length];
 622   if (cut_point<1) mfe1=mfe2=energy;
 623   return energy;
 624 }
 625
 626 PRIVATE void backtrack_co(const char *string, int s, int b /* b=0: start new structure, b \ne 0: add to existing structure */) {
 627
 628   /*------------------------------------------------------------------
 629     trace back through the "c", "fc", "f5" and "fML" arrays to get the
 630     base pairing list. No search for equivalent structures is done.
 631     This is fast, since only few structure elements are recalculated.
 632     ------------------------------------------------------------------*/
 633
 634   int   i, j, k, length, energy, new;
 635   int   no_close, type, type_2, tt;
 636   int   bonus;
 637   int   dangle_model  = P->model_details.dangles;
 638   int   noGUclosure   = P->model_details.noGUclosure;
 639   int   noLP          = P->model_details.noLP;
 640
 641   /* int   b=0;*/
 642
 643   length = strlen(string);
 644   if (s==0) {
 645     sector[++s].i = 1;
 646     sector[s].j   = length;
 647     sector[s].ml  = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
 648   }
 649   while (s>0) {
 650     int ml, fij, fi, cij, traced, i1, j1, mm, p, q, jj=0, gq=0;
 651     int canonical = 1;     /* (i,j) closes a canonical structure */
 652     i  = sector[s].i;
 653     j  = sector[s].j;
 654     ml = sector[s--].ml;   /* ml is a flag indicating if backtracking is to
 655                               occur in the fML- (1) or in the f-array (0) */
 656     if (ml==2) {
 657       base_pair2[++b].i = i;
 658       base_pair2[b].j   = j;
 659       goto repeat1;
 660     }
 661
 662     if (j < i+TURN+1) continue; /* no more pairs in this interval */
 663
 664
 665     if (ml==0) {fij = f5[j]; fi = f5[j-1];}
 666     else if (ml==1) {fij = fML[indx[j]+i]; fi = fML[indx[j-1]+i]+P->MLbase;}
 667     else /* 3 or 4 */ {
 668       fij = fc[j];
 669       fi = (ml==3) ? INF : fc[j-1];
 670     }
 671     if (fij == fi) {  /* 3' end is unpaired */
 672       sector[++s].i = i;
 673       sector[s].j   = j-1;
 674       sector[s].ml  = ml;
 675       continue;
 676     }
 677
 678     if (ml==0 || ml==4) { /* backtrack in f5 or fc[i=cut,j>cut] */
 679       int *ff;
 680       ff = (ml==4) ? fc : f5;
 681       switch(dangle_model){
 682         case 0:   /* j or j-1 is paired. Find pairing partner */
 683                   for (k=j-TURN-1,traced=0; k>=i; k--) {
 684                     int cc;
 685
 686                     if(with_gquad){
 687                       if(fij == ff[k-1] + ggg[indx[j]+k]){
 688                         /* found the decomposition */
 689                         traced = j; jj = k - 1; gq = 1;
 690                         break;
 691                       }
 692                     }
 693
 694                     type = ptype[indx[j]+k];
 695                     if(type){
 696                       cc = c[indx[j]+k];
 697                       if(!SAME_STRAND(k,j)) cc += P->DuplexInit;
 698                       if(fij == ff[k-1] + cc + E_ExtLoop(type, -1, -1, P)){
 699                         traced = j; jj = k-1;
 700                       }
 701                     }
 702                     if(traced) break;
 703                   }
 704
 705                   break;
 706
 707         case 2:   /* j or j-1 is paired. Find pairing partner */
 708                   for (k=j-TURN-1,traced=0; k>=i; k--) {
 709                     int cc;
 710
 711                     if(with_gquad){
 712                       if(fij == ff[k-1] + ggg[indx[j]+k]){
 713                         /* found the decomposition */
 714                         traced = j; jj = k - 1; gq = 1;
 715                         break;
 716                       }
 717                     }
 718
 719                     type = ptype[indx[j]+k];
 720                     if(type){
 721                       cc = c[indx[j]+k];
 722                       if(!SAME_STRAND(k,j)) cc += P->DuplexInit;
 723                       if(fij == ff[k-1] + cc + E_ExtLoop(type, (k>1) && SAME_STRAND(k-1,k) ? S1[k-1] : -1, (j<length) && SAME_STRAND(j,j+1) ? S1[j+1] : -1, P)){
 724                         traced = j; jj = k-1;
 725                       }
 726                     }
 727                     if(traced) break;
 728                   }
 729                   break;
 730
 731         default:  for(k=j-TURN-1,traced=0; k>=i; k--){
 732                     int cc;
 733                     type = ptype[indx[j]+k];
 734
 735                     if(with_gquad){
 736                       if(fij == ff[k-1] + ggg[indx[j]+k]){
 737                         /* found the decomposition */
 738                         traced = j; jj = k - 1; gq = 1;
 739                         break;
 740                       }
 741                     }
 742
 743                     if(type){
 744                       cc = c[indx[j]+k];
 745                       if(!SAME_STRAND(k,j)) cc += P->DuplexInit;
 746                       if(fij == ff[k-1] + cc + E_ExtLoop(type, -1, -1, P)){
 747                         traced = j; jj = k-1; break;
 748                       }
 749                       if((k>1) && SAME_STRAND(k-1,k))
 750                         if(fij == ff[k-2] + cc + E_ExtLoop(type, S1[k-1], -1, P)){
 751                               traced=j; jj=k-2; break;
 752                         }
 753                     }
 754
 755                     type = ptype[indx[j-1]+k];
 756                     if(type && SAME_STRAND(j-1,j)){
 757                       cc = c[indx[j-1]+k];
 758                       if (!SAME_STRAND(k,j-1)) cc += P->DuplexInit; /*???*/
 759                       if (fij == cc + ff[k-1] + E_ExtLoop(type, -1, S1[j], P)){
 760                             traced=j-1; jj = k-1; break;
 761                       }
 762                       if(k>i){
 763                         if (fij == ff[k-2] + cc + E_ExtLoop(type, SAME_STRAND(k-1,k) ? S1[k-1] : -1, S1[j], P)){
 764                           traced=j-1; jj=k-2; break;
 765                         }
 766                       }
 767                     }
 768                   }
 769
 770                   break;
 771       }
 772       if (!traced) nrerror("backtrack failed in f5 (or fc)");
 773       sector[++s].i = i;
 774       sector[s].j   = jj;
 775       sector[s].ml  = ml;
 776
 777       i=k; j=traced;
 778
 779       if(with_gquad && gq){
 780         /* goto backtrace of gquadruplex */
 781         goto repeat_gquad;
 782       }
 783
 784       base_pair2[++b].i = i;
 785       base_pair2[b].j   = j;
 786       goto repeat1;
 787     }
 788     else if (ml==3) { /* backtrack in fc[i<cut,j=cut-1] */
 789       if (fc[i] == fc[i+1]) { /* 5' end is unpaired */
 790         sector[++s].i = i+1;
 791         sector[s].j   = j;
 792         sector[s].ml  = ml;
 793         continue;
 794       }
 795       /* i or i+1 is paired. Find pairing partner */
 796       switch(dangle_model){
 797         case 0:   for (k=i+TURN+1, traced=0; k<=j; k++){
 798                     jj=k+1;
 799                     type = ptype[indx[k]+i];
 800                     if (type) {
 801                       if(fc[i] == fc[k+1] + c[indx[k]+i] + E_ExtLoop(type, -1, -1, P)){
 802                         traced = i;
 803                       }
 804                     }
 805                     if (traced) break;
 806                   }
 807                   break;
 808         case 2:   for (k=i+TURN+1, traced=0; k<=j; k++){
 809                     jj=k+1;
 810                     type = ptype[indx[k]+i];
 811                     if(type){
 812                       if(fc[i] == fc[k+1] + c[indx[k]+i] + E_ExtLoop(type,(i>1 && SAME_STRAND(i-1,i)) ? S1[i-1] : -1,  SAME_STRAND(k,k+1) ? S1[k+1] : -1, P)){
 813                         traced = i;
 814                       }
 815                     }
 816                     if (traced) break;
 817                   }
 818                   break;
 819         default:  for(k=i+TURN+1, traced=0; k<=j; k++){
 820                     jj=k+1;
 821                     type = ptype[indx[k]+i];
 822                     if(type){
 823                       if(fc[i] == fc[k+1] + c[indx[k]+i] + E_ExtLoop(type, -1, -1, P)){
 824                         traced = i; break;
 825                       }
 826                       else if(fc[i] == fc[k+2] + c[indx[k]+i] + E_ExtLoop(type, -1, SAME_STRAND(k,k+1) ? S1[k+1] : -1, P)){
 827                         traced = i; jj=k+2; break;
 828                       }
 829                     }
 830                     type = ptype[indx[k]+i+1];
 831                     if(type){
 832                       if(fc[i] == fc[k+1] + c[indx[k]+i+1] + E_ExtLoop(type, SAME_STRAND(i, i+1) ? S1[i] : -1, -1, P)){
 833                         traced = i+1; break;
 834                       }
 835                       if(k<j){
 836                         if(fc[i] == fc[k+2] + c[indx[k]+i+1] + E_ExtLoop(type, SAME_STRAND(i, i+1) ? S1[i] : -1, SAME_STRAND(k, k+1) ? S1[k+1] : -1, P)){
 837                           traced = i+1; jj=k+2; break;
 838                         }
 839                       }
 840                     }
 841                   }
 842                   break;
 843       }
 844
 845       if (!traced) nrerror("backtrack failed in fc[] 5' of cut");
 846
 847       sector[++s].i = jj;
 848       sector[s].j   = j;
 849       sector[s].ml  = ml;
 850
 851       j=k; i=traced;
 852       base_pair2[++b].i = i;
 853       base_pair2[b].j   = j;
 854       goto repeat1;
 855     }
 856
 857     else { /* true multi-loop backtrack in fML */
 858       if (fML[indx[j]+i+1]+P->MLbase == fij) { /* 5' end is unpaired */
 859         sector[++s].i = i+1;
 860         sector[s].j   = j;
 861         sector[s].ml  = ml;
 862         continue;
 863       }
 864
 865       if(with_gquad){
 866         if(fij == ggg[indx[j]+i] + E_MLstem(0, -1, -1, P)){
 867           /* go to backtracing of quadruplex */
 868           goto repeat_gquad;
 869         }
 870       }
 871
 872       tt  = ptype[indx[j]+i];
 873       cij = c[indx[j]+i];
 874       switch(dangle_model){
 875         case 0:   if(fij == cij + E_MLstem(tt, -1, -1, P)){
 876                     base_pair2[++b].i  = i;
 877                     base_pair2[b].j    = j;
 878                     goto repeat1;
 879                   }
 880                   break;
 881         case 2:   if(fij == cij + E_MLstem(tt, (i>1) ? S1[i-1] : -1, (j<length) ? S1[j+1] : -1, P)){
 882                     base_pair2[++b].i  = i;
 883                     base_pair2[b].j    = j;
 884                     goto repeat1;
 885                   }
 886                   break;
 887         default:  if(fij == cij + E_MLstem(tt, -1, -1, P)){
 888                     base_pair2[++b].i  = i;
 889                     base_pair2[b].j    = j;
 890                     goto repeat1;
 891                   }
 892                   tt = ptype[indx[j]+i+1];
 893                   if(fij == c[indx[j]+i+1] + P->MLbase + E_MLstem(tt, S1[i], -1, P)){
 894                     i++;
 895                     base_pair2[++b].i  = i;
 896                     base_pair2[b].j    = j;
 897                     goto repeat1;
 898                   }
 899                   tt = ptype[indx[j-1]+i];
 900                   if(fij == c[indx[j-1]+i] + P->MLbase + E_MLstem(tt, -1, S1[j], P)){
 901                     j--;
 902                     base_pair2[++b].i  = i;
 903                     base_pair2[b].j    = j;
 904                     goto repeat1;
 905                   }
 906                   tt = ptype[indx[j-1]+i+1];
 907                   if(fij == c[indx[j-1]+i+1] + 2*P->MLbase + E_MLstem(tt, S1[i], S1[j], P)){
 908                     i++; j--;
 909                     base_pair2[++b].i  = i;
 910                     base_pair2[b].j    = j;
 911                     goto repeat1;
 912                   }
 913                   break;
 914       }
 915
 916       /* find next component of multiloop */
 917       for (k = i+1+TURN; k <= j-2-TURN; k++)
 918         if (fij == (fML[indx[k]+i]+fML[indx[j]+k+1]))
 919           break;
 920
 921       if ((dangle_model==3)&&(k>j-2-TURN)) { /* must be coax stack */
 922         ml = 2;
 923         for (k = i+1+TURN; k <= j-2-TURN; k++) {
 924           type = ptype[indx[k]+i];  type= rtype[type];
 925           type_2 = ptype[indx[j]+k+1]; type_2= rtype[type_2];
 926           if (type && type_2)
 927             if (fij == c[indx[k]+i]+c[indx[j]+k+1]+P->stack[type][type_2]+
 928                 2*P->MLintern[1])
 929               break;
 930         }
 931       }
 932
 933       sector[++s].i = i;
 934       sector[s].j   = k;
 935       sector[s].ml  = ml;
 936       sector[++s].i = k+1;
 937       sector[s].j   = j;
 938       sector[s].ml  = ml;
 939
 940       if (k>j-2-TURN) nrerror("backtrack failed in fML");
 941       continue;
 942     }
 943
 944   repeat1:
 945
 946     /*----- begin of "repeat:" -----*/
 947     if (canonical)  cij = c[indx[j]+i];
 948     type = ptype[indx[j]+i];
 949
 950     bonus = 0;
 951
 952     if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
 953     if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
 954
 955     if (noLP)
 956       if (cij == c[indx[j]+i]) {
 957         /* (i.j) closes canonical structures, thus
 958            (i+1.j-1) must be a pair                */
 959         type_2 = ptype[indx[j-1]+i+1]; type_2 = rtype[type_2];
 960         cij -= P->stack[type][type_2] + bonus;
 961         base_pair2[++b].i = i+1;
 962         base_pair2[b].j   = j-1;
 963         i++; j--;
 964         canonical=0;
 965         goto repeat1;
 966       }
 967     canonical = 1;
 968
 969
 970     no_close = (((type==3)||(type==4))&&noGUclosure&&(bonus==0));
 971     if (SAME_STRAND(i,j)) {
 972       if (no_close) {
 973         if (cij == FORBIDDEN) continue;
 974       } else
 975         if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1, P)+bonus)
 976           continue;
 977     }
 978     else {
 979       if(dangle_model){
 980         if(cij == E_ExtLoop(rtype[type], SAME_STRAND(j-1,j) ? S1[j-1] : -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P)) continue;
 981       }
 982       else if(cij == E_ExtLoop(rtype[type], -1, -1, P)) continue;
 983     }
 984
 985     for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
 986       int minq;
 987       minq = j-i+p-MAXLOOP-2;
 988       if (minq<p+1+TURN) minq = p+1+TURN;
 989       for (q = j-1; q >= minq; q--) {
 990
 991         type_2 = ptype[indx[q]+p];
 992         if (type_2==0) continue;
 993         type_2 = rtype[type_2];
 994         if (noGUclosure)
 995           if (no_close||(type_2==3)||(type_2==4))
 996             if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
 997
 998         /* energy = oldLoopEnergy(i, j, p, q, type, type_2); */
 999         if (SAME_STRAND(i,p) && SAME_STRAND(q,j))
1000           energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
1001                               S1[i+1], S1[j-1], S1[p-1], S1[q+1], P);
1002         else {
1003           energy = E_IntLoop_Co(rtype[type], rtype[type_2], i, j, p, q, cut_point, S1[i+1], S1[j-1], S1[p-1], S1[q+1], dangle_model, P);
1004         }
1005
1006         new = energy+c[indx[q]+p]+bonus;
1007         traced = (cij == new);
1008         if (traced) {
1009           base_pair2[++b].i = p;
1010           base_pair2[b].j   = q;
1011           i = p, j = q;
1012           goto repeat1;
1013         }
1014       }
1015     }
1016
1017     /* end of repeat: --------------------------------------------------*/
1018
1019     /* (i.j) must close a fake or true multi-loop */
1020     tt = rtype[type];
1021     i1 = i+1;
1022     j1 = j-1;
1023
1024     if(with_gquad){
1025       /*
1026         The case that is handled here actually resembles something like
1027         an interior loop where the enclosing base pair is of regular
1028         kind and the enclosed pair is not a canonical one but a g-quadruplex
1029         that should then be decomposed further...
1030       */
1031       if(SAME_STRAND(i,j)){
1032         if(backtrack_GQuad_IntLoop(cij, i, j, type, S, ggg, indx, &p, &q, P)){
1033           i = p; j = q;
1034           goto repeat_gquad;
1035         }
1036       }
1037     }
1038
1039     /* fake multi-loop */
1040     if(!SAME_STRAND(i,j)){
1041       int ii, jj, decomp;
1042       ii = jj = 0;
1043       decomp = fc[i1] + fc[j1];
1044       switch(dangle_model){
1045         case 0:   if(cij == decomp + E_ExtLoop(tt, -1, -1, P)){
1046                     ii=i1, jj=j1;
1047                   }
1048                   break;
1049         case 2:   if(cij == decomp + E_ExtLoop(tt, SAME_STRAND(j-1,j) ? S1[j-1] : -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P)){
1050                     ii=i1, jj=j1;
1051                   }
1052
1053                   break;
1054         default:  if(cij == decomp + E_ExtLoop(tt, -1, -1, P)){
1055                     ii=i1, jj=j1;
1056                   }
1057                   else if(cij == fc[i+2] + fc[j-1] + E_ExtLoop(tt, -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P)){
1058                     ii = i+2; jj = j1;
1059                   }
1060                   else if(cij == fc[i+1] + fc[j-2] + E_ExtLoop(tt, SAME_STRAND(j-1,j) ? S1[j-1] : -1, -1, P)){
1061                     ii = i1; jj = j-2;
1062                   }
1063                   else if(cij == fc[i+2] + fc[j-2] + E_ExtLoop(tt, SAME_STRAND(j-1,j) ? S1[j-1] : -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P)){
1064                     ii = i+2; jj = j-2;
1065                   }
1066                   break;
1067       }
1068       if(ii){
1069         sector[++s].i = ii;
1070         sector[s].j   = cut_point-1;
1071         sector[s].ml  = 3;
1072         sector[++s].i = cut_point;
1073         sector[s].j   = jj;
1074         sector[s].ml  = 4;
1075         continue;
1076       }
1077     }
1078
1079     /* true multi-loop */
1080     mm = bonus + P->MLclosing;
1081     sector[s+1].ml  = sector[s+2].ml = 1;
1082     int ml0   = E_MLstem(tt, -1, -1, P);
1083     int ml5   = E_MLstem(tt, SAME_STRAND(j-1,j) ? S1[j-1] : -1, -1, P);
1084     int ml3   = E_MLstem(tt, -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P);
1085     int ml53  = E_MLstem(tt, SAME_STRAND(j-1,j) ? S1[j-1] : -1, SAME_STRAND(i,i+1) ? S1[i+1] : -1, P);
1086     for (traced = 0, k = i+2+TURN; k < j-2-TURN; k++) {
1087       switch(dangle_model){
1088         case 0:   /* no dangles */
1089                   if(cij == mm + fML[indx[k]+i+1] + fML[indx[j-1]+k+1] + ml0)
1090                     traced = i+1;
1091                   break;
1092         case 2:   /*double dangles */
1093                   if(cij == mm + fML[indx[k]+i+1] + fML[indx[j-1]+k+1] + ml53)
1094                     traced = i+1;
1095                   break;
1096         default:  /* normal dangles */
1097                   if(cij == mm + fML[indx[k]+i+1] + fML[indx[j-1]+k+1] + ml0){
1098                     traced = i+1;
1099                     break;
1100                   }
1101                   else if (cij == fML[indx[k]+i+2] + fML[indx[j-1]+k+1] + ml3 + mm + P->MLbase){
1102                     traced = i1 = i+2;
1103                     break;
1104                   }
1105                   else if (cij == fML[indx[k]+i+1] + fML[indx[j-2]+k+1] + ml5 + mm + P->MLbase){
1106                     traced = i1 = i+1; j1 = j-2;
1107                     break;
1108                   }
1109                   else if (cij == fML[indx[k]+i+2] + fML[indx[j-2]+k+1] + ml53 + mm + 2*P->MLbase){
1110                     traced = i1 = i+2; j1 = j-2;
1111                     break;
1112                   }
1113                   break;
1114       }
1115       if(traced) break;
1116       /* coaxial stacking of (i.j) with (i+1.k) or (k.j-1) */
1117       /* use MLintern[1] since coax stacked pairs don't get TerminalAU */
1118       if (dangle_model==3) {
1119         int en;
1120         type_2 = ptype[indx[k]+i+1]; type_2 = rtype[type_2];
1121         if (type_2) {
1122           en = c[indx[k]+i+1]+P->stack[type][type_2]+fML[indx[j-1]+k+1];
1123           if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1124             ml = 2;
1125             sector[s+1].ml  = 2;
1126             break;
1127           }
1128         }
1129         type_2 = ptype[indx[j-1]+k+1]; type_2 = rtype[type_2];
1130         if (type_2) {
1131           en = c[indx[j-1]+k+1]+P->stack[type][type_2]+fML[indx[k]+i+1];
1132           if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1133             sector[s+2].ml = 2;
1134             break;
1135           }
1136         }
1137       }
1138     }
1139     if (k<=j-3-TURN) { /* found the decomposition */
1140       sector[++s].i = i1;
1141       sector[s].j   = k;
1142       sector[++s].i = k+1;
1143       sector[s].j   = j1;
1144     } else {
1145 #if 0
1146       /* Y shaped ML loops don't work yet */
1147       if (dangle_model==3) {
1148         /* (i,j) must close a Y shaped ML loop with coax stacking */
1149         if (cij == fML[indx[j-2]+i+2] + mm + d3 + d5 + P->MLbase + P->MLbase) {
1150           i1 = i+2;
1151           j1 = j-2;
1152         } else if (cij ==  fML[indx[j-2]+i+1] + mm + d5 + P->MLbase)
1153           j1 = j-2;
1154         else if (cij ==  fML[indx[j-1]+i+2] + mm + d3 + P->MLbase)
1155           i1 = i+2;
1156         else /* last chance */
1157           if (cij != fML[indx[j-1]+i+1] + mm + P->MLbase)
1158             fprintf(stderr,  "backtracking failed in repeat");
1159         /* if we arrive here we can express cij via fML[i1,j1]+dangles */
1160         sector[++s].i = i1;
1161         sector[s].j   = j1;
1162       }
1163       else
1164 #endif
1165         nrerror("backtracking failed in repeat");
1166     }
1167
1168     continue; /* this is a workarround to not accidentally proceed in the following block */
1169
1170   repeat_gquad:
1171     /*
1172       now we do some fancy stuff to backtrace the stacksize and linker lengths
1173       of the g-quadruplex that should reside within position i,j
1174     */
1175     {
1176       int l[3], L, a;
1177       L = -1;
1178
1179       get_gquad_pattern_mfe(S, i, j, P, &L, l);
1180       if(L != -1){
1181         /* fill the G's of the quadruplex into base_pair2 */
1182         for(a=0;a<L;a++){
1183           base_pair2[++b].i = i+a;
1184           base_pair2[b].j   = i+a;
1185           base_pair2[++b].i = i+L+l[0]+a;
1186           base_pair2[b].j   = i+L+l[0]+a;
1187           base_pair2[++b].i = i+L+l[0]+L+l[1]+a;
1188           base_pair2[b].j   = i+L+l[0]+L+l[1]+a;
1189           base_pair2[++b].i = i+L+l[0]+L+l[1]+L+l[2]+a;
1190           base_pair2[b].j   = i+L+l[0]+L+l[1]+L+l[2]+a;
1191         }
1192         goto repeat_gquad_exit;
1193       }
1194       nrerror("backtracking failed in repeat_gquad");
1195     }
1196   repeat_gquad_exit:
1197     asm("nop");
1198
1199   } /* end >> while (s>0) << */
1200
1201   base_pair2[0].i = b;    /* save the total number of base pairs */
1202 }
1203
1204 PRIVATE void free_end(int *array, int i, int start) {
1205   int inc, type, energy, length, j, left, right;
1206   int dangle_model = P->model_details.dangles;
1207
1208   inc = (i>start)? 1:-1;
1209   length = S[0];
1210
1211   if (i==start) array[i]=0;
1212   else array[i] = array[i-inc];
1213   if (inc>0) {
1214     left = start; right=i;
1215   } else {
1216     left = i; right = start;
1217   }
1218
1219   for (j=start; inc*(i-j)>TURN; j+=inc) {
1220     int ii, jj;
1221     short si, sj;
1222     if (i>j) { ii = j; jj = i;} /* inc>0 */
1223     else     { ii = i; jj = j;} /* inc<0 */
1224     type = ptype[indx[jj]+ii];
1225     if (type) {  /* i is paired with j */
1226       si = (ii>1)       && SAME_STRAND(ii-1,ii) ? S1[ii-1] : -1;
1227       sj = (jj<length)  && SAME_STRAND(jj,jj+1) ? S1[jj+1] : -1;
1228       energy = c[indx[jj]+ii];
1229       switch(dangle_model){
1230         case 0:
1231                   array[i] = MIN2(array[i], array[j-inc] + energy + E_ExtLoop(type, -1, -1, P));
1232                   break;
1233         case 2:
1234                   array[i] = MIN2(array[i], array[j-inc] + energy + E_ExtLoop(type, si, sj, P));
1235                   break;
1236         default:
1237                   array[i] = MIN2(array[i], array[j-inc] + energy + E_ExtLoop(type, -1, -1, P));
1238                   if(inc > 0){
1239                     if(j > left)
1240                     array[i] = MIN2(array[i], array[j-2] + energy + E_ExtLoop(type, si, -1, P));
1241                   }
1242                   else if(j < right)
1243                     array[i] = MIN2(array[i], array[j+2] + energy + E_ExtLoop(type, -1, sj, P));
1244                   break;
1245       }
1246     }
1247
1248     if(with_gquad){
1249       if(SAME_STRAND(ii, jj))
1250         array[i] = MIN2(array[i], array[j-inc] + ggg[indx[jj]+ii]);
1251     }
1252
1253     if (dangle_model%2==1) {
1254       /* interval ends in a dangle (i.e. i-inc is paired) */
1255       if (i>j) { ii = j; jj = i-1;} /* inc>0 */
1256       else     { ii = i+1; jj = j;} /* inc<0 */
1257       type = ptype[indx[jj]+ii];
1258       if (!type) continue;
1259
1260       si = (ii > left)  && SAME_STRAND(ii-1,ii) ? S1[ii-1] : -1;
1261       sj = (jj < right) && SAME_STRAND(jj,jj+1) ? S1[jj+1] : -1;
1262       energy = c[indx[jj]+ii];
1263       if(inc>0)
1264         array[i] = MIN2(array[i], array[j - inc] + energy + E_ExtLoop(type, -1, sj, P));
1265       else
1266         array[i] = MIN2(array[i], array[j - inc] + energy + E_ExtLoop(type, si, -1, P));
1267       if(j!= start){ /* dangle_model on both sides */
1268         array[i] = MIN2(array[i], array[j-2*inc] + energy + E_ExtLoop(type, si, sj, P));
1269       }
1270     }
1271   }
1272 }
1273
1274 PUBLIC void update_cofold_params(void){
1275   update_cofold_params_par(NULL);
1276 }
1277
1278 PUBLIC void update_cofold_params_par(paramT *parameters){
1279   if(P) free(P);
1280
1281   if(parameters){
1282     P = get_parameter_copy(parameters);
1283   } else {
1284     model_detailsT md;
1285     set_model_details(&md);
1286     P = get_scaled_parameters(temperature, md);
1287   }
1288   make_pair_matrix();
1289   if (init_length < 0) init_length=0;
1290 }
1291
1292 /*---------------------------------------------------------------------------*/
1293
1294 PRIVATE void make_ptypes(const short *S, const char *structure) {
1295   int n,i,j,k,l;
1296   int noLP = P->model_details.noLP;
1297
1298   n=S[0];
1299   for (k=1; k<n-TURN; k++)
1300     for (l=1; l<=2; l++) {
1301       int type,ntype=0,otype=0;
1302       i=k; j = i+TURN+l; if (j>n) continue;
1303       type = pair[S[i]][S[j]];
1304       while ((i>=1)&&(j<=n)) {
1305         if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
1306         if (noLP && (!otype) && (!ntype))
1307           type = 0; /* i.j can only form isolated pairs */
1308         ptype[indx[j]+i] = (char) type;
1309         otype =  type;
1310         type  = ntype;
1311         i--; j++;
1312       }
1313     }
1314
1315   if (struct_constrained && (structure != NULL))
1316     constrain_ptypes(structure, (unsigned int)n, ptype, BP, TURN, 0);
1317 }
1318
1319 PUBLIC void get_monomere_mfes(float *e1, float *e2) {
1320   /*exports monomere free energies*/
1321   *e1 = mfe1;
1322   *e2 = mfe2;
1323 }
1324
1325 PRIVATE void backtrack(const char *sequence) {
1326   /*routine to call backtrack_co from 1 to n, backtrack type??*/
1327   backtrack_co(sequence, 0,0);
1328 }
1329
1330 PRIVATE int comp_pair(const void *A, const void *B) {
1331   bondT *x,*y;
1332   int ex, ey;
1333   x = (bondT *) A;
1334   y = (bondT *) B;
1335   ex = c[indx[x->j]+x->i]+c[indx[x->i+length]+x->j];
1336   ey = c[indx[y->j]+y->i]+c[indx[y->i+length]+y->j];
1337   if (ex>ey) return 1;
1338   if (ex<ey) return -1;
1339   return (indx[x->j]+x->i - indx[y->j]+y->i);
1340 }
1341
1342 PUBLIC SOLUTION *zukersubopt(const char *string) {
1343   return zukersubopt_par(string, NULL);
1344 }
1345
1346 PUBLIC SOLUTION *zukersubopt_par(const char *string, paramT *parameters){
1347
1348 /* Compute zuker suboptimal. Here, we're abusing the cofold() code
1349    "double" sequence, compute dimerarray entries, track back every base pair.
1350    This is slightly wasteful compared to the normal solution */
1351
1352   char      *doubleseq, *structure, *mfestructure, **todo;
1353   int       i, j, counter, num_pairs, psize, p;
1354   float     energy;
1355   SOLUTION  *zukresults;
1356   bondT     *pairlist;
1357
1358   num_pairs       = counter = 0;
1359   zuker           = 1;
1360   length          = (int)strlen(string);
1361   doubleseq       = (char *)space((2*length+1)*sizeof(char));
1362   mfestructure    = (char *) space((unsigned) 2*length+1);
1363   structure       = (char *) space((unsigned) 2*length+1);
1364   zukresults      = (SOLUTION *)space(((length*(length-1))/2)*sizeof(SOLUTION));
1365   mfestructure[0] = '\0';
1366   BP              = (int *)space(sizeof(int)*(2*length+2));
1367
1368   /* double the sequence */
1369   strcpy(doubleseq,string);
1370   strcat(doubleseq,string);
1371   cut_point = length + 1;
1372
1373   /* get mfe and do forward recursion */
1374 #ifdef _OPENMP
1375   /* always init everything since all global static variables are uninitialized when entering a thread */
1376   init_cofold(2 * length, parameters);
1377 #else
1378   if(parameters) init_cofold(2 * length, parameters);
1379   else if ((2 * length) > init_length) init_cofold(2 * length, parameters);
1380   else if (fabs(P->temperature - temperature)>1e-6) update_cofold_params_par(parameters);
1381 #endif
1382
1383
1384   S     = encode_sequence(doubleseq, 0);
1385   S1    = encode_sequence(doubleseq, 1);
1386   S1[0] = S[0]; /* store length at pos. 0 */
1387   make_ptypes(S, NULL); /* no constraint folding possible (yet?) with zukersubopt */
1388
1389   (void)fill_arrays(doubleseq);
1390
1391   psize     = length;
1392   pairlist  = (bondT *) space(sizeof(bondT)*(psize+1));
1393   todo      = (char **) space(sizeof(char *)*(length+1));
1394   for (i=1; i<length; i++) {
1395     todo[i] = (char *) space(sizeof(char)*(length+1));
1396   }
1397
1398   /* Make a list of all base pairs */
1399   for (i=1; i<length; i++) {
1400     for (j=i+TURN2+1/*??*/; j<=length; j++) {
1401       if (ptype[indx[j]+i]==0) continue;
1402       if (num_pairs>=psize) {
1403         psize = 1.2*psize + 32;
1404         pairlist = xrealloc(pairlist, sizeof(bondT)*(psize+1));
1405       }
1406       pairlist[num_pairs].i   = i;
1407       pairlist[num_pairs++].j = j;
1408       todo[i][j]=1;
1409     }
1410   }
1411   qsort(pairlist, num_pairs, sizeof(bondT), comp_pair);
1412
1413   for (p=0; p<num_pairs; p++) {
1414     i=pairlist[p].i;
1415     j=pairlist[p].j;
1416     if (todo[i][j]) {
1417       int k;
1418       sector[1].i   = i;
1419       sector[1].j   = j;
1420       sector[1].ml  = 2;
1421       backtrack_co(doubleseq, 1,0);
1422       sector[1].i   = j;
1423       sector[1].j   = i + length;
1424       sector[1].ml  = 2;
1425       backtrack_co(doubleseq, 1,base_pair2[0].i);
1426       energy = c[indx[j]+i]+c[indx[i+length]+j];
1427       parenthesis_zuker(structure, base_pair2, length);
1428       zukresults[counter].energy      = energy;
1429       zukresults[counter++].structure = strdup(structure);
1430       for (k = 1; k <= base_pair2[0].i; k++) { /* mark all pairs in structure as done */
1431         int x,y;
1432         x=base_pair2[k].i;
1433         y=base_pair2[k].j;
1434         if (x>length) x-=length;
1435         if (y>length) y-=length;
1436         if (x>y) {
1437           int temp;
1438           temp=x; x=y; y=temp;
1439         }
1440         todo[x][y] = 0;
1441       }
1442     }
1443   }
1444   /*free zeugs*/
1445   free(pairlist);
1446   for (i=1; i<length; i++)
1447     free(todo[i]);
1448   free(todo);
1449   free(structure);
1450   free(mfestructure);
1451   free(doubleseq);
1452   zuker=0;
1453   free(S); free(S1); free(BP);
1454   return zukresults;
1455 }
1456
1457
1458 /*###########################################*/
1459 /*# deprecated functions below              #*/
1460 /*###########################################*/
1461
1462 PUBLIC void initialize_cofold(int length){ /* DO NOTHING */ }