1 /* Last changed Time-stamp: <2009-02-18 14:19:51 ivo> */
3 local pair probabilities for RNA secondary structures
5 Stephan Bernhart, Ivo L Hofacker
9 todo: compute energy z-score for each window
17 #include <float.h> /* #defines FLT_MAX ... */
19 #include "energy_par.h"
20 #include "fold_vars.h"
23 #include "part_func.h"
25 #include "loop_energies.h"
35 PRIVATE char rcsid[] UNUSED = "$Id: LPfold.c,v 1.8 2009/02/18 20:34:38 ivo Exp $";
37 #define ISOLATED 256.0
40 #################################
42 #################################
46 #################################
48 #################################
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;
66 PRIVATE double alpha = 1.0;
70 /* NOTE: all variables are assumed to be uninitialized if they are declared as threadprivate
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)
80 #################################
81 # PRIVATE FUNCTION DECLARATIONS #
82 #################################
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);*/
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);
104 #################################
105 # BEGIN OF FUNCTION DEFINITIONS #
106 #################################
109 PRIVATE void init_partfunc_L(int length, pf_paramT *parameters){
110 if (length<1) nrerror("init_partfunc_L: length must be greater 0");
112 /* Explicitly turn off dynamic threads */
114 free_pf_arrays_L(); /* free previous allocation */
116 if (init_length>0) free_pf_arrays_L(); /* free previous allocation */
120 nonstandard_arithmetic();
127 get_arrays_L((unsigned) length);
128 scale_pf_params((unsigned) length, parameters);
131 init_length = length;
135 PRIVATE void get_arrays_L(unsigned int length){
136 /*arrays in 2 dimensions*/
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));
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 *));
162 my_iindx = get_iindx(length);
163 iindx = get_iindx(length); /* for backward compatibility and Perl wrapper */
164 jindx = get_indx(length);
167 PRIVATE void free_pf_arrays_L(void){
179 if(prm_l) free(prm_l);
180 if(prm_l1) free(prm_l1);
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);
191 if(pf_params) free(pf_params);
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;
201 standard_arithmetic();
213 PUBLIC void update_pf_paramsLP(int length){
214 update_pf_paramsLP_par(length, NULL);
217 PUBLIC void update_pf_paramsLP_par(int length, pf_paramT *parameters){
219 init_partfunc_L(length, parameters);
221 if(parameters) init_partfunc_L(length, parameters);
222 else if (length > init_length) init_partfunc_L(length, parameters);
224 /* make_pair_matrix();*/
225 scale_pf_params((unsigned) length, parameters);
230 PUBLIC plist *pfl_fold( char *sequence,
238 return pfl_fold_par(sequence, winSize, pairSize, cutoffb, pU, dpp2, pUfp, spup, NULL);
241 PUBLIC plist *pfl_fold_par( char *sequence,
249 pf_paramT *parameters){
251 int n, m, i, j, k, l, u, u1, ii, type, type_2, tt, ov, do_dpp, simply_putout, noGUclosure;
253 FLT_OR_DBL temp, Qmax, prm_MLb, prmt, prmt1, qbt1, *tmp, expMLclosing;
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");
275 if(dpp !=NULL) do_dpp=1;
277 n = (int) strlen(sequence);
278 if (n<TURN+2) return 0;
281 /* always init everything since all global static variables are uninitialized when entering a thread */
282 init_partfunc_L(n, parameters);
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);
289 expMLclosing = pf_params->expMLclosing;
290 noGUclosure = pf_params->model_details.noGUclosure;
293 max_real = (sizeof(FLT_OR_DBL) == sizeof(float)) ? FLT_MAX : DBL_MAX;
295 S = encode_sequence(sequence, 0);
296 S1 = encode_sequence(sequence, 1);
298 /* make_ptypes(S, structure); das machmadochlieber lokal, ey!*/
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*/
305 for (i=1; i<=ulength; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
308 for (i=1; i<=n; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
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 */
315 pl = (struct plist *)space(1000*sizeof(struct plist));
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)];
324 for (j=TURN+2;j<=n+winSize; j++) {
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) */
336 /*hairpin contribution*/
337 if (((type==3)||(type==4))&&noGUclosure) qbt1 = 0;
339 qbt1 = exp_E_Hairpin(u, type, S1[i+1], S1[j-1], sequence+i-1, pf_params) * scale[u+2];
341 /* interior loops with interior pair k,l */
342 for (k=i+1; k<=MIN2(i+MAXLOOP+1,j-TURN-2); k++) {
344 for (l=MAX2(k+TURN+1,j-1-MAXLOOP+u1); l<j; l++) {
345 type_2 = ptype[k][l];
347 type_2 = rtype[type_2];
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];
354 /*multiple stem loop contribution*/
355 ii = my_iindx[i+1]; /* ii-k=[i+1,k-1] */
357 for (k=i+2; k<=j-1; k++) temp += qm[i+1][k-1]*qqm1[k];
359 qbt1 += temp * expMLclosing * exp_E_MLstem(tt, S1[j-1], S1[i+1], pf_params) * scale[2];
362 } /* end if (type!=0) */
365 /* construction of qqm matrix containing final stem
366 contributions to multiple loop partition function
368 qqm[i] = qqm1[i]*expMLbase[1];
370 qbt1 = qb[i][j] * exp_E_MLstem(type, (i>1) ? S1[i-1] : -1, (j<n) ? S1[j+1] : -1, pf_params);
374 /*construction of qm matrix containing multiple loop
375 partition function contributions from segment i,j */
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]);
384 /*auxiliary matrix qq for cubic order q calculation below */
387 qbt1 *= exp_E_ExtLoop(type, (i>1) ? S1[i-1] : -1, (j < n) ? S1[j+1] : -1, pf_params);
389 qq[i] = qq1[i]*scale[1] + qbt1;
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];
398 if (Qmax>max_real/10.)
399 fprintf(stderr, "Q close to overflow: %d %d %g\n", i,j,temp);
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);
408 tmp = qq1; qq1 =qq; qq =tmp;
409 tmp = qqm1; qqm1=qqm; qqm=tmp;
412 /* just as a general service, I save here the free energy of the windows
413 no output is generated, however,...
415 if ((j>=winSize) && (j<=n) && (ulength) && !(pUoutput)) {
417 Fwindow=(-log(q[j-winSize+1][j])-winSize*log(pf_params->pf_scale))*pf_params->kT/1000.0;
421 if (ulength>=winSize)
422 pU[j][winSize]=scale[winSize]/q[j-winSize+1][j];
428 /* initialize multiloopfs */
429 for (k=j-winSize; k<=MIN2(n,j); k++) {
432 /* prm_l1[k]=0; others stay*/
436 for (l=k+TURN+1; l<=MIN2(n,k+winSize-1); l++) {
438 pR[k][l] = 0; /* set zero at start */
440 if (qb[k][l]==0) continue;
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];
446 pR[k][l]+=1./q[k][l];
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];
453 pR[k][l] *= exp_E_ExtLoop(type, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params);
455 type_2 = ptype[k][l];
456 type_2 = rtype[type_2];
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++) {
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];
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++) {
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];
480 /* 3. bonding k,l as substem of multi-loop enclosed by i,m */
482 if(k>1) /*sonst nix!*/
483 for (l=MIN2(n-1,k+winSize-2); l>=k+TURN+1; l--) { /* opposite direction */
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];
493 prmt *= expMLclosing;
495 prm_l[m] = prm_l1[m]*expMLbase[1]+prmt1;
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];
501 prml[m] = prml[ m] + prm_l[m];
503 if (qb[k][l] == 0.) continue;
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;
517 for (m=MIN2(k+winSize-2,n);m>=l+2; m--)
518 temp += prml[m]*qm[l+1][m-1];
520 temp *= exp_E_MLstem(tt, (k>1) ? S1[k-1] : -1, (l<n) ? S1[l+1] : -1, pf_params) * scale[2];
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]);
529 if (pR[k][l]>=max_real) {
534 } /* end for (l=..) */
535 tmp = prm_l1; prm_l1=prm_l; prm_l=tmp;
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 ;)*/
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 *));
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));
557 //dont want to have that yet?
559 for (i=1; i<=n; i++) pU[i]=(double *)space((MAX2(MAXLOOP,ulength)+2)*sizeof(double));
561 compute_pU_splitup(k-MAXLOOP-1,ulength,pU,pUO,pUH, pUI, pUM, winSize, n, sequence);
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');
570 compute_pU(k-MAXLOOP-1,ulength,pU, winSize, n, sequence);
572 /* here, we put out and free pUs not in use any more (hopefully) */
574 putoutpU(pU,k-MAXLOOP-1, ulength, pUfp);
578 if (j-(2*winSize+MAXLOOP+1)>0) {
579 printpbar(pR,winSize,j-(2*winSize+MAXLOOP+1),n);
581 print_plist(n, j-(2*winSize+MAXLOOP+1), pR, winSize, spup);
584 pl=get_plistW(pl, n, j-(2*winSize+MAXLOOP+1), pR, winSize);
586 if (do_dpp)dpp=get_deppp(dpp,j-(2*winSize-MAXLOOP),pairSize, n);
587 FreeOldArrays(j-(2*winSize+MAXLOOP+1));
589 } /* end if (do_backtrack)*/
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);
600 for (j=MAX2(n-winSize-MAXLOOP,1); j<=n; j++) {
601 printpbar(pR,winSize,j,n);
603 print_plist(n, j, pR, winSize, spup);
606 pl=get_plistW(pl, n, j, pR, winSize);
608 if ((do_dpp)&&j<n) dpp=get_deppp(dpp,j,pairSize, n);
611 /* free_pf_arrays_L(); */
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);
623 PRIVATE void scale_pf_params(unsigned int length, pf_paramT *parameters){
625 double kT, scaling_factor;
627 if(pf_params) free(pf_params);
630 pf_params = get_boltzmann_factor_copy(parameters);
633 set_model_details(&md);
634 pf_params = get_boltzmann_factors(temperature, alpha, md, pf_scale);
637 scaling_factor = pf_params->pf_scale;
638 kT = pf_params->kT; /* kT in cal/mol */
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;
647 scale[1] = 1./scaling_factor;
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];
656 PRIVATE void printpbar(FLT_OR_DBL **prb,int winSize, int i, int n) {
658 int howoften=0; /* how many samples do we have for this pair */
661 for (j=i+TURN; j<MIN2(i+winSize,n+1); j++) {
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;
672 PRIVATE void FreeOldArrays(int i) {
673 /*free arrays no longer needed*/
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));
692 q[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
694 qb[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
696 qm[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
699 qm2[j]=(FLT_OR_DBL *)space((winSize+1)*sizeof(FLT_OR_DBL));
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));
705 ptype[j]=(char *)space((winSize+1)*sizeof(char));
711 PRIVATE void GetPtype(int i, int winSize,const short *S,int n) {
712 /*make new entries in ptype array*/
715 for (j=i; j<=MIN2(i+winSize,n); j++) {
716 type = pair[S[i]][S[j]];
717 ptype[i][j] = (char) type;
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 */
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) {
735 pl=(plist *)xrealloc(pl,max_p*sizeof(plist));
739 pl[num_p++].p=Tpr[start][j];
742 /* mark end of data with zeroes */
746 /* pl=(struct plist *)xrealloc(pl,(count)*sizeof(struct plist)); */
751 PRIVATE plist *get_deppp(plist *pl, int start, int pairsize, int length) {
752 /* compute dependent pair probabilities */
756 temp=(plist *)space(pairsize*sizeof(plist)); /* holds temporary deppp */
757 for (j=start+TURN; j<MIN2(start+pairsize,length); j++) {
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];
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++) {
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 */
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]);
795 /* mark end of data with zeroes */
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
810 FLT_OR_DBL expMLclosing = pf_params->expMLclosing;
812 QBE=(double *) space((MAX2(ulength,MAXLOOP)+2)*sizeof(double));
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++) {
821 printf("i%d j%d ",i5,j3);
824 if (ptype[i5][j3]!=0) {/**/
825 /* (.. >-----|..........)
828 temp = (i5<k) ? qm2[i5+1][k] * expMLbase[j3-k-1] : 0.; /* (..{}{}-----|......) */
831 temp += qm2[k+ulength+1][j3-1] * expMLbase[k+ulength-i5]; /* (..|-----|{}{}) */
833 if((i5<k)&&(j3-1>k+ulength))
834 temp += qm[i5+1][k] * qm[k+ulength+1][j3-1] * expMLbase[ulength]; /* ({}|-----|{}) */
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;
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*/
842 pU[k+ulength][ulength] += temp;
847 /* code doubling to avoid if within loop */
849 /*initialization for interior loops,
850 it is not recomended to have verysmall ulengths!!*/
851 if (ulength<MAXLOOP) {
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:
861 i5<k5<k-TURN k5<=i5+l3+2+MAXLOOP-j3
864 for (i5=MAX2(k+ulength-winSize,1);i5<k-TURN-1;i5++) {
866 for (j3=k+ulength+1; j3<=MIN2(n,MIN2(i5+winSize-1,k+MAXLOOP+1)); j3++) {
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++){
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]);
877 pU[k+ulength][ulength]+= temp;
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
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++) {
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++) {
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]);
898 pU[k+ulength][ulength]+= temp;
902 /* Add up Is QI5[l][m-l-1] QI3 */
903 /* Add up Interior loop terms */
906 for (len=winSize; len>=ulength; len--) temp+=QI3[k][len];
907 for (;len>0; len--) {
913 for (len=winSize; len>=MAX2(ulength,MAXLOOP); len--) temp+=QI5[k][len];
914 for (;len>0; len--) {
916 QBE[len] += temp; /* replace QBE with QI */
918 /* Add Hairpinenergy to QBE */
920 for(obp = MIN2(n, k + winSize - 1); obp > k + ulength; 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--){
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) */
928 /* doubling the code to get the if out of the loop */
930 /* Add up Multiloopterms qmb[l][m]+=prml[m]*dang;
931 q2l[l][m]+=(prml[m]-prm_l[m])*dang; */
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 */
940 for(len = 1; len < ulength; len++){
941 for(obp = k + len + TURN; obp <= MIN2(n, k + winSize - 1); obp++){
943 QBE[len] += qmb[k][obp-k-1] * qm[k+len+1/*2*/][obp-1] * expMLbase[len];
946 for (len=1; len<ulength; len++) {
947 for (obp=k+len+TURN+TURN; obp<=MIN2(n,k+winSize-1); 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 (___()()) */
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];
962 if ((ulength>=winSize)&&(k>=ulength)) {
963 pU[k][winSize]=scale[winSize]/q[k-winSize+1][k];
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--){
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];
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];
978 /* Divide by number of possible windows */
979 pU[k][startu] += temp;
981 int leftmost, rightmost;
983 leftmost = MAX2(1, k - winSize + 1);
984 rightmost = MIN2(n - winSize + 1, k - startu + 1);
985 pU[k][startu] /= (rightmost - leftmost + 1);
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!*/
997 fprintf(fp,"%d\t",k);
998 for (i=1; i<=MIN2(ulength,k); i++) {
999 fprintf(fp,"%.5g\t",pUx[k][i]);
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!*/
1008 fprintf(fp,"%d\t",k);
1009 for (i=1; i<=MIN2(ulength,k); i++) {
1010 fprintf(fp,"%.5g\t",pUx[k][i]);
1012 fprintf(fp,"\t%s\n",ident);
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);
1021 PUBLIC void putoutpU_prob_par(double **pU,int length, int ulength, FILE *fp, int energies, pf_paramT *parameters){
1022 /*put out unpaireds */
1024 double kT = parameters->kT/1000.0;
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);
1033 for (k=1; k<=length; k++){
1034 fprintf(fp,"%d\t",k);
1035 for (i=1; i<=ulength; i++) {
1040 if (energies) temp=-log(pU[k][i])*kT;
1042 fprintf(fp,"%.7g\t",temp);
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);
1054 PUBLIC void putoutpU_prob_bin_par(double **pU,int length, int ulength, FILE *fp, int energies, pf_paramT *parameters) {
1056 /*put out unpaireds */
1058 double kT= parameters->kT/1000.0;
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 */
1069 fwrite(p,sizeof(int),1,fp);
1072 for (i=1; i<=ulength; i++) {
1073 for (k=1; k<=11; k++){/* write first ten entries to 1000000 */
1075 fwrite(p,sizeof(int),1,fp);
1077 for (k=1; k<=length; k++){/* write data now */
1079 p[0]=1000000; /* check if u > pos */
1080 fwrite(p,sizeof(int),1,fp);
1084 p[0]= (int) rint(100 *(-log(pU[k][i])*kT));
1085 fwrite(p,sizeof(int),1,fp);
1088 for (k=1; k<=9; k++){/* finish by writing the last 10 entries */
1090 fwrite(p,sizeof(int),1,fp);
1093 /* free pU array; */
1094 for (k=1; k<=length; k++){
1103 Here: Space for questions...
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
1119 FLT_OR_DBL expMLclosing = pf_params->expMLclosing;
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));
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++) {
1133 printf("i%d j%d ",i5,j3);
1136 if (ptype[i5][j3]!=0) {/**/
1137 /* (.. >-----|..........)
1138 i5 j j+ulength j3 */
1140 temp = (i5<k) ? qm2[i5+1][k] * expMLbase[j3-k-1] : 0.; /* (..{}{}-----|......) */
1143 temp += qm2[k+ulength+1][j3-1] * expMLbase[k+ulength-i5]; /* (..|-----|{}{}) */
1145 if((i5<k)&&(j3-1>k+ulength))
1146 temp += qm[i5+1][k] * qm[k+ulength+1][j3-1] * expMLbase[ulength]; /* ({}|-----|{}) */
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;
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*/
1154 pU[k+ulength][ulength] += temp;
1159 /* code doubling to avoid if within loop */
1161 for (len=winSize; len>=MAX2(ulength,MAXLOOP); len--) temp+=QI5[k][len];
1162 for (;len>0; len--) {
1163 temp += QI5[k][len];
1165 QBE[len] += temp; /* replace QBE with QI */
1167 /* Add Hairpinenergy to QBE */
1169 for(obp = MIN2(n, k + winSize - 1); obp > k + ulength; 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--){
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) */
1178 /* doubling the code to get the if out of the loop */
1180 /* Add up Multiloopterms qmb[l][m]+=prml[m]*dang;
1181 q2l[l][m]+=(prml[m]-prm_l[m])*dang; */
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];
1189 QBE[len] += temp; /* add (()()____) type cont. to I3 */
1191 for(len = 1; len < ulength; len++){
1192 for(obp = k + len + TURN; obp <= MIN2(n, k + winSize - 1); obp++){
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];
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 (___()()) */
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];
1219 if ((ulength>=winSize)&&(k>=ulength)) {
1220 pUO[k][winSize]=scale[winSize]/q[k-winSize+1][k];
1223 if ((ulength>=winSize)&&(k>=ulength)) {
1224 pU[k][winSize]=scale[winSize]/q[k-winSize+1][k];
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--){
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];
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];
1239 /* Divide by number of possible windows */
1240 pU[k][startu] += temp;
1241 pUO[k][startu] += temp;
1244 int leftmost, rightmost;
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);
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 */
1265 double kT= (temperature+K0)*GASCONST/1000.0;
1267 if (energies) fprintf(fp,"#opening energies\n #i$\tl=");
1268 else fprintf(fp,"#unpaired probabilities\n #i$\tl=");
1270 fprintf(fp,"Total\n");
1271 for (i=1; i<=ulength; i++) {
1272 fprintf(fp,"%d\t", i);
1276 for (k=1; k<=length; k++){
1277 fprintf(fp,"%d\t",k);
1278 for (i=1; i<=ulength; i++) {
1283 if (energies) temp=-log(pU[k][i])*kT;
1285 fprintf(fp,"%.7g\t",temp);
1287 fprintf(fp,"\tT\n");
1290 fprintf(fp,"\n###################################################################\nHairpin\n");
1291 for (i=1; i<=ulength; i++) {
1292 fprintf(fp,"%d\t", i);
1296 for (k=1; k<=length; k++){
1297 fprintf(fp,"%d\t",k);
1298 for (i=1; i<=ulength; i++) {
1303 if (energies) temp=-log(pUH[k][i])*kT;
1304 else temp=pUH[k][i];
1305 fprintf(fp,"%.7g\t",temp);
1307 fprintf(fp,"\tH\n");
1310 fprintf(fp,"\n###################################################################\nInterior\n");
1311 for (i=1; i<=ulength; i++) {
1312 fprintf(fp,"%d\t", i);
1316 for (k=1; k<=length; k++){
1317 fprintf(fp,"%d\t",k);
1318 for (i=1; i<=ulength; i++) {
1323 if (energies) temp=-log(pUI[k][i])*kT;
1324 else temp=pUI[k][i];
1325 fprintf(fp,"%.7g\t",temp);
1327 fprintf(fp,"\tI\n");
1330 fprintf(fp,"\n###################################################################\nMultiloop\n");
1331 for (i=1; i<=ulength; i++) {
1332 fprintf(fp,"%d\t", i);
1336 for (k=1; k<=length; k++){
1337 fprintf(fp,"%d\t",k);
1338 for (i=1; i<=ulength; i++) {
1343 if (energies) temp=-log(pUM[k][i])*kT;
1344 else temp=pUM[k][i];
1345 fprintf(fp,"%.7g\t",temp);
1347 fprintf(fp,"\tM\n");
1350 fprintf(fp,"\n###################################################################\nExterior\n");
1351 for (i=1; i<=ulength; i++) {
1352 fprintf(fp,"%d\t", i);
1354 fprintf(fp,"\t E\n");
1356 for (k=1; k<=length; k++){
1357 fprintf(fp,"%d\t",k);
1358 for (i=1; i<=ulength; i++) {
1363 if (energies) temp=-log(pUO[k][i])*kT;
1364 else temp=pUO[k][i];
1365 fprintf(fp,"%.7g\t",temp);
1374 /*###########################################*/
1375 /*# deprecated functions below #*/
1376 /*###########################################*/
1378 PUBLIC void init_pf_foldLP(int length){ /* DO NOTHING */}