binaries/src/ViennaRNA/lib/fold.c

   1 /** \file **/
   2
   3 /*
   4                   minimum free energy
   5                   RNA secondary structure prediction
   6
   7                   c Ivo Hofacker, Chrisoph Flamm
   8                   original implementation by
   9                   Walter Fontana
  10                   g-quadruplex support and threadsafety
  11                   by Ronny Lorenz
  12
  13                   Vienna RNA package
  14 */
  15
  16 #include <config.h>
  17 #include <stdio.h>
  18 #include <stdlib.h>
  19 #include <math.h>
  20 #include <ctype.h>
  21 #include <string.h>
  22 #include <limits.h>
  23
  24 #include "utils.h"
  25 #include "energy_par.h"
  26 #include "fold_vars.h"
  27 #include "pair_mat.h"
  28 #include "params.h"
  29 #include "loop_energies.h"
  30 #include "data_structures.h"
  31 #include "gquad.h"
  32 #include "fold.h"
  33
  34 #ifdef _OPENMP
  35 #include <omp.h>
  36 #endif
  37
  38
  39 #define PAREN
  40 #define STACK_BULGE1      1       /* stacking energies for bulges of size 1 */
  41 #define NEW_NINIO         1       /* new asymetry penalty */
  42 #define MAXSECTORS        500     /* dimension for a backtrack array */
  43 #define LOCALITY          0.      /* locality parameter for base-pairs */
  44
  45 #define SAME_STRAND(I,J)  (((I)>=cut_point)||((J)<cut_point))
  46
  47 /*
  48 #################################
  49 # GLOBAL VARIABLES              #
  50 #################################
  51 */
  52 PUBLIC  int logML     = 0;  /* if nonzero use logarithmic ML energy in energy_of_struct */
  53 PUBLIC  int uniq_ML   = 0;  /* do ML decomposition uniquely (for subopt) */
  54 PUBLIC  int cut_point = -1; /* set to first pos of second seq for cofolding */
  55 PUBLIC  int eos_debug = 0;  /* verbose info from energy_of_struct */
  56
  57 /*
  58 #################################
  59 # PRIVATE VARIABLES             #
  60 #################################
  61 */
  62 PRIVATE int     *indx     = NULL; /* index for moving in the triangle matrices c[] and fMl[]*/
  63 PRIVATE int     *c        = NULL; /* energy array, given that i-j pair */
  64 PRIVATE int     *cc       = NULL; /* linear array for calculating canonical structures */
  65 PRIVATE int     *cc1      = NULL; /*   "     "        */
  66 PRIVATE int     *f5       = NULL; /* energy of 5' end */
  67 PRIVATE int     *f53      = NULL; /* energy of 5' end with 3' nucleotide not available for mismatches */
  68 PRIVATE int     *fML      = NULL; /* multi-loop auxiliary energy array */
  69 PRIVATE int     *fM1      = NULL; /* second ML array, only for subopt */
  70 PRIVATE int     *fM2      = NULL; /* fM2 = multiloop region with exactly two stems, extending to 3' end        */
  71 PRIVATE int     *Fmi      = NULL; /* holds row i of fML (avoids jumps in memory) */
  72 PRIVATE int     *DMLi     = NULL; /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
  73 PRIVATE int     *DMLi1    = NULL; /*             MIN(fML[i+1,k]+fML[k+1,j])  */
  74 PRIVATE int     *DMLi2    = NULL; /*             MIN(fML[i+2,k]+fML[k+1,j])  */
  75 PRIVATE int     *DMLi_a   = NULL; /* DMLi_a[j] holds min energy for at least two multiloop stems in [i,j], where j is available for dangling onto a surrounding stem */
  76 PRIVATE int     *DMLi_o   = NULL; /* DMLi_o[j] holds min energy for at least two multiloop stems in [i,j], where j is unavailable for dangling onto a surrounding stem */
  77 PRIVATE int     *DMLi1_a  = NULL;
  78 PRIVATE int     *DMLi1_o  = NULL;
  79 PRIVATE int     *DMLi2_a  = NULL;
  80 PRIVATE int     *DMLi2_o  = NULL;
  81 PRIVATE int     Fc, FcH, FcI, FcM;  /* parts of the exterior loop energies */
  82 PRIVATE sect    sector[MAXSECTORS]; /* stack of partial structures for backtracking */
  83 PRIVATE char    *ptype = NULL;      /* precomputed array of pair types */
  84 PRIVATE short   *S = NULL, *S1 = NULL;
  85 PRIVATE paramT  *P          = NULL;
  86 PRIVATE int     init_length = -1;
  87 PRIVATE int     *BP = NULL; /* contains the structure constrainsts: BP[i]
  88                         -1: | = base must be paired
  89                         -2: < = base must be paired with j<i
  90                         -3: > = base must be paired with j>i
  91                         -4: x = base must not pair
  92                         positive int: base is paired with int      */
  93 PRIVATE short   *pair_table         = NULL; /* needed by energy of struct */
  94 PRIVATE bondT   *base_pair2         = NULL; /* this replaces base_pair from fold_vars.c */
  95 PRIVATE int     circular            = 0;
  96 PRIVATE int     struct_constrained  = 0;
  97 PRIVATE int     with_gquad          = 0;
  98
  99 PRIVATE int     *ggg = NULL;    /* minimum free energies of the gquadruplexes */
 100
 101 #ifdef _OPENMP
 102
 103 #pragma omp threadprivate(indx, c, cc, cc1, f5, f53, fML, fM1, fM2, Fmi,\
 104                           DMLi, DMLi1, DMLi2, DMLi_a, DMLi_o, DMLi1_a, DMLi1_o, DMLi2_a, DMLi2_o,\
 105                           Fc, FcH, FcI, FcM,\
 106                           sector, ptype, S, S1, P, init_length, BP, pair_table, base_pair2, circular, struct_constrained,\
 107                           ggg, with_gquad)
 108
 109 #endif
 110
 111 /*
 112 #################################
 113 # PRIVATE FUNCTION DECLARATIONS #
 114 #################################
 115 */
 116 PRIVATE void  get_arrays(unsigned int size);
 117 PRIVATE int   stack_energy(int i, const char *string, int verbostiy_level);
 118 PRIVATE int   energy_of_extLoop_pt(int i, short *pair_table);
 119 PRIVATE int   energy_of_ml_pt(int i, short *pt);
 120 PRIVATE int   ML_Energy(int i, int is_extloop);
 121 PRIVATE void  make_ptypes(const short *S, const char *structure);
 122 PRIVATE void  backtrack(const char *sequence, int s);
 123 PRIVATE int   fill_arrays(const char *sequence);
 124 PRIVATE void  fill_arrays_circ(const char *string, int *bt);
 125 PRIVATE void  init_fold(int length, paramT *parameters);
 126 /* needed by cofold/eval */
 127 PRIVATE int   cut_in_loop(int i);
 128
 129 /* deprecated functions */
 130 /*@unused@*/
 131 int oldLoopEnergy(int i, int j, int p, int q, int type, int type_2);
 132 int LoopEnergy(int n1, int n2, int type, int type_2, int si1, int sj1, int sp1, int sq1);
 133 int HairpinE(int size, int type, int si1, int sj1, const char *string);
 134
 135
 136 /*
 137 #################################
 138 # BEGIN OF FUNCTION DEFINITIONS #
 139 #################################
 140 */
 141
 142 /* allocate memory for folding process */
 143 PRIVATE void init_fold(int length, paramT *parameters){
 144
 145 #ifdef _OPENMP
 146 /* Explicitly turn off dynamic threads */
 147   omp_set_dynamic(0);
 148 #endif
 149
 150   if (length<1) nrerror("initialize_fold: argument must be greater 0");
 151   free_arrays();
 152   get_arrays((unsigned) length);
 153   init_length=length;
 154
 155   indx = get_indx((unsigned)length);
 156
 157   update_fold_params_par(parameters);
 158 }
 159
 160 /*--------------------------------------------------------------------------*/
 161
 162 PRIVATE void get_arrays(unsigned int size){
 163   if(size >= (unsigned int)sqrt((double)INT_MAX))
 164     nrerror("get_arrays@fold.c: sequence length exceeds addressable range");
 165
 166   c     = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 167   fML   = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 168   if (uniq_ML)
 169     fM1 = (int *) space(sizeof(int)*((size*(size+1))/2+2));
 170
 171   ptype = (char *)space(sizeof(char)*((size*(size+1))/2+2));
 172   f5    = (int *) space(sizeof(int)*(size+2));
 173   f53   = (int *) space(sizeof(int)*(size+2));
 174   cc    = (int *) space(sizeof(int)*(size+2));
 175   cc1   = (int *) space(sizeof(int)*(size+2));
 176   Fmi   = (int *) space(sizeof(int)*(size+1));
 177   DMLi  = (int *) space(sizeof(int)*(size+1));
 178   DMLi1 = (int *) space(sizeof(int)*(size+1));
 179   DMLi2 = (int *) space(sizeof(int)*(size+1));
 180
 181   DMLi_a  = (int *) space(sizeof(int)*(size+1));
 182   DMLi_o  = (int *) space(sizeof(int)*(size+1));
 183   DMLi1_a = (int *) space(sizeof(int)*(size+1));
 184   DMLi1_o = (int *) space(sizeof(int)*(size+1));
 185   DMLi2_a = (int *) space(sizeof(int)*(size+1));
 186   DMLi2_o = (int *) space(sizeof(int)*(size+1));
 187
 188   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 */
 189
 190   /* extra array(s) for circfold() */
 191   if(circular) fM2 =  (int *) space(sizeof(int)*(size+2));
 192 }
 193
 194 /*--------------------------------------------------------------------------*/
 195
 196 PUBLIC void free_arrays(void){
 197   if(indx)      free(indx);
 198   if(c)         free(c);
 199   if(fML)       free(fML);
 200   if(f5)        free(f5);
 201   if(f53)       free(f53);
 202   if(cc)        free(cc);
 203   if(cc1)       free(cc1);
 204   if(ptype)     free(ptype);
 205   if(fM1)       free(fM1);
 206   if(fM2)       free(fM2);
 207   if(base_pair2) free(base_pair2);
 208   if(Fmi)       free(Fmi);
 209   if(DMLi)      free(DMLi);
 210   if(DMLi1)     free(DMLi1);
 211   if(DMLi2)     free(DMLi2);
 212   if(DMLi_a)    free(DMLi_a);
 213   if(DMLi_o)    free(DMLi_o);
 214   if(DMLi1_a)   free(DMLi1_a);
 215   if(DMLi1_o)   free(DMLi1_o);
 216   if(DMLi2_a)   free(DMLi2_a);
 217   if(DMLi2_o)   free(DMLi2_o);
 218   if(P)         free(P);
 219   if(ggg)       free(ggg);
 220
 221   indx = c = fML = f5 = f53 = cc = cc1 = fM1 = fM2 = Fmi = DMLi = DMLi1 = DMLi2 = ggg = NULL;
 222   DMLi_a = DMLi_o = DMLi1_a = DMLi1_o = DMLi2_a = DMLi2_o = NULL;
 223   ptype       = NULL;
 224   base_pair   = NULL;
 225   base_pair2  = NULL;
 226   P           = NULL;
 227   init_length = 0;
 228 }
 229
 230 /*--------------------------------------------------------------------------*/
 231
 232 PUBLIC void export_fold_arrays( int **f5_p,
 233                                 int **c_p,
 234                                 int **fML_p,
 235                                 int **fM1_p,
 236                                 int **indx_p,
 237                                 char **ptype_p){
 238   /* make the DP arrays available to routines such as subopt() */
 239   *f5_p     = f5;
 240   *c_p      = c;
 241   *fML_p    = fML;
 242   *fM1_p    = fM1;
 243   *indx_p   = indx;
 244   *ptype_p  = ptype;
 245 }
 246
 247 PUBLIC void export_fold_arrays_par( int **f5_p,
 248                                     int **c_p,
 249                                     int **fML_p,
 250                                     int **fM1_p,
 251                                     int **indx_p,
 252                                     char **ptype_p,
 253                                     paramT **P_p){
 254   export_fold_arrays(f5_p, c_p, fML_p, fM1_p, indx_p,ptype_p);
 255   *P_p = P;
 256 }
 257
 258 PUBLIC void export_circfold_arrays( int *Fc_p,
 259                                     int *FcH_p,
 260                                     int *FcI_p,
 261                                     int *FcM_p,
 262                                     int **fM2_p,
 263                                     int **f5_p,
 264                                     int **c_p,
 265                                     int **fML_p,
 266                                     int **fM1_p,
 267                                     int **indx_p,
 268                                     char **ptype_p){
 269   /* make the DP arrays available to routines such as subopt() */
 270   *f5_p     = f5;
 271   *c_p      = c;
 272   *fML_p    = fML;
 273   *fM1_p    = fM1;
 274   *fM2_p    = fM2;
 275   *Fc_p     = Fc;
 276   *FcH_p    = FcH;
 277   *FcI_p    = FcI;
 278   *FcM_p    = FcM;
 279   *indx_p   = indx;
 280   *ptype_p  = ptype;
 281 }
 282
 283 PUBLIC void export_circfold_arrays_par( int *Fc_p,
 284                                     int *FcH_p,
 285                                     int *FcI_p,
 286                                     int *FcM_p,
 287                                     int **fM2_p,
 288                                     int **f5_p,
 289                                     int **c_p,
 290                                     int **fML_p,
 291                                     int **fM1_p,
 292                                     int **indx_p,
 293                                     char **ptype_p,
 294                                     paramT **P_p){
 295   export_circfold_arrays(Fc_p, FcH_p, FcI_p, FcM_p, fM2_p, f5_p, c_p, fML_p, fM1_p, indx_p, ptype_p);
 296   *P_p = P;
 297 }
 298 /*--------------------------------------------------------------------------*/
 299
 300
 301 PUBLIC float fold(const char *string, char *structure){
 302   return fold_par(string, structure, NULL, fold_constrained, 0);
 303 }
 304
 305 PUBLIC float circfold(const char *string, char *structure){
 306   return fold_par(string, structure, NULL, fold_constrained, 1);
 307 }
 308
 309 PUBLIC float fold_par(const char *string,
 310                       char *structure,
 311                       paramT *parameters,
 312                       int is_constrained,
 313                       int is_circular){
 314
 315   int i, length, energy, bonus, bonus_cnt, s;
 316
 317   bonus               = 0;
 318   bonus_cnt           = 0;
 319   s                   = 0;
 320   circular            = is_circular;
 321   struct_constrained  = is_constrained;
 322   length              = (int) strlen(string);
 323
 324 #ifdef _OPENMP
 325   init_fold(length, parameters);
 326 #else
 327   if (parameters) init_fold(length, parameters);
 328   else if (length>init_length) init_fold(length, parameters);
 329   else if (fabs(P->temperature - temperature)>1e-6) update_fold_params();
 330 #endif
 331
 332   with_gquad  = P->model_details.gquad;
 333   S           = encode_sequence(string, 0);
 334   S1          = encode_sequence(string, 1);
 335   BP          = (int *)space(sizeof(int)*(length+2));
 336
 337   make_ptypes(S, structure);
 338
 339   energy = fill_arrays(string);
 340
 341   if(circular){
 342     fill_arrays_circ(string, &s);
 343     energy = Fc;
 344   }
 345   backtrack(string, s);
 346
 347 #ifdef PAREN
 348   parenthesis_structure(structure, base_pair2, length);
 349 #else
 350   letter_structure(structure, base_pair2, length);
 351 #endif
 352
 353   /*
 354   *  Backward compatibility:
 355   *  This block may be removed if deprecated functions
 356   *  relying on the global variable "base_pair" vanishs from within the package!
 357   */
 358   base_pair = base_pair2;
 359   /*
 360   {
 361     if(base_pair) free(base_pair);
 362     base_pair = (bondT *)space(sizeof(bondT) * (1+length/2));
 363     memcpy(base_pair, base_pair2, sizeof(bondT) * (1+length/2));
 364   }
 365   */
 366
 367   /* check constraints */
 368   for(i=1;i<=length;i++) {
 369     if((BP[i]<0)&&(BP[i]>-4)) {
 370       bonus_cnt++;
 371       if((BP[i]==-3)&&(structure[i-1]==')')) bonus++;
 372       if((BP[i]==-2)&&(structure[i-1]=='(')) bonus++;
 373       if((BP[i]==-1)&&(structure[i-1]!='.')) bonus++;
 374     }
 375
 376     if(BP[i]>i) {
 377       int l;
 378       bonus_cnt++;
 379       for(l=1; l<=base_pair2[0].i; l++)
 380         if(base_pair2[l].i != base_pair2[l].j)
 381           if((i==base_pair2[l].i)&&(BP[i]==base_pair2[l].j)) bonus++;
 382     }
 383   }
 384
 385   if (bonus_cnt>bonus) fprintf(stderr,"\ncould not enforce all constraints\n");
 386   bonus*=BONUS;
 387
 388   free(S); free(S1); free(BP);
 389
 390   energy += bonus;      /*remove bonus energies from result */
 391
 392   if (backtrack_type=='C')
 393     return (float) c[indx[length]+1]/100.;
 394   else if (backtrack_type=='M')
 395     return (float) fML[indx[length]+1]/100.;
 396   else
 397     return (float) energy/100.;
 398 }
 399
 400 /**
 401 *** fill "c", "fML" and "f5" arrays and return  optimal energy
 402 **/
 403 PRIVATE int fill_arrays(const char *string) {
 404
 405   int   i, j, k, length, energy, en, mm5, mm3;
 406   int   decomp, new_fML, max_separation;
 407   int   no_close, type, type_2, tt;
 408   int   bonus=0;
 409
 410   int   dangle_model, noGUclosure, with_gquads;
 411
 412   dangle_model  = P->model_details.dangles;
 413   noGUclosure   = P->model_details.noGUclosure;
 414
 415   length = (int) strlen(string);
 416
 417   max_separation = (int) ((1.-LOCALITY)*(double)(length-2)); /* not in use */
 418
 419   if(with_gquad)
 420     ggg = get_gquad_matrix(S, P);
 421
 422
 423   for (j=1; j<=length; j++) {
 424     Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
 425   }
 426
 427   for (j = 1; j<=length; j++)
 428     for (i=(j>TURN?(j-TURN):1); i<j; i++) {
 429       c[indx[j]+i] = fML[indx[j]+i] = INF;
 430       if (uniq_ML) fM1[indx[j]+i] = INF;
 431     }
 432
 433   if (length <= TURN) return 0;
 434
 435   for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
 436
 437     for (j = i+TURN+1; j <= length; j++) {
 438       int p, q, ij, jj, ee;
 439       int minq, maxq, l1, up, c0, c1, c2, c3;
 440       int MLenergy;
 441       ij = indx[j]+i;
 442       bonus = 0;
 443       type = ptype[ij];
 444       energy = INF;
 445       /* enforcing structure constraints */
 446       if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
 447       if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
 448       if ((BP[i]==-4)||(BP[j]==-4)) type=0;
 449
 450       no_close = (((type==3)||(type==4))&&noGUclosure&&(bonus==0));
 451
 452       if (j-i-1 > max_separation) type = 0;  /* forces locality degree */
 453
 454       if (type) {   /* we have a pair */
 455         int new_c=0, stackEnergy=INF;
 456         /* hairpin ----------------------------------------------*/
 457
 458         new_c = (no_close) ? FORBIDDEN : E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+i-1, P);
 459
 460         /*--------------------------------------------------------
 461           check for elementary structures involving more than one
 462           closing pair.
 463           --------------------------------------------------------*/
 464
 465         for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++) {
 466           minq = j-i+p-MAXLOOP-2;
 467           if (minq<p+1+TURN) minq = p+1+TURN;
 468           for (q = minq; q < j; q++) {
 469             type_2 = ptype[indx[q]+p];
 470
 471             if (type_2==0) continue;
 472             type_2 = rtype[type_2];
 473
 474             if (noGUclosure)
 475               if (no_close||(type_2==3)||(type_2==4))
 476                 if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
 477
 478             energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
 479                                 S1[i+1], S1[j-1], S1[p-1], S1[q+1], P);
 480
 481             ee = energy+c[indx[q]+p];
 482             new_c = MIN2(new_c, ee);
 483             if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
 484
 485           } /* end q-loop */
 486         } /* end p-loop */
 487         /* multi-loop decomposition ------------------------*/
 488
 489
 490         if (!no_close) {
 491           decomp = DMLi1[j-1];
 492           tt = rtype[type];
 493           switch(dangle_model){
 494             /* no dangles */
 495             case 0:   decomp += E_MLstem(tt, -1, -1, P);
 496                       break;
 497
 498             /* double dangles */
 499             case 2:   decomp += E_MLstem(tt, S1[j-1], S1[i+1], P);
 500                       break;
 501
 502             /* normal dangles, aka dangles = 1 || 3 */
 503             default:  decomp += E_MLstem(tt, -1, -1, P);
 504                       decomp = MIN2(decomp, DMLi2[j-1] + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase);
 505                       decomp = MIN2(decomp, DMLi2[j-2] + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase);
 506                       decomp = MIN2(decomp, DMLi1[j-2] + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase);
 507                       break;
 508           }
 509           MLenergy = decomp + P->MLclosing;
 510           new_c = MIN2(new_c, MLenergy);
 511         }
 512
 513         /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
 514
 515         if (dangle_model==3) {
 516           decomp = INF;
 517           for (k = i+2+TURN; k < j-2-TURN; k++) {
 518             type_2 = rtype[ptype[indx[k]+i+1]];
 519             if (type_2)
 520               decomp = MIN2(decomp, c[indx[k]+i+1]+P->stack[type][type_2]+fML[indx[j-1]+k+1]);
 521             type_2 = rtype[ptype[indx[j-1]+k+1]];
 522             if (type_2)
 523               decomp = MIN2(decomp, c[indx[j-1]+k+1]+P->stack[type][type_2]+fML[indx[k]+i+1]);
 524           }
 525           /* no TermAU penalty if coax stack */
 526           decomp += 2*P->MLintern[1] + P->MLclosing;
 527           new_c = MIN2(new_c, decomp);
 528         }
 529
 530         if(with_gquad){
 531           /* include all cases where a g-quadruplex may be enclosed by base pair (i,j) */
 532           if (!no_close) {
 533             tt = rtype[type];
 534             energy = E_GQuad_IntLoop(i, j, type, S1, ggg, indx, P);
 535             new_c = MIN2(new_c, energy);
 536           }
 537         }
 538
 539         new_c = MIN2(new_c, cc1[j-1]+stackEnergy);
 540         cc[j] = new_c + bonus;
 541         if (noLonelyPairs)
 542           c[ij] = cc1[j-1]+stackEnergy+bonus;
 543         else
 544           c[ij] = cc[j];
 545
 546       } /* end >> if (pair) << */
 547
 548       else c[ij] = INF;
 549
 550       /* done with c[i,j], now compute fML[i,j] and fM1[i,j] */
 551
 552       /* (i,j) + MLstem ? */
 553       new_fML = INF;
 554       if(type){
 555         new_fML = c[ij];
 556         switch(dangle_model){
 557           case 2:   new_fML += E_MLstem(type, (i==1) ? S1[length] : S1[i-1], S1[j+1], P);
 558                     break;
 559           default:  new_fML += E_MLstem(type, -1, -1, P);
 560                     break;
 561         }
 562       }
 563
 564       if(with_gquad){
 565         new_fML = MIN2(new_fML, ggg[indx[j] + i] + E_MLstem(0, -1, -1, P));
 566       }
 567
 568       if (uniq_ML){
 569         fM1[ij] = MIN2(fM1[indx[j-1]+i] + P->MLbase, new_fML);
 570       }
 571
 572       /* free ends ? -----------------------------------------*/
 573       /*  we must not just extend 3'/5' end by unpaired nucleotides if
 574       *   dangle_model == 1, this could lead to d5+d3 contributions were
 575       *   mismatch must be taken!
 576       */
 577       switch(dangle_model){
 578         /* no dangles */
 579         case 0:   new_fML = MIN2(new_fML, fML[ij+1]+P->MLbase);
 580                   new_fML = MIN2(fML[indx[j-1]+i]+P->MLbase, new_fML);
 581                   break;
 582
 583         /* double dangles */
 584         case 2:   new_fML = MIN2(new_fML, fML[ij+1]+P->MLbase);
 585                   new_fML = MIN2(fML[indx[j-1]+i]+P->MLbase, new_fML);
 586                   break;
 587
 588         /* normal dangles, aka dangle_model = 1 || 3 */
 589         default:  mm5 = ((i>1) || circular) ? S1[i] : -1;
 590                   mm3 = ((j<length) || circular) ? S1[j] : -1;
 591                   new_fML = MIN2(new_fML, fML[ij+1] + P->MLbase);
 592                   new_fML = MIN2(new_fML, fML[indx[j-1]+i] + P->MLbase);
 593                   tt = ptype[ij+1];
 594                   if(tt) new_fML = MIN2(new_fML, c[ij+1] + E_MLstem(tt, mm5, -1, P) + P->MLbase);
 595                   tt = ptype[indx[j-1]+i];
 596                   if(tt) new_fML = MIN2(new_fML, c[indx[j-1]+i] + E_MLstem(tt, -1, mm3, P) + P->MLbase);
 597                   tt = ptype[indx[j-1]+i+1];
 598                   if(tt) new_fML = MIN2(new_fML, c[indx[j-1]+i+1] + E_MLstem(tt, mm5, mm3, P) + 2*P->MLbase);
 599                   break;
 600       }
 601
 602       /* modular decomposition -------------------------------*/
 603       for (decomp = INF, k = i + 1 + TURN; k <= j - 2 - TURN; k++)
 604         decomp = MIN2(decomp, Fmi[k]+fML[indx[j]+k+1]);
 605       DMLi[j] = decomp;               /* store for use in ML decompositon */
 606       new_fML = MIN2(new_fML,decomp);
 607
 608       /* coaxial stacking */
 609       if (dangle_model==3) {
 610         /* additional ML decomposition as two coaxially stacked helices */
 611         for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++) {
 612           type = ptype[indx[k]+i]; type = rtype[type];
 613           type_2 = ptype[indx[j]+k+1]; type_2 = rtype[type_2];
 614           if (type && type_2)
 615             decomp = MIN2(decomp,
 616                           c[indx[k]+i]+c[indx[j]+k+1]+P->stack[type][type_2]);
 617         }
 618
 619         decomp += 2*P->MLintern[1];        /* no TermAU penalty if coax stack */
 620 #if 0
 621         /* This is needed for Y shaped ML loops with coax stacking of
 622            interior pairts, but backtracking will fail if activated */
 623         DMLi[j] = MIN2(DMLi[j], decomp);
 624         DMLi[j] = MIN2(DMLi[j], DMLi[j-1]+P->MLbase);
 625         DMLi[j] = MIN2(DMLi[j], DMLi1[j]+P->MLbase);
 626         new_fML = MIN2(new_fML, DMLi[j]);
 627 #endif
 628         new_fML = MIN2(new_fML, decomp);
 629       }
 630       fML[ij] = Fmi[j] = new_fML;     /* substring energy */
 631
 632     }
 633
 634     {
 635       int *FF; /* rotate the auxilliary arrays */
 636       FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
 637       FF = cc1; cc1=cc; cc=FF;
 638       for (j=1; j<=length; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
 639     }
 640   }
 641
 642   /* calculate energies of 5' and 3' fragments */
 643
 644   f5[TURN+1]= 0;
 645   /* duplicated code may be faster than conditions inside loop ;) */
 646   switch(dangle_model){
 647     /* dont use dangling end and mismatch contributions at all */
 648     case 0:   for(j=TURN+2; j<=length; j++){
 649                 f5[j] = f5[j-1];
 650                 for (i=j-TURN-1; i>1; i--){
 651
 652                   if(with_gquad){
 653                     f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
 654                   }
 655
 656                   type = ptype[indx[j]+i];
 657                   if(!type) continue;
 658                   en = c[indx[j]+i];
 659                   f5[j] = MIN2(f5[j], f5[i-1] + en + E_ExtLoop(type, -1, -1, P));
 660                 }
 661
 662                 if(with_gquad){
 663                   f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
 664                 }
 665
 666                 type=ptype[indx[j]+1];
 667                 if(!type) continue;
 668                 en = c[indx[j]+1];
 669                 f5[j] = MIN2(f5[j], en + E_ExtLoop(type, -1, -1, P));
 670               }
 671               break;
 672
 673     /* always use dangles on both sides */
 674     case 2:   for(j=TURN+2; j<length; j++){
 675                 f5[j] = f5[j-1];
 676                 for (i=j-TURN-1; i>1; i--){
 677
 678                   if(with_gquad){
 679                     f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
 680                   }
 681
 682                   type = ptype[indx[j]+i];
 683                   if(!type) continue;
 684                   en = c[indx[j]+i];
 685                   f5[j] = MIN2(f5[j], f5[i-1] + en + E_ExtLoop(type, S1[i-1], S1[j+1], P));
 686                 }
 687
 688                 if(with_gquad){
 689                   f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
 690                 }
 691
 692                 type=ptype[indx[j]+1];
 693                 if(!type) continue;
 694                 en = c[indx[j]+1];
 695                 f5[j] = MIN2(f5[j], en + E_ExtLoop(type, -1, S1[j+1], P));
 696               }
 697               f5[length] = f5[length-1];
 698               for (i=length-TURN-1; i>1; i--){
 699
 700                 if(with_gquad){
 701                   f5[length] = MIN2(f5[length], f5[i-1] + ggg[indx[length]+i]);
 702                 }
 703
 704                 type = ptype[indx[length]+i];
 705                 if(!type) continue;
 706                 en = c[indx[length]+i];
 707                 f5[length] = MIN2(f5[length], f5[i-1] + en + E_ExtLoop(type, S1[i-1], -1, P));
 708               }
 709
 710               if(with_gquad){
 711                 f5[length] = MIN2(f5[length], ggg[indx[length]+1]);
 712               }
 713
 714               type=ptype[indx[length]+1];
 715               if(!type) break;
 716               en = c[indx[length]+1];
 717               f5[length] = MIN2(f5[length], en + E_ExtLoop(type, -1, -1, P));
 718
 719
 720               break;
 721
 722     /* normal dangles, aka dangle_model = 1 || 3 */
 723     default:  for(j=TURN+2; j<=length; j++){
 724                 f5[j] = f5[j-1];
 725                 for (i=j-TURN-1; i>1; i--){
 726
 727                   if(with_gquad){
 728                     f5[j] = MIN2(f5[j], f5[i-1] + ggg[indx[j]+i]);
 729                   }
 730
 731                   type = ptype[indx[j]+i];
 732                   if(type){
 733                     en = c[indx[j]+i];
 734                     f5[j] = MIN2(f5[j], f5[i-1] + en + E_ExtLoop(type, -1, -1, P));
 735                     f5[j] = MIN2(f5[j], f5[i-2] + en + E_ExtLoop(type, S1[i-1], -1, P));
 736                   }
 737                   type = ptype[indx[j-1]+i];
 738                   if(type){
 739                     en = c[indx[j-1]+i];
 740                     f5[j] = MIN2(f5[j], f5[i-1] + en + E_ExtLoop(type, -1, S1[j], P));
 741                     f5[j] = MIN2(f5[j], f5[i-2] + en + E_ExtLoop(type, S1[i-1], S1[j], P));
 742                   }
 743                 }
 744
 745                 if(with_gquad){
 746                   f5[j] = MIN2(f5[j], ggg[indx[j]+1]);
 747                 }
 748
 749                 type = ptype[indx[j]+1];
 750                 if(type) f5[j] = MIN2(f5[j], c[indx[j]+1] + E_ExtLoop(type, -1, -1, P));
 751                 type = ptype[indx[j-1]+1];
 752                 if(type) f5[j] = MIN2(f5[j], c[indx[j-1]+1] + E_ExtLoop(type, -1, S1[j], P));
 753               }
 754   }
 755   return f5[length];
 756 }
 757
 758 #include "circfold.inc"
 759
 760 /**
 761 *** trace back through the "c", "f5" and "fML" arrays to get the
 762 *** base pairing list. No search for equivalent structures is done.
 763 *** This is fast, since only few structure elements are recalculated.
 764 ***
 765 *** normally s=0.
 766 *** If s>0 then s items have been already pushed onto the sector stack
 767 **/
 768 PRIVATE void backtrack(const char *string, int s) {
 769   int   i, j, ij, k, l1, mm5, mm3, length, energy, en, new;
 770   int   no_close, type, type_2, tt, minq, maxq, c0, c1, c2, c3;
 771   int   bonus;
 772   int   b=0;
 773   int   dangle_model = P->model_details.dangles;
 774
 775   length = strlen(string);
 776   if (s==0) {
 777     sector[++s].i = 1;
 778     sector[s].j = length;
 779     sector[s].ml = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')? 2: 0);
 780   }
 781   while (s>0) {
 782     int ml, fij, fi, cij, traced, i1, j1, p, q, jj=0, gq=0;
 783     int canonical = 1;     /* (i,j) closes a canonical structure */
 784     i  = sector[s].i;
 785     j  = sector[s].j;
 786     ml = sector[s--].ml;   /* ml is a flag indicating if backtracking is to
 787                               occur in the fML- (1) or in the f-array (0) */
 788     if (ml==2) {
 789       base_pair2[++b].i = i;
 790       base_pair2[b].j   = j;
 791       goto repeat1;
 792     }
 793
 794     else if(ml==7) { /* indicates that i,j are enclosing a gquadruplex */
 795       /* actually, do something here */
 796     }
 797
 798     if (j < i+TURN+1) continue; /* no more pairs in this interval */
 799
 800     fij = (ml == 1)? fML[indx[j]+i] : f5[j];
 801     fi  = (ml == 1)?(fML[indx[j-1]+i]+P->MLbase): f5[j-1];
 802
 803     if (fij == fi) {  /* 3' end is unpaired */
 804       sector[++s].i = i;
 805       sector[s].j   = j-1;
 806       sector[s].ml  = ml;
 807       continue;
 808     }
 809
 810     if (ml == 0) { /* backtrack in f5 */
 811       switch(dangle_model){
 812         case 0:   /* j is paired. Find pairing partner */
 813                   for(k=j-TURN-1,traced=0; k>=1; k--){
 814
 815                     if(with_gquad){
 816                       if(fij == f5[k-1] + ggg[indx[j]+k]){
 817                         /* found the decomposition */
 818                         traced = j; jj = k - 1; gq = 1;
 819                         break;
 820                       }
 821                     }
 822
 823                     type = ptype[indx[j]+k];
 824                     if(type)
 825                       if(fij == E_ExtLoop(type, -1, -1, P) + c[indx[j]+k] + f5[k-1]){
 826                         traced=j; jj = k-1;
 827                         break;
 828                       }
 829                   }
 830                   break;
 831
 832         case 2:   mm3 = (j<length) ? S1[j+1] : -1;
 833                   for(k=j-TURN-1,traced=0; k>=1; k--){
 834
 835                     if(with_gquad){
 836                       if(fij == f5[k-1] + ggg[indx[j]+k]){
 837                         /* found the decomposition */
 838                         traced = j; jj = k - 1; gq = 1;
 839                         break;
 840                       }
 841                     }
 842
 843                     type = ptype[indx[j]+k];
 844                     if(type)
 845                       if(fij == E_ExtLoop(type, (k>1) ? S1[k-1] : -1, mm3, P) + c[indx[j]+k] + f5[k-1]){
 846                         traced=j; jj = k-1;
 847                         break;
 848                       }
 849                   }
 850                   break;
 851
 852         default:  for(traced = 0, k=j-TURN-1; k>1; k--){
 853
 854                     if(with_gquad){
 855                       if(fij == f5[k-1] + ggg[indx[j]+k]){
 856                         /* found the decomposition */
 857                         traced = j; jj = k - 1; gq = 1;
 858                         break;
 859                       }
 860                     }
 861
 862                     type = ptype[indx[j] + k];
 863                     if(type){
 864                       en = c[indx[j] + k];
 865                       if(fij == f5[k-1] + en + E_ExtLoop(type, -1, -1, P)){
 866                         traced = j;
 867                         jj = k-1;
 868                         break;
 869                       }
 870                       if(fij == f5[k-2] + en + E_ExtLoop(type, S1[k-1], -1, P)){
 871                         traced = j;
 872                         jj = k-2;
 873                         break;
 874                       }
 875                     }
 876                     type = ptype[indx[j-1] + k];
 877                     if(type){
 878                       en = c[indx[j-1] + k];
 879                       if(fij == f5[k-1] + en + E_ExtLoop(type, -1, S1[j], P)){
 880                         traced = j-1;
 881                         jj = k-1;
 882                         break;
 883                       }
 884                       if(fij == f5[k-2] + en + E_ExtLoop(type, S1[k-1], S1[j], P)){
 885                         traced = j-1;
 886                         jj = k-2;
 887                         break;
 888                       }
 889                     }
 890                   }
 891                   if(!traced){
 892
 893                     if(with_gquad){
 894                       if(fij == ggg[indx[j]+1]){
 895                         /* found the decomposition */
 896                         traced = j; jj = 0; gq = 1;
 897                         break;
 898                       }
 899                     }
 900
 901                     type = ptype[indx[j]+1];
 902                     if(type){
 903                       if(fij == c[indx[j]+1] + E_ExtLoop(type, -1, -1, P)){
 904                         traced = j;
 905                         jj = 0;
 906                         break;
 907                       }
 908                     }
 909                     type = ptype[indx[j-1]+1];
 910                     if(type){
 911                       if(fij == c[indx[j-1]+1] + E_ExtLoop(type, -1, S1[j], P)){
 912                         traced = j-1;
 913                         jj = 0;
 914                         break;
 915                       }
 916                     }
 917                   }
 918                   break;
 919       }
 920
 921       if (!traced){
 922         fprintf(stderr, "%s\n", string);
 923         nrerror("backtrack failed in f5");
 924       }
 925       /* push back the remaining f5 portion */
 926       sector[++s].i = 1;
 927       sector[s].j   = jj;
 928       sector[s].ml  = ml;
 929
 930       /* trace back the base pair found */
 931       i=k; j=traced;
 932
 933       if(with_gquad && gq){
 934         /* goto backtrace of gquadruplex */
 935         goto repeat_gquad;
 936       }
 937
 938       base_pair2[++b].i = i;
 939       base_pair2[b].j   = j;
 940       goto repeat1;
 941     }
 942     else { /* trace back in fML array */
 943       if (fML[indx[j]+i+1]+P->MLbase == fij) { /* 5' end is unpaired */
 944         sector[++s].i = i+1;
 945         sector[s].j   = j;
 946         sector[s].ml  = ml;
 947         continue;
 948       }
 949
 950       ij  = indx[j]+i;
 951
 952       if(with_gquad){
 953         if(fij == ggg[ij] + E_MLstem(0, -1, -1, P)){
 954           /* go to backtracing of quadruplex */
 955           goto repeat_gquad;
 956         }
 957       }
 958
 959       tt  = ptype[ij];
 960       en  = c[ij];
 961       switch(dangle_model){
 962         case 0:   if(fij == en + E_MLstem(tt, -1, -1, P)){
 963                     base_pair2[++b].i = i;
 964                     base_pair2[b].j   = j;
 965                     goto repeat1;
 966                   }
 967                   break;
 968
 969         case 2:   if(fij == en + E_MLstem(tt, S1[i-1], S1[j+1], P)){
 970                     base_pair2[++b].i = i;
 971                     base_pair2[b].j   = j;
 972                     goto repeat1;
 973                   }
 974                   break;
 975
 976         default:  if(fij == en + E_MLstem(tt, -1, -1, P)){
 977                     base_pair2[++b].i = i;
 978                     base_pair2[b].j   = j;
 979                     goto repeat1;
 980                   }
 981                   tt = ptype[ij+1];
 982                   if(fij == c[ij+1] + E_MLstem(tt, S1[i], -1, P) + P->MLbase){
 983                     base_pair2[++b].i = ++i;
 984                     base_pair2[b].j   = j;
 985                     goto repeat1;
 986                   }
 987                   tt = ptype[indx[j-1]+i];
 988                   if(fij == c[indx[j-1]+i] + E_MLstem(tt, -1, S1[j], P) + P->MLbase){
 989                     base_pair2[++b].i = i;
 990                     base_pair2[b].j   = --j;
 991                     goto repeat1;
 992                   }
 993                   tt = ptype[indx[j-1]+i+1];
 994                   if(fij == c[indx[j-1]+i+1] + E_MLstem(tt, S1[i], S1[j], P) + 2*P->MLbase){
 995                     base_pair2[++b].i = ++i;
 996                     base_pair2[b].j   = --j;
 997                     goto repeat1;
 998                   }
 999                   break;
1000       }
1001
1002       for(k = i + 1 + TURN; k <= j - 2 - TURN; k++)
1003         if(fij == (fML[indx[k]+i]+fML[indx[j]+k+1]))
1004           break;
1005
1006       if ((dangle_model==3)&&(k > j - 2 - TURN)) { /* must be coax stack */
1007         ml = 2;
1008         for (k = i+1+TURN; k <= j - 2 - TURN; k++) {
1009           type    = rtype[ptype[indx[k]+i]];
1010           type_2  = rtype[ptype[indx[j]+k+1]];
1011           if (type && type_2)
1012             if (fij == c[indx[k]+i]+c[indx[j]+k+1]+P->stack[type][type_2]+
1013                        2*P->MLintern[1])
1014               break;
1015         }
1016       }
1017       sector[++s].i = i;
1018       sector[s].j   = k;
1019       sector[s].ml  = ml;
1020       sector[++s].i = k+1;
1021       sector[s].j   = j;
1022       sector[s].ml  = ml;
1023
1024       if (k>j-2-TURN) nrerror("backtrack failed in fML");
1025       continue;
1026     }
1027
1028   repeat1:
1029
1030     /*----- begin of "repeat:" -----*/
1031     ij = indx[j]+i;
1032     if (canonical)  cij = c[ij];
1033
1034     type = ptype[ij];
1035
1036     bonus = 0;
1037     if (struct_constrained) {
1038       if ((BP[i]==j)||(BP[i]==-1)||(BP[i]==-2)) bonus -= BONUS;
1039       if ((BP[j]==-1)||(BP[j]==-3)) bonus -= BONUS;
1040     }
1041     if (noLonelyPairs)
1042       if (cij == c[ij]){
1043         /* (i.j) closes canonical structures, thus
1044            (i+1.j-1) must be a pair                */
1045         type_2 = ptype[indx[j-1]+i+1]; type_2 = rtype[type_2];
1046         cij -= P->stack[type][type_2] + bonus;
1047         base_pair2[++b].i = i+1;
1048         base_pair2[b].j   = j-1;
1049         i++; j--;
1050         canonical=0;
1051         goto repeat1;
1052       }
1053     canonical = 1;
1054
1055
1056     no_close = (((type==3)||(type==4))&&no_closingGU&&(bonus==0));
1057     if (no_close) {
1058       if (cij == FORBIDDEN) continue;
1059     } else
1060       if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1, P)+bonus)
1061         continue;
1062
1063     for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
1064       minq = j-i+p-MAXLOOP-2;
1065       if (minq<p+1+TURN) minq = p+1+TURN;
1066       for (q = j-1; q >= minq; q--) {
1067
1068         type_2 = ptype[indx[q]+p];
1069         if (type_2==0) continue;
1070         type_2 = rtype[type_2];
1071         if (no_closingGU)
1072           if (no_close||(type_2==3)||(type_2==4))
1073             if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
1074
1075         /* energy = oldLoopEnergy(i, j, p, q, type, type_2); */
1076         energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
1077                             S1[i+1], S1[j-1], S1[p-1], S1[q+1], P);
1078
1079         new = energy+c[indx[q]+p]+bonus;
1080         traced = (cij == new);
1081         if (traced) {
1082           base_pair2[++b].i = p;
1083           base_pair2[b].j   = q;
1084           i = p, j = q;
1085           goto repeat1;
1086         }
1087       }
1088     }
1089
1090     /* end of repeat: --------------------------------------------------*/
1091
1092     /* (i.j) must close a multi-loop */
1093     tt = rtype[type];
1094     i1 = i+1; j1 = j-1;
1095
1096     if(with_gquad){
1097       /*
1098         The case that is handled here actually resembles something like
1099         an interior loop where the enclosing base pair is of regular
1100         kind and the enclosed pair is not a canonical one but a g-quadruplex
1101         that should then be decomposed further...
1102       */
1103       if(backtrack_GQuad_IntLoop(cij, i, j, type, S, ggg, indx, &p, &q, P)){
1104         i = p; j = q;
1105         goto repeat_gquad;
1106       }
1107     }
1108
1109     sector[s+1].ml  = sector[s+2].ml = 1;
1110
1111     switch(dangle_model){
1112       case 0:   en = cij - E_MLstem(tt, -1, -1, P) - P->MLclosing - bonus;
1113                 for(k = i+2+TURN; k < j-2-TURN; k++){
1114                   if(en == fML[indx[k]+i+1] + fML[indx[j-1]+k+1])
1115                     break;
1116                 }
1117                 break;
1118
1119       case 2:   en = cij - E_MLstem(tt, S1[j-1], S1[i+1], P) - P->MLclosing - bonus;
1120                 for(k = i+2+TURN; k < j-2-TURN; k++){
1121                     if(en == fML[indx[k]+i+1] + fML[indx[j-1]+k+1])
1122                       break;
1123                 }
1124                 break;
1125
1126       default:  for(k = i+2+TURN; k < j-2-TURN; k++){
1127                   en = cij - P->MLclosing - bonus;
1128                   if(en == fML[indx[k]+i+1] + fML[indx[j-1]+k+1] + E_MLstem(tt, -1, -1, P)){
1129                     break;
1130                   }
1131                   else if(en == fML[indx[k]+i+2] + fML[indx[j-1]+k+1] + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase){
1132                     i1 = i+2;
1133                     break;
1134                   }
1135                   else if(en == fML[indx[k]+i+1] + fML[indx[j-2]+k+1] + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase){
1136                     j1 = j-2;
1137                     break;
1138                   }
1139                   else if(en == fML[indx[k]+i+2] + fML[indx[j-2]+k+1] + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase){
1140                     i1 = i+2;
1141                     j1 = j-2;
1142                     break;
1143                   }
1144                   /* coaxial stacking of (i.j) with (i+1.k) or (k.j-1) */
1145                   /* use MLintern[1] since coax stacked pairs don't get TerminalAU */
1146                   if(dangle_model == 3){
1147                     type_2 = rtype[ptype[indx[k]+i+1]];
1148                     if (type_2) {
1149                       en = c[indx[k]+i+1]+P->stack[type][type_2]+fML[indx[j-1]+k+1];
1150                       if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1151                         ml = 2;
1152                         sector[s+1].ml  = 2;
1153                         traced = 1;
1154                         break;
1155                       }
1156                     }
1157                     type_2 = rtype[ptype[indx[j-1]+k+1]];
1158                     if (type_2) {
1159                       en = c[indx[j-1]+k+1]+P->stack[type][type_2]+fML[indx[k]+i+1];
1160                       if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1161                         sector[s+2].ml = 2;
1162                         traced = 1;
1163                         break;
1164                       }
1165                     }
1166                   }
1167                 }
1168                 break;
1169     }
1170
1171     if (k<=j-3-TURN) { /* found the decomposition */
1172       sector[++s].i = i1;
1173       sector[s].j   = k;
1174       sector[++s].i = k+1;
1175       sector[s].j   = j1;
1176     } else {
1177 #if 0
1178       /* Y shaped ML loops fon't work yet */
1179       if (dangle_model==3) {
1180         d5 = P->dangle5[tt][S1[j-1]];
1181         d3 = P->dangle3[tt][S1[i+1]];
1182         /* (i,j) must close a Y shaped ML loop with coax stacking */
1183         if (cij ==  fML[indx[j-2]+i+2] + mm + d3 + d5 + P->MLbase + P->MLbase) {
1184           i1 = i+2;
1185           j1 = j-2;
1186         } else if (cij ==  fML[indx[j-2]+i+1] + mm + d5 + P->MLbase)
1187           j1 = j-2;
1188         else if (cij ==  fML[indx[j-1]+i+2] + mm + d3 + P->MLbase)
1189           i1 = i+2;
1190         else /* last chance */
1191           if (cij != fML[indx[j-1]+i+1] + mm + P->MLbase)
1192             fprintf(stderr,  "backtracking failed in repeat");
1193         /* if we arrive here we can express cij via fML[i1,j1]+dangles */
1194         sector[++s].i = i1;
1195         sector[s].j   = j1;
1196       }
1197       else
1198 #endif
1199         nrerror("backtracking failed in repeat");
1200     }
1201
1202     continue; /* this is a workarround to not accidentally proceed in the following block */
1203
1204   repeat_gquad:
1205     /*
1206       now we do some fancy stuff to backtrace the stacksize and linker lengths
1207       of the g-quadruplex that should reside within position i,j
1208     */
1209     {
1210       int l[3], L, a;
1211       L = -1;
1212
1213       get_gquad_pattern_mfe(S, i, j, P, &L, l);
1214       if(L != -1){
1215         /* fill the G's of the quadruplex into base_pair2 */
1216         for(a=0;a<L;a++){
1217           base_pair2[++b].i = i+a;
1218           base_pair2[b].j   = i+a;
1219           base_pair2[++b].i = i+L+l[0]+a;
1220           base_pair2[b].j   = i+L+l[0]+a;
1221           base_pair2[++b].i = i+L+l[0]+L+l[1]+a;
1222           base_pair2[b].j   = i+L+l[0]+L+l[1]+a;
1223           base_pair2[++b].i = i+L+l[0]+L+l[1]+L+l[2]+a;
1224           base_pair2[b].j   = i+L+l[0]+L+l[1]+L+l[2]+a;
1225         }
1226         goto repeat_gquad_exit;
1227       }
1228       nrerror("backtracking failed in repeat_gquad");
1229     }
1230   repeat_gquad_exit:
1231     asm("nop");
1232
1233   } /* end of infinite while loop */
1234
1235   base_pair2[0].i = b;    /* save the total number of base pairs */
1236 }
1237
1238 PUBLIC char *backtrack_fold_from_pair(char *sequence, int i, int j) {
1239   char *structure;
1240   sector[1].i  = i;
1241   sector[1].j  = j;
1242   sector[1].ml = 2;
1243   base_pair2[0].i=0;
1244   S   = encode_sequence(sequence, 0);
1245   S1  = encode_sequence(sequence, 1);
1246   backtrack(sequence, 1);
1247   structure = (char *) space((strlen(sequence)+1)*sizeof(char));
1248   parenthesis_structure(structure, base_pair2, strlen(sequence));
1249   free(S);free(S1);
1250   return structure;
1251 }
1252
1253 /*---------------------------------------------------------------------------*/
1254
1255 PUBLIC void letter_structure(char *structure, bondT *bp, int length){
1256   int   n, k, x, y;
1257   char  alpha[] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz";
1258
1259   for (n = 0; n < length; structure[n++] = ' ');
1260   structure[length] = '\0';
1261
1262   for (n = 0, k = 1; k <= bp[0].i; k++) {
1263     y = bp[k].j;
1264     x = bp[k].i;
1265     if (x-1 > 0 && y+1 <= length) {
1266       if (structure[x-2] != ' ' && structure[y] == structure[x-2]) {
1267         structure[x-1] = structure[x-2];
1268         structure[y-1] = structure[x-1];
1269         continue;
1270       }
1271     }
1272     if (structure[x] != ' ' && structure[y-2] == structure[x]) {
1273       structure[x-1] = structure[x];
1274       structure[y-1] = structure[x-1];
1275       continue;
1276     }
1277     n++;
1278     structure[x-1] = alpha[n-1];
1279     structure[y-1] = alpha[n-1];
1280   }
1281 }
1282
1283 /*---------------------------------------------------------------------------*/
1284
1285 PUBLIC void parenthesis_structure(char *structure, bondT *bp, int length){
1286   int n, k;
1287
1288   for (n = 0; n < length; structure[n++] = '.');
1289   structure[length] = '\0';
1290
1291   for (k = 1; k <= bp[0].i; k++){
1292
1293     if(bp[k].i == bp[k].j){ /* Gquad bonds are marked as bp[i].i == bp[i].j */
1294       structure[bp[k].i-1] = '+';
1295     } else { /* the following ones are regular base pairs */
1296       structure[bp[k].i-1] = '(';
1297       structure[bp[k].j-1] = ')';
1298     }
1299   }
1300 }
1301
1302 PUBLIC void parenthesis_zuker(char *structure, bondT *bp, int length){
1303   int k, i, j, temp;
1304
1305   for (k = 0; k < length; structure[k++] = '.');
1306   structure[length] = '\0';
1307
1308   for (k = 1; k <= bp[0].i; k++) {
1309     i=bp[k].i;
1310     j=bp[k].j;
1311     if (i>length) i-=length;
1312     if (j>length) j-=length;
1313     if (i>j) {
1314       temp=i; i=j; j=temp;
1315     }
1316     if(i == j){ /* Gquad bonds are marked as bp[i].i == bp[i].j */
1317       structure[i-1] = '+';
1318     } else { /* the following ones are regular base pairs */
1319       structure[i-1] = '(';
1320       structure[j-1] = ')';
1321     }
1322   }
1323 }
1324
1325
1326 /*---------------------------------------------------------------------------*/
1327
1328 PUBLIC void update_fold_params(void){
1329   update_fold_params_par(NULL);
1330 }
1331
1332 PUBLIC void update_fold_params_par(paramT *parameters){
1333   if(P) free(P);
1334   if(parameters){
1335     P = get_parameter_copy(parameters);
1336   } else {
1337     model_detailsT md;
1338     set_model_details(&md);
1339     P = get_scaled_parameters(temperature, md);
1340   }
1341   make_pair_matrix();
1342   if (init_length < 0) init_length=0;
1343 }
1344
1345 /*---------------------------------------------------------------------------*/
1346 PUBLIC float energy_of_structure(const char *string, const char *structure, int verbosity_level){
1347   return energy_of_struct_par(string, structure, NULL, verbosity_level);
1348 }
1349
1350 PUBLIC float energy_of_struct_par(const char *string,
1351                                   const char *structure,
1352                                   paramT *parameters,
1353                                   int verbosity_level){
1354   int   energy;
1355   short *ss, *ss1;
1356
1357   update_fold_params_par(parameters);
1358
1359   if (strlen(structure)!=strlen(string))
1360     nrerror("energy_of_struct: string and structure have unequal length");
1361
1362   /* save the S and S1 pointers in case they were already in use */
1363   ss = S; ss1 = S1;
1364   S   = encode_sequence(string, 0);
1365   S1  = encode_sequence(string, 1);
1366
1367   pair_table = make_pair_table(structure);
1368
1369   energy = energy_of_structure_pt(string, pair_table, S, S1, verbosity_level);
1370
1371   free(pair_table);
1372   free(S); free(S1);
1373   S=ss; S1=ss1;
1374   return  (float) energy/100.;
1375 }
1376
1377 /*  returns a correction term that may be added to the energy retrieved
1378     from energy_of_struct_par() to correct misinterpreted loops. This
1379     correction is necessary since energy_of_struct_par() will forget
1380     about the existance of gquadruplexes and just treat them as unpaired
1381     regions.
1382
1383     recursive variant
1384 */
1385 PRIVATE int en_corr_of_loop_gquad(int i,
1386                                   int j,
1387                                   const char *string,
1388                                   const char *structure,
1389                                   short *pt,
1390                                   int *loop_idx,
1391                                   const short *s1){
1392
1393   int pos, energy, p, q, r, s, u, type, type2;
1394   int L, l[3];
1395
1396   energy = 0;
1397   q = i;
1398   while((pos = parse_gquad(structure + q-1, &L, l)) > 0){
1399     q += pos-1;
1400     p = q - 4*L - l[0] - l[1] - l[2] + 1;
1401     if(q > j) break;
1402     /* we've found the first g-quadruplex at position [p,q] */
1403     energy += E_gquad(L, l, P);
1404     /* check if it's enclosed in a base pair */
1405     if(loop_idx[p] == 0){ q++; continue; /* g-quad in exterior loop */}
1406     else{
1407       energy += E_MLstem(0, -1, -1, P);
1408       /*  find its enclosing pair */
1409       int num_elem, num_g, elem_i, elem_j, up_mis;
1410       num_elem  = 0;
1411       num_g     = 1;
1412       r         = p - 1;
1413       up_mis    = q - p + 1;
1414
1415       /* seek for first pairing base located 5' of the g-quad */
1416       for(r = p - 1; !pt[r] && (r >= i); r--);
1417       if(r < i) nrerror("this should not happen");
1418
1419       if(r < pt[r]){ /* found the enclosing pair */
1420         s = pt[r];
1421       } else {
1422         num_elem++;
1423         elem_i = pt[r];
1424         elem_j = r;
1425         r = pt[r]-1 ;
1426         /* seek for next pairing base 5' of r */
1427         for(; !pt[r] && (r >= i); r--);
1428         if(r < i) nrerror("so nich");
1429         if(r < pt[r]){ /* found the enclosing pair */
1430           s = pt[r];
1431         } else {
1432           /* hop over stems and unpaired nucleotides */
1433           while((r > pt[r]) && (r >= i)){
1434             if(pt[r]){ r = pt[r]; num_elem++;}
1435             r--;
1436           }
1437           if(r < i) nrerror("so nich");
1438           s = pt[r]; /* found the enclosing pair */
1439         }
1440       }
1441       /* now we have the enclosing pair (r,s) */
1442
1443       u = q+1;
1444       /* we know everything about the 5' part of this loop so check the 3' part */
1445       while(u<s){
1446         if(structure[u-1] == '.') u++;
1447         else if (structure[u-1] == '+'){ /* found another gquad */
1448           pos = parse_gquad(structure + u - 1, &L, l);
1449           if(pos > 0){
1450             energy += E_gquad(L, l, P) + E_MLstem(0, -1, -1, P);
1451             up_mis += pos;
1452             u += pos;
1453             num_g++;
1454           }
1455         } else { /* we must have found a stem */
1456           if(!(u < pt[u])) nrerror("wtf!");
1457           num_elem++; elem_i = u; elem_j = pt[u];
1458           energy += en_corr_of_loop_gquad(u, pt[u], string, structure, pt, loop_idx, s1);
1459           u = pt[u] + 1;
1460         }
1461       }
1462       if(u!=s) nrerror("what the hell");
1463       else{ /* we are done since we've found no other 3' structure element */
1464         switch(num_elem){
1465           /* g-quad was misinterpreted as hairpin closed by (r,s) */
1466           case 0:   /* if(num_g == 1)
1467                       if((p-r-1 == 0) || (s-q-1 == 0))
1468                         nrerror("too few unpaired bases");
1469                     */
1470                     type = pair[s1[r]][s1[s]];
1471                     if(dangles == 2)
1472                       energy += P->mismatchI[type][s1[r+1]][s1[s-1]];
1473                     if(type > 2)
1474                       energy += P->TerminalAU;
1475                     energy += P->internal_loop[s - r - 1 - up_mis];
1476                     energy -= E_Hairpin(s - r - 1,
1477                                         type,
1478                                         s1[r + 1],
1479                                         s1[s - 1],
1480                                         string + r - 1,
1481                                         P);
1482                     break;
1483           /* g-quad was misinterpreted as interior loop closed by (r,s) with enclosed pair (elem_i, elem_j) */
1484           case 1:   type = pair[s1[r]][s1[s]];
1485                     type2 = pair[s1[elem_i]][s1[elem_j]];
1486                     energy += P->MLclosing
1487                               + E_MLstem(rtype[type], s1[s-1], s1[r+1], P)
1488                               + (elem_i - r - 1 + s - elem_j - 1 - up_mis) * P->MLbase
1489                               + E_MLstem(type2, s1[elem_i-1], s1[elem_j+1], P);
1490                     energy -= E_IntLoop(elem_i - r - 1,
1491                                         s - elem_j - 1,
1492                                         type,
1493                                         rtype[type2],
1494                                         s1[r + 1],
1495                                         s1[s - 1],
1496                                         s1[elem_i - 1],
1497                                         s1[elem_j + 1],
1498                                         P);
1499                     break;
1500           /* gquad was misinterpreted as unpaired nucleotides in a multiloop */
1501           default:  energy -= (up_mis) * P->MLbase;
1502                     break;
1503         }
1504       }
1505       q = s+1;
1506     }
1507   }
1508   return energy;
1509 }
1510
1511 PUBLIC float energy_of_gquad_structure( const char *string,
1512                                         const char *structure,
1513                                         int verbosity_level){
1514
1515   int   energy, gge, *loop_idx;
1516   short *ss, *ss1;
1517
1518 #ifdef _OPENMP
1519   if(P == NULL) update_fold_params();
1520 #else
1521   if((init_length<0)||(P==NULL)) update_fold_params();
1522 #endif
1523
1524   if (fabs(P->temperature - temperature)>1e-6) update_fold_params();
1525
1526   if (strlen(structure)!=strlen(string))
1527     nrerror("energy_of_struct: string and structure have unequal length");
1528
1529   /* save the S and S1 pointers in case they were already in use */
1530   ss = S; ss1 = S1;
1531   S   = encode_sequence(string, 0);
1532   S1  = encode_sequence(string, 1);
1533
1534   /* the pair_table looses every information about the gquad position
1535      thus we have to find add the energy contributions for each loop
1536      that contains a gquad by ourself, substract all miscalculated
1537      contributions, i.e. loops that actually contain a gquad, from
1538      energy_of_structure_pt()
1539   */
1540   pair_table  = make_pair_table(structure);
1541   energy      = energy_of_structure_pt(string, pair_table, S, S1, verbosity_level);
1542
1543   loop_idx    = make_loop_index_pt(pair_table);
1544   gge         = en_corr_of_loop_gquad(1, S[0], string, structure, pair_table, loop_idx, S1);
1545   energy     += gge;
1546
1547   free(pair_table);
1548   free(loop_idx);
1549   free(S); free(S1);
1550   S=ss; S1=ss1;
1551   return  (float) energy/100.;
1552 }
1553
1554 PUBLIC int energy_of_structure_pt(const char *string,
1555                                   short *ptable,
1556                                   short *s,
1557                                   short *s1,
1558                                   int verbosity_level){
1559   return energy_of_struct_pt_par(string, ptable, s, s1, NULL, verbosity_level);
1560 }
1561
1562 PUBLIC int energy_of_struct_pt_par( const char *string,
1563                                     short *ptable,
1564                                     short *s,
1565                                     short *s1,
1566                                     paramT *parameters,
1567                                     int verbosity_level){
1568   /* auxiliary function for kinfold,
1569      for most purposes call energy_of_struct instead */
1570
1571   int   i, length, energy;
1572   short *ss, *ss1;
1573
1574   update_fold_params_par(parameters);
1575
1576   pair_table = ptable;
1577   ss  = S;
1578   ss1 = S1;
1579   S = s;
1580   S1 = s1;
1581
1582   length = S[0];
1583 /*   energy =  backtrack_type=='M' ? ML_Energy(0, 0) : ML_Energy(0, 1); */
1584     energy =  backtrack_type=='M' ? energy_of_ml_pt(0, ptable) : energy_of_extLoop_pt(0, ptable);
1585   if (verbosity_level>0)
1586     printf("External loop                           : %5d\n", energy);
1587   for (i=1; i<=length; i++) {
1588     if (pair_table[i]==0) continue;
1589     energy += stack_energy(i, string, verbosity_level);
1590     i=pair_table[i];
1591   }
1592   for (i=1; !SAME_STRAND(i,length); i++) {
1593     if (!SAME_STRAND(i,pair_table[i])) {
1594       energy+=P->DuplexInit;
1595       break;
1596     }
1597   }
1598   S   = ss;
1599   S1  = ss1;
1600   return energy;
1601 }
1602
1603 PUBLIC float energy_of_circ_structure(const char *string,
1604                                       const char *structure,
1605                                       int verbosity_level){
1606   return energy_of_circ_struct_par(string, structure, NULL, verbosity_level);
1607 }
1608
1609 PUBLIC float energy_of_circ_struct_par( const char *string,
1610                                         const char *structure,
1611                                         paramT *parameters,
1612                                         int verbosity_level){
1613
1614   int   i, j, length, energy=0, en0, degree=0, type;
1615   short *ss, *ss1;
1616
1617   update_fold_params_par(parameters);
1618
1619   int dangle_model = P->model_details.dangles;
1620
1621   if (strlen(structure)!=strlen(string))
1622     nrerror("energy_of_struct: string and structure have unequal length");
1623
1624   /* save the S and S1 pointers in case they were already in use */
1625   ss = S; ss1 = S1;
1626   S   = encode_sequence(string, 0);
1627   S1  = encode_sequence(string, 1);
1628
1629   pair_table = make_pair_table(structure);
1630
1631   length = S[0];
1632
1633   for (i=1; i<=length; i++) {
1634     if (pair_table[i]==0) continue;
1635     degree++;
1636     energy += stack_energy(i, string, verbosity_level);
1637     i=pair_table[i];
1638   }
1639
1640   if (degree==0) return 0.;
1641   for (i=1; pair_table[i]==0; i++);
1642   j = pair_table[i];
1643   type=pair[S[j]][S[i]];
1644   if (type==0) type=7;
1645   if (degree==1) {
1646     char loopseq[10];
1647     int u, si1, sj1;
1648     for (i=1; pair_table[i]==0; i++);
1649     u = length-j + i-1;
1650     if (u<7) {
1651       strcpy(loopseq , string+j-1);
1652       strncat(loopseq, string, i);
1653     }
1654     si1 = (i==1)?S1[length] : S1[i-1];
1655     sj1 = (j==length)?S1[1] : S1[j+1];
1656     en0 = E_Hairpin(u, type, sj1, si1, loopseq, P);
1657   } else
1658     if (degree==2) {
1659       int p,q, u1,u2, si1, sq1, type_2;
1660       for (p=j+1; pair_table[p]==0; p++);
1661       q=pair_table[p];
1662       u1 = p-j-1;
1663       u2 = i-1 + length-q;
1664       type_2 = pair[S[q]][S[p]];
1665       if (type_2==0) type_2=7;
1666       si1 = (i==1)? S1[length] : S1[i-1];
1667       sq1 = (q==length)? S1[1] : S1[q+1];
1668       en0 = E_IntLoop(u1, u2, type, type_2,
1669                        S1[j+1], si1, S1[p-1], sq1,P);
1670     } else { /* degree > 2 */
1671       en0 = ML_Energy(0, 0) - P->MLintern[0];
1672       if (dangle_model) {
1673         int d5, d3;
1674         if (pair_table[1]) {
1675           j = pair_table[1];
1676           type = pair[S[1]][S[j]];
1677           if (dangle_model==2)
1678             en0 += P->dangle5[type][S1[length]];
1679           else { /* dangle_model==1 */
1680             if (pair_table[length]==0) {
1681               d5 = P->dangle5[type][S1[length]];
1682               if (pair_table[length-1]!=0) {
1683                 int tt;
1684                 tt = pair[S[pair_table[length-1]]][S[length-1]];
1685                 d3 = P->dangle3[tt][S1[length]];
1686                 if (d3<d5) d5 = 0;
1687                 else d5 -= d3;
1688               }
1689               en0 += d5;
1690             }
1691           }
1692         }
1693         if (pair_table[length]) {
1694           i = pair_table[length];
1695           type = pair[S[i]][S[length]];
1696           if (dangle_model==2)
1697             en0 += P->dangle3[type][S1[1]];
1698           else { /* dangle_model==1 */
1699             if (pair_table[1]==0) {
1700               d3 = P->dangle3[type][S1[1]];
1701               if (pair_table[2]) {
1702                 int tt;
1703                 tt = pair[S[2]][S[pair_table[2]]];
1704                 d5 = P->dangle5[tt][1];
1705                 if (d5<d3) d3=0;
1706                 else d3 -= d5;
1707               }
1708               en0 += d3;
1709             }
1710           }
1711         }
1712       }
1713     }
1714
1715   if (verbosity_level>0)
1716     printf("External loop                           : %5d\n", en0);
1717   energy += en0;
1718   /* fprintf(stderr, "ext loop degree %d tot %d\n", degree, energy); */
1719   free(S); free(S1);
1720   S=ss; S1=ss1;
1721   return  (float) energy/100.0;
1722 }
1723
1724 /*---------------------------------------------------------------------------*/
1725 PRIVATE int stack_energy(int i, const char *string, int verbosity_level)
1726 {
1727   /* calculate energy of substructure enclosed by (i,j) */
1728   int ee, energy = 0;
1729   int j, p, q, type;
1730
1731   j=pair_table[i];
1732   type = pair[S[i]][S[j]];
1733   if (type==0) {
1734     type=7;
1735     if (verbosity_level>=0)
1736       fprintf(stderr,"WARNING: bases %d and %d (%c%c) can't pair!\n", i, j,
1737               string[i-1],string[j-1]);
1738   }
1739
1740   p=i; q=j;
1741   while (p<q) { /* process all stacks and interior loops */
1742     int type_2;
1743     while (pair_table[++p]==0);
1744     while (pair_table[--q]==0);
1745     if ((pair_table[q]!=(short)p)||(p>q)) break;
1746     type_2 = pair[S[q]][S[p]];
1747     if (type_2==0) {
1748       type_2=7;
1749       if (verbosity_level>=0)
1750         fprintf(stderr,"WARNING: bases %d and %d (%c%c) can't pair!\n", p, q,
1751                 string[p-1],string[q-1]);
1752     }
1753     /* energy += LoopEnergy(i, j, p, q, type, type_2); */
1754     if ( SAME_STRAND(i,p) && SAME_STRAND(q,j) )
1755       ee = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
1756     else
1757       ee = energy_of_extLoop_pt(cut_in_loop(i), pair_table);
1758     if (verbosity_level>0)
1759       printf("Interior loop (%3d,%3d) %c%c; (%3d,%3d) %c%c: %5d\n",
1760              i,j,string[i-1],string[j-1],p,q,string[p-1],string[q-1], ee);
1761     energy += ee;
1762     i=p; j=q; type = rtype[type_2];
1763   } /* end while */
1764
1765   /* p,q don't pair must have found hairpin or multiloop */
1766
1767   if (p>q) {                       /* hair pin */
1768     if (SAME_STRAND(i,j))
1769       ee = E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+i-1, P);
1770     else
1771       ee = energy_of_extLoop_pt(cut_in_loop(i), pair_table);
1772     energy += ee;
1773     if (verbosity_level>0)
1774       printf("Hairpin  loop (%3d,%3d) %c%c              : %5d\n",
1775              i, j, string[i-1],string[j-1], ee);
1776
1777     return energy;
1778   }
1779
1780   /* (i,j) is exterior pair of multiloop */
1781   while (p<j) {
1782     /* add up the contributions of the substructures of the ML */
1783     energy += stack_energy(p, string, verbosity_level);
1784     p = pair_table[p];
1785     /* search for next base pair in multiloop */
1786     while (pair_table[++p]==0);
1787   }
1788   {
1789     int ii;
1790     ii = cut_in_loop(i);
1791     ee = (ii==0) ? energy_of_ml_pt(i, pair_table) : energy_of_extLoop_pt(ii, pair_table);
1792   }
1793   energy += ee;
1794   if (verbosity_level>0)
1795     printf("Multi    loop (%3d,%3d) %c%c              : %5d\n",
1796            i,j,string[i-1],string[j-1],ee);
1797
1798   return energy;
1799 }
1800
1801 /*---------------------------------------------------------------------------*/
1802
1803
1804
1805 /**
1806 *** Calculate the energy contribution of
1807 *** stabilizing dangling-ends/mismatches
1808 *** for all stems branching off the exterior
1809 *** loop
1810 **/
1811 PRIVATE int energy_of_extLoop_pt(int i, short *pair_table) {
1812   int energy, mm5, mm3;
1813   int p, q, q_prev;
1814   int length = (int)pair_table[0];
1815
1816   /* helper variables for dangles == 1 case */
1817   int E3_available;  /* energy of 5' part where 5' mismatch is available for current stem */
1818   int E3_occupied;   /* energy of 5' part where 5' mismatch is unavailable for current stem */
1819
1820   int dangle_model = P->model_details.dangles;
1821
1822   /* initialize vars */
1823   energy      = 0;
1824   p           = (i==0) ? 1 : i;
1825   q_prev      = -1;
1826
1827   if(dangle_model%2 == 1){
1828     E3_available = INF;
1829     E3_occupied  = 0;
1830   }
1831
1832   /* seek to opening base of first stem */
1833   while(p <= length && !pair_table[p]) p++;
1834
1835   while(p < length){
1836     int tt;
1837     /* p must have a pairing partner */
1838     q  = (int)pair_table[p];
1839     /* get type of base pair (p,q) */
1840     tt = pair[S[p]][S[q]];
1841     if(tt==0) tt=7;
1842
1843     switch(dangle_model){
1844       /* no dangles */
1845       case 0:   energy += E_ExtLoop(tt, -1, -1, P);
1846                 break;
1847       /* the beloved double dangles */
1848       case 2:   mm5 = ((SAME_STRAND(p-1,p)) && (p>1))       ? S1[p-1] : -1;
1849                 mm3 = ((SAME_STRAND(q,q+1)) && (q<length))  ? S1[q+1] : -1;
1850                 energy += E_ExtLoop(tt, mm5, mm3, P);
1851                 break;
1852
1853       default:  {
1854                   int tmp;
1855                   if(q_prev + 2 < p){
1856                     E3_available = MIN2(E3_available, E3_occupied);
1857                     E3_occupied  = E3_available;
1858                   }
1859                   mm5 = ((SAME_STRAND(p-1,p)) && (p>1) && !pair_table[p-1])       ? S1[p-1] : -1;
1860                   mm3 = ((SAME_STRAND(q,q+1)) && (q<length) && !pair_table[q+1])  ? S1[q+1] : -1;
1861                   tmp = MIN2(
1862                                                 E3_occupied  + E_ExtLoop(tt, -1, mm3, P),
1863                                                 E3_available + E_ExtLoop(tt, mm5, mm3, P)
1864                                               );
1865                   E3_available =       MIN2(
1866                                                 E3_occupied  + E_ExtLoop(tt, -1, -1, P),
1867                                                 E3_available + E_ExtLoop(tt, mm5, -1, P)
1868                                               );
1869                   E3_occupied = tmp;
1870                 }
1871                 break;
1872
1873     } /* end switch dangle_model */
1874     /* seek to the next stem */
1875     p = q + 1;
1876     q_prev = q;
1877     while (p <= length && !pair_table[p]) p++;
1878     if(p==i) break; /* cut was in loop */
1879   }
1880
1881   if(dangle_model%2 == 1)
1882     energy = MIN2(E3_occupied, E3_available);
1883
1884   return energy;
1885 }
1886
1887 /**
1888 *** i is the 5'-base of the closing pair
1889 ***
1890 *** since each helix can coaxially stack with at most one of its
1891 *** neighbors we need an auxiliarry variable  cx_energy
1892 *** which contains the best energy given that the last two pairs stack.
1893 *** energy  holds the best energy given the previous two pairs do not
1894 *** stack (i.e. the two current helices may stack)
1895 *** We don't allow the last helix to stack with the first, thus we have to
1896 *** walk around the Loop twice with two starting points and take the minimum
1897 ***/
1898 PRIVATE int energy_of_ml_pt(int i, short *pt){
1899
1900   int energy, cx_energy, tmp, tmp2, best_energy=INF;
1901   int i1, j, p, q, q_prev, q_prev2, u, x, type, count, mm5, mm3, tt, ld5, new_cx, dang5, dang3, dang;
1902   int mlintern[NBPAIRS+1];
1903
1904   /* helper variables for dangles == 1|5 case */
1905   int E_mm5_available;  /* energy of 5' part where 5' mismatch of current stem is available */
1906   int E_mm5_occupied;   /* energy of 5' part where 5' mismatch of current stem is unavailable */
1907   int E2_mm5_available; /* energy of 5' part where 5' mismatch of current stem is available with possible 3' dangle for enclosing pair (i,j) */
1908   int E2_mm5_occupied;  /* energy of 5' part where 5' mismatch of current stem is unavailable with possible 3' dangle for enclosing pair (i,j) */
1909   int dangle_model = P->model_details.dangles;
1910
1911   if(i >= pt[i])
1912     nrerror("energy_of_ml_pt: i is not 5' base of a closing pair!");
1913
1914   j = (int)pt[i];
1915
1916   /* init the variables */
1917   energy      = 0;
1918   p           = i+1;
1919   q_prev      = i-1;
1920   q_prev2     = i;
1921
1922   for (x = 0; x <= NBPAIRS; x++) mlintern[x] = P->MLintern[x];
1923
1924   /* seek to opening base of first stem */
1925   while(p <= j && !pair_table[p]) p++;
1926   u = p - i - 1;
1927
1928   switch(dangle_model){
1929     case 0:   while(p < j){
1930                 /* p must have a pairing partner */
1931                 q  = (int)pair_table[p];
1932                 /* get type of base pair (p,q) */
1933                 tt = pair[S[p]][S[q]];
1934                 if(tt==0) tt=7;
1935                 energy += E_MLstem(tt, -1, -1, P);
1936                 /* seek to the next stem */
1937                 p = q + 1;
1938                 q_prev = q_prev2 = q;
1939                 while (p <= j && !pair_table[p]) p++;
1940                 u += p - q - 1; /* add unpaired nucleotides */
1941               }
1942               /* now lets get the energy of the enclosing stem */
1943               type = pair[S[j]][S[i]]; if (type==0) type=7;
1944               energy += E_MLstem(type, -1, -1, P);
1945               break;
1946
1947     case 2:   while(p < j){
1948                 /* p must have a pairing partner */
1949                 q  = (int)pair_table[p];
1950                 /* get type of base pair (p,q) */
1951                 tt = pair[S[p]][S[q]];
1952                 if(tt==0) tt=7;
1953                 mm5 = (SAME_STRAND(p-1,p))  ? S1[p-1] : -1;
1954                 mm3 = (SAME_STRAND(q,q+1))  ? S1[q+1] : -1;
1955                 energy += E_MLstem(tt, mm5, mm3, P);
1956                 /* seek to the next stem */
1957                 p = q + 1;
1958                 q_prev = q_prev2 = q;
1959                 while (p <= j && !pair_table[p]) p++;
1960                 u += p - q - 1; /* add unpaired nucleotides */
1961               }
1962               type = pair[S[j]][S[i]]; if (type==0) type=7;
1963               mm5 = ((SAME_STRAND(j-1,j)) && !pair_table[j-1])  ? S1[j-1] : -1;
1964               mm3 = ((SAME_STRAND(i,i+1)) && !pair_table[i+1])  ? S1[i+1] : -1;
1965               energy += E_MLstem(type, S1[j-1], S1[i+1], P);
1966               break;
1967
1968     case 3:   /* we treat helix stacking different */
1969               for (count=0; count<2; count++) { /* do it twice */
1970                 ld5 = 0; /* 5' dangle energy on prev pair (type) */
1971                 if ( i==0 ) {
1972                   j = (unsigned int)pair_table[0]+1;
1973                   type = 0;  /* no pair */
1974                 }
1975                 else {
1976                   j = (unsigned int)pair_table[i];
1977                   type = pair[S[j]][S[i]]; if (type==0) type=7;
1978                   /* prime the ld5 variable */
1979                   if (SAME_STRAND(j-1,j)) {
1980                     ld5 = P->dangle5[type][S1[j-1]];
1981                     if ((p=(unsigned int)pair_table[j-2]) && SAME_STRAND(j-2, j-1))
1982                     if (P->dangle3[pair[S[p]][S[j-2]]][S1[j-1]]<ld5) ld5 = 0;
1983                   }
1984                 }
1985                 i1=i; p = i+1; u=0;
1986                 energy = 0; cx_energy=INF;
1987                 do { /* walk around the multi-loop */
1988                   new_cx = INF;
1989
1990                   /* hop over unpaired positions */
1991                   while (p <= (unsigned int)pair_table[0] && pair_table[p]==0) p++;
1992
1993                   /* memorize number of unpaired positions */
1994                   u += p-i1-1;
1995                   /* get position of pairing partner */
1996                   if ( p == (unsigned int)pair_table[0]+1 ){
1997                     q = 0;tt = 0; /* virtual root pair */
1998                   } else {
1999                     q  = (unsigned int)pair_table[p];
2000                     /* get type of base pair P->q */
2001                     tt = pair[S[p]][S[q]]; if (tt==0) tt=7;
2002                   }
2003
2004                   energy += mlintern[tt];
2005                   cx_energy += mlintern[tt];
2006
2007                   dang5=dang3=0;
2008                   if ((SAME_STRAND(p-1,p))&&(p>1))
2009                     dang5=P->dangle5[tt][S1[p-1]];      /* 5'dangle of pq pair */
2010                   if ((SAME_STRAND(i1,i1+1))&&(i1<(unsigned int)S[0]))
2011                     dang3 = P->dangle3[type][S1[i1+1]]; /* 3'dangle of previous pair */
2012
2013                   switch (p-i1-1) {
2014                     case 0:   /* adjacent helices */
2015                               if (i1!=0){
2016                                 if (SAME_STRAND(i1,p)) {
2017                                   new_cx = energy + P->stack[rtype[type]][rtype[tt]];
2018                                   /* subtract 5'dangle and TerminalAU penalty */
2019                                   new_cx += -ld5 - mlintern[tt]-mlintern[type]+2*mlintern[1];
2020                                 }
2021                                 ld5=0;
2022                                 energy = MIN2(energy, cx_energy);
2023                               }
2024                               break;
2025                     case 1:   /* 1 unpaired base between helices */
2026                               dang = MIN2(dang3, dang5);
2027                               energy = energy +dang; ld5 = dang - dang3;
2028                               /* may be problem here: Suppose
2029                                 cx_energy>energy, cx_energy+dang5<energy
2030                                 and the following helices are also stacked (i.e.
2031                                 we'll subtract the dang5 again */
2032                               if (cx_energy+dang5 < energy) {
2033                                 energy = cx_energy+dang5;
2034                                 ld5 = dang5;
2035                               }
2036                               new_cx = INF;  /* no coax stacking with mismatch for now */
2037                               break;
2038                     default:  /* many unpaired base between helices */
2039                               energy += dang5 +dang3;
2040                               energy = MIN2(energy, cx_energy + dang5);
2041                               new_cx = INF;  /* no coax stacking possible */
2042                               ld5 = dang5;
2043                               break;
2044                   }
2045                   type = tt;
2046                   cx_energy = new_cx;
2047                   i1 = q; p=q+1;
2048                 } while (q!=i);
2049                 best_energy = MIN2(energy, best_energy); /* don't use cx_energy here */
2050                 /* fprintf(stderr, "%6.2d\t", energy); */
2051                 /* skip a helix and start again */
2052                 while (pair_table[p]==0) p++;
2053                 if (i == (unsigned int)pair_table[p]) break;
2054                 i = (unsigned int)pair_table[p];
2055               } /* end doing it twice */
2056               energy = best_energy;
2057               break;
2058
2059     default:  E_mm5_available = E2_mm5_available  = INF;
2060               E_mm5_occupied  = E2_mm5_occupied   = 0;
2061               while(p < j){
2062                 /* p must have a pairing partner */
2063                 q  = (int)pair_table[p];
2064                 /* get type of base pair (p,q) */
2065                 tt = pair[S[p]][S[q]];
2066                 if(tt==0) tt=7;
2067                 if(q_prev + 2 < p){
2068                   E_mm5_available = MIN2(E_mm5_available, E_mm5_occupied);
2069                   E_mm5_occupied  = E_mm5_available;
2070                 }
2071                 if(q_prev2 + 2 < p){
2072                   E2_mm5_available  = MIN2(E2_mm5_available, E2_mm5_occupied);
2073                   E2_mm5_occupied   = E2_mm5_available;
2074                 }
2075                 mm5 = ((SAME_STRAND(p-1,p)) && !pair_table[p-1])  ? S1[p-1] : -1;
2076                 mm3 = ((SAME_STRAND(q,q+1)) && !pair_table[q+1])  ? S1[q+1] : -1;
2077                 tmp =                   MIN2(
2078                                               E_mm5_occupied  + E_MLstem(tt, -1, mm3, P),
2079                                               E_mm5_available + E_MLstem(tt, mm5, mm3, P)
2080                                             );
2081                 tmp   =                 MIN2(tmp, E_mm5_available + E_MLstem(tt, -1, mm3, P));
2082                 tmp2  =                 MIN2(
2083                                               E_mm5_occupied  + E_MLstem(tt, -1, -1, P),
2084                                               E_mm5_available + E_MLstem(tt, mm5, -1, P)
2085                                             );
2086                 E_mm5_available =       MIN2(tmp2, E_mm5_available  + E_MLstem(tt, -1, -1, P));
2087                 E_mm5_occupied  = tmp;
2088
2089                 tmp =                  MIN2(
2090                                               E2_mm5_occupied  + E_MLstem(tt, -1, mm3, P),
2091                                               E2_mm5_available + E_MLstem(tt, mm5, mm3, P)
2092                                             );
2093                 tmp =                   MIN2(tmp, E2_mm5_available + E_MLstem(tt, -1, mm3, P));
2094                 tmp2 =                  MIN2(
2095                                               E2_mm5_occupied  + E_MLstem(tt, -1, -1, P),
2096                                               E2_mm5_available + E_MLstem(tt, mm5, -1, P)
2097                                             );
2098                 E2_mm5_available =      MIN2(tmp2, E2_mm5_available + E_MLstem(tt, -1, -1, P));
2099                 E2_mm5_occupied = tmp;
2100                 /* printf("(%d,%d): \n E_o = %d, E_a = %d, E2_o = %d, E2_a = %d\n", p, q, E_mm5_occupied,E_mm5_available,E2_mm5_occupied,E2_mm5_available); */
2101                 /* seek to the next stem */
2102                 p = q + 1;
2103                 q_prev = q_prev2 = q;
2104                 while (p <= j && !pair_table[p]) p++;
2105                 u += p - q - 1; /* add unpaired nucleotides */
2106               }
2107               /* now lets see how we get the minimum including the enclosing stem */
2108               type = pair[S[j]][S[i]]; if (type==0) type=7;
2109               mm5 = ((SAME_STRAND(j-1,j)) && !pair_table[j-1])  ? S1[j-1] : -1;
2110               mm3 = ((SAME_STRAND(i,i+1)) && !pair_table[i+1])  ? S1[i+1] : -1;
2111               if(q_prev + 2 < p){
2112                 E_mm5_available = MIN2(E_mm5_available, E_mm5_occupied);
2113                 E_mm5_occupied  = E_mm5_available;
2114               }
2115               if(q_prev2 + 2 < p){
2116                 E2_mm5_available  = MIN2(E2_mm5_available, E2_mm5_occupied);
2117                 E2_mm5_occupied   = E2_mm5_available;
2118               }
2119               energy = MIN2(E_mm5_occupied  + E_MLstem(type, -1, -1, P),
2120                             E_mm5_available + E_MLstem(type, mm5, -1, P)
2121                           );
2122               energy = MIN2(energy, E_mm5_available   + E_MLstem(type, -1, -1, P));
2123               energy = MIN2(energy, E2_mm5_occupied   + E_MLstem(type, -1, mm3, P));
2124               energy = MIN2(energy, E2_mm5_occupied   + E_MLstem(type, -1, -1, P));
2125               energy = MIN2(energy, E2_mm5_available  + E_MLstem(type, mm5, mm3, P));
2126               energy = MIN2(energy, E2_mm5_available  + E_MLstem(type, -1, mm3, P));
2127               energy = MIN2(energy, E2_mm5_available  + E_MLstem(type, mm5, -1, P));
2128               energy = MIN2(energy, E2_mm5_available  + E_MLstem(type, -1, -1, P));
2129               break;
2130   }/* end switch dangle_model */
2131
2132   energy += P->MLclosing;
2133   /* logarithmic ML loop energy if logML */
2134   if(logML && (u>6))
2135     energy += 6*P->MLbase+(int)(P->lxc*log((double)u/6.));
2136   else
2137     energy += (u*P->MLbase);
2138
2139   return energy;
2140 }
2141
2142 /*---------------------------------------------------------------------------*/
2143
2144 PUBLIC int loop_energy(short * ptable, short *s, short *s1, int i) {
2145   /* compute energy of a single loop closed by base pair (i,j) */
2146   int j, type, p,q, energy;
2147   short *Sold, *S1old, *ptold;
2148
2149   ptold=pair_table;   Sold = S;   S1old = S1;
2150   pair_table = ptable;   S = s;   S1 = s1;
2151
2152   if (i==0) { /* evaluate exterior loop */
2153     energy = energy_of_extLoop_pt(0,pair_table);
2154     pair_table=ptold; S=Sold; S1=S1old;
2155     return energy;
2156   }
2157   j = pair_table[i];
2158   if (j<i) nrerror("i is unpaired in loop_energy()");
2159   type = pair[S[i]][S[j]];
2160   if (type==0) {
2161     type=7;
2162     if (eos_debug>=0)
2163       fprintf(stderr,"WARNING: bases %d and %d (%c%c) can't pair!\n", i, j,
2164               Law_and_Order[S[i]],Law_and_Order[S[j]]);
2165   }
2166   p=i; q=j;
2167
2168
2169   while (pair_table[++p]==0);
2170   while (pair_table[--q]==0);
2171   if (p>q) { /* Hairpin */
2172     char loopseq[8] = "";
2173     if (SAME_STRAND(i,j)) {
2174       if (j-i-1<7) {
2175         int u;
2176         for (u=0; i+u<=j; u++) loopseq[u] = Law_and_Order[S[i+u]];
2177         loopseq[u] = '\0';
2178       }
2179       energy = E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], loopseq, P);
2180     } else {
2181       energy = energy_of_extLoop_pt(cut_in_loop(i), pair_table);
2182     }
2183   }
2184   else if (pair_table[q]!=(short)p) { /* multi-loop */
2185     int ii;
2186     ii = cut_in_loop(i);
2187     energy = (ii==0) ? energy_of_ml_pt(i, pair_table) : energy_of_extLoop_pt(ii, pair_table);
2188   }
2189   else { /* found interior loop */
2190     int type_2;
2191     type_2 = pair[S[q]][S[p]];
2192     if (type_2==0) {
2193       type_2=7;
2194       if (eos_debug>=0)
2195         fprintf(stderr,"WARNING: bases %d and %d (%c%c) can't pair!\n", p, q,
2196                 Law_and_Order[S[p]],Law_and_Order[S[q]]);
2197     }
2198     /* energy += LoopEnergy(i, j, p, q, type, type_2); */
2199     if ( SAME_STRAND(i,p) && SAME_STRAND(q,j) )
2200       energy = E_IntLoop(p-i-1, j-q-1, type, type_2,
2201                           S1[i+1], S1[j-1], S1[p-1], S1[q+1], P);
2202     else
2203       energy = energy_of_extLoop_pt(cut_in_loop(i), pair_table);
2204   }
2205
2206   pair_table=ptold; S=Sold; S1=S1old;
2207   return energy;
2208 }
2209
2210 /*---------------------------------------------------------------------------*/
2211
2212
2213 PUBLIC float energy_of_move(const char *string, const char *structure, int m1, int m2) {
2214   int   energy;
2215   short *ss, *ss1;
2216
2217 #ifdef _OPENMP
2218   if(P == NULL) update_fold_params();
2219 #else
2220   if((init_length<0)||(P==NULL)) update_fold_params();
2221 #endif
2222
2223   if (fabs(P->temperature - temperature)>1e-6) update_fold_params();
2224
2225   if (strlen(structure)!=strlen(string))
2226     nrerror("energy_of_struct: string and structure have unequal length");
2227
2228   /* save the S and S1 pointers in case they were already in use */
2229   ss = S; ss1 = S1;
2230   S   = encode_sequence(string, 0);
2231   S1  = encode_sequence(string, 1);
2232
2233   pair_table = make_pair_table(structure);
2234
2235   energy = energy_of_move_pt(pair_table, S, S1, m1, m2);
2236
2237   free(pair_table);
2238   free(S); free(S1);
2239   S=ss; S1=ss1;
2240   return  (float) energy/100.;
2241 }
2242
2243 /*---------------------------------------------------------------------------*/
2244
2245 PUBLIC int energy_of_move_pt(short *pt, short *s, short *s1, int m1, int m2) {
2246   /*compute change in energy given by move (m1,m2)*/
2247   int en_post, en_pre, i,j,k,l, len;
2248
2249   len = pt[0];
2250   k = (m1>0)?m1:-m1;
2251   l = (m2>0)?m2:-m2;
2252   /* first find the enclosing pair i<k<l<j */
2253   for (j=l+1; j<=len; j++) {
2254     if (pt[j]<=0) continue; /* unpaired */
2255     if (pt[j]<k) break;   /* found it */
2256     if (pt[j]>j) j=pt[j]; /* skip substructure */
2257     else {
2258       fprintf(stderr, "%d %d %d %d ", m1, m2, j, pt[j]);
2259       nrerror("illegal move or broken pair table in energy_of_move()");
2260     }
2261   }
2262   i = (j<=len) ? pt[j] : 0;
2263   en_pre = loop_energy(pt, s, s1, i);
2264   en_post = 0;
2265   if (m1<0) { /*it's a delete move */
2266     en_pre += loop_energy(pt, s, s1, k);
2267     pt[k]=0;
2268     pt[l]=0;
2269   } else { /* insert move */
2270     pt[k]=l;
2271     pt[l]=k;
2272     en_post += loop_energy(pt, s, s1, k);
2273   }
2274   en_post += loop_energy(pt, s, s1, i);
2275   /*  restore pair table */
2276   if (m1<0) {
2277     pt[k]=l;
2278     pt[l]=k;
2279   } else {
2280     pt[k]=0;
2281     pt[l]=0;
2282   }
2283   return (en_post - en_pre);
2284 }
2285
2286
2287
2288 PRIVATE int cut_in_loop(int i) {
2289   /* walk around the loop;  return j pos of pair after cut if
2290      cut_point in loop else 0 */
2291   int  p, j;
2292   p = j = pair_table[i];
2293   do {
2294     i  = pair_table[p];  p = i+1;
2295     while ( pair_table[p]==0 ) p++;
2296   } while (p!=j && SAME_STRAND(i,p));
2297   return SAME_STRAND(i,p) ? 0 : j;
2298 }
2299
2300 /*---------------------------------------------------------------------------*/
2301
2302 PRIVATE void make_ptypes(const short *S, const char *structure) {
2303   int n,i,j,k,l;
2304
2305   n=S[0];
2306   for (k=1; k<n-TURN; k++)
2307     for (l=1; l<=2; l++) {
2308       int type,ntype=0,otype=0;
2309       i=k; j = i+TURN+l; if (j>n) continue;
2310       type = pair[S[i]][S[j]];
2311       while ((i>=1)&&(j<=n)) {
2312         if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
2313         if (noLonelyPairs && (!otype) && (!ntype))
2314           type = 0; /* i.j can only form isolated pairs */
2315         ptype[indx[j]+i] = (char) type;
2316         otype =  type;
2317         type  = ntype;
2318         i--; j++;
2319       }
2320     }
2321
2322   if (struct_constrained && (structure != NULL))
2323     constrain_ptypes(structure, (unsigned int)n, ptype, BP, TURN, 0);
2324 }
2325
2326 PUBLIC void assign_plist_from_db(plist **pl, const char *struc, float pr){
2327   /* convert bracket string to plist */
2328   short *pt;
2329   int i, k = 0, size, n;
2330   plist *gpl, *ptr;
2331
2332   size  = strlen(struc);
2333   n     = 2;
2334
2335   pt  = make_pair_table(struc);
2336   *pl = (plist *)space(n*size*sizeof(plist));
2337   for(i = 1; i < size; i++){
2338     if(pt[i]>i){
2339       (*pl)[k].i      = i;
2340       (*pl)[k].j      = pt[i];
2341       (*pl)[k].p      = pr;
2342       (*pl)[k++].type = 0;
2343     }
2344   }
2345
2346   gpl = get_plist_gquad_from_db(struc, pr);
2347   for(ptr = gpl; ptr->i != 0; ptr++){
2348     if (k == n * size - 1){
2349       n *= 2;
2350       *pl = (plist *)xrealloc(*pl, n * size * sizeof(plist));
2351     }
2352     (*pl)[k].i      = ptr->i;
2353     (*pl)[k].j      = ptr->j;
2354     (*pl)[k].p       = ptr->p;
2355     (*pl)[k++].type = ptr->type;
2356   }
2357   free(gpl);
2358
2359   (*pl)[k].i      = 0;
2360   (*pl)[k].j      = 0;
2361   (*pl)[k].p      = 0.;
2362   (*pl)[k++].type = 0.;
2363   free(pt);
2364   *pl = (plist *)xrealloc(*pl, k * sizeof(plist));
2365 }
2366
2367
2368 /*###########################################*/
2369 /*# deprecated functions below              #*/
2370 /*###########################################*/
2371
2372 PUBLIC int HairpinE(int size, int type, int si1, int sj1, const char *string) {
2373   int energy;
2374
2375   energy = (size <= 30) ? P->hairpin[size] :
2376     P->hairpin[30]+(int)(P->lxc*log((size)/30.));
2377
2378   if (tetra_loop){
2379     if (size == 4) { /* check for tetraloop bonus */
2380       char tl[7]={0}, *ts;
2381       strncpy(tl, string, 6);
2382       if ((ts=strstr(P->Tetraloops, tl)))
2383         return (P->Tetraloop_E[(ts - P->Tetraloops)/7]);
2384     }
2385     if (size == 6) {
2386       char tl[9]={0}, *ts;
2387       strncpy(tl, string, 8);
2388       if ((ts=strstr(P->Hexaloops, tl)))
2389         return (energy = P->Hexaloop_E[(ts - P->Hexaloops)/9]);
2390     }
2391     if (size == 3) {
2392       char tl[6]={0,0,0,0,0,0}, *ts;
2393       strncpy(tl, string, 5);
2394       if ((ts=strstr(P->Triloops, tl))) {
2395         return (P->Triloop_E[(ts - P->Triloops)/6]);
2396       }
2397       if (type>2)  /* neither CG nor GC */
2398         energy += P->TerminalAU; /* penalty for closing AU GU pair IVOO??
2399                                     sind dass jetzt beaunuesse oder mahlnuesse (vorzeichen?)*/
2400       return energy;
2401     }
2402    }
2403    energy += P->mismatchH[type][si1][sj1];
2404
2405   return energy;
2406 }
2407
2408 /*---------------------------------------------------------------------------*/
2409
2410 PUBLIC int oldLoopEnergy(int i, int j, int p, int q, int type, int type_2) {
2411   /* compute energy of degree 2 loop (stack bulge or interior) */
2412   int n1, n2, m, energy;
2413   n1 = p-i-1;
2414   n2 = j-q-1;
2415
2416   if (n1>n2) { m=n1; n1=n2; n2=m; } /* so that n2>=n1 */
2417
2418   if (n2 == 0)
2419     energy = P->stack[type][type_2];   /* stack */
2420
2421   else if (n1==0) {                  /* bulge */
2422     energy = (n2<=MAXLOOP)?P->bulge[n2]:
2423       (P->bulge[30]+(int)(P->lxc*log(n2/30.)));
2424
2425 #if STACK_BULGE1
2426     if (n2==1) energy+=P->stack[type][type_2];
2427 #endif
2428   } else {                           /* interior loop */
2429
2430     if ((n1+n2==2)&&(james_rule))
2431       /* special case for loop size 2 */
2432       energy = P->int11[type][type_2][S1[i+1]][S1[j-1]];
2433     else {
2434       energy = (n1+n2<=MAXLOOP)?(P->internal_loop[n1+n2]):
2435         (P->internal_loop[30]+(int)(P->lxc*log((n1+n2)/30.)));
2436
2437 #if NEW_NINIO
2438       energy += MIN2(MAX_NINIO, (n2-n1)*P->ninio[2]);
2439 #else
2440       m       = MIN2(4, n1);
2441       energy += MIN2(MAX_NINIO,((n2-n1)*P->ninio[m]));
2442 #endif
2443       energy += P->mismatchI[type][S1[i+1]][S1[j-1]]+
2444         P->mismatchI[type_2][S1[q+1]][S1[p-1]];
2445     }
2446   }
2447   return energy;
2448 }
2449
2450 /*--------------------------------------------------------------------------*/
2451
2452 PUBLIC int LoopEnergy(int n1, int n2, int type, int type_2,
2453                       int si1, int sj1, int sp1, int sq1) {
2454   /* compute energy of degree 2 loop (stack bulge or interior) */
2455   int nl, ns, energy;
2456
2457   if (n1>n2) { nl=n1; ns=n2;}
2458   else {nl=n2; ns=n1;}
2459
2460   if (nl == 0)
2461     return P->stack[type][type_2];    /* stack */
2462
2463   if (ns==0) {                       /* bulge */
2464     energy = (nl<=MAXLOOP)?P->bulge[nl]:
2465       (P->bulge[30]+(int)(P->lxc*log(nl/30.)));
2466     if (nl==1) energy += P->stack[type][type_2];
2467     else {
2468       if (type>2) energy += P->TerminalAU;
2469       if (type_2>2) energy += P->TerminalAU;
2470     }
2471     return energy;
2472   }
2473   else {                             /* interior loop */
2474     if (ns==1) {
2475       if (nl==1)                     /* 1x1 loop */
2476         return P->int11[type][type_2][si1][sj1];
2477       if (nl==2) {                   /* 2x1 loop */
2478         if (n1==1)
2479           energy = P->int21[type][type_2][si1][sq1][sj1];
2480         else
2481           energy = P->int21[type_2][type][sq1][si1][sp1];
2482         return energy;
2483       }
2484         else {  /* 1xn loop */
2485         energy = (nl+1<=MAXLOOP)?(P->internal_loop[nl+1]):
2486         (P->internal_loop[30]+(int)(P->lxc*log((nl+1)/30.)));
2487         energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
2488         energy += P->mismatch1nI[type][si1][sj1]+
2489         P->mismatch1nI[type_2][sq1][sp1];
2490         return energy;
2491         }
2492     }
2493     else if (ns==2) {
2494       if(nl==2)      {   /* 2x2 loop */
2495         return P->int22[type][type_2][si1][sp1][sq1][sj1];}
2496       else if (nl==3)  { /* 2x3 loop */
2497         energy = P->internal_loop[5]+P->ninio[2];
2498         energy += P->mismatch23I[type][si1][sj1]+
2499           P->mismatch23I[type_2][sq1][sp1];
2500         return energy;
2501       }
2502
2503     }
2504     { /* generic interior loop (no else here!)*/
2505       energy = (n1+n2<=MAXLOOP)?(P->internal_loop[n1+n2]):
2506         (P->internal_loop[30]+(int)(P->lxc*log((n1+n2)/30.)));
2507
2508       energy += MIN2(MAX_NINIO, (nl-ns)*P->ninio[2]);
2509
2510       energy += P->mismatchI[type][si1][sj1]+
2511         P->mismatchI[type_2][sq1][sp1];
2512     }
2513   }
2514   return energy;
2515 }
2516
2517 PRIVATE int ML_Energy(int i, int is_extloop) {
2518   /* i is the 5'-base of the closing pair (or 0 for exterior loop)
2519      loop is scored as ML if extloop==0 else as exterior loop
2520
2521      since each helix can coaxially stack with at most one of its
2522      neighbors we need an auxiliarry variable  cx_energy
2523      which contains the best energy given that the last two pairs stack.
2524      energy  holds the best energy given the previous two pairs do not
2525      stack (i.e. the two current helices may stack)
2526      We don't allow the last helix to stack with the first, thus we have to
2527      walk around the Loop twice with two starting points and take the minimum
2528   */
2529
2530   int energy, cx_energy, best_energy=INF;
2531   int i1, j, p, q, u, x, type, count;
2532   int mlintern[NBPAIRS+1], mlclosing, mlbase;
2533   int dangle_model = P->model_details.dangles;
2534
2535   if (is_extloop) {
2536     for (x = 0; x <= NBPAIRS; x++)
2537       mlintern[x] = P->MLintern[x]-P->MLintern[1]; /* 0 or TerminalAU */
2538     mlclosing = mlbase = 0;
2539   } else {
2540     for (x = 0; x <= NBPAIRS; x++) mlintern[x] = P->MLintern[x];
2541     mlclosing = P->MLclosing; mlbase = P->MLbase;
2542   }
2543
2544   /*  as we do not only have dangling end but also mismatch contributions,
2545   **  we do this a bit different to previous implementations
2546   */
2547   if(is_extloop){
2548     energy = 0;
2549     i1  = i;
2550     p   = i+1;
2551
2552     int E_mm5_available, E_mm5_occupied;
2553     /* find out if we may have 5' mismatch for the next stem */
2554     while (p <= (int)pair_table[0] && pair_table[p]==0) p++;
2555     /* get position of pairing partner */
2556     if(p < (int)pair_table[0]){
2557         E_mm5_occupied  = (p - i - 1 > 0) ? INF : 0;
2558         E_mm5_available = (p - i - 1 > 0) ? 0 : INF;
2559     }
2560
2561     if(p < (int)pair_table[0])
2562       do{
2563         int tt;
2564         /* p must have a pairing partner */
2565         q  = (int)pair_table[p];
2566         /* get type of base pair (p,q) */
2567         tt = pair[S[p]][S[q]];
2568         if(tt==0) tt=7;
2569
2570         int mm5 = ((SAME_STRAND(p-1,p)) && (p>1)) ? S1[p-1]: -1;
2571         int mm3 = ((SAME_STRAND(q,q+1)) && (q<(unsigned int)pair_table[0])) ? S1[q+1]: -1;
2572
2573         switch(dangle_model){
2574           /* dangle_model == 0 */
2575           case 0: energy += E_ExtLoop(tt, -1, -1, P);
2576                   break;
2577           /* dangle_model == 1 */
2578           case 1: {
2579                     /* check for unpaired nucleotide 3' to the current stem */
2580                     int u3 = ((q < pair_table[0]) && (pair_table[q+1] == 0)) ? 1 : 0;
2581                     if(pair_table[p-1] != 0) mm5 = -1;
2582
2583                     if(!u3){
2584                       mm3 = -1;
2585                       E_mm5_occupied  = MIN2(
2586                                               E_mm5_occupied  + E_ExtLoop(tt, -1, -1, P),
2587                                               E_mm5_available + E_ExtLoop(tt, mm5, -1, P)
2588                                             );
2589                       E_mm5_available = E_mm5_occupied;
2590                     }
2591                     else{
2592                       E_mm5_occupied  = MIN2(
2593                                               E_mm5_occupied  + E_ExtLoop(tt, -1, mm3, P),
2594                                               E_mm5_available + E_ExtLoop(tt, mm5, mm3, P)
2595                                             );
2596                       E_mm5_available = MIN2(
2597                                               E_mm5_occupied  + E_ExtLoop(tt, -1, -1, P),
2598                                               E_mm5_available + E_ExtLoop(tt, mm5, -1, P)
2599                                             );
2600                     }
2601                   }
2602                   break;
2603
2604           /* the beloved case dangle_model == 2 */
2605           case 2: energy += E_ExtLoop(tt, mm5, mm3, P);
2606                   break;
2607
2608           /* dangle_model == 3 a.k.a. helix stacking */
2609           case 3: break;
2610
2611         } /* end switch dangle_model */
2612
2613         /* seek to the next stem */
2614         p = q + 1;
2615         while (p <= (int)pair_table[0] && pair_table[p]==0) p++;
2616         if(p == (int)pair_table[0] + 1){
2617           if(dangle_model == 1)
2618             energy = (p > q + 1) ? E_mm5_occupied : E_mm5_available;
2619           q = 0;
2620           break;
2621         }
2622       } while(q != i);
2623   }
2624   /* not exterior loop */
2625   else{
2626     for (count=0; count<2; count++) { /* do it twice */
2627       int ld5 = 0; /* 5' dangle energy on prev pair (type) */
2628       if ( i==0 ) {
2629         j = (unsigned int)pair_table[0]+1;
2630         type = 0;  /* no pair */
2631       }
2632       else {
2633         j = (unsigned int)pair_table[i];
2634         type = pair[S[j]][S[i]]; if (type==0) type=7;
2635
2636         if (dangle_model==3) { /* prime the ld5 variable */
2637           if (SAME_STRAND(j-1,j)) {
2638             ld5 = P->dangle5[type][S1[j-1]];
2639             if ((p=(unsigned int)pair_table[j-2]) && SAME_STRAND(j-2, j-1))
2640                 if (P->dangle3[pair[S[p]][S[j-2]]][S1[j-1]]<ld5) ld5 = 0;
2641           }
2642         }
2643       }
2644       i1=i; p = i+1; u=0;
2645       energy = 0; cx_energy=INF;
2646       do { /* walk around the multi-loop */
2647         int tt, new_cx = INF;
2648
2649         /* hop over unpaired positions */
2650         while (p <= (unsigned int)pair_table[0] && pair_table[p]==0) p++;
2651
2652         /* memorize number of unpaired positions */
2653         u += p-i1-1;
2654         /* get position of pairing partner */
2655         if ( p == (unsigned int)pair_table[0]+1 ){
2656           q = 0;tt = 0; /* virtual root pair */
2657         } else {
2658         q  = (unsigned int)pair_table[p];
2659           /* get type of base pair P->q */
2660         tt = pair[S[p]][S[q]]; if (tt==0) tt=7;
2661         }
2662
2663         energy += mlintern[tt];
2664         cx_energy += mlintern[tt];
2665
2666         if (dangle_model) {
2667           int dang5=0, dang3=0, dang;
2668           if ((SAME_STRAND(p-1,p))&&(p>1))
2669             dang5=P->dangle5[tt][S1[p-1]];      /* 5'dangle of pq pair */
2670           if ((SAME_STRAND(i1,i1+1))&&(i1<(unsigned int)S[0]))
2671             dang3 = P->dangle3[type][S1[i1+1]]; /* 3'dangle of previous pair */
2672
2673           switch (p-i1-1) {
2674           case 0: /* adjacent helices */
2675             if (dangle_model==2)
2676               energy += dang3+dang5;
2677             else if (dangle_model==3 && i1!=0) {
2678               if (SAME_STRAND(i1,p)) {
2679                 new_cx = energy + P->stack[rtype[type]][rtype[tt]];
2680                 /* subtract 5'dangle and TerminalAU penalty */
2681                 new_cx += -ld5 - mlintern[tt]-mlintern[type]+2*mlintern[1];
2682               }
2683               ld5=0;
2684               energy = MIN2(energy, cx_energy);
2685             }
2686             break;
2687           case 1: /* 1 unpaired base between helices */
2688             dang = (dangle_model==2)?(dang3+dang5):MIN2(dang3, dang5);
2689             if (dangle_model==3) {
2690               energy = energy +dang; ld5 = dang - dang3;
2691               /* may be problem here: Suppose
2692                  cx_energy>energy, cx_energy+dang5<energy
2693                  and the following helices are also stacked (i.e.
2694                  we'll subtract the dang5 again */
2695               if (cx_energy+dang5 < energy) {
2696                 energy = cx_energy+dang5;
2697                 ld5 = dang5;
2698               }
2699               new_cx = INF;  /* no coax stacking with mismatch for now */
2700             } else
2701               energy += dang;
2702             break;
2703           default: /* many unpaired base between helices */
2704             energy += dang5 +dang3;
2705             if (dangle_model==3) {
2706               energy = MIN2(energy, cx_energy + dang5);
2707               new_cx = INF;  /* no coax stacking possible */
2708               ld5 = dang5;
2709             }
2710           }
2711           type = tt;
2712         }
2713         if (dangle_model==3) cx_energy = new_cx;
2714         i1 = q; p=q+1;
2715       } while (q!=i);
2716       best_energy = MIN2(energy, best_energy); /* don't use cx_energy here */
2717       /* fprintf(stderr, "%6.2d\t", energy); */
2718       if (dangle_model!=3 || is_extloop) break;  /* may break cofold with co-ax */
2719       /* skip a helix and start again */
2720       while (pair_table[p]==0) p++;
2721       if (i == (unsigned int)pair_table[p]) break;
2722       i = (unsigned int)pair_table[p];
2723     }
2724     energy = best_energy;
2725     energy += mlclosing;
2726     /* logarithmic ML loop energy if logML */
2727     if ( (!is_extloop) && logML && (u>6) )
2728       energy += 6*mlbase+(int)(P->lxc*log((double)u/6.));
2729     else
2730       energy += mlbase*u;
2731     /* fprintf(stderr, "\n"); */
2732   }
2733   return energy;
2734 }
2735
2736 PUBLIC void initialize_fold(int length){
2737   /* DO NOTHING */
2738 }
2739
2740 PUBLIC float energy_of_struct(const char *string, const char *structure){
2741   return energy_of_structure(string, structure, eos_debug);
2742 }
2743
2744 PUBLIC int energy_of_struct_pt(const char *string, short * ptable, short *s, short *s1){
2745   return energy_of_structure_pt(string, ptable, s, s1, eos_debug);
2746 }
2747
2748 PUBLIC float energy_of_circ_struct(const char *string, const char *structure){
2749   return energy_of_circ_structure(string, structure, eos_debug);
2750 }
2751