API tools faq
paste
Advertisement
Guest User

part_func.c

a guest
Nov 18th, 2013
318
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 54.41 KB | None | 0 0
raw download clone embed print report
  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.     }
  1273.   }
  1274.   /* mark the end of pl */
  1275.   (*pl)[count].i   = 0;
  1276.   (*pl)[count].j   = 0;
  1277.   (*pl)[count++].p = 0.;
  1278.   /* shrink memory to actual size needed */
  1279.   *pl = (plist *)xrealloc(*pl, count * sizeof(plist));
  1280.   free(index);
  1281. }
  1282.  
  1283. /* this doesn't work if free_pf_arrays() is called before */
  1284. PUBLIC void assign_plist_gquad_from_pr( plist **pl,
  1285.                                         int length,
  1286.                                         double cut_off){
  1287.  
  1288.   int i, j, k, n, count, *index;
  1289.   count = 0;
  1290.   n     = 2;
  1291.  
  1292.   if(!probs){ *pl = NULL; return;}
  1293.  
  1294.   index = get_iindx(length);
  1295.  
  1296.   /* first guess of the size needed for pl */
  1297.   *pl = (plist *)space(n*length*sizeof(plist));
  1298.  
  1299.   for (i=1; i<length; i++) {
  1300.     for (j=i+1; j<=length; j++) {
  1301.       /* skip all entries below the cutoff */
  1302.       if (probs[index[i]-j] < cut_off) continue;
  1303.  
  1304.       /* do we need to allocate more memory? */
  1305.       if (count == n * length - 1){
  1306.         n *= 2;
  1307.         *pl = (plist *)xrealloc(*pl, n * length * sizeof(plist));
  1308.       }
  1309.  
  1310.       /* check for presence of gquadruplex */
  1311.       if((S[i] == 3) && (S[j] == 3)){
  1312.           /* add probability of a gquadruplex at position (i,j)
  1313.              for dot_plot
  1314.           */
  1315.           (*pl)[count].i      = i;
  1316.           (*pl)[count].j      = j;
  1317.           (*pl)[count].p      = probs[index[i] - j];
  1318.           (*pl)[count++].type = 1;
  1319.           /* now add the probabilies of it's actual pairing patterns */
  1320.           plist *inner, *ptr;
  1321.           inner = get_plist_gquad_from_pr(S, i, j, G, probs, scale, pf_params);
  1322.           for(ptr=inner; ptr->i != 0; ptr++){
  1323.               if (count == n * length - 1){
  1324.                 n *= 2;
  1325.                 *pl = (plist *)xrealloc(*pl, n * length * sizeof(plist));
  1326.               }
  1327.               /* check if we've already seen this pair */
  1328.               for(k = 0; k < count; k++)
  1329.                 if(((*pl)[k].i == ptr->i) && ((*pl)[k].j == ptr->j))
  1330.                   break;
  1331.               (*pl)[k].i      = ptr->i;
  1332.               (*pl)[k].j      = ptr->j;
  1333.               (*pl)[k].type = 0;
  1334.               if(k == count){
  1335.                 (*pl)[k].p  = ptr->p;
  1336.                 count++;
  1337.               } else
  1338.                 (*pl)[k].p  += ptr->p;
  1339.           }
  1340.       } else {
  1341.           (*pl)[count].i      = i;
  1342.           (*pl)[count].j      = j;
  1343.           (*pl)[count].p      = probs[index[i] - j];
  1344.           (*pl)[count++].type = 0;
  1345.       }
  1346.     }
  1347.   }
  1348.   /* mark the end of pl */
  1349.   (*pl)[count].i   = 0;
  1350.   (*pl)[count].j   = 0;
  1351.   (*pl)[count++].p = 0.;
  1352.   /* shrink memory to actual size needed */
  1353.   *pl = (plist *)xrealloc(*pl, count * sizeof(plist));
  1354.   free(index);
  1355. }
  1356.  
  1357. /* this doesn't work if free_pf_arrays() is called before */
  1358. PUBLIC char *get_centroid_struct_gquad_pr( int length,
  1359.                                           double *dist){
  1360.  
  1361.   /* compute the centroid structure of the ensemble, i.e. the strutcure
  1362.      with the minimal average distance to all other structures
  1363.      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
  1364.      Thus, the centroid is simply the structure containing all pairs with
  1365.      p_ij>0.5 */
  1366.   int i,j, k;
  1367.   double p;
  1368.   char  *centroid;
  1369.   int   *my_iindx = get_iindx(length);
  1370.  
  1371.   if (probs == NULL)
  1372.     nrerror("get_centroid_struct_pr: probs==NULL!");
  1373.  
  1374.   *dist = 0.;
  1375.   centroid = (char *) space((length+1)*sizeof(char));
  1376.   for (i=0; i<length; i++) centroid[i]='.';
  1377.   for (i=1; i<=length; i++)
  1378.     for (j=i+TURN+1; j<=length; j++) {
  1379.       if ((p=probs[my_iindx[i]-j])>0.5) {
  1380.         /* check for presence of gquadruplex */
  1381.         if((S[i] == 3) && (S[j] == 3)){
  1382.           int L, l[3];
  1383.           get_gquad_pattern_pf(S, i, j, pf_params, &L, l);
  1384.           for(k=0;k<L;k++){
  1385.             centroid[i+k-1]\
  1386.             = centroid[i+k+L+l[0]-1]\
  1387.             = centroid[i+k+2*L+l[0]+l[1]-1]\
  1388.             = centroid[i+k+3*L+l[0]+l[1]+l[2]-1]\
  1389.             = '+';
  1390.           }
  1391.           /* skip everything within the gquad */
  1392.           i = j; j = j+TURN+1;
  1393.           *dist += (1-p); /* right? */
  1394.           break;
  1395.         } else {
  1396.             centroid[i-1] = '(';
  1397.             centroid[j-1] = ')';
  1398.         }
  1399.         *dist += (1-p);
  1400.       } else
  1401.         *dist += p;
  1402.     }
  1403.   free(my_iindx);
  1404.   centroid[length] = '\0';
  1405.   return centroid;
  1406. }
  1407.  
  1408. /* this function is a threadsafe replacement for centroid() */
  1409. PUBLIC char *get_centroid_struct_pl(int length, double *dist, plist *pl){
  1410.   /* compute the centroid structure of the ensemble, i.e. the strutcure
  1411.      with the minimal average distance to all other structures
  1412.      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
  1413.      Thus, the centroid is simply the structure containing all pairs with
  1414.      p_ij>0.5 */
  1415.   int i;
  1416.   char *centroid;
  1417.  
  1418.   if (pl==NULL)
  1419.     nrerror("get_centroid_struct: pl==NULL!");
  1420.  
  1421.   *dist = 0.;
  1422.   centroid = (char *) space((length+1)*sizeof(char));
  1423.   for (i=0; i<length; i++) centroid[i]='.';
  1424.   for (i=0; pl[i].i>0; i++){
  1425.     if ((pl[i].p)>0.5) {
  1426.       centroid[pl[i].i-1] = '(';
  1427.       centroid[pl[i].j-1] = ')';
  1428.       *dist += (1-pl[i].p);
  1429.     } else
  1430.       *dist += pl[i].p;
  1431.   }
  1432.   centroid[length] = '\0';
  1433.   return centroid;
  1434. }
  1435.  
  1436. /* this function is a threadsafe replacement for centroid() */
  1437. PUBLIC char *get_centroid_struct_pr(int length, double *dist, FLT_OR_DBL *probs){
  1438.   /* compute the centroid structure of the ensemble, i.e. the strutcure
  1439.      with the minimal average distance to all other structures
  1440.      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
  1441.      Thus, the centroid is simply the structure containing all pairs with
  1442.      p_ij>0.5 */
  1443.   int i,j;
  1444.   double p;
  1445.   char  *centroid;
  1446.   int   *index = get_iindx(length);
  1447.  
  1448.   if (probs == NULL)
  1449.     nrerror("get_centroid_struct_pr: probs==NULL!");
  1450.  
  1451.   *dist = 0.;
  1452.   centroid = (char *) space((length+1)*sizeof(char));
  1453.   for (i=0; i<length; i++) centroid[i]='.';
  1454.   for (i=1; i<=length; i++)
  1455.     for (j=i+TURN+1; j<=length; j++) {
  1456.       if ((p=probs[index[i]-j])>0.5) {
  1457.         centroid[i-1] = '(';
  1458.         centroid[j-1] = ')';
  1459.         *dist += (1-p);
  1460.       } else
  1461.         *dist += p;
  1462.     }
  1463.   free(index);
  1464.   centroid[length] = '\0';
  1465.   return centroid;
  1466. }
  1467.  
  1468. PUBLIC plist *stackProb(double cutoff){
  1469.   plist *pl;
  1470.   int i,j,plsize=256;
  1471.   int num = 0;
  1472.  
  1473.   if (probs==NULL)
  1474.     nrerror("probs==NULL. You need to call pf_fold() before stackProb()");
  1475.  
  1476.   int length  = S[0];
  1477.   int *index  = get_iindx(length);
  1478.  
  1479.   pl = (plist *) space(plsize*sizeof(plist));
  1480.  
  1481.   for (i=1; i<length; i++)
  1482.     for (j=i+TURN+3; j<=length; j++) {
  1483.       double p;
  1484.       if((p=probs[index[i]-j]) < cutoff) continue;
  1485.       if (qb[index[i+1]-(j-1)]<FLT_MIN) continue;
  1486.       p *= qb[index[i+1]-(j-1)]/qb[index[i]-j];
  1487.       p *= exp_E_IntLoop(0,0,ptype[index[i]-j],rtype[ptype[index[i+1]-(j-1)]],
  1488.                          0,0,0,0, pf_params)*scale[2];/* add *scale[u1+u2+2] */
  1489.       if (p>cutoff) {
  1490.         pl[num].i = i;
  1491.         pl[num].j = j;
  1492.         pl[num++].p = p;
  1493.         if (num>=plsize) {
  1494.           plsize *= 2;
  1495.           pl = xrealloc(pl, plsize*sizeof(plist));
  1496.         }
  1497.       }
  1498.     }
  1499.   pl[num].i=0;
  1500.   free(index);
  1501.   return pl;
  1502. }
  1503.  
  1504. /*-------------------------------------------------------------------------*/
  1505. /* make arrays used for pf_fold available to other routines */
  1506. PUBLIC int get_pf_arrays( short **S_p,
  1507.                           short **S1_p,
  1508.                           char **ptype_p,
  1509.                           FLT_OR_DBL **qb_p,
  1510.                           FLT_OR_DBL **qm_p,
  1511.                           FLT_OR_DBL **q1k_p,
  1512.                           FLT_OR_DBL **qln_p){
  1513.  
  1514.   if(qb == NULL) return(0); /* check if pf_fold() has been called */
  1515.   *S_p = S; *S1_p = S1; *ptype_p = ptype;
  1516.   *qb_p = qb; *qm_p = qm;
  1517.   *q1k_p = q1k; *qln_p = qln;
  1518.   return(1); /* success */
  1519. }
  1520.  
  1521. /* get the free energy of a subsequence from the q[] array */
  1522. PUBLIC double get_subseq_F(int i, int j){
  1523.   if (!q)
  1524.     nrerror("call pf_fold() to fill q[] array before calling get_subseq_F()");
  1525.   return ((-log(q[my_iindx[i]-j])-(j-i+1)*log(pf_params->pf_scale))*pf_params->kT/1000.0);
  1526. }
  1527.  
  1528.  
  1529. PUBLIC double mean_bp_distance(int length){
  1530.   return mean_bp_distance_pr(length, probs);
  1531. }
  1532.  
  1533. PUBLIC double mean_bp_distance_pr(int length, FLT_OR_DBL *p){
  1534.   /* compute the mean base pair distance in the thermodynamic ensemble */
  1535.   /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
  1536.      this can be computed from the pair probs p_ij as
  1537.      <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
  1538.   int i,j;
  1539.   double d=0;
  1540.   int *index = get_iindx((unsigned int) length);
  1541.  
  1542.   if (p==NULL)
  1543.     nrerror("p==NULL. You need to supply a valid probability matrix for mean_bp_distance_pr()");
  1544.  
  1545.   for (i=1; i<=length; i++)
  1546.     for (j=i+TURN+1; j<=length; j++)
  1547.       d += p[index[i]-j] * (1-p[index[i]-j]);
  1548.  
  1549.   free(index);
  1550.   return 2*d;
  1551. }
  1552.  
  1553. PUBLIC FLT_OR_DBL *export_bppm(void){
  1554.   return probs;
  1555. }
  1556.  
  1557. /*###########################################*/
  1558. /*# deprecated functions below              #*/
  1559. /*###########################################*/
  1560.  
  1561. /* this function is deprecated since it is not threadsafe */
  1562. PUBLIC char *centroid(int length, double *dist) {
  1563.   /* compute the centroid structure of the ensemble, i.e. the strutcure
  1564.      with the minimal average distance to all other structures
  1565.      <d(S)> = \sum_{(i,j) \in S} (1-p_{ij}) + \sum_{(i,j) \notin S} p_{ij}
  1566.      Thus, the centroid is simply the structure containing all pairs with
  1567.      p_ij>0.5 */
  1568.   int i,j;
  1569.   double p;
  1570.   char *centroid;
  1571.  
  1572.   if (pr==NULL)
  1573.     nrerror("pr==NULL. You need to call pf_fold() before centroid()");
  1574.  
  1575.   *dist = 0.;
  1576.   centroid = (char *) space((length+1)*sizeof(char));
  1577.   for (i=0; i<length; i++) centroid[i]='.';
  1578.   for (i=1; i<=length; i++)
  1579.     for (j=i+TURN+1; j<=length; j++) {
  1580.       if ((p=pr[my_iindx[i]-j])>0.5) {
  1581.         centroid[i-1] = '(';
  1582.         centroid[j-1] = ')';
  1583.         *dist += (1-p);
  1584.       } else
  1585.         *dist += p;
  1586.     }
  1587.   return centroid;
  1588. }
  1589.  
  1590.  
  1591. /* This function is deprecated since it uses the global array pr for calculations */
  1592. PUBLIC double mean_bp_dist(int length) {
  1593.   /* compute the mean base pair distance in the thermodynamic ensemble */
  1594.   /* <d> = \sum_{a,b} p_a p_b d(S_a,S_b)
  1595.      this can be computed from the pair probs p_ij as
  1596.      <d> = \sum_{ij} p_{ij}(1-p_{ij}) */
  1597.   int i,j;
  1598.   double d=0;
  1599.  
  1600.   if (pr==NULL)
  1601.     nrerror("pr==NULL. You need to call pf_fold() before mean_bp_dist()");
  1602.  
  1603.   for (i=1; i<=length; i++)
  1604.     for (j=i+TURN+1; j<=length; j++)
  1605.       d += pr[my_iindx[i]-j] * (1-pr[my_iindx[i]-j]);
  1606.   return 2*d;
  1607. }
  1608.  
  1609. /*----------------------------------------------------------------------*/
  1610. PUBLIC double expHairpinEnergy(int u, int type, short si1, short sj1,
  1611.                                 const char *string) {
  1612. /* compute Boltzmann weight of a hairpin loop, multiply by scale[u+2] */
  1613.   double q, kT;
  1614.   kT = pf_params->kT;   /* kT in cal/mol  */
  1615.   if(u <= 30)
  1616.     q = pf_params->exphairpin[u];
  1617.   else
  1618.     q = pf_params->exphairpin[30] * exp( -(pf_params->lxc*log( u/30.))*10./kT);
  1619.   if ((tetra_loop)&&(u==4)) {
  1620.     char tl[7]={0}, *ts;
  1621.     strncpy(tl, string, 6);
  1622.     if ((ts=strstr(pf_params->Tetraloops, tl)))
  1623.       return (pf_params->exptetra[(ts-pf_params->Tetraloops)/7]);
  1624.   }
  1625.   if ((tetra_loop)&&(u==6)) {
  1626.     char tl[9]={0}, *ts;
  1627.     strncpy(tl, string, 6);
  1628.     if ((ts=strstr(pf_params->Hexaloops, tl)))
  1629.       return  (pf_params->exphex[(ts-pf_params->Hexaloops)/9]);
  1630.   }
  1631.   if (u==3) {
  1632.     char tl[6]={0}, *ts;
  1633.     strncpy(tl, string, 5);
  1634.     if ((ts=strstr(pf_params->Triloops, tl)))
  1635.       return (pf_params->exptri[(ts-pf_params->Triloops)/6]);
  1636.     if (type>2)
  1637.       q *= pf_params->expTermAU;
  1638.   }
  1639.   else /* no mismatches for tri-loops */
  1640.     q *= pf_params->expmismatchH[type][si1][sj1];
  1641.  
  1642.   return q;
  1643. }
  1644. PUBLIC double expLoopEnergy(int u1, int u2, int type, int type2,
  1645.                              short si1, short sj1, short sp1, short sq1) {
  1646. /* compute Boltzmann weight of interior loop,
  1647.    multiply by scale[u1+u2+2] for scaling */
  1648.   double z=0;
  1649.   int no_close = 0;
  1650.  
  1651.   if ((no_closingGU) && ((type2==3)||(type2==4)||(type==2)||(type==4)))
  1652.     no_close = 1;
  1653.  
  1654.   if ((u1==0) && (u2==0)) /* stack */
  1655.     z = pf_params->expstack[type][type2];
  1656.   else if (no_close==0) {
  1657.     if ((u1==0)||(u2==0)) { /* bulge */
  1658.       int u;
  1659.       u = (u1==0)?u2:u1;
  1660.       z = pf_params->expbulge[u];
  1661.       if (u2+u1==1) z *= pf_params->expstack[type][type2];
  1662.       else {
  1663.         if (type>2) z *= pf_params->expTermAU;
  1664.         if (type2>2) z *= pf_params->expTermAU;
  1665.       }
  1666.     }
  1667.     else {     /* interior loop */
  1668.       if (u1+u2==2) /* size 2 is special */
  1669.         z = pf_params->expint11[type][type2][si1][sj1];
  1670.       else if ((u1==1) && (u2==2))
  1671.         z = pf_params->expint21[type][type2][si1][sq1][sj1];
  1672.       else if ((u1==2) && (u2==1))
  1673.         z = pf_params->expint21[type2][type][sq1][si1][sp1];
  1674.       else if ((u1==2) && (u2==2))
  1675.         z = pf_params->expint22[type][type2][si1][sp1][sq1][sj1];
  1676.       else if (((u1==2)&&(u2==3))||((u1==3)&&(u2==2))){ /*2-3 is special*/
  1677.         z = pf_params->expinternal[5]*
  1678.           pf_params->expmismatch23I[type][si1][sj1]*
  1679.           pf_params->expmismatch23I[type2][sq1][sp1];
  1680.         z *= pf_params->expninio[2][1];
  1681.       }
  1682.       else if ((u1==1)||(u2==1)) {  /*1-n is special*/
  1683.         z = pf_params->expinternal[u1+u2]*
  1684.           pf_params->expmismatch1nI[type][si1][sj1]*
  1685.           pf_params->expmismatch1nI[type2][sq1][sp1];
  1686.         z *= pf_params->expninio[2][abs(u1-u2)];
  1687.       }
  1688.       else {
  1689.         z = pf_params->expinternal[u1+u2]*
  1690.           pf_params->expmismatchI[type][si1][sj1]*
  1691.           pf_params->expmismatchI[type2][sq1][sp1];
  1692.         z *= pf_params->expninio[2][abs(u1-u2)];
  1693.       }
  1694.     }
  1695.   }
  1696.   return z;
  1697. }
  1698.  
  1699. PUBLIC void init_pf_circ_fold(int length){
  1700. /* DO NOTHING */
  1701. }
  1702.  
  1703. PUBLIC void init_pf_fold(int length){
  1704. /* DO NOTHING */
  1705. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!