binaries/src/ViennaRNA/lib/LPfold.c

   1 /* Last changed Time-stamp: <2009-02-18 14:19:51 ivo> */
   2 /*
   3   local pair probabilities for RNA secondary structures
   4
   5   Stephan Bernhart, Ivo L Hofacker
   6   Vienna RNA package
   7 */
   8 /*
   9   todo: compute energy z-score for each window
  10
  11 */
  12 #include <config.h>
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include <math.h>
  17 #include <float.h>    /* #defines FLT_MAX ... */
  18 #include "utils.h"
  19 #include "energy_par.h"
  20 #include "fold_vars.h"
  21 #include "pair_mat.h"
  22 #include "PS_dot.h"
  23 #include "part_func.h"
  24 #include "params.h"
  25 #include "loop_energies.h"
  26 #include "LPfold.h"
  27 #include "Lfold.h"
  28
  29 #ifdef _OPENMP
  30 #include <omp.h>
  31 #endif
  32
  33
  34 /*@unused@*/
  35 PRIVATE char rcsid[] UNUSED = "$Id: LPfold.c,v 1.8 2009/02/18 20:34:38 ivo Exp $";
  36
  37 #define ISOLATED  256.0
  38
  39 /*
  40 #################################
  41 # GLOBAL VARIABLES              #
  42 #################################
  43 */
  44
  45 /*
  46 #################################
  47 # PRIVATE VARIABLES             #
  48 #################################
  49 */
  50
  51 PRIVATE float       cutoff;
  52 PRIVATE int         num_p=0; /* for counting basepairs in pairlist pl, can actually be moved into pfl_fold */
  53 PRIVATE FLT_OR_DBL  *expMLbase=NULL;
  54 PRIVATE FLT_OR_DBL  **q=NULL, **qb=NULL, **qm=NULL, *qqm=NULL, *qqm1=NULL, *qq=NULL, *qq1=NULL, **pR=NULL, **qm2=NULL, **QI5=NULL,  **q2l=NULL, **qmb=NULL;/*,**QI3,*/
  55 PRIVATE FLT_OR_DBL  *prml=NULL, *prm_l=NULL, *prm_l1=NULL, *q1k=NULL, *qln=NULL;
  56 PRIVATE FLT_OR_DBL  *scale=NULL;
  57 PRIVATE char        **ptype=NULL; /* precomputed array of pair types */
  58 PRIVATE int         *jindx=NULL;
  59 PRIVATE int         *my_iindx=NULL;
  60 PRIVATE int         init_length = 0;  /* length in last call to init_pf_fold() */
  61 PRIVATE pf_paramT   *pf_params=NULL;
  62 PRIVATE short       *S=NULL, *S1=NULL;
  63 PRIVATE int         unpaired;
  64 PRIVATE int         ulength;
  65 PRIVATE int         pUoutput;
  66 PRIVATE double      alpha = 1.0;
  67
  68 #ifdef _OPENMP
  69
  70 /* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
  71 */
  72 #pragma omp threadprivate(cutoff, num_p, scale, ptype, jindx, my_iindx, init_length, pf_params,\
  73                           expMLbase, q, qb, qm, qqm, qqm1, qq, qq1, pR, qm2, QI5, q2l, qmb,\
  74                           prml, prm_l, prm_l1, q1k, qln,\
  75                           S, S1, unpaired, ulength, pUoutput, alpha)
  76
  77 #endif
  78
  79 /*
  80 #################################
  81 # PRIVATE FUNCTION DECLARATIONS #
  82 #################################
  83 */
  84
  85 PRIVATE void  init_partfunc_L(int length, pf_paramT *parameters);
  86 PRIVATE void  get_arrays_L(unsigned int length);
  87 PRIVATE void  free_pf_arrays_L(void);
  88 PRIVATE void  scale_pf_params(unsigned int length, pf_paramT *parameters);
  89 PRIVATE void  GetPtype(int j, int pairsize, const short *S, int n);
  90 PRIVATE void  FreeOldArrays(int i);
  91 PRIVATE void  GetNewArrays(int j, int winSize);
  92 PRIVATE void  printpbar(FLT_OR_DBL **prb,int winSize, int i, int n);
  93 PRIVATE plist *get_deppp(struct plist *pl, int start, int pairsize, int length);
  94 PRIVATE plist *get_plistW(struct plist *pl, int length, int start, FLT_OR_DBL **Tpr, int winSize);
  95 PRIVATE void  print_plist(int length, int start, FLT_OR_DBL **Tpr, int winSize, FILE *fp);
  96 PRIVATE void  compute_pU(int k, int ulength, double **pU, int winSize, int n, char *sequence);
  97 PRIVATE void  putoutpU(double **pU,int k, int ulength, FILE *fp);
  98 /*PRIVATE void make_ptypes(const short *S, const char *structure);*/
  99
 100 PRIVATE void putoutpU_splitup(double **pUx, int k, int ulength, FILE *fp, char ident);
 101 PRIVATE void compute_pU_splitup(int k, int ulength, double **pU,  double **pUO, double **pUH, double **pUI, double **pUM, int winSize,int n, char *sequence);
 102
 103 /*
 104 #################################
 105 # BEGIN OF FUNCTION DEFINITIONS #
 106 #################################
 107 */
 108
 109 PRIVATE void init_partfunc_L(int length, pf_paramT *parameters){
 110   if (length<1) nrerror("init_partfunc_L: length must be greater 0");
 111 #ifdef _OPENMP
 112 /* Explicitly turn off dynamic threads */
 113   omp_set_dynamic(0);
 114   free_pf_arrays_L(); /* free previous allocation */
 115 #else
 116   if (init_length>0) free_pf_arrays_L(); /* free previous allocation */
 117 #endif
 118
 119 #ifdef SUN4
 120   nonstandard_arithmetic();
 121 #else
 122 #ifdef HP9
 123   fpsetfastmode(1);
 124 #endif
 125 #endif
 126   make_pair_matrix();
 127   get_arrays_L((unsigned) length);
 128   scale_pf_params((unsigned) length, parameters);
 129
 130 #ifndef _OPENMP
 131   init_length = length;
 132 #endif
 133 }
 134
 135 PRIVATE void get_arrays_L(unsigned int length){
 136   /*arrays in 2 dimensions*/
 137
 138   q         = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *)*(length+1));
 139   qb        = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *)*(length+1));
 140   qm        = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *)*(length+1));
 141   pR        = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *)*(length+1));
 142   q1k       = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+1));
 143   qln       = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 144   qq        = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 145   qq1       = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 146   qqm       = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 147   qqm1      = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 148   prm_l     = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 149   prm_l1    = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 150   prml      = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+2));
 151   expMLbase = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+1));
 152   scale     = (FLT_OR_DBL *)  space(sizeof(FLT_OR_DBL)  *(length+1));
 153   ptype     = (char **)       space(sizeof(char *)      *(length+2));
 154
 155   if (ulength>0) {
 156     /* QI3 = (FLT_OR_DBL **) space((length+1)*sizeof(FLT_OR_DBL *));*/
 157     QI5 = (FLT_OR_DBL **) space((length+1)*sizeof(FLT_OR_DBL *));
 158     qmb = (FLT_OR_DBL **) space((length+1)*sizeof(FLT_OR_DBL *));
 159     qm2 = (FLT_OR_DBL **) space((length+1)*sizeof(FLT_OR_DBL *));
 160     q2l = (FLT_OR_DBL **) space((length+1)*sizeof(FLT_OR_DBL *));
 161   }
 162   my_iindx  = get_iindx(length);
 163   iindx     = get_iindx(length); /* for backward compatibility and Perl wrapper */
 164   jindx     = get_indx(length);
 165 }
 166
 167 PRIVATE void free_pf_arrays_L(void){
 168   if(q)         free(q);
 169   if(qb)        free(qb);
 170   if(qm)        free(qm);
 171   if(pR)        free(pR);
 172   if(qm2)       free(qm2);
 173   if(qq)        free(qq);
 174   if(qq1)       free(qq1);
 175   if(qqm)       free(qqm);
 176   if(qqm1)      free(qqm1);
 177   if(q1k)       free(q1k);
 178   if(qln)       free(qln);
 179   if(prm_l)     free(prm_l);
 180   if(prm_l1)    free(prm_l1);
 181   if(prml)      free(prml);
 182   if(expMLbase) free(expMLbase);
 183   if(scale)     free(scale);
 184   if(my_iindx)  free(my_iindx);
 185   if(iindx)     free(iindx); /* for backward compatibility and Perl wrapper */
 186   if(jindx)     free(jindx);
 187   if(ptype)     free(ptype);
 188   if(QI5)       free(QI5);
 189   if(qmb)       free(qmb);
 190   if(q2l)       free(q2l);
 191   if(pf_params) free(pf_params);
 192
 193   q = qb = qm = pR = QI5 = qmb = qm2 = q2l = NULL;
 194   qq = qq1 = qqm = qqm1 = q1k = qln = prml = prm_l = prm_l1 = expMLbase = NULL;
 195   my_iindx = jindx = iindx = NULL;
 196   pf_params = NULL;
 197   ptype     = NULL;
 198   scale = NULL;
 199
 200 #ifdef SUN4
 201   standard_arithmetic();
 202 #else
 203 #ifdef HP9
 204   fpsetfastmode(0);
 205 #endif
 206 #endif
 207
 208 #ifndef _OPENMP
 209   init_length=0;
 210 #endif
 211 }
 212
 213 PUBLIC void update_pf_paramsLP(int length){
 214   update_pf_paramsLP_par(length, NULL);
 215 }
 216
 217 PUBLIC void update_pf_paramsLP_par(int length, pf_paramT *parameters){
 218 #ifdef _OPENMP
 219   init_partfunc_L(length, parameters);
 220 #else
 221   if(parameters) init_partfunc_L(length, parameters);
 222   else if (length > init_length) init_partfunc_L(length, parameters);
 223   else {
 224     /*   make_pair_matrix();*/
 225     scale_pf_params((unsigned) length, parameters);
 226   }
 227 #endif
 228 }
 229
 230 PUBLIC plist *pfl_fold( char *sequence,
 231                         int winSize,
 232                         int pairSize,
 233                         float cutoffb,
 234                         double **pU,
 235                         struct plist **dpp2,
 236                         FILE *pUfp,
 237                         FILE *spup){
 238   return pfl_fold_par(sequence, winSize, pairSize, cutoffb, pU, dpp2, pUfp, spup, NULL);
 239 }
 240
 241 PUBLIC plist *pfl_fold_par( char *sequence,
 242                             int winSize,
 243                             int pairSize,
 244                             float cutoffb,
 245                             double **pU,
 246                             struct plist **dpp2,
 247                             FILE *pUfp,
 248                             FILE *spup,
 249                             pf_paramT *parameters){
 250
 251   int         n, m, i, j, k, l, u, u1, ii, type, type_2, tt, ov, do_dpp, simply_putout, noGUclosure;
 252   double      max_real;
 253   FLT_OR_DBL  temp, Qmax, prm_MLb, prmt, prmt1, qbt1, *tmp, expMLclosing;
 254   plist       *dpp, *pl;
 255   int split=0;
 256
 257   ov            = 0;
 258   Qmax          = 0;
 259   do_dpp        = 0;
 260   simply_putout = 0;
 261   dpp           = NULL;
 262   pl            = NULL;
 263   pUoutput      = 0;
 264   ulength       = 0;
 265   cutoff        = cutoffb;
 266
 267   if(pU != NULL)  ulength       = (int)pU[0][0]+0.49;
 268   if(spup !=NULL) simply_putout = 1; /*can't have one without the other*/
 269   if(pUfp!=NULL)  pUoutput      = 1;
 270   else if((pUoutput)&&(ulength!=0)){
 271     fprintf(stderr, "There was a problem with non existing File Pointer for unpaireds, terminating process\n");
 272     return pl;
 273   }
 274   dpp = *dpp2;
 275   if(dpp !=NULL)  do_dpp=1;
 276
 277   n = (int) strlen(sequence);
 278   if (n<TURN+2) return 0;
 279
 280 #ifdef _OPENMP
 281   /* always init everything since all global static variables are uninitialized when entering a thread */
 282   init_partfunc_L(n, parameters);
 283 #else
 284   if(parameters) init_partfunc_L(n, parameters);
 285   else if (n > init_length) init_partfunc_L(n, parameters);
 286   else if (fabs(pf_params->temperature - temperature)>1e-6) update_pf_paramsLP_par(n, parameters);
 287 #endif
 288
 289   expMLclosing  = pf_params->expMLclosing;
 290   noGUclosure   = pf_params->model_details.noGUclosure;
 291
 292
 293   max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
 294
 295   S   = encode_sequence(sequence, 0);
 296   S1  = encode_sequence(sequence, 1);
 297
 298   /*  make_ptypes(S, structure); das machmadochlieber lokal, ey!*/
 299
 300   /*here, I allocate memory for pU, if has to be saved, I allocate all in one go,
 301     if pU is put out and freed, I only allocate what I really need*/
 302
 303   if (ulength>0){
 304     if (pUoutput) {
 305       for (i=1; i<=ulength; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 306     }
 307     else {
 308       for (i=1; i<=n; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 309      }
 310   }
 311
 312   /*array initialization ; qb,qm,q
 313     qb,qm,q (i,j) are stored as ((n+1-i)*(n-i) div 2 + n+1-j */
 314   num_p = 0;
 315   pl    = (struct plist *)space(1000*sizeof(struct plist));
 316
 317
 318   /*ALWAYS q[i][j] => i>j!!*/
 319   for (j=1; j<MIN2(TURN+2,n); j++) { /*allocate start*/
 320     GetNewArrays(j, winSize);
 321     GetPtype(j,pairSize,S,n);
 322     for (i=1; i<=j; i++) q[i][j]=scale[(j-i+1)];
 323   }
 324   for (j=TURN+2;j<=n+winSize; j++) {
 325     if (j<=n) {
 326       GetNewArrays(j, winSize);
 327       GetPtype(j,pairSize,S,n);
 328       for (i=MAX2(1,j-winSize); i<=j/*-TURN*/; i++)
 329         q[i][j]=scale[(j-i+1)];
 330       for (i=j-TURN-1;i>=MAX2(1,(j-winSize+1)); i--) {
 331         /* construction of partition function of segment i,j*/
 332         /*firstly that given i bound to j : qb(i,j) */
 333         u = j-i-1;
 334         type = ptype[i][j];
 335         if (type!=0) {
 336           /*hairpin contribution*/
 337           if (((type==3)||(type==4))&&noGUclosure) qbt1 = 0;
 338           else
 339             qbt1 = exp_E_Hairpin(u, type, S1[i+1], S1[j-1], sequence+i-1, pf_params) * scale[u+2];
 340
 341           /* interior loops with interior pair k,l */
 342           for (k=i+1; k<=MIN2(i+MAXLOOP+1,j-TURN-2); k++) {
 343             u1 = k-i-1;
 344             for (l=MAX2(k+TURN+1,j-1-MAXLOOP+u1); l<j; l++) {
 345               type_2 = ptype[k][l];
 346               if (type_2) {
 347                 type_2 = rtype[type_2];
 348                 qbt1 += qb[k][l] *
 349                   exp_E_IntLoop(u1, j-l-1, type, type_2,
 350                                 S1[i+1], S1[j-1], S1[k-1], S1[l+1], pf_params) * scale[k-i+j-l];
 351               }
 352             }
 353           }
 354           /*multiple stem loop contribution*/
 355           ii = my_iindx[i+1]; /* ii-k=[i+1,k-1] */
 356           temp = 0.0;
 357           for (k=i+2; k<=j-1; k++) temp += qm[i+1][k-1]*qqm1[k];
 358           tt = rtype[type];
 359           qbt1 += temp * expMLclosing * exp_E_MLstem(tt, S1[j-1], S1[i+1], pf_params) * scale[2];
 360
 361           qb[i][j] = qbt1;
 362         } /* end if (type!=0) */
 363         else qb[i][j] = 0.0;
 364
 365         /* construction of qqm matrix containing final stem
 366            contributions to multiple loop partition function
 367            from segment i,j */
 368         qqm[i] = qqm1[i]*expMLbase[1];
 369         if (type) {
 370           qbt1 = qb[i][j] * exp_E_MLstem(type, (i>1) ? S1[i-1] : -1, (j<n) ? S1[j+1] : -1, pf_params);
 371           qqm[i] += qbt1;
 372         }
 373
 374         /*construction of qm matrix containing multiple loop
 375           partition function contributions from segment i,j */
 376         temp = 0.0;
 377         /*ii = my_iindx[i];   ii-k=[i,k-1] */
 378         /*new qm2 computation done here*/
 379         for (k=i+1; k<=j; k++) temp += (qm[i][k-1])*qqm[k];
 380         if (ulength>0) qm2[i][j]=temp;/*new qm2 computation done here*/
 381         for (k=i+1; k<=j; k++) temp += expMLbase[k-i] * qqm[k];
 382         qm[i][j] = (temp + qqm[i]);
 383
 384         /*auxiliary matrix qq for cubic order q calculation below */
 385         qbt1 = qb[i][j];
 386         if (type) {
 387           qbt1 *= exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j < n) ? S1[j+1] : -1, pf_params);
 388         }
 389         qq[i] = qq1[i]*scale[1] + qbt1;
 390
 391         /*construction of partition function for segment i,j */
 392         temp = 1.0*scale[1+j-i] + qq[i];
 393         for (k=i; k<=j-1; k++) temp += q[i][k]*qq[k+1];
 394         q[i][j] = temp;
 395
 396         if (temp>Qmax) {
 397           Qmax = temp;
 398           if (Qmax>max_real/10.)
 399             fprintf(stderr, "Q close to overflow: %d %d %g\n", i,j,temp);
 400         }
 401         if (temp>=max_real) {
 402           PRIVATE char msg[128];
 403           snprintf(msg, 128, "overflow in pf_fold while calculating q[%d,%d]\n"
 404                   "use larger pf_scale", i,j);
 405           nrerror(msg);
 406         }
 407       } /*end for i*/
 408       tmp = qq1;  qq1 =qq;  qq =tmp;
 409       tmp = qqm1; qqm1=qqm; qqm=tmp;
 410     }
 411
 412     /* just as a general service, I save here the free energy of the windows
 413        no output is generated, however,...
 414     */
 415     if ((j>=winSize) && (j<=n) && (ulength) && !(pUoutput)) {
 416       double Fwindow=0.;
 417       Fwindow=(-log(q[j-winSize+1][j])-winSize*log(pf_params->pf_scale))*pf_params->kT/1000.0;
 418
 419       pU[j][0]=Fwindow;
 420       /*
 421       if (ulength>=winSize)
 422         pU[j][winSize]=scale[winSize]/q[j-winSize+1][j];
 423       */
 424     }
 425     if (j>winSize) {
 426       Qmax=0;
 427       /* i=j-winSize; */
 428       /* initialize multiloopfs */
 429       for (k=j-winSize; k<=MIN2(n,j); k++) {
 430         prml[k]=0;
 431         prm_l[k]=0;
 432         /*        prm_l1[k]=0;  others stay*/
 433       }
 434       prm_l1[j-winSize]=0;
 435       k=j-winSize;
 436       for (l=k+TURN+1; l<=MIN2(n,k+winSize-1); l++) {
 437         int a;
 438         pR[k][l] = 0; /* set zero at start */
 439         type=ptype[k][l];
 440         if (qb[k][l]==0) continue;
 441
 442         for (a=MAX2(1,l-winSize+2); a<MIN2(k,n-winSize+2);a++)
 443           pR[k][l]+=q[a][k-1]*q[l+1][a+winSize-1]/q[a][a+winSize-1];
 444
 445         if (l-k+1==winSize)
 446           pR[k][l]+=1./q[k][l];
 447         else {
 448           if (k+winSize-1<=n)          /* k outermost */
 449             pR[k][l]+=q[l+1][k+winSize-1]/q[k][k+winSize-1];
 450           if (l-winSize+1>=1)  /*l outermost*/
 451             pR[k][l]+=q[l-winSize+1][k-1]/q[l-winSize+1][l];
 452         }
 453         pR[k][l] *= exp_E_ExtLoop(type, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params);
 454
 455         type_2 = ptype[k][l];
 456         type_2 = rtype[type_2];
 457
 458         for (i=MAX2(MAX2(l-winSize+1,k-MAXLOOP-1),1); i<=k-1; i++) {
 459           for (m=l+1; m<=MIN2(MIN2(l+ MAXLOOP -k+i+2,i+winSize-1),n); m++) {
 460             type = ptype[i][m];
 461             if ((pR[i][m]>0))
 462               pR[k][l] += pR[i][m]*exp_E_IntLoop(k-i-1, m-l-1, type, type_2,
 463                                                  S1[i+1], S1[m-1], S1[k-1], S1[l+1], pf_params) * scale[k-i+m-l];
 464           }
 465         }
 466         if (ulength) { /* NOT IF WITHIN INNER LOOP */
 467           for (i=MAX2(MAX2(l-winSize+1,k-MAXLOOP-1),1); i<=k-1; i++) {
 468             for (m=l+1; m<=MIN2(MIN2(l+ MAXLOOP -k+i+2,i+winSize-1),n); m++) {
 469               type = ptype[i][m];
 470               if ((pR[i][m]>0)){
 471                 temp=pR[i][m]*qb[k][l]*exp_E_IntLoop(k-i-1, m-l-1, type, type_2,
 472                                                      S1[i+1], S1[m-1], S1[k-1], S1[l+1], pf_params) * scale[k-i+m-l];
 473                 QI5[l][m-l-1]+=temp;
 474                 QI5[i][k-i-1]+=temp;
 475               }
 476             }
 477            }
 478         }
 479       }
 480       /* 3. bonding k,l as substem of multi-loop enclosed by i,m */
 481       prm_MLb = 0.;
 482       if(k>1) /*sonst nix!*/
 483         for (l=MIN2(n-1,k+winSize-2); l>=k+TURN+1; l--) { /* opposite direction */
 484           m=l+1;
 485           prmt = prmt1 = 0.0;
 486           tt = ptype[k-1][m]; tt=rtype[tt];
 487           prmt1 = pR[k-1][m] * expMLclosing * exp_E_MLstem(tt, S1[l], S1[k], pf_params);
 488           for (i=MAX2(1,l-winSize+2); i<k-1/*TURN*/; i++) {
 489             tt = ptype[i][m]; tt = rtype[tt];
 490             prmt += pR[i][m] * exp_E_MLstem(tt, S1[m-1], S1[i+1], pf_params) * qm[i+1][k-1];
 491           }
 492           tt = ptype[k][l];
 493           prmt *= expMLclosing;
 494           prml[ m] = prmt;
 495           prm_l[m] = prm_l1[m]*expMLbase[1]+prmt1;
 496
 497           prm_MLb = prm_MLb*expMLbase[1] + prml[m];
 498           /* same as:    prm_MLb = 0;
 499              for (i=n; i>k; i--)  prm_MLb += prml[i]*expMLbase[k-i-1];
 500           */
 501           prml[m] = prml[ m] + prm_l[m];
 502
 503           if (qb[k][l] == 0.) continue;
 504
 505           temp = prm_MLb;
 506
 507           if (ulength) {
 508             double dang;
 509             /* coefficient for computations of unpairedarrays */
 510             dang  =   qb[k][l] * exp_E_MLstem(tt, S1[k-1], S1[l+1], pf_params) * scale[2];
 511             for (m=MIN2(k+winSize-2,n);m>=l+2; m--){
 512               qmb[l][m-l-1] +=  prml[m]*dang;
 513               q2l[l][m-l-1] +=  (prml[m]-prm_l[m])*dang;
 514             }
 515           }
 516
 517           for (m=MIN2(k+winSize-2,n);m>=l+2; m--)
 518             temp += prml[m]*qm[l+1][m-1];
 519
 520           temp      *= exp_E_MLstem(tt, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params) * scale[2];
 521           pR[k][l]  += temp;
 522
 523           if (pR[k][l]>Qmax) {
 524             Qmax = pR[k][l];
 525             if (Qmax>max_real/10.)
 526               fprintf(stderr, "P close to overflow: %d %d %g %g\n",
 527                       i, m, pR[k][l], qb[k][l]);
 528           }
 529           if (pR[k][l]>=max_real) {
 530             ov++;
 531             pR[k][l]=FLT_MAX;
 532           }
 533
 534         } /* end for (l=..) */
 535       tmp = prm_l1; prm_l1=prm_l; prm_l=tmp;
 536
 537       /* end for (l=..)   */
 538       if ((ulength)&&(k-MAXLOOP-1>0)){
 539         /* if (pUoutput) pU[k-MAXLOOP-1]=(double *)space((ulength+2)*sizeof(double)); */
 540         if(split){ /*generate the new arrays, if you want them somewhere else, you have to generate them and overgive them ;)*/
 541           double **pUO;
 542           double **pUI;
 543           double **pUM;
 544           double **pUH;
 545           pUO= (double **)  space((n+1)*sizeof(double *));
 546           pUI= (double **)  space((n+1)*sizeof(double *));
 547           pUM= (double **)  space((n+1)*sizeof(double *));
 548           pUH= (double **)  space((n+1)*sizeof(double *));
 549           if (pUoutput) {
 550             for (i=1; i<=ulength; i++) {
 551               pUH[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 552               pUI[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 553               pUO[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 554               pUM[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 555             }
 556           }
 557           //dont want to have that yet?
 558           /*  else {
 559             for (i=1; i<=n; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
 560             }*/
 561           compute_pU_splitup(k-MAXLOOP-1,ulength,pU,pUO,pUH, pUI, pUM, winSize, n, sequence);
 562           if (pUoutput) {
 563             putoutpU_splitup(pUO,k-MAXLOOP-1, ulength, pUfp,'E');
 564             putoutpU_splitup(pUH,k-MAXLOOP-1, ulength, pUfp,'H');
 565             putoutpU_splitup(pUI,k-MAXLOOP-1, ulength, pUfp,'I');
 566             putoutpU_splitup(pUM,k-MAXLOOP-1, ulength, pUfp,'M');
 567           }
 568         }
 569         else {
 570         compute_pU(k-MAXLOOP-1,ulength,pU, winSize, n, sequence);
 571
 572         /* here, we put out and free pUs not in use any more (hopefully) */
 573         if (pUoutput)
 574           putoutpU(pU,k-MAXLOOP-1, ulength, pUfp);
 575       }
 576       }
 577
 578       if (j-(2*winSize+MAXLOOP+1)>0) {
 579         printpbar(pR,winSize,j-(2*winSize+MAXLOOP+1),n);
 580         if (simply_putout) {
 581           print_plist(n, j-(2*winSize+MAXLOOP+1), pR, winSize, spup);
 582         }
 583         else{
 584           pl=get_plistW(pl, n, j-(2*winSize+MAXLOOP+1), pR, winSize);
 585         }
 586         if (do_dpp)dpp=get_deppp(dpp,j-(2*winSize-MAXLOOP),pairSize, n);
 587         FreeOldArrays(j-(2*winSize+MAXLOOP+1));
 588       }
 589     }   /* end if (do_backtrack)*/
 590
 591   }/* end for j */
 592
 593   /* finish output and free */
 594   for (j=MAX2(1,n-MAXLOOP); j<=n;j++) {
 595     /* if (pUoutput) pU[j]=(double *)space((ulength+2)*sizeof(double)); */
 596     if (ulength) compute_pU(j,ulength,pU, winSize, n, sequence);
 597     /*here, we put out and free pUs not in use any more (hopefully)*/
 598     if (pUoutput) putoutpU(pU,j, ulength, pUfp);
 599   }
 600   for (j=MAX2(n-winSize-MAXLOOP,1); j<=n; j++) {
 601     printpbar(pR,winSize,j,n);
 602     if (simply_putout) {
 603       print_plist(n, j, pR, winSize, spup);
 604     }
 605     else {
 606       pl=get_plistW(pl, n, j, pR, winSize);
 607     }
 608     if ((do_dpp)&&j<n) dpp=get_deppp(dpp,j,pairSize, n);
 609     FreeOldArrays(j);
 610   }
 611   /* free_pf_arrays_L(); */
 612   free(S);
 613   free(S1);
 614   S = S1 = NULL;
 615   if (ov>0) fprintf(stderr, "%d overflows occurred while backtracking;\n"
 616                     "you might try a smaller pf_scale than %g\n",
 617                     ov, pf_params->pf_scale);
 618   *dpp2=dpp;
 619
 620   return pl;
 621 }
 622
 623 PRIVATE void scale_pf_params(unsigned int length, pf_paramT *parameters){
 624   unsigned int i;
 625   double  kT, scaling_factor;
 626
 627   if(pf_params) free(pf_params);
 628
 629   if(parameters){
 630     pf_params = get_boltzmann_factor_copy(parameters);
 631   } else {
 632     model_detailsT  md;
 633     set_model_details(&md);
 634     pf_params = get_boltzmann_factors(temperature, alpha, md, pf_scale);
 635   }
 636
 637   scaling_factor = pf_params->pf_scale;
 638   kT = pf_params->kT;   /* kT in cal/mol  */
 639
 640    /* scaling factors (to avoid overflows) */
 641   if (scaling_factor == -1) { /* mean energy for random sequences: 184.3*length cal */
 642     scaling_factor = exp(-(-185+(pf_params->temperature-37.)*7.27)/kT);
 643     if (scaling_factor<1) scaling_factor=1;
 644     pf_params->pf_scale = scaling_factor;
 645   }
 646   scale[0] = 1.;
 647   scale[1] = 1./scaling_factor;
 648   expMLbase[0] = 1;
 649   expMLbase[1] = pf_params->expMLbase/scaling_factor;
 650   for (i=2; i<=length; i++) {
 651     scale[i] = scale[i/2]*scale[i-(i/2)];
 652     expMLbase[i] = pow(pf_params->expMLbase, (double)i) * scale[i];
 653   }
 654 }
 655
 656 PRIVATE void printpbar(FLT_OR_DBL **prb,int winSize, int i, int n) {
 657   int j;
 658   int howoften=0; /* how many samples do we have for this pair */
 659   int pairdist;
 660
 661   for (j=i+TURN; j<MIN2(i+winSize,n+1); j++) {
 662     pairdist=(j-i+1);
 663     /*4cases*/
 664     howoften=MIN2(winSize-pairdist+1,i); /*pairdist,start*/
 665     howoften=MIN2(howoften,n-j+1);       /*end*/
 666     howoften=MIN2(howoften,n-winSize+1); /*windowsize*/
 667     prb[i][j] *= qb[i][j]/howoften;
 668   }
 669   return;
 670 }
 671
 672 PRIVATE void FreeOldArrays(int i) {
 673   /*free arrays no longer needed*/
 674   free(pR[i]+i);
 675   free(q[i]+i);
 676   free(qb[i]+i);
 677   free(qm[i]+i);
 678   if (ulength!=0) {
 679     free(qm2[i]+i);
 680     free(QI5[i]);
 681     free(qmb[i]);
 682     free(q2l[i]);
 683   }
 684   free(ptype[i]+i);
 685   return;
 686 }
 687
 688 PRIVATE void GetNewArrays(int j, int winSize) {
 689   /*allocate new part of arrays*/
 690   pR[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 691   pR[j]-=j;
 692   q[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 693   q[j]-=j;
 694   qb[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 695   qb[j]-=j;
 696   qm[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 697   qm[j]-=j;
 698   if (ulength!=0) {
 699     qm2[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 700     qm2[j]-=j;
 701     QI5[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 702     qmb[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 703     q2l[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
 704   }
 705   ptype[j]=(char *)space((winSize+1)*sizeof(char));
 706   ptype[j]-=j;
 707   return;
 708 }
 709
 710
 711 PRIVATE void GetPtype(int i, int winSize,const short *S,int n) {
 712   /*make new entries in ptype array*/
 713   int j;
 714   int type;
 715   for (j=i; j<=MIN2(i+winSize,n); j++) {
 716     type = pair[S[i]][S[j]];
 717     ptype[i][j] = (char) type;
 718   }
 719   return;
 720 }
 721
 722
 723 PRIVATE plist *get_plistW(plist *pl, int length,
 724                                  int start, FLT_OR_DBL **Tpr, int winSize) {
 725   /* get pair probibilities out of pr array */
 726   int  j,  max_p;
 727   max_p=1000;
 728   while (max_p<num_p)
 729     max_p*=2;
 730
 731   for (j=start+1; j<=MIN2(start+winSize, length); j++) {
 732     if (Tpr[start][j]<cutoff) continue;
 733     if (num_p==max_p-1) {
 734       max_p*=2;
 735       pl=(plist *)xrealloc(pl,max_p*sizeof(plist));
 736     }
 737     pl[num_p].i=start;
 738     pl[num_p].j=j;
 739     pl[num_p++].p=Tpr[start][j];
 740   }
 741
 742   /* mark end of data with zeroes */
 743   pl[num_p].i=0;
 744   pl[num_p].j=0;
 745   pl[num_p].p=0.;
 746   /* pl=(struct plist *)xrealloc(pl,(count)*sizeof(struct plist)); */
 747   return pl;
 748 }
 749
 750
 751 PRIVATE plist *get_deppp(plist *pl, int start, int pairsize, int length) {
 752   /* compute dependent pair probabilities */
 753   int i, j, count=0;
 754   double tmp;
 755   plist *temp;
 756   temp=(plist *)space(pairsize*sizeof(plist)); /* holds temporary deppp */
 757   for (j=start+TURN; j<MIN2(start+pairsize,length); j++) {
 758
 759     if ((qb[start][j]*qb[start-1][(j+1)])>10e-200) {
 760       int type=ptype[start-1][j+1];
 761       int type_2=rtype[ptype[start][j]];
 762       tmp=qb[start][j]/qb[start-1][(j+1)]*exp_E_IntLoop(0, 0, type, type_2,
 763                                                         S1[start], S1[j], S1[start-1], S1[j+1], pf_params) * scale[2];
 764        temp[count].i=start;
 765       temp[count].j=j;
 766       temp[count++].p=tmp;
 767     }
 768   }
 769   /* write it to list of deppps */
 770   for (i=0; pl[i].i!=0; i++);
 771   pl=(plist *)xrealloc(pl,(i+count+1)*sizeof(plist));
 772   for (j=0; j<count; j++) {
 773     pl[i+j].i=temp[j].i;
 774     pl[i+j].j=temp[j].j;
 775     pl[i+j].p=temp[j].p;
 776   }
 777   pl[i+count].i=0;
 778   pl[i+count].j=0;
 779   pl[i+count].p=0;
 780   free(temp);
 781   return pl;
 782 }
 783
 784
 785 PRIVATE void print_plist(int length,int start, FLT_OR_DBL **Tpr, int winSize, FILE *fp) {
 786   /* print out of pr array, do not save */
 787   int  j;
 788
 789
 790   for (j=start+1; j<=MIN2(start+winSize, length); j++) {
 791     if (Tpr[start][j]<cutoff) continue;
 792     fprintf(fp,"%d  %d  %g\n",start,j,Tpr[start][j]);
 793   }
 794
 795   /* mark end of data with zeroes */
 796
 797   return ;
 798 }
 799
 800 PRIVATE void compute_pU(int k, int ulength, double **pU, int winSize,int n, char *sequence) {
 801 /*  here, we try to add a function computing all unpaired probabilities starting at some i,
 802     going down to $unpaired, to be unpaired, i.e. a list with entries from 1 to unpaired for
 803     every i, with the probability of a stretch of length x, starting at i-x+1, to be unpaired
 804 */
 805   int startu;
 806   int i5;
 807   int j3, len, obp;
 808   double temp;
 809   double *QBE;
 810   FLT_OR_DBL  expMLclosing      = pf_params->expMLclosing;
 811
 812   QBE=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
 813
 814   /* first, we will */
 815   /* for k<=ulength, pU[k][k]=0, because no bp can enclose it */
 816   if (pUoutput&&k+ulength<=n)  pU[k+ulength]=(double *)space((ulength+2)*sizeof(double));
 817   /*compute pu[k+ulength][ulength] */
 818    for (i5=MAX2(k+ulength-winSize+1,1);i5<=k;i5++) {
 819     for (j3=k+ulength+1; j3<=MIN2(n,i5+winSize-1); j3++) {
 820       /*  if (k>400) {
 821         printf("i%d j%d  ",i5,j3);
 822         fflush(stdout);
 823         } */
 824       if (ptype[i5][j3]!=0) {/**/
 825         /* (.. >-----|..........)
 826           i5  j     j+ulength  j3              */
 827         /*Multiloops*/
 828         temp = (i5<k) ? qm2[i5+1][k] * expMLbase[j3-k-1] : 0.; /* (..{}{}-----|......) */
 829
 830         if(j3-1>k+ulength)
 831           temp  +=  qm2[k+ulength+1][j3-1] * expMLbase[k+ulength-i5]; /* (..|-----|{}{}) */
 832
 833         if((i5<k)&&(j3-1>k+ulength))
 834           temp  +=  qm[i5+1][k] * qm[k+ulength+1][j3-1] * expMLbase[ulength]; /* ({}|-----|{}) */
 835
 836         /* add dangles, multloopclosing etc. */
 837         temp  *=  exp_E_MLstem(rtype[ptype[i5][j3]], S1[j3-1], S1[i5+1], pf_params) * scale[2] * expMLclosing;
 838         /*add hairpins*/
 839         temp  +=  exp_E_Hairpin(j3-i5-1, ptype[i5][j3], S1[i5+1], S1[j3-1], sequence+i5-1, pf_params) * scale[j3-i5+1];
 840         /*add outer probability*/
 841         temp *= pR[i5][j3];
 842         pU[k+ulength][ulength] += temp;
 843
 844       }
 845     }
 846    }
 847    /* code doubling to avoid if within loop */
 848 #if 0
 849   /*initialization for interior loops,
 850     it is not recomended to have verysmall ulengths!!*/
 851   if (ulength<MAXLOOP) {
 852     int k5;
 853     int l3;
 854     int outype;
 855     /* kl bp is 5' */
 856     /* MAXLOOP>((l5-k5-1)+(j3-l3-1)
 857       k-winSize+ulength<i5<k-TURN-1;
 858       k+ulength<j3<=k+MAXLOOP+1
 859       if i then use l3, it is easier by far:
 860       j3-MAXLOOP<=l3<=k
 861       i5<k5<k-TURN k5<=i5+l3+2+MAXLOOP-j3
 862       k5+TURN<l3<=k
 863     */
 864     for (i5=MAX2(k+ulength-winSize,1);i5<k-TURN-1;i5++) {
 865
 866       for (j3=k+ulength+1; j3<=MIN2(n,MIN2(i5+winSize-1,k+MAXLOOP+1)); j3++) {
 867         double temp=0;
 868         if (outype=ptype[i5][j3]>0) /* oder so halt */
 869           for (l3=MAX2(i5+TURN+1,j3-MAXLOOP-1); l3<=k; l3++){
 870             for (k5=i5+1; k5<=MIN2(l3-TURN-1,MAXLOOP+i5+l3+2-j3); k5++){
 871               if (ptype[k5][l3]) {
 872                 temp+= qb[k5][l3]*expLoopEnergy(k5-i5-1, j3-l3-1, outype, rtype[ptype[k5][l3]], S1[i5+1], S1[j3-1], S1[k5-1], S1[l3+1]);
 873               }
 874             }
 875           }
 876         temp*=pR[i5][j3];
 877         pU[k+ulength][ulength]+= temp;
 878       }
 879     }
 880     /* kl bp is 3' */
 881     /*
 882       k+ulength-MAXLOOP<=i5<=k
 883       k+ulength+1+TURN<j3<i5+winSize
 884       k+ulength+1<=k5<i5+MAXLOOP+2 || k5<j3-TURN
 885       k5<l3<j3 || j3-k5-i5-2-ML<=l3<j3
 886     */
 887     for (i5=MAX2(1,MAX2(k+ulength-winSize,k+ulength-MAXLOOP));i5<=k; i5++){
 888       for (j3=k+ulength+TURN+2; j3<MIN2(n+1,i5+winSize); j3++) {
 889         double temp = 0;
 890         if (outype=ptype[i5][j3]>0) /* oder so halt */
 891           for (k5=k+ulength+1; k5<MIN2(j3-TURN-1,i5+MAXLOOP+2); k5++) {
 892             for (l3=MAX2(k5+TURN+1,j3+k5-i5-2-MAXLOOP); l3<j3; l3++) {
 893               if (ptype[k5][l3])
 894                 temp += qb[k5][l3]*expLoopEnergy(k5-i5-1, j3-l3-1, outype, rtype[ptype[k5][l3]], S1[i5+1], S1[j3-1], S1[k5-1], S1[l3+1]);
 895             }
 896           }
 897         temp*=pR[i5][j3];
 898         pU[k+ulength][ulength]+= temp;
 899       }
 900     }
 901   }
 902   /* Add up Is QI5[l][m-l-1] QI3 */
 903   /* Add up Interior loop terms */
 904   temp=0.;
 905
 906   for (len=winSize; len>=ulength; len--) temp+=QI3[k][len];
 907   for (;len>0; len--) {
 908     temp += QI3[k][len];
 909     QBE[len] += temp;
 910   }
 911 #endif
 912   temp=0.;
 913   for (len=winSize; len>=MAX2(ulength,MAXLOOP); len--) temp+=QI5[k][len];
 914   for (;len>0; len--) {
 915     temp += QI5[k][len];
 916     QBE[len] += temp;  /* replace QBE with QI */
 917   }
 918   /* Add Hairpinenergy to QBE */
 919   temp=0.;
 920   for(obp = MIN2(n, k + winSize - 1); obp > k + ulength; obp--)
 921     if(ptype[k][obp])
 922       temp += pR[k][obp] * exp_E_Hairpin(obp-k-1, ptype[k][obp], S1[k+1], S1[obp-1], sequence+k-1, pf_params) * scale[obp-k+1];
 923   for(obp = MIN2(n, MIN2(k + winSize - 1, k + ulength)); obp > k + 1; obp--){
 924     if (ptype[k][obp])
 925       temp += pR[k][obp] * exp_E_Hairpin(obp-k-1, ptype[k][obp], S1[k+1], S1[obp-1], sequence+k-1, pf_params) * scale[obp-k+1];
 926     QBE[obp-k-1] += temp;  /* add hairpins to QBE (all in one array) */
 927   }
 928   /* doubling the code to get the if out of the loop */
 929
 930   /* Add up Multiloopterms  qmb[l][m]+=prml[m]*dang;
 931     q2l[l][m]+=(prml[m]-prm_l[m])*dang; */
 932
 933   temp=0.;
 934   for(len = winSize; len >= ulength; len--)
 935     temp += q2l[k][len] * expMLbase[len];
 936   for( ; len > 0; len--){
 937     temp += q2l[k][len] * expMLbase[len];
 938     QBE[len] += temp; /* add (()()____) type cont. to I3 */
 939   }
 940   for(len = 1; len < ulength; len++){
 941     for(obp = k + len + TURN; obp <= MIN2(n, k + winSize - 1); obp++){
 942       /* add (()___()) */
 943       QBE[len] += qmb[k][obp-k-1] * qm[k+len+1/*2*/][obp-1] * expMLbase[len];
 944     }
 945   }
 946   for (len=1; len<ulength; len++) {
 947     for (obp=k+len+TURN+TURN; obp<=MIN2(n,k+winSize-1); obp++) {
 948       if (ptype[k][obp]) {
 949         temp      =   exp_E_MLstem(rtype[ptype[k][obp]], S1[obp-1], S1[k+1], pf_params) * scale[2] * expMLbase[len] * expMLclosing; /* k:obp */
 950         QBE[len]  +=  pR[k][obp] * temp * qm2[k+len+1][obp-1]; /* add (___()()) */
 951       }
 952     }
 953   }
 954   /* After computing all these contributions in QBE[len], that k is paired
 955     and the unpaired stretch is AT LEAST len long, we start to add that to
 956     the old unpaired thingies; */
 957   for(len = 1; len < MIN2(MAX2(ulength, MAXLOOP), n - k); len++){
 958     pU[k+len][len] += pU[k+len][len+1] + QBE[len];
 959   }
 960
 961   /*open chain*/
 962   if ((ulength>=winSize)&&(k>=ulength)) {
 963     pU[k][winSize]=scale[winSize]/q[k-winSize+1][k];
 964   }
 965   /* now the not enclosed by any base pair terms for whatever it is we do not need anymore...
 966     ... which should be e.g; k, again */
 967   for(startu = MIN2(ulength, k); startu > 0; startu--){
 968     temp=0.;
 969     for(i5 = MAX2(1, k - winSize + 2); i5 <= MIN2(k - startu, n - winSize + 1); i5++){
 970       temp += q[i5][k - startu] * q[k + 1][i5 + winSize - 1] * scale[startu]/q[i5][i5 + winSize - 1];
 971     }
 972     /* the 2 Cases where the borders are on the edge of the interval */
 973     if((k >= winSize) && (startu + 1 <= winSize))
 974       temp += q[k - winSize + 1][k - startu]*scale[startu]/q[k - winSize + 1][k];
 975     if((k <= n - winSize+ startu) && (k - startu >= 0) && (k < n) && (startu + 1 <= winSize))
 976       temp += q[k + 1][k - startu + winSize] * scale[startu] / q[k - startu + 1][k - startu + winSize];
 977
 978     /* Divide by number of possible windows */
 979     pU[k][startu] += temp;
 980     {
 981       int leftmost, rightmost;
 982
 983       leftmost      = MAX2(1, k - winSize + 1);
 984       rightmost     = MIN2(n - winSize + 1, k - startu + 1);
 985       pU[k][startu] /= (rightmost - leftmost + 1);
 986     }
 987   }
 988   free(QBE);
 989   return;
 990 }
 991
 992
 993 PRIVATE void putoutpU(double **pUx, int k, int ulength, FILE *fp) {
 994   /*put out unpaireds for k, and free pU[k], make sure we don't need pU[k] any more!!*/
 995   /*could use that for hairpins, also!*/
 996   int i;
 997   fprintf(fp,"%d\t",k);
 998   for (i=1; i<=MIN2(ulength,k); i++) {
 999     fprintf(fp,"%.5g\t",pUx[k][i]);
1000   }
1001   fprintf(fp,"\n");
1002   free(pUx[k]);
1003 }
1004 PRIVATE void putoutpU_splitup(double **pUx, int k, int ulength, FILE *fp, char ident) {
1005   /*put out unpaireds for k, and free pU[k], make sure we don't need pU[k] any more!!*/
1006   /*could use that for hairpins, also!*/
1007   int i;
1008   fprintf(fp,"%d\t",k);
1009   for (i=1; i<=MIN2(ulength,k); i++) {
1010     fprintf(fp,"%.5g\t",pUx[k][i]);
1011   }
1012   fprintf(fp,"\t%s\n",ident);
1013   free(pUx[k]);
1014 }
1015
1016 PUBLIC void putoutpU_prob(double **pU,int length, int ulength, FILE *fp, int energies) {
1017   putoutpU_prob_par(pU, length, ulength, fp, energies, pf_params);
1018 }
1019
1020
1021 PUBLIC void putoutpU_prob_par(double **pU,int length, int ulength, FILE *fp, int energies, pf_paramT *parameters){
1022   /*put out unpaireds */
1023   int i,k;
1024   double kT = parameters->kT/1000.0;
1025   double temp;
1026   if (energies) fprintf(fp,"#opening energies\n #i$\tl=");
1027   else  fprintf(fp,"#unpaired probabilities\n #i$\tl=");
1028   for (i=1; i<=ulength; i++) {
1029     fprintf(fp,"%d\t", i);
1030   }
1031   fprintf(fp,"\n");
1032
1033   for (k=1; k<=length; k++){
1034     fprintf(fp,"%d\t",k);
1035     for (i=1; i<=ulength; i++) {
1036       if (i>k) {
1037         fprintf(fp,"NA\t");
1038         continue;
1039       }
1040       if (energies) temp=-log(pU[k][i])*kT;
1041       else temp=pU[k][i];
1042       fprintf(fp,"%.7g\t",temp);
1043     }
1044     fprintf(fp,"\n");
1045     free(pU[k]);
1046   }
1047   fflush(fp);
1048 }
1049
1050 PUBLIC void putoutpU_prob_bin(double **pU,int length, int ulength, FILE *fp, int energies) {
1051   putoutpU_prob_bin_par(pU, length, ulength, fp, energies, pf_params);
1052 }
1053
1054 PUBLIC void putoutpU_prob_bin_par(double **pU,int length, int ulength, FILE *fp, int energies, pf_paramT *parameters) {
1055
1056   /*put out unpaireds */
1057   int i,k;
1058   double kT= parameters->kT/1000.0;
1059   double temp;
1060   int *p;
1061   p = (int*) space(sizeof(int)*1);
1062   /* write first line */
1063   p[0]=ulength; /* u length */
1064   fwrite(p,sizeof(int),1,fp);
1065   p[0]=length; /* seq length */
1066   fwrite(p,sizeof(int),1,fp);
1067   for (k=3; k<=(length+20); k++){ /* all the other lines are set to 1000000 because we are at ulength=0 */
1068     p[0]=1000000;
1069     fwrite(p,sizeof(int),1,fp);
1070   }
1071   /* data */
1072   for (i=1; i<=ulength; i++) {
1073     for (k=1; k<=11; k++){/* write first ten entries to 1000000 */
1074       p[0]=1000000;
1075       fwrite(p,sizeof(int),1,fp);
1076     }
1077     for (k=1; k<=length; k++){/* write data now */
1078       if (i>k) {
1079         p[0]=1000000;         /* check if u > pos */
1080         fwrite(p,sizeof(int),1,fp);
1081         continue;
1082       }
1083       else{
1084         p[0]= (int) rint(100 *(-log(pU[k][i])*kT));
1085         fwrite(p,sizeof(int),1,fp);
1086       }
1087     }
1088     for (k=1; k<=9; k++){/* finish by writing the last 10 entries */
1089       p[0]=1000000;
1090       fwrite(p,sizeof(int),1,fp);
1091     }
1092   }
1093   /* free pU array; */
1094   for (k=1; k<=length; k++){
1095     free(pU[k]);
1096   }
1097   free(p);
1098   fflush(fp);
1099 }
1100
1101
1102 /*
1103  Here: Space for questions...
1104 */
1105 PRIVATE void compute_pU_splitup(int k, int ulength, double **pU,  double **pUO, double **pUH, double **pUI, double **pUM, int winSize,int n, char *sequence) {
1106 /*  here, we try to add a function computing all unpaired probabilities starting at some i,
1107     going down to $unpaired, to be unpaired, i.e. a list with entries from 1 to unpaired for
1108     every i, with the probability of a stretch of length x, starting at i-x+1, to be unpaired
1109 */
1110   int startu;
1111   int i5;
1112   int j3, len, obp;
1113   double temp;
1114   double *QBE;
1115   double *QBI;
1116   double *QBM;
1117   double *QBH;
1118
1119   FLT_OR_DBL  expMLclosing      = pf_params->expMLclosing;
1120
1121   QBE=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
1122   QBM=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
1123   QBI=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
1124   QBH=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
1125
1126   /* first, we will */
1127   /* for k<=ulength, pU[k][k]=0, because no bp can enclose it */
1128   if (pUoutput&&k+ulength<=n)  pU[k+ulength]=(double *)space((ulength+2)*sizeof(double));
1129   /*compute pu[k+ulength][ulength] */
1130    for (i5=MAX2(k+ulength-winSize+1,1);i5<=k;i5++) {
1131     for (j3=k+ulength+1; j3<=MIN2(n,i5+winSize-1); j3++) {
1132       /*  if (k>400) {
1133         printf("i%d j%d  ",i5,j3);
1134         fflush(stdout);
1135         } */
1136       if (ptype[i5][j3]!=0) {/**/
1137         /* (.. >-----|..........)
1138           i5  j     j+ulength  j3              */
1139         /*Multiloops*/
1140         temp = (i5<k) ? qm2[i5+1][k] * expMLbase[j3-k-1] : 0.; /* (..{}{}-----|......) */
1141
1142         if(j3-1>k+ulength)
1143           temp  +=  qm2[k+ulength+1][j3-1] * expMLbase[k+ulength-i5]; /* (..|-----|{}{}) */
1144
1145         if((i5<k)&&(j3-1>k+ulength))
1146           temp  +=  qm[i5+1][k] * qm[k+ulength+1][j3-1] * expMLbase[ulength]; /* ({}|-----|{}) */
1147
1148         /* add dangles, multloopclosing etc. */
1149         temp  *=  exp_E_MLstem(rtype[ptype[i5][j3]], S1[j3-1], S1[i5+1], pf_params) * scale[2] * expMLclosing;
1150         /*add hairpins*/
1151         temp  +=  exp_E_Hairpin(j3-i5-1, ptype[i5][j3], S1[i5+1], S1[j3-1], sequence+i5-1, pf_params) * scale[j3-i5+1];
1152         /*add outer probability*/
1153         temp *= pR[i5][j3];
1154         pU[k+ulength][ulength] += temp;
1155
1156       }
1157     }
1158    }
1159    /* code doubling to avoid if within loop */
1160   temp=0.;
1161   for (len=winSize; len>=MAX2(ulength,MAXLOOP); len--) temp+=QI5[k][len];
1162   for (;len>0; len--) {
1163     temp += QI5[k][len];
1164     QBI[len] += temp;
1165     QBE[len] += temp;  /* replace QBE with QI */
1166   }
1167   /* Add Hairpinenergy to QBE */
1168   temp=0.;
1169   for(obp = MIN2(n, k + winSize - 1); obp > k + ulength; obp--)
1170     if(ptype[k][obp])
1171       temp += pR[k][obp] * exp_E_Hairpin(obp-k-1, ptype[k][obp], S1[k+1], S1[obp-1], sequence+k-1, pf_params) * scale[obp-k+1];
1172   for(obp = MIN2(n, MIN2(k + winSize - 1, k + ulength)); obp > k + 1; obp--){
1173     if (ptype[k][obp])
1174       temp += pR[k][obp] * exp_E_Hairpin(obp-k-1, ptype[k][obp], S1[k+1], S1[obp-1], sequence+k-1, pf_params) * scale[obp-k+1];
1175     QBH[obp-k-1] += temp;
1176     QBE[obp-k-1] += temp;  /* add hairpins to QBE (all in one array) */
1177   }
1178   /* doubling the code to get the if out of the loop */
1179
1180   /* Add up Multiloopterms  qmb[l][m]+=prml[m]*dang;
1181     q2l[l][m]+=(prml[m]-prm_l[m])*dang; */
1182
1183   temp=0.;
1184   for(len = winSize; len >= ulength; len--)
1185     temp += q2l[k][len] * expMLbase[len];
1186   for( ; len > 0; len--){
1187     temp += q2l[k][len] * expMLbase[len];
1188     QBM[len] += temp;
1189     QBE[len] += temp; /* add (()()____) type cont. to I3 */
1190   }
1191   for(len = 1; len < ulength; len++){
1192     for(obp = k + len + TURN; obp <= MIN2(n, k + winSize - 1); obp++){
1193       /* add (()___()) */
1194       QBM[len] += qmb[k][obp-k-1] * qm[k+len+1/*2*/][obp-1] * expMLbase[len];
1195       QBE[len] += qmb[k][obp-k-1] * qm[k+len+1/*2*/][obp-1] * expMLbase[len];
1196     }
1197   }
1198   for (len=1; len<ulength; len++) {
1199     for (obp=k+len+TURN+TURN; obp<=MIN2(n,k+winSize-1); obp++) {
1200       if (ptype[k][obp]) {
1201         temp      =   exp_E_MLstem(rtype[ptype[k][obp]], S1[obp-1], S1[k+1], pf_params) * scale[2] * expMLbase[len] * expMLclosing; /* k:obp */
1202         QBE[len]  +=  pR[k][obp] * temp * qm2[k+len+1][obp-1]; /* add (___()()) */
1203         QBM[len]  +=  pR[k][obp] * temp * qm2[k+len+1][obp-1]; /* add (___()()) */
1204       }
1205     }
1206   }
1207   /* After computing all these contributions in QBE[len], that k is paired
1208     and the unpaired stretch is AT LEAST len long, we start to add that to
1209     the old unpaired thingies; */
1210   for(len = 1; len < MIN2(MAX2(ulength, MAXLOOP), n - k); len++){
1211     pU[k+len][len] += pU[k+len][len+1] + QBE[len];
1212     pUH[k+len][len] += pUH[k+len][len+1] + QBH[len];
1213     pUM[k+len][len] += pUM[k+len][len+1] + QBM[len];
1214     pUI[k+len][len] += pUI[k+len][len+1] + QBI[len];
1215
1216   }
1217
1218   /* open chain */
1219   if ((ulength>=winSize)&&(k>=ulength)) {
1220     pUO[k][winSize]=scale[winSize]/q[k-winSize+1][k];
1221   }
1222   /*open chain*/
1223   if ((ulength>=winSize)&&(k>=ulength)) {
1224     pU[k][winSize]=scale[winSize]/q[k-winSize+1][k];
1225   }
1226   /* now the not enclosed by any base pair terms for whatever it is we do not need anymore...
1227     ... which should be e.g; k, again */
1228   for(startu = MIN2(ulength, k); startu > 0; startu--){
1229     temp=0.;
1230     for(i5 = MAX2(1, k - winSize + 2); i5 <= MIN2(k - startu, n - winSize + 1); i5++){
1231       temp += q[i5][k - startu] * q[k + 1][i5 + winSize - 1] * scale[startu]/q[i5][i5 + winSize - 1];
1232     }
1233     /* the 2 Cases where the borders are on the edge of the interval */
1234     if((k >= winSize) && (startu + 1 <= winSize))
1235       temp += q[k - winSize + 1][k - startu]*scale[startu]/q[k - winSize + 1][k];
1236     if((k <= n - winSize+ startu) && (k - startu >= 0) && (k < n) && (startu + 1 <= winSize))
1237       temp += q[k + 1][k - startu + winSize] * scale[startu] / q[k - startu + 1][k - startu + winSize];
1238
1239     /* Divide by number of possible windows */
1240     pU[k][startu] += temp;
1241     pUO[k][startu] += temp;
1242
1243     {
1244       int leftmost, rightmost;
1245
1246       leftmost      = MAX2(1, k - winSize + 1);
1247       rightmost     = MIN2(n - winSize + 1, k - startu + 1);
1248       pU[k][startu] /= (rightmost - leftmost + 1);
1249       /*Do we want to make a distinction between those?*/
1250       pUH[k][startu] /= (rightmost - leftmost + 1);
1251       pUO[k][startu] /= (rightmost - leftmost + 1);
1252       pUI[k][startu] /= (rightmost - leftmost + 1);
1253       pUM[k][startu] /= (rightmost - leftmost + 1);
1254     }
1255   }
1256   free(QBE);
1257   free(QBI);
1258   free(QBH);
1259   free(QBM);
1260   return;
1261 }
1262 PUBLIC void putoutpU_prob_splitup(double **pU, double **pUO, double **pUH, double **pUI, double **pUM, int length, int ulength, FILE *fp, int energies) {
1263   /*put out unpaireds */
1264   int i,k;
1265   double kT= (temperature+K0)*GASCONST/1000.0;
1266   double temp;
1267   if (energies) fprintf(fp,"#opening energies\n #i$\tl=");
1268   else  fprintf(fp,"#unpaired probabilities\n #i$\tl=");
1269
1270   fprintf(fp,"Total\n");
1271   for (i=1; i<=ulength; i++) {
1272     fprintf(fp,"%d\t", i);
1273   }
1274   fprintf(fp,"\n");
1275
1276   for (k=1; k<=length; k++){
1277     fprintf(fp,"%d\t",k);
1278     for (i=1; i<=ulength; i++) {
1279       if (i>k) {
1280         fprintf(fp,"NA\t");
1281         continue;
1282       }
1283       if (energies) temp=-log(pU[k][i])*kT;
1284       else temp=pU[k][i];
1285       fprintf(fp,"%.7g\t",temp);
1286     }
1287     fprintf(fp,"\tT\n");
1288     free(pU[k]);
1289   }
1290   fprintf(fp,"\n###################################################################\nHairpin\n");
1291   for (i=1; i<=ulength; i++) {
1292     fprintf(fp,"%d\t", i);
1293   }
1294   fprintf(fp,"\n");
1295
1296   for (k=1; k<=length; k++){
1297     fprintf(fp,"%d\t",k);
1298     for (i=1; i<=ulength; i++) {
1299       if (i>k) {
1300         fprintf(fp,"NA\t");
1301         continue;
1302       }
1303       if (energies) temp=-log(pUH[k][i])*kT;
1304       else temp=pUH[k][i];
1305       fprintf(fp,"%.7g\t",temp);
1306     }
1307     fprintf(fp,"\tH\n");
1308     free(pUH[k]);
1309   }
1310   fprintf(fp,"\n###################################################################\nInterior\n");
1311   for (i=1; i<=ulength; i++) {
1312     fprintf(fp,"%d\t", i);
1313   }
1314   fprintf(fp,"\n");
1315
1316   for (k=1; k<=length; k++){
1317     fprintf(fp,"%d\t",k);
1318     for (i=1; i<=ulength; i++) {
1319       if (i>k) {
1320         fprintf(fp,"NA\t");
1321         continue;
1322       }
1323       if (energies) temp=-log(pUI[k][i])*kT;
1324       else temp=pUI[k][i];
1325       fprintf(fp,"%.7g\t",temp);
1326     }
1327     fprintf(fp,"\tI\n");
1328     free(pUI[k]);
1329   }
1330   fprintf(fp,"\n###################################################################\nMultiloop\n");
1331   for (i=1; i<=ulength; i++) {
1332     fprintf(fp,"%d\t", i);
1333   }
1334   fprintf(fp,"\n");
1335
1336   for (k=1; k<=length; k++){
1337     fprintf(fp,"%d\t",k);
1338     for (i=1; i<=ulength; i++) {
1339       if (i>k) {
1340         fprintf(fp,"NA\t");
1341         continue;
1342       }
1343       if (energies) temp=-log(pUM[k][i])*kT;
1344       else temp=pUM[k][i];
1345       fprintf(fp,"%.7g\t",temp);
1346     }
1347     fprintf(fp,"\tM\n");
1348     free(pUM[k]);
1349   }
1350   fprintf(fp,"\n###################################################################\nExterior\n");
1351   for (i=1; i<=ulength; i++) {
1352     fprintf(fp,"%d\t", i);
1353   }
1354   fprintf(fp,"\t E\n");
1355
1356   for (k=1; k<=length; k++){
1357     fprintf(fp,"%d\t",k);
1358     for (i=1; i<=ulength; i++) {
1359       if (i>k) {
1360         fprintf(fp,"NA\t");
1361         continue;
1362       }
1363       if (energies) temp=-log(pUO[k][i])*kT;
1364       else temp=pUO[k][i];
1365       fprintf(fp,"%.7g\t",temp);
1366     }
1367     fprintf(fp,"\n");
1368     free(pU[k]);
1369   }
1370   fflush(fp);
1371 }
1372
1373
1374 /*###########################################*/
1375 /*# deprecated functions below              #*/
1376 /*###########################################*/
1377
1378 PUBLIC void init_pf_foldLP(int length){ /* DO NOTHING */}
1379