2 read energy parameters from a file
4 Stephan Kopp, Ivo Hofacker
16 #include "energy_const.h"
17 #include "energy_par.h"
18 #include "read_epars.h"
21 #define PRIVATE static
29 PRIVATE void display_array(int *p, int size, int line, FILE *fp);
30 PRIVATE char *get_array1(int *arr, int size);
31 PRIVATE void ignore_comment(char *line);
32 PRIVATE void check_symmetry(void);
33 PRIVATE void update_nst(int array[NBPAIRS+1][NBPAIRS+1][5][5][5][5]);
36 *** read a 1dimensional array from file
37 *** \param array a pointer to the first element in the array
38 *** \param dim the size of the array
39 *** \param shift the first position the new values will be written in
41 PRIVATE void rd_1dim(int *array, int dim, int shift);
42 PRIVATE void rd_1dim_slice(int *array, int dim, int shift, int post);
43 PRIVATE void rd_2dim(int *array,
45 int shift1, int shift2);
46 PRIVATE void rd_2dim_slice(int *array,
48 int shift1, int shift2,
49 int post1, int post2);
50 PRIVATE void rd_3dim(int *array,
51 int dim1, int dim2, int dim3,
52 int shift1, int shift2, int shift3);
53 PRIVATE void rd_3dim_slice(int *array,
54 int dim1, int dim2, int dim3,
55 int shift1, int shift2, int shift3,
56 int post1, int post2, int post3);
57 PRIVATE void rd_4dim(int *array,
58 int dim1, int dim2, int dim3, int dim4,
59 int shift1, int shift2, int shift3, int shift4);
60 PRIVATE void rd_4dim_slice(int *array,
61 int dim1, int dim2, int dim3, int dim4,
62 int shift1, int shift2, int shift3, int shift4,
63 int post1, int post2, int post3, int post4);
64 PRIVATE void rd_5dim(int *array,
65 int dim1, int dim2, int dim3, int dim4, int dim5,
66 int shift1, int shift2, int shift3, int shift4, int shift5);
67 PRIVATE void rd_5dim_slice(int *array,
68 int dim1, int dim2, int dim3, int dim4, int dim5,
69 int shift1, int shift2, int shift3, int shift4, int shift5,
70 int post1, int post2, int post3, int post4, int post5);
71 PRIVATE void rd_6dim(int *array,
72 int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
73 int shift1, int shift2, int shift3, int shift4, int shift5, int shift6);
74 PRIVATE void rd_6dim_slice(int *array,
75 int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
76 int shift1, int shift2, int shift3, int shift4, int shift5, int shift6,
77 int post1, int post2, int post3, int post4, int post5, int post6);
78 PRIVATE void rd_Tetraloop37(void);
79 PRIVATE void rd_Triloop37(void);
80 PRIVATE void rd_Hexaloop37(void);
82 /*------------------------------------------------------------*/
83 PUBLIC void read_parameter_file(const char fname[]){
84 char *line, ident[256];
88 if (!(fp=fopen(fname,"r"))) {
90 "\nread_parameter_file:\n"
91 "\t\tcan't open file %s\n"
92 "\t\tusing default parameters instead.\n", fname);
96 if (!(line = get_line(fp))) {
97 fprintf(stderr," File %s is inproper.\n", fname);
102 if (strncmp(line,"## RNAfold parameter file v2.0",30)!=0) {
104 "Missing header line in file.\n"
105 "May be this file has not v2.0 format.\n"
106 "Use INTERRUPT-key to stop.\n");
110 while((line=get_line(fp))) {
112 r = sscanf(line, "# %255s", ident);
114 type = gettype(ident);
117 case S: rd_2dim(&(stack37[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
119 case S_H: rd_2dim(&(stackdH[0][0]), NBPAIRS+1, NBPAIRS+1, 1, 1);
121 case HP: rd_1dim(&(hairpin37[0]), 31, 0);
123 case HP_H: rd_1dim(&(hairpindH[0]), 31, 0);
125 case B: rd_1dim(&(bulge37[0]), 31, 0);
127 case B_H: rd_1dim(&(bulgedH[0]), 31, 0);
129 case IL: rd_1dim(&(internal_loop37[0]), 31, 0);
131 case IL_H: rd_1dim(&(internal_loopdH[0]), 31, 0);
133 case MME: rd_3dim(&(mismatchExt37[0][0][0]),
137 case MME_H: rd_3dim(&(mismatchExtdH[0][0][0]),
141 case MMH: rd_3dim(&(mismatchH37[0][0][0]),
145 case MMH_H: rd_3dim(&(mismatchHdH[0][0][0]),
149 case MMI: rd_3dim(&(mismatchI37[0][0][0]),
153 case MMI_H: rd_3dim(&(mismatchIdH[0][0][0]),
157 case MMI1N: rd_3dim(&(mismatch1nI37[0][0][0]),
161 case MMI1N_H: rd_3dim(&(mismatch1nIdH[0][0][0]),
165 case MMI23: rd_3dim(&(mismatch23I37[0][0][0]),
169 case MMI23_H: rd_3dim(&(mismatch23IdH[0][0][0]),
173 case MMM: rd_3dim(&(mismatchM37[0][0][0]),
177 case MMM_H: rd_3dim(&(mismatchMdH[0][0][0]),
181 case INT11: rd_4dim(&(int11_37[0][0][0][0]),
182 NBPAIRS+1, NBPAIRS+1, 5, 5,
185 case INT11_H: rd_4dim(&(int11_dH[0][0][0][0]),
186 NBPAIRS+1, NBPAIRS+1, 5, 5,
189 case INT21: rd_5dim(&(int21_37[0][0][0][0][0]),
190 NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
193 case INT21_H: rd_5dim(&(int21_dH[0][0][0][0][0]),
194 NBPAIRS+1, NBPAIRS+1, 5, 5, 5,
197 case INT22: rd_6dim_slice(&(int22_37[0][0][0][0][0][0]),
198 NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
201 update_nst(int22_37);
203 case INT22_H: rd_6dim_slice(&(int22_dH[0][0][0][0][0][0]),
204 NBPAIRS+1, NBPAIRS+1, 5, 5, 5, 5,
207 update_nst(int22_dH);
209 case D5: rd_2dim(&(dangle5_37[0][0]), NBPAIRS+1, 5, 1, 0);
211 case D5_H: rd_2dim(&(dangle5_dH[0][0]), NBPAIRS+1, 5, 1, 0);
213 case D3: rd_2dim(&(dangle3_37[0][0]), NBPAIRS+1, 5, 1, 0);
215 case D3_H: rd_2dim(&(dangle3_dH[0][0]), NBPAIRS+1, 5, 1, 0);
219 rd_1dim(&values[0], 6, 0);
220 ML_BASE37 = values[0];
221 ML_BASEdH = values[1];
222 ML_closing37 = values[2];
223 ML_closingdH = values[3];
224 ML_intern37 = values[4];
225 ML_interndH = values[5];
230 rd_1dim(&values[0], 3, 0);
233 MAX_NINIO = values[2];
238 rd_1dim(&values[0], 4, 0);
239 DuplexInit37 = values[0];
240 DuplexInitdH = values[1];
241 TerminalAU37 = values[2];
242 TerminalAUdH = values[3];
245 case TL: rd_Tetraloop37();
247 case TRI: rd_Triloop37();
249 case HEX: rd_Hexaloop37();
251 default: /* do nothing but complain */
252 fprintf(stderr,"read_epars: Unknown field identifier in `%s'\n", line);
254 } /* else ignore line */
263 /*------------------------------------------------------------*/
265 PRIVATE void display_array(int *p, int size, int nl, FILE *fp)
269 for (i=1; i<=size; i++, p++) {
272 case INF: fprintf(fp," INF"); break;
273 case -INF: fprintf(fp," -INf"); break;
274 case DEF: fprintf(fp," DEF"); break;
275 default: fprintf(fp,"%6d", *p); break;
277 if ((i%nl)==0) fprintf(fp,"\n");
279 if (size%nl) fprintf(fp,"\n");
283 /*------------------------------------------------------------*/
285 PRIVATE char *get_array1(int *arr, int size)
287 int i, p, pos, pp, r, last;
294 if (!line) nrerror("unexpected end of file in get_array1");
295 ignore_comment(line);
297 while ((i<size)&&(sscanf(line+pos,"%15s%n", buf, &pp)==1)) {
299 if (buf[0]=='*') {i++; continue;}
300 else if (buf[0]=='x') { /* should only be used for loop parameters */
301 if (i==0) nrerror("can't extrapolate first value");
302 p = arr[last] + (int) (0.5+ lxc37*log(((double) i)/(double)(last)));
304 else if (strcmp(buf,"DEF") == 0) p = DEF;
305 else if (strcmp(buf,"INF") == 0) p = INF;
306 else if (strcmp(buf,"NST") == 0) p = NST;
308 r=sscanf(buf,"%d", &p);
311 fprintf(stderr, "can't interpret `%s' in get_array1\n", buf);
324 PRIVATE void rd_1dim(int *array, int dim, int shift){
325 rd_1dim_slice(array, dim, shift, 0);
328 PRIVATE void rd_1dim_slice(int *array, int dim, int shift, int post){
330 cp = get_array1(array+shift, dim-shift-post);
333 fprintf(stderr,"\nrd_1dim: %s\n", cp);
339 PRIVATE void rd_2dim(int *array, int dim1, int dim2, int shift1, int shift2){
340 rd_2dim_slice(array, dim1, dim2, shift1, shift2, 0, 0);
343 PRIVATE void rd_2dim_slice(int *array,
345 int shift1, int shift2,
346 int post1, int post2){
348 int delta_pre = shift1 + shift2;
349 int delta_post = post1 + post2;
351 if(delta_pre + delta_post == 0){
352 rd_1dim(array, dim1 * dim2, 0);
355 for (i=shift1; i<dim1 - post1; i++)
356 rd_1dim_slice(array + (i*dim2), dim2, shift2, post2);
360 PRIVATE void rd_3dim(int *array, int dim1, int dim2, int dim3, int shift1, int shift2, int shift3){
363 shift1, shift2, shift3,
367 PRIVATE void rd_3dim_slice(int *array,
368 int dim1, int dim2, int dim3,
369 int shift1, int shift2, int shift3,
370 int post1, int post2, int post3){
372 int delta_pre = shift1 + shift2 + shift3;
373 int delta_post = post1 + post2 + post3;
375 if(delta_pre + delta_post == 0){
376 rd_1dim(array, dim1 * dim2 * dim3, 0);
379 for (i=shift1; i<dim1 - post1; i++){
380 rd_2dim_slice(array + (i * dim2 * dim3),
388 PRIVATE void rd_4dim(int *array,
389 int dim1, int dim2, int dim3, int dim4,
390 int shift1, int shift2, int shift3, int shift4){
392 dim1, dim2, dim3, dim4,
393 shift1, shift2, shift3, shift4,
397 PRIVATE void rd_4dim_slice(int *array,
398 int dim1, int dim2, int dim3, int dim4,
399 int shift1, int shift2, int shift3, int shift4,
400 int post1, int post2, int post3, int post4){
402 int delta_pre = shift1 + shift2 + shift3 + shift4;
403 int delta_post = post1 + post2 + post3 + post4;
405 if(delta_pre + delta_post == 0){
406 rd_1dim(array, dim1 * dim2 * dim3 * dim4, 0);
409 for(i=shift1; i<dim1 - post1; i++){
410 rd_3dim_slice(array + (i * dim2 * dim3 * dim4),
412 shift2, shift3, shift4,
413 post2, post3, post4);
418 PRIVATE void rd_5dim(int *array,
419 int dim1, int dim2, int dim3, int dim4, int dim5,
420 int shift1, int shift2, int shift3, int shift4, int shift5){
422 dim1, dim2, dim3, dim4, dim5,
423 shift1, shift2, shift3, shift4, shift5,
427 PRIVATE void rd_5dim_slice(int *array,
428 int dim1, int dim2, int dim3, int dim4, int dim5,
429 int shift1, int shift2, int shift3, int shift4, int shift5,
430 int post1, int post2, int post3, int post4, int post5){
432 int delta_pre = shift1 + shift2 + shift3 + shift4 + shift5;
433 int delta_post = post1 + post2 + post3 + post4 + post5;
435 if(delta_pre + delta_post == 0){
436 rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5, 0);
439 for(i=shift1; i<dim1 - post1; i++)
440 rd_4dim_slice(array + (i * dim2 * dim3 * dim4 * dim5),
441 dim2, dim3, dim4, dim5,
442 shift2, shift3, shift4, shift5,
443 post2, post3, post4, post5);
448 *** \param dim1 The size of the first dimension
449 *** \param shift1 The pre shift for the first dimension
451 PRIVATE void rd_6dim(int *array,
452 int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
453 int shift1, int shift2, int shift3, int shift4, int shift5, int shift6){
455 dim1, dim2, dim3, dim4, dim5, dim6,
456 shift1, shift2, shift3, shift4, shift5, shift6,
460 PRIVATE void rd_6dim_slice(int *array,
461 int dim1, int dim2, int dim3, int dim4, int dim5, int dim6,
462 int shift1, int shift2, int shift3, int shift4, int shift5, int shift6,
463 int post1, int post2, int post3, int post4, int post5, int post6){
465 int delta_pre = shift1 + shift2 + shift3 + shift4 + shift5 + shift6;
466 int delta_post = post1 + post2 + post3 + post4 + post5 + post6;
468 if(delta_pre + delta_post == 0){
469 rd_1dim(array, dim1 * dim2 * dim3 * dim4 * dim5 * dim6, 0);
472 for(i=shift1; i<dim1 - post1; i++)
473 rd_5dim_slice(array + (i * dim2 * dim3 * dim4 * dim5 * dim6),
474 dim2, dim3, dim4, dim5, dim6,
475 shift2, shift3, shift4, shift5, shift6,
476 post2, post3, post4, post5, post6);
481 /*------------------------------------------------------------*/
482 PRIVATE void rd_Tetraloop37(void)
488 /* erase old tetraloop entries */
489 memset(&Tetraloops, 0, 281);
490 memset(&Tetraloop37, 0, sizeof(int)*40);
491 memset(&TetraloopdH, 0, sizeof(int)*40);
494 if (buf==NULL) break;
495 r = sscanf(buf,"%6s %d %d", &Tetraloops[7*i], &Tetraloop37[i], &TetraloopdH[i]);
496 strcat(Tetraloops, " ");
499 } while((r==3)&&(i<40));
503 /*------------------------------------------------------------*/
504 PRIVATE void rd_Hexaloop37(void)
510 /* erase old hexaloop entries */
511 memset(&Hexaloops, 0, 361);
512 memset(&Hexaloop37, 0, sizeof(int)*40);
513 memset(&HexaloopdH, 0, sizeof(int)*40);
516 if (buf==NULL) break;
517 r = sscanf(buf,"%8s %d %d", &Hexaloops[9*i], &Hexaloop37[i], &HexaloopdH[i]);
518 strcat(Hexaloops, " ");
521 } while((r==3)&&(i<40));
525 /*------------------------------------------------------------*/
526 PRIVATE void rd_Triloop37(void)
532 /* erase old hexaloop entries */
533 memset(&Triloops, 0, 241);
534 memset(&Triloop37, 0, sizeof(int)*40);
535 memset(&TriloopdH, 0, sizeof(int)*40);
538 if (buf==NULL) break;
539 r = sscanf(buf,"%5s %d %d", &Triloops[6*i], &Triloop37[i], &TriloopdH[i]);
540 strcat(Triloops, " ");
543 } while((r==3)&&(i<40));
547 /*------------------------------------------------------------*/
550 PRIVATE void ignore_comment(char * line)
552 /* excise C style comments */
553 /* only one comment per line, no multiline comments */
556 if ((cp1=strstr(line, "/*"))) {
557 cp2 = strstr(cp1, "*/");
559 nrerror("unclosed comment in parameter file");
560 /* can't use strcpy for overlapping strings */
561 for (cp2+=2; *cp2!='\0'; cp2++, cp1++)
568 /*------------------------------------------------------------*/
570 PUBLIC char *settype(enum parset s){
572 case S: return "stack";
573 case S_H: return "stack_enthalpies";
574 case HP: return "hairpin";
575 case HP_H: return "hairpin_enthalpies";
576 case B: return "bulge";
577 case B_H: return "bulge_enthalpies";
578 case IL: return "interior";
579 case IL_H: return "interior_enthalpies";
580 case MME: return "mismatch_exterior";
581 case MME_H: return "mismatch_exterior_enthalpies";
582 case MMH: return "mismatch_hairpin";
583 case MMH_H: return "mismatch_hairpin_enthalpies";
584 case MMI: return "mismatch_interior";
585 case MMI_H: return "mismatch_interior_enthalpies";
586 case MMI1N: return "mismatch_interior_1n";
587 case MMI1N_H: return "mismatch_interior_1n_enthalpies";
588 case MMI23: return "mismatch_interior_23";
589 case MMI23_H: return "mismatch_interior_23_enthalpies";
590 case MMM: return "mismatch_multi";
591 case MMM_H: return "mismatch_multi_enthalpies";
592 case D5: return "dangle5";
593 case D5_H: return "dangle5_enthalpies";
594 case D3: return "dangle3";
595 case D3_H: return "dangle3_enthalpies";
596 case INT11: return "int11";
597 case INT11_H: return "int11_enthalpies";
598 case INT21: return "int21";
599 case INT21_H: return "int21_enthalpies";
600 case INT22: return "int22";
601 case INT22_H: return "int22_enthalpies";
602 case ML: return "ML_params";
603 case NIN: return "NINIO";
604 case TRI: return "Triloops";
605 case TL: return "Tetraloops";
606 case HEX: return "Hexaloops";
607 case QUIT: return "END";
608 case MISC: return "Misc";
609 default: nrerror("\nThe answer is: 42\n");
613 /*------------------------------------------------------------*/
615 PUBLIC enum parset gettype(const char *ident){
616 if (strcmp(ident,"stack") == 0) return S;
617 else if (strcmp(ident,"stack_enthalpies") == 0) return S_H;
618 else if (strcmp(ident,"hairpin") == 0) return HP;
619 else if (strcmp(ident,"hairpin_enthalpies") == 0) return HP_H;
620 else if (strcmp(ident,"bulge") == 0) return B;
621 else if (strcmp(ident,"bulge_enthalpies") == 0) return B_H;
622 else if (strcmp(ident,"interior") == 0) return IL;
623 else if (strcmp(ident,"interior_enthalpies") == 0) return IL_H;
624 else if (strcmp(ident,"mismatch_exterior") == 0) return MME;
625 else if (strcmp(ident,"mismatch_exterior_enthalpies") == 0) return MME_H;
626 else if (strcmp(ident,"mismatch_hairpin") == 0) return MMH;
627 else if (strcmp(ident,"mismatch_hairpin_enthalpies") == 0) return MMH_H;
628 else if (strcmp(ident,"mismatch_interior") == 0) return MMI;
629 else if (strcmp(ident,"mismatch_interior_enthalpies") == 0) return MMI_H;
630 else if (strcmp(ident,"mismatch_interior_1n") == 0) return MMI1N;
631 else if (strcmp(ident,"mismatch_interior_1n_enthalpies") == 0) return MMI1N_H;
632 else if (strcmp(ident,"mismatch_interior_23") == 0) return MMI23;
633 else if (strcmp(ident,"mismatch_interior_23_enthalpies") == 0) return MMI23_H;
634 else if (strcmp(ident,"mismatch_multi") == 0) return MMM;
635 else if (strcmp(ident,"mismatch_multi_enthalpies") == 0) return MMM_H;
636 else if (strcmp(ident,"int11") == 0) return INT11;
637 else if (strcmp(ident,"int11_enthalpies") == 0) return INT11_H;
638 else if (strcmp(ident,"int21") == 0) return INT21;
639 else if (strcmp(ident,"int21_enthalpies") == 0) return INT21_H;
640 else if (strcmp(ident,"int22") == 0) return INT22;
641 else if (strcmp(ident,"int22_enthalpies") == 0) return INT22_H;
642 else if (strcmp(ident,"dangle5")== 0) return D5;
643 else if (strcmp(ident,"dangle5_enthalpies")== 0) return D5_H;
644 else if (strcmp(ident,"dangle3")== 0) return D3;
645 else if (strcmp(ident,"dangle3_enthalpies")== 0) return D3_H;
646 else if (strcmp(ident,"ML_params")== 0) return ML;
647 else if (strcmp(ident,"NINIO") == 0) return NIN;
648 else if (strcmp(ident,"Triloops") == 0) return TRI;
649 else if (strcmp(ident,"Tetraloops") == 0) return TL;
650 else if (strcmp(ident,"Hexaloops") == 0) return HEX;
651 else if (strcmp(ident,"Misc") == 0) return MISC;
652 else if (strcmp(ident,"END") == 0) return QUIT;
656 /*---------------------------------------------------------------*/
658 PUBLIC void write_parameter_file(const char fname[]){
661 char *pnames[] = {"NP", "CG", "GC", "GU", "UG", "AU", "UA", " @"};
662 char bnames[] = "@ACGU";
663 outfp = fopen(fname, "w");
665 fprintf(stderr, "can't open file %s\n", fname);
668 fprintf(outfp,"## RNAfold parameter file v2.0\n");
670 fprintf(outfp,"\n# %s\n", settype(S));
671 fprintf(outfp,"/* CG GC GU UG AU UA @ */\n");
672 for (c=1; c<NBPAIRS+1; c++)
673 display_array(stack37[c]+1,NBPAIRS,NBPAIRS, outfp);
675 fprintf(outfp,"\n# %s\n", settype(S_H));
676 fprintf(outfp,"/* CG GC GU UG AU UA @ */\n");
677 for (c=1; c<NBPAIRS+1; c++)
678 display_array(stackdH[c]+1,NBPAIRS,NBPAIRS, outfp);
680 fprintf(outfp,"\n# %s\n", settype(MMH));
682 for (k=1; k<NBPAIRS+1; k++)
684 display_array(mismatchH37[k][i],5,5, outfp);
687 fprintf(outfp,"\n# %s\n", settype(MMH_H));
689 for (k=1; k<NBPAIRS+1; k++)
691 display_array(mismatchHdH[k][i],5,5, outfp);
695 fprintf(outfp,"\n# %s\n", settype(MMI));
697 for (k=1; k<NBPAIRS+1; k++)
699 display_array(mismatchI37[k][i],5,5, outfp);
702 fprintf(outfp,"\n# %s\n", settype(MMI_H));
704 for (k=1; k<NBPAIRS+1; k++)
706 display_array(mismatchIdH[k][i],5,5, outfp);
709 fprintf(outfp,"\n# %s\n", settype(MMI1N));
711 for (k=1; k<NBPAIRS+1; k++)
713 display_array(mismatch1nI37[k][i],5,5, outfp);
716 fprintf(outfp,"\n# %s\n", settype(MMI1N_H));
718 for (k=1; k<NBPAIRS+1; k++)
720 display_array(mismatch1nIdH[k][i],5,5, outfp);
723 fprintf(outfp,"\n# %s\n", settype(MMI23));
725 for (k=1; k<NBPAIRS+1; k++)
727 display_array(mismatch23I37[k][i],5,5, outfp);
730 fprintf(outfp,"\n# %s\n", settype(MMI23_H));
732 for (k=1; k<NBPAIRS+1; k++)
734 display_array(mismatch23IdH[k][i],5,5, outfp);
737 fprintf(outfp,"\n# %s\n", settype(MMM));
739 for (k=1; k<NBPAIRS+1; k++)
741 display_array(mismatchM37[k][i],5,5, outfp);
744 fprintf(outfp,"\n# %s\n", settype(MMM_H));
746 for (k=1; k<NBPAIRS+1; k++)
748 display_array(mismatchMdH[k][i],5,5, outfp);
751 fprintf(outfp,"\n# %s\n", settype(MME));
753 for (k=1; k<NBPAIRS+1; k++)
755 display_array(mismatchExt37[k][i],5,5, outfp);
758 fprintf(outfp,"\n# %s\n", settype(MME_H));
760 for (k=1; k<NBPAIRS+1; k++)
762 display_array(mismatchExtdH[k][i],5,5, outfp);
765 fprintf(outfp,"\n# %s\n", settype(D5));
766 fprintf(outfp,"/* @ A C G U */\n");
767 for (c=1; c<NBPAIRS+1; c++)
768 display_array(dangle5_37[c], 5, 5, outfp);
770 fprintf(outfp,"\n# %s\n", settype(D5_H));
771 fprintf(outfp,"/* @ A C G U */\n");
772 for (c=1; c<NBPAIRS+1; c++)
773 display_array(dangle5_dH[c], 5, 5, outfp);
775 fprintf(outfp,"\n# %s\n", settype(D3));
776 fprintf(outfp,"/* @ A C G U */\n");
777 for (c=1; c<NBPAIRS+1; c++)
778 display_array(dangle3_37[c], 5, 5, outfp);
780 fprintf(outfp,"\n# %s\n", settype(D3_H));
781 fprintf(outfp,"/* @ A C G U */\n");
782 for (c=1; c<NBPAIRS+1; c++)
783 display_array(dangle3_dH[c], 5, 5, outfp);
786 /* dont print "no pair" entries for interior loop arrays */
787 fprintf(outfp,"\n# %s\n", settype(INT11));
789 for (k=1; k<NBPAIRS+1; k++)
790 for (l=1; l<NBPAIRS+1; l++) {
791 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
793 display_array(int11_37[k][l][i], 5, 5, outfp);
797 fprintf(outfp,"\n# %s\n", settype(INT11_H));
799 for (k=1; k<NBPAIRS+1; k++)
800 for (l=1; l<NBPAIRS+1; l++) {
801 fprintf(outfp, "/* %2s..%2s */\n", pnames[k], pnames[l]);
803 display_array(int11_dH[k][l][i],5,5, outfp);
807 fprintf(outfp,"\n# %s\n", settype(INT21));
809 for (p1=1; p1<NBPAIRS+1; p1++)
810 for (p2=1; p2<NBPAIRS+1; p2++)
811 for (i=0; i<5; i++) {
812 fprintf(outfp, "/* %2s.%c..%2s */\n",
813 pnames[p1], bnames[i], pnames[p2]);
815 display_array(int21_37[p1][p2][i][j],5,5, outfp);
819 fprintf(outfp,"\n# %s\n", settype(INT21_H));
821 for (p1=1; p1<NBPAIRS+1; p1++)
822 for (p2=1; p2<NBPAIRS+1; p2++)
823 for (i=0; i<5; i++) {
824 fprintf(outfp, "/* %2s.%c..%2s */\n",
825 pnames[p1], bnames[i], pnames[p2]);
827 display_array(int21_dH[p1][p2][i][j],5,5, outfp);
831 fprintf(outfp,"\n# %s\n", settype(INT22));
832 { int p1, p2, i, j, k;
833 for (p1=1; p1<NBPAIRS; p1++)
834 for (p2=1; p2<NBPAIRS; p2++)
836 for (j=1; j<5; j++) {
837 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
838 pnames[p1], bnames[i], bnames[j], pnames[p2]);
840 display_array(int22_37[p1][p2][i][j][k]+1,4,5, outfp);
844 fprintf(outfp,"\n# %s\n", settype(INT22_H));
845 { int p1, p2, i, j, k;
846 for (p1=1; p1<NBPAIRS; p1++)
847 for (p2=1; p2<NBPAIRS; p2++)
849 for (j=1; j<5; j++) {
850 fprintf(outfp, "/* %2s.%c%c..%2s */\n",
851 pnames[p1], bnames[i], bnames[j], pnames[p2]);
853 display_array(int22_dH[p1][p2][i][j][k]+1,4,5, outfp);
857 fprintf(outfp,"\n# %s\n", settype(HP));
858 display_array(hairpin37, 31, 10, outfp);
860 fprintf(outfp,"\n# %s\n", settype(HP_H));
861 display_array(hairpindH, 31, 10, outfp);
863 fprintf(outfp,"\n# %s\n", settype(B));
864 display_array(bulge37, 31, 10, outfp);
866 fprintf(outfp,"\n# %s\n", settype(B_H));
867 display_array(bulgedH, 31, 10, outfp);
869 fprintf(outfp,"\n# %s\n", settype(IL));
870 display_array(internal_loop37, 31, 10, outfp);
872 fprintf(outfp,"\n# %s\n", settype(IL_H));
873 display_array(internal_loopdH, 31, 10, outfp);
875 fprintf(outfp,"\n# %s\n", settype(ML));
876 fprintf(outfp,"/* F = cu*n_unpaired + cc + ci*loop_degree (+TermAU) */\n");
877 fprintf(outfp,"/*\t cu\t cu_dH\t cc\t cc_dH\t ci\t ci_dH */\n");
878 fprintf(outfp,"\t%6d\t%6d\t%6d\t%6d\t%6d\t%6d\n", ML_BASE37, ML_BASEdH, ML_closing37, ML_closingdH, ML_intern37, ML_interndH);
880 fprintf(outfp,"\n# %s\n", settype(NIN));
881 fprintf(outfp,"/* Ninio = MIN(max, m*|n1-n2| */\n"
882 "/*\t m\t m_dH max */\n"
883 "\t%6d\t%6d\t%6d\n", ninio37, niniodH, MAX_NINIO);
885 fprintf(outfp,"\n# %s\n", settype(MISC));
886 fprintf(outfp,"/* all parameters are pairs of 'energy enthalpy' */\n");
887 fprintf(outfp,"/* DuplexInit TerminalAU LXC */\n");
888 fprintf(outfp," %6d %6d %6d %6d %3.6f %6d\n", DuplexInit37, DuplexInitdH, TerminalAU37, TerminalAUdH, lxc37, 0);
890 fprintf(outfp,"\n# %s\n", settype(HEX));
891 for (c=0; c< strlen(Hexaloops)/9; c++)
892 fprintf(outfp,"\t%.8s %6d %6d\n", Hexaloops+c*9, Hexaloop37[c], HexaloopdH[c]);
894 fprintf(outfp,"\n# %s\n", settype(TL));
895 for (c=0; c< strlen(Tetraloops)/7; c++)
896 fprintf(outfp,"\t%.6s %6d %6d\n", Tetraloops+c*7, Tetraloop37[c], TetraloopdH[c]);
898 fprintf(outfp,"\n# %s\n", settype(TRI));
899 for (c=0; c< strlen(Triloops)/6; c++)
900 fprintf(outfp,"\t%.5s %6d %6d\n", Triloops+c*6, Triloop37[c], TriloopdH[c]);
902 fprintf(outfp,"\n# %s\n", settype(QUIT));
906 PRIVATE void check_symmetry(void) {
909 for (i=0; i<=NBPAIRS; i++)
910 for (j=0; j<=NBPAIRS; j++)
911 if (stack37[i][j] != stack37[j][i])
912 fprintf(stderr, "WARNING: stacking energies not symmetric\n");
914 for (i=0; i<=NBPAIRS; i++)
915 for (j=0; j<=NBPAIRS; j++)
916 if (stackdH[i][j] != stackdH[j][i])
917 fprintf(stderr, "WARNING: stacking enthalpies not symmetric\n");
920 /* interior 1x1 loops */
921 for (i=0; i<=NBPAIRS; i++)
922 for (j=0; j<=NBPAIRS; j++)
925 if (int11_37[i][j][k][l] != int11_37[j][i][l][k])
926 fprintf(stderr, "WARNING: int11 energies not symmetric (%d,%d,%d,%d) (%d vs. %d)\n", i, j, k, l, int11_37[i][j][k][l], int11_37[j][i][l][k]);
928 for (i=0; i<=NBPAIRS; i++)
929 for (j=0; j<=NBPAIRS; j++)
932 if (int11_dH[i][j][k][l] != int11_dH[j][i][l][k])
933 fprintf(stderr, "WARNING: int11 enthalpies not symmetric\n");
935 /* interior 2x2 loops */
936 for (i=0; i<=NBPAIRS; i++)
937 for (j=0; j<=NBPAIRS; j++)
939 for (l=0; l<5; l++) {
943 if (int22_37[i][j][k][l][m][n] != int22_37[j][i][m][n][k][l])
944 fprintf(stderr, "WARNING: int22 energies not symmetric\n");
947 for (i=0; i<=NBPAIRS; i++)
948 for (j=0; j<=NBPAIRS; j++)
950 for (l=0; l<5; l++) {
954 if (int22_dH[i][j][k][l][m][n] != int22_dH[j][i][m][n][k][l])
955 fprintf(stderr, "WARNING: int22 enthalpies not symmetric: %d %d %d %d %d %d\n", i,j,k,l,m,n);
959 /* update nonstandard nucleotide/basepair involved contributions for int22 */
960 PRIVATE void update_nst(int array[NBPAIRS+1][NBPAIRS+1][5][5][5][5]){
961 int i, j, k, l, m, n;
962 int max, max2, max3, max4, max5, max6;
964 /* get maxima for one nonstandard nucleotide */
965 for (i=1; i<NBPAIRS; i++){
966 for (j=1; j<NBPAIRS; j++){
970 max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
972 max = MAX2(max, array[i][j][k][l][m][n]);
973 max2 = MAX2(max2, array[i][j][k][l][n][m]);
974 max3 = MAX2(max3, array[i][j][k][n][l][m]);
975 max4 = MAX2(max4, array[i][j][n][k][l][m]);
977 array[i][j][k][l][m][0] = max;
978 array[i][j][k][l][0][m] = max2;
979 array[i][j][k][0][l][m] = max3;
980 array[i][j][0][k][l][m] = max4;
986 /* get maxima for two nonstandard nucleotides */
987 for (i=1; i<NBPAIRS; i++){
988 for (j=1; j<NBPAIRS; j++){
991 max = max2 = max3 = max4 = max5 = max6 = -INF; /* max of {CGAU} */
993 max = MAX2(max, array[i][j][k][l][m][0]);
994 max2 = MAX2(max2, array[i][j][k][m][0][l]);
995 max3 = MAX2(max3, array[i][j][m][0][k][l]);
996 max4 = MAX2(max4, array[i][j][0][k][l][m]);
997 max5 = MAX2(max5, array[i][j][0][k][m][l]);
998 max6 = MAX2(max6, array[i][j][k][0][l][m]);
1000 array[i][j][k][l][0][0] = max;
1001 array[i][j][k][0][0][l] = max2;
1002 array[i][j][0][0][k][l] = max3;
1003 array[i][j][k][0][l][0] = max6;
1004 array[i][j][0][k][0][l] = max5;
1005 array[i][j][0][k][l][0] = max4;
1010 /* get maxima for three nonstandard nucleotides */
1011 for (i=1; i<NBPAIRS; i++){
1012 for (j=1; j<NBPAIRS; j++){
1013 for (k=1; k<5; k++){
1014 max = max2 = max3 = max4 = -INF; /* max of {CGAU} */
1015 for (l=1; l<5; l++){
1016 /* should be arbitrary where index l resides in last 3 possible locations */
1017 max = MAX2(max, array[i][j][k][l][0][0]);
1018 max2 = MAX2(max2, array[i][j][0][k][l][0]);
1019 max3 = MAX2(max3, array[i][j][0][0][k][l]);
1020 max4 = MAX2(max4, array[i][j][0][0][l][k]);
1022 array[i][j][k][0][0][0] = max;
1023 array[i][j][0][k][0][0] = max2;
1024 array[i][j][0][0][k][0] = max3;
1025 array[i][j][0][0][0][k] = max4;
1029 /* get maxima for 4 nonstandard nucleotides */
1030 for (i=1; i<NBPAIRS; i++){
1031 for (j=1; j<NBPAIRS; j++){
1032 max = -INF; /* max of {CGAU} */
1033 for (k=1; k<5; k++){
1034 max = MAX2(max, array[i][j][k][0][0][0]);
1036 array[i][j][0][0][0][0] = max;
1040 /* now compute contributions for nonstandard base pairs ... */
1041 /* first, 1 nonstandard bp */
1042 for (i=1; i<NBPAIRS; i++){
1043 for (k=0; k<5; k++){
1044 for (l=0; l<5; l++){
1045 for (m=0; m<5; m++){
1048 for(j=1;j<NBPAIRS;j++){
1049 max = MAX2(max, array[i][j][k][l][m][n]);
1050 max2 = MAX2(max2, array[j][i][k][l][m][n]);
1052 array[i][NBPAIRS][k][l][m][n] = max;
1053 array[NBPAIRS][i][k][l][m][n] = max2;
1060 /* now 2 nst base pairs */
1061 for (k=0; k<5; k++){
1062 for (l=0; l<5; l++){
1063 for (m=0; m<5; m++){
1066 for(j=1;j<NBPAIRS;j++){
1067 max = MAX2(max, array[NBPAIRS][j][k][l][m][n]);
1069 array[NBPAIRS][NBPAIRS][k][l][m][n] = max;