binaries/src/ViennaRNA/lib/part_func_up.c

   1 /* Last changed Time-stamp: <2008-07-04 15:57:03 ulim> */
   2 /*
   3                   partiton function for RNA secondary structures
   4
   5                   Ivo L Hofacker
   6
   7                   Vienna RNA package
   8 */
   9 /*
  10   $Log: part_func_up.c,v $
  11   Revision 1.4  2008/07/04 14:27:36  ivo
  12   Modify output (again)
  13
  14   Revision 1.3  2008/05/08 14:11:55  ivo
  15   minor output changes
  16
  17   Revision 1.2  2007/12/13 10:19:54  ivo
  18   major RNAup update from Ulli
  19
  20   Revision 1.1  2007/04/30 15:13:13  ivo
  21   merge RNAup into package
  22
  23   Revision 1.11  2006/07/17 11:11:43  ulim
  24   removed all globals from fold_vars.h,c, cleaned code
  25
  26   Revision 1.10  2006/07/12 09:19:29  ulim
  27   global variables w, incr3 and incr5 are now local
  28
  29   Revision 1.9  2006/07/11 12:45:02  ulim
  30   remove redundancy in function pf_interact(...)
  31
  32   Revision 1.8  2006/03/08 15:26:37  ulim
  33   modified -o[1|2], added meaningful default
  34
  35   Revision 1.5  2006/01/23 11:27:04  ulim
  36   include file into new package version. cleaned it
  37
  38   Revision 1.2  2005/07/29 15:13:37  ulim
  39   put the function, calculating the probability of an unpaired region in
  40   an RNA and the function calculating the prob. of interaction between 2 RNAs
  41   in a seperate file (pf_two.c)
  42
  43   Revision 1.1  2005/07/26 13:27:12  ulim
  44   Initial revision
  45
  46   Revision 1.2  2005/07/01 13:14:57  ulim
  47   fixed error in scaling, included new commandline options -incr5, -incr3 to
  48   allow a variable number of unpaired positions 5' and 3' of the site of
  49   interaction between the two RNAs
  50
  51   Revision 1.1  2005/04/19 08:16:38  ulim
  52   Initial revision
  53 */
  54
  55 #include <config.h>
  56 #include <stdio.h>
  57 #include <stdlib.h>
  58 #include <string.h>
  59 #include <math.h>
  60 #include <float.h>    /* #defines FLT_MAX ... */
  61 #include <unistd.h>
  62 #include "fold.h"
  63 #include "utils.h"
  64 #include "energy_par.h"
  65 #include "fold_vars.h"
  66 #include "pair_mat.h"
  67 #include "params.h"
  68 #include "part_func.h"
  69 #include "loop_energies.h"
  70 #include "part_func_up.h"
  71 #include "duplex.h"
  72 /*@unused@*/
  73 static char rcsid[] UNUSED = "$Id: part_func_up.c,v 1.4 2008/07/04 14:27:36 ivo Exp $";
  74
  75 #define CO_TURN 0
  76 #define ZERO(A) (fabs(A) < DBL_EPSILON)
  77 #define EQUAL(A,B) (fabs((A)-(B)) < 1000*DBL_EPSILON)
  78 #define ISOLATED  256.0
  79 /* #define NUMERIC 1 */
  80
  81 /*
  82 #################################
  83 # GLOBAL VARIABLES              #
  84 #################################
  85 */
  86
  87 /*
  88 #################################
  89 # PRIVATE VARIABLES             #
  90 #################################
  91 */
  92 PRIVATE short       *S=NULL, *S1=NULL, *SS=NULL, *SS2=NULL;
  93 PRIVATE pf_paramT   *Pf = NULL;/* use this structure for all the exp-arrays*/
  94 PRIVATE FLT_OR_DBL  *qb=NULL, *qm=NULL, *prpr=NULL; /* add arrays for pf_unpaired()*/
  95 PRIVATE FLT_OR_DBL  *probs=NULL;
  96 PRIVATE FLT_OR_DBL  *q1k=NULL, *qln=NULL;
  97 PRIVATE double      *qqm2=NULL, *qq_1m2=NULL, *qqm=NULL, *qqm1=NULL;
  98 PRIVATE FLT_OR_DBL  *scale=NULL, *expMLbase=NULL;
  99 PRIVATE char        *ptype=NULL; /* precomputed array of pair types */
 100 PRIVATE int         init_length;  /* length in last call to init_pf_fold()*/
 101 PRIVATE double      init_temp; /* temperature in last call to scale_pf_params */
 102 PRIVATE int         *my_iindx = NULL;
 103 /* make iptypes array for intermolecular constrains (ipidx for indexing)*/
 104
 105
 106 /*
 107 #################################
 108 # PRIVATE FUNCTION DECLARATIONS #
 109 #################################
 110 */
 111 PRIVATE pu_out      *get_u_vals(pu_contrib *p_c,
 112                                 int **unpaired_values,
 113                                 char *select_contrib);
 114
 115 PRIVATE int         plot_free_pu_out( pu_out* res,
 116                                       interact *pint,
 117                                       char *ofile,
 118                                       char *head);
 119
 120 PRIVATE void        scale_stru_pf_params(unsigned int length);
 121
 122 PRIVATE void        init_pf_two(int length);
 123
 124 PRIVATE void        scale_int(const char *s,
 125                               const char *sl,
 126                               double *sc_int);
 127
 128 PRIVATE void        encode_seq( const char *s1,
 129                                 const char *s2);
 130
 131 PRIVATE constrain   *get_ptypes(char *S,
 132                                 const char *structure);
 133
 134 PRIVATE void        get_up_arrays(unsigned int length);
 135
 136 PRIVATE void        free_up_arrays(void);
 137
 138 PRIVATE void        set_encoded_seq(const char *sequence,
 139                                     short **S,
 140                                     short **S1);
 141
 142 PRIVATE void        get_interact_arrays(unsigned int n1,
 143                                         unsigned int n2,
 144                                         pu_contrib *p_c,
 145                                         pu_contrib *p_c2,
 146                                         int w,
 147                                         int incr5,
 148                                         int incr3,
 149                                         double ***p_c_S,
 150                                         double ***p_c2_S);
 151
 152 /*
 153 #################################
 154 # BEGIN OF FUNCTION DEFINITIONS #
 155 #################################
 156 */
 157
 158 PUBLIC pu_contrib *get_pu_contrib_struct(unsigned int n, unsigned int w){
 159   unsigned int i;
 160   pu_contrib  *pu = (pu_contrib *)space(sizeof(pu_contrib));
 161   pu->length      = n;
 162   pu->w           = w;
 163   /* contributions to probability of being unpaired witihin a(n)
 164    H hairpin,
 165    I interior loop,
 166    M muliloop,
 167    E exterior loop*/
 168   /* pu_test->X[i][j] where i <= j and i [1...n], j = [1...w[ */
 169   pu->H           = (double **)space(sizeof(double *) * (n + 1));
 170   pu->I           = (double **)space(sizeof(double *) * (n + 1));
 171   pu->M           = (double **)space(sizeof(double *) * (n + 1));
 172   pu->E           = (double **)space(sizeof(double *) * (n + 1));
 173   for(i=0;i<=n;i++){
 174     pu->H[i]  = (double *)space(sizeof(double) * (w + 1));
 175     pu->I[i]  = (double *)space(sizeof(double) * (w + 1));
 176     pu->M[i]  = (double *)space(sizeof(double) * (w + 1));
 177     pu->E[i]  = (double *)space(sizeof(double) * (w + 1));
 178   }
 179   return pu;
 180 }
 181
 182 PUBLIC  void  free_pu_contrib_struct(pu_contrib *pu){
 183   unsigned int i;
 184   if(pu != NULL){
 185     for(i=0;i<=pu->length;i++){
 186       free(pu->H[i]);
 187       free(pu->I[i]);
 188       free(pu->M[i]);
 189       free(pu->E[i]);
 190     }
 191     free(pu->H);
 192     free(pu->I);
 193     free(pu->M);
 194     free(pu->E);
 195     free(pu);
 196   }
 197 }
 198
 199 /* you have to call pf_fold(sequence, structure); befor pf_unstru */
 200 PUBLIC pu_contrib *pf_unstru(char *sequence, int w){
 201   int           n, i, j, v, k, l, o, p, ij, kl, po, u, u1, d, type, type_2, tt, uu;
 202   unsigned int  size;
 203   double        temp, tqm2;
 204   double        qbt1, *tmp, Zuij, sum_l, *sum_M;
 205   double        *store_H, *store_Io, **store_I2o; /* hairp., interior contribs */
 206   double        *store_M_qm_o,*store_M_mlbase;    /* multiloop contributions */
 207   pu_contrib    *pu_test;
 208
 209   sum_l           = 0.0;
 210   temp            = 0;
 211   n               = (int) strlen(sequence);
 212   sum_M           = (double *)  space((n+1) * sizeof(double));
 213   pu_test         = get_pu_contrib_struct((unsigned)n, (unsigned)w);
 214   size            = ((n+1)*(n+2))>>1;
 215
 216   get_up_arrays((unsigned) n);
 217   init_pf_two(n);
 218
 219   /* init everything */
 220   for (d=0; d<=TURN; d++)
 221     for (i=1; i<=n-d; i++){
 222       j=i+d;
 223       ij = my_iindx[i]-j;
 224       if(d < w) {
 225         pu_test->H[i][d]=pu_test->I[i][d]=pu_test->M[i][d]=pu_test->E[i][d]=0.;
 226       }
 227     }
 228
 229
 230   for (i=0; i<size; i++)
 231     prpr[i]= probs[i];
 232
 233   sum_M[0] = 0.;
 234   for (i=1; i<=n; i++){
 235     /* set auxillary arrays to 0, reuse qqm and qqm1, reuse qqm2 and qq_1m2*/
 236     sum_M[i] = qqm[i] = qqm1[i] = qqm2[i] = qq_1m2[i] = 0;
 237     for (j=i+TURN+1; j<=n; j++){
 238       ij = my_iindx[i]-j;
 239       /* i need the part_func of all structures outside bp[ij] */
 240       if(qb[ij] > 0.0) prpr[ij]= (probs[ij]/qb[ij]);
 241     }
 242   }
 243
 244   /* alloc even more memory */
 245   store_I2o = (double **)space(sizeof(double *) * (n + 1)); /* for p,k */
 246   for(i=0;i<=n;i++)
 247     store_I2o[i] = (double *)space(sizeof(double) * (MAXLOOP + 2));
 248
 249   /* expMLbase[i-p]*dangles_po */
 250   store_M_mlbase = (double *)space(sizeof(double) * (size + 1));
 251
 252   /* 2. exterior bp (p,o) encloses unpaired region [i,i+w[*/
 253   for (o=TURN+2;o<=n; o++) {
 254     double sum_h;
 255     /*allocate space for arrays to store different contributions to H, I & M */
 256     store_H       = (double *)space(sizeof(double) * (o+2));
 257     /* unpaired between ]l,o[ */
 258     store_Io      = (double *)space(sizeof(double) * (o+2));
 259     /* qm[p+1,i-1]*dangles_po */
 260     store_M_qm_o  = (double *)space(sizeof(double) * (n+1));
 261
 262     for (p=o-TURN-1; p>=1; p--) {
 263       /* construction of partition function of segment [p,o], given that
 264          an unpaired region [i,i+w[ exists within [p,o] */
 265       u = o-p-1;
 266       po = my_iindx[p]-o;
 267       type = ptype[po];
 268       if(type){
 269
 270         /*hairpin contribution*/
 271         if (((type==3)||(type==4))&&no_closingGU)
 272           temp = 0.;
 273         else
 274           temp = prpr[po] * exp_E_Hairpin(u, type, S1[p+1], S1[o-1], sequence+p-1, Pf) * scale[u+2];
 275         /* all H contribs are collect for the longest unpaired region */
 276         store_H[p+1] = temp;
 277
 278         /* interior loops with interior pair k,l and an unpaired region of
 279          length w between p and k || l and o*/
 280         for (k=p+1; k<=MIN2(p+MAXLOOP+1,o-TURN-2); k++) {
 281           u1    = k-p-1;
 282           sum_l = 0.;
 283           for (l=MAX2(k+TURN+1,o-1-MAXLOOP+u1); l<o; l++) {
 284             kl      = my_iindx[k]-l;
 285             type_2  = ptype[kl];
 286             if((l+1) < o) store_Io[l+1] += sum_l;
 287
 288             temp=0.;
 289             if (type_2){
 290               type_2 = rtype[type_2];
 291               temp = prpr[po] * qb[kl] * exp_E_IntLoop(u1, o-l-1, type, type_2, S1[p+1], S1[o-1], S1[k-1], S1[l+1], Pf) *scale[u1+o-l+1];
 292               if((l+1) < o) store_Io[l+1] += temp; /* unpaired region between ]l,o[ */
 293               sum_l += temp;
 294             } /* end of if pair(k,l) */
 295           } /* end of l */
 296           /* unpaired in region ]p,k[  */
 297           for(i=p+1;i <= k-1;i++)
 298             store_I2o[i][MIN2(w-1,k-i-1)] += sum_l;
 299         } /* end of k */
 300       } /*end of if(type) test for bp (p,o) */
 301
 302       /* multiple stem loop contribution
 303          calculate qm2[my_iindx[i]-j] in the course of the calculation
 304          of the multiple stem loop contribution:
 305          advantage: you save memory:
 306          instead of a (n+1)*n array for qqm2 you only need 2*n arrays
 307          disadvantage: you have to use two times the op-loop for the full
 308          multiloop contribution
 309          first op-loop: index o goes from 1...n and
 310                         index p from o-TURN-1 ... 1
 311          second op-loop: index o goes from n...1 and
 312                          index p from o+TURN+1 ... n !!
 313          HERE index o goes from 1...n and index p o-TURN-1 ... 1 ,
 314          we calculate the contributions to multiple stem loop
 315          where exp(i+w-1-p)*(qqm2 values between i+w and o-1)
 316          AND qm[iindex[p+1]-(i-1)]*exp(beta*w)*qm[iindex[i+w]-(o-1)]
 317          you have to recalculate of qqm matrix containing final stem
 318          contributions to multiple loop partition function
 319          from segment p,o */
 320
 321       /* recalculate qqm[]
 322          qqm[p] := (contribution with exact one loop in region (p,o)*/
 323       qqm[p]  = qqm1[p] * expMLbase[1];
 324       if(type){
 325         qbt1    =   qb[po] * exp_E_MLstem(type, (p>1) ? S1[p-1] : -1, (o<n) ? S1[o+1] : -1, Pf);
 326         qqm[p]  +=  qbt1;
 327         /* reverse dangles for prpr[po]*... */
 328         temp    =   0.;
 329         tt      =   rtype[type];
 330         temp    =   prpr[po] * exp_E_MLstem(tt, S1[o-1], S1[p+1], Pf) * scale[2] * Pf->expMLclosing;
 331         for(i=p+1; i < o; i++) {
 332           int p1i = (p+1) < (i-1)  ? my_iindx[p+1]-(i-1)  : 0;
 333           /*unpaired region expMLbase[i-p] left of structured
 334             region qq_1m2[i+1]*/
 335           /* @expMLbase:  note distance of i-p == i-(p+1)+1 */
 336           store_M_mlbase[my_iindx[p+1]-i] += expMLbase[i-p] * temp * qq_1m2[i+1];
 337           /* structured region qm[p1i] left of unpaired region */
 338           /* contribition for unpaired region is added after the p-loop */
 339           store_M_qm_o[i] += qm[p1i] * temp;
 340         } /*end of for i ... */
 341       }
 342
 343       for(tqm2 = 0., i=p+1; i < o; i++)
 344         tqm2  +=  qm[my_iindx[p]-i] * qqm[i+1];
 345
 346       /* qqm2[p] contrib with at least 2 loops in region (p,o) */
 347       qqm2[p] = tqm2;
 348     } /* end for (p=..) */
 349
 350     for(sum_h = 0., i=1; i < o; i++) {
 351       int max_v, vo;
 352       sum_h +=  store_H[i];
 353       max_v =   MIN2(w-1,o-i-1);
 354       for(v=max_v; v >= 0; v--){
 355         /* Hairpins */
 356         pu_test->H[i][v] += sum_h;/* store_H[i][v] + store_H[i][max_v]; */
 357         /* Interior loops: unpaired region between  ]l,o[ calculated here !*/
 358         /* unpaired region between ]p,k[ collected after after o-loop */
 359         if(v <= MIN2(max_v,MAXLOOP)) {
 360           pu_test->I[i][v] += store_Io[i]; /* ]l,o[ */
 361         }
 362         /* Multiloops:*/
 363         /* unpaired region [i,v] between structured regions ]p,i[ and ]v,o[. */
 364         /* store_M_qm_o[i] = part. funct over all structured regions ]p,i[ */
 365         vo = (i+v+1) <= (o-1) ? my_iindx[i+v+1]-(o-1): 0;
 366         pu_test->M[i][v] += store_M_qm_o[i]*expMLbase[v+1]*qm[vo];
 367       }
 368     }
 369     tmp = qqm1; qqm1=qqm; qqm=tmp;
 370     tmp = qqm2; qqm2=qq_1m2; qq_1m2=tmp;
 371
 372     free(store_Io);
 373     free(store_H);
 374     free(store_M_qm_o);
 375   }/* end for (o=..) */
 376
 377   for(i=1; i < n; i++) {
 378     int     max_v;
 379     double  sum_iv;
 380     sum_iv  = 0.;
 381     max_v   = MIN2(w-1,n-i);
 382     for(v=n; v >=0; v--) {
 383       if(v <= MIN2(max_v,MAXLOOP)) {
 384         /* all unpaired regions [i,v] between p and k in interior loops */
 385         /* notice v runs from max_v -> 0, sum_iv sums all int. l. contribs */
 386         /* for each x, v < x =< max_v, since they contribute to [i,v] */
 387         sum_iv            += store_I2o[i][v];
 388         pu_test->I[i][v]  += sum_iv;
 389       }
 390       /* all unpaired region [i,v] for a fixed v, given that */
 391       /* region ]v,o[ contains at least 2 structures qq_1m2[v+1]; */
 392       if(v >= i) {
 393         sum_M[v] += store_M_mlbase[my_iindx[i]-v];
 394         if(v-i<=max_v) {
 395           pu_test->M[i][v-i] += sum_M[v];
 396         }
 397       }
 398     }
 399   }
 400
 401   for(i=0;i<=n;i++) {
 402     free(store_I2o[i]);
 403   }
 404   free(store_I2o);
 405
 406   for (i=1; i<=n; i++) {
 407     /* set auxillary arrays to 0 */
 408     qqm[i] = qqm1[i] = qqm2[i] = qq_1m2[i] = 0;
 409   }
 410
 411   /* 2. exterior bp (p,o) encloses unpaired region [i,j]
 412      HERE index o goes from n...1 and index p from o+TURN+1 ... n,
 413      that is, we add the one multiloop contribution that we
 414      could not calculate before  */
 415
 416 /* is free'ing plus allocating faster than looping over all entries an setting them to 0? */
 417 #if 0
 418   free(store_M_mlbase);
 419   store_M_mlbase = (double *) space(sizeof(double) * (size + 1));
 420 #else
 421   /* this should be the fastest way to set everything to 0 */
 422   memset(store_M_mlbase, 0, sizeof(double) * (size + 1));
 423 #endif
 424
 425   for (o=n-TURN-1;o>=1; o--) {
 426     for (p=o+TURN+1; p<=n; p++) {
 427       po    = my_iindx[o]-p;
 428       type  = ptype[po];
 429       /* recalculate of qqm matrix containing final stem
 430          contributions to multiple loop partition function
 431          from segment [o,p] */
 432       qqm[p] = qqm1[p] * expMLbase[1];
 433       if (type) {
 434         qbt1 = qb[po];
 435         qbt1 *= exp_E_MLstem(type, (o>1) ? S1[o-1] : -1, (p<n) ? S1[p+1] : -1, Pf);
 436         qqm[p] += qbt1;
 437         /* revers dangles for prpr[po]...  */
 438         temp=0.;
 439         tt=rtype[type];
 440         temp = prpr[po]*exp_E_MLstem(tt, S1[p-1], S1[o+1], Pf) * Pf->expMLclosing * scale[2];
 441       }
 442       tqm2=0.;
 443       for(i=o+1; i < p; i++) {
 444         uu=0;
 445         tqm2+=qqm[i]*qm[my_iindx[i+1]-p];
 446
 447         if(type !=0) {
 448           double temp2;
 449           temp2=0;
 450           /* structured region qq_1m2[i-1] left of unpaired r. expMLbase[p-i]*/
 451           /* @expMLbase:  note distance of p-i == p+1-i+1 */
 452            store_M_mlbase[my_iindx[i]-p+1] +=  qq_1m2[i-1]*expMLbase[p-i]*temp;
 453         }
 454       }/*end of for i ....*/
 455       qqm2[p] = tqm2;
 456     }/* end for (p=..) */
 457     tmp = qqm1; qqm1=qqm; qqm=tmp;
 458     tmp = qqm2; qqm2=qq_1m2; qq_1m2=tmp;
 459   }/* end for (o=..) */
 460   /* now collect the missing multiloop contributions */
 461   for(i=0;i<=n;i++) { sum_M[i]=0.; }
 462   for(i=1; i<=n;i++) {
 463     int v_max = MIN2(w-1,n-i);
 464     for(v=n; v>=i; v--){
 465       sum_M[i]  += store_M_mlbase[my_iindx[i]-v];
 466       if ((v-i <= v_max) ) {
 467         pu_test->M[i][v-i] += sum_M[i];
 468       }
 469     }
 470   }
 471
 472   /* 1. region [i,j] exterior to all loops */
 473   Zuij=0.;
 474   for (i=1; i<=n; i++) {
 475     uu=0;
 476     for(j=i; j<MIN2(i+w,n+1);j++){
 477       ij=my_iindx[i]-j;
 478       temp=q1k[i-1]*1*scale[j-i+1]*qln[j+1]/q1k[n];
 479       pu_test->E[i][j-i]+=temp;
 480
 481     }
 482   }
 483
 484   free(sum_M);
 485   free(store_M_mlbase);
 486   free_up_arrays();
 487   return pu_test;
 488 }
 489
 490
 491 PRIVATE void  get_interact_arrays(unsigned int n1,
 492                                   unsigned int n2,
 493                                   pu_contrib *p_c,
 494                                   pu_contrib *p_c2,
 495                                   int w,
 496                                   int incr5,
 497                                   int incr3,
 498                                   double ***p_c_S,
 499                                   double ***p_c2_S){
 500
 501   unsigned int i;
 502   int pc_size, j;
 503   *p_c_S = (double **)space(sizeof(double *)*(n1+1));
 504
 505   for (i=1; i<=n1; i++){
 506     pc_size = MIN2((w + incr5 + incr3), (int)n1);
 507     (*p_c_S)[i] = (double *)space(sizeof(double) * (pc_size + 1));
 508     for (j=0; j < pc_size; j++)
 509       (*p_c_S)[i][j] = p_c->H[i][j] + p_c->I[i][j] + p_c->M[i][j] + p_c->E[i][j];
 510   }
 511
 512   if(p_c2 != NULL){
 513     (*p_c2_S) = (double **)space(sizeof(double *) * (n2 + 1));
 514     for (i=1; i<=n2; i++){
 515       pc_size = MIN2(w, (int)n2);
 516       (*p_c2_S)[i]  = (double *)space(sizeof(double) * (pc_size + 2));
 517       for (j=0; j < pc_size; j++)
 518         (*p_c2_S)[i][j] = p_c2->H[i][j] + p_c2->I[i][j] + p_c2->M[i][j] + p_c2->E[i][j];
 519     }
 520   }
 521 }
 522
 523 /*------------------------------------------------------------------------*/
 524 /* s1 is the longer seq */
 525 PUBLIC interact *pf_interact( const char *s1,
 526                               const char *s2,
 527                               pu_contrib *p_c,
 528                               pu_contrib *p_c2,
 529                               int w,
 530                               char *cstruc,
 531                               int incr3,
 532                               int incr5){
 533
 534   int         i, j, k,l,n1,n2,add_i5,add_i3,i_max,k_max, pc_size;
 535   double      temp, Z, rev_d, E, Z2,**p_c_S, **p_c2_S, int_scale;
 536   FLT_OR_DBL  ****qint_4, **qint_ik;
 537   /* PRIVATE double **pint; array for pf_up() output */
 538   interact    *Int;
 539   double      G_min, G_is,Gi_min;
 540   int         gi,gj,gk,gl,ci,cj,ck,cl,prev_k,prev_l;
 541   FLT_OR_DBL  **int_ik;
 542   double      Z_int, temp_int, temppfs;
 543   double      const_scale, const_T;
 544   constrain   *cc = NULL;  /* constrains for cofolding */
 545   char        *Seq, *i_long,*i_short,*pos=NULL; /* short seq appended to long one */
 546   /* int ***pu_jl; */ /* positions of interaction in the short RNA */
 547
 548   G_min = G_is = Gi_min = 100.0;
 549   gi = gj = gk = gl = ci = cj = ck = cl = 0;
 550
 551   n1      = (int) strlen(s1);
 552   n2      = (int) strlen(s2);
 553   prev_k  = 1;
 554   prev_l  = n2;
 555
 556   i_long  = (char *) space (sizeof(char)*(n1+1));
 557   i_short = (char *) space (sizeof(char)*(n2+1));
 558   Seq     = (char *) space (sizeof(char)*(n1+n2+2));
 559
 560   strcpy(Seq,s1);
 561   strcat(Seq,s2);
 562
 563   set_encoded_seq(s1, &S, &S1);
 564   set_encoded_seq(s2, &SS, &SS2);
 565
 566   cc = get_ptypes(Seq,cstruc);
 567
 568   get_interact_arrays(n1, n2, p_c, p_c2, w, incr5, incr3, &p_c_S, &p_c2_S);
 569
 570   /*array for pf_up() output */
 571   Int = (interact *) space(sizeof(interact)*1);
 572   Int->Pi = (double *) space(sizeof(double)*(n1+2));
 573   Int->Gi = (double *) space(sizeof(double)*(n1+2));
 574
 575   /* use a different scaling for pf_interact*/
 576   scale_int(s2, s1, &int_scale);
 577
 578   /* set the global scale array and the global variable pf_scale to the
 579      values used to scale the interaction, keep their former values !! */
 580   temppfs = pf_scale;
 581   pf_scale = int_scale;
 582
 583   /* in order to scale expLoopEnergy correctly call*/
 584   scale_stru_pf_params((unsigned) n1);
 585
 586   qint_ik = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *) * (n1+1));
 587   for (i=1; i<=n1; i++) {
 588     qint_ik[i] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (n1+1));
 589   }
 590 /* int_ik */
 591   int_ik = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL *) * (n1+1));
 592   for (i=1; i<=n1; i++) {
 593     int_ik[i] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (n1+1));
 594   }
 595   Z_int=0.;
 596   /*  Gint = ( -log(int_ik[gk][gi])-( ((int) w/2)*log(pf_scale)) )*((Pf->temperature+K0)*GASCONST/1000.0); */
 597   const_scale = ((int) w/2)*log(pf_scale);
 598   const_T = (Pf->kT/1000.0);
 599   encode_seq(s1, s2);
 600   /* static  short *S~S1, *S1~SS1, *SS~S2, *SS2; */
 601   for (i=0; i<=n1; i++) {
 602     Int->Pi[i]=Int->Gi[i]=0.;
 603   }
 604   E=0.;
 605   Z=0.;
 606
 607   if ( fold_constrained && cstruc != NULL) {
 608     pos = strchr(cstruc,'|');
 609     if(pos) {
 610       ci=ck=cl=cj=0;
 611       /* long seq              & short seq
 612          .........||..|||||....&....||||...  w = maximal interaction length
 613                  ck       ci       cj  cl    */
 614       strncpy(i_long,cstruc,n1);
 615       i_long[n1] = '\0';
 616       strncpy(i_short,&cstruc[n1],n2);
 617       i_short[n2] ='\0';
 618       pos = strchr(i_long,'|');
 619       if(pos) ck = (int) (pos-i_long)+1; /* k */
 620       pos = strrchr(i_long,'|');
 621       if(pos) ci = (int) (pos-i_long)+1; /* i */
 622       pos = strrchr(i_short,'|');
 623       if(pos) cl = (int) (pos-i_short)+1; /* l */
 624       pos = strchr(i_short,'|');
 625       if(pos) cj = (int) (pos-i_short)+1; /* j */
 626
 627       if(ck > 0 && ci > 0 && ci-ck+1 > w) {
 628         fprintf(stderr, "distance between constrains in longer seq, %d, larger than -w = %d",ci-ck+1,w);
 629         nrerror("pf_interact: could not satisfy all constraints");
 630       }
 631       if(cj > 0 && cl > 0 && cl-cj+1 > w) {
 632         fprintf(stderr, "distance between constrains in shorter seq, %d, larger than -w = %d",cl-cj+1,w);
 633         nrerror("pf_interact: could not satisfy all constraints");
 634       }
 635     }
 636
 637   } else if ( fold_constrained && cstruc == NULL) {
 638     nrerror("option -C selected, but no constrained structure given\n");
 639   }
 640   if(fold_constrained) pos = strchr(cstruc,'|');
 641
 642   /*  qint_4[i][j][k][l] contribution that region (k-i) in seq1 (l=n1)
 643       is paired to region (l-j) in seq 2(l=n2) that is
 644       a region closed by bp k-l  and bp i-j */
 645   qint_4 = (FLT_OR_DBL ****) space(sizeof(FLT_OR_DBL ***) * (n1+1));
 646
 647   /* qint_4[i][j][k][l] */
 648   for (i=1; i<=n1; i++) {
 649     int end_k;
 650     end_k = i-w;
 651     if(fold_constrained && pos && ci) end_k= MAX2(i-w, ci-w);
 652     /* '|' constrains for long sequence: index i from 1 to n1 (5' to 3')*/
 653     /* interaction has to include 3' most '|' constrain, ci */
 654     if(fold_constrained && pos && ci && i==1 && i<ci)
 655       i= ci-w+1 > 1 ? ci-w+1 : 1;
 656     /* interaction has to include 5' most '|' constrain, ck*/
 657     if(fold_constrained && pos && ck && i > ck+w-1) break;
 658
 659     /* note: qint_4[i] will be freed before we allocate qint_4[i+1] */
 660     qint_4[i] = (FLT_OR_DBL ***) space(sizeof(FLT_OR_DBL **) * (n2+1));
 661     for (j=n2; j>0; j--) {
 662       qint_4[i][j] = (FLT_OR_DBL **) space(sizeof(FLT_OR_DBL*) * (w+1));
 663       for (k=0; k<=w; k++) {
 664         qint_4[i][j][k] = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL) * (w+1));
 665       }
 666     }
 667
 668      prev_k=1;
 669     for (j=n2; j>0; j--) {
 670       int type, type2,end_l;
 671       end_l = j+w;
 672       if(fold_constrained && pos && ci) end_l= MIN2(cj+w,j+w);
 673       /* '|' constrains for short sequence: index j from n2 to 1 (3' to 5')*/
 674       /* interaction has to include 5' most '|' constrain, cj */
 675       if(fold_constrained && pos && cj && j==n2 && j>cj)
 676         j = cj+w-1 > n2 ? n2 : cj+w-1;
 677       /* interaction has to include 3' most '|' constrain, cl*/
 678       if(fold_constrained && pos && cl && j < cl-w+1) break;
 679       type = cc->ptype[cc->indx[i]-(n1+j)];
 680       qint_4[i][j][0][0] = type ? Pf->expDuplexInit : 0;
 681
 682       if (!type) continue;
 683       qint_4[i][j][0][0] *= exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j<n2) ? SS2[j+1] : -1, Pf);
 684
 685       rev_d = exp_E_ExtLoop(rtype[type], (j>1) ? SS2[j-1] : -1, (i<n1) ? S1[i+1] : -1, Pf);
 686
 687       /* add inc5 and incr3 */
 688       if((i-incr5) > 0 ) add_i5=i-incr5;
 689       else add_i5=1;
 690       add_i3=incr3;
 691       pc_size = MIN2((w+incr3+incr5),n1);
 692       if(incr3 < pc_size) add_i3=incr3;
 693       else add_i3=pc_size-1;
 694
 695       /* only one bp (no interior loop) */
 696       if(p_c2 == NULL) {/* consider only structure of longer seq. */
 697         qint_ik[i][i]+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*scale[((int) w/2)];
 698         Z+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*scale[((int) w/2)];
 699       } else {/* consider structures of both seqs. */
 700         qint_ik[i][i]+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*p_c2_S[j][0]*scale[((int) w/2)];
 701         Z+=qint_4[i][j][0][0]*rev_d*p_c_S[add_i5][add_i3]*p_c2_S[j][0]*scale[((int) w/2)];
 702       }
 703
 704 /* int_ik */
 705       /* check deltaG_ges = deltaG_int + deltaG_unstr; */
 706       int_ik[i][i]+=qint_4[i][j][0][0]*rev_d*scale[((int) w/2)];
 707       Z_int+=qint_4[i][j][0][0]*rev_d*scale[((int) w/2)];
 708       temp_int=0.;
 709
 710       temp=0.;
 711       prev_l = n2;
 712       for (k=i-1; k>end_k && k>0; k--) {
 713         if (fold_constrained && pos && cstruc[k-1] == '|' && k > prev_k)
 714           prev_k=k;
 715         for (l=j+1; l< end_l && l<=n2; l++) {
 716           int a,b,ia,ib,isw;
 717           double scalew, tt, intt;
 718
 719           type2 = cc->ptype[cc->indx[k]-(n1+l)];
 720           /* '|' : l HAS TO be paired: not pair (k,x) where x>l allowed */
 721           if(fold_constrained && pos && cstruc[n1+l-1] == '|' && l < prev_l)
 722             prev_l=l; /*break*/
 723           if(fold_constrained && pos && (k<=ck || i>=ci) && !type2) continue;
 724           if(fold_constrained && pos && ((cstruc[k-1] == '|') || (cstruc[n1+l-1] == '|')) && !type2) break;
 725
 726           if (!type2) continue;
 727           /* to save memory keep only qint_4[i-w...i][][][] in memory
 728              use indices qint_4[i][j][a={0,1,...,w-1}][b={0,1,...,w-1}] */
 729           a=i-k;/* k -> a from 1...w-1*/
 730           b=l-j;/* l -> b from 1...w-1 */
 731
 732           /* scale everything to w/2 */
 733           isw = ((int) w/2);
 734           if ((a+b) < isw ){
 735             scalew = ( scale[isw - (a+b)] );
 736           } else if ( (a+b) > isw ) {
 737             scalew = 1/( scale[(a+b) - isw] );
 738           } else {
 739             scalew = 1;
 740           }
 741
 742           if (i-k+l-j-2<=MAXLOOP) {
 743             if(k >= prev_k && l <= prev_l) { /* don't violate constrains */
 744               E = exp_E_IntLoop(i-k-1,l-j-1, type2, rtype[type],
 745                                 S1[k+1], SS2[l-1], S1[i-1], SS2[j+1], Pf) *
 746                                 scale[i-k+l-j]; /* add *scale[u1+u2+2] */
 747
 748               qint_4[i][j][a][b] += ( qint_4[k][l][0][0]*E);
 749
 750               /* use ia and ib to go from a....w-1 and from b....w-1  */
 751               ia=ib=1;
 752               while((a+ia)<w && i-(a+ia)>=1 && (b+ib)<w && (j+b+ib)<=n2) {
 753                 int iaa,ibb;
 754
 755                 qint_4[i][j][a+ia][b+ib] += qint_4[k][l][ia][ib]*E;
 756
 757                 iaa=ia+1;
 758                 while(a+iaa<w && i-(a+iaa)>=1) {
 759                   qint_4[i][j][a+iaa][b+ib] += qint_4[k][l][iaa][ib]*E;
 760                   ++iaa;
 761                 }
 762
 763                 ibb=ib+1;
 764                 while( (b+ibb)<w && (j+b+ibb)<=n2 ) {
 765                   qint_4[i][j][a+ia][b+ibb] += qint_4[k][l][ia][ibb]*E;
 766                   ++ibb;
 767                 }
 768                 ++ia;
 769                 ++ib;
 770               }
 771             }
 772           }
 773           /* '|' constrain in long sequence */
 774           /* collect interactions starting before 5' most '|' constrain */
 775           if ( fold_constrained && pos && ci && i < ci) continue;
 776           /* collect interactions ending after 3' most '|' constrain*/
 777           if ( fold_constrained && pos && ck &&  k > ck) continue;
 778           /* '|' constrain in short sequence */
 779           /* collect interactions starting before 5' most '|' constrain */
 780           if ( fold_constrained && pos && cj && j > cj) continue;
 781           /* collect interactions ending after 3' most '|' constrain*/
 782           if ( fold_constrained && pos && cl && l < cl) continue;
 783
 784           /* scale everything to w/2*/
 785           /* qint_ik[k][i] all interactions where k and i both are paired */
 786           /* substract incr5 from k */
 787           if(k-incr5 > 0) add_i5=k-incr5;
 788           else add_i5=1;
 789           /* add incr3 to i */
 790           pc_size = MIN2((w+incr3+incr5),n1);
 791           if(i-k+incr3 < pc_size) add_i3=i-k+incr3;
 792           else add_i3=pc_size-1;
 793
 794           if(p_c2 == NULL) {/* consider only structure of longer seq. */
 795             tt = qint_4[i][j][a][b]*p_c_S[add_i5][add_i3]*scalew*rev_d;
 796           } else { /* consider structures of both seqs. */
 797             tt = qint_4[i][j][a][b]*p_c_S[add_i5][add_i3]*p_c2_S[j][b]*scalew*rev_d;
 798           }
 799           temp+= tt;
 800           qint_ik[k][i]+= tt;
 801           /* int_ik */
 802           /* check deltaG_ges = deltaG_int + deltaG_unstr; */
 803           intt = qint_4[i][j][a][b]*scalew*rev_d;
 804           temp_int += intt;
 805           int_ik[k][i]+= intt;
 806           G_is = (-log(tt)-const_scale)*(const_T);
 807           if (G_is < G_min || EQUAL(G_is,G_min)) {
 808             G_min = G_is;
 809             Gi_min =(-log(intt)-const_scale)*(const_T);
 810             gi=i;
 811             gj=j;
 812             gk=k;
 813             gl=l;
 814           }
 815         }
 816       }
 817       Z+=temp;
 818       /* int_ik */
 819       Z_int+=temp_int;
 820     }
 821
 822     /* free qint_4 values not needed any more */
 823     if(i > w) {
 824       int bla;
 825       bla=i-w;
 826       if (fold_constrained && pos && ci && i-w < ci-w+1) continue;
 827       if (fold_constrained && pos && ci) bla = MAX2(ci-w+1,i-w);
 828       for (j=n2; j>0; j--) {
 829         for (k=0; k<=w; k++){
 830           free(qint_4[bla][j][k]);
 831         }
 832         free(qint_4[bla][j]);
 833       }
 834       free(qint_4[bla]);
 835       qint_4[bla] = NULL;
 836     }
 837   }
 838
 839
 840   Z2=0.0;
 841   i_max = k_max = 0;
 842   for (i=1; i<=n1; i++) {
 843     for (k=i; k<=n1 && k<i+w; k++) {
 844       Z2+=qint_ik[i][k];
 845       for(l=i;l<=k;l++) {
 846         /* Int->Pi[l]: prob that position l is within a paired region */
 847         /* qint_ik[i][k] as well as Z are scaled to scale[((int) w/2) */
 848         Int->Pi[l]+=qint_ik[i][k]/Z;
 849         /* Int->Gi[l]: minimal delta G at position [l] */
 850         Int->Gi[l]=MIN2(Int->Gi[l],
 851                        ( -log(qint_ik[i][k])-( ((int) w/2)*log(pf_scale)) )*
 852                        (Pf->kT/1000.0) );
 853       }
 854     }
 855   }
 856   if(n1 > w){
 857     int start_i,end_i;
 858     start_i = n1-w+1;
 859     end_i=n1;
 860     if (fold_constrained && pos && ci) {
 861       /* a break in the k loop might result in unfreed values */
 862       start_i = ci-w+1 < n1-w+1 ? ci-w+1 : n1-w+1;
 863       start_i = start_i > 0 ? start_i : 1;
 864       /* start_i = ck; */
 865       end_i = ck+w-1 > n1 ? n1 : ck+w-1;
 866     }
 867     for (i=start_i; i<=end_i; i++) {
 868       if(qint_4[i] == NULL ) continue;
 869       for (j=n2; j>0; j--) {
 870         for (k=0; k<=w; k++) {
 871           free(qint_4[i][j][k]);
 872         }
 873         free(qint_4[i][j]);
 874       }
 875       free(qint_4[i]);
 876     }
 877     free(qint_4);
 878   } else {
 879     int start_i,end_i;
 880     start_i = 1;
 881     end_i=n1;
 882     if (fold_constrained && pos) {
 883       start_i = ci-w+1 > 0 ? ci-w+1 : 1;
 884       end_i = ck+w-1 > n1 ? n1 : ck+w-1;
 885     }
 886
 887     for (i=start_i; i<=end_i; i++) {
 888       for (j=n2; j>0; j--) {
 889         for (k=0; k<=w; k++) {
 890           free(qint_4[i][j][k]);
 891         }
 892         free(qint_4[i][j]);
 893       }
 894       free(qint_4[i]);
 895     }
 896     free(qint_4);
 897   }
 898   if(fold_constrained && (gi==0 || gk==0 ||  gl==0 || gj==0)) {
 899     nrerror("pf_interact: could not satisfy all constraints");
 900   }
 901   /* fill structure interact */
 902   Int->length = n1;
 903   Int->i = gi;
 904   Int->j = gj;
 905   Int->k = gk;
 906   Int->l = gl;
 907   Int->Gikjl = G_min;
 908   Int->Gikjl_wo = Gi_min;
 909
 910   free(i_long);
 911   free(i_short);
 912
 913   for (i=1; i<=n1; i++) {
 914     free(int_ik[i]);
 915   }
 916   free(int_ik);
 917   for (i=1; i<=n1; i++) {
 918     free(qint_ik[i]);
 919   }
 920   free(qint_ik);
 921
 922   /* reset the global variables pf_scale and scale to their original values */
 923   pf_scale = temppfs;/* reset pf_scale */
 924   scale_stru_pf_params((unsigned) n1);/* reset the scale array */
 925   free_pf_arrays(); /* for arrays for pf_fold(...) */
 926
 927   if(expMLbase != NULL) {
 928     free(expMLbase);
 929     expMLbase = NULL;
 930   }
 931   if(scale != NULL) {
 932     free(scale);
 933     scale = NULL;
 934   }
 935   for (i=1; i<=n1; i++) {
 936     free(p_c_S[i]);
 937   }
 938   free(p_c_S);
 939   if(p_c2 != NULL) {
 940     for (i=1; i<=n2; i++) {
 941       free(p_c2_S[i]);
 942     }
 943     free(p_c2_S);
 944   }
 945   free(Seq);
 946   free(cc->indx);
 947   free(cc->ptype);
 948   free(cc);
 949   return(Int);
 950 }
 951 /*------------------------------------------------------------------------*/
 952 /* use an extra scale for pf_interact, here sl is the longer sequence */
 953 PRIVATE void scale_int(const char *s, const char *sl, double *sc_int){
 954   int       n,nl,l_scales;
 955   duplexT   mfe;
 956   double    kT;
 957
 958   n         = strlen(s);
 959   nl        = strlen(sl);
 960
 961   expMLbase = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*(nl+1));
 962   scale     = (FLT_OR_DBL *) space(sizeof(FLT_OR_DBL)*((nl+1)*2));
 963
 964   /* use RNA duplex to get a realistic estimate for the best possible
 965      interaction energy between the short RNA s and its target sl */
 966   mfe = duplexfold(s,sl);
 967
 968   kT = Pf->kT/1000.0;   /* in Kcal */
 969
 970   /* sc_int is similar to pf_scale: i.e. one time the scale */
 971   *sc_int = exp(-(mfe.energy)/kT/n);
 972
 973   /* free the structure returned by duplexfold */
 974   free(mfe.structure);
 975 }
 976
 977 /*----------------------------------------------------------------------*/
 978 /* init_pf_two(n) :gets the arrays, that you need, from part_func.c */
 979 /* get_pf_arrays(&S, &S1, &ptype, &qb, &qm, &q1k, &qln);*/
 980 /* init_pf_fold(), update_pf_params, encode_char(), make_ptypes() are called by pf_fold() */
 981 PRIVATE void init_pf_two(int length){
 982 #ifdef SUN4
 983   nonstandard_arithmetic();
 984 #else
 985 #ifdef HP9
 986   fpsetfastmode(1);
 987 #endif
 988 #endif
 989   make_pair_matrix();
 990
 991   /* gets the arrays, that we need, from part_func.c */
 992   if(!get_pf_arrays(&S, &S1, &ptype, &qb, &qm, &q1k, &qln))
 993     nrerror("init_pf_two: pf_fold() has to be called before calling pf_unstru()\n");
 994   /* get a pointer to the base pair probs */
 995   probs = export_bppm();
 996
 997   scale_stru_pf_params((unsigned) length);
 998
 999   init_length=length;
1000   if(init_temp != Pf->temperature)
1001     nrerror("init_pf_two: inconsistency with temperature");
1002 }
1003
1004 PRIVATE void  get_up_arrays(unsigned int length){
1005   unsigned int l1 = length + 1;
1006   unsigned int l2 = length + 2;
1007   prpr      = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL)  * ((l1*l2)>>1));
1008   expMLbase = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL)  * l1);
1009   scale     = (FLT_OR_DBL *)space(sizeof(FLT_OR_DBL)  * l1);
1010   qqm2      = (double *)    space(sizeof(double)      * l2);
1011   qq_1m2    = (double *)    space(sizeof(double)      * l2);
1012   qqm       = (double *)    space(sizeof(double)      * l2);
1013   qqm1      = (double *)    space(sizeof(double)      * l2);
1014   my_iindx  = get_iindx(length);
1015 }
1016
1017 PRIVATE void  free_up_arrays(void){
1018   if(prpr       != NULL){ free(prpr);       prpr      = NULL;}
1019   if(expMLbase  != NULL){ free(expMLbase);  expMLbase = NULL;}
1020   if(scale      != NULL){ free(scale);      scale     = NULL;}
1021   if(qqm        != NULL){ free(qqm);        qqm       = NULL;}
1022   if(qqm1       != NULL){ free(qqm1);       qqm1      = NULL;}
1023   if(qqm2       != NULL){ free(qqm2);       qqm2      = NULL;}
1024   if(qq_1m2     != NULL){ free(qq_1m2);     qq_1m2    = NULL;}
1025   if(my_iindx   != NULL){ free(my_iindx);   my_iindx  = NULL;}
1026 }
1027
1028 PUBLIC void free_interact(interact *pin) {
1029   if(S != NULL && pin != NULL){
1030     free(S);
1031     S=NULL;
1032   }
1033   if(S1 != NULL && pin != NULL){
1034     free(S1);
1035     S1=NULL;
1036   }
1037   if(pin != NULL){
1038     free(pin->Pi);
1039     free(pin->Gi);
1040     free(pin);
1041     pin=NULL;
1042   }
1043 }
1044 /*---------------------------------------------------------------------------*/
1045
1046 PRIVATE void encode_seq(const char *s1, const char *s2) {
1047   unsigned int i,l;
1048
1049   l = strlen(s1);
1050   /* S and S1 are freed by free_pf_arrays(); ! */
1051   S = (short *) space(sizeof(short)*(l+1));
1052   S1= (short *) space(sizeof(short)*(l+1));
1053   /* S1 exists only for the special X K and I bases and energy_set!=0 */
1054   S[0] = l;
1055   for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
1056     S[i]= (short) encode_char(toupper(s1[i-1]));
1057     S1[i] = alias[S[i]];   /* for mismatches of nostandard bases */
1058   }
1059   if(s2 != NULL) {
1060     l = strlen(s2);
1061     /* SS2 exists only for the special X K and I bases and energy_set!=0 */
1062     SS[0] = l;
1063     for (i=1; i<=l; i++) { /* make numerical encoding of sequence */
1064       SS[i]= (short) encode_char(toupper(s2[i-1]));
1065       SS2[i] = alias[SS[i]];   /* for mismatches of nostandard bases */
1066     }
1067   }
1068 }
1069
1070 /*-------------------------------------------------------------------------*/
1071  /* scale energy parameters and pre-calculate Boltzmann weights:
1072   most of this is done in structure Pf see params.c,h (function:
1073   get_scaled_pf_parameters(), only arrays scale and expMLbase are handled here*/
1074 PRIVATE void scale_stru_pf_params(unsigned int length)
1075 {
1076   unsigned int i;
1077   double  kT;
1078
1079
1080   /* Do this only at the first call for get_scaled_pf_parameters()
1081      and/or if temperature has changed*/
1082   if(init_temp != temperature) {
1083     if(Pf) free(Pf);
1084     Pf=get_scaled_pf_parameters();
1085   }
1086
1087   init_temp = Pf->temperature;
1088
1089   kT = Pf->kT;   /* kT in cal/mol  */
1090
1091    /* scaling factors (to avoid overflows) */
1092   if (pf_scale == -1) { /* mean energy for random sequences: 184.3*length cal */
1093     pf_scale = exp(-(-185+(Pf->temperature-37.)*7.27)/kT);
1094     if (pf_scale<1) pf_scale=1;
1095   }
1096   scale[0] = 1.;
1097   scale[1] = 1./pf_scale;
1098   expMLbase[0] = 1;
1099   expMLbase[1] = Pf->expMLbase/pf_scale;
1100   for (i=2; i<=length; i++) {
1101     scale[i] = scale[i/2]*scale[i-(i/2)];
1102     expMLbase[i] = pow(Pf->expMLbase, (double)i) * scale[i];
1103   }
1104 }
1105 /*-------------------------------------------------------------------------*/
1106 /* make a results structure containing all u-values & the header */
1107 PUBLIC pu_out *get_u_vals(pu_contrib *p_c, int **unpaired_values, char *select_contrib) {
1108   int i, j, k, l, g,ws,num_u_vals,unstr,count,contribs,size,w,len;
1109   int S,E,H,I,M;
1110   int off_S, off_E, off_H, off_I, off_M;
1111   /* double **p_cont,**p_cont_sh, dG_u; p_u AND its contributions */
1112   pu_out* u_results;
1113
1114   len = p_c->length;
1115
1116   /* number of different -u values */
1117   for (num_u_vals = 0, i = 1; i <= unpaired_values[0][0]; i++) {
1118     j = unpaired_values[i][0];
1119     do num_u_vals++; while(++j <= unpaired_values[i][1]);
1120   }
1121   /* check which contributions ([-c "SHIME"] ) are desired by the user,
1122      set the offset for each contribution */
1123   contribs = 0;
1124   S = E = H = I = M = 0;
1125   off_S = off_E = off_H = off_I = off_M = 0;
1126   if(strchr(select_contrib, 'S')) {
1127     S=1;
1128     off_S = contribs;
1129     ++contribs;
1130   }
1131   if(strchr(select_contrib, 'E')) {
1132     E=1;
1133     off_E = contribs;
1134     ++contribs;
1135   }
1136   if(strchr(select_contrib, 'H')) {
1137     H=1;
1138     off_H = contribs;
1139     ++contribs;
1140   }
1141   if(strchr(select_contrib, 'I')) {
1142     I=1;
1143     off_I = contribs;
1144     ++contribs;
1145   }
1146   if(strchr(select_contrib, 'M')) {
1147     M=1;
1148     off_M = contribs;
1149     ++contribs;
1150   }
1151
1152   if(contribs > 5) {
1153     nrerror("get_u_vals: error with contribs!");
1154   }
1155   /* allocate the results structure */
1156   u_results = (pu_out *) space(1*sizeof(pu_out));
1157   u_results->len = len; /* sequence length */
1158   /*num_u_vals differnet -u values, contribs [-c "SHIME"] */
1159   u_results->u_vals = num_u_vals;
1160   u_results->contribs = contribs;
1161   /* add 1 column for position within the sequence and
1162      add 1 column for the free energy of interaction values */
1163   /* header e.g. u3I (contribution for u3 interior loops */
1164   size = 1 + (num_u_vals*contribs) + 1;
1165   u_results->header = (char **) space((size+1)*sizeof(char*));
1166   for(i=0;i<(size+1);i++){
1167     u_results->header[i] = (char *) space(10*sizeof(char));
1168   }
1169   /* different free energies for all  -u and -c combinations */
1170   u_results->u_values = (double**) space((size+1) *sizeof(double*));
1171   for(i=0;i<(size+1);i++){
1172     /* position within the sequence  */
1173     u_results->u_values[i] = (double*) space((len+3)*sizeof(double));
1174   }
1175   /* write the position within the sequence in the u_results array
1176      at column zerro */
1177   sprintf(u_results->header[0],"pos");
1178   for(i=0;i<=len;i++){
1179     /* add the position*/
1180     u_results->u_values[0][i] = i;
1181   }
1182   /* go over the different -u values, u_vals[] listy of different -u values*/
1183   for (count = k = 1; k <= unpaired_values[0][0]; k++) {
1184     l = unpaired_values[k][0];
1185     do{
1186       int offset; /* offset for the respective -u value (depents on the number
1187                    of the -u value and on the numbers of contribs */
1188
1189       offset = ((count - 1) * contribs) + 1; /* first colum is the position */
1190       /* set the current value of -u : here we call it w */
1191       w = l; /* set w to the actual -u value */
1192       if(w > len) break; /* corr caro */
1193       /* make the header - look which contribitions are wanted */
1194       if(S) sprintf(u_results->header[offset+off_S],"u%dS",w);
1195       if(E) sprintf(u_results->header[offset+off_E],"u%dE",w);
1196       if(H) sprintf(u_results->header[offset+off_H],"u%dH",w);
1197       if(I) sprintf(u_results->header[offset+off_I],"u%dI",w);
1198       if(M) sprintf(u_results->header[offset+off_M],"u%dM",w);
1199
1200       if(p_c != NULL) {
1201         for (i=1; i<=len; i++) { /* for each position */
1202           /* w goes form j to i (intervall end at i) */
1203           for (j=i; j < MIN2((i+w),len+1); j++) { /* for each -u value < w
1204                                                 this is not necessay ->
1205                                                 calculate j from i and w
1206                                                 : (j-i+1) == w */
1207             double blubb;
1208             /* if (j-i+1) == w we have the -u = w value wanted */
1209             if( (j-i+1) == w && i+w-1 <= len) {
1210               blubb = p_c->H[i][j-i]+p_c->I[i][j-i]+p_c->M[i][j-i]+p_c->E[i][j-i];
1211
1212               /* printf("len %d  blubb %.3f \n",len, blubb); */
1213               if(S) u_results->u_values[offset+off_S][i+w-1]+=blubb;
1214               if(E) u_results->u_values[offset+off_E][i+w-1]+=p_c->E[i][j-i];
1215               if(H) u_results->u_values[offset+off_H][i+w-1]+=p_c->H[i][j-i];
1216               if(I) u_results->u_values[offset+off_I][i+w-1]+=p_c->I[i][j-i];
1217               if(M) u_results->u_values[offset+off_M][i+w-1]+=p_c->M[i][j-i];
1218
1219             }
1220             if(i<w && (j-i+1) != w && i+w-1 > len &&  i+w-1 < len+3) {
1221               if(S) u_results->u_values[offset+off_S][i+w-1]=-1;
1222               if(E) u_results->u_values[offset+off_E][i+w-1]=-1;
1223               if(H) u_results->u_values[offset+off_H][i+w-1]=-1;
1224               if(I) u_results->u_values[offset+off_I][i+w-1]=-1;
1225               if(M) u_results->u_values[offset+off_M][i+w-1]=-1;
1226             }
1227           }
1228         }
1229       } else return(NULL); /* error */
1230       count++;
1231     } while(++l <= unpaired_values[k][1]);
1232   }
1233   return(u_results); /*success*/
1234 }
1235 /* plot the results structure */
1236 /* when plotting the results for the target seq we add a header */
1237 /* when plotting the results for the interaction partner u want no header,
1238    set s1 to NULL to avoid plotting the header */
1239 /* currently we plot the free energies to a file: the probability of
1240    being unpaired for region [i,j], p_u[i,j], is related to the free
1241    energy to open region [i,j], dG_u[i,j] by:
1242    dG_u[i,j] = -log(p_u[i,j])*(temperature+K0)*GASCONST/1000.0; */
1243 PUBLIC int plot_free_pu_out(pu_out* res, interact *pint, char *ofile, char *head) {
1244   int size,s,i,len;
1245   double dG_u;
1246   char nan[4], *time, dg[11];
1247   FILE *wastl;
1248   double  kT = Pf->kT;
1249   wastl = fopen(ofile,"a");
1250   if (wastl==NULL) {
1251     fprintf(stderr, "p_cont: can't open %s for Up_plot\n", ofile);
1252     return(0);
1253   }
1254   sprintf(dg,"dG");
1255
1256   /* printf("T=%.16f \n(temperature+K0)*GASCONST/1000.0 = %.16f\n",temperature,(temperature+K0)*GASCONST/1000.0); */
1257
1258   /* write the header of the output file:  */
1259   /*  # timestamp commandlineaufruf   */
1260   /*  # length and name of first sequence (target)
1261   /*  # first seq */
1262   /*  # length and name of second sequence (interaction partner) */
1263   /*  # second seq */
1264   /* the next line is the output for the target: colums
1265      position in target | dG_unpaired values for target | interaction energy */
1266   /*  # pos   u1S   u1H  dg */
1267   /*  values for target */
1268   /* if -b was choosen: the next lines are the dG_unpaired values for
1269      the interaction partner */
1270   /*  # pos   u1S   u1H  */
1271   /*  values for the interaction partner */
1272
1273   /* print header, if nh is zerro */
1274   if(head){
1275     time = time_stamp();
1276     fprintf(wastl,"# %s\n", time);
1277     fprintf(wastl,"%s\n",head);
1278   }
1279   fprintf(wastl,"# ");
1280   /* }  else { fprintf(wastl," "); } close if before  */
1281   len  = res->len;
1282   size = res->u_vals * res->contribs;
1283
1284   sprintf(nan,"NA");
1285   nan[2] = '\0';
1286
1287   for(i=0;i<=len; i++) {
1288     for(s=0;s<=size+1;s++) { /* that is for different contribution */
1289       if ( i== 0 && s > size && pint != NULL)
1290         fprintf(wastl,"%8s  ",dg);
1291       if(i != 0) {
1292         if(s>0 && s<=size) {
1293           if(res->u_values[s][i] > 0.0) {
1294             dG_u = -log(res->u_values[s][i])*kT/1000.0;
1295             fprintf(wastl,"%8.3f  ",dG_u);
1296           } else { /* no p_u value was defined print nan*/
1297             fprintf(wastl,"%8s  ",nan);
1298           }
1299
1300         } else if (s > size && pint != NULL) {
1301           fprintf(wastl,"%8.3f  ",pint->Gi[i]);
1302         } else if (s == 0) {
1303           fprintf(wastl,"%8.0f  ",res->u_values[s][i]);
1304         }
1305       } else {
1306         if(s>1) {
1307           fprintf(wastl,"%8s  ",res->header[s]);
1308         } else {
1309           fprintf(wastl,"%7s  ",res->header[s]);
1310         }
1311       }
1312     }
1313     fprintf(wastl,"\n");
1314   }
1315   fclose(wastl);
1316   /*free pu_out* res */
1317   if(res != NULL) {
1318     for(i=0;i<=(size+2);i++) {
1319       free(res->u_values[i]);
1320       free(res->header[i]);
1321     }
1322     free(res->u_values);
1323     free(res->header);
1324     free(res);
1325     res = NULL;
1326   }
1327
1328   return(1); /* success */
1329 }
1330
1331 PUBLIC int Up_plot(pu_contrib *p_c, pu_contrib *p_c_sh, interact *pint, char *ofile, int **unpaired_values, char *select_contrib, char *head, unsigned int mode) {
1332   pu_out *dada;
1333   int ret;
1334   /* check what case we have */
1335
1336   /* upmode = 1 only one seq */
1337   /* if(p_c != NULL && pint == NULL) { */
1338   if(mode & RNA_UP_MODE_1){
1339     dada = get_u_vals(p_c,unpaired_values,select_contrib);
1340     ret = plot_free_pu_out(dada,NULL,ofile,head);
1341
1342   /* upmode > 1 cofolding */
1343   /* } else if (p_c != NULL && pint != NULL) { */
1344   } else if(mode & RNA_UP_MODE_2) {
1345     dada = get_u_vals(p_c,unpaired_values,select_contrib);
1346     ret = plot_free_pu_out(dada,pint,ofile,head);
1347
1348   /* upmode = 3  cofolding*/
1349   /* } else if (p_c == NULL && p_c_sh != NULL) { */
1350   }
1351   if(mode & RNA_UP_MODE_3) {
1352     dada  = get_u_vals(p_c,unpaired_values, select_contrib);
1353     ret   = plot_free_pu_out(dada, pint, ofile, head);
1354
1355     dada = get_u_vals(p_c_sh, unpaired_values, select_contrib);
1356     ret = plot_free_pu_out(dada,NULL,ofile, NULL);
1357   }
1358   return(ret);
1359 }
1360
1361 /*-------------------------------------------------------------------------*/
1362 /* copy from part_func_co.c */
1363 PRIVATE constrain *get_ptypes(char *Seq, const char *structure) {
1364   int n,i,j,k,l, length;
1365   constrain *con;
1366   short *s, *s1;
1367
1368   length = strlen(Seq);
1369   make_pair_matrix();
1370   con = (constrain *) space(sizeof(constrain));
1371   con->indx = (int *) space(sizeof(int)*(length+1));
1372   for (i=1; i<=length; i++) {
1373     con->indx[i] = ((length+1-i)*(length-i))/2 +length+1;
1374   }
1375   con->ptype = (char *) space(sizeof(char)*((length+1)*(length+2)/2));
1376
1377   set_encoded_seq((const char *)Seq, &s, &s1);
1378
1379   n=s[0];
1380   for (k=1; k<=n-CO_TURN-1; k++)
1381     for (l=1; l<=2; l++) {
1382       int type,ntype=0,otype=0;
1383       i=k; j = i+CO_TURN+l; if (j>n) continue;
1384       type = pair[s[i]][s[j]];
1385       while ((i>=1)&&(j<=n)) {
1386         if ((i>1)&&(j<n)) ntype = pair[s[i-1]][s[j+1]];
1387         if (noLonelyPairs && (!otype) && (!ntype))
1388           type = 0; /* i.j can only form isolated pairs */
1389         con->ptype[con->indx[i]-j] = (char) type;
1390         otype =  type;
1391         type  = ntype;
1392         i--; j++;
1393       }
1394     }
1395
1396   if (fold_constrained&&(structure!=NULL)) {
1397     int hx, *stack;
1398     char type;
1399     stack = (int *) space(sizeof(int)*(n+1));
1400     for(hx=0, j=1; j<=n; j++) {
1401       switch (structure[j-1]) {
1402       case 'x': /* can't pair */
1403         for (l=1; l<j-CO_TURN; l++) con->ptype[con->indx[l]-j] = 0;
1404         for (l=j+CO_TURN+1; l<=n; l++) con->ptype[con->indx[j]-l] = 0;
1405         break;
1406       case '(':
1407         stack[hx++]=j;
1408         /* fallthrough */
1409       case '<': /* pairs upstream */
1410         break;
1411       case ')':
1412         if (hx<=0) {
1413           fprintf(stderr, "%s\n", structure);
1414           nrerror("1. unbalanced brackets in constraints");
1415         }
1416         i = stack[--hx];
1417         type = con->ptype[con->indx[i]-j];
1418         /* don't allow pairs i<k<j<l */
1419         for (k=i; k<=j; k++)
1420           for (l=j; l<=n; l++) con->ptype[con->indx[k]-l] = 0;
1421         /* don't allow pairs k<i<l<j */
1422         for (k=1; k<=i; k++)
1423           for (l=i; l<=j; l++) con->ptype[con->indx[k]-l] = 0;
1424         con->ptype[con->indx[i]-j] = (type==0)?7:type;
1425       case '>': /* pairs downstream */
1426         break;
1427       }
1428     }
1429     if (hx!=0) {
1430       fprintf(stderr, "%s\n", structure);
1431       nrerror("2. unbalanced brackets in constraint string");
1432     }
1433     free(stack);
1434   }
1435   free(s);
1436   free(s1);
1437   return con;
1438 }
1439 PRIVATE  void  set_encoded_seq(const char *sequence, short **S, short **S1){
1440   unsigned int i,l;
1441   l = strlen(sequence);
1442   if(S!= NULL){
1443     *S  = (short *)space(sizeof(short) * (l + 2));
1444     for(i=1; i<=l; i++) /* make numerical encoding of sequence */
1445       (*S)[i]= (short) encode_char(toupper(sequence[i-1]));
1446     (*S)[l+1] = (*S)[1];
1447     (*S)[0]   = (short) l;
1448   }
1449   /* S1 exists only for the special X K and I bases and energy_set!=0 */
1450   if(S1 != NULL){
1451     *S1 = (short *)space(sizeof(short) * (l + 2));
1452     for(i=1; i<=l; i++) /* make numerical encoding of sequence */
1453       (*S1)[i]  = alias[(short) encode_char(toupper(sequence[i-1]))]; /* for mismatches of nostandard bases */
1454     /* for circular folding add first base at position n+1 and last base at position 0 in S1 */
1455     (*S1)[l+1]  = (*S1)[1];
1456     (*S1)[0]    = (*S1)[l];
1457   }
1458 }