binaries/src/ViennaRNA/lib/Lfold.c

   1 /* Last changed Time-stamp: <2007-09-04 11:16:01 ivo> */
   2 /*
   3                   minimum free energy
   4                   RNA secondary structure prediction
   5                   with maximum distance base pairs
   6
   7                   c Ivo Hofacker, Peter Stadler
   8
   9                   Vienna RNA package
  10 */
  11
  12 #include <config.h>
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 #include <ctype.h>
  17 #include <string.h>
  18 #include "utils.h"
  19 #include "energy_par.h"
  20 #include "fold_vars.h"
  21 #include "pair_mat.h"
  22 #include "params.h"
  23 #include "loop_energies.h"
  24 #include "gquad.h"
  25 #include "Lfold.h"
  26
  27 #ifdef USE_SVM
  28 #include "svm.h"
  29 #include "svm_utils.h"
  30 #endif
  31
  32 #ifdef _OPENMP
  33 #include <omp.h>
  34 #endif
  35
  36 #define PAREN
  37
  38 #define STACK_BULGE1      1   /* stacking energies for bulges of size 1 */
  39 #define NEW_NINIO         1   /* new asymetry penalty */
  40 /*@unused@*/
  41 #define MAXSECTORS        500     /* dimension for a backtrack array */
  42 #define LOCALITY          0.      /* locality parameter for base-pairs */
  43
  44 /*
  45 #################################
  46 # GLOBAL VARIABLES              #
  47 #################################
  48 */
  49
  50 /*
  51 #################################
  52 # PRIVATE VARIABLES             #
  53 #################################
  54 */
  55 PRIVATE paramT        *P = NULL;
  56 PRIVATE int           **c = NULL;        /* energy array, given that i-j pair */
  57 PRIVATE int           *cc = NULL;        /* linear array for calculating canonical structures */
  58 PRIVATE int           *cc1 = NULL;       /*   "     "        */
  59 PRIVATE int           *f3 = NULL;        /* energy of 5' end */
  60 PRIVATE int           **fML = NULL;      /* multi-loop auxiliary energy array */
  61 PRIVATE int           *Fmi = NULL;       /* holds row i of fML (avoids jumps in memory) */
  62 PRIVATE int           *DMLi = NULL;      /* DMLi[j] holds MIN(fML[i,k]+fML[k+1,j])  */
  63 PRIVATE int           *DMLi1 = NULL;     /*             MIN(fML[i+1,k]+fML[k+1,j])  */
  64 PRIVATE int           *DMLi2 = NULL;     /*             MIN(fML[i+2,k]+fML[k+1,j])  */
  65 PRIVATE char          **ptype = NULL;    /* precomputed array of pair types */
  66 PRIVATE short         *S = NULL, *S1 = NULL;
  67 PRIVATE unsigned int  length;
  68 PRIVATE char          *prev = NULL;
  69 #ifdef USE_SVM
  70 PRIVATE struct svm_model  *avg_model = NULL;
  71 PRIVATE struct svm_model  *sd_model = NULL;
  72 #endif
  73
  74 PRIVATE int           with_gquad  = 0;
  75 PRIVATE int           **ggg       = NULL;
  76
  77
  78 #ifdef _OPENMP
  79
  80 #ifdef USE_SVM
  81 #pragma omp threadprivate(P, c, cc, cc1, f3, fML, Fmi, DMLi, DMLi1, DMLi2, ptype, S, S1, length, sd_model, avg_model, ggg, with_gquad)
  82 #else
  83 #pragma omp threadprivate(P, c, cc, cc1, f3, fML, Fmi, DMLi, DMLi1, DMLi2, ptype, S, S1, length, ggg, with_gquad)
  84 #endif
  85
  86 #endif
  87
  88 /*
  89 #################################
  90 # PRIVATE FUNCTION DECLARATIONS #
  91 #################################
  92 */
  93 PRIVATE void  initialize_Lfold(int length, int maxdist);
  94 PRIVATE void  update_fold_params(void);
  95 PRIVATE void  get_arrays(unsigned int size, int maxdist);
  96 PRIVATE void  free_arrays(int maxdist);
  97 PRIVATE void  make_ptypes(const short *S, int i, int maxdist, int n);
  98 PRIVATE char  *backtrack(const char *sequence, int start, int maxdist);
  99 PRIVATE int   fill_arrays(const char *sequence, int maxdist, int zsc, double min_z);
 100
 101 /*
 102 #################################
 103 # BEGIN OF FUNCTION DEFINITIONS #
 104 #################################
 105 */
 106
 107 /*--------------------------------------------------------------------------*/
 108 PRIVATE void initialize_Lfold(int length, int maxdist){
 109
 110   if (length<1) nrerror("initialize_Lfold: argument must be greater 0");
 111   get_arrays((unsigned) length, maxdist);
 112   update_fold_params();
 113 }
 114
 115 /*--------------------------------------------------------------------------*/
 116 PRIVATE void get_arrays(unsigned int size, int maxdist){
 117   int i;
 118   c     = (int **)  space(sizeof(int *) *(size+1));
 119   fML   = (int **)  space(sizeof(int *) *(size+1));
 120   ptype = (char **) space(sizeof(char *)*(size+1));
 121   f3    = (int *)   space(sizeof(int)   *(size+2));  /* has to be one longer */
 122   cc    = (int *)   space(sizeof(int)   *(maxdist+5));
 123   cc1   = (int *)   space(sizeof(int)   *(maxdist+5));
 124   Fmi   = (int *)   space(sizeof(int)   *(maxdist+5));
 125   DMLi  = (int *)   space(sizeof(int)   *(maxdist+5));
 126   DMLi1 = (int *)   space(sizeof(int)   *(maxdist+5));
 127   DMLi2 = (int *)   space(sizeof(int)   *(maxdist+5));
 128
 129   for (i=size; (i>(int)size-maxdist-5) && (i>=0); i--) {
 130     c[i] = (int *) space(sizeof(int)*(maxdist+5));
 131   }
 132   for (i=size; (i>(int)size-maxdist-5) && (i>=0); i--) {
 133     fML[i] = (int *) space(sizeof(int)*(maxdist+5));
 134   }
 135   for (i=size; (i>(int)size-maxdist-5) && (i>=0); i--) {
 136     ptype[i] = (char *) space(sizeof(char)*(maxdist+5));
 137   }
 138 }
 139
 140 /*--------------------------------------------------------------------------*/
 141
 142 PRIVATE void free_arrays(int maxdist){
 143   int i;
 144   for (i=0; (i<maxdist+5) && (i<=length); i++){
 145     free(c[i]);
 146     free(fML[i]);
 147     free(ptype[i]);
 148   }
 149   free(c);
 150   free(fML);
 151   free(ptype);
 152   free(f3);
 153   free(cc);
 154   free(cc1);
 155   free(Fmi);
 156   free(DMLi);
 157   free(DMLi1);
 158   free(DMLi2);
 159
 160   if(ggg){
 161     for (i=0; (i<maxdist+5) && (i<=length); i++){
 162       free(ggg[i]);
 163     }
 164     free(ggg);
 165     ggg = NULL;
 166   }
 167
 168   f3    = cc = cc1 = Fmi = DMLi = DMLi1 = DMLi2 = NULL;
 169   c     = fML = NULL;
 170   ptype = NULL;
 171 }
 172
 173 /*--------------------------------------------------------------------------*/
 174
 175 PUBLIC  float Lfold(const char *string, char *structure, int maxdist){
 176   return Lfoldz(string, structure, maxdist, 0, 0.0);
 177 }
 178
 179 PUBLIC  float Lfoldz(const char *string, char *structure, int maxdist, int zsc, double min_z){
 180   int i, energy;
 181
 182   length = (int) strlen(string);
 183   if (maxdist>length) maxdist = length;
 184   initialize_Lfold(length, maxdist);
 185   if (fabs(P->temperature - temperature)>1e-6) update_fold_params();
 186
 187   with_gquad  = P->model_details.gquad;
 188   S           = encode_sequence(string, 0);
 189   S1          = encode_sequence(string, 1);
 190
 191   for (i=length; i>=(int)length-(int)maxdist-4 && i>0; i--)
 192     make_ptypes(S, i, maxdist, length);
 193
 194 #ifdef USE_SVM  /*svm*/
 195   if(zsc){
 196     avg_model = svm_load_model_string(avg_model_string);
 197     sd_model  = svm_load_model_string(sd_model_string);
 198   }
 199 #endif
 200
 201   energy = fill_arrays(string, maxdist, zsc, min_z);
 202
 203 #ifdef USE_SVM  /*svm*/
 204   if(zsc){
 205     svm_destroy_model(avg_model);
 206     svm_destroy_model(sd_model);
 207   }
 208 #endif
 209
 210   free(S); free(S1);
 211   free_arrays(maxdist);
 212
 213   return (float) energy/100.;
 214 }
 215
 216 PRIVATE int fill_arrays(const char *string, int maxdist, int zsc, double min_z) {
 217   /* fill "c", "fML" and "f3" arrays and return  optimal energy */
 218
 219   int   i, j, k, length, energy;
 220   int   decomp, new_fML;
 221   int   no_close, type, type_2, tt;
 222   int   fij;
 223   int   lind;
 224
 225   length = (int) strlen(string);
 226   prev=NULL;
 227   for (j=0; j<maxdist+5; j++)
 228     Fmi[j]=DMLi[j]=DMLi1[j]=DMLi2[j]=INF;
 229   for (j=length; j>length-maxdist-4; j--) {
 230     for (i=(length-maxdist-4>0)?length-maxdist-4:1 ; i<j; i++)
 231       c[i][j-i] = fML[i][j-i] = INF;
 232   }
 233
 234   if(with_gquad){
 235     ggg = NULL;
 236     ggg = get_gquad_L_matrix(S, length - maxdist - 4, maxdist, ggg, P);
 237   }
 238
 239   for (i = length-TURN-1; i >= 1; i--) { /* i,j in [1..length] */
 240     for (j = i+TURN+1; j <= length && j <= i+maxdist; j++) {
 241       int p, q;
 242       type = ptype[i][j-i];
 243
 244       no_close = (((type==3)||(type==4))&&no_closingGU);
 245
 246       if (type) {   /* we have a pair */
 247         int new_c=0, stackEnergy=INF;
 248         /* hairpin ----------------------------------------------*/
 249
 250         new_c = (no_close) ? FORBIDDEN : E_Hairpin(j-i-1, type, S1[i+1], S1[j-1], string+i-1, P);
 251
 252         /*--------------------------------------------------------
 253           check for elementary structures involving more than one
 254           closing pair.
 255           --------------------------------------------------------*/
 256
 257         for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1) ; p++){
 258           int minq = j-i+p-MAXLOOP-2;
 259           if (minq<p+1+TURN) minq = p+1+TURN;
 260           for (q = minq; q < j; q++) {
 261             type_2 = ptype[p][q-p];
 262
 263             if (type_2==0) continue;
 264             type_2 = rtype[type_2];
 265
 266             if (no_closingGU)
 267               if (no_close||(type_2==3)||(type_2==4))
 268                 if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
 269
 270             energy = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
 271             new_c = MIN2(new_c, energy + c[p][q-p]);
 272             if ((p==i+1)&&(j==q+1)) stackEnergy = energy; /* remember stack energy */
 273           } /* end q-loop */
 274         } /* end p-loop */
 275
 276         /* multi-loop decomposition ------------------------*/
 277         if (!no_close) {
 278           decomp  = DMLi1[j-1-(i+1)];
 279           tt      = rtype[type];
 280           switch(dangles){
 281             /* no dangles */
 282             case 0:   decomp += E_MLstem(tt, -1, -1, P);
 283                       break;
 284             /* double dangles */
 285             case 2:   decomp += E_MLstem(tt, S1[j-1], S1[i+1], P);
 286                       break;
 287             /* normal dangles, aka dangles = 1 */
 288             default:  decomp += E_MLstem(tt, -1, -1, P);
 289                       decomp = MIN2(decomp, DMLi2[j-1-(i+2)] + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase);
 290                       decomp = MIN2(decomp, DMLi2[j-2-(i+2)] + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase);
 291                       decomp = MIN2(decomp, DMLi1[j-2-(i+1)] + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase);
 292                       break;
 293           }
 294           new_c = MIN2(new_c, decomp + P->MLclosing);
 295         }
 296
 297         /* coaxial stacking of (i.j) with (i+1.k) or (k+1.j-1) */
 298
 299         if (dangles==3) {
 300           decomp = INF;
 301           for (k = i+2+TURN; k < j-2-TURN; k++) {
 302             type_2 = ptype[i+1][k-i-1]; type_2 = rtype[type_2];
 303             if (type_2)
 304               decomp = MIN2(decomp, c[i+1][k-i-1]+P->stack[type][type_2]+
 305                             fML[k+1][j-1-k-1]);
 306             type_2 = ptype[k+1][j-1-k-1]; type_2 = rtype[type_2];
 307             if (type_2)
 308               decomp = MIN2(decomp, c[k+1][j-1-k-1]+P->stack[type][type_2]+
 309                             fML[i+1][k-i-1]);
 310           }
 311           /* no TermAU penalty if coax stack */
 312           decomp += 2*P->MLintern[1] + P->MLclosing;
 313           new_c = MIN2(new_c, decomp);
 314         }
 315
 316         if(with_gquad){
 317           /* include all cases where a g-quadruplex may be enclosed by base pair (i,j) */
 318           if (!no_close) {
 319             tt = rtype[type];
 320             energy = E_GQuad_IntLoop_L(i, j, type, S1, ggg, maxdist, P);
 321             new_c = MIN2(new_c, energy);
 322           }
 323         }
 324
 325         new_c = MIN2(new_c, cc1[j-1-(i+1)]+stackEnergy);
 326         cc[j-i] = new_c;
 327         if (noLonelyPairs)
 328           c[i][j-i] = cc1[j-1-(i+1)]+stackEnergy;
 329         else
 330           c[i][j-i] = cc[j-i];
 331
 332       } /* end >> if (pair) << */
 333
 334       else c[i][j-i] = INF;
 335
 336       /* done with c[i,j], now compute fML[i,j] */
 337       /* free ends ? -----------------------------------------*/
 338       new_fML = INF;
 339       switch(dangles){
 340         /* no dangles */
 341         case 0:   new_fML = fML[i+1][j-i-1] + P->MLbase;
 342                   new_fML = MIN2(new_fML, fML[i][j-1-i] + P->MLbase);
 343                   new_fML = MIN2(new_fML, c[i][j-i] + E_MLstem(type, -1, -1, P));
 344                   break;
 345         /* double dangles */
 346         case 2:   new_fML = fML[i+1][j-i-1] + P->MLbase;
 347                   new_fML = MIN2(fML[i][j-1-i] + P->MLbase, new_fML);
 348                   new_fML = MIN2(new_fML,  c[i][j-i] + E_MLstem(type, (i>1) ? S1[i-1] : -1, (j<length) ? S1[j+1] : -1, P));
 349                   break;
 350         /* normal dangles, aka dangles = 1 */
 351         default:  /* i unpaired */
 352                   new_fML = fML[i+1][j-i-1] + P->MLbase;
 353                   /* j unpaired */
 354                   new_fML = MIN2(new_fML, fML[i][j-1-i] + P->MLbase);
 355                   /* i,j */
 356                   if(type) new_fML = MIN2(new_fML, c[i][j-i] + E_MLstem(type, -1, -1, P));
 357                   /* i+1,j */
 358                   tt = ptype[i+1][j-i-1];
 359                   if(tt) new_fML = MIN2(new_fML, c[i+1][j-i-1] + E_MLstem(tt, S1[i], -1, P) + P->MLbase);
 360                   /* i, j-1 */
 361                   tt = ptype[i][j-1-i];
 362                   if(tt) new_fML = MIN2(new_fML, c[i][j-1-i] + E_MLstem(tt, -1, S1[j], P) + P->MLbase);
 363                   /* i+1,j-1 */
 364                   tt = ptype[i+1][j-1-i-1];
 365                   if(tt) new_fML = MIN2(new_fML, c[i+1][j-1-i-1] + E_MLstem(tt, S1[i], S1[j], P) + 2*P->MLbase);
 366                   break;
 367       }
 368
 369       if(with_gquad){
 370         new_fML = MIN2(new_fML, ggg[i][j - i] + E_MLstem(0, -1, -1, P));
 371       }
 372
 373       /* modular decomposition -------------------------------*/
 374       for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++)
 375         decomp = MIN2(decomp, Fmi[k-i]+fML[k+1][j-k-1]);
 376
 377       DMLi[j-i] = decomp;               /* store for use in ML decompositon */
 378       new_fML   = MIN2(new_fML, decomp);
 379
 380       /* coaxial stacking */
 381       if (dangles==3) {
 382         /* additional ML decomposition as two coaxially stacked helices */
 383         for (decomp = INF, k = i+1+TURN; k <= j-2-TURN; k++) {
 384           type = ptype[i][k-i]; type = rtype[type];
 385           type_2 = ptype[k+1][j-k-1]; type_2 = rtype[type_2];
 386           if (type && type_2)
 387             decomp = MIN2(decomp,
 388                           c[i][k-i]+c[k+1][j-k-1]+P->stack[type][type_2]);
 389         }
 390
 391         decomp += 2*P->MLintern[1];          /* no TermAU penalty if coax stack */
 392 #if 0
 393         /* This is needed for Y shaped ML loops with coax stacking of
 394            interior pairts, but backtracking will fail if activated */
 395         DMLi[j-i] = MIN2(DMLi[j-i], decomp);
 396         DMLi[j-i] = MIN2(DMLi[j-i], DMLi[j-1-i]+P->MLbase);
 397         DMLi[j-i] = MIN2(DMLi[j-i], DMLi1[j-(i+1)]+P->MLbase);
 398         new_fML = MIN2(new_fML, DMLi[j-i]);
 399 #endif
 400         new_fML = MIN2(new_fML, decomp);
 401       }
 402       fML[i][j-i] = Fmi[j-i] = new_fML;     /* substring energy */
 403     } /* for (j...) */
 404
 405     /* calculate energies of 5' and 3' fragments */
 406     {
 407       static int do_backtrack = 0, prev_i=0;
 408       char *ss=NULL;
 409       double prevz;
 410       f3[i] = f3[i+1];
 411       switch(dangles){
 412         /* dont use dangling end and mismatch contributions at all */
 413         case 0:   for(j=i+TURN+1; j<length && j<=i+maxdist; j++){
 414                     type = ptype[i][j-i];
 415
 416                     if(with_gquad){
 417                       f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
 418                     }
 419
 420                     if(type)
 421                       f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, -1, -1, P));
 422                   }
 423                   if(length<=i+maxdist){
 424                     j=length;
 425
 426                     if(with_gquad){
 427                       f3[i] = MIN2(f3[i], ggg[i][j-i]);
 428                     }
 429
 430                     type = ptype[i][j-i];
 431                     if(type)
 432                       f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, -1, -1, P));
 433                   }
 434                   break;
 435         /* always use dangles on both sides */
 436         case 2:   for(j=i+TURN+1; j<length && j<=i+maxdist; j++){
 437                     type = ptype[i][j-i];
 438
 439                     if(with_gquad){
 440                       if(ggg[i][j-i] != INF)
 441                         f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
 442                     }
 443
 444                     if(type)
 445                       f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, S1[j+1], P));
 446                   }
 447                   if(length<=i+maxdist){
 448                     j=length;
 449
 450                     if(with_gquad){
 451                       f3[i] = MIN2(f3[i], ggg[i][j-i]);
 452                     }
 453
 454                     type = ptype[i][j-i];
 455                     if(type)
 456                       f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, -1, P));
 457                   }
 458                   break;
 459         /* normal dangles, aka dangles = 1 */
 460         default:  for(j=i+TURN+1; j<length && j<=i+maxdist; j++){
 461                     type = ptype[i][j-i];
 462
 463                     if(with_gquad){
 464                       f3[i] = MIN2(f3[i], f3[j+1] + ggg[i][j-i]);
 465                     }
 466
 467                     if(type){
 468                       f3[i] = MIN2(f3[i], f3[j+1] + c[i][j-i] + E_ExtLoop(type, -1, -1, P));
 469                       f3[i] = MIN2(f3[i], ((j+2<=length) ? f3[j+2] : 0) + c[i][j-i] + E_ExtLoop(type, -1, S1[j+1], P));
 470                     }
 471                     type = ptype[i+1][j-i-1];
 472                     if(type){
 473                       f3[i] = MIN2(f3[i], f3[j+1] + c[i+1][j-i-1] + E_ExtLoop(type, S1[i], -1, P));
 474                       f3[i] = MIN2(f3[i], ((j + 1 < length) ? f3[j+2] : 0) + c[i+1][j-i-1] + E_ExtLoop(type, S1[i], S1[j+1], P));
 475                     }
 476                   }
 477                   if(length<=i+maxdist){
 478                     j     = length;
 479
 480                     if(with_gquad){
 481                       f3[i] = MIN2(f3[i], ggg[i][j-i]);
 482                     }
 483
 484                     type  = ptype[i][j-i];
 485                     if(type)
 486                       f3[i] = MIN2(f3[i], c[i][j-i] + E_ExtLoop(type, -1, -1, P));
 487                     type  = ptype[i+1][j-i-1];
 488                     if(type)
 489                       f3[i] = MIN2(f3[i], c[i+1][j-i-1] + E_ExtLoop(type, S1[i], -1, P));
 490                   }
 491                   break;
 492       } /* switch(dangles)... */
 493
 494       /* backtrack partial structure */
 495       if (f3[i] < f3[i+1]) do_backtrack=1;
 496       else if (do_backtrack) {
 497         int pairpartner; /*i+1?? is paired with pairpartner*/
 498         int cc;
 499         int traced2=0;
 500         fij = f3[i+1];
 501         lind=i+1;
 502
 503         /*start "short" backtrack*/
 504
 505         /*get paired base*/
 506         while(fij==f3[lind+1]) lind++;
 507
 508         /*get pairpartner*/
 509         for (pairpartner = lind + TURN; pairpartner <= lind + maxdist; pairpartner++){
 510           type = ptype[lind][pairpartner-lind];
 511           switch(dangles){
 512             case 0:   if(type){
 513                         cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
 514                         if(fij == cc + f3[pairpartner + 1])
 515                           traced2 = 1;
 516                       }
 517                       else if(with_gquad) {
 518                         cc = ggg[lind][pairpartner-lind];
 519                         if(fij == cc + f3[pairpartner + 1])
 520                           traced2 = 1;
 521                       }
 522
 523                       break;
 524             case 2:   if(type){
 525                         cc = c[lind][pairpartner-lind] + E_ExtLoop(type, (lind > 1) ? S1[lind-1] : -1, (pairpartner < length) ? S1[pairpartner+1] : -1, P);
 526                         if(fij == cc + f3[pairpartner + 1])
 527                           traced2 = 1;
 528                       }
 529                       else if(with_gquad){
 530                         cc = ggg[lind][pairpartner-lind];
 531                         if(fij == cc + f3[pairpartner + 1])
 532                           traced2 = 1;
 533                       }
 534
 535                       break;
 536             default:  if(type){
 537                         cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
 538                         if(fij == cc + f3[pairpartner + 1]){
 539                           traced2 = 1;
 540                           break;
 541                         }
 542                         else if(pairpartner < length){
 543                           cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, S1[pairpartner+1], P);
 544                           if(fij == cc + f3[pairpartner + 2]){
 545                             traced2 = 1;
 546                             break;
 547                           }
 548                         }
 549                       }
 550                       else if(with_gquad){
 551                         cc = ggg[lind][pairpartner-lind];
 552                         if(fij == cc + f3[pairpartner + 1])
 553                           traced2 = 1;
 554                       }
 555
 556                       type = ptype[lind+1][pairpartner-lind-1];
 557                       if(type){
 558                         cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], -1, P);
 559                         if(fij == cc + f3[pairpartner+1]){
 560                           traced2 = 1;
 561                           break;
 562                         }
 563                         else if(pairpartner < length){
 564                           cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], S1[pairpartner+1], P);
 565                           if(fij == cc + f3[pairpartner+2])
 566                             traced2 = 1;
 567                         }
 568                       }
 569                       break;
 570           }
 571           if(traced2) break;
 572         }
 573         if (!traced2) nrerror("backtrack failed in short backtrack 1");
 574         if (zsc){
 575 #ifdef USE_SVM
 576           int info_avg;
 577           double average_free_energy;
 578           double sd_free_energy;
 579           double my_z;
 580           int *AUGC = get_seq_composition(S, lind-1, MIN2((pairpartner+1),length));
 581           /*\svm*/
 582           average_free_energy = avg_regression(AUGC[0], AUGC[1], AUGC[2], AUGC[3], AUGC[4], avg_model, &info_avg);
 583           if (info_avg == 0)  {
 584             double difference;
 585             double min_sd = minimal_sd(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4]);
 586             difference=(fij-f3[pairpartner+1])/100.-average_free_energy;
 587             if ( difference - ( min_z * min_sd ) <= 0.0001 ) {
 588               sd_free_energy = sd_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],sd_model);
 589               my_z=difference/sd_free_energy;
 590               if (my_z<=min_z){
 591                 ss =  backtrack(string, lind , pairpartner+1);
 592                 if (prev) {
 593                   if ((i+strlen(ss)<prev_i+strlen(prev)) ||
 594                       strncmp(ss+prev_i-i,prev,strlen(prev))) { /* ss does not contain prev */
 595                     if (dangles==2)
 596                       printf(".%s (%6.2f) %4d z= %.3f\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1, prevz);
 597                     else
 598                       printf("%s (%6.2f) %4d z=%.3f\n ", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i, prevz);
 599                   }
 600                   free(prev);
 601                 }
 602                 prev=ss; prev_i = lind; prevz=my_z;
 603               }
 604             }
 605
 606           }
 607           free(AUGC);
 608           do_backtrack=0;
 609 #endif
 610         }
 611         else {
 612           /* original code for Lfold*/
 613           ss =  backtrack(string, lind , pairpartner+1);
 614           if (prev) {
 615             if ((i+strlen(ss)<prev_i+strlen(prev)) || strncmp(ss+prev_i-i,prev,strlen(prev))){
 616               /* ss does not contain prev */
 617               if (dangles==2)
 618                 printf(".%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1);
 619               else
 620                 printf("%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i);
 621             }
 622             free(prev);
 623           }
 624           prev=ss;
 625           prev_i = lind;
 626           do_backtrack=0;
 627         }
 628       }
 629       if (i==1) {
 630         if (prev) {
 631           if(zsc) {
 632             if (dangles==2)
 633               printf(".%s (%6.2f) %4d z= %.2f\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1, prevz);
 634            else
 635               printf("%s (%6.2f) %4dz= %.2f \n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i, prevz);
 636             free(prev); prev=NULL;
 637           }
 638           else {
 639             if (dangles==2)
 640               printf(".%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)-1])/100., prev_i-1);
 641             else
 642               printf("%s (%6.2f) %4d\n", prev, (f3[prev_i]-f3[prev_i+strlen(prev)])/100., prev_i);
 643           }
 644         } else if (f3[i]<0)do_backtrack=1;
 645
 646         if (do_backtrack) {
 647           int pairpartner; /*i+1?? is paired with pairpartner*/
 648           int cc;
 649           double average_free_energy;
 650           double sd_free_energy;
 651           int info_avg;
 652           double my_z;
 653           int traced2 = 0;
 654           fij = f3[i];
 655           lind=i;
 656           while(fij==f3[lind+1]) lind++;
 657           /*get pairpartner*/
 658           for(pairpartner = lind + TURN; pairpartner <= lind + maxdist; pairpartner++){
 659             type = ptype[lind][pairpartner-lind];
 660             switch(dangles){
 661               case 0:   if(type){
 662                           cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
 663                           if(fij == cc + f3[pairpartner + 1])
 664                             traced2 = 1;
 665                         }
 666                         else if(with_gquad){
 667                           cc = ggg[lind][pairpartner-lind];
 668                           if(fij == cc + f3[pairpartner + 1])
 669                             traced2 = 1;
 670                         }
 671
 672                         break;
 673               case 2:   if(type){
 674                           cc = c[lind][pairpartner-lind] + E_ExtLoop(type, (lind > 1) ? S1[lind-1] : -1, (pairpartner < length) ? S1[pairpartner+1] : -1, P);
 675                           if(fij == cc + f3[pairpartner + 1])
 676                             traced2 = 1;
 677                         }
 678                         else if(with_gquad){
 679                           cc = ggg[lind][pairpartner-lind];
 680                           if(fij == cc + f3[pairpartner + 1])
 681                             traced2 = 1;
 682                         }
 683
 684                         break;
 685               default:  if(type){
 686                           cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, -1, P);
 687                           if(fij == cc + f3[pairpartner + 1]){
 688                             traced2 = 1;
 689                             break;
 690                           }
 691                           else if(pairpartner < length){
 692                             cc = c[lind][pairpartner-lind] + E_ExtLoop(type, -1, S1[pairpartner + 1], P);
 693                             if(fij == cc + f3[pairpartner + 1]){
 694                               traced2 = 1;
 695                               break;
 696                             }
 697                           }
 698                         }
 699                         else if(with_gquad){
 700                           cc = ggg[lind][pairpartner-lind];
 701                           if(fij == cc + f3[pairpartner + 1])
 702                             traced2 = 1;
 703                         }
 704
 705                         type = ptype[lind+1][pairpartner-lind-1];
 706                         if(type){
 707                           cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], -1, P);
 708                           if(fij == cc + f3[pairpartner+1]){
 709                             traced2 = 1;
 710                             break;
 711                           }
 712                           else if (pairpartner < length){
 713                             cc = c[lind+1][pairpartner-(lind+1)] + E_ExtLoop(type, S1[lind], S1[pairpartner+1], P);
 714                             if(fij == cc + f3[pairpartner + 2]){
 715                               traced2 =1;
 716                               break;
 717                             }
 718                           }
 719                         }
 720             }
 721             if(traced2) break;
 722           }
 723           if (!traced2) nrerror("backtrack failed in short backtrack 2");
 724
 725           if(zsc){
 726 #ifdef USE_SVM
 727             int *AUGC = get_seq_composition(S, lind-1, MIN2((pairpartner+1),length));
 728             average_free_energy = avg_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],avg_model,&info_avg);
 729             if (info_avg == 0)  {
 730               double difference;
 731               double min_sd = minimal_sd(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4]);
 732               difference=(fij-f3[pairpartner+1])/100.-average_free_energy;
 733               if ( difference - ( min_z * min_sd ) <= 0.0001 ) {
 734                 sd_free_energy = sd_regression(AUGC[0],AUGC[1],AUGC[2],AUGC[3],AUGC[4],sd_model);
 735                 my_z=difference/sd_free_energy;
 736                 if (my_z<=min_z){
 737                   ss =  backtrack(string, lind , pairpartner+1);
 738                   printf("%s (%6.2f) %4d z= %.2f\n", ss, (f3[lind]-f3[lind+strlen(ss)-1])/100., lind, my_z);
 739                 }
 740               }
 741             }
 742             free(AUGC);
 743 #endif
 744           }
 745           else {
 746             ss =  backtrack(string, lind , pairpartner+1);
 747             if (dangles==2)
 748               printf("%s (%6.2f) %4d\n", ss, (f3[lind]-f3[lind+strlen(ss)-1])/100., 1);
 749             else
 750               printf("%s (%6.2f) %4d\n", ss, (f3[lind]-f3[lind+strlen(ss)])/100., 1);
 751             free(ss);
 752           }
 753         }
 754         do_backtrack=0;
 755       }
 756     }
 757     {
 758       int ii, *FF; /* rotate the auxilliary arrays */
 759       FF = DMLi2; DMLi2 = DMLi1; DMLi1 = DMLi; DMLi = FF;
 760       FF = cc1; cc1=cc; cc=FF;
 761       for (j=0; j< maxdist+5; j++) {cc[j]=Fmi[j]=DMLi[j]=INF; }
 762       if (i+maxdist+4<=length) {
 763         c[i-1] = c[i+maxdist+4]; c[i+maxdist+4] = NULL;
 764         fML[i-1] = fML[i+maxdist+4]; fML[i+maxdist+4]=NULL;
 765         ptype[i-1] = ptype[i+maxdist+4]; ptype[i+maxdist+4] = NULL;
 766         if (i>1){
 767           make_ptypes(S, i-1, maxdist, length);
 768
 769           if(with_gquad){
 770             ggg = get_gquad_L_matrix(S, i - 1, maxdist, ggg, P);
 771           }
 772         }
 773         for (ii=0; ii<maxdist+5; ii++) {
 774           c[i-1][ii] = INF;
 775           fML[i-1][ii] = INF;
 776         }
 777       }
 778     }
 779   }
 780
 781   return f3[1];
 782 }
 783
 784 PRIVATE char *backtrack(const char *string, int start, int maxdist){
 785   /*------------------------------------------------------------------
 786     trace back through the "c", "f3" and "fML" arrays to get the
 787     base pairing list. No search for equivalent structures is done.
 788     This is fast, since only few structure elements are recalculated.
 789     ------------------------------------------------------------------*/
 790   sect  sector[MAXSECTORS];   /* backtracking sectors */
 791   int   i, j, k, energy, new, no_close, type, type_2, tt, s=0;
 792   char  *structure;
 793
 794   /* length = strlen(string); */
 795   sector[++s].i = start;
 796   sector[s].j   = MIN2(length, maxdist+1);
 797   sector[s].ml  = (backtrack_type=='M') ? 1 : ((backtrack_type=='C')?2:0);
 798
 799   structure = (char *) space((MIN2(length-start, maxdist)+3)*sizeof(char));
 800   for (i=0; i<=MIN2(length-start, maxdist); i++) structure[i] = '-';
 801
 802   while (s>0) {
 803     int ml, fij, cij, traced, i1, j1, d3, d5, mm, mm5, mm3, mm53, p, q, jj=0, gq=0;
 804     int canonical = 1;     /* (i,j) closes a canonical structure */
 805     i  = sector[s].i;
 806     j  = sector[s].j;
 807     ml = sector[s--].ml;   /* ml is a flag indicating if backtracking is to
 808                               occur in the fML- (1) or in the f-array (0) */
 809     if (ml==2) {
 810       structure[i-start] = '(';
 811       structure[j-start] = ')';
 812       goto repeat1;
 813     }
 814
 815     if (j < i+TURN+1) continue; /* no more pairs in this interval */
 816
 817     fij = (ml)? fML[i][j-i] : f3[i];
 818
 819     if (ml == 0) { /* backtrack in f3 */
 820
 821       if (fij == f3[i+1]) {
 822         sector[++s].i = i+1;
 823         sector[s].j   = j;
 824         sector[s].ml  = ml;
 825         continue;
 826       }
 827       /* i or i+1 is paired. Find pairing partner */
 828       switch(dangles){
 829         case 0:   for(traced = 0, k=j; k>i+TURN; k--){
 830
 831                     if(with_gquad){
 832                       if(fij == ggg[i][k-i] + f3[k+1]){
 833                         /* found the decomposition */
 834                         traced = i; jj = k + 1; gq = 1;
 835                         break;
 836                       }
 837                     }
 838
 839                     jj    = k+1;
 840                     type  = ptype[i][k-i];
 841                     if(type)
 842                       if(fij == c[i][k-i] + E_ExtLoop(type, -1, -1, P) + f3[k+1]){
 843                         traced = i;
 844                         break;
 845                       }
 846                   }
 847                   break;
 848         case 2:   for(traced = 0, k=j; k>i+TURN; k--){
 849
 850                     if(with_gquad){
 851                       if(fij == ggg[i][k-i] + f3[k+1]){
 852                         /* found the decomposition */
 853                         traced = i; jj = k + 1; gq = 1;
 854                         break;
 855                       }
 856                     }
 857
 858                     jj    = k+1;
 859                     type  = ptype[i][k-i];
 860                     if(type)
 861                       if(fij == c[i][k-i] + E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (k<length) ? S1[k+1] : -1, P) + f3[k+1]){
 862                         traced = i;
 863                         break;
 864                       }
 865                   }
 866                   break;
 867         default:  for(traced = 0,k=j; k>i+TURN; k--){
 868
 869                     if(with_gquad){
 870                       if(fij == ggg[i][k-i] + f3[k+1]){
 871                         /* found the decomposition */
 872                         traced = i; jj = k + 1; gq = 1;
 873                         break;
 874                       }
 875                     }
 876
 877                     jj = k+1;
 878                     type = ptype[i+1][k-(i+1)];
 879                     if(type){
 880                       if(fij == c[i+1][k-(i+1)] + E_ExtLoop(type, S1[i], -1, P) + f3[k+1]){
 881                         traced=i+1;
 882                       }
 883                       if(k < length){
 884                         if(fij == c[i+1][k-(i+1)] + E_ExtLoop(type, S1[i], S1[k+1], P) + f3[k+2]){
 885                           traced  = i+1;
 886                           jj      = k+2;
 887                         }
 888                       }
 889                     }
 890                     type = ptype[i][k-i];
 891                     if(type){
 892                       if(fij == c[i][k-i] + E_ExtLoop(type, -1, -1, P) + f3[k+1]){
 893                         traced = i;
 894                       }
 895                       if(k<length){
 896                         if(fij == c[i][k-i] + E_ExtLoop(type, -1, S1[k+1], P) + f3[k+2]){
 897                           traced  = i;
 898                           jj      = k+2;
 899                         }
 900                       }
 901                     }
 902                     if(traced) break;
 903                   }
 904                   break;
 905       } /* switch(dangles)...*/
 906
 907       if (!traced) nrerror("backtrack failed in f3");
 908       if (j==length) { /* backtrack only one component, unless j==length */
 909         sector[++s].i = jj;
 910         sector[s].j   = j;
 911         sector[s].ml  = ml;
 912       }
 913       i=traced; j=k;
 914
 915       if(with_gquad && gq){
 916         /* goto backtrace of gquadruplex */
 917         goto repeat_gquad;
 918       }
 919
 920       structure[i-start] = '('; structure[j-start] = ')';
 921       if (((jj==j+2) || (dangles==2)) && (j < length)) structure[j+1-start] = '.';
 922       goto repeat1;
 923     }
 924     else { /* trace back in fML array */
 925       if (fML[i][j-1-i]+P->MLbase == fij) {  /* 3' end is unpaired */
 926         sector[++s].i = i;
 927         sector[s].j   = j-1;
 928         sector[s].ml  = ml;
 929         continue;
 930       }
 931       if (fML[i+1][j-(i+1)]+P->MLbase == fij) { /* 5' end is unpaired */
 932         sector[++s].i = i+1;
 933         sector[s].j   = j;
 934         sector[s].ml  = ml;
 935         continue;
 936       }
 937
 938       if(with_gquad){
 939         if(fij == ggg[i][j-i] + E_MLstem(0, -1, -1, P)){
 940           /* go to backtracing of quadruplex */
 941           goto repeat_gquad;
 942         }
 943       }
 944
 945       switch(dangles){
 946         case 0:   tt = ptype[i][j-i];
 947                   if(fij == c[i][j-i] + E_MLstem(tt, -1, -1, P)){
 948                     structure[i-start] = '(';
 949                     structure[j-start] = ')';
 950                     goto repeat1;
 951                   }
 952                   break;
 953         case 2:   tt = ptype[i][j-i];
 954                   if(fij == c[i][j-i] + E_MLstem(tt, (i>1) ? S1[i-1] : -1, (j < length) ? S1[j+1] : -1, P)){
 955                     structure[i-start] = '(';
 956                     structure[j-start] = ')';
 957                     goto repeat1;
 958                   }
 959                   break;
 960         default:  tt = ptype[i][j-i];
 961                   if(fij == c[i][j-i] + E_MLstem(tt, -1, -1, P)){
 962                     structure[i-start] = '(';
 963                     structure[j-start] = ')';
 964                     goto repeat1;
 965                   }
 966                   tt = ptype[i+1][j-(i+1)];
 967                   if(fij == c[i+1][j-(i+1)] + E_MLstem(tt, S1[i], -1, P) + P->MLbase){
 968                     structure[++i-start] = '(';
 969                     structure[j-start] = ')';
 970                     goto repeat1;
 971                   }
 972                   tt = ptype[i][j-1-i];
 973                   if(fij == c[i][j-1-i] + E_MLstem(tt, -1, S1[j], P) + P->MLbase){
 974                     structure[i-start] = '(';
 975                     structure[--j-start] = ')';
 976                     goto repeat1;
 977                   }
 978                   tt = ptype[i+1][j-1-(i+1)];
 979                   if(fij == c[i+1][j-1-(i+1)] + E_MLstem(tt, S1[i], S1[j], P) + 2*P->MLbase){
 980                     structure[++i-start] = '(';
 981                     structure[--j-start] = ')';
 982                     goto repeat1;
 983                   }
 984                   break;
 985       } /* switch(dangles)... */
 986
 987       /* modular decomposition */
 988       for (k = i+1+TURN; k <= j-2-TURN; k++)
 989         if (fij == (fML[i][k-i]+fML[k+1][j-(k+1)]))
 990           break;
 991
 992       if ((dangles==3)&&(k>j-2-TURN)) { /* must be coax stack */
 993         ml = 2;
 994         for (k = i+1+TURN; k <= j-2-TURN; k++) {
 995           type = ptype[i][k-i];  type= rtype[type];
 996           type_2 = ptype[k+1][j-(k+1)]; type_2= rtype[type_2];
 997           if (type && type_2)
 998             if (fij == c[i][k-i]+c[k+1][j-(k+1)]+P->stack[type][type_2]+
 999                        2*P->MLintern[1])
1000               break;
1001         }
1002       }
1003
1004       sector[++s].i = i;
1005       sector[s].j   = k;
1006       sector[s].ml  = ml;
1007       sector[++s].i = k+1;
1008       sector[s].j   = j;
1009       sector[s].ml  = ml;
1010
1011       if (k>j-2-TURN) nrerror("backtrack failed in fML");
1012       continue;
1013     }
1014
1015   repeat1:
1016
1017     /*----- begin of "repeat:" -----*/
1018     if (canonical)  cij = c[i][j-i];
1019
1020     type = ptype[i][j-i];
1021
1022
1023     if (noLonelyPairs)
1024       if (cij == c[i][j-i]) {
1025         /* (i.j) closes canonical structures, thus
1026            (i+1.j-1) must be a pair                */
1027         type_2 = ptype[i+1][j-1-(i+1)]; type_2 = rtype[type_2];
1028         cij -= P->stack[type][type_2];
1029         structure[i+1-start] = '('; structure[j-1-start] = ')';
1030         i++; j--;
1031         canonical=0;
1032         goto repeat1;
1033       }
1034     canonical = 1;
1035
1036
1037     no_close = (((type==3)||(type==4))&&no_closingGU);
1038     if (no_close) {
1039       if (cij == FORBIDDEN) continue;
1040     } else
1041       if (cij == E_Hairpin(j-i-1, type, S1[i+1], S1[j-1],string+i-1, P))
1042         continue;
1043
1044     for (p = i+1; p <= MIN2(j-2-TURN,i+MAXLOOP+1); p++) {
1045       int minq;
1046       minq = j-i+p-MAXLOOP-2;
1047       if (minq<p+1+TURN) minq = p+1+TURN;
1048       for (q = j-1; q >= minq; q--) {
1049
1050         type_2 = ptype[p][q-p];
1051         if (type_2==0) continue;
1052         type_2 = rtype[type_2];
1053         if (no_closingGU)
1054           if (no_close||(type_2==3)||(type_2==4))
1055             if ((p>i+1)||(q<j-1)) continue;  /* continue unless stack */
1056
1057         /* energy = oldLoopEnergy(i, j, p, q, type, type_2); */
1058         energy = E_IntLoop(p-i-1, j-q-1, type, type_2, S1[i+1], S1[j-1], S1[p-1], S1[q+1],P);
1059
1060         new = energy+c[p][q-p];
1061         traced = (cij == new);
1062         if (traced) {
1063           structure[p-start] = '(';
1064           structure[q-start] = ')';
1065           i = p, j = q;
1066           goto repeat1;
1067         }
1068       }
1069     }
1070
1071     /* end of repeat: --------------------------------------------------*/
1072
1073     /* (i.j) must close a multi-loop */
1074     tt = rtype[type];
1075     i1 = i+1; j1 = j-1;
1076
1077     if(with_gquad){
1078       /*
1079         The case that is handled here actually resembles something like
1080         an interior loop where the enclosing base pair is of regular
1081         kind and the enclosed pair is not a canonical one but a g-quadruplex
1082         that should then be decomposed further...
1083       */
1084       if(backtrack_GQuad_IntLoop_L(cij, i, j, type, S, ggg, maxdist, &p, &q, P)){
1085         i = p; j = q;
1086         goto repeat_gquad;
1087       }
1088     }
1089
1090     sector[s+1].ml  = sector[s+2].ml = 1;
1091
1092     switch(dangles){
1093       case 0:   mm = P->MLclosing + E_MLstem(tt, -1, -1, P);
1094                 for(k = i+2+TURN; k < j-2-TURN; k++){
1095                   if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
1096                     break;
1097                 }
1098                 break;
1099       case 2:   mm = P->MLclosing + E_MLstem(tt, S1[j-1], S1[i+1], P);
1100                 for(k = i+2+TURN; k < j-2-TURN; k++){
1101                   if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
1102                     break;
1103                 }
1104                 break;
1105       default:  mm    = P->MLclosing + E_MLstem(tt, -1, -1, P);
1106                 mm5   = P->MLclosing + E_MLstem(tt, S1[j-1], -1, P) + P->MLbase;
1107                 mm3   = P->MLclosing + E_MLstem(tt, -1, S1[i+1], P) + P->MLbase;
1108                 mm53  = P->MLclosing + E_MLstem(tt, S1[j-1], S1[i+1], P) + 2*P->MLbase;
1109                 for(k = i+2+TURN; k < j-2-TURN; k++){
1110                   if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-1-(k+1)] + mm)
1111                     break;
1112                   else if(cij == fML[i+2][k-(i+2)] + fML[k+1][j-1-(k+1)] + mm3){
1113                     i1 = i+2;
1114                     break;
1115                   }
1116                   else if(cij == fML[i+1][k-(i+1)] + fML[k+1][j-2-(k+1)] + mm5){
1117                     j1 = j-2;
1118                     break;
1119                   }
1120                   else if(cij == fML[i+2][k-(i+2)] + fML[k+1][j-2-(k+1)] + mm53){
1121                     i1 = i+2;
1122                     j1 = j-2;
1123                     break;
1124                   }
1125                   /* coaxial stacking of (i.j) with (i+1.k) or (k.j-1) */
1126                   /* use MLintern[1] since coax stacked pairs don't get TerminalAU */
1127                   if (dangles==3) {
1128                     int en;
1129                     type_2 = ptype[i+1][k-(i+1)]; type_2 = rtype[type_2];
1130                     if (type_2) {
1131                       en = c[i+1][k-(i+1)]+P->stack[type][type_2]+fML[k+1][j-1-(k+1)];
1132                       if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1133                         ml = 2;
1134                         sector[s+1].ml  = 2;
1135                         break;
1136                       }
1137                     }
1138                     type_2 = ptype[k+1][j-1-(k+1)]; type_2 = rtype[type_2];
1139                     if (type_2) {
1140                       en = c[k+1][j-1-(k+1)]+P->stack[type][type_2]+fML[i+1][k-(i+1)];
1141                       if (cij == en+2*P->MLintern[1]+P->MLclosing) {
1142                         sector[s+2].ml = 2;
1143                         break;
1144                       }
1145                     }
1146                   }
1147                 }
1148                 break;
1149     } /* switch(dangles)... */
1150
1151     if (k<=j-3-TURN) { /* found the decomposition */
1152       sector[++s].i = i1;
1153       sector[s].j   = k;
1154       sector[++s].i = k+1;
1155       sector[s].j   = j1;
1156     } else {
1157 #if 0
1158       /* Y shaped ML loops fon't work yet */
1159       if (dangles==3) {
1160         /* (i,j) must close a Y shaped ML loop with coax stacking */
1161         if (cij ==  fML[i+1][j-2-(i+2)] + mm + d3 + d5 + P->MLbase + P->MLbase) {
1162           i1 = i+2;
1163           j1 = j-2;
1164         } else if (cij ==  fML[i+1][j-2-(i+1)] + mm + d5 + P->MLbase)
1165           j1 = j-2;
1166         else if (cij ==  fML[i+2][j-1-(i+2)] + mm + d3 + P->MLbase)
1167           i1 = i+2;
1168         else /* last chance */
1169           if (cij != fML[i+1][j-1-(i+1)] + mm + P->MLbase)
1170             fprintf(stderr,  "backtracking failed in repeat");
1171         /* if we arrive here we can express cij via fML[i1,j1]+dangles */
1172         sector[++s].i = i1;
1173         sector[s].j   = j1;
1174       }
1175       else
1176 #endif
1177         nrerror("backtracking failed in repeat");
1178     }
1179
1180     continue; /* this is a workarround to not accidentally proceed in the following block */
1181
1182   repeat_gquad:
1183     /*
1184       now we do some fancy stuff to backtrace the stacksize and linker lengths
1185       of the g-quadruplex that should reside within position i,j
1186     */
1187     {
1188       int l[3], L, a;
1189       L = -1;
1190
1191       get_gquad_pattern_mfe(S, i, j, P, &L, l);
1192       if(L != -1){
1193         /* fill the G's of the quadruplex into the structure string */
1194         for(a=0;a<L;a++){
1195           structure[i+a-start] = '+';
1196           structure[i+L+l[0]+a-start] = '+';
1197           structure[i+L+l[0]+L+l[1]+a-start] = '+';
1198           structure[i+L+l[0]+L+l[1]+L+l[2]+a-start] = '+';
1199         }
1200         goto repeat_gquad_exit;
1201       }
1202       nrerror("backtracking failed in repeat_gquad");
1203     }
1204   repeat_gquad_exit:
1205     asm("nop");
1206
1207   }
1208
1209   for (i=strlen(structure)-1; i>0 && structure[i] == '-'; i--)
1210     structure[i] = '\0';
1211   for (;i>=0; i--)
1212    if (structure[i]=='-') structure[i]='.';
1213
1214   return structure;
1215 }
1216
1217
1218 PRIVATE void update_fold_params(void){
1219   if(P) free(P);
1220   P = scale_parameters();
1221   make_pair_matrix();
1222 }
1223
1224 /*---------------------------------------------------------------------------*/
1225
1226 PRIVATE void make_ptypes(const short *S, int i, int maxdist, int n) {
1227   int j,k, type;
1228
1229   for (k=TURN+1; k<maxdist; k++) {
1230     j = i+k;
1231     if (j>n) continue;
1232     type = pair[S[i]][S[j]];
1233     if (noLonelyPairs && type) {
1234       if (!ptype[i+1][j-1-i-1])
1235         if (j==n || i==1 || (!pair[S[i-1]][S[j+1]])) type=0;
1236     }
1237     ptype[i][j-i]=type;
1238   }
1239 }