binaries/src/ViennaRNA/lib/part_func.c

   1 /*
   2                   partiton function for RNA secondary structures
   3
   4                   Ivo L Hofacker
   5                   Vienna RNA package
   6 */
   7 /*
   8   $Log: part_func.c,v $
   9   Revision 1.29  2008/02/23 10:10:49  ivo
  10   list returned from StackProb was sometimes not terminated correctly
  11
  12   Revision 1.28  2008/01/08 15:08:10  ivo
  13   circular fold would fail for open chain
  14
  15   Revision 1.27  2007/12/05 13:04:04  ivo
  16   add various circfold variants from Ronny
  17
  18   Revision 1.26  2007/09/19 12:41:56  ivo
  19   add computation of centroid() structure for RNAfold -p
  20
  21   Revision 1.25  2007/04/30 15:12:00  ivo
  22   merge RNAup into package
  23
  24   Revision 1.24  2007/03/03 17:57:44  ivo
  25   make sure entries in scale[] decrease to 0
  26
  27   Revision 1.23  2006/12/01 12:40:23  ivo
  28   undo Ulli's accidental commit
  29
  30   Revision 1.21  2006/08/04 15:39:06  ivo
  31   new function stackProb returns probability for stacks
  32   p[(i,j)(i+1,j-1)]
  33
  34   Revision 1.20  2004/08/12 12:14:46  ivo
  35   update
  36
  37   Revision 1.19  2004/05/14 16:28:05  ivo
  38   fix the bug in make_ptype here too (fixed in 1.27 of fold.c)
  39
  40   Revision 1.18  2004/02/17 10:46:52  ivo
  41   make sure init_pf_fold is called before scale_parameters
  42
  43   Revision 1.17  2004/02/09 18:37:59  ivo
  44   new mean_bp_dist() function to compute ensemble diversity
  45
  46   Revision 1.16  2003/08/04 09:14:09  ivo
  47   finish up stochastic backtracking
  48
  49   Revision 1.15  2002/03/19 16:51:12  ivo
  50   more on stochastic backtracking (still incomplete)
  51
  52   Revision 1.14  2002/02/08 17:37:23  ivo
  53   set freed S,S1 pointers to NULL
  54
  55   Revision 1.13  2001/11/16 17:30:04  ivo
  56   add stochastic backtracking (still incomplete)
  57 */
  58 #include <config.h>
  59 #include <stdio.h>
  60 #include <stdlib.h>
  61 #include <string.h>
  62 #include <math.h>
  63 #include <float.h>    /* #defines FLT_MAX ... */
  64 #include <limits.h>
  65
  66 #include "utils.h"
  67 #include "energy_par.h"
  68 #include "fold_vars.h"
  69 #include "pair_mat.h"
  70 #include "params.h"
  71 #include "loop_energies.h"
  72 #include "gquad.h"
  73 #include "part_func.h"
  74
  75 #ifdef _OPENMP
  76 #include <omp.h>
  77 #endif
  78
  79 #define ISOLATED  256.0
  80
  81 /*
  82 #################################
  83 # GLOBAL VARIABLES              #
  84 #################################
  85 */
  86 PUBLIC  int         st_back = 0;
  87
  88 /*
  89 #################################
  90 # PRIVATE VARIABLES             #
  91 #################################
  92 */
  93 PRIVATE FLT_OR_DBL  *q = NULL, *qb=NULL, *qm = NULL, *qm1 = NULL, *qqm = NULL, *qqm1 = NULL, *qq = NULL, *qq1 = NULL;
  94 PRIVATE FLT_OR_DBL  *probs=NULL, *prml=NULL, *prm_l=NULL, *prm_l1=NULL, *q1k=NULL, *qln=NULL;
  95 PRIVATE FLT_OR_DBL  *scale=NULL;
  96 PRIVATE FLT_OR_DBL  *expMLbase=NULL;
  97 PRIVATE FLT_OR_DBL  qo=0., qho=0., qio=0., qmo=0., *qm2=NULL;
  98 PRIVATE int         *jindx=NULL;
  99 PRIVATE int         *my_iindx=NULL;
 100 PRIVATE int         init_length = -1;   /* length in last call to init_pf_fold() */
 101 PRIVATE int         circular=0;
 102 PRIVATE int         do_bppm = 1;             /* do backtracking per default */
 103 PRIVATE int         struct_constrained = 0;
 104 PRIVATE char        *pstruc=NULL;
 105 PRIVATE char        *sequence=NULL;
 106 PRIVATE char        *ptype=NULL;        /* precomputed array of pair types */
 107 PRIVATE pf_paramT   *pf_params=NULL;    /* the precomputed Boltzmann weights */
 108 PRIVATE short       *S=NULL, *S1=NULL;
 109 PRIVATE int         with_gquad = 0;
 110
 111 PRIVATE FLT_OR_DBL  *G = NULL, *Gj = NULL, *Gj1 = NULL;
 112
 113 #ifdef _OPENMP
 114
 115 #pragma omp threadprivate(q, qb, qm, qm1, qqm, qqm1, qq, qq1, prml, prm_l, prm_l1, q1k, qln,\
 116                           probs, scale, expMLbase, qo, qho, qio, qmo, qm2, jindx, my_iindx, init_length,\
 117                           circular, pstruc, sequence, ptype, pf_params, S, S1, do_bppm, struct_constrained,\
 118                           G, Gj, Gj1, with_gquad)
 119
 120 #endif
 121
 122 /*
 123 #################################
 124 # PRIVATE FUNCTION DECLARATIONS #
 125 #################################
 126 */
 127 PRIVATE void  init_partfunc(int length, pf_paramT *parameters);
 128 PRIVATE void  scale_pf_params(unsigned int length, pf_paramT *parameters);
 129 PRIVATE void  get_arrays(unsigned int length);
 130 PRIVATE void  make_ptypes(const short *S, const char *structure);
 131 PRIVATE void  pf_circ(const char *sequence, char *structure);
 132 PRIVATE void  pf_linear(const char *sequence, char *structure);
 133 PRIVATE void  pf_create_bppm(const char *sequence, char *structure);
 134 PRIVATE void  backtrack(int i, int j);
 135 PRIVATE void  backtrack_qm(int i, int j);
 136 PRIVATE void  backtrack_qm1(int i,int j);
 137 PRIVATE void  backtrack_qm2(int u, int n);
 138
 139 /*
 140 #################################
 141 # BEGIN OF FUNCTION DEFINITIONS #
 142 #################################
 143 */
 144
 145 PRIVATE void init_partfunc(int length, pf_paramT *parameters){
 146   if (length<1) nrerror("init_pf_fold: length must be greater 0");
 147
 148 #ifdef _OPENMP
 149 /* Explicitly turn off dynamic threads */
 150   omp_set_dynamic(0);
 151   free_pf_arrays(); /* free previous allocation */
 152 #else
 153   if (init_length>0) free_pf_arrays(); /* free previous allocation */
 154 #endif
 155
 156 #ifdef SUN4
 157   nonstandard_arithmetic();
 158 #else
 159 #ifdef HP9
 160   fpsetfastmode(1);
 161 #endif
 162 #endif
 163   make_pair_matrix();
 164   get_arrays((unsigned) length);
 165   scale_pf_params((unsigned) length, parameters);
 166
 167   init_length = length;
 168 }
 169
 170 PRIVATE void get_arrays(unsigned int length){
 171   unsigned int size;
 172
 173   if((length +1) >= (unsigned int)sqrt((double)INT_MAX))
 174     nrerror("get_arrays@part_func.c: sequence length exceeds addressable range");
 175
 176   size  = sizeof(FLT_OR_DBL) * ((length+1)*(length+2)/2);
 177
 178   q     = (FLT_OR_DBL *) space(size);
 179   qb    = (FLT_OR_DBL *) space(size);
 180   qm    = (FLT_OR_DBL *) space(size);
 181   qm1   = (st_back || circular) ? (FLT_OR_DBL *) space(size) : NULL;
 182   qm2   = (circular) ? (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2)) : NULL;
 183   probs = (do_bppm) ? (FLT_OR_DBL *) space(size) : NULL;
 184
 185   ptype     = (char *) space(sizeof(char)*((length+1)*(length+2)/2));
 186   q1k       = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+1));
 187   qln       = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 188   qq        = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 189   qq1       = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 190   qqm       = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 191   qqm1      = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 192   prm_l     = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 193   prm_l1    = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 194   prml      = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 195   expMLbase = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+1));
 196   scale     = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+1));
 197
 198   Gj        = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 199   Gj1       = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(length+2));
 200
 201   my_iindx  = get_iindx(length);
 202   iindx     = get_iindx(length); /* for backward compatibility and Perl wrapper */
 203   jindx     = get_indx(length);
 204 }
 205
 206 /**
 207 *** Allocate memory for all matrices and other stuff
 208 **/
 209 PUBLIC void free_pf_arrays(void){
 210   if(q)         free(q);
 211   if(qb)        free(qb);
 212   if(qm)        free(qm);
 213   if(qm1)       free(qm1);
 214   if(qm2)       free(qm2);
 215   if(ptype)     free(ptype);
 216   if(qq)        free(qq);
 217   if(qq1)       free(qq1);
 218   if(qqm)       free(qqm);
 219   if(qqm1)      free(qqm1);
 220   if(q1k)       free(q1k);
 221   if(qln)       free(qln);
 222   if(probs)     free(probs);
 223   if(prm_l)     free(prm_l);
 224   if(prm_l1)    free(prm_l1);
 225   if(prml)      free(prml);
 226   if(expMLbase) free(expMLbase);
 227   if(scale)     free(scale);
 228   if(my_iindx)  free(my_iindx);
 229   if(iindx)     free(iindx); /* for backward compatibility and Perl wrapper */
 230   if(jindx)     free(jindx);
 231   if(S)         free(S);
 232   if(S1)        free(S1);
 233   if(G)         free(G);
 234   if(Gj)        free(Gj);
 235   if(Gj1)       free(Gj1);
 236
 237   S = S1 = NULL;
 238   q = pr = probs = qb = qm = qm1 = qm2 = qq = qq1 = qqm = qqm1 = q1k = qln = prm_l = prm_l1 = prml = expMLbase = scale = G = Gj = Gj1 = NULL;
 239   my_iindx = jindx = iindx = NULL;
 240
 241   ptype = NULL;
 242
 243 #ifdef SUN4
 244   standard_arithmetic();
 245 #else
 246 #ifdef HP9
 247   fpsetfastmode(0);
 248 #endif
 249 #endif
 250
 251   init_length = 0;
 252 }
 253
 254 /*-----------------------------------------------------------------*/
 255 PUBLIC float pf_fold(const char *sequence, char *structure){
 256   return pf_fold_par(sequence, structure, NULL, do_backtrack, fold_constrained, 0);
 257 }
 258
 259 PUBLIC float pf_circ_fold(const char *sequence, char *structure){
 260   return pf_fold_par(sequence, structure, NULL, do_backtrack, fold_constrained, 1);
 261 }
 262
 263 PUBLIC float pf_fold_par( const char *sequence,
 264                           char *structure,
 265                           pf_paramT *parameters,
 266                           int calculate_bppm,
 267                           int is_constrained,
 268                           int is_circular){
 269
 270   FLT_OR_DBL  Q;
 271   double      free_energy;
 272   int         n = (int) strlen(sequence);
 273
 274   circular            = is_circular;
 275   do_bppm             = calculate_bppm;
 276   struct_constrained  = is_constrained;
 277
 278 #ifdef _OPENMP
 279   init_partfunc(n, parameters);
 280 #else
 281   if(parameters) init_partfunc(n, parameters);
 282   else if (n > init_length) init_partfunc(n, parameters);
 283   else if (fabs(pf_params->temperature - temperature)>1e-6) update_pf_params_par(n, parameters);
 284 #endif
 285
 286   with_gquad  = pf_params->model_details.gquad;
 287   S           = encode_sequence(sequence, 0);
 288   S1          = encode_sequence(sequence, 1);
 289
 290   make_ptypes(S, structure);
 291
 292   /* do the linear pf fold and fill all matrices  */
 293   pf_linear(sequence, structure);
 294
 295   if(circular)
 296     pf_circ(sequence, structure); /* do post processing step for circular RNAs */
 297
 298   if (backtrack_type=='C')      Q = qb[my_iindx[1]-n];
 299   else if (backtrack_type=='M') Q = qm[my_iindx[1]-n];
 300   else Q = (circular) ? qo : q[my_iindx[1]-n];
 301
 302   /* ensemble free energy in Kcal/mol              */
 303   if (Q<=FLT_MIN) fprintf(stderr, "pf_scale too large\n");
 304   free_energy = (-log(Q)-n*log(pf_params->pf_scale))*pf_params->kT/1000.0;
 305   /* in case we abort because of floating point errors */
 306   if (n>1600) fprintf(stderr, "free energy = %8.2f\n", free_energy);
 307
 308   /* calculate base pairing probability matrix (bppm)  */
 309   if(do_bppm){
 310     pf_create_bppm(sequence, structure);
 311     /*
 312     *  Backward compatibility:
 313     *  This block may be removed if deprecated functions
 314     *  relying on the global variable "pr" vanish from within the package!
 315     */
 316     pr = probs;
 317     /*
 318      {
 319       if(pr) free(pr);
 320       pr = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2));
 321       memcpy(pr, probs, sizeof(FLT_OR_DBL) * ((n+1)*(n+2)/2));
 322     }
 323     */
 324   }
 325   return free_energy;
 326 }
 327
 328 PRIVATE void pf_linear(const char *sequence, char *structure){
 329
 330   int n, i,j,k,l, ij, u,u1,d,ii, type, type_2, tt, minl, maxl;
 331
 332   int noGUclosure;
 333   FLT_OR_DBL expMLstem = 0.;
 334
 335   FLT_OR_DBL temp, Qmax=0;
 336   FLT_OR_DBL qbt1, *tmp;
 337
 338   FLT_OR_DBL  expMLclosing = pf_params->expMLclosing;
 339   double      max_real;
 340
 341   max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
 342
 343   n = (int) strlen(sequence);
 344
 345
 346   noGUclosure = pf_params->model_details.noGUclosure;
 347
 348   /*array initialization ; qb,qm,q
 349     qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */
 350
 351   if(with_gquad){
 352     expMLstem = exp_E_MLstem(0, -1, -1, pf_params);
 353     G         = get_gquad_pf_matrix(S, scale, pf_params);
 354   }
 355
 356   for (d=0; d<=TURN; d++)
 357     for (i=1; i<=n-d; i++) {
 358       j=i+d;
 359       ij = my_iindx[i]-j;
 360       q[ij]=1.0*scale[d+1];
 361       qb[ij]=qm[ij]=0.0;
 362     }
 363
 364   for (i=1; i<=n; i++)
 365     qq[i]=qq1[i]=qqm[i]=qqm1[i]=0;
 366
 367   for (j=TURN+2;j<=n; j++) {
 368     for (i=j-TURN-1; i>=1; i--) {
 369       /* construction of partition function of segment i,j*/
 370       /*firstly that given i binds j : qb(i,j) */
 371       u = j-i-1; ij = my_iindx[i]-j;
 372       type = ptype[ij];
 373       if (type!=0) {
 374         /*hairpin contribution*/
 375         if (((type==3)||(type==4))&&noGUclosure) qbt1 = 0;
 376         else
 377           qbt1 = exp_E_Hairpin(u, type, S1[i+1], S1[j-1], sequence+i-1, pf_params) * scale[u+2];
 378         /* interior loops with interior pair k,l */
 379         for (k=i+1; k<=MIN2(i+MAXLOOP+1,j-TURN-2); k++) {
 380           u1 = k-i-1;
 381           for (l=MAX2(k+TURN+1,j-1-MAXLOOP+u1); l<j; l++) {
 382             type_2 = ptype[my_iindx[k]-l];
 383             if (type_2) {
 384               type_2 = rtype[type_2];
 385               qbt1 += qb[my_iindx[k]-l] * (scale[u1+j-l+1] *
 386                                         exp_E_IntLoop(u1, j-l-1, type, type_2,
 387                                         S1[i+1], S1[j-1], S1[k-1], S1[l+1], pf_params));
 388             }
 389           }
 390         }
 391         /*multiple stem loop contribution*/
 392         ii = my_iindx[i+1]; /* ii-k=[i+1,k-1] */
 393         temp = 0.0;
 394         for (k=i+2; k<=j-1; k++) temp += qm[ii-(k-1)]*qqm1[k];
 395         tt = rtype[type];
 396         qbt1 += temp * expMLclosing * exp_E_MLstem(tt, S1[j-1], S1[i+1], pf_params) * scale[2];
 397
 398         if(with_gquad){
 399           qbt1 += exp_E_GQuad_IntLoop(i, j, type, S1, G, my_iindx, pf_params) * scale[2];
 400         }
 401
 402         qb[ij] = qbt1;
 403       }
 404       /* end if (type!=0) */
 405       else
 406         qb[ij] = 0.0;
 407
 408       /* construction of qqm matrix containing final stem
 409          contributions to multiple loop partition function
 410          from segment i,j */
 411       qqm[i] = qqm1[i]*expMLbase[1];
 412       if (type) {
 413         qbt1 = qb[ij] * exp_E_MLstem(type, ((i>1) || circular) ? S1[i-1] : -1, ((j<n) || circular) ? S1[j+1] : -1, pf_params);
 414         qqm[i] += qbt1;
 415       }
 416
 417       if(with_gquad){
 418         /*include gquads into qqm*/
 419         qqm[i] += G[my_iindx[i]-j] * expMLstem;
 420       }
 421
 422       if (qm1) qm1[jindx[j]+i] = qqm[i]; /* for stochastic backtracking and circfold */
 423
 424       /*construction of qm matrix containing multiple loop
 425         partition function contributions from segment i,j */
 426       temp = 0.0;
 427       ii = my_iindx[i];  /* ii-k=[i,k-1] */
 428       for (k=j; k>i; k--) temp += (qm[ii-(k-1)] + expMLbase[k-i])*qqm[k];
 429       qm[ij] = (temp + qqm[i]);
 430
 431       /*auxiliary matrix qq for cubic order q calculation below */
 432       qbt1=0.0;
 433       if (type){
 434         qbt1 += qb[ij];
 435         qbt1 *= exp_E_ExtLoop(type, ((i>1) || circular) ? S1[i-1] : -1, ((j<n) || circular) ? S1[j+1] : -1, pf_params);
 436       }
 437
 438       if(with_gquad){
 439         qbt1 += G[ij];
 440       }
 441
 442       qq[i] = qq1[i]*scale[1] + qbt1;
 443
 444       /*construction of partition function for segment i,j */
 445       temp = 1.0*scale[1+j-i] + qq[i];
 446       for (k=i; k<=j-1; k++) temp += q[ii-k]*qq[k+1];
 447       q[ij] = temp;
 448       if (temp>Qmax) {
 449         Qmax = temp;
 450         if (Qmax>max_real/10.)
 451           fprintf(stderr, "Q close to overflow: %d %d %g\n", i,j,temp);
 452       }
 453       if (temp>=max_real) {
 454         PRIVATE char msg[128];
 455         sprintf(msg, "overflow in pf_fold while calculating q[%d,%d]\n"
 456                      "use larger pf_scale", i,j);
 457         nrerror(msg);
 458       }
 459     }
 460     tmp = qq1;  qq1 =qq;  qq =tmp;
 461     tmp = qqm1; qqm1=qqm; qqm=tmp;
 462
 463     if(with_gquad){ /* rotate the auxilary g-quadruplex matrices */
 464       tmp = Gj1; Gj1=Gj; Gj=tmp;
 465     }
 466   }
 467 }
 468
 469 /* calculate partition function for circular case */
 470 /* NOTE: this is the postprocessing step ONLY     */
 471 /* You have to call pf_linear first to calculate  */
 472 /* complete circular case!!!                      */
 473 PRIVATE void pf_circ(const char *sequence, char *structure){
 474
 475   int u, p, q, k, l;
 476   int noGUclosure;
 477   int n = (int) strlen(sequence);
 478
 479   FLT_OR_DBL  qot;
 480   FLT_OR_DBL  expMLclosing  = pf_params->expMLclosing;
 481
 482   noGUclosure = pf_params->model_details.noGUclosure;
 483   qo = qho = qio = qmo = 0.;
 484
 485   /* construct qm2 matrix from qm1 entries  */
 486   for(k=1; k<n-TURN-1; k++){
 487     qot = 0.;
 488     for (u=k+TURN+1; u<n-TURN-1; u++)
 489       qot += qm1[jindx[u]+k]*qm1[jindx[n]+(u+1)];
 490     qm2[k] = qot;
 491    }
 492
 493   for(p = 1; p < n; p++){
 494     for(q = p + TURN + 1; q <= n; q++){
 495       int type;
 496       /* 1. get exterior hairpin contribution  */
 497       u = n-q + p-1;
 498       if (u<TURN) continue;
 499       type = ptype[my_iindx[p]-q];
 500       if (!type) continue;
 501        /* cause we want to calc the exterior loops, we need the reversed pair type from now on  */
 502       type=rtype[type];
 503
 504       char loopseq[10];
 505       if (u<7){
 506         strcpy(loopseq , sequence+q-1);
 507         strncat(loopseq, sequence, p);
 508       }
 509       qho += (((type==3)||(type==4))&&noGUclosure) ? 0. : qb[my_iindx[p]-q] * exp_E_Hairpin(u, type, S1[q+1], S1[p-1],  loopseq, pf_params) * scale[u];
 510
 511       /* 2. exterior interior loops, i "define" the (k,l) pair as "outer pair"  */
 512       /* so "outer type" is rtype[type[k,l]] and inner type is type[p,q]        */
 513       qot = 0.;
 514       for(k=q+1; k < n; k++){
 515         int ln1, lstart;
 516         ln1 = k - q - 1;
 517         if(ln1+p-1>MAXLOOP) break;
 518         lstart = ln1+p-1+n-MAXLOOP;
 519         if(lstart<k+TURN+1) lstart = k + TURN + 1;
 520         for(l=lstart;l <= n; l++){
 521           int ln2, type2;
 522           ln2 = (p - 1) + (n - l);
 523
 524           if((ln1+ln2) > MAXLOOP) continue;
 525
 526           type2 = ptype[my_iindx[k]-l];
 527           if(!type2) continue;
 528           qio += qb[my_iindx[p]-q] * qb[my_iindx[k]-l] * exp_E_IntLoop(ln2, ln1, rtype[type2], type, S1[l+1], S1[k-1], S1[p-1], S1[q+1], pf_params) * scale[ln1+ln2];
 529         }
 530       } /* end of kl double loop */
 531     }
 532   } /* end of pq double loop */
 533
 534   /* 3. Multiloops  */
 535   for(k=TURN+2; k<n-2*TURN-3; k++)
 536     qmo += qm[my_iindx[1]-k] * qm2[k+1] * expMLclosing;
 537
 538   /* add an additional pf of 1.0 to take the open chain into account too           */
 539   qo = qho + qio + qmo + 1.0*scale[n];
 540 }
 541
 542 /* calculate base pairing probs */
 543 PUBLIC void pf_create_bppm(const char *sequence, char *structure){
 544   int n, i,j,k,l, ij, kl, ii, i1, ll, type, type_2, tt, u1, ov=0;
 545   FLT_OR_DBL  temp, Qmax=0, prm_MLb;
 546   FLT_OR_DBL  prmt,prmt1;
 547   FLT_OR_DBL  *tmp;
 548   FLT_OR_DBL  tmp2;
 549   FLT_OR_DBL  expMLclosing = pf_params->expMLclosing;
 550   double      max_real;
 551
 552   FLT_OR_DBL  expMLstem = (with_gquad) ? exp_E_MLstem(0, -1, -1, pf_params) : 0;
 553
 554   max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
 555
 556   if((S != NULL) && (S1 != NULL)){
 557     n = S[0];
 558     Qmax=0;
 559
 560     for (k=1; k<=n; k++) {
 561       q1k[k] = q[my_iindx[1] - k];
 562       qln[k] = q[my_iindx[k] -n];
 563     }
 564     q1k[0] = 1.0;
 565     qln[n+1] = 1.0;
 566
 567 /*  pr = q; */     /* recycling */
 568
 569
 570     /* 1. exterior pair i,j and initialization of pr array */
 571     if(circular){
 572       for (i=1; i<=n; i++) {
 573         for (j=i; j<=MIN2(i+TURN,n); j++)
 574           probs[my_iindx[i]-j] = 0;
 575         for (j=i+TURN+1; j<=n; j++) {
 576           ij = my_iindx[i]-j;
 577           type = ptype[ij];
 578           if (type&&(qb[ij]>0.)) {
 579             probs[ij] = 1./qo;
 580             int rt = rtype[type];
 581
 582             /* 1.1. Exterior Hairpin Contribution */
 583             int u = i + n - j -1;
 584             /* get the loop sequence */
 585             char loopseq[10];
 586             if (u<7){
 587               strcpy(loopseq , sequence+j-1);
 588               strncat(loopseq, sequence, i);
 589             }
 590             tmp2 = exp_E_Hairpin(u, rt, S1[j+1], S1[i-1], loopseq, pf_params) * scale[u];
 591
 592             /* 1.2. Exterior Interior Loop Contribution                    */
 593             /* 1.2.1. i,j  delimtis the "left" part of the interior loop    */
 594             /* (j,i) is "outer pair"                                                */
 595             for(k=1; k < i-TURN-1; k++){
 596               int ln1, lstart;
 597               ln1 = k + n - j - 1;
 598               if(ln1>MAXLOOP) break;
 599               lstart = ln1+i-1-MAXLOOP;
 600               if(lstart<k+TURN+1) lstart = k + TURN + 1;
 601               for(l=lstart; l < i; l++){
 602                 int ln2, type_2;
 603                 type_2 = ptype[my_iindx[k]-l];
 604                 if (type_2==0) continue;
 605                 ln2 = i - l - 1;
 606                 if(ln1+ln2>MAXLOOP) continue;
 607                 tmp2 += qb[my_iindx[k] - l]
 608                         * exp_E_IntLoop(ln1,
 609                                         ln2,
 610                                         rt,
 611                                         rtype[type_2],
 612                                         S1[j+1],
 613                                         S1[i-1],
 614                                         S1[k-1],
 615                                         S1[l+1],
 616                                         pf_params)
 617                         * scale[ln1 + ln2];
 618               }
 619             }
 620             /* 1.2.2. i,j  delimtis the "right" part of the interior loop  */
 621             for(k=j+1; k < n-TURN; k++){
 622               int ln1, lstart;
 623               ln1 = k - j - 1;
 624               if((ln1 + i - 1)>MAXLOOP) break;
 625               lstart = ln1+i-1+n-MAXLOOP;
 626               if(lstart<k+TURN+1) lstart = k + TURN + 1;
 627               for(l=lstart; l <= n; l++){
 628                 int ln2, type_2;
 629                 type_2 = ptype[my_iindx[k]-l];
 630                 if (type_2==0) continue;
 631                 ln2 = i - 1 + n - l;
 632                 if(ln1+ln2>MAXLOOP) continue;
 633                 tmp2 += qb[my_iindx[k] - l]
 634                         * exp_E_IntLoop(ln2,
 635                                         ln1,
 636                                         rtype[type_2],
 637                                         rt,
 638                                         S1[l+1],
 639                                         S1[k-1],
 640                                         S1[i-1],
 641                                         S1[j+1],
 642                                         pf_params)
 643                         * scale[ln1 + ln2];
 644               }
 645             }
 646             /* 1.3 Exterior multiloop decomposition */
 647             /* 1.3.1 Middle part                    */
 648             if((i>TURN+2) && (j<n-TURN-1))
 649               tmp2 += qm[my_iindx[1]-i+1]
 650                       * qm[my_iindx[j+1]-n]
 651                       * expMLclosing
 652                       * exp_E_MLstem(type, S1[i-1], S1[j+1], pf_params);
 653
 654             /* 1.3.2 Left part                      */
 655             for(k=TURN+2; k < i-TURN-2; k++)
 656               tmp2 += qm[my_iindx[1]-k]
 657                       * qm1[jindx[i-1]+k+1]
 658                       * expMLbase[n-j]
 659                       * expMLclosing
 660                       * exp_E_MLstem(type, S1[i-1], S1[j+1], pf_params);
 661
 662             /* 1.3.3 Right part                      */
 663             for(k=j+TURN+2; k < n-TURN-1;k++)
 664               tmp2 += qm[my_iindx[j+1]-k]
 665                       * qm1[jindx[n]+k+1]
 666                       * expMLbase[i-1]
 667                       * expMLclosing
 668                       * exp_E_MLstem(type, S1[i-1], S1[j+1], pf_params);
 669
 670             /* all exterior loop decompositions for pair i,j done  */
 671             probs[ij] *= tmp2;
 672
 673           }
 674           else probs[ij] = 0;
 675         }
 676       }
 677     } /* end if(circular)  */
 678     else {
 679       for (i=1; i<=n; i++) {
 680         for (j=i; j<=MIN2(i+TURN,n); j++) probs[my_iindx[i]-j] = 0;
 681         for (j=i+TURN+1; j<=n; j++) {
 682           ij = my_iindx[i]-j;
 683           type = ptype[ij];
 684           if (type&&(qb[ij]>0.)) {
 685             probs[ij] = q1k[i-1]*qln[j+1]/q1k[n];
 686             probs[ij] *= exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j<n) ? S1[j+1] : -1, pf_params);
 687           }
 688           else
 689             probs[ij] = 0.;
 690         }
 691       }
 692     } /* end if(!circular)  */
 693
 694     for (l=n; l>TURN+1; l--) {
 695
 696       /* 2. bonding k,l as substem of 2:loop enclosed by i,j */
 697       for (k=1; k<l-TURN; k++) {
 698         kl = my_iindx[k]-l;
 699         type_2 = ptype[kl];
 700         if (type_2==0) continue;
 701         type_2 = rtype[type_2];
 702         if (qb[kl]==0.) continue;
 703
 704         tmp2 = 0.;
 705         for (i=MAX2(1,k-MAXLOOP-1); i<=k-1; i++)
 706           for (j=l+1; j<=MIN2(l+ MAXLOOP -k+i+2,n); j++) {
 707             ij = my_iindx[i] - j;
 708             type = ptype[ij];
 709             if (type && (probs[ij]>0.)) {
 710               /* add *scale[u1+u2+2] */
 711               tmp2 +=  probs[ij]
 712                        * (scale[k-i+j-l]
 713                        * exp_E_IntLoop(k - i - 1,
 714                                        j - l - 1,
 715                                        type,
 716                                        type_2,
 717                                        S1[i + 1],
 718                                        S1[j - 1],
 719                                        S1[k - 1],
 720                                        S1[l + 1],
 721                                        pf_params));
 722             }
 723           }
 724         probs[kl] += tmp2;
 725       }
 726
 727       if(with_gquad){
 728         /* 2.5. bonding k,l as gquad enclosed by i,j */
 729         FLT_OR_DBL *expintern = &(pf_params->expinternal[0]);
 730
 731         if(l < n - 3){
 732           for(k = 2; k <= l - VRNA_GQUAD_MIN_BOX_SIZE; k++){
 733             kl = my_iindx[k]-l;
 734             if (G[kl]==0.) continue;
 735             tmp2 = 0.;
 736             i = k - 1;
 737             for(j = MIN2(l + MAXLOOP + 1, n); j > l + 3; j--){
 738               ij = my_iindx[i] - j;
 739               type = ptype[ij];
 740               if(!type) continue;
 741               tmp2 += probs[ij] * expintern[j-l-1] * pf_params->expmismatchI[type][S1[i+1]][S1[j-1]] * scale[2];
 742             }
 743             probs[kl] += tmp2 * G[kl];
 744           }
 745         }
 746
 747         if(l < n){
 748           for(k = 4; k <= l - VRNA_GQUAD_MIN_BOX_SIZE; k++){
 749             kl = my_iindx[k]-l;
 750             if (G[kl]==0.) continue;
 751             tmp2 = 0.;
 752             j = l + 1;
 753             for (i=MAX2(1,k-MAXLOOP-1); i < k - 3; i++){
 754               ij = my_iindx[i] - j;
 755               type = ptype[ij];
 756               if(!type) continue;
 757               tmp2 += probs[ij] * expintern[k - i - 1] * pf_params->expmismatchI[type][S1[i+1]][S1[j-1]] * scale[2];
 758             }
 759             probs[kl] += tmp2 * G[kl];
 760           }
 761         }
 762
 763         if (l < n - 1){
 764           for (k=3; k<=l-VRNA_GQUAD_MIN_BOX_SIZE; k++) {
 765             kl = my_iindx[k]-l;
 766             if (G[kl]==0.) continue;
 767             tmp2 = 0.;
 768             for (i=MAX2(1,k-MAXLOOP-1); i<=k-2; i++){
 769               u1 = k - i - 1;
 770               for (j=l+2; j<=MIN2(l + MAXLOOP - u1 + 1,n); j++) {
 771                 ij = my_iindx[i] - j;
 772                 type = ptype[ij];
 773                 if(!type) continue;
 774                 tmp2 += probs[ij] * expintern[u1+j-l-1] * pf_params->expmismatchI[type][S1[i+1]][S1[j-1]] * scale[2];
 775               }
 776             }
 777             probs[kl] += tmp2 * G[kl];
 778           }
 779         }
 780       }
 781
 782       /* 3. bonding k,l as substem of multi-loop enclosed by i,j */
 783       prm_MLb = 0.;
 784       if (l<n) for (k=2; k<l-TURN; k++) {
 785         i = k-1;
 786         prmt = prmt1 = 0.0;
 787
 788         ii = my_iindx[i];     /* ii-j=[i,j]     */
 789         ll = my_iindx[l+1];   /* ll-j=[l+1,j-1] */
 790         tt = ptype[ii-(l+1)]; tt=rtype[tt];
 791         /* (i, l+1) closes the ML with substem (k,l) */
 792         if(tt)
 793           prmt1 = probs[ii-(l+1)] * expMLclosing * exp_E_MLstem(tt, S1[l], S1[i+1], pf_params);
 794
 795         /* (i,j) with j>l+1 closes the ML with substem (k,l) */
 796         for (j=l+2; j<=n; j++) {
 797           tt = ptype[ii-j]; tt = rtype[tt];
 798           if(tt)
 799             prmt += probs[ii-j] * exp_E_MLstem(tt, S1[j-1], S1[i+1], pf_params) * qm[ll-(j-1)];
 800         }
 801         kl = my_iindx[k]-l;
 802         tt = ptype[kl];
 803         prmt *= expMLclosing;
 804         prml[ i] = prmt;
 805         prm_l[i] = prm_l1[i]*expMLbase[1]+prmt1;
 806
 807         prm_MLb = prm_MLb*expMLbase[1] + prml[i];
 808         /* same as:    prm_MLb = 0;
 809            for (i=1; i<=k-1; i++) prm_MLb += prml[i]*expMLbase[k-i-1]; */
 810
 811         prml[i] = prml[ i] + prm_l[i];
 812
 813         if(with_gquad){
 814           if ((!tt) && (G[kl] == 0.)) continue;
 815         } else {
 816           if (qb[kl] == 0.) continue;
 817         }
 818
 819         temp = prm_MLb;
 820
 821         for (i=1;i<=k-2; i++)
 822           temp += prml[i]*qm[my_iindx[i+1] - (k-1)];
 823
 824         if(with_gquad){
 825           if(tt)
 826             temp    *= exp_E_MLstem(tt, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params) * scale[2];
 827           else
 828             temp    *= G[kl] * expMLstem * scale[2];
 829         } else {
 830           temp    *= exp_E_MLstem(tt, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params) * scale[2];
 831         }
 832
 833         probs[kl]  += temp;
 834
 835         if (probs[kl]>Qmax) {
 836           Qmax = probs[kl];
 837           if (Qmax>max_real/10.)
 838             fprintf(stderr, "P close to overflow: %d %d %g %g\n",
 839               i, j, probs[kl], qb[kl]);
 840         }
 841         if (probs[kl]>=max_real) {
 842           ov++;
 843           probs[kl]=FLT_MAX;
 844         }
 845
 846       } /* end for (k=..) */
 847       tmp = prm_l1; prm_l1=prm_l; prm_l=tmp;
 848
 849     }  /* end for (l=..)   */
 850
 851     for (i=1; i<=n; i++)
 852       for (j=i+TURN+1; j<=n; j++) {
 853         ij = my_iindx[i]-j;
 854
 855         if(with_gquad){
 856           if (qb[ij] > 0.)
 857             probs[ij] *= qb[ij];
 858           if (G[ij] > 0.){
 859             probs[ij] += q1k[i-1] * G[ij] * qln[j+1]/q1k[n];
 860           }
 861         } else {
 862           if (qb[ij] > 0.)
 863             probs[ij] *= qb[ij];
 864         }
 865       }
 866
 867     if (structure!=NULL)
 868       bppm_to_structure(structure, probs, n);
 869     if (ov>0) fprintf(stderr, "%d overflows occurred while backtracking;\n"
 870         "you might try a smaller pf_scale than %g\n",
 871         ov, pf_params->pf_scale);
 872   } /* end if((S != NULL) && (S1 != NULL))  */
 873   else
 874     nrerror("bppm calculations have to be done after calling forward recursion\n");
 875   return;
 876 }
 877
 878 PRIVATE void scale_pf_params(unsigned int length, pf_paramT *parameters){
 879   unsigned int  i;
 880   double        scaling_factor;
 881
 882   if(pf_params) free(pf_params);
 883
 884   if(parameters){
 885     pf_params = get_boltzmann_factor_copy(parameters);
 886   } else {
 887     model_detailsT md;
 888     set_model_details(&md);
 889     pf_params = get_boltzmann_factors(temperature, 1.0, md, pf_scale);
 890   }
 891
 892   scaling_factor = pf_params->pf_scale;
 893
 894   /* scaling factors (to avoid overflows) */
 895   if (scaling_factor == -1) { /* mean energy for random sequences: 184.3*length cal */
 896     scaling_factor = exp(-(-185+(pf_params->temperature-37.)*7.27)/pf_params->kT);
 897     if (scaling_factor<1) scaling_factor=1;
 898     pf_params->pf_scale = scaling_factor;
 899     pf_scale = pf_params->pf_scale; /* compatibility with RNAup, may be removed sometime */
 900   }
 901   scale[0] = 1.;
 902   scale[1] = 1./scaling_factor;
 903   expMLbase[0] = 1;
 904   expMLbase[1] = pf_params->expMLbase/scaling_factor;
 905   for (i=2; i<=length; i++) {
 906     scale[i] = scale[i/2]*scale[i-(i/2)];
 907     expMLbase[i] = pow(pf_params->expMLbase, (double)i) * scale[i];
 908   }
 909 }
 910
 911 /*---------------------------------------------------------------------------*/
 912
 913 PUBLIC void update_pf_params(int length){
 914   update_pf_params_par(length, NULL);
 915 }
 916
 917 PUBLIC void update_pf_params_par(int length, pf_paramT *parameters){
 918 #ifdef _OPENMP
 919   make_pair_matrix(); /* is this really necessary? */
 920   scale_pf_params((unsigned) length, parameters);
 921 #else
 922   if(parameters) init_partfunc(length, parameters);
 923   else if (length>init_length) init_partfunc(length, parameters);  /* init not update */
 924   else {
 925     make_pair_matrix();
 926     scale_pf_params((unsigned) length, parameters);
 927   }
 928 #endif
 929 }
 930
 931 /*---------------------------------------------------------------------------*/
 932
 933 PUBLIC char bppm_symbol(const float *x){
 934 /*  if( ((x[1]-x[2])*(x[1]-x[2]))<0.1&&x[0]<=0.677) return '|'; */
 935   if( x[0] > 0.667 )  return '.';
 936   if( x[1] > 0.667 )  return '(';
 937   if( x[2] > 0.667 )  return ')';
 938   if( (x[1]+x[2]) > x[0] ) {
 939     if( (x[1]/(x[1]+x[2])) > 0.667) return '{';
 940     if( (x[2]/(x[1]+x[2])) > 0.667) return '}';
 941     else return '|';
 942   }
 943   if( x[0] > (x[1]+x[2]) ) return ',';
 944   return ':';
 945 }
 946
 947 PUBLIC void bppm_to_structure(char *structure, FLT_OR_DBL *p, unsigned int length){
 948   int    i, j;
 949   int   *index = get_iindx(length);
 950   float  P[3];   /* P[][0] unpaired, P[][1] upstream p, P[][2] downstream p */
 951
 952   for( j=1; j<=length; j++ ) {
 953     P[0] = 1.0;
 954     P[1] = P[2] = 0.0;
 955     for( i=1; i<j; i++) {
 956       P[2] += p[index[i]-j];    /* j is paired downstream */
 957       P[0] -= p[index[i]-j];    /* j is unpaired */
 958     }
 959     for( i=j+1; i<=length; i++ ) {
 960       P[1] += p[index[j]-i];    /* j is paired upstream */
 961       P[0] -= p[index[j]-i];    /* j is unpaired */
 962     }
 963     structure[j-1] = bppm_symbol(P);
 964   }
 965   structure[length] = '\0';
 966   free(index);
 967 }
 968
 969
 970 /*---------------------------------------------------------------------------*/
 971 PRIVATE void make_ptypes(const short *S, const char *structure){
 972   int n,i,j,k,l, noLP;
 973
 974   noLP = pf_params->model_details.noLP;
 975
 976   n=S[0];
 977   for (k=1; k<n-TURN; k++)
 978     for (l=1; l<=2; l++) {
 979       int type,ntype=0,otype=0;
 980       i=k; j = i+TURN+l; if (j>n) continue;
 981       type = pair[S[i]][S[j]];
 982       while ((i>=1)&&(j<=n)) {
 983         if ((i>1)&&(j<n)) ntype = pair[S[i-1]][S[j+1]];
 984         if (noLP && (!otype) && (!ntype))
 985           type = 0; /* i.j can only form isolated pairs */
 986         qb[my_iindx[i]-j] = 0.;
 987         ptype[my_iindx[i]-j] = (char) type;
 988         otype =  type;
 989         type  = ntype;
 990         i--; j++;
 991       }
 992     }
 993
 994   if (struct_constrained && (structure != NULL))
 995     constrain_ptypes(structure, (unsigned int)n, ptype, NULL, TURN, 1);
 996 }
 997
 998 /*
 999   stochastic backtracking in pf_fold arrays
1000   returns random structure S with Boltzman probabilty
1001   p(S) = exp(-E(S)/kT)/Z
1002 */
1003 char *pbacktrack(char *seq){
1004   double r, qt;
1005   int i,j,n, start;
1006   sequence = seq;
1007   n = strlen(sequence);
1008
1009   if (init_length<1)
1010     nrerror("can't backtrack without pf arrays.\n"
1011             "Call pf_fold() before pbacktrack()");
1012   pstruc = space((n+1)*sizeof(char));
1013
1014   for (i=0; i<n; i++) pstruc[i] = '.';
1015
1016   start = 1;
1017   while (start<n) {
1018   /* find i position of first pair */
1019     for (i=start; i<n; i++) {
1020       r = urn() * qln[i];
1021       if (r > qln[i+1]*scale[1])  break; /* i is paired */
1022     }
1023     if (i>=n) break; /* no more pairs */
1024     /* now find the pairing partner j */
1025     r = urn() * (qln[i] - qln[i+1]*scale[1]);
1026     for (qt=0, j=i+1; j<=n; j++) {
1027       int type;
1028       type = ptype[my_iindx[i]-j];
1029       if (type) {
1030         double qkl;
1031         qkl = qb[my_iindx[i]-j];
1032         if (j<n) qkl *= qln[j+1];
1033         qkl *=  exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j<n) ? S1[j+1] : -1, pf_params);
1034         qt += qkl;
1035         if (qt > r) break; /* j is paired */
1036       }
1037     }
1038     if (j==n+1) nrerror("backtracking failed in ext loop");
1039     start = j+1;
1040     backtrack(i,j);
1041   }
1042
1043   return pstruc;
1044 }
1045 char *pbacktrack_circ(char *seq){
1046   double r, qt;
1047   int i, j, k, l, n;
1048   FLT_OR_DBL  expMLclosing      = pf_params->expMLclosing;
1049
1050   sequence = seq;
1051   n = strlen(sequence);
1052
1053   if (init_length<1)
1054     nrerror("can't backtrack without pf arrays.\n"
1055       "Call pf_circ_fold() before pbacktrack_circ()");
1056   pstruc = space((n+1)*sizeof(char));
1057
1058   /* initialize pstruct with single bases  */
1059   for (i=0; i<n; i++) pstruc[i] = '.';
1060
1061   qt = 1.0*scale[n];
1062   r = urn() * qo;
1063
1064   /* open chain? */
1065   if(qt > r) return pstruc;
1066
1067   for(i=1; (i < n); i++){
1068     for(j=i+TURN+1;(j<=n); j++){
1069
1070       int type, u;
1071       /* 1. first check, wether we can do a hairpin loop  */
1072       u = n-j + i-1;
1073       if (u<TURN) continue;
1074
1075       type = ptype[my_iindx[i]-j];
1076       if (!type) continue;
1077
1078       type=rtype[type];
1079
1080       char loopseq[10];
1081       if (u<7){
1082         strcpy(loopseq , sequence+j-1);
1083         strncat(loopseq, sequence, i);
1084       }
1085
1086       qt += qb[my_iindx[i]-j] * exp_E_Hairpin(u, type, S1[j+1], S1[i-1],  loopseq, pf_params) * scale[u];
1087       /* found a hairpin? so backtrack in the enclosed part and we're done  */
1088       if(qt>r){ backtrack(i,j); return pstruc;}
1089
1090       /* 2. search for (k,l) with which we can close an interior loop  */
1091       for(k=j+1; (k < n); k++){
1092         int ln1, lstart;
1093         ln1 = k - j - 1;
1094         if(ln1+i-1>MAXLOOP) break;
1095
1096         lstart = ln1+i-1+n-MAXLOOP;
1097         if(lstart<k+TURN+1) lstart = k + TURN + 1;
1098         for(l=lstart; (l <= n); l++){
1099             int ln2, type2;
1100             ln2 = (i - 1) + (n - l);
1101             if((ln1+ln2) > MAXLOOP) continue;
1102
1103             type2 = ptype[my_iindx[k]-l];
1104             if(!type) continue;
1105             type2 = rtype[type2];
1106             qt += qb[my_iindx[i]-j] * qb[my_iindx[k]-l] * exp_E_IntLoop(ln2, ln1, type2, type, S1[l+1], S1[k-1], S1[i-1], S1[j+1], pf_params) * scale[ln1 + ln2];
1107             /* found an exterior interior loop? also this time, we can go straight  */
1108             /* forward and backtracking the both enclosed parts and we're done      */
1109             if(qt>r){ backtrack(i,j); backtrack(k,l); return pstruc;}
1110         }
1111       } /* end of kl double loop */
1112     }
1113   } /* end of ij double loop  */
1114   {
1115     /* as we reach this part, we have to search for our barrier between qm and qm2  */
1116     qt = 0.;
1117     r = urn()*qmo;
1118     for(k=TURN+2; k<n-2*TURN-3; k++){
1119       qt += qm[my_iindx[1]-k] * qm2[k+1] * expMLclosing;
1120       /* backtrack in qm and qm2 if we've found a valid barrier k  */
1121       if(qt>r){ backtrack_qm(1,k); backtrack_qm2(k+1,n); return pstruc;}
1122     }
1123   }
1124   /* if we reach the actual end of this function, an error has occured  */
1125   /* cause we HAVE TO find an exterior loop or an open chain!!!         */
1126   nrerror("backtracking failed in exterior loop");
1127   return pstruc;
1128 }
1129
1130 PRIVATE void backtrack_qm(int i, int j){
1131   /* divide multiloop into qm and qm1  */
1132   double qmt, r;
1133   int k;
1134   while(j>i){
1135     /* now backtrack  [i ... j] in qm[] */
1136     r = urn() * qm[my_iindx[i] - j];
1137     qmt = qm1[jindx[j]+i]; k=i;
1138     if(qmt<r)
1139       for(k=i+1; k<=j; k++){
1140         qmt += (qm[my_iindx[i]-(k-1)]+expMLbase[k-i])*qm1[jindx[j]+k];
1141         if(qmt >= r) break;
1142       }
1143     if(k>j) nrerror("backtrack failed in qm");
1144
1145     backtrack_qm1(k,j);
1146
1147     if(k<i+TURN) break; /* no more pairs */
1148     r = urn() * (qm[my_iindx[i]-(k-1)] + expMLbase[k-i]);
1149     if(expMLbase[k-i] >= r) break; /* no more pairs */
1150     j = k-1;
1151   }
1152 }
1153
1154 PRIVATE void backtrack_qm1(int i,int j){
1155   /* i is paired to l, i<l<j; backtrack in qm1 to find l */
1156   int ii, l, type;
1157   double qt, r;
1158   r = urn() * qm1[jindx[j]+i];
1159   ii = my_iindx[i];
1160   for (qt=0., l=i+TURN+1; l<=j; l++) {
1161     type = ptype[ii-l];
1162     if (type)
1163       qt +=  qb[ii-l] * exp_E_MLstem(type, S1[i-1], S1[l+1], pf_params) * expMLbase[j-l];
1164     if (qt>=r) break;
1165   }
1166   if (l>j) nrerror("backtrack failed in qm1");
1167   backtrack(i,l);
1168 }
1169
1170 PRIVATE void backtrack_qm2(int k, int n){
1171   double qom2t, r;
1172   int u;
1173   r= urn()*qm2[k];
1174   /* we have to search for our barrier u between qm1 and qm1  */
1175   for (qom2t = 0.,u=k+TURN+1; u<n-TURN-1; u++){
1176     qom2t += qm1[jindx[u]+k]*qm1[jindx[n]+(u+1)];
1177     if(qom2t > r) break;
1178   }
1179   if(u==n-TURN) nrerror("backtrack failed in qm2");
1180   backtrack_qm1(k,u);
1181   backtrack_qm1(u+1,n);
1182 }
1183
1184 PRIVATE void backtrack(int i, int j){
1185   int noGUclosure = pf_params->model_details.noGUclosure;
1186
1187   do {
1188     double r, qbt1;
1189     int k, l, type, u, u1;
1190
1191     pstruc[i-1] = '('; pstruc[j-1] = ')';
1192
1193     r = urn() * qb[my_iindx[i]-j];
1194     type = ptype[my_iindx[i]-j];
1195     u = j-i-1;
1196     /*hairpin contribution*/
1197     if (((type==3)||(type==4))&&noGUclosure) qbt1 = 0;
1198     else
1199       qbt1 = exp_E_Hairpin(u, type, S1[i+1], S1[j-1], sequence+i-1, pf_params)*
1200         scale[u+2]; /* add scale[u+2] */
1201
1202     if (qbt1>=r) return; /* found the hairpin we're done */
1203
1204     for (k=i+1; k<=MIN2(i+MAXLOOP+1,j-TURN-2); k++) {
1205       u1 = k-i-1;
1206       for (l=MAX2(k+TURN+1,j-1-MAXLOOP+u1); l<j; l++) {
1207         int type_2;
1208         type_2 = ptype[my_iindx[k]-l];
1209         if (type_2) {
1210           type_2 = rtype[type_2];
1211           /* add *scale[u1+u2+2] */
1212           qbt1 += qb[my_iindx[k]-l] * (scale[u1+j-l+1] *
1213             exp_E_IntLoop(u1, j-l-1, type, type_2,
1214                           S1[i+1], S1[j-1], S1[k-1], S1[l+1], pf_params));
1215         }
1216         if (qbt1 > r) break;
1217       }
1218       if (qbt1 > r) break;
1219     }
1220     if (l<j) {
1221       i=k; j=l;
1222     }
1223     else break;
1224   } while (1);
1225
1226   /* backtrack in multi-loop */
1227   {
1228     double r, qt;
1229     int k, ii, jj;
1230
1231     i++; j--;
1232     /* find the first split index */
1233     ii = my_iindx[i]; /* ii-j=[i,j] */
1234     jj = jindx[j]; /* jj+i=[j,i] */
1235     for (qt=0., k=i+1; k<j; k++) qt += qm[ii-(k-1)]*qm1[jj+k];
1236     r = urn() * qt;
1237     for (qt=0., k=i+1; k<j; k++) {
1238       qt += qm[ii-(k-1)]*qm1[jj+k];
1239       if (qt>=r) break;
1240     }
1241     if (k>=j) nrerror("backtrack failed, can't find split index ");
1242
1243     backtrack_qm1(k, j);
1244
1245     j = k-1;
1246     backtrack_qm(i,j);
1247   }
1248 }
1249
1250 PUBLIC void assign_plist_from_pr(plist **pl, FLT_OR_DBL *probs, int length, double cut_off){
1251   int i, j, n, count, *index;
1252   count = 0;
1253   n     = 2;
1254
1255   index = get_iindx(length);
1256
1257   /* first guess of the size needed for pl */
1258   *pl = (plist *)space(n*length*sizeof(plist));
1259
1260   for (i=1; i<length; i++) {
1261     for (j=i+1; j<=length; j++) {
1262       /* skip all entries below the cutoff */
1263       if (probs[index[i]-j] < cut_off) continue;
1264       /* do we need to allocate more memory? */
1265       if (count == n * length - 1){
1266         n *= 2;
1267         *pl = (plist *)xrealloc(*pl, n * length * sizeof(plist));
1268       }
1269       (*pl)[count].i    = i;
1270       (*pl)[count].j    = j;
1271       (*pl)[count].p  = probs[index[i] - j];
1272       (*pl)[count++].type = 0;
1273     }
1274   }
1275   /* mark the end of pl */
1276   (*pl)[count].i   = 0;
1277   (*pl)[count].j   = 0;
1278   (*pl)[count].p = 0.;
1279   (*pl)[count++].type = 0;
1280   /* shrink memory to actual size needed */
1281   *pl = (plist *)xrealloc(*pl, count * sizeof(plist));
1282
1283   free(index);
1284 }
1285
1286 /* this doesn't work if free_pf_arrays() is called before */
1287 PUBLIC void assign_plist_gquad_from_pr( plist **pl,
1288                                         int length,
1289                                         double cut_off){
1290
1291   int i, j, k, n, count, *index;
1292   count = 0;
1293   n     = 2;
1294
1295   if(!probs){ *pl = NULL; return;}
1296
1297   index = get_iindx(length);
1298
1299   /* first guess of the size needed for pl */
1300   *pl = (plist *)space(n*length*sizeof(plist));
1301
1302   for (i=1; i<length; i++) {
1303     for (j=i+1; j<=length; j++) {
1304       /* skip all entries below the cutoff */
1305       if (probs[index[i]-j] < cut_off) continue;
1306
1307       /* do we need to allocate more memory? */
1308       if (count == n * length - 1){
1309         n *= 2;
1310         *pl = (plist *)xrealloc(*pl, n * length * sizeof(plist));
1311       }
1312
1313       /* check for presence of gquadruplex */
1314       if((S[i] == 3) && (S[j] == 3)){
1315           /* add probability of a gquadruplex at position (i,j)
1316              for dot_plot
1317           */
1318           (*pl)[count].i      = i;
1319           (*pl)[count].j      = j;
1320           (*pl)[count].p      = probs[index[i] - j];
1321           (*pl)[count++].type = 1;
1322           /* now add the probabilies of it's actual pairing patterns */
1323           plist *inner, *ptr;
1324           inner = get_plist_gquad_from_pr(S, i, j, G, probs, scale, pf_params);
1325           for(ptr=inner; ptr->i != 0; ptr++){
1326               if (count == n * length - 1){
1327                 n *= 2;
1328                 *pl = (plist *)xrealloc(*pl, n * length * sizeof(plist));
1329               }
1330               /* check if we've already seen this pair */
1331               for(k = 0; k < count; k++)
1332                 if(((*pl)[k].i == ptr->i) && ((*pl)[k].j == ptr->j))
1333                   break;
1334               (*pl)[k].i      = ptr->i;
1335               (*pl)[k].j      = ptr->j;
1336               (*pl)[k].type = 0;
1337               if(k == count){
1338                 (*pl)[k].p  = ptr->p;
1339                 count++;
1340               } else
1341                 (*pl)[k].p  += ptr->p;
1342           }
1343       } else {
1344           (*pl)[count].i      = i;
1345           (*pl)[count].j      = j;
1346           (*pl)[count].p      = probs[index[i] - j];
1347           (*pl)[count++].type = 0;
1348       }
1349     }
1350   }
1351   /* mark the end of pl */
1352   (*pl)[count].i   = 0;
1353   (*pl)[count].j   = 0;
1354   (*pl)[count++].p = 0.;
1355   /* shrink memory to actual size needed */
1356   *pl = (plist *)xrealloc(*pl, count * sizeof(plist));
1357   free(index);
1358 }
1359
1360 /* this doesn't work if free_pf_arrays() is called before */
1361 PUBLIC char *get_centroid_struct_gquad_pr( int length,
1362                                           double *dist){
1363
1364   /* compute the centroid structure of the ensemble, i.e. the strutcure
1365      with the minimal average distance to all other structures
1366      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
1367      Thus, the centroid is simply the structure containing all pairs with
1368      p_ij>0.5 */
1369   int i,j, k;
1370   double p;
1371   char  *centroid;
1372   int   *my_iindx = get_iindx(length);
1373
1374   if (probs == NULL)
1375     nrerror("get_centroid_struct_pr: probs==NULL!");
1376
1377   *dist = 0.;
1378   centroid = (char *) space((length+1)*sizeof(char));
1379   for (i=0; i<length; i++) centroid[i]='.';
1380   for (i=1; i<=length; i++)
1381     for (j=i+TURN+1; j<=length; j++) {
1382       if ((p=probs[my_iindx[i]-j])>0.5) {
1383         /* check for presence of gquadruplex */
1384         if((S[i] == 3) && (S[j] == 3)){
1385           int L, l[3];
1386           get_gquad_pattern_pf(S, i, j, pf_params, &L, l);
1387           for(k=0;k<L;k++){
1388             centroid[i+k-1]\
1389             = centroid[i+k+L+l[0]-1]\
1390             = centroid[i+k+2*L+l[0]+l[1]-1]\
1391             = centroid[i+k+3*L+l[0]+l[1]+l[2]-1]\
1392             = '+';
1393           }
1394           /* skip everything within the gquad */
1395           i = j; j = j+TURN+1;
1396           *dist += (1-p); /* right? */
1397           break;
1398         } else {
1399             centroid[i-1] = '(';
1400             centroid[j-1] = ')';
1401         }
1402         *dist += (1-p);
1403       } else
1404         *dist += p;
1405     }
1406   free(my_iindx);
1407   centroid[length] = '\0';
1408   return centroid;
1409 }
1410
1411 /* this function is a threadsafe replacement for centroid() */
1412 PUBLIC char *get_centroid_struct_pl(int length, double *dist, plist *pl){
1413   /* compute the centroid structure of the ensemble, i.e. the strutcure
1414      with the minimal average distance to all other structures
1415      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
1416      Thus, the centroid is simply the structure containing all pairs with
1417      p_ij>0.5 */
1418   int i;
1419   char *centroid;
1420
1421   if (pl==NULL)
1422     nrerror("get_centroid_struct: pl==NULL!");
1423
1424   *dist = 0.;
1425   centroid = (char *) space((length+1)*sizeof(char));
1426   for (i=0; i<length; i++) centroid[i]='.';
1427   for (i=0; pl[i].i>0; i++){
1428     if ((pl[i].p)>0.5) {
1429       centroid[pl[i].i-1] = '(';
1430       centroid[pl[i].j-1] = ')';
1431       *dist += (1-pl[i].p);
1432     } else
1433       *dist += pl[i].p;
1434   }
1435   centroid[length] = '\0';
1436   return centroid;
1437 }
1438
1439 /* this function is a threadsafe replacement for centroid() */
1440 PUBLIC char *get_centroid_struct_pr(int length, double *dist, FLT_OR_DBL *probs){
1441   /* compute the centroid structure of the ensemble, i.e. the strutcure
1442      with the minimal average distance to all other structures
1443      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
1444      Thus, the centroid is simply the structure containing all pairs with
1445      p_ij>0.5 */
1446   int i,j;
1447   double p;
1448   char  *centroid;
1449   int   *index = get_iindx(length);
1450
1451   if (probs == NULL)
1452     nrerror("get_centroid_struct_pr: probs==NULL!");
1453
1454   *dist = 0.;
1455   centroid = (char *) space((length+1)*sizeof(char));
1456   for (i=0; i<length; i++) centroid[i]='.';
1457   for (i=1; i<=length; i++)
1458     for (j=i+TURN+1; j<=length; j++) {
1459       if ((p=probs[index[i]-j])>0.5) {
1460         centroid[i-1] = '(';
1461         centroid[j-1] = ')';
1462         *dist += (1-p);
1463       } else
1464         *dist += p;
1465     }
1466   free(index);
1467   centroid[length] = '\0';
1468   return centroid;
1469 }
1470
1471 PUBLIC plist *stackProb(double cutoff){
1472   plist *pl;
1473   int i,j,plsize=256;
1474   int num = 0;
1475
1476   if (probs==NULL)
1477     nrerror("probs==NULL. You need to call pf_fold() before stackProb()");
1478
1479   int length  = S[0];
1480   int *index  = get_iindx(length);
1481
1482   pl = (plist *) space(plsize*sizeof(plist));
1483
1484   for (i=1; i<length; i++)
1485     for (j=i+TURN+3; j<=length; j++) {
1486       double p;
1487       if((p=probs[index[i]-j]) < cutoff) continue;
1488       if (qb[index[i+1]-(j-1)]<FLT_MIN) continue;
1489       p *= qb[index[i+1]-(j-1)]/qb[index[i]-j];
1490       p *= exp_E_IntLoop(0,0,ptype[index[i]-j],rtype[ptype[index[i+1]-(j-1)]],
1491                          0,0,0,0, pf_params)*scale[2];/* add *scale[u1+u2+2] */
1492       if (p>cutoff) {
1493         pl[num].i = i;
1494         pl[num].j = j;
1495         pl[num++].p = p;
1496         if (num>=plsize) {
1497           plsize *= 2;
1498           pl = xrealloc(pl, plsize*sizeof(plist));
1499         }
1500       }
1501     }
1502   pl[num].i=0;
1503   free(index);
1504   return pl;
1505 }
1506
1507 /*-------------------------------------------------------------------------*/
1508 /* make arrays used for pf_fold available to other routines */
1509 PUBLIC int get_pf_arrays( short **S_p,
1510                           short **S1_p,
1511                           char **ptype_p,
1512                           FLT_OR_DBL **qb_p,
1513                           FLT_OR_DBL **qm_p,
1514                           FLT_OR_DBL **q1k_p,
1515                           FLT_OR_DBL **qln_p){
1516
1517   if(qb == NULL) return(0); /* check if pf_fold() has been called */
1518   *S_p = S; *S1_p = S1; *ptype_p = ptype;
1519   *qb_p = qb; *qm_p = qm;
1520   *q1k_p = q1k; *qln_p = qln;
1521   return(1); /* success */
1522 }
1523
1524 /* get the free energy of a subsequence from the q[] array */
1525 PUBLIC double get_subseq_F(int i, int j){
1526   if (!q)
1527     nrerror("call pf_fold() to fill q[] array before calling get_subseq_F()");
1528   return ((-log(q[my_iindx[i]-j])-(j-i+1)*log(pf_params->pf_scale))*pf_params->kT/1000.0);
1529 }
1530
1531
1532 PUBLIC double mean_bp_distance(int length){
1533   return mean_bp_distance_pr(length, probs);
1534 }
1535
1536 PUBLIC double mean_bp_distance_pr(int length, FLT_OR_DBL *p){
1537   /* compute the mean base pair distance in the thermodynamic ensemble */
1538   /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
1539      this can be computed from the pair probs p_ij as
1540      <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
1541   int i,j;
1542   double d=0;
1543   int *index = get_iindx((unsigned int) length);
1544
1545   if (p==NULL)
1546     nrerror("p==NULL. You need to supply a valid probability matrix for mean_bp_distance_pr()");
1547
1548   for (i=1; i<=length; i++)
1549     for (j=i+TURN+1; j<=length; j++)
1550       d += p[index[i]-j] * (1-p[index[i]-j]);
1551
1552   free(index);
1553   return 2*d;
1554 }
1555
1556 PUBLIC FLT_OR_DBL *export_bppm(void){
1557   return probs;
1558 }
1559
1560 /*###########################################*/
1561 /*# deprecated functions below              #*/
1562 /*###########################################*/
1563
1564 /* this function is deprecated since it is not threadsafe */
1565 PUBLIC char *centroid(int length, double *dist) {
1566   /* compute the centroid structure of the ensemble, i.e. the strutcure
1567      with the minimal average distance to all other structures
1568      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
1569      Thus, the centroid is simply the structure containing all pairs with
1570      p_ij>0.5 */
1571   int i,j;
1572   double p;
1573   char *centroid;
1574
1575   if (pr==NULL)
1576     nrerror("pr==NULL. You need to call pf_fold() before centroid()");
1577
1578   *dist = 0.;
1579   centroid = (char *) space((length+1)*sizeof(char));
1580   for (i=0; i<length; i++) centroid[i]='.';
1581   for (i=1; i<=length; i++)
1582     for (j=i+TURN+1; j<=length; j++) {
1583       if ((p=pr[my_iindx[i]-j])>0.5) {
1584         centroid[i-1] = '(';
1585         centroid[j-1] = ')';
1586         *dist += (1-p);
1587       } else
1588         *dist += p;
1589     }
1590   return centroid;
1591 }
1592
1593
1594 /* This function is deprecated since it uses the global array pr for calculations */
1595 PUBLIC double mean_bp_dist(int length) {
1596   /* compute the mean base pair distance in the thermodynamic ensemble */
1597   /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
1598      this can be computed from the pair probs p_ij as
1599      <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
1600   int i,j;
1601   double d=0;
1602
1603   if (pr==NULL)
1604     nrerror("pr==NULL. You need to call pf_fold() before mean_bp_dist()");
1605
1606   for (i=1; i<=length; i++)
1607     for (j=i+TURN+1; j<=length; j++)
1608       d += pr[my_iindx[i]-j] * (1-pr[my_iindx[i]-j]);
1609   return 2*d;
1610 }
1611
1612 /*----------------------------------------------------------------------*/
1613 PUBLIC double expHairpinEnergy(int u, int type, short si1, short sj1,
1614                                 const char *string) {
1615 /* compute Boltzmann weight of a hairpin loop, multiply by scale[u+2] */
1616   double q, kT;
1617   kT = pf_params->kT;   /* kT in cal/mol  */
1618   if(u <= 30)
1619     q = pf_params->exphairpin[u];
1620   else
1621     q = pf_params->exphairpin[30] * exp( -(pf_params->lxc*log( u/30.))*10./kT);
1622   if ((tetra_loop)&&(u==4)) {
1623     char tl[7]={0}, *ts;
1624     strncpy(tl, string, 6);
1625     if ((ts=strstr(pf_params->Tetraloops, tl)))
1626       return (pf_params->exptetra[(ts-pf_params->Tetraloops)/7]);
1627   }
1628   if ((tetra_loop)&&(u==6)) {
1629     char tl[9]={0}, *ts;
1630     strncpy(tl, string, 6);
1631     if ((ts=strstr(pf_params->Hexaloops, tl)))
1632       return  (pf_params->exphex[(ts-pf_params->Hexaloops)/9]);
1633   }
1634   if (u==3) {
1635     char tl[6]={0}, *ts;
1636     strncpy(tl, string, 5);
1637     if ((ts=strstr(pf_params->Triloops, tl)))
1638       return (pf_params->exptri[(ts-pf_params->Triloops)/6]);
1639     if (type>2)
1640       q *= pf_params->expTermAU;
1641   }
1642   else /* no mismatches for tri-loops */
1643     q *= pf_params->expmismatchH[type][si1][sj1];
1644
1645   return q;
1646 }
1647 PUBLIC double expLoopEnergy(int u1, int u2, int type, int type2,
1648                              short si1, short sj1, short sp1, short sq1) {
1649 /* compute Boltzmann weight of interior loop,
1650    multiply by scale[u1+u2+2] for scaling */
1651   double z=0;
1652   int no_close = 0;
1653
1654   if ((no_closingGU) && ((type2==3)||(type2==4)||(type==2)||(type==4)))
1655     no_close = 1;
1656
1657   if ((u1==0) && (u2==0)) /* stack */
1658     z = pf_params->expstack[type][type2];
1659   else if (no_close==0) {
1660     if ((u1==0)||(u2==0)) { /* bulge */
1661       int u;
1662       u = (u1==0)?u2:u1;
1663       z = pf_params->expbulge[u];
1664       if (u2+u1==1) z *= pf_params->expstack[type][type2];
1665       else {
1666         if (type>2) z *= pf_params->expTermAU;
1667         if (type2>2) z *= pf_params->expTermAU;
1668       }
1669     }
1670     else {     /* interior loop */
1671       if (u1+u2==2) /* size 2 is special */
1672         z = pf_params->expint11[type][type2][si1][sj1];
1673       else if ((u1==1) && (u2==2))
1674         z = pf_params->expint21[type][type2][si1][sq1][sj1];
1675       else if ((u1==2) && (u2==1))
1676         z = pf_params->expint21[type2][type][sq1][si1][sp1];
1677       else if ((u1==2) && (u2==2))
1678         z = pf_params->expint22[type][type2][si1][sp1][sq1][sj1];
1679       else if (((u1==2)&&(u2==3))||((u1==3)&&(u2==2))){ /*2-3 is special*/
1680         z = pf_params->expinternal[5]*
1681           pf_params->expmismatch23I[type][si1][sj1]*
1682           pf_params->expmismatch23I[type2][sq1][sp1];
1683         z *= pf_params->expninio[2][1];
1684       }
1685       else if ((u1==1)||(u2==1)) {  /*1-n is special*/
1686         z = pf_params->expinternal[u1+u2]*
1687           pf_params->expmismatch1nI[type][si1][sj1]*
1688           pf_params->expmismatch1nI[type2][sq1][sp1];
1689         z *= pf_params->expninio[2][abs(u1-u2)];
1690       }
1691       else {
1692         z = pf_params->expinternal[u1+u2]*
1693           pf_params->expmismatchI[type][si1][sj1]*
1694           pf_params->expmismatchI[type2][sq1][sp1];
1695         z *= pf_params->expninio[2][abs(u1-u2)];
1696       }
1697     }
1698   }
1699   return z;
1700 }
1701
1702 PUBLIC void init_pf_circ_fold(int length){
1703 /* DO NOTHING */
1704 }
1705
1706 PUBLIC void init_pf_fold(int length){
1707 /* DO NOTHING */
1708 }
1709
1710